From a64799f06f33922493dc0dd3eb4aa2e7b5ca3753 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 6 Mar 2024 09:30:46 +0000 Subject: [PATCH 01/57] Prepare next dev version --- CITATION.cff | 4 ++-- DESCRIPTION | 2 +- NEWS.md | 2 ++ README.md | 2 +- 4 files changed, 6 insertions(+), 4 deletions(-) diff --git a/CITATION.cff b/CITATION.cff index 6222c213..417bdfaf 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -8,7 +8,7 @@ message: 'To cite package "multinma" in publications use:' type: software license: GPL-3.0-only title: 'multinma: Bayesian Network Meta-Analysis of Individual and Aggregate Data' -version: 0.6.1 +version: 0.6.1.9000 doi: 10.5281/zenodo.3904454 abstract: Network meta-analysis and network meta-regression models for aggregate data, individual patient data, and mixtures of both individual and aggregate data using @@ -28,7 +28,7 @@ preferred-citation: email: david.phillippo@bristol.ac.uk orcid: https://orcid.org/0000-0003-2672-7841 year: '2023' - notes: R package version 0.6.1 + notes: R package version 0.6.1.9000 url: https://dmphillippo.github.io/multinma/ doi: 10.5281/zenodo.3904454 repository: https://CRAN.R-project.org/package=multinma diff --git a/DESCRIPTION b/DESCRIPTION index f6441312..f6ed702e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: multinma Title: Bayesian Network Meta-Analysis of Individual and Aggregate Data -Version: 0.6.1 +Version: 0.6.1.9000 Authors@R: person(given = c("David", "M."), family = "Phillippo", diff --git a/NEWS.md b/NEWS.md index dbf9fb61..c3edd949 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,5 @@ +# multinma 0.6.1.9000 + # multinma 0.6.1 * Fix: Piecewise exponential hazard models no longer give errors during set-up. diff --git a/README.md b/README.md index 7723194a..700e6d7e 100644 --- a/README.md +++ b/README.md @@ -78,7 +78,7 @@ papers: The `multinma` package can be cited as follows: > Phillippo, D. M. (2024). *multinma: Bayesian Network Meta-Analysis of -> Individual and Aggregate Data*. R package version 0.6.1, doi: +> Individual and Aggregate Data*. R package version 0.6.1.9000, doi: > [10.5281/zenodo.3904454](https://doi.org/10.5281/zenodo.3904454). When fitting ML-NMR models, please cite the methods paper: From 1fd3d71c3148d12ae26d355cf5efff68a077c430 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 6 Mar 2024 09:40:45 +0000 Subject: [PATCH 02/57] Update actions to checkout@v4 (node.js 20) --- .github/workflows/R-CMD-check-HTML5.yaml | 2 +- .github/workflows/R-CMD-check.yaml | 2 +- .github/workflows/commit-commands.yaml | 4 ++-- .github/workflows/pkgdown.yaml | 2 +- .github/workflows/pr-commands.yaml | 4 ++-- 5 files changed, 7 insertions(+), 7 deletions(-) diff --git a/.github/workflows/R-CMD-check-HTML5.yaml b/.github/workflows/R-CMD-check-HTML5.yaml index adb1ec18..115e6039 100644 --- a/.github/workflows/R-CMD-check-HTML5.yaml +++ b/.github/workflows/R-CMD-check-HTML5.yaml @@ -16,7 +16,7 @@ jobs: R_KEEP_PKG_SOURCE: yes _R_CHECK_RD_VALIDATE_RD2HTML_: TRUE steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v4 - name: Install pdflatex run: sudo apt-get install texlive-latex-base texlive-fonts-recommended texlive-fonts-extra texlive-latex-extra diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index bdf5ae43..b70e09bf 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -27,7 +27,7 @@ jobs: _R_CHECK_DONTTEST_EXAMPLES_: false steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-pandoc@v2 diff --git a/.github/workflows/commit-commands.yaml b/.github/workflows/commit-commands.yaml index 29a4a1f6..ac249718 100644 --- a/.github/workflows/commit-commands.yaml +++ b/.github/workflows/commit-commands.yaml @@ -14,7 +14,7 @@ jobs: env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-r@v2 with: @@ -44,7 +44,7 @@ jobs: env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-pandoc@v2 diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index ddec85bf..07b6161d 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -17,7 +17,7 @@ jobs: env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-pandoc@v2 diff --git a/.github/workflows/pr-commands.yaml b/.github/workflows/pr-commands.yaml index a2dafcfd..bae4be1b 100644 --- a/.github/workflows/pr-commands.yaml +++ b/.github/workflows/pr-commands.yaml @@ -14,7 +14,7 @@ jobs: env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/pr-fetch@v2 with: @@ -51,7 +51,7 @@ jobs: env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/pr-fetch@v2 with: From 964fc5b4bd5d1f8dce13cd5ae8ad6022f854e7a3 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 6 Mar 2024 12:29:04 +0000 Subject: [PATCH 03/57] Avoid bayesplot pairs theme warnings --- R/stan_nma-class.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/stan_nma-class.R b/R/stan_nma-class.R index 95bfc346..98dc32d9 100644 --- a/R/stan_nma-class.R +++ b/R/stan_nma-class.R @@ -901,8 +901,8 @@ pairs.stan_nma <- function(x, ..., pars, include = TRUE) { condition = bayesplot::pairs_condition(nuts = "accept_stat__"), .homonyms = "first") - thm <- bayesplot::bayesplot_theme_set(theme_multinma()) + thm <- ggplot2::theme_set(theme_multinma()) out <- do.call(bayesplot::mcmc_pairs, args = args) - bayesplot::bayesplot_theme_set(thm) + ggplot2::theme_set(thm) return(out) } From 936a325d6f3935bd90795bb777ccae63d32da2e9 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Mon, 11 Mar 2024 14:01:57 +0000 Subject: [PATCH 04/57] Fix fallback formatting when crayon not installed --- R/nma_data-class.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/R/nma_data-class.R b/R/nma_data-class.R index fb29f50a..ef1dc755 100644 --- a/R/nma_data-class.R +++ b/R/nma_data-class.R @@ -213,25 +213,25 @@ subtle <- function(...) { if (require_pkg("crayon", error = FALSE)) return(crayon::silver(...)) else - return(...) + return(paste0(...)) } bold <- function(...) { if (require_pkg("crayon", error = FALSE)) return(crayon::bold(...)) else - return(...) + return(paste0(...)) } red <- function(...) { if (require_pkg("crayon", error = FALSE)) return(crayon::red(...)) else - return(...) + return(paste0(...)) } green <- function(...) { if (require_pkg("crayon", error = FALSE)) return(crayon::green(...)) else - return(...) + return(paste0(...)) } #' Convert networks to graph objects From f65d0e1c524046404e06c8b729556e0d22fac529 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Mon, 11 Mar 2024 16:46:29 +0000 Subject: [PATCH 05/57] Update news --- NEWS.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/NEWS.md b/NEWS.md index c3edd949..e9c338cd 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,8 @@ # multinma 0.6.1.9000 +* Fix: Fallback formatting used by print methods when the crayon package is not +installed now works properly, rather than giving errors. + # multinma 0.6.1 * Fix: Piecewise exponential hazard models no longer give errors during set-up. From aab2bf025884095c036b222ac187bf441ec65b74 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Mon, 11 Mar 2024 16:46:36 +0000 Subject: [PATCH 06/57] Bump dev version --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index f6ed702e..67f5d651 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: multinma Title: Bayesian Network Meta-Analysis of Individual and Aggregate Data -Version: 0.6.1.9000 +Version: 0.6.1.9001 Authors@R: person(given = c("David", "M."), family = "Phillippo", From 22b2424e453c303315f5aa1d1b8c72e62e1f4515 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Tue, 16 Apr 2024 09:32:47 +0100 Subject: [PATCH 07/57] Fix typo causing predict to fail for AgD meta-regression newdata --- DESCRIPTION | 2 +- NEWS.md | 2 ++ R/predict.R | 2 +- tests/testthat/test-example_bcg_vaccine.R | 13 +++++++++---- vignettes/example_bcg_vaccine.Rmd | 13 +++++++++---- 5 files changed, 22 insertions(+), 10 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 67f5d651..307340c5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: multinma Title: Bayesian Network Meta-Analysis of Individual and Aggregate Data -Version: 0.6.1.9001 +Version: 0.6.1.9002 Authors@R: person(given = c("David", "M."), family = "Phillippo", diff --git a/NEWS.md b/NEWS.md index e9c338cd..ea17d27a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,8 @@ * Fix: Fallback formatting used by print methods when the crayon package is not installed now works properly, rather than giving errors. +* Fix: Small bug caused `predict()` for AgD meta-regression models with new data +and `baseline_type = "response"` to fail with an error. # multinma 0.6.1 diff --git a/R/predict.R b/R/predict.R index c8fd19b7..adbb6b72 100644 --- a/R/predict.R +++ b/R/predict.R @@ -1354,7 +1354,7 @@ predict.stan_nma <- function(object, ..., # Convert to linear predictor scale if baseline_type = "response" if (any(baseline_type == "response")) { - mu[ , , baseline_type == "response"] <- link_fun(mu[ , , baseline_type = "response"], link = object$link) + mu[ , , baseline_type == "response"] <- link_fun(mu[ , , baseline_type == "response"], link = object$link) } # Convert to samples on network ref trt if trt_ref given diff --git a/tests/testthat/test-example_bcg_vaccine.R b/tests/testthat/test-example_bcg_vaccine.R index 857a2d3d..ef7c091b 100644 --- a/tests/testthat/test-example_bcg_vaccine.R +++ b/tests/testthat/test-example_bcg_vaccine.R @@ -299,16 +299,21 @@ test_that("Relative effects and predict work with data.frame", { new <- tibble::tibble(latitude = seq(10, 50, by = 10), label = paste0(latitude, "\u00B0 latitude")) expect_identical(relative_effects(bcg_fit_lat, newdata = new, study = label), relative_effects(bcg_fit_lat, newdata = as.data.frame(new), study = label)) + # For predict() we need to account for the random baseline sample - # expect_identical(withr::with_seed(predict(bcg_fit_lat, newdata = new, study = label, - # baseline = distr(qnorm, mean = -2, sd = 0.1)), seed = 1234), - # withr::with_seed(predict(bcg_fit_lat, newdata = as.data.frame(new), study = label, - # baseline = distr(qnorm, mean = -2, sd = 0.1)), seed = 1234)) qcons <- function(p, cons = 0) {cons} + expect_identical(predict(bcg_fit_lat, newdata = new, study = label, baseline = distr(qcons, -2)), predict(bcg_fit_lat, newdata = as.data.frame(new), study = label, baseline = distr(qcons, -2))) + + expect_identical(predict(bcg_fit_lat, newdata = new, study = label, + baseline = distr(qcons, 0.5), + baseline_type = "response"), + predict(bcg_fit_lat, newdata = as.data.frame(new), study = label, + baseline = distr(qcons, 0.5), + baseline_type = "response")) }) test_that("Predictions using network baselines are correct", { diff --git a/vignettes/example_bcg_vaccine.Rmd b/vignettes/example_bcg_vaccine.Rmd index df6d3d64..e562b837 100644 --- a/vignettes/example_bcg_vaccine.Rmd +++ b/vignettes/example_bcg_vaccine.Rmd @@ -360,16 +360,21 @@ test_that("Relative effects and predict work with data.frame", { new <- tibble::tibble(latitude = seq(10, 50, by = 10), label = paste0(latitude, "\u00B0 latitude")) expect_identical(relative_effects(bcg_fit_lat, newdata = new, study = label), relative_effects(bcg_fit_lat, newdata = as.data.frame(new), study = label)) + # For predict() we need to account for the random baseline sample - # expect_identical(withr::with_seed(predict(bcg_fit_lat, newdata = new, study = label, - # baseline = distr(qnorm, mean = -2, sd = 0.1)), seed = 1234), - # withr::with_seed(predict(bcg_fit_lat, newdata = as.data.frame(new), study = label, - # baseline = distr(qnorm, mean = -2, sd = 0.1)), seed = 1234)) qcons <- function(p, cons = 0) {cons} + expect_identical(predict(bcg_fit_lat, newdata = new, study = label, baseline = distr(qcons, -2)), predict(bcg_fit_lat, newdata = as.data.frame(new), study = label, baseline = distr(qcons, -2))) + + expect_identical(predict(bcg_fit_lat, newdata = new, study = label, + baseline = distr(qcons, 0.5), + baseline_type = "response"), + predict(bcg_fit_lat, newdata = as.data.frame(new), study = label, + baseline = distr(qcons, 0.5), + baseline_type = "response")) }) test_that("Predictions using network baselines are correct", { From c399888dd901cc96663ed1f8f81c646b09a52b37 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Tue, 13 Feb 2024 10:39:01 +0000 Subject: [PATCH 08/57] Start marginal_effects() --- R/marginal_effects.R | 87 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) create mode 100644 R/marginal_effects.R diff --git a/R/marginal_effects.R b/R/marginal_effects.R new file mode 100644 index 00000000..d2ee2127 --- /dev/null +++ b/R/marginal_effects.R @@ -0,0 +1,87 @@ +#' Marginal treatment effects +#' +#' Generate population-average marginal treatment effects. These are formed from +#' population-average absolute predictions, so this function is a wrapper around +#' [predict.stan_nma()]. +#' +#' @param object A `stan_nma` object created by [nma()]. +#' @param ... Arguments passed to [predict.stan_nma()], for example to specify +#' the covariate distribution and baseline risk for a target population. +#' @param mtype The type of marginal effect to construct from the average +#' absolute effects, either `"difference"` (the default) for a difference of +#' absolute effects such as a risk difference, `"ratio"` for a ratio of +#' absolute effects such as a risk ratio, or `"link"` for a difference on the +#' scale of the link function used in fitting the model such as a marginal log +#' odds ratio. +#' @param all_contrasts Logical, generate estimates for all contrasts (`TRUE`), +#' or just the "basic" contrasts against the network reference treatment +#' (`FALSE`)? Default `FALSE`. +#' @param trt_ref Reference treatment to construct relative effects against, if +#' `all_contrasts = FALSE`. By default, relative effects will be against the +#' network reference treatment. Coerced to character string. +#' @param probs Numeric vector of quantiles of interest to present in computed +#' summary, default `c(0.025, 0.25, 0.5, 0.75, 0.975)` +#' @param predictive_distribution Logical, when a random effects model has been +#' fitted, should the predictive distribution for marginal effects in a new +#' study be returned? Default `FALSE`. +#' @param summary Logical, calculate posterior summaries? Default `TRUE`. +#' +#' @return +#' @export +#' +#' @examples +marginal_effects <- function(object, + ..., + mtype = c("difference", "ratio", "link"), + all_contrasts = FALSE, trt_ref = NULL, + probs = c(0.025, 0.25, 0.5, 0.75, 0.975), + predictive_distribution = FALSE, + summary = TRUE) { + + # Checks + if (!inherits(object, "stan_nma")) abort("Expecting a `stan_nma` object, as returned by nma().") + + mtype <- rlang::arg_match(mtype) + + if (!rlang::is_bool(all_contrasts)) + abort("`all_contrasts` should be TRUE or FALSE.") + + if (!is.null(trt_ref)) { + if (all_contrasts) { + warn("Ignoring `trt_ref` when all_contrasts = TRUE.") + trt_ref <- NULL + } else { + if (length(trt_ref) > 1) abort("`trt_ref` must be length 1.") + trt_ref <- as.character(trt_ref) + lvls_trt <- levels(object$network$treatments) + if (! trt_ref %in% lvls_trt) + abort(sprintf("`trt_ref` does not match a treatment in the network.\nSuitable values are: %s", + ifelse(length(lvls_trt) <= 5, + paste0(lvls_trt, collapse = ", "), + paste0(paste0(lvls_trt[1:5], collapse = ", "), ", ...")))) + } + } + + if (!rlang::is_bool(summary)) + abort("`summary` should be TRUE or FALSE.") + + check_probs(probs) + + if(!rlang::is_bool(predictive_distribution)) + abort("`predictive_distribution` should be TRUE or FALSE") + if (predictive_distribution && x$trt_effects != "random") predictive_distribution <- FALSE + + # Cannot produce marginal effects for inconsistency models + if (x$consistency != "consistency") + abort(glue::glue("Cannot produce marginal effects under inconsistency '{x$consistency}' model.")) + + # Get network reference treatment + nrt <- levels(x$network$treatments)[1] + + # Call predict to get absolute predictions + pred <- predict(object, + type = "response", level = "aggregate", + summary = FALSE, + predictive_distribution = predictive_distribution, + ...) +} From c581a698a5be0afcfb77ab2a9b6f66ce03b2e507 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Tue, 13 Feb 2024 11:17:40 +0000 Subject: [PATCH 09/57] Always output predict metadata, even if summary=FALSE --- R/predict.R | 134 +++++++++++++++++++++++----------------------------- 1 file changed, 60 insertions(+), 74 deletions(-) diff --git a/R/predict.R b/R/predict.R index adbb6b72..12a516f2 100644 --- a/R/predict.R +++ b/R/predict.R @@ -761,65 +761,54 @@ predict.stan_nma <- function(object, ..., } # Produce nma_summary - if (summary) { - pred_summary <- summary_mcmc_array(pred_array, probs) - - if (object$likelihood == "ordered") { - pred_summary <- tibble::add_column(pred_summary, - .trt = rep(preddat$.trt, each = n_cc), - .category = rep(l_cc, times = nrow(preddat)), - .before = 1) - - } else if (object$likelihood %in% valid_lhood$survival) { - if (type %in% c("survival", "hazard", "cumhaz")) { - preddat <- tidyr::unnest(preddat, cols = ".time") - pred_summary <- tibble::add_column(pred_summary, - .trt = preddat$.trt, - .time = preddat$.time, - .before = 1) - } else if (type == "rmst") { - pred_summary <- tibble::add_column(pred_summary, - .trt = preddat$.trt, - .time = preddat$.time, - .before = 1) - } else if (type == "quantile") { - pred_summary <- tibble::add_column(pred_summary, - .trt = rep(preddat$.trt, each = length(quantiles)), - .quantile = rep(quantiles, times = nrow(preddat)), - .before = 1) - } else { - pred_summary <- tibble::add_column(pred_summary, - .trt = preddat$.trt, - .before = 1) - } + if (object$likelihood == "ordered") { + pred_meta <- tibble::tibble(.trt = rep(preddat$.trt, each = n_cc), + .category = rep(l_cc, times = nrow(preddat))) + } else if (object$likelihood %in% valid_lhood$survival) { + if (type %in% c("survival", "hazard", "cumhaz")) { + preddat <- tidyr::unnest(preddat, cols = ".time") + pred_meta <- tibble::tibble(.trt = preddat$.trt, + .time = preddat$.time) + } else if (type == "rmst") { + pred_meta <- tibble::tibble(.trt = preddat$.trt, + .time = preddat$.time) + } else if (type == "quantile") { + pred_meta <- tibble::tibble(.trt = rep(preddat$.trt, each = length(quantiles)), + .quantile = rep(quantiles, times = nrow(preddat))) } else { - pred_summary <- tibble::add_column(pred_summary, - .trt = preddat$.trt, - .before = 1) + pred_meta <- tibble::tibble(.trt = preddat$.trt) } - if (is.null(baseline)) { - if (object$likelihood == "ordered") { - pred_summary <- tibble::add_column(pred_summary, - .study = rep(preddat$.study, each = n_cc), - .before = 1) - } else if (type == "quantile") { - pred_summary <- tibble::add_column(pred_summary, - .study = rep(preddat$.study, each = length(quantiles)), - .before = 1) - } else { - pred_summary <- tibble::add_column(pred_summary, - .study = preddat$.study, - .before = 1) - } + } else { + pred_meta <- tibble::tibble(.trt = preddat$.trt) + } + + if (is.null(baseline)) { + if (object$likelihood == "ordered") { + pred_meta <- tibble::add_column(pred_meta, + .study = rep(preddat$.study, each = n_cc), + .before = 1) + } else if (type == "quantile") { + pred_meta <- tibble::add_column(pred_meta, + .study = rep(preddat$.study, each = length(quantiles)), + .before = 1) + } else { + pred_meta <- tibble::add_column(pred_meta, + .study = preddat$.study, + .before = 1) } + } - out <- list(summary = pred_summary, sims = pred_array) + if (summary) { + pred_summary <- summary_mcmc_array(pred_array, probs) + pred_summary <- dplyr::bind_cols(pred_meta, pred_summary) } else { - out <- list(sims = pred_array) + pred_summary <- pred_meta } + out <- list(summary = pred_summary, sims = pred_array) + # With regression model ------------------------------------------------------ } else { @@ -1922,35 +1911,32 @@ predict.stan_nma <- function(object, ..., } # Produce nma_summary + if (object$likelihood == "ordered") { + pred_meta <- tibble::tibble(.study = rep(preddat$.study, each = n_cc), + .trt = rep(preddat$.trt, each = n_cc), + .category = rep(l_cc, times = nrow(preddat))) + } else if (is_surv) { + pred_meta <- tibble::tibble(.study = outdat$.study, + .trt = outdat$.trt) + if (type %in% c("survival", "hazard", "cumhaz", "rmst")) { + pred_meta <- tibble::add_column(pred_meta, .time = outdat$.time, .after = ".trt") + } else if (type == "quantile") { + pred_meta <- tibble::add_column(pred_meta, .quantile = outdat$.quantile, .after = ".trt") + } + } else { + pred_meta <- tibble::tibble(.study = preddat$.study, + .trt = preddat$.trt) + } + if (summary) { pred_summary <- summary_mcmc_array(pred_array, probs) - if (object$likelihood == "ordered") { - pred_summary <- tibble::add_column(pred_summary, - .study = rep(preddat$.study, each = n_cc), - .trt = rep(preddat$.trt, each = n_cc), - .category = rep(l_cc, times = nrow(preddat)), - .before = 1) - } else if (is_surv) { - pred_summary <- tibble::add_column(pred_summary, - .study = outdat$.study, - .trt = outdat$.trt, - .before = 1) - if (type %in% c("survival", "hazard", "cumhaz", "rmst")) { - pred_summary <- tibble::add_column(pred_summary, .time = outdat$.time, .after = ".trt") - } else if (type == "quantile") { - pred_summary <- tibble::add_column(pred_summary, .quantile = outdat$.quantile, .after = ".trt") - } - } else { - pred_summary <- tibble::add_column(pred_summary, - .study = preddat$.study, - .trt = preddat$.trt, - .before = 1) - } - out <- list(summary = pred_summary, sims = pred_array) + pred_summary <- dplyr::bind_cols(pred_meta, pred_summary) } else { - out <- list(sims = pred_array) + pred_summary <- pred_meta } + out <- list(summary = pred_summary, sims = pred_array) + } if (summary) { From d1aa78ec46a396f37749ae9214b655b7910825c4 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Tue, 13 Feb 2024 17:48:38 +0000 Subject: [PATCH 10/57] Flesh out marginal_effects --- R/marginal_effects.R | 158 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 147 insertions(+), 11 deletions(-) diff --git a/R/marginal_effects.R b/R/marginal_effects.R index d2ee2127..6cd2d10e 100644 --- a/R/marginal_effects.R +++ b/R/marginal_effects.R @@ -60,6 +60,8 @@ marginal_effects <- function(object, paste0(lvls_trt, collapse = ", "), paste0(paste0(lvls_trt[1:5], collapse = ", "), ", ...")))) } + } else { + trt_ref <- levels(object$network$treatments)[1] } if (!rlang::is_bool(summary)) @@ -69,19 +71,153 @@ marginal_effects <- function(object, if(!rlang::is_bool(predictive_distribution)) abort("`predictive_distribution` should be TRUE or FALSE") - if (predictive_distribution && x$trt_effects != "random") predictive_distribution <- FALSE + if (predictive_distribution && object$trt_effects != "random") predictive_distribution <- FALSE # Cannot produce marginal effects for inconsistency models - if (x$consistency != "consistency") - abort(glue::glue("Cannot produce marginal effects under inconsistency '{x$consistency}' model.")) - - # Get network reference treatment - nrt <- levels(x$network$treatments)[1] + if (object$consistency != "consistency") + abort(glue::glue("Cannot produce marginal effects under inconsistency '{object$consistency}' model.")) # Call predict to get absolute predictions - pred <- predict(object, - type = "response", level = "aggregate", - summary = FALSE, - predictive_distribution = predictive_distribution, - ...) + if (object$likelihood %in% valid_lhood$survival) { + pred <- predict(object, + level = "aggregate", + summary = FALSE, + predictive_distribution = predictive_distribution, + ...) + } else { + pred <- predict(object, + type = "response", level = "aggregate", + summary = FALSE, + predictive_distribution = predictive_distribution, + ...) + } + + # Set up output metadata + pred_meta <- pred$summary + out_meta <- pred_meta + + vars <- intersect(c(".study", ".time", ".quantile", ".category"), colnames(pred_meta)) + + if (!all_contrasts) { + out_meta <- dplyr::filter(out_meta, .data$.trt != trt_ref) %>% + dplyr::mutate(.trtb = .data$.trt, .trta = trt_ref) %>% + dplyr::select(-".trt", ".study", ".trtb", ".trta", dplyr::everything()) + } else { + contrs <- utils::combn(object$network$treatments, 2) + trtb <- contrs[2, ] + trta <- contrs[1, ] + out_meta <- dplyr::distinct(out_meta, dplyr::pick(dplyr::all_of(vars))) %>% + dplyr::mutate(.trtb = list(trtb), .trta = list(trta)) %>% + tidyr::unnest(cols = c(.data$.trtb, .data$.trta)) %>% + dplyr::select(".study", ".trtb", ".trta", dplyr::everything()) %>% + dplyr::arrange(".study", ".trta", ".trtb", dplyr::pick(dplyr::everything())) + } + + pred_meta$id <- 1:nrow(pred_meta) + out_meta$id <- 1:nrow(out_meta) + + # Output parameter names + pnames <- paste0("D[", + if (dplyr::n_distinct(out_meta$.study) > 1) paste0(out_meta$.study, ": ") else character(), + out_meta$.trtb, + if (all_contrasts) paste0(" vs. ", out_meta$.trta) else character(), + if (".time" %in% vars) paste0(", ", rep(1:dplyr::n_distinct(out_meta$.time), times = dplyr::n_distinct(out_meta$.study, out_meta$.trtb, out_meta$.trta))) + else if (".quantile" %in% vars) paste0(", ", out_meta$.quantile) + else if (".category" %in% vars) paste0(", ", out_meta$.category) + else character() + , "]") + + # Set up output array + d_out <- dim(pred$sims) + d_out[[3]] <- nrow(out_meta) + dn_out <- list(iterations = NULL, chains = NULL, parameters = pnames) + out_array <- array(NA_real_, dim = d_out, dimnames = dn_out) + + # Set contrast function + if (mtype == "difference") { + cf <- function(b, a) b - a + } else if (mtype == "ratio") { + cf <- function(b, a) b / a + } else if (mtype == "link") { + cf <- function(b, a) link_fun(b, link = object$link) - link_fun(a, link = object$link) + } + + cfv <- function(x) cf(x[2], x[1]) + mk_contr <- function(x) combn(x, 2, FUN = cfv) + + # Calculate contrasts, looping over groups + grps <- dplyr::distinct(out_meta, dplyr::pick(dplyr::all_of(vars))) + + for (i in 1:nrow(grps)) { + grpi <- grps[i, ] + predi <- dplyr::semi_join(pred_meta, grpi, by = vars)$id + outi <- dplyr::semi_join(out_meta, grpi, by = vars)$id + + if (all_contrasts) { + out_array[ , , outi] <- aperm(apply(pred$sims[ , , predi], MARGIN = 1:2, FUN = mk_contr), c(2, 3, 1)) + } else { + predi_ref <- dplyr::semi_join(pred_meta, dplyr::mutate(grpi, .trt = trt_ref), by = c(vars, ".trt"))$id + predi_nonref <- setdiff(predi, predi_ref) + out_array[ , , outi] <- sweep(pred$sims[ , , predi_nonref, drop = FALSE], MARGIN = 1:2, STATS = pred$sims[ , , predi_ref], FUN = cf) + } + } + + out_meta <- dplyr::select(out_meta, -"id") + + # Summarise + if (summary) { + out_summary <- summary_mcmc_array(out_array, probs) + out_summary <- dplyr::bind_cols(out_meta, out_summary) + } else { + out_summary <- out_meta + } + + out <- list(summary = out_summary, sims = out_array) + + # Set attributes + if (summary) { + attr(out, "xlab") <- "Treatment" + + if (object$likelihood %in% valid_lhood$survival) { + dots <- rlang::enquos(...) + if (rlang::has_name(dots, "type")) type <- rlang::eval_tidy(dots$type) + else type <- "survival" + } else if (mtype == "link") { + type <- "link" + } else { + type <- "response" + } + + if (mtype == "difference") { + attr(out, "ylab") <- paste0("Marginal ", + get_scale_name(likelihood = object$likelihood, + link = object$link, + measure = "absolute", + type = type), + " Difference") + } else if (mtype == "ratio") { + attr(out, "ylab") <- paste0("Marginal ", + get_scale_name(likelihood = object$likelihood, + link = object$link, + measure = "absolute", + type = type), + " Ratio") + } else if (mtype == "link") { + attr(out, "ylab") <- paste0("Marginal ", + get_scale_name(likelihood = object$likelihood, + link = object$link, + measure = "relative", + type = type)) + } + + if (object$likelihood == "ordered") { + class(out) <- c("ordered_nma_summary", "nma_summary") + } else if (type %in% c("survival", "hazard", "cumhaz")) { + class(out) <- c("surv_nma_summary", "nma_summary") + attr(out, "xlab") <- "Time" + } else { + class(out) <- "nma_summary" + } + } + return(out) } From e9a879c7eb0e40fef96e65209ae15a724639df45 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Tue, 13 Feb 2024 17:52:17 +0000 Subject: [PATCH 11/57] Set times_seq automatically within plot for marginal_effects too --- R/predict.R | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/R/predict.R b/R/predict.R index 12a516f2..7e849f89 100644 --- a/R/predict.R +++ b/R/predict.R @@ -2012,7 +2012,9 @@ predict.stan_nma_surv <- function(object, times = NULL, abort("`times_seq` must be a single positive integer.") # Set times_seq by default if called within plot() - if (is.null(times_seq) && type %in% c("survival", "hazard", "cumhaz") && deparse(sys.call(-2)[[1]]) == "plot") + if (is.null(times_seq) && type %in% c("survival", "hazard", "cumhaz") && + (deparse(sys.call(-2)[[1]]) == "plot" || + deparse(sys.call(-2)[[1]]) == "marginal_effects" && deparse(sys.call(-3)[[1]]) == "plot")) times_seq <- 50 # Other checks (including times, aux) in predict.stan_nma() From fa8f9243b42346fbda65728d35e7b505da70ba89 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 14 Feb 2024 14:44:48 +0000 Subject: [PATCH 12/57] Fix time labeling --- R/marginal_effects.R | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/R/marginal_effects.R b/R/marginal_effects.R index 6cd2d10e..8ec53023 100644 --- a/R/marginal_effects.R +++ b/R/marginal_effects.R @@ -116,12 +116,18 @@ marginal_effects <- function(object, pred_meta$id <- 1:nrow(pred_meta) out_meta$id <- 1:nrow(out_meta) + if (".time" %in% vars) { + time_id <- dplyr::group_by(out_meta, .data$.study, .data$.trta, .data$.trtb) %>% + dplyr::mutate(time_id = 1:dplyr::n()) %>% + dplyr::pull("time_id") + } + # Output parameter names pnames <- paste0("D[", if (dplyr::n_distinct(out_meta$.study) > 1) paste0(out_meta$.study, ": ") else character(), out_meta$.trtb, if (all_contrasts) paste0(" vs. ", out_meta$.trta) else character(), - if (".time" %in% vars) paste0(", ", rep(1:dplyr::n_distinct(out_meta$.time), times = dplyr::n_distinct(out_meta$.study, out_meta$.trtb, out_meta$.trta))) + if (".time" %in% vars) paste0(", ", time_id) else if (".quantile" %in% vars) paste0(", ", out_meta$.quantile) else if (".category" %in% vars) paste0(", ", out_meta$.category) else character() From f1024a63aa2357e73f470a4d5ba24f0f411f478b Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 14 Feb 2024 14:46:03 +0000 Subject: [PATCH 13/57] Avoid namespace warnings --- R/marginal_effects.R | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/R/marginal_effects.R b/R/marginal_effects.R index 8ec53023..3dd88341 100644 --- a/R/marginal_effects.R +++ b/R/marginal_effects.R @@ -79,17 +79,17 @@ marginal_effects <- function(object, # Call predict to get absolute predictions if (object$likelihood %in% valid_lhood$survival) { - pred <- predict(object, - level = "aggregate", - summary = FALSE, - predictive_distribution = predictive_distribution, - ...) + pred <- stats::predict(object, + level = "aggregate", + summary = FALSE, + predictive_distribution = predictive_distribution, + ...) } else { - pred <- predict(object, - type = "response", level = "aggregate", - summary = FALSE, - predictive_distribution = predictive_distribution, - ...) + pred <- stats::predict(object, + type = "response", level = "aggregate", + summary = FALSE, + predictive_distribution = predictive_distribution, + ...) } # Set up output metadata @@ -149,7 +149,7 @@ marginal_effects <- function(object, } cfv <- function(x) cf(x[2], x[1]) - mk_contr <- function(x) combn(x, 2, FUN = cfv) + mk_contr <- function(x) utils::combn(x, 2, FUN = cfv) # Calculate contrasts, looping over groups grps <- dplyr::distinct(out_meta, dplyr::pick(dplyr::all_of(vars))) From 21a2d178119fc21028cf7f68f59c001743d87e0e Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Fri, 16 Feb 2024 15:03:57 +0000 Subject: [PATCH 14/57] Match predict in metadata format --- R/marginal_effects.R | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/R/marginal_effects.R b/R/marginal_effects.R index 3dd88341..958e7879 100644 --- a/R/marginal_effects.R +++ b/R/marginal_effects.R @@ -99,9 +99,7 @@ marginal_effects <- function(object, vars <- intersect(c(".study", ".time", ".quantile", ".category"), colnames(pred_meta)) if (!all_contrasts) { - out_meta <- dplyr::filter(out_meta, .data$.trt != trt_ref) %>% - dplyr::mutate(.trtb = .data$.trt, .trta = trt_ref) %>% - dplyr::select(-".trt", ".study", ".trtb", ".trta", dplyr::everything()) + out_meta <- dplyr::filter(out_meta, .data$.trt != trt_ref) } else { contrs <- utils::combn(object$network$treatments, 2) trtb <- contrs[2, ] @@ -109,23 +107,23 @@ marginal_effects <- function(object, out_meta <- dplyr::distinct(out_meta, dplyr::pick(dplyr::all_of(vars))) %>% dplyr::mutate(.trtb = list(trtb), .trta = list(trta)) %>% tidyr::unnest(cols = c(.data$.trtb, .data$.trta)) %>% - dplyr::select(".study", ".trtb", ".trta", dplyr::everything()) %>% - dplyr::arrange(".study", ".trta", ".trtb", dplyr::pick(dplyr::everything())) + dplyr::select(dplyr::any_of(".study"), ".trtb", ".trta", dplyr::everything()) %>% + dplyr::arrange(dplyr::pick(dplyr::any_of(".study")), ".trta", ".trtb", dplyr::pick(dplyr::everything())) } pred_meta$id <- 1:nrow(pred_meta) out_meta$id <- 1:nrow(out_meta) if (".time" %in% vars) { - time_id <- dplyr::group_by(out_meta, .data$.study, .data$.trta, .data$.trtb) %>% + time_id <- dplyr::group_by(out_meta, .data$.study, dplyr::pick(dplyr::any_of(c(".trt", ".trta", ".trtb")))) %>% dplyr::mutate(time_id = 1:dplyr::n()) %>% dplyr::pull("time_id") } # Output parameter names - pnames <- paste0("D[", - if (dplyr::n_distinct(out_meta$.study) > 1) paste0(out_meta$.study, ": ") else character(), - out_meta$.trtb, + pnames <- paste0("marg[", + if (rlang::has_name(out_meta, ".study")) paste0(out_meta$.study, ": ") else character(), + if (all_contrasts) out_meta$.trtb else out_meta$.trt, if (all_contrasts) paste0(" vs. ", out_meta$.trta) else character(), if (".time" %in% vars) paste0(", ", time_id) else if (".quantile" %in% vars) paste0(", ", out_meta$.quantile) From 29855797952c3a272d7525a74ceaa6cd5c4e688e Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Fri, 16 Feb 2024 15:04:10 +0000 Subject: [PATCH 15/57] Add progress bar to predict --- R/predict.R | 40 +++++++++++++++++++++++++++++++++++++--- 1 file changed, 37 insertions(+), 3 deletions(-) diff --git a/R/predict.R b/R/predict.R index 7e849f89..536fa70d 100644 --- a/R/predict.R +++ b/R/predict.R @@ -104,6 +104,11 @@ #' fitted, should the predictive distribution for absolute effects in a new #' study be returned? Default `FALSE`. #' @param summary Logical, calculate posterior summaries? Default `TRUE`. +#' @param progress Logical, display progress for potentially long-running +#' calculations? Population-average predictions from ML-NMR models are +#' computationally intensive, especially for survival outcomes. Currently the +#' default is to display progress only when running interactively and +#' producing predictions for a survival ML-NMR model. #' #' @details #' # Aggregate-level predictions from IPD NMA and ML-NMR models @@ -281,10 +286,13 @@ predict.stan_nma <- function(object, ..., baseline_level = c("individual", "aggregate"), probs = c(0.025, 0.25, 0.5, 0.75, 0.975), predictive_distribution = FALSE, - summary = TRUE) { + summary = TRUE, + progress = FALSE) { # Checks if (!inherits(object, "stan_nma")) abort("Expecting a `stan_nma` object, as returned by nma().") + if (!rlang::is_bool(progress)) abort("`progress` must be a single logical value, TRUE/FALSE.") + is_surv <- inherits(object, "stan_nma_surv") # Survival flag if (!is_surv) type <- rlang::arg_match(type) @@ -1556,6 +1564,16 @@ predict.stan_nma <- function(object, ..., l_cc <- stringr::str_replace(dimnames(cc)[[3]], "^cc\\[(.+)\\]$", "\\1") } + studies <- levels(forcats::fct_drop(preddat$.study)) + n_studies <- length(studies) + treatments <- levels(forcats::fct_drop(preddat$.trt)) + n_trt <- length(treatments) + + if (progress) { + pb <- utils::txtProgressBar(max = n_studies * n_trt, style = 3, width = min(100, getOption("width"))) + on.exit(close(pb)) + } + # Make prediction arrays if (is_surv) { # Handle survival models separately - produce predictions study by study @@ -1631,6 +1649,8 @@ predict.stan_nma <- function(object, ..., } for (trt in 1:n_trt) { + #if (progress) cglue("\rCalculating predictions for {studies[s]}, {treatments[trt]} [{(s-1)*n_trt + trt} / {n_studies * n_trt}]", sep = "") + # Collapse preddat by unique rows for efficiency collapse_by <- setdiff(colnames(preddat), c(".time", ".obs_id", ".sample_size")) @@ -1775,6 +1795,7 @@ predict.stan_nma <- function(object, ..., pred_array[ , , outdat$.study == studies[s] & outdat$.trt == treatments[trt]] <- s_pred_array } + if (progress) utils::setTxtProgressBar(pb, (s-1)*n_trt + trt) } } @@ -1810,6 +1831,8 @@ predict.stan_nma <- function(object, ..., pred_array <- inverse_link(pred_array, link = object$link) } + if (progress) utils::setTxtProgressBar(pb, n_studies * n_trt) + } else { # Predictions aggregated over each population # Produce aggregated predictions study by study - more memory efficient @@ -1843,6 +1866,8 @@ predict.stan_nma <- function(object, ..., for (s in 1:n_studies) { + # if (progress) cglue("\rCalculating predictions for {studies[s]} [{s} / {n_studies}]", sep = "") + # Study select ss <- preddat$.study == studies[s] @@ -1905,11 +1930,15 @@ predict.stan_nma <- function(object, ..., pred_array[ , , outdat$.study == studies[s]] <- tcrossprod_mcmc_array(s_pred_array, X_weighted_mean) + if (progress) utils::setTxtProgressBar(pb, s * n_trt) + } preddat <- dplyr::distinct(preddat, .data$.study, .data$.trt) } + # if (progress) cat("\rDone!") + # Produce nma_summary if (object$likelihood == "ordered") { pred_meta <- tibble::tibble(.study = rep(preddat$.study, each = n_cc), @@ -1954,6 +1983,9 @@ predict.stan_nma <- function(object, ..., class(out) <- "nma_summary" } } + + if (progress) cat("\r") + return(out) } @@ -1998,7 +2030,8 @@ predict.stan_nma_surv <- function(object, times = NULL, times_seq = NULL, probs = c(0.025, 0.25, 0.5, 0.75, 0.975), predictive_distribution = FALSE, - summary = TRUE) { + summary = TRUE, + progress = interactive()) { type <- rlang::arg_match(type) times <- rlang::enquo(times) @@ -2025,7 +2058,8 @@ predict.stan_nma_surv <- function(object, times = NULL, quantiles = quantiles, times_seq = times_seq, baseline_level = "individual", - baseline_type = "link") + baseline_type = "link", + progress = progress) } #' Produce survival predictions from arrays of linear predictors and auxiliary parameters From 96794519f8e312ca72b973bbaf1e444ccc9b0539 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Fri, 16 Feb 2024 15:04:22 +0000 Subject: [PATCH 16/57] Add marginal_effects() examples --- R/marginal_effects.R | 87 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 86 insertions(+), 1 deletion(-) diff --git a/R/marginal_effects.R b/R/marginal_effects.R index 958e7879..c82e97b9 100644 --- a/R/marginal_effects.R +++ b/R/marginal_effects.R @@ -29,7 +29,92 @@ #' @return #' @export #' -#' @examples +#' @examples ## Smoking cessation +#' @template ex_smoking_nma_re_example +#' @examples \donttest{ +#' # Marginal risk difference in each study population in the network +#' marginal_effects(smk_fit_RE, mtype = "difference") +#' +#' # Since there are no covariates in the model, the marginal and conditional +#' # (log) odds ratios here coincide +#' marginal_effects(smk_fit_RE, mtype = "link") +#' relative_effects(smk_fit_RE) +#' +#' # Marginal risk differences in a population with 67 observed events out of +#' # 566 individuals on No Intervention, corresponding to a Beta(67, 566 - 67) +#' # distribution on the baseline probability of response +#' (smk_rd_RE <- marginal_effects(smk_fit_RE, +#' baseline = distr(qbeta, 67, 566 - 67), +#' baseline_type = "response", +#' mtype = "difference")) +#' plot(smk_rd_RE) +#' } +#' +#' ## Plaque psoriasis ML-NMR +#' @template ex_plaque_psoriasis_mlnmr_example +#' @examples \donttest{ +#' # Population-average marginal probit differences in each study in the network +#' (pso_marg <- marginal_effects(pso_fit, mtype = "link")) +#' plot(pso_marg, ref_line = c(0, 1)) +#' +#' # Population-average marginal probit differences in a new target population, +#' # with means and SDs or proportions given by +#' new_agd_int <- data.frame( +#' bsa_mean = 0.6, +#' bsa_sd = 0.3, +#' prevsys = 0.1, +#' psa = 0.2, +#' weight_mean = 10, +#' weight_sd = 1, +#' durnpso_mean = 3, +#' durnpso_sd = 1 +#' ) +#' +#' # We need to add integration points to this data frame of new data +#' # We use the weighted mean correlation matrix computed from the IPD studies +#' new_agd_int <- add_integration(new_agd_int, +#' durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd), +#' prevsys = distr(qbern, prob = prevsys), +#' bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd), +#' weight = distr(qgamma, mean = weight_mean, sd = weight_sd), +#' psa = distr(qbern, prob = psa), +#' cor = pso_net$int_cor, +#' n_int = 64) +#' +#' # Population-average marginal probit differences of achieving PASI 75 in this +#' # target population, given a Normal(-1.75, 0.08^2) distribution on the +#' # baseline probit-probability of response on Placebo (at the reference levels +#' # of the covariates), are given by +#' (pso_marg_new <- marginal_effects(pso_fit, +#' mtype = "link", +#' newdata = new_agd_int, +#' baseline = distr(qnorm, -1.75, 0.08))) +#' plot(pso_marg_new) +#' } +#' +#' ## Progression free survival with newly-diagnosed multiple myeloma +#' @template ex_ndmm_example +#' @examples \donttest{ +#' # We can produce a range of marginal effects from models with survival +#' # outcomes, specified with the mtype and type arguments. For example: +#' +#' # Marginal survival probability difference at 5 years, all contrasts +#' marginal_effects(ndmm_fit, type = "survival", mtype = "difference", +#' times = 5, all_contrasts = TRUE) +#' +#' # Marginal difference in RMST up to 5 years +#' marginal_effects(ndmm_fit, type = "rmst", mtype = "difference", times = 5) +#' +#' # Marginal median survival time ratios +#' marginal_effects(ndmm_fit, type = "median", mtype = "ratio") +#' +#' # Marginal hazard ratios, notice that these are time-varying +#' # Here we specify a vector of times to avoid attempting to plot undefined +#' # hazard ratios for some studies at t=0 +#' plot(marginal_effects(ndmm_fit, type = "hazard", mtype = "ratio", +#' times = seq(0.001, 6, length.out = 50))) +#' +#' } marginal_effects <- function(object, ..., mtype = c("difference", "ratio", "link"), From 68d6309127c92285cfb52425fd5a607c256ae71b Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Fri, 16 Feb 2024 15:05:04 +0000 Subject: [PATCH 17/57] Add return value --- R/marginal_effects.R | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/R/marginal_effects.R b/R/marginal_effects.R index c82e97b9..2122b0c4 100644 --- a/R/marginal_effects.R +++ b/R/marginal_effects.R @@ -26,7 +26,9 @@ #' study be returned? Default `FALSE`. #' @param summary Logical, calculate posterior summaries? Default `TRUE`. #' -#' @return +#' @return A [nma_summary] object if `summary = TRUE`, otherwise a list +#' containing a 3D MCMC array of samples and (for regression models) a data +#' frame of study information. #' @export #' #' @examples ## Smoking cessation From ec7a72f75807db39158786e0e79ce49e5bca073a Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Fri, 16 Feb 2024 15:18:39 +0000 Subject: [PATCH 18/57] Document --- NAMESPACE | 1 + man/marginal_effects.Rd | 154 ++++++++++++++++++++++++++++++++++++++++ man/predict.stan_nma.Rd | 12 +++- 3 files changed, 165 insertions(+), 2 deletions(-) create mode 100644 man/marginal_effects.Rd diff --git a/NAMESPACE b/NAMESPACE index 403893a3..abbd0d11 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -83,6 +83,7 @@ export(is_network_connected) export(log_normal) export(log_student_t) export(make_knots) +export(marginal_effects) export(multi) export(nma) export(normal) diff --git a/man/marginal_effects.Rd b/man/marginal_effects.Rd new file mode 100644 index 00000000..06d40ed2 --- /dev/null +++ b/man/marginal_effects.Rd @@ -0,0 +1,154 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/marginal_effects.R +\name{marginal_effects} +\alias{marginal_effects} +\title{Marginal treatment effects} +\usage{ +marginal_effects( + object, + ..., + mtype = c("difference", "ratio", "link"), + all_contrasts = FALSE, + trt_ref = NULL, + probs = c(0.025, 0.25, 0.5, 0.75, 0.975), + predictive_distribution = FALSE, + summary = TRUE +) +} +\arguments{ +\item{object}{A \code{stan_nma} object created by \code{\link[=nma]{nma()}}.} + +\item{...}{Arguments passed to \code{\link[=predict.stan_nma]{predict.stan_nma()}}, for example to specify +the covariate distribution and baseline risk for a target population.} + +\item{mtype}{The type of marginal effect to construct from the average +absolute effects, either \code{"difference"} (the default) for a difference of +absolute effects such as a risk difference, \code{"ratio"} for a ratio of +absolute effects such as a risk ratio, or \code{"link"} for a difference on the +scale of the link function used in fitting the model such as a marginal log +odds ratio.} + +\item{all_contrasts}{Logical, generate estimates for all contrasts (\code{TRUE}), +or just the "basic" contrasts against the network reference treatment +(\code{FALSE})? Default \code{FALSE}.} + +\item{trt_ref}{Reference treatment to construct relative effects against, if +\code{all_contrasts = FALSE}. By default, relative effects will be against the +network reference treatment. Coerced to character string.} + +\item{probs}{Numeric vector of quantiles of interest to present in computed +summary, default \code{c(0.025, 0.25, 0.5, 0.75, 0.975)}} + +\item{predictive_distribution}{Logical, when a random effects model has been +fitted, should the predictive distribution for marginal effects in a new +study be returned? Default \code{FALSE}.} + +\item{summary}{Logical, calculate posterior summaries? Default \code{TRUE}.} +} +\value{ +A \link{nma_summary} object if \code{summary = TRUE}, otherwise a list +containing a 3D MCMC array of samples and (for regression models) a data +frame of study information. +} +\description{ +Generate population-average marginal treatment effects. These are formed from +population-average absolute predictions, so this function is a wrapper around +\code{\link[=predict.stan_nma]{predict.stan_nma()}}. +} +\examples{ +## Smoking cessation +\donttest{ +# Run smoking RE NMA example if not already available +if (!exists("smk_fit_RE")) example("example_smk_re", run.donttest = TRUE) +} +\donttest{ +# Marginal risk difference in each study population in the network +marginal_effects(smk_fit_RE, mtype = "difference") + +# Since there are no covariates in the model, the marginal and conditional +# (log) odds ratios here coincide +marginal_effects(smk_fit_RE, mtype = "link") +relative_effects(smk_fit_RE) + +# Marginal risk differences in a population with 67 observed events out of +# 566 individuals on No Intervention, corresponding to a Beta(67, 566 - 67) +# distribution on the baseline probability of response +(smk_rd_RE <- marginal_effects(smk_fit_RE, + baseline = distr(qbeta, 67, 566 - 67), + baseline_type = "response", + mtype = "difference")) +plot(smk_rd_RE) +} + +## Plaque psoriasis ML-NMR +\donttest{ +# Run plaque psoriasis ML-NMR example if not already available +if (!exists("pso_fit")) example("example_pso_mlnmr", run.donttest = TRUE) +} +\donttest{ +# Population-average marginal probit differences in each study in the network +(pso_marg <- marginal_effects(pso_fit, mtype = "link")) +plot(pso_marg, ref_line = c(0, 1)) + +# Population-average marginal probit differences in a new target population, +# with means and SDs or proportions given by +new_agd_int <- data.frame( + bsa_mean = 0.6, + bsa_sd = 0.3, + prevsys = 0.1, + psa = 0.2, + weight_mean = 10, + weight_sd = 1, + durnpso_mean = 3, + durnpso_sd = 1 +) + +# We need to add integration points to this data frame of new data +# We use the weighted mean correlation matrix computed from the IPD studies +new_agd_int <- add_integration(new_agd_int, + durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd), + prevsys = distr(qbern, prob = prevsys), + bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd), + weight = distr(qgamma, mean = weight_mean, sd = weight_sd), + psa = distr(qbern, prob = psa), + cor = pso_net$int_cor, + n_int = 64) + +# Population-average marginal probit differences of achieving PASI 75 in this +# target population, given a Normal(-1.75, 0.08^2) distribution on the +# baseline probit-probability of response on Placebo (at the reference levels +# of the covariates), are given by +(pso_marg_new <- marginal_effects(pso_fit, + mtype = "link", + newdata = new_agd_int, + baseline = distr(qnorm, -1.75, 0.08))) +plot(pso_marg_new) +} + +## Progression free survival with newly-diagnosed multiple myeloma +\donttest{ +# Run newly-diagnosed multiple myeloma example if not already available +if (!exists("ndmm_fit")) example("example_ndmm", run.donttest = TRUE) +} +\donttest{ +# We can produce a range of marginal effects from models with survival +# outcomes, specified with the mtype and type arguments. For example: + +# Marginal survival probability difference at 5 years, all contrasts +marginal_effects(ndmm_fit, type = "survival", mtype = "difference", + times = 5, all_contrasts = TRUE) + +# Marginal difference in RMST up to 5 years +marginal_effects(ndmm_fit, type = "rmst", mtype = "difference", times = 5) + +# Marginal median survival time ratios +marginal_effects(ndmm_fit, type = "median", mtype = "ratio") + +# Marginal hazard ratios, notice that these are time-varying +# Here we specify a vector of times to avoid attempting to plot undefined +# hazard ratios for some studies at t=0 +plot(marginal_effects(ndmm_fit, type = "hazard", mtype = "ratio", + times = seq(0.001, 6, length.out = 50))) + +} +} diff --git a/man/predict.stan_nma.Rd b/man/predict.stan_nma.Rd index c3690643..aba401fa 100644 --- a/man/predict.stan_nma.Rd +++ b/man/predict.stan_nma.Rd @@ -18,7 +18,8 @@ baseline_level = c("individual", "aggregate"), probs = c(0.025, 0.25, 0.5, 0.75, 0.975), predictive_distribution = FALSE, - summary = TRUE + summary = TRUE, + progress = FALSE ) \method{predict}{stan_nma_surv}( @@ -36,7 +37,8 @@ times_seq = NULL, probs = c(0.025, 0.25, 0.5, 0.75, 0.975), predictive_distribution = FALSE, - summary = TRUE + summary = TRUE, + progress = interactive() ) } \arguments{ @@ -143,6 +145,12 @@ study be returned? Default \code{FALSE}.} \item{summary}{Logical, calculate posterior summaries? Default \code{TRUE}.} +\item{progress}{Logical, display progress for potentially long-running +calculations? Population-average predictions from ML-NMR models are +computationally intensive, especially for survival outcomes. Currently the +default is to display progress only when running interactively and +producing predictions for a survival ML-NMR model.} + \item{times}{A numeric vector of times to evaluate predictions at. Alternatively, if \code{newdata} is specified, \code{times} can be the name of a column in \code{newdata} which contains the times. If \code{NULL} (the default) then From 01250bf6c293e5c4f78d42e8cf7fc1e73a4c49f0 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Fri, 16 Feb 2024 15:19:52 +0000 Subject: [PATCH 19/57] Include in pkgdown index --- _pkgdown.yml | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/_pkgdown.yml b/_pkgdown.yml index 0d50706f..38bcea20 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -118,10 +118,12 @@ reference: - title: Posterior summaries and working with fitted models desc: > - Producing and plotting relative effects, absolute predictions, posterior ranks - and rank probabilities. Converting to MCMC arrays and matrices. + Producing and plotting relative effects, absolute predictions, marginal + effects, posterior ranks and rank probabilities. Converting to MCMC arrays + and matrices. contents: - relative_effects + - marginal_effects - predict.stan_nma - posterior_ranks - posterior_rank_probs From 7940694a700de2cef26c4aaf164179eb5e7deeb6 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Fri, 16 Feb 2024 18:54:28 +0000 Subject: [PATCH 20/57] Fix all_contrasts labelling --- R/marginal_effects.R | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/R/marginal_effects.R b/R/marginal_effects.R index 2122b0c4..ca202a92 100644 --- a/R/marginal_effects.R +++ b/R/marginal_effects.R @@ -194,8 +194,7 @@ marginal_effects <- function(object, out_meta <- dplyr::distinct(out_meta, dplyr::pick(dplyr::all_of(vars))) %>% dplyr::mutate(.trtb = list(trtb), .trta = list(trta)) %>% tidyr::unnest(cols = c(.data$.trtb, .data$.trta)) %>% - dplyr::select(dplyr::any_of(".study"), ".trtb", ".trta", dplyr::everything()) %>% - dplyr::arrange(dplyr::pick(dplyr::any_of(".study")), ".trta", ".trtb", dplyr::pick(dplyr::everything())) + dplyr::select(dplyr::any_of(".study"), ".trtb", ".trta", dplyr::everything()) } pred_meta$id <- 1:nrow(pred_meta) From 29cdfeb1fdb67ef54f4bf059cf7e13a8f82221f2 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 28 Feb 2024 11:01:43 +0000 Subject: [PATCH 21/57] Rename predict trt_ref argument to baseline_trt --- R/predict.R | 80 +++++++++++++++++++++++++++++------------------------ 1 file changed, 44 insertions(+), 36 deletions(-) diff --git a/R/predict.R b/R/predict.R index 536fa70d..c97c4ae2 100644 --- a/R/predict.R +++ b/R/predict.R @@ -39,8 +39,8 @@ #' always `"link"`, and `baseline_level` is `"individual"` for IPD NMA or #' ML-NMR, and `"aggregate"` for AgD NMA). #' -#' Use the `trt_ref` argument to specify which treatment this distribution -#' applies to. +#' Use the `baseline_trt` argument to specify which treatment this +#' distribution applies to. #' @param newdata Only required if a regression model is fitted and `baseline` #' is specified. A data frame of covariate details, for which to produce #' predictions. Column names must match variables in the regression model. @@ -59,10 +59,6 @@ #' specified: if `newdata` contains integration points produced by #' [add_integration()], studies will be labelled sequentially by row; #' otherwise data will be assumed to come from a single study. -#' @param trt_ref Treatment to which the `baseline` response distribution -#' refers, if `baseline` is specified. By default, the baseline response -#' distribution will refer to the network reference treatment. Coerced to -#' character string. #' @param type Whether to produce predictions on the `"link"` scale (the #' default, e.g. log odds) or `"response"` scale (e.g. probabilities). #' @@ -84,6 +80,10 @@ #' specified, predictions are produced for all IPD studies in the network if #' `level` is `"individual"` or `"aggregate"`, and for all arm-based AgD #' studies in the network if `level` is `"aggregate"`. +#' @param baseline_trt Treatment to which the `baseline` response distribution +#' refers, if `baseline` is specified. By default, the baseline response +#' distribution will refer to the network reference treatment. Coerced to +#' character string. #' @param baseline_type When a `baseline` distribution is given, specifies #' whether this corresponds to the `"link"` scale (the default, e.g. log odds) #' or `"response"` scale (e.g. probabilities). For survival models, `baseline` @@ -109,6 +109,7 @@ #' computationally intensive, especially for survival outcomes. Currently the #' default is to display progress only when running interactively and #' producing predictions for a survival ML-NMR model. +#' @param trt_ref Deprecated, use `baseline_trt` instead. #' #' @details #' # Aggregate-level predictions from IPD NMA and ML-NMR models @@ -279,15 +280,17 @@ #' predict(ndmm_fit, type = "rmst", times = 5) # 5-year RMST #' } predict.stan_nma <- function(object, ..., - baseline = NULL, newdata = NULL, study = NULL, trt_ref = NULL, + baseline = NULL, newdata = NULL, study = NULL, type = c("link", "response"), level = c("aggregate", "individual"), + baseline_trt = NULL, baseline_type = c("link", "response"), baseline_level = c("individual", "aggregate"), probs = c(0.025, 0.25, 0.5, 0.75, 0.975), predictive_distribution = FALSE, summary = TRUE, - progress = FALSE) { + progress = FALSE, + trt_ref = NULL) { # Checks if (!inherits(object, "stan_nma")) abort("Expecting a `stan_nma` object, as returned by nma().") @@ -315,23 +318,26 @@ predict.stan_nma <- function(object, ..., # Get network reference treatment nrt <- levels(object$network$treatments)[1] - if (!is.null(trt_ref)) { + # Deprecated trt_ref in favour of baseline_trt + if (is.null(baseline_trt) && !is.null(trt_ref)) baseline_trt <- trt_ref + + if (!is.null(baseline_trt)) { if (is.null(baseline)) { - # warn("Ignoring `trt_ref` since `baseline` is not given.") - trt_ref <- nrt + # warn("Ignoring `baseline_trt` since `baseline` is not given.") + baseline_trt <- nrt } else { - if (length(trt_ref) > 1) abort("`trt_ref` must be length 1.") - trt_ref <- as.character(trt_ref) + if (length(baseline_trt) > 1) abort("`baseline_trt` must be length 1.") + baseline_trt <- as.character(baseline_trt) lvls_trt <- levels(object$network$treatments) - if (! trt_ref %in% lvls_trt) - abort(sprintf("`trt_ref` does not match a treatment in the network.\nSuitable values are: %s", + if (! baseline_trt %in% lvls_trt) + abort(sprintf("`baseline_trt` does not match a treatment in the network.\nSuitable values are: %s", ifelse(length(lvls_trt) <= 5, paste0(lvls_trt, collapse = ", "), paste0(paste0(lvls_trt[1:5], collapse = ", "), ", ...")))) } } else { - # Set trt_ref to network reference treatment if unset - trt_ref <- nrt + # Set baseline_trt to network reference treatment if unset + baseline_trt <- nrt } # Define auxiliary parameters for survival models @@ -558,9 +564,9 @@ predict.stan_nma <- function(object, ..., mu <- link_fun(mu, link = object$link) } - # Convert to samples on network ref trt if trt_ref given - if (trt_ref != nrt) { - mu <- mu - d[ , , paste0("d[", trt_ref, "]"), drop = FALSE] + # Convert to samples on network ref trt if baseline_trt given + if (baseline_trt != nrt) { + mu <- mu - d[ , , paste0("d[", baseline_trt, "]"), drop = FALSE] } # Combine mu and d @@ -1354,9 +1360,9 @@ predict.stan_nma <- function(object, ..., mu[ , , baseline_type == "response"] <- link_fun(mu[ , , baseline_type == "response"], link = object$link) } - # Convert to samples on network ref trt if trt_ref given - if (trt_ref != nrt) { - mu <- sweep(mu, 1:2, post_temp[ , , paste0("d[", trt_ref, "]"), drop = FALSE], FUN = "-") + # Convert to samples on network ref trt if baseline_trt given + if (baseline_trt != nrt) { + mu <- sweep(mu, 1:2, post_temp[ , , paste0("d[", baseline_trt, "]"), drop = FALSE], FUN = "-") } } else { # ML-NMR or IPD NMR if (any(baseline_level == "individual")) { @@ -1366,9 +1372,9 @@ predict.stan_nma <- function(object, ..., mu[ , , baseline_level == "individual" & baseline_type == "response"] <- link_fun(mu[ , , baseline_level == "individual" & baseline_type == "response"], link = object$link) } - # Convert to samples on network ref trt if trt_ref given - if (trt_ref != nrt) { - mu[ , , baseline_level == "individual" & baseline_level == "individual"] <- sweep(mu[ , , baseline_level == "individual" & baseline_level == "individual", drop = FALSE], 1:2, post_temp[ , , paste0("d[", trt_ref, "]"), drop = FALSE], FUN = "-") + # Convert to samples on network ref trt if baseline_trt given + if (baseline_trt != nrt) { + mu[ , , baseline_level == "individual" & baseline_level == "individual"] <- sweep(mu[ , , baseline_level == "individual" & baseline_level == "individual", drop = FALSE], 1:2, post_temp[ , , paste0("d[", baseline_trt, "]"), drop = FALSE], FUN = "-") } } @@ -1383,21 +1389,21 @@ predict.stan_nma <- function(object, ..., mu0[ , , baseline_level == "aggregate" & baseline_type == "link"] <- inverse_link(mu[ , , baseline_level == "aggregate" & baseline_type == "link"], link = object$link) } - preddat_trt_ref <- dplyr::filter(preddat, .data$.trt == trt_ref) + preddat_trt_ref <- dplyr::filter(preddat, .data$.trt == baseline_trt) - # Get posterior samples of betas and d[trt_ref] + # Get posterior samples of betas and d[baseline_trt] post_beta <- as.array(object, pars = "beta") - if (trt_ref == nrt) { + if (baseline_trt == nrt) { post_d <- 0 } else { - post_d <- as.array(object, pars = paste0("d[", trt_ref, "]")) + post_d <- as.array(object, pars = paste0("d[", baseline_trt, "]")) } - # Get design matrix for regression for trt_ref - X_trt_ref <- X_all[preddat$.trt == trt_ref, , drop = FALSE] + # Get design matrix for regression for baseline_trt + X_trt_ref <- X_all[preddat$.trt == baseline_trt, , drop = FALSE] X_beta_trt_ref <- X_trt_ref[ , !grepl("^(\\.study|\\.trt|\\.contr)[^:]+$", colnames(X_trt_ref)), drop = FALSE] - if (!is.null(offset_all)) offset_trt_ref <- offset_all[preddat$.trt == trt_ref] + if (!is.null(offset_all)) offset_trt_ref <- offset_all[preddat$.trt == baseline_trt] range_mu <- range(as.array(object, pars = "mu")) @@ -1425,7 +1431,7 @@ predict.stan_nma <- function(object, ..., rtsolve <- uniroot(mu_solve, interval = range_mu, extendInt = "yes", ..., mu0 = mu0[i_iter, i_chain, s, drop = TRUE], post_beta = post_beta[i_iter, i_chain, , drop = TRUE], - post_d = if (trt_ref == nrt) 0 else post_d[i_iter, i_chain, , drop = TRUE], + post_d = if (baseline_trt == nrt) 0 else post_d[i_iter, i_chain, , drop = TRUE], X_beta = s_X_beta, offset = if (!is.null(offset_all)) s_offset else 0, link = object$link) @@ -2021,9 +2027,10 @@ predict.stan_nma <- function(object, ..., #' @rdname predict.stan_nma predict.stan_nma_surv <- function(object, times = NULL, ..., + baseline_trt = NULL, baseline = NULL, aux = NULL, - newdata = NULL, study = NULL, trt_ref = NULL, + newdata = NULL, study = NULL, type = c("survival", "hazard", "cumhaz", "mean", "median", "quantile", "rmst", "link"), quantiles = c(0.25, 0.5, 0.75), level = c("aggregate", "individual"), @@ -2031,7 +2038,8 @@ predict.stan_nma_surv <- function(object, times = NULL, probs = c(0.025, 0.25, 0.5, 0.75, 0.975), predictive_distribution = FALSE, summary = TRUE, - progress = interactive()) { + progress = interactive(), + trt_ref = NULL) { type <- rlang::arg_match(type) times <- rlang::enquo(times) From fd21c0fe796870fb06077c5563e27a9d358b4ae1 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 28 Feb 2024 13:04:43 +0000 Subject: [PATCH 22/57] Ignore type and level arguments with friendly warning, rather than error --- R/marginal_effects.R | 29 +++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/R/marginal_effects.R b/R/marginal_effects.R index ca202a92..7be1ac1f 100644 --- a/R/marginal_effects.R +++ b/R/marginal_effects.R @@ -165,18 +165,27 @@ marginal_effects <- function(object, abort(glue::glue("Cannot produce marginal effects under inconsistency '{object$consistency}' model.")) # Call predict to get absolute predictions + if ("level" %in% ...names()) warn('Ignoring `level` argument, this is always "aggregate" for marginal effects.') + if (object$likelihood %in% valid_lhood$survival) { - pred <- stats::predict(object, - level = "aggregate", - summary = FALSE, - predictive_distribution = predictive_distribution, - ...) + pred <- rlang::eval_tidy(rlang::call2( + stats::predict, + !!! rlang::dots_list(object, + level = "aggregate", + summary = FALSE, + predictive_distribution = predictive_distribution, + !!! rlang::enquos(...), + .homonyms = "first"))) } else { - pred <- stats::predict(object, - type = "response", level = "aggregate", - summary = FALSE, - predictive_distribution = predictive_distribution, - ...) + if ("type" %in% ...names()) warn('Ignoring `type` argument, this is always "response" for marginal effects (except for survival models).') + pred <- rlang::eval_tidy(rlang::call2( + stats::predict, + !!! rlang::dots_list(object, + type = "response", level = "aggregate", + summary = FALSE, + predictive_distribution = predictive_distribution, + !!! rlang::enquos(...), + .homonyms = "first"))) } # Set up output metadata From 146b4284b72958b578261ff2a2a4143b31f97683 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 28 Feb 2024 13:05:06 +0000 Subject: [PATCH 23/57] More detail for ... args --- R/marginal_effects.R | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/R/marginal_effects.R b/R/marginal_effects.R index 7be1ac1f..9021da0e 100644 --- a/R/marginal_effects.R +++ b/R/marginal_effects.R @@ -6,7 +6,13 @@ #' #' @param object A `stan_nma` object created by [nma()]. #' @param ... Arguments passed to [predict.stan_nma()], for example to specify -#' the covariate distribution and baseline risk for a target population. +#' the covariate distribution and baseline risk for a target population, e.g. +#' `newdata`, `baseline`, and related arguments. For survival outcomes, `type` +#' can also be specified to determine the quantity from which to form a +#' marginal effect. For example, `type = "hazard"` with `mtype = "ratio"` +#' produces marginal hazard ratios, `type = "median"` with `mtype = +#' "difference"` produces marginal median survival time differences, and so +#' on. #' @param mtype The type of marginal effect to construct from the average #' absolute effects, either `"difference"` (the default) for a difference of #' absolute effects such as a risk difference, `"ratio"` for a ratio of From 47c5961198c2e7728d1facd2490bd37af8e07798 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Thu, 7 Mar 2024 11:54:02 +0000 Subject: [PATCH 24/57] Document --- man/marginal_effects.Rd | 7 ++++++- man/predict.stan_nma.Rd | 26 +++++++++++++++----------- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/man/marginal_effects.Rd b/man/marginal_effects.Rd index 06d40ed2..649dad58 100644 --- a/man/marginal_effects.Rd +++ b/man/marginal_effects.Rd @@ -19,7 +19,12 @@ marginal_effects( \item{object}{A \code{stan_nma} object created by \code{\link[=nma]{nma()}}.} \item{...}{Arguments passed to \code{\link[=predict.stan_nma]{predict.stan_nma()}}, for example to specify -the covariate distribution and baseline risk for a target population.} +the covariate distribution and baseline risk for a target population, e.g. +\code{newdata}, \code{baseline}, and related arguments. For survival outcomes, \code{type} +can also be specified to determine the quantity from which to form a +marginal effect. For example, \code{type = "hazard"} with \code{mtype = "ratio"} +produces marginal hazard ratios, \code{type = "median"} with \code{mtype = "difference"} produces marginal median survival time differences, and so +on.} \item{mtype}{The type of marginal effect to construct from the average absolute effects, either \code{"difference"} (the default) for a difference of diff --git a/man/predict.stan_nma.Rd b/man/predict.stan_nma.Rd index aba401fa..5a1e16b6 100644 --- a/man/predict.stan_nma.Rd +++ b/man/predict.stan_nma.Rd @@ -11,26 +11,27 @@ baseline = NULL, newdata = NULL, study = NULL, - trt_ref = NULL, type = c("link", "response"), level = c("aggregate", "individual"), + baseline_trt = NULL, baseline_type = c("link", "response"), baseline_level = c("individual", "aggregate"), probs = c(0.025, 0.25, 0.5, 0.75, 0.975), predictive_distribution = FALSE, summary = TRUE, - progress = FALSE + progress = FALSE, + trt_ref = NULL ) \method{predict}{stan_nma_surv}( object, times = NULL, ..., + baseline_trt = NULL, baseline = NULL, aux = NULL, newdata = NULL, study = NULL, - trt_ref = NULL, type = c("survival", "hazard", "cumhaz", "mean", "median", "quantile", "rmst", "link"), quantiles = c(0.25, 0.5, 0.75), level = c("aggregate", "individual"), @@ -38,7 +39,8 @@ probs = c(0.025, 0.25, 0.5, 0.75, 0.975), predictive_distribution = FALSE, summary = TRUE, - progress = interactive() + progress = interactive(), + trt_ref = NULL ) } \arguments{ @@ -70,8 +72,8 @@ intercept parameters in the linear predictor (i.e. \code{baseline_type} is always \code{"link"}, and \code{baseline_level} is \code{"individual"} for IPD NMA or ML-NMR, and \code{"aggregate"} for AgD NMA). -Use the \code{trt_ref} argument to specify which treatment this distribution -applies to.} +Use the \code{baseline_trt} argument to specify which treatment this +distribution applies to.} \item{newdata}{Only required if a regression model is fitted and \code{baseline} is specified. A data frame of covariate details, for which to produce @@ -93,11 +95,6 @@ specified: if \code{newdata} contains integration points produced by \code{\link[=add_integration]{add_integration()}}, studies will be labelled sequentially by row; otherwise data will be assumed to come from a single study.} -\item{trt_ref}{Treatment to which the \code{baseline} response distribution -refers, if \code{baseline} is specified. By default, the baseline response -distribution will refer to the network reference treatment. Coerced to -character string.} - \item{type}{Whether to produce predictions on the \code{"link"} scale (the default, e.g. log odds) or \code{"response"} scale (e.g. probabilities). @@ -120,6 +117,11 @@ specified, predictions are produced for all IPD studies in the network if \code{level} is \code{"individual"} or \code{"aggregate"}, and for all arm-based AgD studies in the network if \code{level} is \code{"aggregate"}.} +\item{baseline_trt}{Treatment to which the \code{baseline} response distribution +refers, if \code{baseline} is specified. By default, the baseline response +distribution will refer to the network reference treatment. Coerced to +character string.} + \item{baseline_type}{When a \code{baseline} distribution is given, specifies whether this corresponds to the \code{"link"} scale (the default, e.g. log odds) or \code{"response"} scale (e.g. probabilities). For survival models, \code{baseline} @@ -151,6 +153,8 @@ computationally intensive, especially for survival outcomes. Currently the default is to display progress only when running interactively and producing predictions for a survival ML-NMR model.} +\item{trt_ref}{Deprecated, use \code{baseline_trt} instead.} + \item{times}{A numeric vector of times to evaluate predictions at. Alternatively, if \code{newdata} is specified, \code{times} can be the name of a column in \code{newdata} which contains the times. If \code{NULL} (the default) then From 389351132afb616a8ba71bd921f3af28fc52e848 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Thu, 7 Mar 2024 12:37:27 +0000 Subject: [PATCH 25/57] Always use .trtb and .trta columns in output, matching relative_effects --- R/marginal_effects.R | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/R/marginal_effects.R b/R/marginal_effects.R index 9021da0e..4e200d44 100644 --- a/R/marginal_effects.R +++ b/R/marginal_effects.R @@ -201,7 +201,9 @@ marginal_effects <- function(object, vars <- intersect(c(".study", ".time", ".quantile", ".category"), colnames(pred_meta)) if (!all_contrasts) { - out_meta <- dplyr::filter(out_meta, .data$.trt != trt_ref) + out_meta <- dplyr::filter(out_meta, .data$.trt != trt_ref) %>% + dplyr::mutate(.trtb = .data$.trt, .trta = factor(trt_ref, levels = levels(object$network$treatments))) %>% + dplyr::select(dplyr::any_of(".study"), ".trtb", ".trta", dplyr::everything(), -".trt") } else { contrs <- utils::combn(object$network$treatments, 2) trtb <- contrs[2, ] @@ -216,7 +218,7 @@ marginal_effects <- function(object, out_meta$id <- 1:nrow(out_meta) if (".time" %in% vars) { - time_id <- dplyr::group_by(out_meta, .data$.study, dplyr::pick(dplyr::any_of(c(".trt", ".trta", ".trtb")))) %>% + time_id <- dplyr::group_by(out_meta, .data$.study, .data$.trta, .data$.trtb) %>% dplyr::mutate(time_id = 1:dplyr::n()) %>% dplyr::pull("time_id") } @@ -224,7 +226,7 @@ marginal_effects <- function(object, # Output parameter names pnames <- paste0("marg[", if (rlang::has_name(out_meta, ".study")) paste0(out_meta$.study, ": ") else character(), - if (all_contrasts) out_meta$.trtb else out_meta$.trt, + out_meta$.trtb, if (all_contrasts) paste0(" vs. ", out_meta$.trta) else character(), if (".time" %in% vars) paste0(", ", time_id) else if (".quantile" %in% vars) paste0(", ", out_meta$.quantile) From 9cbf9980478bb91d5ad2d9f53ec53b5c73decdd8 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Thu, 7 Mar 2024 12:37:48 +0000 Subject: [PATCH 26/57] Tests --- tests/testthat/test-marginal_effects.R | 288 +++++++++++++++++++++++++ 1 file changed, 288 insertions(+) create mode 100644 tests/testthat/test-marginal_effects.R diff --git a/tests/testthat/test-marginal_effects.R b/tests/testthat/test-marginal_effects.R new file mode 100644 index 00000000..3e68a9ca --- /dev/null +++ b/tests/testthat/test-marginal_effects.R @@ -0,0 +1,288 @@ + +smk_net <- set_agd_arm(smoking, + study = studyn, + trt = trtc, + r = r, + n = n, + trt_ref = "No intervention") + +# Only test gradients, no sampling +smk_fit_RE <- nma(smk_net, + trt_effects = "random", + prior_intercept = normal(scale = 100), + prior_trt = normal(scale = 100), + prior_het = normal(scale = 5), + test_grad = TRUE) + +test_that("all_contrasts argument", { + m <- "should be TRUE or FALSE" + expect_error(marginal_effects(smk_fit_RE, all_contrasts = "a"), m) + expect_error(marginal_effects(smk_fit_RE, all_contrasts = 1), m) + expect_error(marginal_effects(smk_fit_RE, all_contrasts = list()), m) + expect_error(marginal_effects(smk_fit_RE, all_contrasts = NA), m) + expect_error(marginal_effects(smk_fit_RE, all_contrasts = NULL), m) +}) + +test_that("trt_ref argument", { + m <- "does not match a treatment in the network" + expect_error(marginal_effects(smk_fit_RE, trt_ref = "a"), m) + expect_error(marginal_effects(smk_fit_RE, trt_ref = 1), m) + expect_error(marginal_effects(smk_fit_RE, trt_ref = list("a")), m) + expect_error(marginal_effects(smk_fit_RE, trt_ref = NA), m) + + # expect_warning(marginal_effects(smk_fit_RE, all_contrasts = TRUE, trt_ref = "Self-help"), "Ignoring `trt_ref`") +}) + +test_that("probs argument", { + m <- "numeric vector of probabilities" + expect_error(marginal_effects(smk_fit_RE, probs = "a"), m) + expect_error(marginal_effects(smk_fit_RE, probs = -1), m) + expect_error(marginal_effects(smk_fit_RE, probs = 1.5), m) + expect_error(marginal_effects(smk_fit_RE, probs = Inf), m) + expect_error(marginal_effects(smk_fit_RE, probs = list()), m) + expect_error(marginal_effects(smk_fit_RE, probs = NA), m) + expect_error(marginal_effects(smk_fit_RE, probs = NULL), m) +}) + +test_that("summary argument", { + m <- "should be TRUE or FALSE" + expect_error(marginal_effects(smk_fit_RE, summary = "a"), m) + expect_error(marginal_effects(smk_fit_RE, summary = 1), m) + expect_error(marginal_effects(smk_fit_RE, summary = list()), m) + expect_error(marginal_effects(smk_fit_RE, summary = NA), m) + expect_error(marginal_effects(smk_fit_RE, summary = NULL), m) +}) + +test_that("newdata argument", { + m <- "not a data frame" + expect_error(marginal_effects(smk_fit_RE, newdata = "a"), m) + expect_error(marginal_effects(smk_fit_RE, newdata = 1), m) + expect_error(marginal_effects(smk_fit_RE, newdata = list()), m) + expect_error(marginal_effects(smk_fit_RE, newdata = NA), m) +}) + + +skip_on_cran() # Reduce CRAN check time + +# Only small number of samples to test +smk_fit_RE <- suppressWarnings(nma(smk_net, + trt_effects = "random", + prior_intercept = normal(scale = 100), + prior_trt = normal(scale = 100), + prior_het = normal(scale = 5), + iter = 10)) + +test_that(".trta, .trtb columns are correct", { + re1.1 <- tibble::as_tibble(marginal_effects(smk_fit_RE, mtype = "difference")) + expect_identical(paste0("marg[", re1.1$.study, ": ", re1.1$.trtb, "]"), + re1.1$parameter) + expect_identical(as.character(re1.1$.trta), rep("No intervention", times = nrow(re1.1))) + + re1.2 <- tibble::as_tibble(marginal_effects(smk_fit_RE, mtype = "ratio")) + expect_identical(paste0("marg[", re1.1$.study, ": ", re1.2$.trtb, "]"), + re1.2$parameter) + expect_identical(as.character(re1.2$.trta), rep("No intervention", times = nrow(re1.2))) + + re1.3 <- tibble::as_tibble(marginal_effects(smk_fit_RE, mtype = "link")) + expect_identical(paste0("marg[", re1.1$.study, ": ", re1.3$.trtb, "]"), + re1.3$parameter) + expect_identical(as.character(re1.3$.trta), rep("No intervention", times = nrow(re1.3))) + + + re2.1 <- tibble::as_tibble(marginal_effects(smk_fit_RE, all_contrasts = TRUE, mtype = "difference")) + expect_identical(paste0("marg[", re2.1$.study, ": ", re2.1$.trtb, " vs. ", re2.1$.trta, "]"), + re2.1$parameter) + + re2.2 <- tibble::as_tibble(marginal_effects(smk_fit_RE, all_contrasts = TRUE, mtype = "ratio")) + expect_identical(paste0("marg[", re2.1$.study, ": ", re2.2$.trtb, " vs. ", re2.2$.trta, "]"), + re2.2$parameter) + + re2.3 <- tibble::as_tibble(marginal_effects(smk_fit_RE, all_contrasts = TRUE, mtype = "link")) + expect_identical(paste0("marg[", re2.1$.study, ": ", re2.3$.trtb, " vs. ", re2.3$.trta, "]"), + re2.3$parameter) + + + re3.1 <- tibble::as_tibble(marginal_effects(smk_fit_RE, trt_ref = "Self-help", mtype = "difference")) + expect_identical(paste0("marg[", re1.1$.study, ": ", re3.1$.trtb, "]"), + re3.1$parameter) + expect_identical(as.character(re3.1$.trta), rep("Self-help", times = nrow(re3.1))) + + re3.2 <- tibble::as_tibble(marginal_effects(smk_fit_RE, trt_ref = "Self-help", mtype = "ratio")) + expect_identical(paste0("marg[", re1.1$.study, ": ", re3.2$.trtb, "]"), + re3.2$parameter) + expect_identical(as.character(re3.2$.trta), rep("Self-help", times = nrow(re3.2))) + + re3.3 <- tibble::as_tibble(marginal_effects(smk_fit_RE, trt_ref = "Self-help", mtype = "link")) + expect_identical(paste0("marg[", re3.1$.study, ": ", re3.3$.trtb, "]"), + re3.3$parameter) + expect_identical(as.character(re3.3$.trta), rep("Self-help", times = nrow(re3.3))) +}) + + +test_that(".trta, .trtb columns are correct in new populations", { + re1.1 <- tibble::as_tibble(marginal_effects(smk_fit_RE, baseline = distr(qnorm, 0, 0.1), mtype = "difference")) + expect_identical(paste0("marg[", re1.1$.trtb, "]"), + re1.1$parameter) + expect_identical(as.character(re1.1$.trta), rep("No intervention", times = nrow(re1.1))) + + re1.2 <- tibble::as_tibble(marginal_effects(smk_fit_RE, baseline = distr(qnorm, 0, 0.1), mtype = "ratio")) + expect_identical(paste0("marg[", re1.2$.trtb, "]"), + re1.2$parameter) + expect_identical(as.character(re1.2$.trta), rep("No intervention", times = nrow(re1.2))) + + re1.3 <- tibble::as_tibble(marginal_effects(smk_fit_RE, baseline = distr(qnorm, 0, 0.1), mtype = "link")) + expect_identical(paste0("marg[", re1.3$.trtb, "]"), + re1.3$parameter) + expect_identical(as.character(re1.3$.trta), rep("No intervention", times = nrow(re1.3))) + + + re2.1 <- tibble::as_tibble(marginal_effects(smk_fit_RE, baseline = distr(qnorm, 0, 0.1), all_contrasts = TRUE, mtype = "difference")) + expect_identical(paste0("marg[", re2.1$.trtb, " vs. ", re2.1$.trta, "]"), + re2.1$parameter) + + re2.2 <- tibble::as_tibble(marginal_effects(smk_fit_RE, baseline = distr(qnorm, 0, 0.1), all_contrasts = TRUE, mtype = "ratio")) + expect_identical(paste0("marg[", re2.2$.trtb, " vs. ", re2.2$.trta, "]"), + re2.2$parameter) + + re2.3 <- tibble::as_tibble(marginal_effects(smk_fit_RE, baseline = distr(qnorm, 0, 0.1), all_contrasts = TRUE, mtype = "link")) + expect_identical(paste0("marg[", re2.3$.trtb, " vs. ", re2.3$.trta, "]"), + re2.3$parameter) + + + re3.1 <- tibble::as_tibble(marginal_effects(smk_fit_RE, baseline = distr(qnorm, 0, 0.1), trt_ref = "Self-help", mtype = "difference")) + expect_identical(paste0("marg[", re3.1$.trtb, "]"), + re3.1$parameter) + expect_identical(as.character(re3.1$.trta), rep("Self-help", times = nrow(re3.1))) + + re3.2 <- tibble::as_tibble(marginal_effects(smk_fit_RE, baseline = distr(qnorm, 0, 0.1), trt_ref = "Self-help", mtype = "ratio")) + expect_identical(paste0("marg[", re3.2$.trtb, "]"), + re3.2$parameter) + expect_identical(as.character(re3.2$.trta), rep("Self-help", times = nrow(re3.2))) + + re3.3 <- tibble::as_tibble(marginal_effects(smk_fit_RE, baseline = distr(qnorm, 0, 0.1), trt_ref = "Self-help", mtype = "link")) + expect_identical(paste0("marg[", re3.3$.trtb, "]"), + re3.3$parameter) + expect_identical(as.character(re3.3$.trta), rep("Self-help", times = nrow(re3.3))) +}) + +library(dplyr) + +pso_ipd <- mutate(plaque_psoriasis_ipd, + bsa = bsa / 100, + prevsys = as.numeric(prevsys), + psa = as.numeric(psa), + weight = weight / 10, + durnpso = durnpso / 10, + # Treatment classes + trtclass = case_when(trtn == 1 ~ "Placebo", + trtn %in% c(2, 3, 5, 6) ~ "IL blocker", + trtn == 4 ~ "TNFa blocker"), + # Check complete cases for covariates of interest + complete = complete.cases(durnpso, prevsys, bsa, weight, psa) + ) + +pso_agd <- mutate(plaque_psoriasis_agd, + bsa_mean = bsa_mean / 100, + bsa_sd = bsa_sd / 100, + prevsys = prevsys / 100, + psa = psa / 100, + weight_mean = weight_mean / 10, + weight_sd = weight_sd / 10, + durnpso_mean = durnpso_mean / 10, + durnpso_sd = durnpso_sd / 10, + # Treatment classes + trtclass = case_when(trtn == 1 ~ "Placebo", + trtn %in% c(2, 3, 5, 6) ~ "IL blocker", + trtn == 4 ~ "TNFa blocker") + ) + +pso_net <- combine_network( + set_ipd(filter(pso_ipd, complete), + studyc, trtc, + r = pasi75), + set_agd_arm(pso_agd, + study = studyc, trt = trtc, + r = pasi75_r, n = pasi75_n)) + +pso_net <- add_integration(pso_net, + durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd), + prevsys = distr(qbern, prob = prevsys), + bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd), + weight = distr(qgamma, mean = weight_mean, sd = weight_sd), + psa = distr(qbern, prob = psa), + n_int = 4) + +# Only small number of samples to test +pso_fit <- suppressWarnings(nma(pso_net, + trt_effects = "fixed", + regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt, + prior_intercept = normal(scale = 10), + prior_trt = normal(scale = 10), + prior_reg = normal(scale = 10), + init_r = 0.1, + iter = 10)) + +pso_new <- tibble( + bsa_mean = 0.6, + bsa_sd = 0.3, + prevsys = 0.1, + psa = 0.2, + weight_mean = 10, + weight_sd = 1, + durnpso_mean = 3, + durnpso_sd = 1, + studyc = "NEW" +) + +pso_new <- add_integration(pso_new, + durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd), + prevsys = distr(qbern, prob = prevsys), + bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd), + weight = distr(qgamma, mean = weight_mean, sd = weight_sd), + psa = distr(qbern, prob = psa), + cor = pso_net$int_cor, + n_int = 4) + + +test_that("newdata validation", { + expect_error(marginal_effects(pso_fit, newdata = dplyr::select(pso_new, -".int_durnpso"), baseline = distr(qnorm, 0, 0.1)), + 'Regression variable "durnpso" not found in `newdata`') + expect_error(marginal_effects(pso_fit, newdata = dplyr::select(pso_new, -".int_durnpso", -".int_weight"), baseline = distr(qnorm, 0, 0.1)), + 'Regression variables "durnpso" and "weight" not found in `newdata`') +}) + +test_that("baseline required when newdata specified", { + m <- "Specify both `newdata` and `baseline`, or neither" + expect_error(marginal_effects(pso_fit, newdata = pso_new), m) + expect_error(marginal_effects(pso_fit, baseline = distr(qnorm, 0, 0.1)), m) +}) + +test_that(".study, .trta, .trtb columns are correct", { + re1 <- tibble::as_tibble(marginal_effects(pso_fit)) + expect_identical(paste0("marg[", re1$.study, ": ", re1$.trtb, "]"), + re1$parameter) + expect_identical(as.character(re1$.trta), rep(levels(pso_net$treatments)[1], times = nrow(re1))) + + re2 <- tibble::as_tibble(marginal_effects(pso_fit, all_contrasts = TRUE)) + expect_identical(paste0("marg[", re2$.study, ": ", re2$.trtb, " vs. ", re2$.trta, "]"), + re2$parameter) + + re3 <- tibble::as_tibble(marginal_effects(pso_fit, trt_ref = "ETN")) + expect_identical(paste0("marg[", re3$.study, ": ", re3$.trtb, "]"), + re3$parameter) + expect_identical(as.character(re3$.trta), rep("ETN", times = nrow(re3))) + + re4 <- tibble::as_tibble(marginal_effects(pso_fit, baseline = distr(qnorm, 0, 0.1), newdata = pso_new, study = studyc)) + expect_identical(paste0("marg[", re4$.study, ": ", re4$.trtb, "]"), + re4$parameter) + expect_identical(as.character(re4$.trta), rep(levels(pso_net$treatments)[1], times = nrow(re4))) + + re5 <- tibble::as_tibble(marginal_effects(pso_fit, all_contrasts = TRUE, baseline = distr(qnorm, 0, 0.1), newdata = pso_new, study = studyc)) + expect_identical(paste0("marg[", re5$.study, ": ", re5$.trtb, " vs. ", re5$.trta, "]"), + re5$parameter) + + re6 <- tibble::as_tibble(marginal_effects(pso_fit, trt_ref = "ETN", baseline = distr(qnorm, 0, 0.1), newdata = pso_new, study = studyc)) + expect_identical(paste0("marg[", re6$.study, ": ", re6$.trtb, "]"), + re6$parameter) + expect_identical(as.character(re6$.trta), rep("ETN", times = nrow(re6))) +}) From a8501b1259fe4edc0f7397aa9b66f116cd986da6 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Thu, 7 Mar 2024 12:53:26 +0000 Subject: [PATCH 27/57] Ordered categorical tests --- tests/testthat/test-marginal_effects.R | 126 ++++++++++++++++++++++++- 1 file changed, 123 insertions(+), 3 deletions(-) diff --git a/tests/testthat/test-marginal_effects.R b/tests/testthat/test-marginal_effects.R index 3e68a9ca..bace3565 100644 --- a/tests/testthat/test-marginal_effects.R +++ b/tests/testthat/test-marginal_effects.R @@ -64,6 +64,9 @@ test_that("newdata argument", { skip_on_cran() # Reduce CRAN check time + +# Binary AgD only -------------------------------------------------------------- + # Only small number of samples to test smk_fit_RE <- suppressWarnings(nma(smk_net, trt_effects = "random", @@ -72,7 +75,7 @@ smk_fit_RE <- suppressWarnings(nma(smk_net, prior_het = normal(scale = 5), iter = 10)) -test_that(".trta, .trtb columns are correct", { +test_that("Binary AgD: .trta, .trtb columns are correct", { re1.1 <- tibble::as_tibble(marginal_effects(smk_fit_RE, mtype = "difference")) expect_identical(paste0("marg[", re1.1$.study, ": ", re1.1$.trtb, "]"), re1.1$parameter) @@ -119,7 +122,7 @@ test_that(".trta, .trtb columns are correct", { }) -test_that(".trta, .trtb columns are correct in new populations", { +test_that("Binary AgD: .trta, .trtb columns are correct in new populations", { re1.1 <- tibble::as_tibble(marginal_effects(smk_fit_RE, baseline = distr(qnorm, 0, 0.1), mtype = "difference")) expect_identical(paste0("marg[", re1.1$.trtb, "]"), re1.1$parameter) @@ -165,6 +168,9 @@ test_that(".trta, .trtb columns are correct in new populations", { expect_identical(as.character(re3.3$.trta), rep("Self-help", times = nrow(re3.3))) }) + +# Binary with IPD regression --------------------------------------------------- + library(dplyr) pso_ipd <- mutate(plaque_psoriasis_ipd, @@ -257,7 +263,7 @@ test_that("baseline required when newdata specified", { expect_error(marginal_effects(pso_fit, baseline = distr(qnorm, 0, 0.1)), m) }) -test_that(".study, .trta, .trtb columns are correct", { +test_that("Binary IPD: .study, .trta, .trtb columns are correct", { re1 <- tibble::as_tibble(marginal_effects(pso_fit)) expect_identical(paste0("marg[", re1$.study, ": ", re1$.trtb, "]"), re1$parameter) @@ -286,3 +292,117 @@ test_that(".study, .trta, .trtb columns are correct", { re6$parameter) expect_identical(as.character(re6$.trta), rep("ETN", times = nrow(re6))) }) + + +# Ordered categorical ----------------------------------------------------- + +ord_net <- set_agd_arm(plaque_psoriasis_agd, + studyc, trtc, + r = multi(r0 = pasi75_n, + PASI75 = pasi75_r, + PASI90 = pasi90_r, + PASI100 = pasi100_r, + type = "ordered", inclusive = TRUE)) + + +# Only small number of samples to test +ord_fit <- suppressWarnings(nma(ord_net, + trt_effects = "random", + prior_intercept = normal(scale = 100), + prior_trt = normal(scale = 100), + prior_het = normal(scale = 5), + iter = 10)) + +test_that("Ordered categorical: .study, .trta, .trtb, .category columns are correct", { + re1.1 <- tibble::as_tibble(marginal_effects(ord_fit, mtype = "difference")) + expect_identical(paste0("marg[", re1.1$.study, ": ", re1.1$.trtb, ", ", re1.1$.category, "]"), + re1.1$parameter) + expect_identical(as.character(re1.1$.trta), rep("SEC_300", times = nrow(re1.1))) + + re1.2 <- tibble::as_tibble(marginal_effects(ord_fit, mtype = "ratio")) + expect_identical(paste0("marg[", re1.1$.study, ": ", re1.2$.trtb, ", ", re1.2$.category, "]"), + re1.2$parameter) + expect_identical(as.character(re1.2$.trta), rep("SEC_300", times = nrow(re1.2))) + + re1.3 <- tibble::as_tibble(marginal_effects(ord_fit, mtype = "link")) + expect_identical(paste0("marg[", re1.1$.study, ": ", re1.3$.trtb, ", ", re1.3$.category, "]"), + re1.3$parameter) + expect_identical(as.character(re1.3$.trta), rep("SEC_300", times = nrow(re1.3))) + + + re2.1 <- tibble::as_tibble(marginal_effects(ord_fit, all_contrasts = TRUE, mtype = "difference")) + expect_identical(paste0("marg[", re2.1$.study, ": ", re2.1$.trtb, " vs. ", re2.1$.trta, ", ", re2.1$.category, "]"), + re2.1$parameter) + + re2.2 <- tibble::as_tibble(marginal_effects(ord_fit, all_contrasts = TRUE, mtype = "ratio")) + expect_identical(paste0("marg[", re2.1$.study, ": ", re2.2$.trtb, " vs. ", re2.2$.trta, ", ", re2.2$.category, "]"), + re2.2$parameter) + + re2.3 <- tibble::as_tibble(marginal_effects(ord_fit, all_contrasts = TRUE, mtype = "link")) + expect_identical(paste0("marg[", re2.1$.study, ": ", re2.3$.trtb, " vs. ", re2.3$.trta, ", ", re2.3$.category, "]"), + re2.3$parameter) + + + re3.1 <- tibble::as_tibble(marginal_effects(ord_fit, trt_ref = "PBO", mtype = "difference")) + expect_identical(paste0("marg[", re1.1$.study, ": ", re3.1$.trtb, ", ", re3.1$.category, "]"), + re3.1$parameter) + expect_identical(as.character(re3.1$.trta), rep("PBO", times = nrow(re3.1))) + + re3.2 <- tibble::as_tibble(marginal_effects(ord_fit, trt_ref = "PBO", mtype = "ratio")) + expect_identical(paste0("marg[", re1.1$.study, ": ", re3.2$.trtb, ", ", re3.2$.category, "]"), + re3.2$parameter) + expect_identical(as.character(re3.2$.trta), rep("PBO", times = nrow(re3.2))) + + re3.3 <- tibble::as_tibble(marginal_effects(ord_fit, trt_ref = "PBO", mtype = "link")) + expect_identical(paste0("marg[", re3.1$.study, ": ", re3.3$.trtb, ", ", re3.3$.category, "]"), + re3.3$parameter) + expect_identical(as.character(re3.3$.trta), rep("PBO", times = nrow(re3.3))) +}) + + +test_that("Ordered categorical: .trta, .trtb, .category columns are correct in new populations", { + re1.1 <- tibble::as_tibble(marginal_effects(ord_fit, baseline = distr(qnorm, 0, 0.1), mtype = "difference")) + expect_identical(paste0("marg[", re1.1$.trtb, ", ", re1.1$.category, "]"), + re1.1$parameter) + expect_identical(as.character(re1.1$.trta), rep("SEC_300", times = nrow(re1.1))) + + re1.2 <- tibble::as_tibble(marginal_effects(ord_fit, baseline = distr(qnorm, 0, 0.1), mtype = "ratio")) + expect_identical(paste0("marg[", re1.2$.trtb, ", ", re1.2$.category, "]"), + re1.2$parameter) + expect_identical(as.character(re1.2$.trta), rep("SEC_300", times = nrow(re1.2))) + + re1.3 <- tibble::as_tibble(marginal_effects(ord_fit, baseline = distr(qnorm, 0, 0.1), mtype = "link")) + expect_identical(paste0("marg[", re1.3$.trtb, ", ", re1.3$.category, "]"), + re1.3$parameter) + expect_identical(as.character(re1.3$.trta), rep("SEC_300", times = nrow(re1.3))) + + + re2.1 <- tibble::as_tibble(marginal_effects(ord_fit, baseline = distr(qnorm, 0, 0.1), all_contrasts = TRUE, mtype = "difference")) + expect_identical(paste0("marg[", re2.1$.trtb, " vs. ", re2.1$.trta, ", ", re2.1$.category, "]"), + re2.1$parameter) + + re2.2 <- tibble::as_tibble(marginal_effects(ord_fit, baseline = distr(qnorm, 0, 0.1), all_contrasts = TRUE, mtype = "ratio")) + expect_identical(paste0("marg[", re2.2$.trtb, " vs. ", re2.2$.trta, ", ", re2.2$.category, "]"), + re2.2$parameter) + + re2.3 <- tibble::as_tibble(marginal_effects(ord_fit, baseline = distr(qnorm, 0, 0.1), all_contrasts = TRUE, mtype = "link")) + expect_identical(paste0("marg[", re2.3$.trtb, " vs. ", re2.3$.trta, ", ", re2.3$.category, "]"), + re2.3$parameter) + + + re3.1 <- tibble::as_tibble(marginal_effects(ord_fit, baseline = distr(qnorm, 0, 0.1), trt_ref = "PBO", mtype = "difference")) + expect_identical(paste0("marg[", re3.1$.trtb, ", ", re3.1$.category, "]"), + re3.1$parameter) + expect_identical(as.character(re3.1$.trta), rep("PBO", times = nrow(re3.1))) + + re3.2 <- tibble::as_tibble(marginal_effects(ord_fit, baseline = distr(qnorm, 0, 0.1), trt_ref = "PBO", mtype = "ratio")) + expect_identical(paste0("marg[", re3.2$.trtb, ", ", re3.2$.category, "]"), + re3.2$parameter) + expect_identical(as.character(re3.2$.trta), rep("PBO", times = nrow(re3.2))) + + re3.3 <- tibble::as_tibble(marginal_effects(ord_fit, baseline = distr(qnorm, 0, 0.1), trt_ref = "PBO", mtype = "link")) + expect_identical(paste0("marg[", re3.3$.trtb, ", ", re3.3$.category, "]"), + re3.3$parameter) + expect_identical(as.character(re3.3$.trta), rep("PBO", times = nrow(re3.3))) +}) + From 1c214d939cc65f1a1ec9ddc07b054d1640ee1d9b Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Thu, 7 Mar 2024 14:44:46 +0000 Subject: [PATCH 28/57] Group by time_id rather than time (avoid errors with duplicate times) --- R/marginal_effects.R | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/R/marginal_effects.R b/R/marginal_effects.R index 4e200d44..04942dbc 100644 --- a/R/marginal_effects.R +++ b/R/marginal_effects.R @@ -196,10 +196,18 @@ marginal_effects <- function(object, # Set up output metadata pred_meta <- pred$summary - out_meta <- pred_meta - vars <- intersect(c(".study", ".time", ".quantile", ".category"), colnames(pred_meta)) + if (".time" %in% vars) { + pred_meta <- dplyr::group_by(pred_meta, .data$.study, .data$.trt) %>% + dplyr::mutate(.time_id = 1:dplyr::n()) %>% + dplyr::ungroup() + + vars <- c(vars, ".time_id") + } + + out_meta <- pred_meta + if (!all_contrasts) { out_meta <- dplyr::filter(out_meta, .data$.trt != trt_ref) %>% dplyr::mutate(.trtb = .data$.trt, .trta = factor(trt_ref, levels = levels(object$network$treatments))) %>% @@ -217,18 +225,12 @@ marginal_effects <- function(object, pred_meta$id <- 1:nrow(pred_meta) out_meta$id <- 1:nrow(out_meta) - if (".time" %in% vars) { - time_id <- dplyr::group_by(out_meta, .data$.study, .data$.trta, .data$.trtb) %>% - dplyr::mutate(time_id = 1:dplyr::n()) %>% - dplyr::pull("time_id") - } - # Output parameter names pnames <- paste0("marg[", if (rlang::has_name(out_meta, ".study")) paste0(out_meta$.study, ": ") else character(), out_meta$.trtb, if (all_contrasts) paste0(" vs. ", out_meta$.trta) else character(), - if (".time" %in% vars) paste0(", ", time_id) + if (".time" %in% vars) paste0(", ", out_meta$.time_id) else if (".quantile" %in% vars) paste0(", ", out_meta$.quantile) else if (".category" %in% vars) paste0(", ", out_meta$.category) else character() @@ -270,6 +272,7 @@ marginal_effects <- function(object, } out_meta <- dplyr::select(out_meta, -"id") + if (".time_id" %in% vars) out_meta <- dplyr::select(out_meta, -".time_id") # Summarise if (summary) { From b2207e520b1fa5e0b2015451c8f9219ce431824b Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Thu, 7 Mar 2024 15:17:57 +0000 Subject: [PATCH 29/57] Ensure all_contrasts still works with 2 treatments --- R/marginal_effects.R | 4 ++-- tests/testthat/test-marginal_effects.R | 18 ++++++++++++++++++ 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/R/marginal_effects.R b/R/marginal_effects.R index 04942dbc..a79a1ed4 100644 --- a/R/marginal_effects.R +++ b/R/marginal_effects.R @@ -142,7 +142,7 @@ marginal_effects <- function(object, if (!is.null(trt_ref)) { if (all_contrasts) { warn("Ignoring `trt_ref` when all_contrasts = TRUE.") - trt_ref <- NULL + trt_ref <- levels(object$network$treatments)[1] } else { if (length(trt_ref) > 1) abort("`trt_ref` must be length 1.") trt_ref <- as.character(trt_ref) @@ -262,7 +262,7 @@ marginal_effects <- function(object, predi <- dplyr::semi_join(pred_meta, grpi, by = vars)$id outi <- dplyr::semi_join(out_meta, grpi, by = vars)$id - if (all_contrasts) { + if (all_contrasts && nlevels(object$network$treatments) > 2) { out_array[ , , outi] <- aperm(apply(pred$sims[ , , predi], MARGIN = 1:2, FUN = mk_contr), c(2, 3, 1)) } else { predi_ref <- dplyr::semi_join(pred_meta, dplyr::mutate(grpi, .trt = trt_ref), by = c(vars, ".trt"))$id diff --git a/tests/testthat/test-marginal_effects.R b/tests/testthat/test-marginal_effects.R index bace3565..373a72b3 100644 --- a/tests/testthat/test-marginal_effects.R +++ b/tests/testthat/test-marginal_effects.R @@ -168,6 +168,24 @@ test_that("Binary AgD: .trta, .trtb columns are correct in new populations", { expect_identical(as.character(re3.3$.trta), rep("Self-help", times = nrow(re3.3))) }) +blocker_net <- set_agd_arm(blocker, studyn, trtc, r = r, n = n) +blocker_fit <- suppressWarnings(nma(blocker_net, + prior_intercept = normal(scale = 100), + prior_trt = normal(scale = 100), + iter = 10)) + +test_that("all_contrasts works with only two treatments", { + expect_equal( + dplyr::select(as.data.frame(marginal_effects(blocker_fit)), -"parameter"), + dplyr::select(as.data.frame(marginal_effects(blocker_fit, all_contrasts = TRUE)) ,-"parameter"), + check.attributes = FALSE + ) + expect_identical( + unname(as.array(marginal_effects(blocker_fit))), + unname(as.array(marginal_effects(blocker_fit, all_contrasts = TRUE))) + ) +}) + # Binary with IPD regression --------------------------------------------------- From b28df9374d4b9ee15ab62a2a618d54c48bd2e3d2 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Thu, 7 Mar 2024 15:27:57 +0000 Subject: [PATCH 30/57] Survival tests --- tests/testthat/test-marginal_effects.R | 72 ++++++++++++++++++++++++-- 1 file changed, 68 insertions(+), 4 deletions(-) diff --git a/tests/testthat/test-marginal_effects.R b/tests/testthat/test-marginal_effects.R index 373a72b3..7f6a7c41 100644 --- a/tests/testthat/test-marginal_effects.R +++ b/tests/testthat/test-marginal_effects.R @@ -338,12 +338,12 @@ test_that("Ordered categorical: .study, .trta, .trtb, .category columns are corr expect_identical(as.character(re1.1$.trta), rep("SEC_300", times = nrow(re1.1))) re1.2 <- tibble::as_tibble(marginal_effects(ord_fit, mtype = "ratio")) - expect_identical(paste0("marg[", re1.1$.study, ": ", re1.2$.trtb, ", ", re1.2$.category, "]"), + expect_identical(paste0("marg[", re1.2$.study, ": ", re1.2$.trtb, ", ", re1.2$.category, "]"), re1.2$parameter) expect_identical(as.character(re1.2$.trta), rep("SEC_300", times = nrow(re1.2))) re1.3 <- tibble::as_tibble(marginal_effects(ord_fit, mtype = "link")) - expect_identical(paste0("marg[", re1.1$.study, ": ", re1.3$.trtb, ", ", re1.3$.category, "]"), + expect_identical(paste0("marg[", re1.3$.study, ": ", re1.3$.trtb, ", ", re1.3$.category, "]"), re1.3$parameter) expect_identical(as.character(re1.3$.trta), rep("SEC_300", times = nrow(re1.3))) @@ -353,11 +353,11 @@ test_that("Ordered categorical: .study, .trta, .trtb, .category columns are corr re2.1$parameter) re2.2 <- tibble::as_tibble(marginal_effects(ord_fit, all_contrasts = TRUE, mtype = "ratio")) - expect_identical(paste0("marg[", re2.1$.study, ": ", re2.2$.trtb, " vs. ", re2.2$.trta, ", ", re2.2$.category, "]"), + expect_identical(paste0("marg[", re2.2$.study, ": ", re2.2$.trtb, " vs. ", re2.2$.trta, ", ", re2.2$.category, "]"), re2.2$parameter) re2.3 <- tibble::as_tibble(marginal_effects(ord_fit, all_contrasts = TRUE, mtype = "link")) - expect_identical(paste0("marg[", re2.1$.study, ": ", re2.3$.trtb, " vs. ", re2.3$.trta, ", ", re2.3$.category, "]"), + expect_identical(paste0("marg[", re2.3$.study, ": ", re2.3$.trtb, " vs. ", re2.3$.trta, ", ", re2.3$.category, "]"), re2.3$parameter) @@ -424,3 +424,67 @@ test_that("Ordered categorical: .trta, .trtb, .category columns are correct in n expect_identical(as.character(re3.3$.trta), rep("PBO", times = nrow(re3.3))) }) + +# Survival ---------------------------------------------------------------- + +surv_net <- set_ipd(ndmm_agd, study, trt, Surv = Surv(eventtime, status)) + +# Only small number of samples to test +surv_fit <- suppressWarnings(nma(surv_net, + likelihood = "weibull", + prior_intercept = normal(scale = 100), + prior_trt = normal(scale = 100), + prior_aux = normal(scale = 5), + iter = 10)) + +test_that("Survival: .study, .trta, .trtb columns are correct", { + re1.1 <- tibble::as_tibble(marginal_effects(surv_fit, type = "median", mtype = "difference")) + expect_identical(paste0("marg[", re1.1$.study, ": ", re1.1$.trtb, "]"), + re1.1$parameter) + expect_identical(as.character(re1.1$.trta), rep("Pbo", times = nrow(re1.1))) + + re1.2 <- tibble::as_tibble(marginal_effects(surv_fit, type = "hazard", mtype = "ratio")) %>% + group_by(.study, .trtb, .trta) %>% + mutate(.id = 1:n()) + expect_identical(paste0("marg[", re1.2$.study, ": ", re1.2$.trtb, ", ", re1.2$.id, "]"), + re1.2$parameter) + expect_identical(as.character(re1.2$.trta), rep("Pbo", times = nrow(re1.2))) + + re1.3 <- tibble::as_tibble(marginal_effects(surv_fit, type = "quantile", mtype = "difference")) + expect_identical(paste0("marg[", re1.3$.study, ": ", re1.3$.trtb, ", ", re1.3$.quantile, "]"), + re1.3$parameter) + expect_identical(as.character(re1.3$.trta), rep("Pbo", times = nrow(re1.3))) + + + re2.1 <- tibble::as_tibble(marginal_effects(surv_fit, all_contrasts = TRUE, type = "median", mtype = "difference")) + expect_identical(paste0("marg[", re2.1$.study, ": ", re2.1$.trtb, " vs. ", re2.1$.trta, "]"), + re2.1$parameter) + + re2.2 <- tibble::as_tibble(marginal_effects(surv_fit, all_contrasts = TRUE, type = "hazard", mtype = "ratio")) %>% + group_by(.study, .trtb, .trta) %>% + mutate(.id = 1:n()) + expect_identical(paste0("marg[", re2.2$.study, ": ", re2.2$.trtb, " vs. ", re2.2$.trta, ", ", re2.2$.id, "]"), + re2.2$parameter) + + re2.3 <- tibble::as_tibble(marginal_effects(surv_fit, all_contrasts = TRUE, type = "quantile", mtype = "difference")) + expect_identical(paste0("marg[", re2.3$.study, ": ", re2.3$.trtb, " vs. ", re2.3$.trta, ", ", re2.3$.quantile, "]"), + re2.3$parameter) + + + re3.1 <- tibble::as_tibble(marginal_effects(surv_fit, trt_ref = "Len", type = "median", mtype = "difference")) + expect_identical(paste0("marg[", re3.1$.study, ": ", re3.1$.trtb, "]"), + re3.1$parameter) + expect_identical(as.character(re3.1$.trta), rep("Len", times = nrow(re3.1))) + + re3.2 <- tibble::as_tibble(marginal_effects(surv_fit, trt_ref = "Len", type = "hazard", mtype = "ratio")) %>% + group_by(.study, .trtb, .trta) %>% + mutate(.id = 1:n()) + expect_identical(paste0("marg[", re3.2$.study, ": ", re3.2$.trtb, ", ", re3.2$.id, "]"), + re3.2$parameter) + expect_identical(as.character(re3.2$.trta), rep("Len", times = nrow(re3.2))) + + re3.3 <- tibble::as_tibble(marginal_effects(surv_fit, trt_ref = "Len", type = "quantile", mtype = "difference")) + expect_identical(paste0("marg[", re3.3$.study, ": ", re3.3$.trtb, ", ", re3.3$.quantile, "]"), + re3.3$parameter) + expect_identical(as.character(re3.3$.trta), rep("Len", times = nrow(re3.3))) +}) From c40a1fa2fb666b3f697dc45c5df734e3a89e5687 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Fri, 19 Apr 2024 13:18:56 +0100 Subject: [PATCH 31/57] Fix prior scales typo in dietary fat example (issue #34) --- tests/testthat/test-example_dietary_fat.R | 40 +-- vignettes/example_dietary_fat.Rmd | 8 +- vignettes/example_dietary_fat.html | 406 +++++++++++----------- 3 files changed, 231 insertions(+), 223 deletions(-) diff --git a/tests/testthat/test-example_dietary_fat.R b/tests/testthat/test-example_dietary_fat.R index 809f4aab..32fa64a5 100644 --- a/tests/testthat/test-example_dietary_fat.R +++ b/tests/testthat/test-example_dietary_fat.R @@ -8,17 +8,17 @@ skip_on_cran() params <- list(run_tests = FALSE) -## ----code=readLines("children/knitr_setup.R"), include=FALSE-------------------------------------- +## ---- code=readLines("children/knitr_setup.R"), include=FALSE------------------------------------- -## ----eval = FALSE--------------------------------------------------------------------------------- +## ---- eval = FALSE-------------------------------------------------------------------------------- ## library(multinma) ## options(mc.cores = parallel::detectCores()) ## ----setup, echo = FALSE-------------------------------------------------------------------------- library(multinma) -nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")), - "true" =, "warn" = 2, +nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")), + "true" =, "warn" = 2, parallel::detectCores()) options(mc.cores = nc) @@ -28,10 +28,10 @@ head(dietary_fat) ## ------------------------------------------------------------------------------------------------- -diet_net <- set_agd_arm(dietary_fat, +diet_net <- set_agd_arm(dietary_fat, study = studyc, trt = trtc, - r = r, + r = r, E = E, trt_ref = "Control", sample_size = n) @@ -43,7 +43,7 @@ summary(normal(scale = 100)) ## ------------------------------------------------------------------------------------------------- -diet_fit_FE <- nma(diet_net, +diet_fit_FE <- nma(diet_net, trt_effects = "fixed", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100)) @@ -53,7 +53,7 @@ diet_fit_FE <- nma(diet_net, diet_fit_FE -## ----eval=FALSE----------------------------------------------------------------------------------- +## ---- eval=FALSE---------------------------------------------------------------------------------- ## # Not run ## print(diet_fit_FE, pars = c("d", "mu")) @@ -67,19 +67,19 @@ summary(normal(scale = 100)) summary(half_normal(scale = 5)) -## ----eval=FALSE----------------------------------------------------------------------------------- +## ---- eval=FALSE---------------------------------------------------------------------------------- ## diet_fit_RE <- nma(diet_net, ## trt_effects = "random", -## prior_intercept = normal(scale = 10), -## prior_trt = normal(scale = 10), +## prior_intercept = normal(scale = 100), +## prior_trt = normal(scale = 100), ## prior_het = half_normal(scale = 5)) -## ----echo=FALSE, warning=FALSE-------------------------------------------------------------------- +## ---- echo=FALSE, warning=FALSE------------------------------------------------------------------- diet_fit_RE <- nowarn_on_ci( - nma(diet_net, + nma(diet_net, trt_effects = "random", - prior_intercept = normal(scale = 10), - prior_trt = normal(scale = 10), + prior_intercept = normal(scale = 100), + prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5)) ) @@ -88,7 +88,7 @@ diet_fit_RE <- nowarn_on_ci( diet_fit_RE -## ----eval=FALSE----------------------------------------------------------------------------------- +## ---- eval=FALSE---------------------------------------------------------------------------------- ## # Not run ## print(diet_fit_RE, pars = c("d", "mu", "delta")) @@ -113,15 +113,15 @@ plot(dic_RE) ## ----diet_pred_FE, fig.height = 2----------------------------------------------------------------- -pred_FE <- predict(diet_fit_FE, - baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5), +pred_FE <- predict(diet_fit_FE, + baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5), type = "response") pred_FE plot(pred_FE) ## ----diet_pred_RE, fig.height = 2----------------------------------------------------------------- -pred_RE <- predict(diet_fit_RE, - baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5), +pred_RE <- predict(diet_fit_RE, + baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5), type = "response") pred_RE plot(pred_RE) diff --git a/vignettes/example_dietary_fat.Rmd b/vignettes/example_dietary_fat.Rmd index d7b5929f..da73c62a 100644 --- a/vignettes/example_dietary_fat.Rmd +++ b/vignettes/example_dietary_fat.Rmd @@ -98,16 +98,16 @@ Fitting the RE model ```{r, eval=FALSE} diet_fit_RE <- nma(diet_net, trt_effects = "random", - prior_intercept = normal(scale = 10), - prior_trt = normal(scale = 10), + prior_intercept = normal(scale = 100), + prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5)) ``` ```{r, echo=FALSE, warning=FALSE} diet_fit_RE <- nowarn_on_ci( nma(diet_net, trt_effects = "random", - prior_intercept = normal(scale = 10), - prior_trt = normal(scale = 10), + prior_intercept = normal(scale = 100), + prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5)) ) ``` diff --git a/vignettes/example_dietary_fat.html b/vignettes/example_dietary_fat.html index 030038c9..296ace58 100644 --- a/vignettes/example_dietary_fat.html +++ b/vignettes/example_dietary_fat.html @@ -362,20 +362,28 @@

