In dmphillippo/multinma: Bayesian Network Meta-Analysis of Individual and Aggregate Data

set.seed(18284729)

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 13 trials investigating BCG vaccination vs. no vaccination for prevention of Tuberculosis (TB) [@TSD3;@Berkey1995]. The data are available in this package as bcg_vaccine:

head(bcg_vaccine)

@TSD3 used these data to demonstrate meta-regression models adjusting for the continuous covariate latitude, the absolute degrees latitude at which the study was conducted, which we recreate here.

Setting up the network

We have data giving the number diagnosed with TB during trial follow-up (r) out of the total (n) in each arm, so we use the function set_agd_arm() to set up the network. We set "unvaccinated" as the network reference treatment.

bcg_net <- set_agd_arm(bcg_vaccine, 
                       study = studyn,
                       trt = trtc,
                       r = r, 
                       n = n,
                       trt_ref = "Unvaccinated")
bcg_net

The latitude variable in the bcg_vaccine data frame will automatically be available to use in a meta-regression model.

Meta-analysis models

We fit random effects (RE) models, firstly without any covariates, and then with a meta-regression on the continuous covariate latitude.

RE meta-analysis (no covariate)

We start by fitting a standard RE model without any covariates. We use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effect $d_\mathrm{Vaccine}$ and study-specific intercepts $\mu_j$, and a $\textrm{half-N}(0, 5^2)$ prior distribution 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))

The model is fitted with the nma() function, with a random effects model specified by trt_effects = "random".

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))

bcg_fit_unadj <- nowarn_on_ci(
  nma(bcg_net, 
      seed = 14308133,
      iter = 10000,
      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:

bcg_fit_unadj

By default, summaries of the study-specific intercepts $\mu_j$ and random effects $\delta_j$ are hidden, but could be examined by changing the pars argument:

# Not run
print(bcg_fit_unadj, pars = c("d", "mu", "delta", "tau"))

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

plot_prior_posterior(bcg_fit_unadj, prior = c("trt", "het"))

RE meta-regression with covariate latitude

We now fit a RE meta-regression model, adjusting for latitude. We use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effect $d_\mathrm{Vaccine}$, study-specific intercepts $\mu_j$, and regression coefficient $\beta$. We use a $\text{half-N}(0, 5^2)$ prior distribution 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))

Again, the model is fitted with the nma() function. The regression formula ~ .trt:latitude means that the interaction of latitude with treatment will be included; the .trt special variable indicates treatment, and latitude is in the original data set. We increase adapt_delta to 0.99 to remove a small number of divergent transition errors (the default for RE models is set to 0.95).

bcg_fit_lat <- nma(bcg_net, 
                   trt_effects = "random",
                   regression = ~.trt:latitude,
                   prior_intercept = normal(scale = 100),
                   prior_trt = normal(scale = 100),
                   prior_reg = normal(scale = 100),
                   prior_het = half_normal(scale = 5),
                   adapt_delta = 0.99)

bcg_fit_lat <- nowarn_on_ci(
                 nma(bcg_net, 
                     seed = 1932599147,
                     iter = 5000,
                     trt_effects = "random",
                     regression = ~.trt:latitude,
                     prior_intercept = normal(scale = 100),
                     prior_trt = normal(scale = 100),
                     prior_reg = normal(scale = 100),
                     prior_het = half_normal(scale = 5),
                     adapt_delta = 0.99)
                 )

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

bcg_fit_lat

Note that latitude has automatically been centered at 33.46, the mean value for the studies in the network.

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(bcg_fit_lat, pars = c("d", "beta", "mu", "delta", "tau"))

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

plot_prior_posterior(bcg_fit_lat, prior = c("trt", "reg", "het"))

Model fit and comparison

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

(bcg_dic_unadj <- dic(bcg_fit_unadj))

(bcg_dic_lat <- dic(bcg_fit_lat))

The DIC is very similar between the two models, so we might at first choose the unadjusted model. The posterior mean residual deviance is larger for the model with the covariate, but this model also has a lower effective number of parameters $p_D$ so is allowing for more shrinkage of the random treatment effects. Moreover, the model with the covariate has a much lower estimated heterogeneity standard deviation:

summary(bcg_fit_unadj, pars = "tau")
summary(bcg_fit_lat, pars = "tau")

Adjusting for latitude is explaining a substantial amount of heterogeneity in the data. The 95% Credible Interval for the regression coefficient also excludes zero:

summary(bcg_fit_lat, pars = "beta")

