WINBUGS to R: Replicating a published Bayesian network meta-analysis in R
Bayesian Network Meta-analyses in R
Network meta-analyses (NMA) are a powerful statistical tool frequently employed in medical research. As alluded to in the overview of NMAs, there are two different types, namely frequentist or Bayesian NMAs. Here, I'm going to focus on Bayesian NMAs, where I've written a short introductory guide to using the gemtc package in R to conduct Bayesian-NMAs. Rather than use a well-known or common dataset to showcase this approach I've decided to replicate an previously published analysis.
The recently published study by Langford et al. (2020) uses a Bayesian-NMA to compare the effect of various sodium-glucose co-transporter (SGLT) inhibitors and metformin in patients with type 1 diabetes mellitus (T1DM). In doing so they collate data from nine studies for the change from baseline (CFB) of HbA1c (average level of glycated-hameoglobin), weight, total daily insulin dose, and systolic blood pressure at two time points (24-26 weeks and 52 weeks). To analyse their data using a Bayesian-NMA the authors used the WINBUGS software. Accordingly, here we will replicate their analysis in R using the gemtc package and compare our results to those from the original paper. Rather than focus on all four of the metrics, we are just going to focus on the effect of the various SGLT inhibitor treatments on HbA1c, although the same approach could be used for all other response metrics.
In this tutorial I'm solely going to focus on the code used to conduct this analysis as opposed to the biological interpretation of the results. If you're interested in the medical implications of these results I suggest that you read the original paper!
As I mentioned above, the data we will be using here is for the CFB in HbA1c under each SGLT inhibitor treatment. For these studies, the aim was for any SGLT inhibitor treatment to reduce the level of HbA1c. Hence, if a SGLT inhibitor treatment has a beneficial effect on HbA1c, the effect size should be negative. Here, CFB in HbA1c is measured in percentage (%) and the effect sizes we have calculated are mean differences. In the original study, Langford et al. compare each SGLT inhibitor treatment to either a placebo or to the standard treatment of Metformin (2000mg), accordingly in this tutorial we will replicate both of these analyses.
Within this tutorial I won't go into the theory surronding Bayes theroem, but there are plenty of great resources online. Likewise, I'm going to assume that you have a basic understanding of meta-analyses, although if not have a look at this extensive resource. This tutorial will use R to conduct Bayesian-NMAs, however there are a variety of different programs which can be used. Indeed, the recently published NICE Meta-Analysis Guideline Methodology Document outlines the various different options that could be used for NMAs (as well as those which should be avoided).
Packages and data
To start our analysis we first need to load the various packages we will be using. Additionally, we will set the seed so that our analysis will be exactly reproducible.
#install the various packages needed for the tutorial
#all of these packages are available via CRAN and can be installed via the install.packages() function
library(gemtc)
library(rjags)
library(data.table)
library(dmetar)
library(igraph)
library(viridis)
library(ggplot2)
library(patchwork)
#set the seed to ensure reproducibility of the analyses
set.seed(36361)Langford et al. (2020) make their data openly available in the supplementary information of their paper which I have put into a csv file for this analysis. As a quick check of their data, I also extracted the data from their studies myself and found that it aligned well with the data provided by Langford et al. (2020). From the data provided by Langford et al. (2020), I have calculated effect sizes (mean difference) for each SGLT inhibitor treatment relative to the placebo treatment. Likewise, I have formatted the data into the format required by gemtc.
The dataset I've loaded contains data for both time points (24-26 weeks or 52 weeks). However, we will be conducting seperate meta-analyses for each of these time points and so the data needs to be split at this point. Additionally, we will remove a couple of additional columns in the datasets which we no longer require.
#load in the data for this analysis
df <- fread("C:/Users/Burgess/Documents/Network_MetaAnalysis/Langford_2020_gemtcstyle.csv")
#show the first ten rows of the dataset
head(df, n = 6)## study metric time treatment diff std.err
## 1: Lund2008 HbA1c 52 Metformin_2000 0.13 0.15556349
## 2: Lund2008 HbA1c 52 Placebo NA NA
## 3: REMOVAL HbA1c 24 Metformin_2000 -0.15 0.07071068
## 4: REMOVAL HbA1c 52 Metformin_2000 -0.06 0.06403124
## 5: REMOVAL HbA1c 24 Placebo NA NA
## 6: REMOVAL HbA1c 52 Placebo NA NA#seperate out the dataset depending on whether the data
#corresponds to 24-26 or 52 weeks
df_24 <- subset(df, time == 24)
df_52 <- subset(df, time == 52)
#remove the columns no longer needed
df_24 <- subset(df_24, select = -c(time, metric))
df_52 <- subset(df_52, select = -c(time, metric))At this point we will also load data on the cumulative number of patients for each of the SGLT inhibitor treatments at each time point. We're not going to need this for our analysis, but will be using it to help visualise our data in the network plot below.
#load up data for the cumulative number of patients per treatment
#this is to be used in our network plot
node_weights <- fread("C:/Users/Burgess/Documents/Network_MetaAnalysis/Node_Weights.csv")
#seperate out the dataset depending on whether the data
#corresponds to 24-26 or 52 weeks
node_weights_24 <- subset(node_weights, Time == 24)$Total
node_weights_52 <- subset(node_weights, Time == 52)$TotalHbA1c Analysis: 24-26 week data
We'll start by analysing the data for the 24-26 week time point. The first step we need to do is generate the network for our dataset which is done using the mtc.network function in the gemtc package. The output of this function produces our network which contains a variety of data. At this point we will edit the description of each of our treatments in our network. We'll do this as it makes any output clearer and easier to read.
#create the network for the 24-26 weeks data
network_24 <- mtc.network(data.re = df_24)
#format the treatment names in the description component of the network
description_names <- network_24$treatments$description
description_names <- paste0(description_names, "mg")
description_names <- gsub("_", " ", description_names)
description_names <- gsub("2 5", "2.5", description_names)
description_names <- gsub("Placebomg", "Placebo", description_names)
network_24$treatments$description <- description_names
rm(description_names)At this point we can generate the network plot for our network. Which is done using the below code.
#set colours for the rest of the analyses
#colours have been selected based upon the treatments
#the viridis package has been used to choose colours
vertex_colours_24 <- c(viridis(4)[4],
viridis(4)[4],
viridis(4)[2],
viridis(4)[2],
viridis(4)[2],
viridis(4)[1],
"white",
viridis(4)[3],
viridis(4)[3])
#generate a network plot
plot(network_24,
use.description = TRUE,
vertex.color = vertex_colours_24,
vertex.label.color = "gray10",
vertex.shape = "circle",
vertex.label.family = "Arial",
vertex.size = (node_weights_24)/100, #based on cumulative number of patients per treatment
vertex.label.dist = -2.5,
vertex.label.cex = 0.75,
layout = layout.fruchterman.reingold)
As we can clearly see, our network contains nine nodes which correspond to nine different treatments. Our treatments can be grouped into five different classes shown by the nodes of different colours. We have data for the CFB of HbA1c under five different medications at various different dosages. The treatments are: Metformin (2000 mg); Dapagliflozin (5mg, 10mg); Empagliflozin (2.5mg, 10mg, 25mg); Sotagliflozin (200mg, 400mg); Placebo.
As we discussed earlier, the size of each node corresponds to the cumulative number of patients in that group.
The next stage for our analysis is to conduct the Bayesian-NMA. To do this, we need to specify the various settings we want to use. Importantly, we are going to be running two different Bayesian-NMAs, where we specify either fixed or random effects in our models. The rationale for this is that we can compare the fit of these two models and subsequently choose the model which best fits our data.
We can specify the settings for both our models using the below code.
#generate the settings for the random effects model for the 24-26 weeks data
model_24_random <- mtc.model(network_24,
likelihood = "normal",
link = "identity",
linearModel = "random",
n.chain = 5)
#generate the settings for the fixed effects model for the 24-26 weeks data
model_24_fixed <- mtc.model(network_24,
likelihood = "normal",
link = "identity",
linearModel = "fixed",
n.chain = 5)As we have now specified our settings, we can run both our fixed and random effects models. Here, the main thing to note is that I have specified the burnin (60000) and iterations (140000) for our Bayesian-NMAs. These parameter values were chosen to mirror those selected by Langford et al. in their original analysis. Likewise, I have specified an identity link function and normal likelihood for both models, again mirroring those settings chosen by Langford et al. (2020) in their orignial analysis.
#run the bayesian network meta-analysis for the random effects model
mtc_24_random <- mtc.run(model_24_random, n.adapt = 60000, n.iter = 140000)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 8
## Unobserved stochastic nodes: 24
## Total graph size: 326
##
## Initializing model#run the bayesian network meta-analysis for the random effects model
mtc_24_fixed <- mtc.run(model_24_fixed, n.adapt = 60000, n.iter = 140000)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 8
## Unobserved stochastic nodes: 8
## Total graph size: 248
##
## Initializing model#note that the burnin and number of iterations have been chosen
#to mirror those from Langford et al. (2020)Once our models have run, we can extract the Deviance Information Criterion (DIC) for each model. This is a similar metric to AIC or BIC and can be used to compare the fit of two models. As with AIC or BIC, a lower value of DIC indicates a better fitting model. As such, using the following code, we can select the better fitting model.
#extract the DIC values for the random and fixed effect models
mtc_24_random_DIC <- summary(mtc_24_random)$DIC[3]
mtc_24_fixed_DIC <- summary(mtc_24_fixed)$DIC[3]
mtc_24_random_DIC## DIC
## 22.24051mtc_24_fixed_DIC ## DIC
## 20.6937#select the best fitting model (i.e. lowest DIC values)
if(mtc_24_random_DIC < mtc_24_fixed_DIC){
mtc_24 <- mtc_24_random
}
if(mtc_24_fixed_DIC <= mtc_24_random_DIC){
mtc_24 <- mtc_24_fixed
}For our models, it appears as though the Bayesian-NMA which incorporates fixed effects is the best fitting model, hence we will use the fixed effects model for the rest of the analysis of the 24-26 week HbA1c data. This results matches the findings of Langford et al. who also report that their fixed effects model had a lower DIC compared to their fixed effects model.
Additionally the \(I^{2}\) metric for the fixed effects model is 0% which indicates consistency in effect sizes across the different studies.
Trace and density plots
Before we consider the results of our NMA, we can check several model attributes to ensure that our model is behaving as we expect it to. The first is done using a trace plot. We can use this to check whether our model is coverging as we expect it to and that there no obvious long term trends.
The second attribute we can consider, is a density plot for the estimates of effect sizes which we would expect to be a normal distribution.
We can generate both these plots of all the different comparisons in our dataset (e.g., Metformin 2000mg v Placebo, Dapagliflozin 5mg v Placebo etc.), but in the interests of clarity we will only consider the figures for the comparison of Metformin (2000mg) to the Placebo.
#generate the trace plot for each treatments
#plot(mtc_24)
#Only generate the trace figure for metformin relative to the placebo
plot(mtc_24$samples[,6])
As we can see, both plots appear as we would expect, with there being no obvious trends in the trace plot and gaussian distribution for our density plot.
Potential Scale Reduction Factor
Another attribute of our NMA which we can assess is the Potential Scale Reduction Factor (PSRF) which compares any varitation both within and between our chains. Ideally, this value should be no more than 1.05 and close to 1 for the larger number of iterations. We can assess this using Gelman-Rubin plots. Although we can generate these plots for each comparison (e.g., Metformin 2000mg v Placebo, Dapagliflozin 5mg v Placebo etc.) we will again only show the comparison of Metformin 2000mg v Placebo for clarity.
#generate the Gelman-Rubin plots for each treatment
#gelman.plot(mtc_24)
#Only generate the Gelman-Rubin plot for metformin relative to the placebo
gelman.plot(mtc_24$samples[,6])
We can also extract an overall PSRF for our models using the following code.
#calculate the overall PSRF
gelman.diag(mtc_24)$mpsrf## [1] 1.000044As we can see, the value of this PSRF is very close to 1.
Treatment rankings
Given that we have now assessed various attributes of our model, we can begin to analyse the results. The first result we can consider is the ranking of our different treatments - although this was not previously reported by Langford et al. (2020). To do this we can generate a rankogram using the following code. The interpretation of this figure is relatively straightforward. The treatment with the largest bar for rank 1 is likely to be the most effective treatment, while the largest bar for rank 2 is likely to be represent the second most effective treatment etc.
#Calculate the rank data
mtc_24_rank <- rank.probability(mtc_24, preferredDirection = -1)
#convert our rank data into a plottable format
rank_pos <- c()
rank_prob <- c()
for(i in 1:length(network_24$treatments$id)){
rank_pos <- c(rank_pos, rep(i, length(network_24$treatments$id)))
rank_prob <- c(rank_prob, mtc_24_rank[,i])
}
#generate a dataframe to store the rank data
df_rank_plot_24 <- data.frame(Probability = rank_prob,
treatment_old = network_24$treatments$id,
Treatment = network_24$treatments$description,
Rank = factor(rank_pos))
#change the Treatment to a factor variable
#this is so that we can order the different treatments in the plot in the order we want
#here we are ordering treatments alphabetically and then by increasing dosage
df_rank_plot_24$Treatment <- factor(df_rank_plot_24$Treatment,
levels = df_rank_plot_24$Treatment[c(2, 1, 4, 3, 5:9)])
#generate the rankogram
rank_24_plot <- ggplot(df_rank_plot_24, aes(x=Treatment, y=Probability, fill=Rank))+
geom_bar(stat="identity", color="black", position=position_dodge())+
scale_fill_grey(start=0, end=1) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust=1))+
ylab("")
rank_24_plot
As we can see from our results there are some relatively strong signals apparent. Firstly, the rankogram suggests that Empagliflozin (25mg) is likely to be the most effective treatment, while Empagliflozin (10mg) is likely to be the second most effective treatment. In contrast, we can see that the Placebo and Metformin (2000mg) are likely to be the least effective and second least effective treatments, respectively. However, while the rankogram is interesting, it provides no information on the magnitude of the effect of each treatment. For this information we would need to see the results of the Bayesian-NMA.
SUCRA values
An additional way to rank the various treatments is to use a metric called Surface Under the Cumulative Ranking (SUCRA). As before, this produces a histogram which can be used to rank the various treatments - although this was not previously reported by Langford et al. (2020). The code to generate a plot of SUCRA values for each treatment is as follows.
#calculate the SUCRA values for each of the treatments using the rankings
mtc_24_sucra <- sucra(mtc_24_rank, lower.is.better = FALSE)
#extract the SUCRA data in an easily plotable format
x_24_sucra <- plot(mtc_24_sucra)
#generate a dataframe to store the SUCRA data
df_sucra_plot_24 <- data.frame(SUCRA = NA,
treatment_old = network_24$treatments$id,
Treatment = network_24$treatments$description,
colours = vertex_colours_24)
#extract the SUCRA values for each treatment
for(i in 1:dim(df_sucra_plot_24)[1]){
df_sucra_plot_24$SUCRA[i] <- subset(x_24_sucra$data,
Treatment == df_sucra_plot_24$treatment_old[i])$SUCRA
}
#change the Treatment to a factor variable
#this is so that we can order the different treatments in the plot in the order we want
#here we are ordering treatments alphabetically and then by increasing dosage
df_sucra_plot_24$Treatment <- factor(df_sucra_plot_24$Treatment,
levels = df_sucra_plot_24$Treatment[c(2, 1, 4, 3, 5:9)])
#generate the SUCRA plot using ggplot
sucra_24_plot <- ggplot(df_sucra_plot_24, aes(x=Treatment, y=SUCRA, fill=Treatment))+
geom_bar(stat="identity", color="black")+
scale_fill_manual(values=df_sucra_plot_24$colours)+
theme_minimal() +
theme(legend.position = "None") +
theme(axis.text.x = element_text(angle = 45, hjust=1))
sucra_24_plot
In line with the findings of the previous results, we can see that using SUCRA values that Empagliflozin (25mg) is likely to be the most effective treatment, while Empagliflozin (10mg) is likely to be the second most effective treatment. Again, as before, the lowest SUCRA values are for Metformin (2000mg) and the Placebo, indicating that these are likely to be the least effective treatments.
Forest plots
In their paper, Langford et al. report forest plots showcasing how the various SGLT inhibitor treatments compare either to the placebo treatment or the Metformin (2000mg) treatment. As such, we can likewise generate these forest for our NMA and compare our results to those reported by Langford et al. (2020).
Comparison of treatments to placebo
We'll start by generating the forest plot showing the effect sizes (mean difference) in CFB HbA1c (%) for each SGLT inhibitor treatment and Metformin in comparison to the placebo treatment. The code to generate this figure is a little complex, although it provides a clear figure using ggplot in contrast to the less clear figure generated by the gemtc package.
#extract the summary information where we calculate the mean difference in HbAc1 for each treatment relative to the placebo
summary_placebo_24 <- summary(relative.effect(mtc_24, t1 = "Placebo"))
#generate a dataframe to store the data
df_forest_placebo_24 <- data.frame(Treatment = network_24$treatments$description,
colours = vertex_colours_24,
Mean = NA,
CrI_Lower = NA,
CrI_Upper = NA)
#remove any data corresponding to the placebo treatment from our dataframe
df_forest_placebo_24 <- subset(df_forest_placebo_24, Treatment != "Placebo")
#calculate the dimensions of our dataframe - this is to shorten the code later
dim_df <- dim(df_forest_placebo_24)[1]
#Extract mean values for the estimate of the effect of each treatment
df_forest_placebo_24$Mean <- summary_placebo_24$summaries$statistics[1:dim_df]
#Extract lower bounds for the credible intervals of each treatment effect
df_forest_placebo_24$CrI_Lower <- summary_placebo_24$summaries$quantiles[1:dim_df]
#Extract upper bounds for the credible intervals of each treatment effect
df_forest_placebo_24$CrI_Upper <- summary_placebo_24$summaries$quantiles[((4*dim_df)+1):(5*dim_df)]
#change the Treatment to a factor variable
#this is so that we can order the different treatments in the plot in the order we want
#here we are ordering treatments in same order as Langford et al. (2020)
df_forest_placebo_24$Treatment <- factor(df_forest_placebo_24$Treatment,
levels = df_forest_placebo_24$Treatment[c(6, 2, 1, 7, 8, 4, 3, 5)])
#generate our forest plot for placebo comparisons
forest_placebo_24_plot <- ggplot(df_forest_placebo_24, aes(x=Mean,
y=Treatment,
fill = Treatment)) +
geom_errorbar(aes(y = Treatment, xmin=CrI_Lower, xmax=CrI_Upper, width = 0)) +
geom_point(size = 5, shape = 21, color = "black") +
scale_fill_manual(values = df_forest_placebo_24$colours[c(6, 2, 1, 7, 8, 4, 3, 5)]) +
theme_classic() +
theme(legend.position = "none") +
geom_vline(colour = "gray50", xintercept = 0, linetype="dashed") +
xlab("Mean difference in CFB HbA1c relative to placebo (%)") +
scale_y_discrete(limits=rev) +
ggtitle("A) Week 24-26")
forest_placebo_24_plot
This forest plot shows the estimated effect sizes for the mean difference (%) in CFB HbA1c for each SGLT inhibitor treatment and Metformin relative to the placebo. Here, each effect size estimate has been determined using our fixed effects Bayesian-NMA, with the error bars representing the 95% credible intervals for ech effect size. Our analysis suggests that all treatments, at 24-26 weeks, led to significant reductions in HbA1c compared to the placebo treatment. However, the magnitude of any effect substantially varied between treatments. In line with the results of the rank and SUCRA analysis, the effect size for Metformin (2000mg) is substantially less (compared to the placebo) than all other treatments. Likewise the reduction in HbA1c was greatest for the higher dosage treatments of Empagliflozin (10mg, 25mg).
The figure we have generated here can be directly compared to Figure 3a from Langford et al. (2020). By comparing the figures, we can see the results of our analysis match those reported by Langford et al. (2020) with any minor differences likely to be due to the stochastic nature of the analysis or differences between R and WINBUGS.
Comparison of treatments to Metformin (2000mg)
Next, we'll generate the forest plot showing the effect sizes (mean difference) in CFB HbA1c (%) for each treatment in comparison to the Metformin (2000mg) treatment. This can be done using the following code.
#extract the summary information where we calculate the mean difference in HbAc1 for each treatment relative to the Metformin treatment
summary_metformin_24 <- summary(relative.effect(mtc_24, t1 = "Metformin_2000"))
#generate a dataframe to store the data
df_forest_metformin_24 <- data.frame(Treatment = network_24$treatments$description,
colours = vertex_colours_24,
Mean = NA,
CrI_Lower = NA,
CrI_Upper = NA)
#remove any data corresponding to the Metformin treatment from our dataframe
df_forest_metformin_24 <- subset(df_forest_metformin_24, Treatment != "Metformin 2000mg")
#calculate the dimensions of our dataframe - this is to shorten the code later
dim_df <- dim(df_forest_metformin_24)[1]
#Extract mean values for the estimate of the effect of each treatment
df_forest_metformin_24$Mean <- summary_metformin_24$summaries$statistics[1:dim_df]
#Extract lower bounds for the credible intervals of each treatment effect
df_forest_metformin_24$CrI_Lower <- summary_metformin_24$summaries$quantiles[1:dim_df]
#Extract upper bounds for the credible intervals of each treatment effect
df_forest_metformin_24$CrI_Upper <- summary_metformin_24$summaries$quantiles[((4*dim_df)+1):(5*dim_df)]
#change the Treatment to a factor variable
#this is so that we can order the different treatments in the plot in the order we want
#here we are ordering treatments in same order as Langford et al. (2020)
df_forest_metformin_24$Treatment <- factor(df_forest_metformin_24$Treatment,
levels = df_forest_metformin_24$Treatment[c(6, 2, 1, 7, 8, 4, 3, 5)])
#generate our forest plot for metformin comparisons
forest_metformin_24_plot <- ggplot(df_forest_metformin_24, aes(x=Mean,
y=Treatment,
fill = Treatment)) +
geom_errorbar(aes(y = Treatment, xmin=CrI_Lower, xmax=CrI_Upper, width = 0)) +
geom_point(size = 5, shape = 21, color = "black") +
scale_fill_manual(values = df_forest_metformin_24$colours[c(6, 2, 1, 7, 8, 4, 3, 5)]) +
theme_classic() +
theme(legend.position = "none") +
geom_vline(colour = "gray50", xintercept = 0, linetype="dashed") +
xlab("Mean difference in CFB HbA1c relative to Metformin 2000mg (%)") +
scale_y_discrete(limits=rev) +
ggtitle("A) Week 24-26")
forest_metformin_24_plot
This forest plot shows the estimated effect sizes for the mean difference (%) in CFB HbA1c for each SGLT inhibitor treatment and the placebo relative to Metformin (2000mg). Here, each effect size estimate has been determined using our fixed effects Bayesian-NMA, with the error bars representing the 95% credible intervals for ech effect size. Our analysis suggests that most treatments, at 24-26 weeks, led to significant reductions in HbA1c compared to the Metformin (2000mg) treatment. However, one exception to this was that the Metformin (2000mg) treatment was significantly better at reducing HbA1c compared to the placebo treatment. Likewise, given that the 95% credible interval for the Empagliflozin (2.5mg) treatment overlaps zero, there is no statistically significant difference between the Empagliflozin (2.5mg) and Metformin (2000mg) treatments. Additionally, magnitude of any effect substantially varied between treatments, although there is substantial overlap between the credible intervals of each treatment.
The figure we have generated here can be directly compared to Figure 4a from Langford et al. (2020). By comparing the figures, we can see the results of our analysis match those reported by Langford et al. (2020) with any minor differences likely to be due to the stochastic nature of the analysis or differences between R and WINBUGS.
Importantly, we can interpret these results in line with the network plot we generated earlier. As we can see, with the exception of the placebo treatment, there are no direct comparisons between Metformin (2000mg) and any other SGLT inhibitor treatment. Accordingly, as noted by Langford et al., this means that all of the reported effect sizes are generated from indirect comparisons and as such are potentially liable to bias. This is an important consideration when interpretting these results.
HbA1c Analysis: 52 week data
Alongside the analysis of the 24-26 week data, Langford et al. (2020) also undertake a Bayesian-NMA for the data at the 52 week time point. This analysis is detailed below. However, for a comprehensive explanation of each step, refer to the 24-26 week analysis.
Firstly, we create the network for the 52 week data, provide thorough descriptions for each treatment, assign vertex colours, and generate the network plot.
#create the network for the 52 week data
network_52 <- mtc.network(data.re = df_52)
#format the treatment names in the description component of the network
description_names <- network_52$treatments$description
description_names <- paste0(description_names, "mg")
description_names <- gsub("_", " ", description_names)
description_names <- gsub("Placebomg", "Placebo", description_names)
network_52$treatments$description <- description_names
rm(description_names)
#set colours for the rest of the analyses
#colours have been selected based upon the treatments
#the viridis package has been used to choose colours
vertex_colours_52 <- c(viridis(4)[4],
viridis(4)[4],
viridis(4)[2],
viridis(4)[2],
viridis(4)[1],
"white",
viridis(4)[3],
viridis(4)[3])
#generate a network plot
plot(network_52,
use.description = TRUE,
vertex.color = vertex_colours_52,
vertex.label.color = "gray10",
vertex.shape = "circle",
vertex.label.family = "Arial",
vertex.size = (node_weights_52)/100, #based on cumulative number of patients per treatment
vertex.label.dist = -2.5,
vertex.label.cex = 0.75,
layout = layout.fruchterman.reingold)
We can see that the structure of the network is similar to that for the 24-26 week data. However, the major difference is that there is no data at 52 weeks for the Empagliflozin (2.5mg) treatment. As such, the subsequent Bayesian-NMA differs to the above analysis in that respect.
The next step is to specify the settings for the Bayesian-NMAs we are going to conduct. As before, we will be conducting two NMAs, with contain either a fixed or random effect. As before, we will use the DIC of each model to choose the Bayesian-NMA which best fits our data.
#generate the settings for the random effects model for the 24-26 weeks data
model_52_random <- mtc.model(network_52,
likelihood = "normal",
link = "identity",
linearModel = "random",
n.chain = 5)
#generate the settings for the fixed effects model for the 24-26 weeks data
model_52_fixed <- mtc.model(network_52,
likelihood = "normal",
link = "identity",
linearModel = "fixed",
n.chain = 5)We can run the Bayesian-NMAs using the following code, which again uses the same burnin period, number of iterations, link function, and likelihood as that specfied by Langford et al. (2020).
#run the bayesian network meta-analysis for the random effects model
mtc_52_random <- mtc.run(model_52_random, n.adapt = 60000, n.iter = 140000)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 7
## Unobserved stochastic nodes: 20
## Total graph size: 263
##
## Initializing model#run the bayesian network meta-analysis for the random effects model
mtc_52_fixed <- mtc.run(model_52_fixed, n.adapt = 60000, n.iter = 140000)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 7
## Unobserved stochastic nodes: 7
## Total graph size: 203
##
## Initializing model#note that the burnin and number of iterations have been chosen
#to mirror those from Langford et al. (2020)Here, we will extract the DIC values and use them to determine the best fitting NMA.
#extract the DIC values for the random and fixed effect models
mtc_52_random_DIC <- summary(mtc_52_random)$DIC[3]
mtc_52_fixed_DIC <- summary(mtc_52_fixed)$DIC[3]
mtc_52_random_DIC## DIC
## 19.04535mtc_52_fixed_DIC ## DIC
## 17.27502#select the best fitting model (i.e. lowest DIC values)
if(mtc_52_random_DIC < mtc_52_fixed_DIC){
mtc_52 <- mtc_52_random
}
if(mtc_52_fixed_DIC < mtc_52_random_DIC){
mtc_52 <- mtc_52_fixed
}As with the 24-26 week data, the best fitting NMA was the fixed effects model which had a lower DIC than the random effects NMA. Again, this mirrors the findings of Langford et al. (2020).
Additionally the \(I^{2}\) metric for the fixed effects model is 0% which indicates consistency in effect sizes across the different studies.
Trace and density plots
As before we can assess various model attributes. Here, we can assess the trace and density plots for the various comparisons in our dataset. As before, in the interests of clarity, only the trace and density plots for the comparison of Metformin (2000mg) to the placebo are shown.
#generate the trace plot for each treatments
#plot(mtc_52)
#Only generate the trace figure for metformin relative to the placebo
plot(mtc_52$samples[,5])
As before, there are no trends in the trace plot and our density plot shows a gaussian distribution. These attributes indicate that our NMA is acting as expected.
Potential Scale Reduction Factor
As before, we can also generate Gelman-Rubin plots for our NMA, though only the Gelman-Rubin plot for Metformin v the placebo is shown.
#generate the Gelman-Rubin plots for each treatment
#gelman.plot(mtc_52)
#Only generate the Gelman-Rubin plot for metformin relative to the placebo
gelman.plot(mtc_52$samples[,5])
We can also extract an overall PSRF for our models using the following code.
#calculate the overall PSRF
gelman.diag(mtc_52)$mpsrf## [1] 1.000031As we can see, the value of this PSRF is very close to 1.
Treatment rankings
As before, we can calculate the rankings for each treatment and illustrate this in a rankogram. This can be done using the following code.
mtc_52_rank <- rank.probability(mtc_52, preferredDirection = -1)
rank_pos <- c()
rank_prob <- c()
for(i in 1:length(network_52$treatments$id)){
rank_pos <- c(rank_pos, rep(i, length(network_52$treatments$id)))
rank_prob <- c(rank_prob, mtc_52_rank[,i])
}
df_rank_plot_52 <- data.frame(Probability = rank_prob,
treatment_old = network_52$treatments$id,
Treatment = network_52$treatments$description,
Rank = factor(rank_pos))
df_rank_plot_52$Treatment <- factor(df_rank_plot_52$Treatment,
levels = df_rank_plot_52$Treatment[c(2, 1, 3:8)])
rank_52_plot <- ggplot(df_rank_plot_52, aes(x=Treatment, y=Probability, fill=Rank))+
geom_bar(stat="identity", color="black", position=position_dodge())+
scale_fill_grey(start=0, end=1) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust=1))+
ylab("")
rank_52_plot
The results of this rankogram closely match those of the 24-26 week NMA. Again we see that Empagliflozin treatments are likely to be the most effective (25mg) and second most effective (10mg) treatments at reducing HbA1c. Likewise we also see that the Placebo and Metformin (2000mg) treatments are likely to be the least and second least effective treatments.
SUCRA values
As before, we can also calculate SUCRA values to compliment our rankogram. This can be conducted using the following code.
#calculate the SUCRA values for each of the treatments using the rankings
mtc_52_sucra <- dmetar::sucra(mtc_52_rank, lower.is.better = FALSE)
#extract the SUCRA data in an easily plotable format
x_52_sucra <- plot(mtc_52_sucra)
#generate a dataframe to store the SUCRA data
df_sucra_plot_52 <- data.frame(SUCRA = NA,
treatment_old = network_52$treatments$id,
Treatment = network_52$treatments$description,
colours = vertex_colours_52)
#extract the SUCRA values for each treatment
for(i in 1:dim(df_sucra_plot_52)[1]){
df_sucra_plot_52$SUCRA[i] <- subset(x_52_sucra$data,
Treatment == df_sucra_plot_52$treatment_old[i])$SUCRA
}
#change the Treatment to a factor variable
#this is so that we can order the different treatments in the plot in the order we want
#here we are ordering treatments alphabetically and then by increasing dosage
df_sucra_plot_52$Treatment <- factor(df_sucra_plot_52$Treatment,
levels = df_sucra_plot_52$Treatment[c(2, 1, 3:8)])
#generate the SUCRA plot using ggplot
sucra_52_plot <- ggplot(df_sucra_plot_52, aes(x=Treatment, y=SUCRA, fill=Treatment))+
geom_bar(stat="identity", color="black")+
scale_fill_manual(values=df_sucra_plot_52$colours)+
theme_minimal() +
theme(legend.position = "None") +
theme(axis.text.x = element_text(angle = 45, hjust=1))
sucra_52_plot
The results of this figure closely align with the conclusions of the rankogram, with the same treatments being identified as being either the most (Empagliflozin 25mg, 10mg) or least (Placebo, Metformin 2000mg) effective at reducing HbA1c.
Forest plots
In their paper, Langford et al. report forest plots showcasing how the various SGLT inhibitor treatments compare either to the placebo treatment or the Metformin (2000mg) treatment at 52 weeks. As such, we can likewise generate these forest for our model and compare our results to those reported by Langford et al. (2020).
Comparison of treatments to placebo
As we did previously, we'll start by generating the forest plot showing the effect sizes (mean difference) in CFB HbA1c (%) for each treatment in comparison to the placebo treatment. This can be done using the following code.
#extract the summary information where we calculate the mean difference in HbAc1 for each treatment relative to the placebo
summary_placebo_52 <- summary(relative.effect(mtc_52, t1 = "Placebo"))
#generate a dataframe to store the data
df_forest_placebo_52 <- data.frame(Treatment = network_52$treatments$description,
colours = vertex_colours_52,
Mean = NA,
CrI_Lower = NA,
CrI_Upper = NA)
#remove any data corresponding to the placebo treatment from our dataframe
df_forest_placebo_52 <- subset(df_forest_placebo_52, Treatment != "Placebo")
#calculate the dimensions of our dataframe - this is to shorten the code later
dim_df <- dim(df_forest_placebo_52)[1]
#Extract mean values for the estimate of the effect of each treatment
df_forest_placebo_52$Mean <- summary_placebo_52$summaries$statistics[1:dim_df]
#Extract lower bounds for the credible intervals of each treatment effect
df_forest_placebo_52$CrI_Lower <- summary_placebo_52$summaries$quantiles[1:dim_df]
#Extract upper bounds for the credible intervals of each treatment effect
df_forest_placebo_52$CrI_Upper <- summary_placebo_52$summaries$quantiles[((4*dim_df)+1):(5*dim_df)]
#change the Treatment to a factor variable
#this is so that we can order the different treatments in the plot in the order we want
#here we are ordering treatments in same order as Langford et al. (2020)
df_forest_placebo_52$Treatment <- factor(df_forest_placebo_52$Treatment,
levels = df_forest_placebo_52$Treatment[c(5, 2, 1, 6, 7, 3, 4)])
#generate our forest plot for placebo comparisons
forest_placebo_52_plot <- ggplot(df_forest_placebo_52, aes(x=Mean,
y=Treatment,
fill = Treatment)) +
geom_errorbar(aes(y = Treatment, xmin=CrI_Lower, xmax=CrI_Upper, width = 0)) +
geom_point(size = 5, shape = 21, color = "black") +
scale_fill_manual(values = df_forest_placebo_52$colours[c(5, 2, 1, 6, 7, 3, 4)]) +
theme_classic() +
theme(legend.position = "none") +
geom_vline(colour = "gray50", xintercept = 0, linetype="dashed") +
xlab("Mean difference in CFB HbA1c relative to placebo (%)") +
scale_y_discrete(limits=rev) +
ggtitle("B) Week 52")
forest_placebo_52_plot
This forest plot shows the estimated effect sizes, for the 52 week data, for the mean difference (%) in CFB HbA1c for each treatment relative to the placebo. Here, each effect size estimate has been determined using our fixed effects Bayesian-NMA, with the error bars representing the 95% credible intervals for each effect size. Our analysis suggests that all treatments (except for Metformin 2000mg), at 52 weeks, led to significant reductions in HbA1c compared to the placebo treatment. Firslty, our results suggest that there is no significant difference between the Metformin (2000mg) and placebo groups, a finding which contrasts the results for the 24-26 week data. Interestingly, the magnitude of the effects for each treatment at 52 weeks appear to be less than those at 24-26 weeks. Additionally, there is substantial overlap between the 95% credible intervals for each treatment (aside from Metformin 2000mg).
Again, the figure we have generated here can be directly compared to Figure 3a from Langford et al. (2020). By comparing the figures, we can see the results of our analysis match those reported by Langford et al. (2020) with any minor differences likely to be due to the stochastic nature of the analysis or differences between R and WINBUGS.
Comparison of treatments to Metformin (2000mg)
As before, we can also generate a forest plot showing the effect sizes (mean difference) in CFB HbA1c (%) for each SGLT inhibitor treatment and the placebo in comparison to the Metformin (2000mg) treatment. This can be done using the following code.
#extract the summary information where we calculate the mean difference in HbAc1 for each treatment relative to the Metformin treatment
summary_metformin_52 <- summary(relative.effect(mtc_52, t1 = "Metformin_2000"))
#generate a dataframe to store the data
df_forest_metformin_52 <- data.frame(Treatment = network_52$treatments$description,
colours = vertex_colours_52,
Mean = NA,
CrI_Lower = NA,
CrI_Upper = NA)
#remove any data corresponding to the Metformin treatment from our dataframe
df_forest_metformin_52 <- subset(df_forest_metformin_52, Treatment != "Metformin 2000mg")
#calculate the dimensions of our dataframe - this is to shorten the code later
dim_df <- dim(df_forest_metformin_52)[1]
#Extract mean values for the estimate of the effect of each treatment
df_forest_metformin_52$Mean <- summary_metformin_52$summaries$statistics[1:dim_df]
#Extract lower bounds for the credible intervals of each treatment effect
df_forest_metformin_52$CrI_Lower <- summary_metformin_52$summaries$quantiles[1:dim_df]
#Extract upper bounds for the credible intervals of each treatment effect
df_forest_metformin_52$CrI_Upper <- summary_metformin_52$summaries$quantiles[((4*dim_df)+1):(5*dim_df)]
#change the Treatment to a factor variable
#this is so that we can order the different treatments in the plot in the order we want
#here we are ordering treatments in same order as Langford et al. (2020)
df_forest_metformin_52$Treatment <- factor(df_forest_metformin_52$Treatment,
levels = df_forest_metformin_52$Treatment[c(5, 2, 1, 6, 7, 3, 4)])
#generate our forest plot for metformin comparisons
forest_metformin_52_plot <- ggplot(df_forest_metformin_52, aes(x=Mean,
y=Treatment,
fill = Treatment)) +
geom_errorbar(aes(y = Treatment, xmin=CrI_Lower, xmax=CrI_Upper, width = 0)) +
geom_point(size = 5, shape = 21, color = "black") +
scale_fill_manual(values = df_forest_metformin_52$colours[c(5, 2, 1, 6, 7, 3, 4)]) +
theme_classic() +
theme(legend.position = "none") +
geom_vline(colour = "gray50", xintercept = 0, linetype="dashed") +
xlab("Mean difference in CFB HbA1c relative to Metformin 2000mg (%)") +
scale_y_discrete(limits=rev) +
ggtitle("B) Week 52")
forest_metformin_52_plot
This forest plot shows the estimated effect sizes for the mean difference (%) in CFB HbA1c for each SGLT inhibitor treatment and the placebo relative to Metformin (2000mg) at 52 weeks. Here, each effect size estimate has been determined using our fixed effects Bayesian-NMA, with the error bars representing the 95% credible intervals for each effect size. Our analysis suggests that all treatments (bar one) led to significant reductions in HbA1c compared to the Metformin (2000mg) treatment. As with the 24-26 week data, the single exception to this was that there was no significant difference between the placebo treatment and Metformin (2000mg) treatment illustrated by 95% credible interval for the placebo treatment overlapping zero. The magnitude of any effect substantially varied between treatments, although there is substantial overlap between the credible intervals of each treatment.
The figure we have generated here can be directly compared to Figure 4a from Langford et al. (2020). By comparing the figures, we can see the results of our analysis match those reported by Langford et al. (2020) with any minor differences likely to be due to the stochastic nature of the analysis or differences between R and WINBUGS.
Importantly, all comparisons with the Metformin (2000mg) treatment (aside from the placebo) are indirect (in line with the 24-26 week analysis). As noted by Langford et al. (2020), all of these effect sizes are derived solely from indirect comparisons to the Metformin (2000mg) treatment and hence potentially liable to bias.
Combined HbA1c Forest Plots
Figures 3a and 4a in Langford et al. (2020) combine the forest plots of the 24-26 week and 52 week Bayesian-NMAs. Accordingly, this can be done for our analysis using the below code.
Comparison of treatments to the placebo
The code to combine the 24-26 week and 52 week forest plots which compare each treatment to the placebo is as follows.
#standardise x-axis and remove the x-axis title
forest_placebo_24_plot <- forest_placebo_24_plot +
xlim(c(-0.7, 0.3)) +
xlab("")
#standardise x-axis
forest_placebo_52_plot <- forest_placebo_52_plot +
xlim(c(-0.7, 0.3))
#generate the combined plot
combined_forest_placebo_plot <- (forest_placebo_24_plot / forest_placebo_52_plot)
combined_forest_placebo_plot
As before, the results shown in this figure match those shown in Figure 3a from Langford et al. (2020).
Comparison of treatments to the metformin
The code to combine the 24-26 week and 52 week forest plots which compare each treatment to the Metformin (2000mg) treatment is as follows.
#standardise x-axis and remove the x-axis title
forest_metformin_24_plot <- forest_metformin_24_plot +
xlim(c(-0.6, 0.6)) +
xlab("")
#standardise x-axis
forest_metformin_52_plot <- forest_metformin_52_plot +
xlim(c(-0.6, 0.6))
#generate the combined plot
combined_forest_metformin_plot <- (forest_metformin_24_plot / forest_metformin_52_plot)
combined_forest_metformin_plot
As before, the results shown in this figure match those shown in Figure 4a from Langford et al. (2020).
Concluding remarks
In this tutorial, we have provided an introduction to Bayesian-NMAs using a real-world dataset from a recently published paper. The dataset in question considered how various SGLT inhibitor treatments (both different medications and dosages) affected various metrics in patients with T1DM. The original analysis was conducted in WINBUGS, but here we have been able to replicate the analysis of Langford et al. (2020) in R. Importantly, the results we obtained here almost exactly match those in the original paper, with any differences likely to be caused by mechanisitic differences between WINBUGS and R or the stochastic nature of Bayesian-NMAs.
Here, we have solely focussed on conducting a Bayesian-NMA for data (at 24-26 and 52 weeks) on HbA1c (a measurement of glycated haemoglobin). Based on our findings, it appears as though SGLT inhibitors may reduce HbA1c to a greater extent than the placebo or Metformin treatments. Indeed, higher doses of Empagliflozin (10mg / 25mg) appear to be particularly effective at reducing HbA1c (relative to either the placebo or Metformin treatments) at both time points. However, as noted above and by Langford et al. (2020) only indirect comparisons can be made between the Metformin and SGLT inhibitor treatments meaning that these comparisons are potentially subjected to bias.
Within this tutorial we have solely considered one metric (CFB in HbA1c). However, Langford et al. (2020) consider several other primary metrics (weight, total daily insulin dose, and systolic blood pressure). Indeed, each of the considered SGLT inhibitor treatments is unlikely to affect a single metric (e.g., HbA1c) in isolation, therefore it is imperative that the effect of each treatment on multiple metrics is considered. Additionally, consideration of various safety outcomes for each treatment is also important. For instance, a SGLT inhibitor treatment may effectively reduce a metric of interest but also induce unwanted or potentially dangerous side effects. Langford et al. (2020) additionally analyse the safety outcomes for a range of adverse events including hypoglycaemia and diabetic ketoacidosis. Indeed, in their analysis of safety outcomes, Langford et al. (2020) highlight that (depsite its apparent efficacy in reducing HbA1c) Sotagliflozin use may be associated with an increased prevalence of diabetic ketoacidosis, a finding supported by previous analyses.
From a meta-analytical perspective, there are additional statistical attributes that could be explored. For instance, inconsistency in an NMA could be considered. Here, the relatively simple network structure with few links between treatments means that there are no comparisons which can be assessed for inconsistency. Indeed, any attempts to model inconsistency using the mtc.nodesplit function in gemtc fail (for our NMA) for this reason. However, if there is the potential for inconsistency to occur in a network then analysis of inconsistency should be conducted.
Overall, this tutorial has attempted to provide an introduction to conducting Bayesian-NMAs in R using a real-world medical dataset. Hopefully this has been informative and outlined the various steps and approaches that need to be taken to conduct a statistically robust NMA in R!
References
Langford, B. E., Evans, M., Haskins-Coulter, T., O'Connor, M., Cant, H. E., Eddowes, L. A., ... & Tank, A. (2020). Systematic literature review and network meta-analysis of sodium-glucose co-transporter inhibitors vs metformin as add-on to insulin in type 1 diabetes. Diabetes, Obesity and Metabolism, 22(1), 39-50.