Example: Dietary fat

-
library(multinma)
-options(mc.cores = parallel::detectCores())
+
library(multinma)
+options(mc.cores = parallel::detectCores())
+
#> For execution on a local, multicore CPU with excess RAM we recommend calling
+#> options(mc.cores = parallel::detectCores())
+#> 
+#> Attaching package: 'multinma'
+#> The following objects are masked from 'package:stats':
+#> 
+#>     dgamma, pgamma, qgamma

This vignette describes the analysis of 10 trials comparing reduced fat diets to control (non-reduced fat diets) for preventing mortality -(Hooper et al. 2000; Dias et al. 2011). The data are +(Hooper et al. +2000; Dias et al. 2011). The data are available in this package as dietary_fat:

-
head(dietary_fat)
-#>   studyn            studyc trtn        trtc   r    n      E
-#> 1      1              DART    1     Control 113 1015 1917.0
-#> 2      1              DART    2 Reduced Fat 111 1018 1925.0
-#> 3      2 London Corn/Olive    1     Control   1   26   43.6
-#> 4      2 London Corn/Olive    2 Reduced Fat   5   28   41.3
-#> 5      2 London Corn/Olive    2 Reduced Fat   3   26   38.0
-#> 6      3    London Low Fat    1     Control  24  129  393.5
+
head(dietary_fat)
+#>   studyn            studyc trtn        trtc   r    n      E
+#> 1      1              DART    1     Control 113 1015 1917.0
+#> 2      1              DART    2 Reduced Fat 111 1018 1925.0
+#> 3      2 London Corn/Olive    1     Control   1   26   43.6
+#> 4      2 London Corn/Olive    2 Reduced Fat   5   28   41.3
+#> 5      2 London Corn/Olive    2 Reduced Fat   3   26   38.0
+#> 6      3    London Low Fat    1     Control  24  129  393.5