plot(bcg_fit_lat, 
     pars = "beta", 
     ref_line = 0,
     stat = "halfeye")

Altogether, we might prefer the model with the adjustment for latitude. When considering covariates in random effects models it is important not to just look at the DIC [@TSD3]. We should also consider any reductions in heterogeneity, and the estimated regression coefficients and their standard error. The DIC is not sensitive to changes in the heterogeneity, as RE models are very flexible and can fit the data well whatever the level of heterogeneity.

Further results

We can produce estimates of the relative effect of vaccination at any latitude using the relative_effects() function. The newdata argument specifies a data frame containing the values of the covariate latitude that we are interested in, and the study argument is used to specify a column of newdata for an informative label.

bcg_releff_lat <- relative_effects(bcg_fit_lat,
                                   newdata = tibble::tibble(latitude = seq(10, 50, by = 10),
                                                            label = paste0(latitude, "\u00B0 latitude")),
                                   study = label)

bcg_releff_lat

The plot() method may be used to visually compare these estimates:

plot(bcg_releff_lat, 
     ref_line = 0)

A more sophisticated plot shows the regression line and confidence band for the effect of latitude, overlaid on the observed log odds ratios in each study:

library(dplyr)
library(ggplot2)

# Get data for regression line
lat_range <- range(bcg_vaccine$latitude)
lat_dat <- tibble(latitude = seq(lat_range[1], lat_range[2], by = 1))

bcg_lat_reg <- relative_effects(bcg_fit_lat, 
                                newdata = lat_dat) %>% 
  as_tibble() %>% 
  bind_cols(lat_dat)

# Get study log odds ratios
bcg_lor <- bcg_vaccine %>% 
  group_by(studyn) %>% 
  mutate(lor = log(r / (n - r)) - log(first(r) / (first(n) - first(r))),
         sample_size = sum(n)) %>% 
  slice(-1)

# Plot
ggplot(aes(x = latitude), data = bcg_lor) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), data = bcg_lat_reg,
              fill = "darkred", alpha = 0.3) +
  geom_line(aes(y = mean), data = bcg_lat_reg,
            colour = "darkred") +
  geom_point(aes(y = lor, size = sample_size), alpha = 0.6) +
  coord_cartesian(xlim = c(0, 60)) +
  xlab("Degrees Latitude") + ylab("log Odds Ratio") +
  scale_size("Sample Size") +
  theme_multinma()

In the presence of heterogeneity, it has been argued that decision makers should consider the predictive distribution of relative effects in a new study, instead of the posterior distribution of the mean treatment effects, as this reflects uncertainty due to heterogeneity and may better represent uncertainty about a future roll-out of a treatment [see @TSD3]. We can produce predictive distributions using the predictive_distribution = TRUE argument to relative_effects(). Dias et al. [-@Dias2018, section 8.3.2] consider the predictive distributions for this BCG vaccine analysis.

In the unadjusted analysis, whilst there is substantial evidence that vaccination is effective on average and essentially zero probability of harm based on the mean effect, the predictive distribution for effectiveness in a new study is wide and covers a range of harmful effects:

(bcg_predeff_unadj <- relative_effects(bcg_fit_unadj, predictive_distribution = TRUE))

The predictive probability of a new trial showing a harmful effect is:

mean(as.matrix(bcg_predeff_unadj) > 0)

For the analysis adjusting for latitude, the predictive distribution of relative effects now depends on latitude; here we calculate these in increments of 10 degrees from the equator:

bcg_predeff_lat <- relative_effects(bcg_fit_lat,
                                   newdata = tibble::tibble(latitude = seq(0, 50, by = 10),
                                                            label = paste0(latitude, "\u00B0 latitude")),
                                   study = label,
                                   predictive_distribution = TRUE)

bcg_predeff_lat

The predictive probabilities of a new trial carried out at a given latitude showing a harmful effect can be calculated as:

colMeans(as.matrix(bcg_predeff_lat) > 0)

So the predictive probability that a new trial carried out at the equator shows a harmful effect is around 80%, whereas at 50 degrees latitude the predictive probability is only 0.7%.

References

#--- Test against TSD 3 results ---
library(testthat)
library(dplyr)

tol <- 0.05
tol_dic <- 0.1

# Relative effects
bcg_unadj_releff <- as.data.frame(summary(bcg_fit_unadj, pars = "d"))

