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 van Valkenhoef, G et al. (2012). 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 (Plummer 2017). 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 (Lunn et al. 2012).
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 (Lu and Ades 2006). 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.