Setting up the network

We begin by setting up the network - here just a pairwise @@ -383,35 +391,35 @@

Setting up the network

(r) and the person-years at risk (E) in each arm, so we use the function set_agd_arm(). We set “Control” as the reference treatment.

-
diet_net <- set_agd_arm(dietary_fat, 
-                        study = studyc,
-                        trt = trtc,
-                        r = r, 
-                        E = E,
-                        trt_ref = "Control",
-                        sample_size = n)
-diet_net
-#> A network with 10 AgD studies (arm-based).
-#> 
-#> ------------------------------------------------------- AgD studies (arm-based) ---- 
-#>  Study                   Treatment arms                        
-#>  DART                    2: Control | Reduced Fat              
-#>  London Corn/Olive       3: Control | Reduced Fat | Reduced Fat
-#>  London Low Fat          2: Control | Reduced Fat              
-#>  Minnesota Coronary      2: Control | Reduced Fat              
-#>  MRC Soya                2: Control | Reduced Fat              
-#>  Oslo Diet-Heart         2: Control | Reduced Fat              
-#>  STARS                   2: Control | Reduced Fat              
-#>  Sydney Diet-Heart       2: Control | Reduced Fat              
-#>  Veterans Administration 2: Control | Reduced Fat              
-#>  Veterans Diet & Skin CA 2: Control | Reduced Fat              
-#> 
-#>  Outcome type: rate
-#> ------------------------------------------------------------------------------------
-#> Total number of treatments: 2
-#> Total number of studies: 10
-#> Reference treatment is: Control
-#> Network is connected
+
diet_net <- set_agd_arm(dietary_fat, 
+                        study = studyc,
+                        trt = trtc,
+                        r = r, 
+                        E = E,
+                        trt_ref = "Control",
+                        sample_size = n)
+diet_net
+#> A network with 10 AgD studies (arm-based).
+#> 
+#> ------------------------------------------------------- AgD studies (arm-based) ---- 
+#>  Study                   Treatment arms                        
+#>  DART                    2: Control | Reduced Fat              
+#>  London Corn/Olive       3: Control | Reduced Fat | Reduced Fat
+#>  London Low Fat          2: Control | Reduced Fat              
+#>  Minnesota Coronary      2: Control | Reduced Fat              
+#>  MRC Soya                2: Control | Reduced Fat              
+#>  Oslo Diet-Heart         2: Control | Reduced Fat              
+#>  STARS                   2: Control | Reduced Fat              
+#>  Sydney Diet-Heart       2: Control | Reduced Fat              
+#>  Veterans Administration 2: Control | Reduced Fat              
+#>  Veterans Diet & Skin CA 2: Control | Reduced Fat              
+#> 
+#>  Outcome type: rate
+#> ------------------------------------------------------------------------------------
+#> Total number of treatments: 2
+#> Total number of studies: 10
+#> Reference treatment is: Control
+#> Network is connected