test_that("Unadjusted relative effects", {
  expect_equivalent(bcg_unadj_releff$mean, -0.762, tolerance = tol)
  expect_equivalent(bcg_unadj_releff$sd, 0.22, tolerance = tol)
  expect_equivalent(bcg_unadj_releff$`2.5%`, -1.21, tolerance = tol)
  expect_equivalent(bcg_unadj_releff$`97.5%`, -0.34, tolerance = tol)
})

test_that("Unadjusted predictive distribution", {
  bcg_unadj_releff_pred <- as.data.frame(relative_effects(bcg_fit_unadj, predictive_distribution = TRUE))
  expect_equivalent(bcg_unadj_releff_pred$mean, -0.762, tolerance = tol)
  expect_equivalent(bcg_unadj_releff_pred$`2.5%`, -2.27, tolerance = tol)
  expect_equivalent(bcg_unadj_releff_pred$`97.5%`, 0.72, tolerance = tol)
})

bcg_lat_releff <- as.data.frame(summary(bcg_fit_lat, pars = "d"))

test_that("Regression relative effects", {
  expect_equivalent(bcg_lat_releff$mean, -0.763, tolerance = tol)
  expect_equivalent(bcg_lat_releff$sd, 0.126, tolerance = tol)
  expect_equivalent(bcg_lat_releff$`2.5%`, -1.04, tolerance = tol)
  expect_equivalent(bcg_lat_releff$`97.5%`, -0.52, tolerance = tol)
})


test_that("Regression predictive distribution", {
  bcg_lat_releff_pred <- relative_effects(bcg_fit_lat,
                                          newdata = data.frame(latitude = c(0, 13, 50)),
                                          predictive_distribution = TRUE)
  expect_equivalent(colMeans(as.matrix(bcg_lat_releff_pred) > 0), c(0.8, 0.35, 0.006), tolerance = tol)
})

# Regression coefficients
bcg_lat_beta <- as.data.frame(summary(bcg_fit_lat, pars = "beta"))

test_that("Regression beta", {
  expect_equivalent(bcg_lat_beta$mean, -0.032, tolerance = tol)
  expect_equivalent(bcg_lat_beta$sd, 0.009, tolerance = tol)
  expect_equivalent(bcg_lat_beta$`2.5%`, -0.05, tolerance = tol)
  expect_equivalent(bcg_lat_beta$`97.5%`, -0.01, tolerance = tol)
})

# RE heterogeneity SD
bcg_unadj_sd <- as.data.frame(summary(bcg_fit_unadj, pars = "tau"))

test_that("Unadjusted heterogeneity SD", {
  expect_equivalent(bcg_unadj_sd$`50%`, 0.649, tolerance = tol)
  expect_equivalent(bcg_unadj_sd$sd, 0.202, tolerance = tol)
  expect_equivalent(bcg_unadj_sd$`2.5%`, 0.39, tolerance = tol)
  expect_equivalent(bcg_unadj_sd$`97.5%`, 1.17, tolerance = tol)
})

bcg_lat_sd <- as.data.frame(summary(bcg_fit_lat, pars = "tau"))

test_that("Regression heterogeneity SD", {
  expect_equivalent(bcg_lat_sd$`50%`, 0.272, tolerance = tol)
  expect_equivalent(bcg_lat_sd$sd, 0.188, tolerance = tol)
  expect_equivalent(bcg_lat_sd$`2.5%`, 0.03, tolerance = tol)
  expect_equivalent(bcg_lat_sd$`97.5%`, 0.75, tolerance = tol)
})

# DIC
test_that("Unadjusted DIC", {
  expect_equivalent(bcg_dic_unadj$resdev, 26.1, tolerance = tol_dic)
  expect_equivalent(bcg_dic_unadj$pd, 23.5, tolerance = tol_dic)
  expect_equivalent(bcg_dic_unadj$dic, 49.6, tolerance = tol_dic)
})

test_that("Regression DIC", {
  expect_equivalent(bcg_dic_lat$resdev, 30.4, tolerance = tol_dic)
  expect_equivalent(bcg_dic_lat$pd, 21.1, tolerance = tol_dic)
  expect_equivalent(bcg_dic_lat$dic, 51.5, tolerance = tol_dic)
})

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
  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", {
  pred_all <- predict(bcg_fit_lat, type = "response")
  pred_12 <- predict(bcg_fit_lat, type = "response", 
                     baseline = list("1" = "1", "2" = "2"),
                     newdata = subset(bcg_vaccine, studyn %in% 1:2),
                     study = studyn)

  expect_equal(unclass(as.array(pred_12)),
               as.array(pred_all)[ , , 1:4])
})

dmphillippo/multinma documentation built on April 12, 2025, 11:41 a.m.