set.seed(2684319)
library(multinma) options(mc.cores = parallel::detectCores())
library(multinma) nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")), "true" =, "warn" = 2, parallel::detectCores()) options(mc.cores = nc)
This vignette describes the analysis of 6 trials comparing transfusion of granulocytes (white blood cells) to control for preventing mortality in patients with neutropenia or neutrophil dysfunction [@Stanworth2005; @Turner2012].
The data are available in this package as transfusion:
head(transfusion)
@Turner2012 previously used this dataset to demonstrate the application of informative priors for heterogeneity, an analysis which we recreate here.
Setting up the network
We begin by setting up the network - here just a pairwise meta-analysis.
We have arm-level count data giving the number of deaths (r) out of the total (n) in each arm, so we use the function set_agd_arm().
We set "Control" as the reference treatment.
tr_net <- set_agd_arm(transfusion, study = studyc, trt = trtc, r = r, n = n, trt_ref = "Control") tr_net
Meta-analysis models
We fit two random effects models, first with a non-informative prior for the heterogeneity, then using the informative prior described by @Turner2012.
Random effects meta-analysis with non-informative heterogeneity prior
We fit a random effects model using the nma() function with trt_effects = "random".
We use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$ and study-specific intercepts $\mu_j$, and a non-informative $\textrm{half-N}(5^2)$ prior for the 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)) summary(half_normal(scale = 5))
Fitting the RE model
tr_fit_RE_noninf <- nma(tr_net, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5))
tr_fit_RE_noninf <- nma(tr_net, seed = 857369814, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5))
tr_fit_RE_noninf <- suppressWarnings(nma(tr_net, seed = 857369814, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5), iter = 10000, save_warmup = FALSE))
Basic parameter summaries are given by the print() method:
tr_fit_RE_noninf
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(tr_fit_RE_noninf, pars = c("d", "mu", "delta"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(tr_fit_RE_noninf, prior = "het")
The posterior distribution for the heterogeneity variance $\tau^2$ is summarised by
noninf_tau <- as.array(tr_fit_RE_noninf, pars = "tau") noninf_tausq <- noninf_tau^2 names(noninf_tausq) <- "tausq" summary(noninf_tausq)
Random effects meta-analysis with informative heterogeneity prior
Keeping the rest of the model setup the same, we now use an informative $\textrm{log-N}(-3.93, 1.51^2)$ prior for the heterogeneity variance $\tau^2$.
We can examine the range of parameter values implied by this prior distribution with the summary() method:
summary(log_normal(-3.93, 1.51))
Fitting the RE model, we specify the log_normal prior distribution in the prior_het argument, and set prior_het_type = "var" to indicate that this prior distribution is on the variance scale (instead of the standard deviation, the default).
tr_fit_RE_inf <- nma(tr_net, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = log_normal(-3.93, 1.51), prior_het_type = "var")
tr_fit_RE_inf <- nma(tr_net, seed = 1803772660, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = log_normal(-3.93, 1.51), prior_het_type = "var")
tr_fit_RE_inf <- suppressWarnings(nma(tr_net, seed = 1803772660, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = log_normal(-3.93, 1.51), prior_het_type = "var", iter = 10000, save_warmup = FALSE))
Basic parameter summaries are given by the print() method:
tr_fit_RE_inf
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(tr_fit_RE_inf, pars = c("d", "mu", "delta"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(tr_fit_RE_inf, prior = "het")
Note: The heterogeneity variance $\tau^2$ is plotted here since the prior was specified on $\tau^2$.
The posterior distribution for the heterogeneity variance $\tau^2$ is summarised by
inf_tau <- as.array(tr_fit_RE_inf, pars = "tau") inf_tausq <- inf_tau^2 names(inf_tausq) <- "tausq" summary(inf_tausq)
References
#--- Test against TSD 2 results --- library(testthat) tol <- 0.05 # Non-informative prior tr_RE_noninf_var <- as.data.frame(summary(noninf_tausq)) test_that("Non-informative RE heterogeneity variance", { skip("Non-informative priors not identical") expect_equivalent(tr_RE_noninf_var$`50%`, 2.74, tolerance = tol) expect_equivalent(tr_RE_noninf_var$`2.5%`, 0.34, tolerance = tol) expect_equivalent(tr_RE_noninf_var$`97.5%`, 18.1, tolerance = tol) }) # Informative prior tr_RE_inf_var <- as.data.frame(summary(inf_tausq)) test_that("Informative RE heterogeneity variance", { expect_equivalent(tr_RE_inf_var$`50%`, 0.18, tolerance = tol) expect_equivalent(tr_RE_inf_var$`2.5%`, 0.003, tolerance = tol) skip_on_ci() expect_equivalent(tr_RE_inf_var$`97.5%`, 1.84, tolerance = tol) })