We also specify the optional sample_size argument, although it is not strictly necessary here. In this case sample_size would only be required to produce a network @@ -433,41 +441,41 @@

Fixed effect meta-analysis

study-specific intercepts \(\mu_j\). We can examine the range of parameter values implied by these prior distributions with the summary() method:

-
summary(normal(scale = 100))
-#> A Normal prior distribution: location = 0, scale = 100.
-#> 50% of the prior density lies between -67.45 and 67.45.
-#> 95% of the prior density lies between -196 and 196.
+
summary(normal(scale = 100))
+#> A Normal prior distribution: location = 0, scale = 100.
+#> 50% of the prior density lies between -67.45 and 67.45.
+#> 95% of the prior density lies between -196 and 196.

The model is fitted using the nma() function. By default, this will use a Poisson likelihood with a log link function, auto-detected from the data.

-
diet_fit_FE <- nma(diet_net, 
-                   trt_effects = "fixed",
-                   prior_intercept = normal(scale = 100),
-                   prior_trt = normal(scale = 100))
+
diet_fit_FE <- nma(diet_net, 
+                   trt_effects = "fixed",
+                   prior_intercept = normal(scale = 100),
+                   prior_trt = normal(scale = 100))

Basic parameter summaries are given by the print() method:

-
diet_fit_FE
-#> A fixed effects NMA with a poisson likelihood (log link).
-#> Inference for Stan model: poisson.
-#> 4 chains, each with iter=2000; warmup=1000; thin=1; 
-#> post-warmup draws per chain=1000, total post-warmup draws=4000.
-#> 
-#>                   mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
-#> d[Reduced Fat]   -0.01    0.00 0.05   -0.11   -0.04   -0.01    0.03    0.10  3448    1
-#> lp__           5325.54    0.06 2.28 5320.39 5324.19 5325.82 5327.18 5329.17  1560    1
-#> 
-#> Samples were drawn using NUTS(diag_e) at Tue Jan  9 17:56:03 2024.
-#> For each parameter, n_eff is a crude measure of effective sample size,
-#> and Rhat is the potential scale reduction factor on split chains (at 
-#> convergence, Rhat=1).
+
diet_fit_FE
+#> A fixed effects NMA with a poisson likelihood (log link).
+#> Inference for Stan model: poisson.
+#> 4 chains, each with iter=2000; warmup=1000; thin=1; 
+#> post-warmup draws per chain=1000, total post-warmup draws=4000.
+#> 
+#>                   mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
+#> d[Reduced Fat]   -0.01    0.00 0.05   -0.11   -0.05   -0.01    0.03    0.09  3426    1
+#> lp__           5386.35    0.06 2.36 5380.89 5384.99 5386.66 5388.07 5389.89  1528    1
+#> 
+#> Samples were drawn using NUTS(diag_e) at Fri Apr 19 13:11:58 2024.
+#> For each parameter, n_eff is a crude measure of effective sample size,
+#> and Rhat is the potential scale reduction factor on split chains (at 
+#> convergence, Rhat=1).

