In audrey-b/BUGSnet: Bayesian network meta-analyses in compliance with best practice and reporting guidelines

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(BUGSnet)

data(afib)
dataprep <- data.prep(arm.data = afib,
                      varname.t = "treatment",
                      varname.s = "study")

Feasibility Assessment

Network Plot

 net.plot(dataprep, node.scale = 1, 
          edge.scale = 1,
          label.offset1 = 2)

Generate Network Characteristics via net.tab()

network.char <- net.tab(data = dataprep,
                        outcome = "events",
                        N = "sampleSize",
                        type.outcome = "binomial")

Network Characteristics

network.char$network generates characteristics about the network, such as connectedness, number of treatments in the network, and in the case of a binomial outcome, the number of events in the network.

knitr::kable(network.char$network)

Intervention Characteristics

network.char$intervention generates outcome and sample size data broken down by treatment.

knitr::kable(network.char$intervention)

Comparison Characteristics

network.char$comparison generates outcome and sample size data broken down by treatment comparison.

knitr::kable(network.char$comparison)

Main analysis

nma.model() creates BUGS code and that will be put into nma.run() and analysed through JAGS [@JAGS]. The reference parameter indicates the name of the treatment that will be seen as the 'referent' comparator, this is often a placebo of some sort. In our case, it is treatment 02. Since our outcome is dichotomous, and we are not interested in event rates, we are using the "binomial" family. In our case, we want to compare odds ratios, so we are using the $logit$ link.

fixed_effects_model <- nma.model(data = dataprep,
                                  outcome = "events",
                                  N = "sampleSize",
                                  reference = "02",
                                  family = "binomial",
                                  link = "logit",
                                  effects = "fixed",
                                  covariate = "stroke",
                                  prior.beta = "EXCHANGEABLE")

random_effects_model <- nma.model(data = dataprep,
                                  outcome = "events",
                                  N = "sampleSize",
                                  reference = "02",
                                  family = "binomial",
                                  link = "logit",
                                  effects = "random",
                                  covariate = "stroke",
                                  prior.beta = "EXCHANGEABLE")

If you want to review the BUGS code, you can review it by outputting cat(random_effects_model$bugs).

The next step is to run the NMA model using nma.run(). Since we are working in a Bayesian framework, we need to specify the number of adaptations, burn-ins, and iterations. A description of Bayesian MCMC is omitted here, we direct the reader to any introductory text on Bayesian Modelling [@lunn2012bugs].

fixed_effects_results <- nma.run(fixed_effects_model,
                           n.adapt = 1000,
                           n.burnin = 1000,
                           n.iter = 5000)

random_effects_results <- nma.run(random_effects_model,
                           n.adapt = 1000,
                           n.burnin = 1000,
                           n.iter = 5000)

Model Choice

Compare fixed vs random effects models by comparing leverage plots and the DIC

par(mfrow = c(1,2))
nma.fit(fixed_effects_results, main = "Fixed Effects Model" )
nma.fit(random_effects_results, main = "Random Effects Model")

The random effects model seems to be more appropriate here due to a lower DIC, and slightly better leverage plot.

Results

Regression plot

We can see the effect of the covariate on the log odds ratio on the following plot

nma.regplot(random_effects_results)

The lines are not parallel because we specified EXCHANGEABLE priors on the regression coefficients.

Let's looks at the relative treatment differences and ranking of treaments when the proportion of patients with a prior stroke is equal to 0.1.

Sucra Plot

sucra.out <- nma.rank(random_effects_results, largerbetter = FALSE, cov.value = 0.1, sucra.palette = "Set1")
sucra.out$sucraplot

League Plot

league.out <- nma.league(random_effects_results, 
                             central.tdcy = "median",
                             order = as.vector(t(dataprep$treatments)),
                             cov.value = 0.1, 
                             log.scale = FALSE)
league.out$heatplot

Forest Plot

nma.forest(random_effects_results, 
           comparator = "02", 
           central.tdcy = "median",
           cov.value = 0.1,
           x.trans = "log")

References

audrey-b/BUGSnet documentation built on Feb. 2, 2025, 5:10 p.m.