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)
The purpose of this example is to demonstrate how to use BUGSnet to perform network meta-analysis on time-to-event data, where each study has a different followup time. The data was taken from the NICE DSU TSD2 [@TSD2], page 66. The results of our package are concordant with those reported in the TSD.
The Diabetes data has 22 studies and 6 treatments, each study has a follow up time between 1 and 6.1 years. The outcome for this data (diabetes.sim) is the number of new diabetes cases observed during the trial period. An age variable, including its standard deviation, was simulated for demonstration purposes. This data contains multi-arm trials, which can easily be accomodated using BUGSnet. The data is shown here:
data(diabetes.sim)
knitr::kable(diabetes.sim)
Before any analysis is done, we need to put our data in a form that BUGSnet can easily handle. This is done by simple running the data.prep() function and specifying the name of the treatment, study, and sample size variables:
library(BUGSnet)
data(diabetes.sim)
rate2.slr <- data.prep(arm.data = diabetes.sim,
varname.t = "Treatment",
varname.s = "Study")
We will use this data.prep object throughout the rest of our example.
Network of Evidence
To generate a network plot, simply input the data.prep object into the net.plot() function as shown:
par(mar=c(1,1,1,1)) #reduce margins for nicer plotting
net.plot(rate2.slr, node.scale = 3, edge.scale=1.5)
Our network has 6 treatments and a high degree of connectedness. If a user wants to highlight a specific node and its direct evidence, this can be done using the flag option.
par(mar=c(1,1,1,1)) #reduce margins for nicer plotting
net.plot(rate2.slr, node.scale = 3, edge.scale=1.5, flag = "ARB")
This gives us the names of the 5 trials connecting the ARB node to 4 treatments.
Generate Network Characteristics via net.tab()
In addition to a plot of the network, BUGSnet can summarize network, intervention, and comparison statistics that can be used to better understand our network of evidence.
network.char <- net.tab(data = rate2.slr,
outcome = "diabetes",
N = "n",
type.outcome = "rate2",
time = "followup")
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)
Patient Characteristics
Within an RCT, patient characteristics between arms should be well balanced. But between RCTs, patient characteristics can vary quite a bit. All potential effect modifiers should be assessed prior to running network meta-analysis. BUGSnet allows for a visual display of potential effect modifiers using data.plot(). Here we can see that age is well balanced between trials and treatments. Note that an error will be produced due to the fact that standard deviation was not reported in 2 studies (hence the missing error bars in the plot).
data.plot(data = rate2.slr,
covariate = "age",
half.length = "age_SD",
by = "treatment",
avg.hline = TRUE, # add overall average line?
text.size = 12) +
ggplot2::ylab("Age (years)") # customize y-axis label
If there are any imbalances in effect modifiers between trials, subgroup analysis or meta-regression can be employed. For this example, we will assume that there is no significant imbalance in effect modifiers.
For a guide on how to perform meta-regression in BUGSnet, consult BUGSnet's documentation and other vignettes.
Main analysis
Model Choice
In the case where study followup time is reported along with the number of events we need to employ the "binomial" family along with the complementary log-log (cloglog) link function. We will run both random and fixed effects models and compare their fit statistics using BUGSnet. Note that both of these models assume a constant hazard ratio for the duration of follow up. A mathematical description of the models 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, where treatment 1 is chosen to represent the reference treatment. Let $R_{ik}$ be the measured aggregate response in arm $k$ of study $i$. Let $\mathbf R = (R_{11}, ..., R_{(1a_1)}, ..., R_{M1}, ..., R_{(Ma_M)})\top$ be the vector of responses across all study arms, where $a_1,...,a_M$ represent the number of arms in studies $1,...,M$. The data is modeled using the binomial family: $R_{ik} \sim Binomial(n_{ik},\phi_{ik})$, where $\phi_{ik}$ is the probability of experiencing the event in arm $k$ of study $i$ and $n_{ik}$ is the sample size in arm $k$ of study $i$. The latent parameters $\phi_{ik}$ are transformed using a cloglog link function and modeled with the following linear model:
$$cloglog(p_{ik})=log(f_i) +\mu_i +\delta_{ik},$$
where $f_i$ is the follow-up time in study $i$, $\mu_i$ represents the fixed effect of the treatment from arm 1 in study i (a control treatment) and $\delta_{ik}$ represents the (fixed or random) effect of the treatment from arm $k$ of study $i$ relative to the treatment in arm 1. Also, $\delta_{i1}=0$ for $i=1,...,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 cloglog scale 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)}$ and $\mu_1,...,\mu_M$. 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 thought the transitivity relation $d_{(t_{i1},t_{ik})}=d_{(1,t_{ik})}-d_{(1,t_{i1})}$.
Turning back to the analysis, the first step is to use nma.model() to build a model using BUGS code based on the user's specifications. 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, we have set the referent treatment to "Diuretic" to be consistent with the TSD and literature.
random_effects_model <- nma.model(data = rate2.slr,
outcome = "diabetes",
N = "n",
reference = "Diuretic",
family = "binomial",
link = "cloglog",
time = "followup",
effects = "random")
fixed_effects_model <- nma.model(data = rate2.slr,
outcome = "diabetes",
N = "n",
reference = "Diuretic",
family = "binomial",
link = "cloglog",
time = "followup",
effects = "fixed")
If you want to review the BUGS code that was created, you can use the commands random_effects_model$model and fixed_effects_model$model (see appendix for an example). Priors are defined on $d_{(1,2)},...,d_{(1,T)}$ and $\mu_1,...,\mu_M$ 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.mu, prior.dand prior.sigma in the nma.model() function.
Once the NMA models have been specified, the next step is to run the models using nma.run() which will analyse the data with the BUGS model through the JAGS software [@JAGS]. 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(20190828)
fixed_effects_results <- nma.run(fixed_effects_model,
n.adapt = 1000,
n.burnin = 1000,
n.iter = 10000)
random_effects_results <- nma.run(random_effects_model,
n.adapt = 1000,
n.burnin = 1000,
n.iter = 10000)
Compare fixed vs random effects models by comparing leverage plots and the DIC. A plot of the $leverage_{ik}$ vs $w_{ik}$ (Bayesian deviance residual) for each of the $k$ arms in all $i$ trials can help highlight potential outliers when fitting our model. If a data point lies outside the purple arc ($x^2 + y =3$), then this data point can be said to be contributing to the model's poor fit.
par(mfrow = c(1,2))
fe_model_fit <- nma.fit(fixed_effects_results, main = "Fixed Effects Model" )
re_model_fit <- nma.fit(random_effects_results, main = "Random Effects Model")
The random effects model is obviously the better choice here. The DIC is considerably lower in the random effects model. The fixed effects model shows that 8 points are largely contributing to the model's poor fit. The random effects model appears to have only 1 outlier, which should be investigated. We can see which arm/study is reponsible for this datapoint by using the command re_model_fit$w. From this we can see that Placebo arm from the "MRC-E" study is contributing to the model's poor fit. This may call for a reconsideration of including this study in our evidence base, or to investigate potential effect modifiers. But for the sake of this example, we will assume that this study arm should remain in the evidence base and that there are no imbalanced effect modifiers.
Check inconsistency
Next, we will assess consistency in the network by fitting a random effects inconsistency model and comparing it to our random effects consistency model [@lu2006assessing]. If our inconsistency model shows a better fit than the consistency model, then it is likely that there is inconsistency in the network.
re_inconsistency_model <- nma.model(data = rate2.slr,
outcome = "diabetes",
N = "n",
reference = "Diuretic",
family = "binomial",
link = "cloglog",
time = "followup",
type = "inconsistency", #specifies inconsistency model
effects = "random")
re_inconsistency_results <- nma.run(re_inconsistency_model,
n.adapt = 1000,
n.burnin = 1000,
n.iter = 10000)
Again, rainbow plots and DIC calculations can highlight outliers and can compare model fit between the two models.
par(mfrow = c(1,2))
re_model_fit <- nma.fit(random_effects_results, main = "Consistency Model" )
inconsist_model_fit <- nma.fit(re_inconsistency_results, main = "Inconsistency Model")
When assessing the fit of both models, the consistency model has a smaller DIC, but the inconsistency model seems to have a neater leverage plot. You'll notice that there are a couple of other metrics that are outputted with the leverage plots, namely, the leverage ($pD$) and the posterior mean of the residual deviance ($\bar{D}{res}$). The leverage can be thought of as the "effective number of parameters", while the residual deviance can be thought of as the "model fit" (the lower the better). We can see that the inconsistency model produces slightly a smaller $\bar{D}{res}$, but at the expense of increased model complexity (larger $pD$). In general, the consistency model seems to be favourable here, due to it's adequate fit and parsimony.
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.
nma.compare(re_model_fit, inconsist_model_fit)
With the exception of 1 or 2 points, the data lies on or near the $y=x$ line, indicating a general agreement between the two models. This suggests that we proceed with the more parsimonious (consistency) model.
Results
Treatment Rankings
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.
The ranking probabilities can be displayed in the form of a rankogram.
sucra.out <- nma.rank(random_effects_results, largerbetter=FALSE, sucra.palette = "Set1")
sucra.out$rankogram
From this rankogram, we can see that there is a $\sim$ 77% probability that ARB's are the best treatments, and a $\sim$ 80% probability that diuretics are the worst treatment. If you want the exact ranking percentages, these results are also available in tabular form via sucra.out$ranktable.
knitr::kable(sucra.out$ranktable)
A surface under the cumulative ranking curve (SUCRA) plot can also be displayed via sucra.out$sucraplot.
sucra.out$sucraplot
If you want the exact percentages, these results are also available in tabular form via sucra.out$sucratable. The last row contains the SUCRA percentages calculated as the average of the cumulative percentages for ranks 1 to 5.
knitr::kable(sucra.out$sucratable)
Since those plots were produced using ggplot2, one can easily customize the plots using ggplot2 commands, for instance via sucra.out$sucraplot + theme_gray():
sucra.out$sucraplot + ggplot2::theme_gray()
League Table and Heatplot
League tables are another great 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(random_effects_results,
central.tdcy = "median",
order = sucra.out$order,
log.scale = FALSE,
low.colour = "springgreen4",
mid.colour = "white",
high.colour = "red",
digits = 2)
league.out$heatplot
Because we used the option log.scale = FALSE, the values reported are on the hazard ratio scale as opposed to the log hazard ratio scale. In this example, a green cell indicates that a treatment performed better than its comparator (estimate smaller than 1), while a red cell indicates that the treatment performed worst than its comparator (estimate greater than 1). The symbols (**) are used to highlight credible intervals that do not contain the neutral value 1, meaning that there is evidence of a statistically significant difference between the treatment and its comparator at the 95% confidence level. The values can also be outputted in a table via the command league.out$table.
knitr::kable(league.out$table)
Forest Plot
Forest plots are another great way to summarize the results of an NMA with respect to a particular comparator.
nma.forest(random_effects_results,
central.tdcy = "median",
comparator = "Placebo",
log.scale = FALSE)
Conclusion
A random effects model showed to be adequate in modelling the treatment effects across the evidence network. Both the SUCRA plot and league table suggest that all treatments are effective when compared diuretic, with the exception of Beta Blockers. ARB's seem to be the most effective treatment, followed by ACE inhibitors and Placebo.
Appendix
BUGS code for final model
The following BUGS code is produced by the nma.model() function.
cat(random_effects_model$bugs)
Assessment of convergence of final model
After a model is fit to the evidence, check the trace plots of the model using nma.diag(). Look for spikes or general irregularities in the posterior distributions of the $d$'s. Irregularities in these distributions suggest that the model never converged, and will likely require you to run more have a longer burn-in period, and more iterations. If you wish to know which treatments are causing irregularities, you can use the command random_effects_model$trt.key to see which treatment correspond to which integer number on the traceplot. Here, the model shows strong signs of convergence.
nma.diag(random_effects_results, plot_prompt = FALSE)
References
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 = "#>" )
library(BUGSnet)
The purpose of this example is to demonstrate how to use BUGSnet to perform network meta-analysis on time-to-event data, where each study has a different followup time. The data was taken from the NICE DSU TSD2 [@TSD2], page 66. The results of our package are concordant with those reported in the TSD.
The Diabetes data has 22 studies and 6 treatments, each study has a follow up time between 1 and 6.1 years. The outcome for this data (diabetes.sim) is the number of new diabetes cases observed during the trial period. An age variable, including its standard deviation, was simulated for demonstration purposes. This data contains multi-arm trials, which can easily be accomodated using BUGSnet. The data is shown here:
data(diabetes.sim) knitr::kable(diabetes.sim)
Before any analysis is done, we need to put our data in a form that BUGSnet can easily handle. This is done by simple running the data.prep() function and specifying the name of the treatment, study, and sample size variables:
library(BUGSnet) data(diabetes.sim) rate2.slr <- data.prep(arm.data = diabetes.sim, varname.t = "Treatment", varname.s = "Study")
We will use this data.prep object throughout the rest of our example.
Network of Evidence
To generate a network plot, simply input the data.prep object into the net.plot() function as shown:
par(mar=c(1,1,1,1)) #reduce margins for nicer plotting net.plot(rate2.slr, node.scale = 3, edge.scale=1.5)
Our network has 6 treatments and a high degree of connectedness. If a user wants to highlight a specific node and its direct evidence, this can be done using the flag option.
par(mar=c(1,1,1,1)) #reduce margins for nicer plotting net.plot(rate2.slr, node.scale = 3, edge.scale=1.5, flag = "ARB")
This gives us the names of the 5 trials connecting the ARB node to 4 treatments.
Generate Network Characteristics via net.tab()
In addition to a plot of the network, BUGSnet can summarize network, intervention, and comparison statistics that can be used to better understand our network of evidence.
network.char <- net.tab(data = rate2.slr, outcome = "diabetes", N = "n", type.outcome = "rate2", time = "followup")
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)
Patient Characteristics
Within an RCT, patient characteristics between arms should be well balanced. But between RCTs, patient characteristics can vary quite a bit. All potential effect modifiers should be assessed prior to running network meta-analysis. BUGSnet allows for a visual display of potential effect modifiers using data.plot(). Here we can see that age is well balanced between trials and treatments. Note that an error will be produced due to the fact that standard deviation was not reported in 2 studies (hence the missing error bars in the plot).
data.plot(data = rate2.slr, covariate = "age", half.length = "age_SD", by = "treatment", avg.hline = TRUE, # add overall average line? text.size = 12) + ggplot2::ylab("Age (years)") # customize y-axis label
If there are any imbalances in effect modifiers between trials, subgroup analysis or meta-regression can be employed. For this example, we will assume that there is no significant imbalance in effect modifiers.
For a guide on how to perform meta-regression in BUGSnet, consult BUGSnet's documentation and other vignettes.
Main analysis
Model Choice
In the case where study followup time is reported along with the number of events we need to employ the "binomial" family along with the complementary log-log (cloglog) link function. We will run both random and fixed effects models and compare their fit statistics using BUGSnet. Note that both of these models assume a constant hazard ratio for the duration of follow up. A mathematical description of the models 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, where treatment 1 is chosen to represent the reference treatment. Let $R_{ik}$ be the measured aggregate response in arm $k$ of study $i$. Let $\mathbf R = (R_{11}, ..., R_{(1a_1)}, ..., R_{M1}, ..., R_{(Ma_M)})\top$ be the vector of responses across all study arms, where $a_1,...,a_M$ represent the number of arms in studies $1,...,M$. The data is modeled using the binomial family: $R_{ik} \sim Binomial(n_{ik},\phi_{ik})$, where $\phi_{ik}$ is the probability of experiencing the event in arm $k$ of study $i$ and $n_{ik}$ is the sample size in arm $k$ of study $i$. The latent parameters $\phi_{ik}$ are transformed using a cloglog link function and modeled with the following linear model:
$$cloglog(p_{ik})=log(f_i) +\mu_i +\delta_{ik},$$
where $f_i$ is the follow-up time in study $i$, $\mu_i$ represents the fixed effect of the treatment from arm 1 in study i (a control treatment) and $\delta_{ik}$ represents the (fixed or random) effect of the treatment from arm $k$ of study $i$ relative to the treatment in arm 1. Also, $\delta_{i1}=0$ for $i=1,...,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 cloglog scale 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)}$ and $\mu_1,...,\mu_M$. 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 thought the transitivity relation $d_{(t_{i1},t_{ik})}=d_{(1,t_{ik})}-d_{(1,t_{i1})}$.
Turning back to the analysis, the first step is to use nma.model() to build a model using BUGS code based on the user's specifications. 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, we have set the referent treatment to "Diuretic" to be consistent with the TSD and literature.
random_effects_model <- nma.model(data = rate2.slr, outcome = "diabetes", N = "n", reference = "Diuretic", family = "binomial", link = "cloglog", time = "followup", effects = "random") fixed_effects_model <- nma.model(data = rate2.slr, outcome = "diabetes", N = "n", reference = "Diuretic", family = "binomial", link = "cloglog", time = "followup", effects = "fixed")
If you want to review the BUGS code that was created, you can use the commands random_effects_model$model and fixed_effects_model$model (see appendix for an example). Priors are defined on $d_{(1,2)},...,d_{(1,T)}$ and $\mu_1,...,\mu_M$ 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.mu, prior.dand prior.sigma in the nma.model() function.
Once the NMA models have been specified, the next step is to run the models using nma.run() which will analyse the data with the BUGS model through the JAGS software [@JAGS]. 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(20190828) fixed_effects_results <- nma.run(fixed_effects_model, n.adapt = 1000, n.burnin = 1000, n.iter = 10000) random_effects_results <- nma.run(random_effects_model, n.adapt = 1000, n.burnin = 1000, n.iter = 10000)
Compare fixed vs random effects models by comparing leverage plots and the DIC. A plot of the $leverage_{ik}$ vs $w_{ik}$ (Bayesian deviance residual) for each of the $k$ arms in all $i$ trials can help highlight potential outliers when fitting our model. If a data point lies outside the purple arc ($x^2 + y =3$), then this data point can be said to be contributing to the model's poor fit.
par(mfrow = c(1,2)) fe_model_fit <- nma.fit(fixed_effects_results, main = "Fixed Effects Model" ) re_model_fit <- nma.fit(random_effects_results, main = "Random Effects Model")
The random effects model is obviously the better choice here. The DIC is considerably lower in the random effects model. The fixed effects model shows that 8 points are largely contributing to the model's poor fit. The random effects model appears to have only 1 outlier, which should be investigated. We can see which arm/study is reponsible for this datapoint by using the command re_model_fit$w. From this we can see that Placebo arm from the "MRC-E" study is contributing to the model's poor fit. This may call for a reconsideration of including this study in our evidence base, or to investigate potential effect modifiers. But for the sake of this example, we will assume that this study arm should remain in the evidence base and that there are no imbalanced effect modifiers.
Check inconsistency
Next, we will assess consistency in the network by fitting a random effects inconsistency model and comparing it to our random effects consistency model [@lu2006assessing]. If our inconsistency model shows a better fit than the consistency model, then it is likely that there is inconsistency in the network.
re_inconsistency_model <- nma.model(data = rate2.slr, outcome = "diabetes", N = "n", reference = "Diuretic", family = "binomial", link = "cloglog", time = "followup", type = "inconsistency", #specifies inconsistency model effects = "random") re_inconsistency_results <- nma.run(re_inconsistency_model, n.adapt = 1000, n.burnin = 1000, n.iter = 10000)
Again, rainbow plots and DIC calculations can highlight outliers and can compare model fit between the two models.
par(mfrow = c(1,2)) re_model_fit <- nma.fit(random_effects_results, main = "Consistency Model" ) inconsist_model_fit <- nma.fit(re_inconsistency_results, main = "Inconsistency Model")
When assessing the fit of both models, the consistency model has a smaller DIC, but the inconsistency model seems to have a neater leverage plot. You'll notice that there are a couple of other metrics that are outputted with the leverage plots, namely, the leverage ($pD$) and the posterior mean of the residual deviance ($\bar{D}{res}$). The leverage can be thought of as the "effective number of parameters", while the residual deviance can be thought of as the "model fit" (the lower the better). We can see that the inconsistency model produces slightly a smaller $\bar{D}{res}$, but at the expense of increased model complexity (larger $pD$). In general, the consistency model seems to be favourable here, due to it's adequate fit and parsimony.
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.
nma.compare(re_model_fit, inconsist_model_fit)
With the exception of 1 or 2 points, the data lies on or near the $y=x$ line, indicating a general agreement between the two models. This suggests that we proceed with the more parsimonious (consistency) model.
Results
Treatment Rankings
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.
The ranking probabilities can be displayed in the form of a rankogram.
sucra.out <- nma.rank(random_effects_results, largerbetter=FALSE, sucra.palette = "Set1") sucra.out$rankogram
From this rankogram, we can see that there is a $\sim$ 77% probability that ARB's are the best treatments, and a $\sim$ 80% probability that diuretics are the worst treatment. If you want the exact ranking percentages, these results are also available in tabular form via sucra.out$ranktable.
knitr::kable(sucra.out$ranktable)
A surface under the cumulative ranking curve (SUCRA) plot can also be displayed via sucra.out$sucraplot.
sucra.out$sucraplot
If you want the exact percentages, these results are also available in tabular form via sucra.out$sucratable. The last row contains the SUCRA percentages calculated as the average of the cumulative percentages for ranks 1 to 5.
knitr::kable(sucra.out$sucratable)
Since those plots were produced using ggplot2, one can easily customize the plots using ggplot2 commands, for instance via sucra.out$sucraplot + theme_gray():
sucra.out$sucraplot + ggplot2::theme_gray()
League Table and Heatplot
League tables are another great 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(random_effects_results, central.tdcy = "median", order = sucra.out$order, log.scale = FALSE, low.colour = "springgreen4", mid.colour = "white", high.colour = "red", digits = 2) league.out$heatplot
Because we used the option log.scale = FALSE, the values reported are on the hazard ratio scale as opposed to the log hazard ratio scale. In this example, a green cell indicates that a treatment performed better than its comparator (estimate smaller than 1), while a red cell indicates that the treatment performed worst than its comparator (estimate greater than 1). The symbols (**) are used to highlight credible intervals that do not contain the neutral value 1, meaning that there is evidence of a statistically significant difference between the treatment and its comparator at the 95% confidence level. The values can also be outputted in a table via the command league.out$table.
knitr::kable(league.out$table)
Forest Plot
Forest plots are another great way to summarize the results of an NMA with respect to a particular comparator.
nma.forest(random_effects_results, central.tdcy = "median", comparator = "Placebo", log.scale = FALSE)
Conclusion
A random effects model showed to be adequate in modelling the treatment effects across the evidence network. Both the SUCRA plot and league table suggest that all treatments are effective when compared diuretic, with the exception of Beta Blockers. ARB's seem to be the most effective treatment, followed by ACE inhibitors and Placebo.
Appendix
BUGS code for final model
The following BUGS code is produced by the nma.model() function.
cat(random_effects_model$bugs)
Assessment of convergence of final model
After a model is fit to the evidence, check the trace plots of the model using nma.diag(). Look for spikes or general irregularities in the posterior distributions of the $d$'s. Irregularities in these distributions suggest that the model never converged, and will likely require you to run more have a longer burn-in period, and more iterations. If you wish to know which treatments are causing irregularities, you can use the command random_effects_model$trt.key to see which treatment correspond to which integer number on the traceplot. Here, the model shows strong signs of convergence.
nma.diag(random_effects_results, plot_prompt = FALSE)
References
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.