By default, summaries of the study-specific intercepts \(\mu_j\) are hidden, but could be examined by changing the pars argument:

-
# Not run
-print(diet_fit_FE, pars = c("d", "mu"))
+
# Not run
+print(diet_fit_FE, pars = c("d", "mu"))

The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:

-
plot_prior_posterior(diet_fit_FE)
-

+
plot_prior_posterior(diet_fit_FE)
+

Random effects meta-analysis

@@ -479,169 +487,169 @@

Random effects meta-analysis

heterogeneity standard deviation \(\tau\). We can examine the range of parameter values implied by these prior distributions with the summary() method:

-
summary(normal(scale = 100))
-#> A Normal prior distribution: location = 0, scale = 100.
-#> 50% of the prior density lies between -67.45 and 67.45.
-#> 95% of the prior density lies between -196 and 196.
-summary(half_normal(scale = 5))
-#> A half-Normal prior distribution: location = 0, scale = 5.
-#> 50% of the prior density lies between 0 and 3.37.
-#> 95% of the prior density lies between 0 and 9.8.
+
summary(normal(scale = 100))
+#> A Normal prior distribution: location = 0, scale = 100.
+#> 50% of the prior density lies between -67.45 and 67.45.
+#> 95% of the prior density lies between -196 and 196.
+summary(half_normal(scale = 5))
+#> A half-Normal prior distribution: location = 0, scale = 5.
+#> 50% of the prior density lies between 0 and 3.37.
+#> 95% of the prior density lies between 0 and 9.8.

