In audrey-b/BUGSnet: Bayesian network meta-analyses in compliance with best practice and reporting guidelines
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
The purpose of this vignette is to show how BUGSnet can be used to perform NMA using contrast-based data. The data we will be using is parkinsons in BUGSnet [@TSD2]. Treatment 1 corresponds to a placebo and treatments 2-5 correspond to active drugs. This dataset was used in [@TSD2] to demonstrate contrast-based NMA, and our results are very similar.
Data Preparation
The format for the contrast-based summary data object that gets passed to data.prep is important. The object should have one row per arm. The contrast-based summary data goes in one column. In this dataset, this column is called differences. In this column, the rows corresponding to arm 1 for each study has NA. The standard errors of the contrast-based summary goes in another column. In this dataset, that column is called se.diffs. Similarly to the previous column, rows corresponding to arm 1 of a study should have NA. If the network includes multi-arm trials, we require an additional column used to specify the variance of the response in arm 1 for each trial, which has a variance in the rows corresponding to a first arm and NAs elsewhere. In this dataset, this column is called var.arm1. Additionally, there must be a column which specifies which study each arm corresponds to (in this example, study) and a column which specifies which treatment is used in each arm (in this example, treatment). The parkinsons data is shown below.
library(BUGSnet)
data(parkinsons)
knitr::kable(parkinsons)
The data is prepared for analysis in BUGSnet by running it through data.prep().
library(BUGSnet)
data(parkinsons)
contrast_prep <- data.prep(arm.data = parkinsons,
varname.t = "treatment",
varname.s = "study")
Network of Evidence
To generate the network plot, input the data.prep object into the net.plot() function.
net.plot(contrast_prep, node.scale = 1.5, edge.scale=0.5)
The network has 5 treatments and is somewhat connected.
Main Analysis
When contrast-based summary data is used, the Normal distribution is used with the identity link for all types of contrast data. Some examples of contrast-based data summaries are log odds ratios, mean differences, log hazard ratios, and risk differences. In the parkinsons dataset, the data are mean differences. A mathematical description of the model follows.
Let $i=1,...,M$ be labels for the $M$ studies in the data and let $t_{ik}\in{1,…,T}$ denote the treatment used in arm $k$ of study $i$. The set ${1,…,T}$ represents the set of treatments that were assessed across the $M$ studies. Let $y_{i,k}$ be the contrast (mean difference, log-odds, ratio, etc) of arm $k$ compared to the reference arm (arm 1) of study $i$. Then we assume $y_{ik}\sim \text{Normal}(\theta_{ik}, V_{ik})$. An identity link is used such that
$$\theta_{ik} = \delta_{ik}$$
where $\delta_{ik}$ represents the fixed or random effect of the treatment from arm $k$ compared to the treatment in arm 1. Then we also have $\delta_{i1} = 0$ for $i = 1,\dots,M$.
In the random effects model, the $\boldsymbol \delta_i’s=(\delta_{i2},...,\delta_{ia_i} )^\top$ are conditionally independent with distributions
$$[\boldsymbol\delta_i│\boldsymbol d_i,\Sigma] \sim MVNormal(\boldsymbol d_i,\Sigma),$$
where $\boldsymbol d_i=(d_{(t_{i1},t_{i2})},...,d_{(t_{i1},t_{ia_i)}})^\top$ and $d_{(t_{i1},t_{ik})}=d_{(1,t_{ik})}-d_{(1,t_{i1})}$ is the difference in the treatment effect of treatments $t_{i1}$ and $t_{ik}$ on the scale of the data and $d_{(1,1)}=0$. For $\Sigma$, a compound symmetry structure is specified following (16), with variances $\sigma^2$ and covariances $0.5\sigma^2$, where $\sigma^2$ represents the between-trial variability in treatment effects (heterogeneity).
In the fixed effect model, the $\delta_{ik}$’s are treated as fixed (from a frequentist perspective) and are defined as $\delta_{ik}=d_{(t_{i1},t_{ik})}=d_{(1,t_{ik})}-d_{(1,t_{i1})}$ with $d_{(1,1)}=0$.
Independent priors are specified on $d_{(1,2)},...,d_{(1,T)}$. Note that in both fixed and random-effects models, the quantities of interest are $d_{(1,2)},...,d_{(1,T)}$ which allow us to conclude on all treatment contrasts though the transitivity relation, $d_{(t_{i1},t_{ik})}=d_{(1,t_{ik})}-d_{(1,t_{i1})}$.
Returning to the analysis, nma.model.contrast() creates BUGS code for contrast-based data 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, and corresponds to the treatment labelled as 1. In our case, we will use treatment 1 (placebo). The scale parameter is used to specify what the scale of the contrasts is, for example, Log-odds ratio or Mean difference, and is used in subsequent plotting functions.
We will fit a fixed effects model and a random effects model.
fixed_effects_model <- nma.model.contrast(data=contrast_prep,
differences="differences",
se.diffs = "se.diffs",
var.arm1 = "var.arm1",
reference = "1",
effects = "fixed",
scale = "Mean differences")
random_effects_model <- nma.model.contrast(data=contrast_prep,
differences="differences",
se.diffs = "se.diffs",
var.arm1 = "var.arm1",
reference = "1",
effects = "random",
scale = "Mean differences")
If you want to review or modify the BUGS code, you can review it by outputting cat(fixed_effects_model$bugs) and cat(random_effects_model$bugs). Priors are defined on $d_{(1,2)},...,d_{(1,T)}$ in both fixed and random effect models, and an additional prior is defined on $\sigma$ in the random effect model. By default, BUGSnet implements the vague priors proposed in @gemtc. You may specify your own priors instead by specifying the options prior.d and prior.sigma in the nma.model.contrast() function.
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].
set.seed(2021)
fixed_effects_results <- nma.run(fixed_effects_model,
n.adapt=1000,
n.burnin=2000,
n.iter=12000)
random_effects_results <- nma.run(random_effects_model,
n.adapt=1000,
n.burnin=4000,
n.iter=15000)
Assess model fit
Compare the fit of the fixed and random effects models by comparing the leverage plots and 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 DIC is less for the fixed effect model. The leverage plot for the fixed effect model also looks a bit better, although the differences are small. Overall, we will continue with the fixed effects model.
Check Inconsistency
Next, we will assess consistency in the network by fitting a fixed effects inconsistency model. We can then compare the fit of this model to the consistency model. If the inconsistency model has a better fit, then it is likely that there is inconsistency in the network [@lu2006assessing].
fe_inconsistency_model <- nma.model.contrast(data=contrast_prep,
differences="differences",
se.diffs = "se.diffs",
var.arm1 = "var.arm1",
reference = "1",
effects = "fixed",
type = "inconsistency",
scale = "Mean differences")
fe_inconsistency_results <- nma.run(fe_inconsistency_model,
n.adapt=1000,
n.burnin=1000,
n.iter=10000)
We use rainbow plots and DIC to assess fit.
par(mfrow = c(1,2))
fe_model_fit <- nma.fit(fixed_effects_results, main = "Consistency Model" )
inconsist_model_fit <- nma.fit(fe_inconsistency_results, main= "Inconsistency Model")
Here, the consistency model has a lower DIC.
Furthermore, a plot of the posterior mean deviance of the individual data points in the inconsistency model against their posterior mean deviance in the consistency model can highlight descrepancies between the two models.
par(mfrow=c(1,1))
nma.compare(fe_model_fit, inconsist_model_fit)
The data lies near the line $y=x$, indicating agreement between the models. This suggests we can proceed with the more parsimonious model, the consistency model.
Results
Simulaneous comparison of every treatment based on the results of the NMA analysis can be achieved by comparing the posterior probabilities of being the best, second best,…, worst treatment. In BUGSnet, simply input the results of your model into the nma.rank() function, and specify largerbetter=TRUE if a larger outcome is associated with better treatments, and FALSE otherwise. In our case, a more negative outcome is associated with better treatments, so we set largerbetter=FALSE.
Sucra Plot
sucra.out <- nma.rank(fixed_effects_results, largerbetter=FALSE, sucra.palette= "Set1")
sucra.out$sucraplot
League Heat Plot
League tables are another way to summarize the results of an NMA. League tables contain all information about relative effectiveness for all possible pairs of interventions [@rouse2017network]. BUGSnet includes 95% credible intervals. You can also plot the league table as a heatplot using the following code:
league.out <- nma.league(fixed_effects_results,
central.tdcy="median",
order = sucra.out$order,
low.colour = "springgreen4",
mid.colour = "white",
high.colour = "red")
league.out$heatplot
Forest Plot
Forest plots are another great way to summarize the results of an NMA with respect to a particular comparator. The x-axis label is based on the scale specified when nma.model.contrast was run.
nma.forest(fixed_effects_results,
central.tdcy="median",
comparator = "1")
Appendix
The traceplots and densities for each chain can be used to assess convergence of the MCMC chains.
nma.diag(fixed_effects_results, plot_prompt = F)
nma.diag(random_effects_results, plot_prompt = F)
audrey-b/BUGSnet documentation built on Feb. 2, 2025, 5:10 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
The purpose of this vignette is to show how BUGSnet can be used to perform NMA using contrast-based data. The data we will be using is parkinsons in BUGSnet [@TSD2]. Treatment 1 corresponds to a placebo and treatments 2-5 correspond to active drugs. This dataset was used in [@TSD2] to demonstrate contrast-based NMA, and our results are very similar.
Data Preparation
The format for the contrast-based summary data object that gets passed to data.prep is important. The object should have one row per arm. The contrast-based summary data goes in one column. In this dataset, this column is called differences. In this column, the rows corresponding to arm 1 for each study has NA. The standard errors of the contrast-based summary goes in another column. In this dataset, that column is called se.diffs. Similarly to the previous column, rows corresponding to arm 1 of a study should have NA. If the network includes multi-arm trials, we require an additional column used to specify the variance of the response in arm 1 for each trial, which has a variance in the rows corresponding to a first arm and NAs elsewhere. In this dataset, this column is called var.arm1. Additionally, there must be a column which specifies which study each arm corresponds to (in this example, study) and a column which specifies which treatment is used in each arm (in this example, treatment). The parkinsons data is shown below.
library(BUGSnet) data(parkinsons) knitr::kable(parkinsons)
The data is prepared for analysis in BUGSnet by running it through data.prep().
library(BUGSnet) data(parkinsons) contrast_prep <- data.prep(arm.data = parkinsons, varname.t = "treatment", varname.s = "study")
Network of Evidence
To generate the network plot, input the data.prep object into the net.plot() function.
net.plot(contrast_prep, node.scale = 1.5, edge.scale=0.5)
The network has 5 treatments and is somewhat connected.
Main Analysis
When contrast-based summary data is used, the Normal distribution is used with the identity link for all types of contrast data. Some examples of contrast-based data summaries are log odds ratios, mean differences, log hazard ratios, and risk differences. In the parkinsons dataset, the data are mean differences. A mathematical description of the model follows.
Let $i=1,...,M$ be labels for the $M$ studies in the data and let $t_{ik}\in{1,…,T}$ denote the treatment used in arm $k$ of study $i$. The set ${1,…,T}$ represents the set of treatments that were assessed across the $M$ studies. Let $y_{i,k}$ be the contrast (mean difference, log-odds, ratio, etc) of arm $k$ compared to the reference arm (arm 1) of study $i$. Then we assume $y_{ik}\sim \text{Normal}(\theta_{ik}, V_{ik})$. An identity link is used such that
$$\theta_{ik} = \delta_{ik}$$
where $\delta_{ik}$ represents the fixed or random effect of the treatment from arm $k$ compared to the treatment in arm 1. Then we also have $\delta_{i1} = 0$ for $i = 1,\dots,M$. In the random effects model, the $\boldsymbol \delta_i’s=(\delta_{i2},...,\delta_{ia_i} )^\top$ are conditionally independent with distributions
$$[\boldsymbol\delta_i│\boldsymbol d_i,\Sigma] \sim MVNormal(\boldsymbol d_i,\Sigma),$$ where $\boldsymbol d_i=(d_{(t_{i1},t_{i2})},...,d_{(t_{i1},t_{ia_i)}})^\top$ and $d_{(t_{i1},t_{ik})}=d_{(1,t_{ik})}-d_{(1,t_{i1})}$ is the difference in the treatment effect of treatments $t_{i1}$ and $t_{ik}$ on the scale of the data and $d_{(1,1)}=0$. For $\Sigma$, a compound symmetry structure is specified following (16), with variances $\sigma^2$ and covariances $0.5\sigma^2$, where $\sigma^2$ represents the between-trial variability in treatment effects (heterogeneity).
In the fixed effect model, the $\delta_{ik}$’s are treated as fixed (from a frequentist perspective) and are defined as $\delta_{ik}=d_{(t_{i1},t_{ik})}=d_{(1,t_{ik})}-d_{(1,t_{i1})}$ with $d_{(1,1)}=0$. Independent priors are specified on $d_{(1,2)},...,d_{(1,T)}$. Note that in both fixed and random-effects models, the quantities of interest are $d_{(1,2)},...,d_{(1,T)}$ which allow us to conclude on all treatment contrasts though the transitivity relation, $d_{(t_{i1},t_{ik})}=d_{(1,t_{ik})}-d_{(1,t_{i1})}$.
Returning to the analysis, nma.model.contrast() creates BUGS code for contrast-based data 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, and corresponds to the treatment labelled as 1. In our case, we will use treatment 1 (placebo). The scale parameter is used to specify what the scale of the contrasts is, for example, Log-odds ratio or Mean difference, and is used in subsequent plotting functions.
We will fit a fixed effects model and a random effects model.
fixed_effects_model <- nma.model.contrast(data=contrast_prep, differences="differences", se.diffs = "se.diffs", var.arm1 = "var.arm1", reference = "1", effects = "fixed", scale = "Mean differences") random_effects_model <- nma.model.contrast(data=contrast_prep, differences="differences", se.diffs = "se.diffs", var.arm1 = "var.arm1", reference = "1", effects = "random", scale = "Mean differences")
If you want to review or modify the BUGS code, you can review it by outputting cat(fixed_effects_model$bugs) and cat(random_effects_model$bugs). Priors are defined on $d_{(1,2)},...,d_{(1,T)}$ in both fixed and random effect models, and an additional prior is defined on $\sigma$ in the random effect model. By default, BUGSnet implements the vague priors proposed in @gemtc. You may specify your own priors instead by specifying the options prior.d and prior.sigma in the nma.model.contrast() function.
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].
set.seed(2021) fixed_effects_results <- nma.run(fixed_effects_model, n.adapt=1000, n.burnin=2000, n.iter=12000) random_effects_results <- nma.run(random_effects_model, n.adapt=1000, n.burnin=4000, n.iter=15000)
Assess model fit
Compare the fit of the fixed and random effects models by comparing the leverage plots and 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 DIC is less for the fixed effect model. The leverage plot for the fixed effect model also looks a bit better, although the differences are small. Overall, we will continue with the fixed effects model.
Check Inconsistency
Next, we will assess consistency in the network by fitting a fixed effects inconsistency model. We can then compare the fit of this model to the consistency model. If the inconsistency model has a better fit, then it is likely that there is inconsistency in the network [@lu2006assessing].
fe_inconsistency_model <- nma.model.contrast(data=contrast_prep, differences="differences", se.diffs = "se.diffs", var.arm1 = "var.arm1", reference = "1", effects = "fixed", type = "inconsistency", scale = "Mean differences") fe_inconsistency_results <- nma.run(fe_inconsistency_model, n.adapt=1000, n.burnin=1000, n.iter=10000)
We use rainbow plots and DIC to assess fit.
par(mfrow = c(1,2)) fe_model_fit <- nma.fit(fixed_effects_results, main = "Consistency Model" ) inconsist_model_fit <- nma.fit(fe_inconsistency_results, main= "Inconsistency Model")
Here, the consistency model has a lower DIC.
Furthermore, a plot of the posterior mean deviance of the individual data points in the inconsistency model against their posterior mean deviance in the consistency model can highlight descrepancies between the two models.
par(mfrow=c(1,1)) nma.compare(fe_model_fit, inconsist_model_fit)
The data lies near the line $y=x$, indicating agreement between the models. This suggests we can proceed with the more parsimonious model, the consistency model.
Results
Simulaneous comparison of every treatment based on the results of the NMA analysis can be achieved by comparing the posterior probabilities of being the best, second best,…, worst treatment. In BUGSnet, simply input the results of your model into the nma.rank() function, and specify largerbetter=TRUE if a larger outcome is associated with better treatments, and FALSE otherwise. In our case, a more negative outcome is associated with better treatments, so we set largerbetter=FALSE.
Sucra Plot
sucra.out <- nma.rank(fixed_effects_results, largerbetter=FALSE, sucra.palette= "Set1") sucra.out$sucraplot
League Heat Plot
League tables are another way to summarize the results of an NMA. League tables contain all information about relative effectiveness for all possible pairs of interventions [@rouse2017network]. BUGSnet includes 95% credible intervals. You can also plot the league table as a heatplot using the following code:
league.out <- nma.league(fixed_effects_results, central.tdcy="median", order = sucra.out$order, low.colour = "springgreen4", mid.colour = "white", high.colour = "red") league.out$heatplot
Forest Plot
Forest plots are another great way to summarize the results of an NMA with respect to a particular comparator. The x-axis label is based on the scale specified when nma.model.contrast was run.
nma.forest(fixed_effects_results, central.tdcy="median", comparator = "1")
Appendix
The traceplots and densities for each chain can be used to assess convergence of the MCMC chains.
nma.diag(fixed_effects_results, plot_prompt = F) nma.diag(random_effects_results, plot_prompt = F)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.