Fitting the RE model

-
diet_fit_RE <- nma(diet_net, 
-                   trt_effects = "random",
-                   prior_intercept = normal(scale = 10),
-                   prior_trt = normal(scale = 10),
-                   prior_het = half_normal(scale = 5))
+
diet_fit_RE <- nma(diet_net, 
+                   trt_effects = "random",
+                   prior_intercept = normal(scale = 100),
+                   prior_trt = normal(scale = 100),
+                   prior_het = half_normal(scale = 5))

Basic parameter summaries are given by the print() method:

-
diet_fit_RE
-#> A random effects NMA with a poisson likelihood (log link).
-#> Inference for Stan model: poisson.
-#> 4 chains, each with iter=2000; warmup=1000; thin=1; 
-#> post-warmup draws per chain=1000, total post-warmup draws=4000.
-#> 
-#>                   mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
-#> d[Reduced Fat]   -0.02    0.00 0.09   -0.20   -0.07   -0.02    0.03    0.15  2045    1
-#> lp__           5340.70    0.13 3.95 5332.26 5338.20 5340.96 5343.47 5347.61   971    1
-#> tau               0.13    0.00 0.11    0.01    0.05    0.10    0.18    0.41  1172    1
-#> 
-#> Samples were drawn using NUTS(diag_e) at Tue Jan  9 17:56:14 2024.
-#> For each parameter, n_eff is a crude measure of effective sample size,
-#> and Rhat is the potential scale reduction factor on split chains (at 
-#> convergence, Rhat=1).
+
diet_fit_RE
+#> A random effects NMA with a poisson likelihood (log link).
+#> Inference for Stan model: poisson.
+#> 4 chains, each with iter=2000; warmup=1000; thin=1; 
+#> post-warmup draws per chain=1000, total post-warmup draws=4000.
+#> 
+#>                   mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
+#> d[Reduced Fat]   -0.01    0.00 0.09   -0.20   -0.06   -0.01    0.04    0.18  1222    1
+#> lp__           5379.44    0.12 3.85 5371.02 5377.13 5379.69 5382.13 5386.15  1062    1
+#> tau               0.14    0.00 0.13    0.00    0.05    0.11    0.19    0.47   865    1
+#> 
+#> Samples were drawn using NUTS(diag_e) at Fri Apr 19 13:13:08 2024.
+#> For each parameter, n_eff is a crude measure of effective sample size,
+#> and Rhat is the potential scale reduction factor on split chains (at 
+#> convergence, Rhat=1).

By default, summaries of the study-specific intercepts \(\mu_j\) and study-specific relative effects \(\delta_{jk}\) are hidden, but could be examined by changing the pars argument:

-
# Not run
-print(diet_fit_RE, pars = c("d", "mu", "delta"))
+
# Not run
+print(diet_fit_RE, pars = c("d", "mu", "delta"))

The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:

-
plot_prior_posterior(diet_fit_RE, prior = c("trt", "het"))
-

+
plot_prior_posterior(diet_fit_RE, prior = c("trt", "het"))
+

Model comparison

Model fit can be checked using the dic() function:

-
(dic_FE <- dic(diet_fit_FE))
-#> Residual deviance: 22.2 (on 21 data points)
-#>                pD: 11
-#>               DIC: 33.2
-
(dic_RE <- dic(diet_fit_RE))
-#> Residual deviance: 21.4 (on 21 data points)
-#>                pD: 13.6
-#>               DIC: 35
+
(dic_FE <- dic(diet_fit_FE))
+#> Residual deviance: 22.1 (on 21 data points)
+#>                pD: 10.9
+#>               DIC: 33
+
(dic_RE <- dic(diet_fit_RE))
+#> Residual deviance: 21.2 (on 21 data points)
+#>                pD: 13.6
+#>               DIC: 34.9

Both models appear to fit the data well, as the residual deviance is close to the number of data points. The DIC is very similar between models, so the FE model may be preferred for parsimony.

We can also examine the residual deviance contributions with the corresponding plot() method.

-
plot(dic_FE)
-

-
plot(dic_RE)
-

+
plot(dic_FE)
+

+
plot(dic_RE)
+

Further results

-

Dias et al. (2011) produce absolute predictions of -the mortality rates on reduced fat and control diets, assuming a Normal +

Dias et al. (2011) produce absolute predictions of the +mortality rates on reduced fat and control diets, assuming a Normal distribution on the baseline log rate of mortality with mean \(-3\) and precision \(1.77\). We can replicate these results using the predict() method. The baseline argument takes a distr() distribution object, with which we specify the corresponding Normal distribution. We set type = "response" to produce predicted rates (type = "link" would produce predicted log rates).

-
pred_FE <- predict(diet_fit_FE, 
-                   baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5), 
-                   type = "response")
-pred_FE
-#>                   mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[Control]     0.07 0.06 0.01 0.03 0.05 0.08  0.21     4202     3972    1
-#> pred[Reduced Fat] 0.06 0.06 0.01 0.03 0.05 0.08  0.21     4224     3978    1
-plot(pred_FE)
-

-
pred_RE <- predict(diet_fit_RE, 
-                   baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5), 
-                   type = "response")
-pred_RE
-#>                   mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[Control]     0.07 0.06 0.01 0.03 0.05 0.08  0.22     4043     3537    1
-#> pred[Reduced Fat] 0.07 0.06 0.01 0.03 0.05 0.08  0.22     4020     3720    1
-plot(pred_RE)
-

+
pred_FE <- predict(diet_fit_FE, 
+                   baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5), 
+                   type = "response")
+pred_FE
+#>                   mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[Control]     0.07 0.06 0.01 0.03 0.05 0.08  0.23     3991     3835    1
+#> pred[Reduced Fat] 0.07 0.06 0.01 0.03 0.05 0.08  0.23     3922     3757    1
+plot(pred_FE)
+

+
pred_RE <- predict(diet_fit_RE, 
+                   baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5), 
+                   type = "response")
+pred_RE
+#>                   mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[Control]     0.06 0.05 0.01 0.03 0.05 0.08  0.21     4081     4077    1
+#> pred[Reduced Fat] 0.06 0.05 0.01 0.03 0.05 0.08  0.21     4057     3912    1
+plot(pred_RE)
+

If the baseline argument is omitted, predicted rates will be produced for every study in the network based on their estimated baseline log rate \(\mu_j\):

-
pred_FE_studies <- predict(diet_fit_FE, type = "response")
-pred_FE_studies
-#> ------------------------------------------------------------------- Study: DART ---- 
-#> 
-#>                         mean sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[DART: Control]     0.06  0 0.05 0.06 0.06 0.06  0.07     5543     2859    1
-#> pred[DART: Reduced Fat] 0.06  0 0.05 0.06 0.06 0.06  0.07     6370     3212    1
-#> 
-#> ------------------------------------------------------ Study: London Corn/Olive ---- 
-#> 
-#>                                      mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[London Corn/Olive: Control]     0.07 0.03 0.03 0.06 0.07 0.09  0.13     7039     2639    1
-#> pred[London Corn/Olive: Reduced Fat] 0.07 0.02 0.03 0.05 0.07 0.09  0.13     7254     2690    1
-#> 
-#> --------------------------------------------------------- Study: London Low Fat ---- 
-#> 
-#>                                   mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[London Low Fat: Control]     0.06 0.01 0.04 0.05 0.06 0.06  0.08     8037     2986    1
-#> pred[London Low Fat: Reduced Fat] 0.06 0.01 0.04 0.05 0.06 0.06  0.08     9224     3205    1
-#> 
-#> ----------------------------------------------------- Study: Minnesota Coronary ---- 
-#> 
-#>                                       mean sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[Minnesota Coronary: Control]     0.05  0 0.05 0.05 0.05 0.06  0.06     4832     3586    1
-#> pred[Minnesota Coronary: Reduced Fat] 0.05  0 0.05 0.05 0.05 0.06  0.06     6395     3547    1
-#> 
-#> --------------------------------------------------------------- Study: MRC Soya ---- 
-#> 
-#>                             mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[MRC Soya: Control]     0.04 0.01 0.03 0.04 0.04 0.04  0.05     6442     2770    1
-#> pred[MRC Soya: Reduced Fat] 0.04 0.01 0.03 0.04 0.04 0.04  0.05     7057     2832    1
-#> 
-#> -------------------------------------------------------- Study: Oslo Diet-Heart ---- 
-#> 
-#>                                    mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[Oslo Diet-Heart: Control]     0.06 0.01 0.05 0.06 0.06 0.07  0.08     6028     2675    1
-#> pred[Oslo Diet-Heart: Reduced Fat] 0.06 0.01 0.05 0.06 0.06 0.07  0.08     8013     3279    1
-#> 
-#> ------------------------------------------------------------------ Study: STARS ---- 
-#> 
-#>                          mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[STARS: Control]     0.02 0.01 0.01 0.01 0.02 0.03  0.05     7028     2787    1
-#> pred[STARS: Reduced Fat] 0.02 0.01 0.01 0.01 0.02 0.03  0.05     6997     2750    1
-#> 
-#> ------------------------------------------------------ Study: Sydney Diet-Heart ---- 
-#> 
-#>                                      mean sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
-#> pred[Sydney Diet-Heart: Control]     0.03  0 0.03 0.03 0.03 0.04  0.04     6388     3364    1
-#> pred[Sydney Diet-Heart: Reduced Fat] 0.03  0 0.03 0.03 0.03 0.04  0.04     7326     3336    1
-#> 
-#> ------------------------------------------------ Study: Veterans Administration ---- 
-#> 
-#>                                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS
-#> pred[Veterans Administration: Control]     0.11 0.01  0.1 0.11 0.11 0.12  0.13     4882     3030
-#> pred[Veterans Administration: Reduced Fat] 0.11 0.01  0.1 0.11 0.11 0.12  0.13     7066     3282
-#>                                            Rhat
-#> pred[Veterans Administration: Control]        1
-#> pred[Veterans Administration: Reduced Fat]    1
-#> 
-#> ------------------------------------------------ Study: Veterans Diet & Skin CA ---- 
-#> 
-#>                                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS
-#> pred[Veterans Diet & Skin CA: Control]     0.01 0.01    0 0.01 0.01 0.02  0.03     6770     2728
-#> pred[Veterans Diet & Skin CA: Reduced Fat] 0.01 0.01    0 0.01 0.01 0.02  0.03     6895     2652
-#>                                            Rhat
-#> pred[Veterans Diet & Skin CA: Control]        1
-#> pred[Veterans Diet & Skin CA: Reduced Fat]    1
-plot(pred_FE_studies) + ggplot2::facet_grid(Study~., labeller = ggplot2::label_wrap_gen(width = 10))
-

+
pred_FE_studies <- predict(diet_fit_FE, type = "response")
+pred_FE_studies
+#> ------------------------------------------------------------------- Study: DART ---- 
+#> 
+#>                         mean sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[DART: Control]     0.06  0 0.05 0.06 0.06 0.06  0.07     5601     3098    1
+#> pred[DART: Reduced Fat] 0.06  0 0.05 0.06 0.06 0.06  0.07     7752     3176    1
+#> 
+#> ------------------------------------------------------ Study: London Corn/Olive ---- 
+#> 
+#>                                      mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[London Corn/Olive: Control]     0.07 0.02 0.03 0.06 0.07 0.09  0.13     7183     2988 1.01
+#> pred[London Corn/Olive: Reduced Fat] 0.07 0.02 0.03 0.06 0.07 0.09  0.13     7121     3035 1.01
+#> 
+#> --------------------------------------------------------- Study: London Low Fat ---- 
+#> 
+#>                                   mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[London Low Fat: Control]     0.06 0.01 0.04 0.05 0.06 0.06  0.08     6555     2868    1
+#> pred[London Low Fat: Reduced Fat] 0.06 0.01 0.04 0.05 0.06 0.06  0.08     7683     3049    1
+#> 
+#> ----------------------------------------------------- Study: Minnesota Coronary ---- 
+#> 
+#>                                       mean sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[Minnesota Coronary: Control]     0.05  0 0.05 0.05 0.05 0.06  0.06     4396     3458    1
+#> pred[Minnesota Coronary: Reduced Fat] 0.05  0 0.05 0.05 0.05 0.06  0.06     8137     3772    1
+#> 
+#> --------------------------------------------------------------- Study: MRC Soya ---- 
+#> 
+#>                             mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[MRC Soya: Control]     0.04 0.01 0.03 0.04 0.04 0.04  0.05     6224     3055    1
+#> pred[MRC Soya: Reduced Fat] 0.04 0.01 0.03 0.04 0.04 0.04  0.05     7525     3247    1
+#> 
+#> -------------------------------------------------------- Study: Oslo Diet-Heart ---- 
+#> 
+#>                                    mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[Oslo Diet-Heart: Control]     0.06 0.01 0.05 0.06 0.06 0.07  0.08     6372     3078    1
+#> pred[Oslo Diet-Heart: Reduced Fat] 0.06 0.01 0.05 0.06 0.06 0.07  0.08     7681     3100    1
+#> 
+#> ------------------------------------------------------------------ Study: STARS ---- 
+#> 
+#>                          mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[STARS: Control]     0.02 0.01 0.01 0.01 0.02 0.03  0.05     7184     2361    1
+#> pred[STARS: Reduced Fat] 0.02 0.01 0.01 0.01 0.02 0.03  0.05     7356     2306    1
+#> 
+#> ------------------------------------------------------ Study: Sydney Diet-Heart ---- 
+#> 
+#>                                      mean sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
+#> pred[Sydney Diet-Heart: Control]     0.03  0 0.03 0.03 0.03 0.04  0.04     6935     3061    1
+#> pred[Sydney Diet-Heart: Reduced Fat] 0.03  0 0.03 0.03 0.03 0.04  0.04     7792     3483    1
+#> 
+#> ------------------------------------------------ Study: Veterans Administration ---- 
+#> 
+#>                                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS
+#> pred[Veterans Administration: Control]     0.11 0.01  0.1 0.11 0.11 0.12  0.13     5552     3154
+#> pred[Veterans Administration: Reduced Fat] 0.11 0.01  0.1 0.11 0.11 0.12  0.12     7289     3471
+#>                                            Rhat
+#> pred[Veterans Administration: Control]        1
+#> pred[Veterans Administration: Reduced Fat]    1
+#> 
+#> ------------------------------------------------ Study: Veterans Diet & Skin CA ---- 
+#> 
+#>                                            mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS
+#> pred[Veterans Diet & Skin CA: Control]     0.01 0.01    0 0.01 0.01 0.02  0.03     6465     2903
+#> pred[Veterans Diet & Skin CA: Reduced Fat] 0.01 0.01    0 0.01 0.01 0.02  0.03     6548     2765
+#>                                            Rhat
+#> pred[Veterans Diet & Skin CA: Control]        1
+#> pred[Veterans Diet & Skin CA: Reduced Fat]    1
+plot(pred_FE_studies) + ggplot2::facet_grid(Study~., labeller = ggplot2::label_wrap_gen(width = 10))
+

References

From 8bed26f5555895aa560c217412144187e1821232 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 24 Apr 2024 11:56:58 +0100 Subject: [PATCH 32/57] Add visual tests for network plots --- DESCRIPTION | 1 + .../atrial-fibrillation-network-star.svg | 107 ++++++ .../atrial-fibrillation-network.svg | 125 +++++++ .../plot_nma_data/bcg-vaccine-network.svg | 41 +++ .../_snaps/plot_nma_data/blocker-network.svg | 41 +++ .../_snaps/plot_nma_data/diabetes-network.svg | 81 ++++ .../plot_nma_data/dietary-fat-network.svg | 41 +++ .../duplicated-agd-arms-network.svg | 41 +++ .../duplicated-agd-contrasts-network.svg | 41 +++ .../duplicated-ipd-arms-network.svg | 41 +++ .../plot_nma_data/hta-psoriasis-network.svg | 75 ++++ .../_snaps/plot_nma_data/ndmm-network.svg | 61 +++ .../plot_nma_data/parkinsons-arm-network.svg | 65 ++++ .../parkinsons-contrast-network.svg | 65 ++++ .../parkinsons-mixed-network.svg | 65 ++++ .../plaque-psoriasis-network.svg | 93 +++++ .../_snaps/plot_nma_data/smoking-network.svg | 61 +++ .../_snaps/plot_nma_data/statins-network.svg | 41 +++ .../thrombolytics-dias-network.svg | 83 +++++ .../plot_nma_data/thrombolytics-network.svg | 83 +++++ .../plot_nma_data/transfusion-network.svg | 41 +++ tests/testthat/test-plot_nma_data.R | 348 ++++++++++++++++++ 22 files changed, 1641 insertions(+) create mode 100644 tests/testthat/_snaps/plot_nma_data/atrial-fibrillation-network-star.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/atrial-fibrillation-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/bcg-vaccine-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/blocker-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/diabetes-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/dietary-fat-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/duplicated-agd-arms-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/duplicated-agd-contrasts-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/duplicated-ipd-arms-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/hta-psoriasis-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/ndmm-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/parkinsons-arm-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/parkinsons-contrast-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/parkinsons-mixed-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/plaque-psoriasis-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/smoking-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/statins-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/thrombolytics-dias-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/thrombolytics-network.svg create mode 100644 tests/testthat/_snaps/plot_nma_data/transfusion-network.svg create mode 100644 tests/testthat/test-plot_nma_data.R diff --git a/DESCRIPTION b/DESCRIPTION index 307340c5..befaf63f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -57,6 +57,7 @@ RoxygenNote: 7.2.3 Roxygen: list(markdown = TRUE) Suggests: testthat (>= 2.1.0), + vdiffr, withr, knitr, rmarkdown, diff --git a/tests/testthat/_snaps/plot_nma_data/atrial-fibrillation-network-star.svg b/tests/testthat/_snaps/plot_nma_data/atrial-fibrillation-network-star.svg new file mode 100644 index 00000000..cc7cfa72 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/atrial-fibrillation-network-star.svg @@ -0,0 +1,107 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Standard adjusted dose anti-coagulant + +Acenocoumarol + +Alternate day aspirin + +Dipyridamole + +Fixed dose warfarin + +Fixed dose warfarin + low dose aspirin + +Fixed dose warfarin + medium dose aspirin + +High dose aspirin + +Indobufen + +Low adjusted dose anti-coagulant + +Low dose aspirin + +Low dose aspirin + copidogrel + +Low dose aspirin + dipyridamole + +Medium dose aspirin + +Placebo/Standard care + +Triflusal + +Ximelagatran +Number of studies + + + + + +1 +2 +3 +4 +5 +Atrial fibrillation network star + + diff --git a/tests/testthat/_snaps/plot_nma_data/atrial-fibrillation-network.svg b/tests/testthat/_snaps/plot_nma_data/atrial-fibrillation-network.svg new file mode 100644 index 00000000..5167ca83 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/atrial-fibrillation-network.svg @@ -0,0 +1,125 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Standard adjusted dose anti-coagulant +Acenocoumarol +Alternate day aspirin +Dipyridamole +Fixed dose warfarin +Fixed dose warfarin + low dose aspirin +Fixed dose warfarin + medium dose aspirin +High dose aspirin +Indobufen +Low adjusted dose anti-coagulant +Low dose aspirin +Low dose aspirin + copidogrel +Low dose aspirin + dipyridamole +Medium dose aspirin +Placebo/Standard care +Triflusal +Ximelagatran +Treatment Class + + + + +Anti-coagulant +Anti-platelet +Control +Mixed +Number of studies + + + + + +1 +2 +3 +4 +5 +Total sample size + + + + +2500 +5000 +7500 +10000 +Atrial fibrillation network + + diff --git a/tests/testthat/_snaps/plot_nma_data/bcg-vaccine-network.svg b/tests/testthat/_snaps/plot_nma_data/bcg-vaccine-network.svg new file mode 100644 index 00000000..1b7c705d --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/bcg-vaccine-network.svg @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + +Unvaccinated + +Vaccinated +Number of studies + +13 +BCG vaccine network + + diff --git a/tests/testthat/_snaps/plot_nma_data/blocker-network.svg b/tests/testthat/_snaps/plot_nma_data/blocker-network.svg new file mode 100644 index 00000000..7ed0e20c --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/blocker-network.svg @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + +Control + +Beta Blocker +Number of studies + +22 +Blocker network + + diff --git a/tests/testthat/_snaps/plot_nma_data/diabetes-network.svg b/tests/testthat/_snaps/plot_nma_data/diabetes-network.svg new file mode 100644 index 00000000..0b903c13 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/diabetes-network.svg @@ -0,0 +1,81 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Beta Blocker +ACE Inhibitor +ARB +CCB +Diuretic +Placebo +Total sample size + + + + + +15000 +20000 +25000 +30000 +35000 +Number of studies + + + + + +1 +2 +3 +4 +5 +Diabetes network + + diff --git a/tests/testthat/_snaps/plot_nma_data/dietary-fat-network.svg b/tests/testthat/_snaps/plot_nma_data/dietary-fat-network.svg new file mode 100644 index 00000000..3257d1f8 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/dietary-fat-network.svg @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + +Control + +Reduced Fat +Number of studies + +11 +Dietary fat network + + diff --git a/tests/testthat/_snaps/plot_nma_data/duplicated-agd-arms-network.svg b/tests/testthat/_snaps/plot_nma_data/duplicated-agd-arms-network.svg new file mode 100644 index 00000000..5240a4c5 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/duplicated-agd-arms-network.svg @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + +2 + +1 +Number of studies + +3 +Duplicated AgD arms network + + diff --git a/tests/testthat/_snaps/plot_nma_data/duplicated-agd-contrasts-network.svg b/tests/testthat/_snaps/plot_nma_data/duplicated-agd-contrasts-network.svg new file mode 100644 index 00000000..5c8083a9 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/duplicated-agd-contrasts-network.svg @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + +2 + +1 +Number of studies + +3 +Duplicated AgD contrasts network + + diff --git a/tests/testthat/_snaps/plot_nma_data/duplicated-ipd-arms-network.svg b/tests/testthat/_snaps/plot_nma_data/duplicated-ipd-arms-network.svg new file mode 100644 index 00000000..4854038e --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/duplicated-ipd-arms-network.svg @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + +1 + +2 +Number of studies + +1 +Duplicated IPD arms network + + diff --git a/tests/testthat/_snaps/plot_nma_data/hta-psoriasis-network.svg b/tests/testthat/_snaps/plot_nma_data/hta-psoriasis-network.svg new file mode 100644 index 00000000..223a1a89 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/hta-psoriasis-network.svg @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Supportive care +Ciclosporin +Efalizumab +Etanercept 25 mg +Etanercept 50 mg +Fumaderm +Infliximab +Methotrexate +Total sample size + + + +500 +1000 +1500 +Number of studies + + + + + +1 +2 +3 +4 +5 +HTA psoriasis network + + diff --git a/tests/testthat/_snaps/plot_nma_data/ndmm-network.svg b/tests/testthat/_snaps/plot_nma_data/ndmm-network.svg new file mode 100644 index 00000000..1780a131 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/ndmm-network.svg @@ -0,0 +1,61 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Pbo +Len +Thal +Total sample size + + + +800 +1200 +1600 +Number of studies + + + +1 +2 +3 +Data + + +AgD +IPD +NDMM network + + diff --git a/tests/testthat/_snaps/plot_nma_data/parkinsons-arm-network.svg b/tests/testthat/_snaps/plot_nma_data/parkinsons-arm-network.svg new file mode 100644 index 00000000..5b3e80b2 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/parkinsons-arm-network.svg @@ -0,0 +1,65 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +4 +1 +2 +3 +5 +Total sample size + + + + + +250 +300 +350 +400 +450 +Number of studies + + +1 +2 +Parkinsons arm network + + diff --git a/tests/testthat/_snaps/plot_nma_data/parkinsons-contrast-network.svg b/tests/testthat/_snaps/plot_nma_data/parkinsons-contrast-network.svg new file mode 100644 index 00000000..8ddae391 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/parkinsons-contrast-network.svg @@ -0,0 +1,65 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +4 +1 +2 +3 +5 +Total sample size + + + + + +250 +300 +350 +400 +450 +Number of studies + + +1 +2 +Parkinsons contrast network + + diff --git a/tests/testthat/_snaps/plot_nma_data/parkinsons-mixed-network.svg b/tests/testthat/_snaps/plot_nma_data/parkinsons-mixed-network.svg new file mode 100644 index 00000000..32b32586 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/parkinsons-mixed-network.svg @@ -0,0 +1,65 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +4 +1 +2 +3 +5 +Total sample size + + + + + +250 +300 +350 +400 +450 +Number of studies + + +1 +2 +Parkinsons mixed network + + diff --git a/tests/testthat/_snaps/plot_nma_data/plaque-psoriasis-network.svg b/tests/testthat/_snaps/plot_nma_data/plaque-psoriasis-network.svg new file mode 100644 index 00000000..92352b5b --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/plaque-psoriasis-network.svg @@ -0,0 +1,93 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +PBO +ETN +IXE_Q2W +IXE_Q4W +SEC_150 +SEC_300 +UST +Total sample size + + + + +500 +750 +1000 +1250 +Data + + +AgD +IPD +Number of studies + + + + +1 +2 +3 +4 +Treatment class + + + + +IL-12/23 blocker +IL-17 blocker +Placebo +TNFa blocker +Plaque psoriasis network + + diff --git a/tests/testthat/_snaps/plot_nma_data/smoking-network.svg b/tests/testthat/_snaps/plot_nma_data/smoking-network.svg new file mode 100644 index 00000000..8d421443 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/smoking-network.svg @@ -0,0 +1,61 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +No intervention +Group counselling +Individual counselling +Self-help +Number of studies + + + +2 +9 +16 +Total sample size + + + +2000 +4000 +6000 +Smoking network + + diff --git a/tests/testthat/_snaps/plot_nma_data/statins-network.svg b/tests/testthat/_snaps/plot_nma_data/statins-network.svg new file mode 100644 index 00000000..e8405361 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/statins-network.svg @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + +Placebo + +Statin +Number of studies + +19 +Statins network + + diff --git a/tests/testthat/_snaps/plot_nma_data/thrombolytics-dias-network.svg b/tests/testthat/_snaps/plot_nma_data/thrombolytics-dias-network.svg new file mode 100644 index 00000000..9fe231aa --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/thrombolytics-dias-network.svg @@ -0,0 +1,83 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +SK +t-PA +Acc t-PA +SK + t-PA +r-PA +TNK +PTCA +UK +ASPAC +Number of studies + + + +1 +6 +11 +Total sample size + + + + +10000 +20000 +30000 +40000 +Thrombolytics Dias network + + diff --git a/tests/testthat/_snaps/plot_nma_data/thrombolytics-network.svg b/tests/testthat/_snaps/plot_nma_data/thrombolytics-network.svg new file mode 100644 index 00000000..7f8b31c7 --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/thrombolytics-network.svg @@ -0,0 +1,83 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +SK +Acc t-PA +ASPAC +PTCA +r-PA +SK + t-PA +t-PA +TNK +UK +Number of studies + + + +1 +6 +11 +Total sample size + + + + +10000 +20000 +30000 +40000 +Thrombolytics network + + diff --git a/tests/testthat/_snaps/plot_nma_data/transfusion-network.svg b/tests/testthat/_snaps/plot_nma_data/transfusion-network.svg new file mode 100644 index 00000000..86ed847a --- /dev/null +++ b/tests/testthat/_snaps/plot_nma_data/transfusion-network.svg @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + +Control + +Transfusion +Number of studies + +6 +Transfusion network + + diff --git a/tests/testthat/test-plot_nma_data.R b/tests/testthat/test-plot_nma_data.R new file mode 100644 index 00000000..c0ed8fcd --- /dev/null +++ b/tests/testthat/test-plot_nma_data.R @@ -0,0 +1,348 @@ + +skip_on_cran() + +library(vdiffr) +library(ggplot2) +library(dplyr) + +test_that("Atrial fibrillation", { + af_net <- set_agd_arm(atrial_fibrillation[atrial_fibrillation$studyc != "WASPO", ], + study = studyc, + trt = trtc, + r = r, + n = n, + trt_class = trt_class) + + expect_doppelganger("Atrial fibrillation network", + plot(af_net, weight_nodes = TRUE, weight_edges = TRUE, show_trt_class = TRUE) + + ggplot2::theme(legend.position = "bottom", legend.box = "vertical") + ) + + expect_doppelganger("Atrial fibrillation network star", + plot(af_net, layout = "star") + ) +}) + +test_that("BCG vaccine", { + bcg_net <- set_agd_arm(bcg_vaccine, + study = studyn, + trt = trtc, + r = r, + n = n, + trt_ref = "Unvaccinated") + + expect_doppelganger("BCG vaccine network", plot(bcg_net)) +}) + +test_that("Blocker", { + blocker_net <- set_agd_arm(blocker, + study = studyn, + trt = trtc, + r = r, + n = n, + trt_ref = "Control") + + expect_doppelganger("Blocker network", plot(blocker_net)) +}) + +test_that("Diabetes", { + db_net <- set_agd_arm(diabetes, + study = studyc, + trt = trtc, + r = r, + n = n) + expect_doppelganger("Diabetes network", + plot(db_net, weight_edges = TRUE, weight_nodes = TRUE) + ggplot2::theme(legend.box.margin = ggplot2::unit(c(0, 0, 0, 4), "lines")) + ) +}) + +test_that("Dietary fat", { + diet_net <- set_agd_arm(dietary_fat, + study = studyc, + trt = trtc, + r = r, + E = E, + trt_ref = "Control", + sample_size = n) + + expect_doppelganger("Dietary fat network", plot(diet_net)) +}) + +test_that("HTA psoriasis", { + pso_net <- set_agd_arm(hta_psoriasis, + study = paste(studyc, year), + trt = trtc, + r = multi(r0 = sample_size - rowSums(cbind(PASI50, PASI75, PASI90), na.rm = TRUE), + PASI50, PASI75, PASI90, + inclusive = FALSE, + type = "ordered")) + + expect_doppelganger("HTA psoriasis network", + plot(pso_net, weight_edges = TRUE, weight_nodes = TRUE) + + # Nudge the legend over + ggplot2::theme(legend.box.spacing = ggplot2::unit(0.75, "in"), + plot.margin = ggplot2::margin(0.1, 0, 0.1, 0.75, "in")) + ) +}) + +test_that("Parkinsons", { + + arm_net <- set_agd_arm(parkinsons, + study = studyn, + trt = trtn, + y = y, + se = se, + sample_size = n) + + expect_doppelganger("Parkinsons arm network", + plot(arm_net, weight_edges = TRUE, weight_nodes = TRUE) + ) + + contr_net <- set_agd_contrast(parkinsons, + study = studyn, + trt = trtn, + y = diff, + se = se_diff, + sample_size = n) + + expect_doppelganger("Parkinsons contrast network", + plot(contr_net, weight_edges = TRUE, weight_nodes = TRUE) + ) + + studies <- parkinsons$studyn + (parkinsons_arm <- parkinsons[studies %in% 1:3, ]) + (parkinsons_contr <- parkinsons[studies %in% 4:7, ]) + + mix_arm_net <- set_agd_arm(parkinsons_arm, + study = studyn, + trt = trtn, + y = y, + se = se, + sample_size = n) + + mix_contr_net <- set_agd_contrast(parkinsons_contr, + study = studyn, + trt = trtn, + y = diff, + se = se_diff, + sample_size = n) + + mix_net <- combine_network(mix_arm_net, mix_contr_net) + + expect_doppelganger("Parkinsons mixed network", + plot(mix_net, weight_edges = TRUE, weight_nodes = TRUE) + ) +}) + +test_that("Plaque psoriasis", { + + # IPD studies + pso_ipd <- plaque_psoriasis_ipd %>% + mutate( + # Variable transformations + bsa = bsa / 100, + weight = weight / 10, + durnpso = durnpso / 10, + prevsys = as.numeric(prevsys), + psa = as.numeric(psa), + # Treatment classes + trtclass = case_when(trtn == 1 ~ "Placebo", + trtn %in% c(2, 3, 5, 6) ~ "IL-17 blocker", + trtn == 4 ~ "TNFa blocker", + trtn == 7 ~ "IL-12/23 blocker"), + # Check complete cases for covariates of interest + is_complete = complete.cases(durnpso, prevsys, bsa, weight, psa) + ) %>% + arrange(studyc, trtn) %>% + filter(is_complete) + + # AgD studies + pso_agd <- plaque_psoriasis_agd %>% + mutate( + # Variable transformations + bsa_mean = bsa_mean / 100, + bsa_sd = bsa_sd / 100, + weight_mean = weight_mean / 10, + weight_sd = weight_sd / 10, + durnpso_mean = durnpso_mean / 10, + durnpso_sd = durnpso_sd / 10, + prevsys = prevsys / 100, + psa = psa / 100, + # Treatment classes + trtclass = case_when(trtn == 1 ~ "Placebo", + trtn %in% c(2, 3, 5, 6) ~ "IL-17 blocker", + trtn == 4 ~ "TNFa blocker", + trtn == 7 ~ "IL-12/23 blocker") + ) %>% + arrange(studyc, trtn) + + pso_net <- combine_network( + set_ipd(pso_ipd, + study = studyc, + trt = trtc, + r = multi(r0 = 1, + PASI75 = pasi75, + PASI90 = pasi90, + PASI100 = pasi100, + type = "ordered", inclusive = TRUE), + trt_class = trtclass), + set_agd_arm(pso_agd, + study = studyc, + trt = trtc, + r = multi(r0 = pasi75_n, + PASI75 = pasi75_r, + PASI90 = pasi90_r, + PASI100 = pasi100_r, + type = "ordered", inclusive = TRUE), + trt_class = trtclass) + ) + + + class_pal <- c("#D95F02", "#7570B3", "#E7298A", "#E6AB02") + + + expect_doppelganger("Plaque psoriasis network", + plot(pso_net, weight_nodes = TRUE, weight_edges = TRUE, show_trt_class = TRUE, nudge = 0.1) + + ggraph::scale_edge_colour_manual("Data", + values = c(AgD = "#113259", IPD = "#55A480")) + + scale_fill_manual("Treatment class", + values = class_pal, + aesthetics = c("fill", "colour")) + + guides(edge_colour = guide_legend(override.aes = list(edge_width = 2)), + fill = guide_legend(override.aes = list(size = 2))) + ) +}) + +test_that("Smoking", { + + smknet <- set_agd_arm(smoking, + study = studyn, + trt = trtc, + r = r, + n = n, + trt_ref = "No intervention") + + expect_doppelganger("Smoking network", + plot(smknet, weight_edges = TRUE, weight_nodes = TRUE) + ggplot2::theme(plot.margin = ggplot2::unit(c(1, 1, 1, 6), "lines")) + ) +}) + +test_that("Statins", { + + statin_net <- set_agd_arm(statins, + study = studyc, + trt = trtc, + r = r, + n = n, + trt_ref = "Placebo") + + expect_doppelganger("Statins network", + plot(statin_net) + ) +}) + +test_that("Thrombolytics", { + + thrombo_net <- set_agd_arm(thrombolytics, + study = studyn, + trt = trtc, + r = r, + n = n) + + expect_doppelganger("Thrombolytics network", + plot(thrombo_net, weight_edges = TRUE, weight_nodes = TRUE) + ) + + # Make trtf factor to order treatments in same way as Dias analysis + trts <- dplyr::distinct(thrombolytics, trtn, trtc) + trts <- dplyr::arrange(trts, trtn) + thrombolytics$trtf <- factor(thrombolytics$trtn, levels = trts$trtn, labels = trts$trtc) + thrombo_net2 <- set_agd_arm(thrombolytics, + study = studyn, + trt = trtf, + r = r, + n = n) + + expect_doppelganger("Thrombolytics Dias network", + plot(thrombo_net2, weight_edges = TRUE, weight_nodes = TRUE) + ) +}) + +test_that("Transfusion", { + + tr_net <- set_agd_arm(transfusion, + study = studyc, + trt = trtc, + r = r, + n = n, + trt_ref = "Control") + + expect_doppelganger("Transfusion network", + plot(tr_net) + ) +}) + +test_that("NDMM", { + + ndmm_ipd$trtclass <- case_match(ndmm_ipd$trtf, + "Pbo" ~ "Placebo", + c("Len", "Thal") ~ "Active") + + ndmm_agd$trtclass <- case_match(ndmm_agd$trtf, + "Pbo" ~ "Placebo", + c("Len", "Thal") ~ "Active") + + ndmm_net <- combine_network( + set_ipd(ndmm_ipd, + study = studyf, + trt = trtf, + trt_class = trtclass, + Surv = Surv(eventtime, status)), + set_agd_surv(ndmm_agd, + study = studyf, + trt = trtf, + trt_class = trtclass, + Surv = Surv(eventtime, status), + covariates = ndmm_agd_covs) + ) + + ndmm_net <- add_integration(ndmm_net, + age = distr(qgamma, mean = age_mean, sd = age_sd), + iss_stage3 = distr(qbern, iss_stage3), + response_cr_vgpr = distr(qbern, response_cr_vgpr), + male = distr(qbern, male)) + + expect_doppelganger("NDMM network", + plot(ndmm_net, + weight_nodes = TRUE, + weight_edges = TRUE, + # Nudge treatment labels away from nodes + nudge = 0.1, + # Manual layout + layout = data.frame(x = c(0, -1, 1), + y = c(-0.5, 0, 0))) + + guides(edge_colour = guide_legend(override.aes = list(edge_width = 2))) + + theme(legend.position = "bottom", legend.direction = "vertical") + ) +}) + +test_that("Duplicated study arms plotted correctly", { + + dat <- data.frame(study = "a", + trt = c(1, 2, 2, 2), + r = 1, n = 1) + + net1 <- set_agd_arm(dat, study, trt, r = r, n = n) + expect_doppelganger("Duplicated AgD arms network", plot(net1)) + + + net2 <- set_ipd(dat, study, trt, r = r) + expect_doppelganger("Duplicated IPD arms network", plot(net2)) + + dat_c <- data.frame(study = "a", + trt = c(1, 2, 2, 2), + y = c(NA, 1, 1, 1), + se = c(0.5, 1, 1, 1)) + net3 <- set_agd_contrast(dat_c, study, trt, y = y, se = se) + expect_doppelganger("Duplicated AgD contrasts network", plot(net3)) +}) From db2330412de7c2ed46e26dc18891776c74fb6540 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 24 Apr 2024 11:58:04 +0100 Subject: [PATCH 33/57] Correct nstudy in network plots with duplicated treatment arms --- R/nma_data-class.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/nma_data-class.R b/R/nma_data-class.R index ef1dc755..fc7a9e7b 100644 --- a/R/nma_data-class.R +++ b/R/nma_data-class.R @@ -279,7 +279,7 @@ as.igraph.nma_data <- function(x, ..., collapse = TRUE) { if (collapse) { e_ipd <- e_ipd %>% dplyr::group_by(.data$.trt, .data$.trt_b) %>% - dplyr::summarise(.nstudy = dplyr::n(), .type = "IPD") + dplyr::summarise(.nstudy = dplyr::n_distinct(.data$.study), .type = "IPD") } else { e_ipd$.type <- "IPD" } @@ -300,7 +300,7 @@ as.igraph.nma_data <- function(x, ..., collapse = TRUE) { if (collapse) { e_agd <- e_agd %>% dplyr::group_by(.data$.trt, .data$.trt_b) %>% - dplyr::summarise(.nstudy = dplyr::n(), .type = "AgD") + dplyr::summarise(.nstudy = dplyr::n_distinct(.data$.study), .type = "AgD") } else { e_agd$.type <- "AgD" } From 8e0347412dac667b174f98c817033ae178327e3f Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 24 Apr 2024 12:03:47 +0100 Subject: [PATCH 34/57] Confirm snapshots of corrected plots --- tests/testthat/_snaps/plot_nma_data/dietary-fat-network.svg | 2 +- .../_snaps/plot_nma_data/duplicated-agd-arms-network.svg | 2 +- .../_snaps/plot_nma_data/duplicated-agd-contrasts-network.svg | 2 +- tests/testthat/_snaps/plot_nma_data/hta-psoriasis-network.svg | 4 ++-- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/testthat/_snaps/plot_nma_data/dietary-fat-network.svg b/tests/testthat/_snaps/plot_nma_data/dietary-fat-network.svg index 3257d1f8..059f3ebb 100644 --- a/tests/testthat/_snaps/plot_nma_data/dietary-fat-network.svg +++ b/tests/testthat/_snaps/plot_nma_data/dietary-fat-network.svg @@ -35,7 +35,7 @@ Reduced Fat Number of studies -11 +10 Dietary fat network diff --git a/tests/testthat/_snaps/plot_nma_data/duplicated-agd-arms-network.svg b/tests/testthat/_snaps/plot_nma_data/duplicated-agd-arms-network.svg index 5240a4c5..46efef17 100644 --- a/tests/testthat/_snaps/plot_nma_data/duplicated-agd-arms-network.svg +++ b/tests/testthat/_snaps/plot_nma_data/duplicated-agd-arms-network.svg @@ -35,7 +35,7 @@ 1 Number of studies -3 +1 Duplicated AgD arms network diff --git a/tests/testthat/_snaps/plot_nma_data/duplicated-agd-contrasts-network.svg b/tests/testthat/_snaps/plot_nma_data/duplicated-agd-contrasts-network.svg index 5c8083a9..7f0bed07 100644 --- a/tests/testthat/_snaps/plot_nma_data/duplicated-agd-contrasts-network.svg +++ b/tests/testthat/_snaps/plot_nma_data/duplicated-agd-contrasts-network.svg @@ -35,7 +35,7 @@ 1 Number of studies -3 +1 Duplicated AgD contrasts network diff --git a/tests/testthat/_snaps/plot_nma_data/hta-psoriasis-network.svg b/tests/testthat/_snaps/plot_nma_data/hta-psoriasis-network.svg index 223a1a89..4b254f81 100644 --- a/tests/testthat/_snaps/plot_nma_data/hta-psoriasis-network.svg +++ b/tests/testthat/_snaps/plot_nma_data/hta-psoriasis-network.svg @@ -28,13 +28,13 @@ - + - + From a06a2a69567bd0ca90eb965c61b001e4e6ce779f Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 24 Apr 2024 12:05:04 +0100 Subject: [PATCH 35/57] Avoid duplicated override.aes warning in plaque psoriasis example --- vignettes/example_plaque_psoriasis.Rmd | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/vignettes/example_plaque_psoriasis.Rmd b/vignettes/example_plaque_psoriasis.Rmd index 49893136..3efa73a5 100644 --- a/vignettes/example_plaque_psoriasis.Rmd +++ b/vignettes/example_plaque_psoriasis.Rmd @@ -517,12 +517,12 @@ class_pal <- c("#D95F02", "#7570B3", "#E7298A", "#E6AB02") plot(pso_net, weight_nodes = TRUE, weight_edges = TRUE, show_trt_class = TRUE, nudge = 0.1) + ggraph::scale_edge_colour_manual("Data", - values = c(AgD = "#113259", IPD = "#55A480"), - guide = guide_legend(override.aes = list(edge_width = 2))) + + values = c(AgD = "#113259", IPD = "#55A480")) + scale_fill_manual("Treatment class", values = class_pal, - aesthetics = c("fill", "colour"), - guide = guide_legend(override.aes = list(size = 2))) + aesthetics = c("fill", "colour")) + + guides(edge_colour = guide_legend(override.aes = list(edge_width = 2)), + fill = guide_legend(override.aes = list(size = 2))) ``` ### Numerical integration for ML-NMR models From 142927c1722b1932041487323bdcd93ef0c3f679 Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 24 Apr 2024 12:07:07 +0100 Subject: [PATCH 36/57] Rebuild vignettes [ci build-vignettes] From 750a368af7a8f1171c8c190bab71d4b1d3c540fe Mon Sep 17 00:00:00 2001 From: David Phillippo Date: Wed, 24 Apr 2024 13:36:16 +0100 Subject: [PATCH 37/57] Fix actions commits [ci build-vignettes] --- .github/workflows/commit-commands.yaml | 8 ++++---- .github/workflows/pr-commands.yaml | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/.github/workflows/commit-commands.yaml b/.github/workflows/commit-commands.yaml index ac249718..97033567 100644 --- a/.github/workflows/commit-commands.yaml +++ b/.github/workflows/commit-commands.yaml @@ -31,8 +31,8 @@ jobs: - name: Commit run: | - git config --local user.name "$GITHUB_ACTOR" - git config --local user.email "$GITHUB_ACTOR@users.noreply.github.com" + git config user.email "41898282+github-actions[bot]@users.noreply.github.com" + git config user.name "github-actions[bot]" git add man/\* NAMESPACE git commit -m 'Document' git push @@ -64,8 +64,8 @@ jobs: - name: Commit run: | - git config --local user.name "$GITHUB_ACTOR" - git config --local user.email "$GITHUB_ACTOR@users.noreply.github.com" + git config user.email "41898282+github-actions[bot]@users.noreply.github.com" + git config user.name "github-actions[bot]" git add vignettes/\* tests/testthat/\*.R git commit -m 'Build pre-compiled vignettes' git push diff --git a/.github/workflows/pr-commands.yaml b/.github/workflows/pr-commands.yaml index bae4be1b..14fa3053 100644 --- a/.github/workflows/pr-commands.yaml +++ b/.github/workflows/pr-commands.yaml @@ -35,8 +35,8 @@ jobs: - name: Commit run: | - git config --local user.name "$GITHUB_ACTOR" - git config --local user.email "$GITHUB_ACTOR@users.noreply.github.com" + git config user.email "41898282+github-actions[bot]@users.noreply.github.com" + git config user.name "github-actions[bot]" git add man/\* NAMESPACE git commit -m 'Document' @@ -75,8 +75,8 @@ jobs: - name: Commit run: | - git config --local user.name "$GITHUB_ACTOR" - git config --local user.email "$GITHUB_ACTOR@users.noreply.github.com" + git config user.email "41898282+github-actions[bot]@users.noreply.github.com" + git config user.name "github-actions[bot]" git add vignettes/\* tests/testthat/\*.R git commit -m 'Build pre-compiled vignettes' From 5062e935ccaaeb6d9bf5cdf23721f4510c32cb06 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Wed, 24 Apr 2024 13:30:14 +0000 Subject: [PATCH 38/57] Build pre-compiled vignettes --- tests/testthat/test-example_bcg_vaccine.R | 60 +- tests/testthat/test-example_dietary_fat.R | 32 +- .../testthat/test-example_plaque_psoriasis.R | 8 +- tests/testthat/test-example_smoking.R | 8 +- vignettes/example_atrial_fibrillation.html | 2102 +++++----- vignettes/example_bcg_vaccine.html | 640 ++- vignettes/example_blocker.html | 313 +- vignettes/example_diabetes.html | 935 ++--- vignettes/example_dietary_fat.html | 394 +- vignettes/example_hta_psoriasis.html | 673 +-- vignettes/example_ndmm.html | 1269 +++--- vignettes/example_parkinsons.html | 1408 +++---- vignettes/example_plaque_psoriasis.html | 3604 +++++++++-------- vignettes/example_smoking.html | 249 +- vignettes/example_statins.html | 333 +- vignettes/example_thrombolytics.html | 1018 +++-- vignettes/example_transfusion.html | 232 +- 17 files changed, 6642 insertions(+), 6636 deletions(-) diff --git a/tests/testthat/test-example_bcg_vaccine.R b/tests/testthat/test-example_bcg_vaccine.R index ef7c091b..608ddc9c 100644 --- a/tests/testthat/test-example_bcg_vaccine.R +++ b/tests/testthat/test-example_bcg_vaccine.R @@ -8,17 +8,17 @@ skip_on_cran() params <- list(run_tests = FALSE) -## ---- code=readLines("children/knitr_setup.R"), include=FALSE----------------- +## ----code=readLines("children/knitr_setup.R"), include=FALSE-------------------------------------- -## ---- include=FALSE----------------------------------------------------------- +## ----include=FALSE-------------------------------------------------------------------------------- set.seed(18284729) -## ---- eval = FALSE------------------------------------------------------------ +## ----eval = FALSE--------------------------------------------------------------------------------- ## library(multinma) ## options(mc.cores = parallel::detectCores()) -## ----setup, echo = FALSE------------------------------------------------------ +## ----setup, echo = FALSE-------------------------------------------------------------------------- library(multinma) nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")), "true" =, "warn" = 2, @@ -26,11 +26,11 @@ nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")), options(mc.cores = nc) -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- head(bcg_vaccine) -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- bcg_net <- set_agd_arm(bcg_vaccine, study = studyn, trt = trtc, @@ -40,19 +40,19 @@ bcg_net <- set_agd_arm(bcg_vaccine, bcg_net -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- summary(normal(scale = 100)) summary(half_normal(scale = 5)) -## ---- eval = FALSE------------------------------------------------------------ +## ----eval = FALSE--------------------------------------------------------------------------------- ## bcg_fit_unadj <- nma(bcg_net, ## trt_effects = "random", ## prior_intercept = normal(scale = 100), ## prior_trt = normal(scale = 100), ## prior_het = half_normal(scale = 5)) -## ---- echo = FALSE------------------------------------------------------------ +## ----echo = FALSE--------------------------------------------------------------------------------- bcg_fit_unadj <- nowarn_on_ci( nma(bcg_net, seed = 14308133, @@ -64,25 +64,25 @@ bcg_fit_unadj <- nowarn_on_ci( ) -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- bcg_fit_unadj -## ---- eval=FALSE-------------------------------------------------------------- +## ----eval=FALSE----------------------------------------------------------------------------------- ## # Not run ## print(bcg_fit_unadj, pars = c("d", "mu", "delta", "tau")) -## ----bcg_unadj_pp_plot-------------------------------------------------------- +## ----bcg_unadj_pp_plot---------------------------------------------------------------------------- plot_prior_posterior(bcg_fit_unadj, prior = c("trt", "het")) -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- summary(normal(scale = 100)) summary(half_normal(scale = 5)) -## ---- eval = FALSE------------------------------------------------------------ +## ----eval = FALSE--------------------------------------------------------------------------------- ## bcg_fit_lat <- nma(bcg_net, ## trt_effects = "random", ## regression = ~.trt:latitude, @@ -92,7 +92,7 @@ summary(half_normal(scale = 5)) ## prior_het = half_normal(scale = 5), ## adapt_delta = 0.99) -## ---- echo = FALSE------------------------------------------------------------ +## ----echo = FALSE--------------------------------------------------------------------------------- bcg_fit_lat <- nowarn_on_ci( nma(bcg_net, seed = 1932599147, @@ -107,32 +107,32 @@ bcg_fit_lat <- nowarn_on_ci( ) -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- bcg_fit_lat -## ---- eval=FALSE-------------------------------------------------------------- +## ----eval=FALSE----------------------------------------------------------------------------------- ## # Not run ## print(bcg_fit_lat, pars = c("d", "beta", "mu", "delta", "tau")) -## ----bcg_lat_pp_plot---------------------------------------------------------- +## ----bcg_lat_pp_plot------------------------------------------------------------------------------ plot_prior_posterior(bcg_fit_lat, prior = c("trt", "reg", "het")) -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- (bcg_dic_unadj <- dic(bcg_fit_unadj)) -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- (bcg_dic_lat <- dic(bcg_fit_lat)) -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- summary(bcg_fit_unadj, pars = "tau") summary(bcg_fit_lat, pars = "tau") -## ----bcg_vaccine_beta_lat, fig.height = 4------------------------------------- +## ----bcg_vaccine_beta_lat, fig.height = 4--------------------------------------------------------- summary(bcg_fit_lat, pars = "beta") plot(bcg_fit_lat, @@ -141,7 +141,7 @@ plot(bcg_fit_lat, stat = "halfeye") -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- bcg_releff_lat <- relative_effects(bcg_fit_lat, newdata = tibble::tibble(latitude = seq(10, 50, by = 10), label = paste0(latitude, "\u00B0 latitude")), @@ -150,12 +150,12 @@ bcg_releff_lat <- relative_effects(bcg_fit_lat, bcg_releff_lat -## ----bcg_vaccine_releff_lat, fig.height = 5----------------------------------- +## ----bcg_vaccine_releff_lat, fig.height = 5------------------------------------------------------- plot(bcg_releff_lat, ref_line = 0) -## ----bcg_vaccine_reg_plot----------------------------------------------------- +## ----bcg_vaccine_reg_plot------------------------------------------------------------------------- library(dplyr) library(ggplot2) @@ -189,15 +189,15 @@ ggplot(aes(x = latitude), data = bcg_lor) + theme_multinma() -## ----bcg_vaccine_predictive_unadj--------------------------------------------- +## ----bcg_vaccine_predictive_unadj----------------------------------------------------------------- (bcg_predeff_unadj <- relative_effects(bcg_fit_unadj, predictive_distribution = TRUE)) -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- mean(as.matrix(bcg_predeff_unadj) > 0) -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- bcg_predeff_lat <- relative_effects(bcg_fit_lat, newdata = tibble::tibble(latitude = seq(0, 50, by = 10), label = paste0(latitude, "\u00B0 latitude")), @@ -207,11 +207,11 @@ bcg_predeff_lat <- relative_effects(bcg_fit_lat, bcg_predeff_lat -## ----------------------------------------------------------------------------- +## ------------------------------------------------------------------------------------------------- colMeans(as.matrix(bcg_predeff_lat) > 0) -## ----bcg_vaccine_tests, include=FALSE, eval=params$run_tests------------------ +## ----bcg_vaccine_tests, include=FALSE, eval=params$run_tests-------------------------------------- #--- Test against TSD 3 results --- library(testthat) library(dplyr) diff --git a/tests/testthat/test-example_dietary_fat.R b/tests/testthat/test-example_dietary_fat.R index 32fa64a5..7abb5a6a 100644 --- a/tests/testthat/test-example_dietary_fat.R +++ b/tests/testthat/test-example_dietary_fat.R @@ -8,17 +8,17 @@ skip_on_cran() params <- list(run_tests = FALSE) -## ---- code=readLines("children/knitr_setup.R"), include=FALSE------------------------------------- +## ----code=readLines("children/knitr_setup.R"), include=FALSE-------------------------------------- -## ---- eval = FALSE-------------------------------------------------------------------------------- +## ----eval = FALSE--------------------------------------------------------------------------------- ## library(multinma) ## options(mc.cores = parallel::detectCores()) ## ----setup, echo = FALSE-------------------------------------------------------------------------- library(multinma) -nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")), - "true" =, "warn" = 2, +nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")), + "true" =, "warn" = 2, parallel::detectCores()) options(mc.cores = nc) @@ -28,10 +28,10 @@ head(dietary_fat) ## ------------------------------------------------------------------------------------------------- -diet_net <- set_agd_arm(dietary_fat, +diet_net <- set_agd_arm(dietary_fat, study = studyc, trt = trtc, - r = r, + r = r, E = E, trt_ref = "Control", sample_size = n) @@ -43,7 +43,7 @@ summary(normal(scale = 100)) ## ------------------------------------------------------------------------------------------------- -diet_fit_FE <- nma(diet_net, +diet_fit_FE <- nma(diet_net, trt_effects = "fixed", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100)) @@ -53,7 +53,7 @@ diet_fit_FE <- nma(diet_net, diet_fit_FE -## ---- eval=FALSE---------------------------------------------------------------------------------- +## ----eval=FALSE----------------------------------------------------------------------------------- ## # Not run ## print(diet_fit_FE, pars = c("d", "mu")) @@ -67,16 +67,16 @@ summary(normal(scale = 100)) summary(half_normal(scale = 5)) -## ---- eval=FALSE---------------------------------------------------------------------------------- +## ----eval=FALSE----------------------------------------------------------------------------------- ## diet_fit_RE <- nma(diet_net, ## trt_effects = "random", ## prior_intercept = normal(scale = 100), ## prior_trt = normal(scale = 100), ## prior_het = half_normal(scale = 5)) -## ---- echo=FALSE, warning=FALSE------------------------------------------------------------------- +## ----echo=FALSE, warning=FALSE-------------------------------------------------------------------- diet_fit_RE <- nowarn_on_ci( - nma(diet_net, + nma(diet_net, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), @@ -88,7 +88,7 @@ diet_fit_RE <- nowarn_on_ci( diet_fit_RE -## ---- eval=FALSE---------------------------------------------------------------------------------- +## ----eval=FALSE----------------------------------------------------------------------------------- ## # Not run ## print(diet_fit_RE, pars = c("d", "mu", "delta")) @@ -113,15 +113,15 @@ plot(dic_RE) ## ----diet_pred_FE, fig.height = 2----------------------------------------------------------------- -pred_FE <- predict(diet_fit_FE, - baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5), +pred_FE <- predict(diet_fit_FE, + baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5), type = "response") pred_FE plot(pred_FE) ## ----diet_pred_RE, fig.height = 2----------------------------------------------------------------- -pred_RE <- predict(diet_fit_RE, - baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5), +pred_RE <- predict(diet_fit_RE, + baseline = distr(qnorm, mean = -3, sd = 1.77^-0.5), type = "response") pred_RE plot(pred_RE) diff --git a/tests/testthat/test-example_plaque_psoriasis.R b/tests/testthat/test-example_plaque_psoriasis.R index eb9b81eb..0cab1659 100644 --- a/tests/testthat/test-example_plaque_psoriasis.R +++ b/tests/testthat/test-example_plaque_psoriasis.R @@ -383,12 +383,12 @@ plot(pso_pred_FE_new, ref_line = c(0, 1)) ## ## plot(pso_net, weight_nodes = TRUE, weight_edges = TRUE, show_trt_class = TRUE, nudge = 0.1) + ## ggraph::scale_edge_colour_manual("Data", -## values = c(AgD = "#113259", IPD = "#55A480"), -## guide = guide_legend(override.aes = list(edge_width = 2))) + +## values = c(AgD = "#113259", IPD = "#55A480")) + ## scale_fill_manual("Treatment class", ## values = class_pal, -## aesthetics = c("fill", "colour"), -## guide = guide_legend(override.aes = list(size = 2))) +## aesthetics = c("fill", "colour")) + +## guides(edge_colour = guide_legend(override.aes = list(edge_width = 2)), +## fill = guide_legend(override.aes = list(size = 2))) ## ----eval=!params$run_tests----------------------------------------------------------------------- diff --git a/tests/testthat/test-example_smoking.R b/tests/testthat/test-example_smoking.R index 74ce8109..0e00b48b 100644 --- a/tests/testthat/test-example_smoking.R +++ b/tests/testthat/test-example_smoking.R @@ -8,10 +8,10 @@ skip_on_cran() params <- list(run_tests = FALSE) -## ---- code=readLines("children/knitr_setup.R"), include=FALSE------------------------------------- +## ----code=readLines("children/knitr_setup.R"), include=FALSE-------------------------------------- -## ---- eval = FALSE-------------------------------------------------------------------------------- +## ----eval = FALSE--------------------------------------------------------------------------------- ## library(multinma) ## options(mc.cores = parallel::detectCores()) @@ -37,7 +37,7 @@ smknet <- set_agd_arm(smoking, smknet -## ---- eval=FALSE---------------------------------------------------------------------------------- +## ----eval=FALSE----------------------------------------------------------------------------------- ## plot(smknet, weight_edges = TRUE, weight_nodes = TRUE) ## ----smoking_network_plot, echo=FALSE, fig.width=8, fig.height=6, out.width="100%"---------------- @@ -63,7 +63,7 @@ smkfit <- nma(smknet, smkfit -## ---- eval=FALSE---------------------------------------------------------------------------------- +## ----eval=FALSE----------------------------------------------------------------------------------- ## # Not run ## print(smkfit, pars = c("d", "tau", "mu", "delta")) diff --git a/vignettes/example_atrial_fibrillation.html b/vignettes/example_atrial_fibrillation.html index 4f5fbb0c..e59a0c44 100644 --- a/vignettes/example_atrial_fibrillation.html +++ b/vignettes/example_atrial_fibrillation.html @@ -49,7 +49,7 @@