In TasnimHamza/crosnma:
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(crosnma)
library(rjags) # only for bias.method='adjust2' - to use the function load.module()
load.module('mix') # only for bias.method='adjust2' - to call the bimodal normal models to JAGS
Introduction
In network meta-analysis we synthesize all relevant available evidence about health outcomes from competing treatments. That evidence might come from different study designs and in different formats: from non-randomized studies (NRS) or randomized controlled trials (RCT) as individual participant data (IPD) or as aggregate data (AD). We set up the package crosnma to synthesize all available evidence.
This document demonstrates how to use crosnma to synthesize cross-design evidence and cross-format data via Bayesian network meta-analysis and meta-regression (NMA and NMR). All models are implemented in JAGS [@plummer_jags].
We describe the workflow within the package using a worked example from a network meta-analysis of studies for treatments in relapsing remitting multiple sclerosis (RRMS). The primary outcome is the occurrence of relapses in two years (binary outcome). In the analysis, the relative effect will be the odds ratio (OR). The aim is to compare the efficacy of four treatments using the data from 6 different studies in different formats and different designs.
The synthesis models
We first introduce the model that synthesizes studies with individual-level (IPD) or/and aggregate data (AD) ignoring their design (naïve synthesis). Then, we present three possible models that account for the different study designs. In the table below we set the notation that will be used in the description of the four models.
| Notation | Description| Argument in crosnma.model()|
|:-------------------|:-----------|:------|
|$i=1, ..., np_j$ | participant id| |
|$j=1, ..., ns$ | study id| |
|$k=1, ..., K$ | treatment index| |
|$ns_{IPD}, ns_{AD}, ns_{RCT}, ns_{NRS}$| the number of studies. The index refers to the design or format of the study| |
|$y_{ijk}$ | binary outcome (0/1)| outcome |
|$p_{ijk}$ | probability of the event to occur| |
|$r_{jk}$ | the number of events per arm| outcome |
|$n_{jk}$ | the sample size per arm| n |
|$b$ |the study-specific reference||
|$u_{j}$ | The treatment effect of the study-specific reference | |
|$\delta_{jk}$|log(OR) of treatment k relative to $b$||
|$x_{ijk}$|the covariate|covariate|
|$\bar{x}{j}$|the mean covariate for study $j$||
|$d{Ak}$| the basic parameters. Here, $d_{AA}=0$ when A set as the reference in the network|use reference to assign the reference treatment|
|$z_j$| study characteristics to estimate the bias probability $\pi_j$| bias.covariate |
|$w$| common inflation factor of variance for the NRS estimates | the element var.infl in run.nrs|
|$\zeta$| common mean shift of the NRS estimates | the element mean.shift in run.nrs|
If the reference in the network ($A$) is available on the study, it is assigned automatically to that reference. If not, it is assigned to the first alphabetically ordered treatment on the study.
Naive synthesis
We synthesize the evidence from RCT and NRS without acknowledging the differences between them. We combine the IPD data from RCT and NRS in one model and we do the same in another model with the AD information. Then, we combine the estimates from both parts as described in Section 2.5.
model IPD only
$$
y_{ijk} \sim Bernoulli(p_{ijk})
$$
\begin{equation}
logit(p_{ijk}) =
\begin{cases}
u_j +\beta_{0,j} x_{ijk} & \text{if $k=b$}\
u_j +\delta_{jk} + (\beta_{0,j}+\beta^w_{1,jk}) x_{ijk}+
(\beta^B_{1,jk}-\beta^w_{1,jk}) \bar{x}_{.j} & \text{if $k≠b$}
\end{cases}
\end{equation}
model AD only
$$
r_{jk} \sim Binomial(p_{jk},n_{jk})
$$
\begin{equation}
logit(p_{jk}) =
\begin{cases}
u_j & \text{if $k=b$}\
u_j +\delta_{jk} +\beta^B_{1,jk} \bar{x}_{j} & \text{if $k≠b$}
\end{cases}
\end{equation}
Using non-randomized studies (NRS) as a prior
First we estimate the relative treatment effects using only the NRS (use run.nrs in crosnma.model() to control this process). Then we use the NRS estimates ($\tilde{d}^{NRS}{Ak},V^{NRS}{Ak}$) as a prior information for the basic parameters of RCT data, $d_{Ak} \sim N(\tilde{d}^{NRS}{Ak},V^{NRS}{Ak})$. To control the NRS influence in the RCT estimates, we can either inflate the prior variance by dividing it by a common inflation factor $w$ (the inflated variances are $V^{NRS}_{Ak}/w$) or shift the NRS means by $\zeta$.
Bias-adjusted model 1
In this and the next model we incorporate the judgments about the study risk of bias (RoB) by extending the method introduced by @dias_2010. In model 1, we either multiply the bias term $\gamma_{1,j}^{R_j}$ or add $\gamma_{2,j} R_j$ to the relative treatment effect on both the AD and IPD parts of the model. We might include only one of the aforementioned bias terms and in this case a single coefficient ($\gamma_j$) will be estimated.
model IPD only
\begin{equation}
logit(p_{ijk}) =
\begin{cases}
u_j +\beta_{0,j} x_{ijk} & \text{if $k=b$}\
u_j +\delta_{jk}*\gamma_{1,j}^{R_j}+\gamma_{2,j} R_j+ (\beta_{0,j}+\beta^w_{1,jk}) x_{ijk}+
(\beta^B_{1,jk}-\beta^w_{1,jk}) \bar{x}_{.j} & \text{if $k≠b$}
\end{cases}
\end{equation}
model AD only
\begin{equation}
logit(p_{jk}) =
\begin{cases}
u_j & \text{if $k=b$}\
u_j +\delta_{jk} *\gamma_{1,j}^{R_j}+\gamma_{2,j} R_j+
\beta^B_{1,jk} \bar{x}_{j} & \text{if $k≠b$}
\end{cases}
\end{equation}
where the bias indicator $R_j$ follows the following distribution
$$
R_j \sim Bernoulli(\pi_j)
$$
The bias probabilities $\pi_j$ are study-specific and can be estimated in two different ways. They are either given informative beta priors (${Beta(a,b)}$) that are set according to the risk of bias. The default beta priors are as follows: high bias RCT pi.high.rct='dbeta(10,1)', low bias RCT pi.low.rct='dbeta(1,10)', high bias NRS pi.high.nrs='dbeta(30,1)' and low bias NRS pi.low.nrs='dbeta(1,30)'. The ratio ${a/b}$ controls the skewness of the beta distribution. The closer to 1 the ratio a/b, the more the mean of probability of bias gets closer to 1 and the study acquires 'major' bias adjustment. Alternatively, we can use the study characteristics $z_j$ to predict $\pi_j$ through a logistic transformation (internally coded).
Bias-adjusted model 2
Another way to incorporate the RoB of the study is by replacing $\delta_{jk}$ by a "bias-adjusted" relative treatment effect $\theta_{jk}$. Then $\theta_{jk}$ is modeled with a bimodal normal distribution as described in Section 2.5. For more details see @verde_2020.
model IPD only
\begin{equation}
logit(p_{ijk}) =
\begin{cases}
u_j +\beta_{0,j} x_{ijk} & \text{if $k=b$}\
u_j +\theta_{jk} + (\beta_{0,j}+\beta^w_{1,jk}) x_{ijk}+
(\beta^B_{1,jk}-\beta^w_{1,jk}) \bar{x}_{j} & \text{if $k≠b$}
\end{cases}
\end{equation}
model AD only
\begin{equation}
logit(p_{jk}) =
\begin{cases}
u_j & \text{if $k=b$}\
u_j +\theta_{jk} +\beta^B_{1,jk} \bar{x}_{j} & \text{if $k≠b$}
\end{cases}
\end{equation}
where the bias-adjusted relative treatment effect ($\theta_{jk}$) are modeled as follows
$$
\theta_{jk} \sim \pi_j N(d_{Ak}-d_{Ab}, \tau_\delta^2) + (1-\pi_j)N(d_{Ak}-d_{Ab}+\gamma_j, \tau_\delta^2+\tau_\gamma^2)
$$
Assumptions about the model parameters
The table below summarizes the different assumptions implemented in the package about combining the parameters in the models described above.
| Parameter | Assumptions| Argument in crosnma.model()|
|:----------------------|:-----------|:------|
|relative treatment effect ($\delta_{jk}$)| Random-effects: $\delta_{jk} \sim N(d_{Ak}-d_{Ab}, \tau_\delta^2)$| trt.effect='random' |
| |Common-effect: $\delta_{jk}=d_{Ak}-d_{Ab}$| trt.effect='common'|
|Covariate effect $\beta_{0, j}$ | Independent effects: $\beta_{0, j} \sim N(0, 10^2)$| reg0.effect='independent' |
| |Random-effects: $\beta_{0,j} \sim N(B_0, \tau_{\beta_0})$| reg0.effect='random'|
|Within-study covariate-treatment interaction ($\beta_{1, jk}^W$)| Random-effects: $\beta_{1, jk}^W \sim N(B_{1, Ak}^W-B_{1, Ab}^W, \tau_{\beta_1^W})$| regw.effect='random' |
| |Common-effect: $\beta_{1,jk}^W = B_{1, Ak}^W-B_{1, Ab}^W$| regw.effect='common'|
|Between-study covariate-treatment interaction ($\beta_{1, jk}^B$)| Random-effects: $\beta_{1,jk}^B \sim N(B_{1, Ak}^B-B_{1, Ab}^B, \tau_{\beta_1^B})$| regb.effect='random' |
| |Common-effect: $\beta_{1,jk}^B = B_{1, Ak}^B-B_{1, Ab}^B$| regb.effect='common'|
|Bias effect ($\gamma_{{1,2},j}$)| Random-effects: $\gamma_{{1,2},j} \sim N(\Gamma_{{1,2}}, \tau_{\gamma})$| bias.effect='random' |
| |Common-effect: $\gamma_{{1,2},j}=\Gamma_{{1,2}}$| bias.effect='common'|
|Bias probability ($\pi_j$)| $\pi_j \sim Beta(a,b)$| |
|| $\pi_j = a+bz_j$| |
Synthesis of studies comparing drugs for relapsing-remitting multiple sclerosis
Description of the data
The data we use are fictitious but have been developed to resample to real RCTs with IPD and aggregate data included in @Tramacere15. The studies provide either aggregate data std.data (2 RCTs) or as individual participant data prt.data (3 RCTs and 1 cohort study). Both datasets compare in total four drugs which are anonymized.
The prt.data contains 2950 rows, each row refers to a participant in the study. We display the first few rows of the data set:
head(prt.data)
For each participant, we have information for the outcome relapse (0=no, 1=yes), the treatment label trt, the age (in years) and sex (0 = Female, 1 = Male) of the participant. The following columns are set on study-level (it is repeated for each participant in each study): the study id, the design of the study (needs to be either rct or nrs), the risk of bias on each study (can be set as low, high or unclear) and the year of publication .
The aggregate data has the standard format for meta-analysis
head(std.data)
Analysis
The network should be checked for its connectivity before running the analysis. This is a vital step as the model will run even if the network is not connected.
netplot(prt.data,std.data)
In the following table, we summarize the number of studies from each design and each data format:
knitr::kable(ns.tab(prt.data,std.data))
There are two steps to run the NMA/NMR model. The first step is to create a JAGS model using crosnma.model() which produces the JAGS code and the data. In the second step, the output of that function will be used in crosnma.run() to run the analysis through JAGS.
Naïve network meta-analysis
We start by indicating the names of the datasets on participant-level (prt.data) and aggregate data (std.data). Then, the name of the variables on each dataset needs to be given respectively in prt.data and std.data. Next, the reference treatment needs to be assigned (we set it to drug A). By choosing trt.effect=random, we are assigning a normal distribution to each relative treatment effect to allow the synthesis across studies, see the table in Section 2.1. Finally, the different designs; RCT and NRS are combined with the information taken at face-value; method.bias = 'naive'.
Optionally, we can specify a prior to the common heterogeneity of the treatment effect across studies. We indicate that distribution in the argument prior as tau.trt='dunif(0,3)', see below.
# jags model: code+data
mod1 <- crosnma.model(prt.data=prt.data,
std.data=std.data,
trt=c('trt','trt'),
study=c('study','study'),
outcome=c('outcome','outcome'),
n='n',
design=c('design','design'),
trt.effect='random',
reference='A',
method.bias = 'naive',
#---------- assign a prior ----------
prior=list(tau.trt='dunif(0,3)')
)
Next, we fit the NMA model using crosnma.run()which requires us to set the number of adaptations, iterations, thinning and chains.
# run jags
jagsfit1 <- crosnma.run(model=mod1,
n.adapt = 50,
n.iter=500,
n.burnin = 200,
thin=1,
n.chains=2)
We summarize the estimated parameters in the following table.
knitr::kable(summary(jagsfit1,expo=F))
The estimated OR of B vs A can be obtained as exp(d.B) and similarly for exp(d.C) and exp(d.D) are the ORs of C and D relative to A, respectively. The value of tau refers to the estimates of the heterogeneity standard deviation in the relative treatment effects across studies.
We need also to assess the convergence of the MCMC chains either by visually inspect the trace plot or checking the Gelman and Rubin statistic, Rhat (it should be approximately 1) in the table above.
coda::traceplot(jagsfit1$samples)
Naïve network meta-regression
In this part, we run NMR model by adding age as a covariate from both datasets. We set a list of elements covariate=list(c('age'),c('age')) representing the names of the covariate in prt.data and std.data, respectively.
# jags model: code+data
mod2 <- crosnma.model(prt.data=prt.data,
std.data=std.data,
trt=c('trt','trt'),
study=c('study','study'),
outcome=c('outcome','outcome'),
n='n',
design=c('design','design'),
reference='A',
trt.effect='random',
#---------- meta-regression ----------
covariate = list(c('age'),c('age')),
split.regcoef = F,
#---------- bias adjustment ----------
method.bias='naive'
)
The MCMC is run under the same set up as in the network meta-analysis.
# run jags
jagsfit2 <- crosnma.run(model=mod2,
n.adapt = 50,
n.iter=500,
n.burnin = 200,
thin=1,
n.chains=2)
and the output table is presented below
knitr::kable(summary(jagsfit2,expo=FALSE))
Now, we additionally estimate b_1 which indicates the mean effect of age and tau.b_1 which refers to the heterogeneity standard deviation in the effect of age across studies. Here, we obtain a single estimate because we choose to not split the within- and between-study age coefficients $(\beta^W_{1,jk} = \beta^B_{1,jk}=\beta_{1,jk})$ to improve the convergence of MCMC.
Again, we check convergence with trace plots
coda::traceplot(jagsfit2$samples)
Using non-randomized studies as a prior in network meta-regression
To run NMA with a prior from NRS, two additional arguments are needed: we indicate using NRS as a prior by setting method.bias='prior'. That means that the model runs internally NMA with only NRS data which are then used to construct informative priors. This requires defining MCMC settings (the number of adaptations, iterations, burn-ins, thinning and chains) in the argument run.nrs.
In this method, the prior for the basic parameters is set to a normal distribution. For basic parameters not examined in the NRS, the code sets a minimally informative prior d~dnorm(0, 1e-4). To account for possible bias, the means of the distribution can be shifted by mean.shift to reflect the potential bias in NRS and/or the variance can be inflated by var.infl to control the influence of NRS on the final estimation. Both should be provided in run.nrs.
# jags model: code+data
mod3 <- crosnma.model(prt.data=prt.data,
std.data=std.data,
trt=c('trt','trt'),
study=c('study','study'),
outcome=c('outcome','outcome'),
n='n',
design=c('design','design'),
reference='A',
trt.effect='random',
#---------- meta-regression ----------
covariate = list(c('age'),c('age')),
split.regcoef = F,
#---------- bias adjustment ----------
method.bias='prior',
run.nrs=list(var.infl=0.6,
n.adapt = 500,
n.iter=10000,
n.burnin = 4000,
thin=1,
n.chains=2))
# run jags
jagsfit3 <- crosnma.run(model=mod3,
n.adapt = 50,
n.iter=500,
n.burnin = 200,
thin=1,
n.chains=2)
knitr::kable(summary(jagsfit3,expo=F))
coda::traceplot(jagsfit3$samples)
Bias-adjusted model 1
In this part, the overall relative treatment effects are estimated from both NRS and RCT with adjustment to study-specific bias.
To fit the model, we set method.bias='adjust1' and we need to provide the name of the bias variable bias=c('bias','bias') in prt.data and std.data, respectively. By default, the effect of bias is assumed to be additive bias.type='add' and equal across studies bias.effect='common'. We also use the year of study publication to estimate the study-probability of bias, bias.covariate = c('year','year').
# jags model: code+data
mod4 <- crosnma.model(prt.data=prt.data,
std.data=std.data,
trt=c('trt','trt'),
study=c('study','study'),
outcome=c('outcome','outcome'),
n='n',
design=c('design','design'),
reference='A',
trt.effect='random',
#---------- bias adjustment ----------
method.bias='adjust1',
bias=c('bias','bias'),
bias.type='add',
bias.effect='common',
bias.covariate = c('year','year')
)
# run jags
jagsfit4 <- crosnma.run(model=mod4,
n.adapt = 50,
n.iter=500,
n.burnin = 200,
thin=1,
n.chains=2)
The results are presented below
knitr::kable(summary(jagsfit4,expo=F))
The parameter g refers to the mean bias effect, common for all studies.
The trace plots are shown below
coda::traceplot(jagsfit4$samples)
Bias-adjusted model 2
The arguments for method.bias='adjust2' are similar to the ones used before in method.bias='adjust1'.
# jags model: code+data
mod5 <- crosnma.model(prt.data=prt.data,
std.data=std.data,
trt=c('trt','trt'),
study=c('study','study'),
outcome=c('outcome','outcome'),
n='n',
design=c('design','design'),
reference='A',
trt.effect='random',
#---------- bias adjustment ----------
method.bias='adjust2',
bias=c('bias','bias'),
bias.type='add',
bias.effect='common',
bias.covariate = c('year','year')
)
# run jags
jagsfit5 <- crosnma.run(model=mod5,
n.adapt = 50,
n.iter=500,
n.burnin = 200,
thin=1,
n.chains=2)
knitr::kable(summary(jagsfit5,expo=F))
coda::traceplot(jagsfit5$samples)
References
TasnimHamza/crosnma documentation built on Dec. 18, 2021, 4:05 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(crosnma) library(rjags) # only for bias.method='adjust2' - to use the function load.module() load.module('mix') # only for bias.method='adjust2' - to call the bimodal normal models to JAGS
Introduction
In network meta-analysis we synthesize all relevant available evidence about health outcomes from competing treatments. That evidence might come from different study designs and in different formats: from non-randomized studies (NRS) or randomized controlled trials (RCT) as individual participant data (IPD) or as aggregate data (AD). We set up the package crosnma to synthesize all available evidence.
This document demonstrates how to use crosnma to synthesize cross-design evidence and cross-format data via Bayesian network meta-analysis and meta-regression (NMA and NMR). All models are implemented in JAGS [@plummer_jags].
We describe the workflow within the package using a worked example from a network meta-analysis of studies for treatments in relapsing remitting multiple sclerosis (RRMS). The primary outcome is the occurrence of relapses in two years (binary outcome). In the analysis, the relative effect will be the odds ratio (OR). The aim is to compare the efficacy of four treatments using the data from 6 different studies in different formats and different designs.
The synthesis models
We first introduce the model that synthesizes studies with individual-level (IPD) or/and aggregate data (AD) ignoring their design (naïve synthesis). Then, we present three possible models that account for the different study designs. In the table below we set the notation that will be used in the description of the four models.
| Notation | Description| Argument in crosnma.model()|
|:-------------------|:-----------|:------|
|$i=1, ..., np_j$ | participant id| |
|$j=1, ..., ns$ | study id| |
|$k=1, ..., K$ | treatment index| |
|$ns_{IPD}, ns_{AD}, ns_{RCT}, ns_{NRS}$| the number of studies. The index refers to the design or format of the study| |
|$y_{ijk}$ | binary outcome (0/1)| outcome |
|$p_{ijk}$ | probability of the event to occur| |
|$r_{jk}$ | the number of events per arm| outcome |
|$n_{jk}$ | the sample size per arm| n |
|$b$ |the study-specific reference||
|$u_{j}$ | The treatment effect of the study-specific reference | |
|$\delta_{jk}$|log(OR) of treatment k relative to $b$||
|$x_{ijk}$|the covariate|covariate|
|$\bar{x}{j}$|the mean covariate for study $j$||
|$d{Ak}$| the basic parameters. Here, $d_{AA}=0$ when A set as the reference in the network|use reference to assign the reference treatment|
|$z_j$| study characteristics to estimate the bias probability $\pi_j$| bias.covariate |
|$w$| common inflation factor of variance for the NRS estimates | the element var.infl in run.nrs|
|$\zeta$| common mean shift of the NRS estimates | the element mean.shift in run.nrs|
If the reference in the network ($A$) is available on the study, it is assigned automatically to that reference. If not, it is assigned to the first alphabetically ordered treatment on the study.
Naive synthesis
We synthesize the evidence from RCT and NRS without acknowledging the differences between them. We combine the IPD data from RCT and NRS in one model and we do the same in another model with the AD information. Then, we combine the estimates from both parts as described in Section 2.5.
model IPD only
$$
y_{ijk} \sim Bernoulli(p_{ijk})
$$
\begin{equation}
logit(p_{ijk}) =
\begin{cases}
u_j +\beta_{0,j} x_{ijk} & \text{if $k=b$}\
u_j +\delta_{jk} + (\beta_{0,j}+\beta^w_{1,jk}) x_{ijk}+
(\beta^B_{1,jk}-\beta^w_{1,jk}) \bar{x}_{.j} & \text{if $k≠b$}
\end{cases}
\end{equation}
model AD only
$$
r_{jk} \sim Binomial(p_{jk},n_{jk})
$$
\begin{equation}
logit(p_{jk}) =
\begin{cases}
u_j & \text{if $k=b$}\
u_j +\delta_{jk} +\beta^B_{1,jk} \bar{x}_{j} & \text{if $k≠b$}
\end{cases}
\end{equation}
Using non-randomized studies (NRS) as a prior
First we estimate the relative treatment effects using only the NRS (use run.nrs in crosnma.model() to control this process). Then we use the NRS estimates ($\tilde{d}^{NRS}{Ak},V^{NRS}{Ak}$) as a prior information for the basic parameters of RCT data, $d_{Ak} \sim N(\tilde{d}^{NRS}{Ak},V^{NRS}{Ak})$. To control the NRS influence in the RCT estimates, we can either inflate the prior variance by dividing it by a common inflation factor $w$ (the inflated variances are $V^{NRS}_{Ak}/w$) or shift the NRS means by $\zeta$.
Bias-adjusted model 1
In this and the next model we incorporate the judgments about the study risk of bias (RoB) by extending the method introduced by @dias_2010. In model 1, we either multiply the bias term $\gamma_{1,j}^{R_j}$ or add $\gamma_{2,j} R_j$ to the relative treatment effect on both the AD and IPD parts of the model. We might include only one of the aforementioned bias terms and in this case a single coefficient ($\gamma_j$) will be estimated.
model IPD only
\begin{equation}
logit(p_{ijk}) =
\begin{cases}
u_j +\beta_{0,j} x_{ijk} & \text{if $k=b$}\
u_j +\delta_{jk}*\gamma_{1,j}^{R_j}+\gamma_{2,j} R_j+ (\beta_{0,j}+\beta^w_{1,jk}) x_{ijk}+
(\beta^B_{1,jk}-\beta^w_{1,jk}) \bar{x}_{.j} & \text{if $k≠b$}
\end{cases}
\end{equation}
model AD only
\begin{equation}
logit(p_{jk}) =
\begin{cases}
u_j & \text{if $k=b$}\
u_j +\delta_{jk} *\gamma_{1,j}^{R_j}+\gamma_{2,j} R_j+
\beta^B_{1,jk} \bar{x}_{j} & \text{if $k≠b$}
\end{cases}
\end{equation}
where the bias indicator $R_j$ follows the following distribution
$$ R_j \sim Bernoulli(\pi_j) $$
The bias probabilities $\pi_j$ are study-specific and can be estimated in two different ways. They are either given informative beta priors (${Beta(a,b)}$) that are set according to the risk of bias. The default beta priors are as follows: high bias RCT pi.high.rct='dbeta(10,1)', low bias RCT pi.low.rct='dbeta(1,10)', high bias NRS pi.high.nrs='dbeta(30,1)' and low bias NRS pi.low.nrs='dbeta(1,30)'. The ratio ${a/b}$ controls the skewness of the beta distribution. The closer to 1 the ratio a/b, the more the mean of probability of bias gets closer to 1 and the study acquires 'major' bias adjustment. Alternatively, we can use the study characteristics $z_j$ to predict $\pi_j$ through a logistic transformation (internally coded).
Bias-adjusted model 2
Another way to incorporate the RoB of the study is by replacing $\delta_{jk}$ by a "bias-adjusted" relative treatment effect $\theta_{jk}$. Then $\theta_{jk}$ is modeled with a bimodal normal distribution as described in Section 2.5. For more details see @verde_2020.
model IPD only
\begin{equation}
logit(p_{ijk}) =
\begin{cases}
u_j +\beta_{0,j} x_{ijk} & \text{if $k=b$}\
u_j +\theta_{jk} + (\beta_{0,j}+\beta^w_{1,jk}) x_{ijk}+
(\beta^B_{1,jk}-\beta^w_{1,jk}) \bar{x}_{j} & \text{if $k≠b$}
\end{cases}
\end{equation}
model AD only
\begin{equation}
logit(p_{jk}) =
\begin{cases}
u_j & \text{if $k=b$}\
u_j +\theta_{jk} +\beta^B_{1,jk} \bar{x}_{j} & \text{if $k≠b$}
\end{cases}
\end{equation}
where the bias-adjusted relative treatment effect ($\theta_{jk}$) are modeled as follows $$ \theta_{jk} \sim \pi_j N(d_{Ak}-d_{Ab}, \tau_\delta^2) + (1-\pi_j)N(d_{Ak}-d_{Ab}+\gamma_j, \tau_\delta^2+\tau_\gamma^2) $$
Assumptions about the model parameters
The table below summarizes the different assumptions implemented in the package about combining the parameters in the models described above.
| Parameter | Assumptions| Argument in crosnma.model()|
|:----------------------|:-----------|:------|
|relative treatment effect ($\delta_{jk}$)| Random-effects: $\delta_{jk} \sim N(d_{Ak}-d_{Ab}, \tau_\delta^2)$| trt.effect='random' |
| |Common-effect: $\delta_{jk}=d_{Ak}-d_{Ab}$| trt.effect='common'|
|Covariate effect $\beta_{0, j}$ | Independent effects: $\beta_{0, j} \sim N(0, 10^2)$| reg0.effect='independent' |
| |Random-effects: $\beta_{0,j} \sim N(B_0, \tau_{\beta_0})$| reg0.effect='random'|
|Within-study covariate-treatment interaction ($\beta_{1, jk}^W$)| Random-effects: $\beta_{1, jk}^W \sim N(B_{1, Ak}^W-B_{1, Ab}^W, \tau_{\beta_1^W})$| regw.effect='random' |
| |Common-effect: $\beta_{1,jk}^W = B_{1, Ak}^W-B_{1, Ab}^W$| regw.effect='common'|
|Between-study covariate-treatment interaction ($\beta_{1, jk}^B$)| Random-effects: $\beta_{1,jk}^B \sim N(B_{1, Ak}^B-B_{1, Ab}^B, \tau_{\beta_1^B})$| regb.effect='random' |
| |Common-effect: $\beta_{1,jk}^B = B_{1, Ak}^B-B_{1, Ab}^B$| regb.effect='common'|
|Bias effect ($\gamma_{{1,2},j}$)| Random-effects: $\gamma_{{1,2},j} \sim N(\Gamma_{{1,2}}, \tau_{\gamma})$| bias.effect='random' |
| |Common-effect: $\gamma_{{1,2},j}=\Gamma_{{1,2}}$| bias.effect='common'|
|Bias probability ($\pi_j$)| $\pi_j \sim Beta(a,b)$| |
|| $\pi_j = a+bz_j$| |
Synthesis of studies comparing drugs for relapsing-remitting multiple sclerosis
Description of the data
The data we use are fictitious but have been developed to resample to real RCTs with IPD and aggregate data included in @Tramacere15. The studies provide either aggregate data std.data (2 RCTs) or as individual participant data prt.data (3 RCTs and 1 cohort study). Both datasets compare in total four drugs which are anonymized.
The prt.data contains 2950 rows, each row refers to a participant in the study. We display the first few rows of the data set:
head(prt.data)
For each participant, we have information for the outcome relapse (0=no, 1=yes), the treatment label trt, the age (in years) and sex (0 = Female, 1 = Male) of the participant. The following columns are set on study-level (it is repeated for each participant in each study): the study id, the design of the study (needs to be either rct or nrs), the risk of bias on each study (can be set as low, high or unclear) and the year of publication .
The aggregate data has the standard format for meta-analysis
head(std.data)
Analysis
The network should be checked for its connectivity before running the analysis. This is a vital step as the model will run even if the network is not connected.
netplot(prt.data,std.data)
In the following table, we summarize the number of studies from each design and each data format:
knitr::kable(ns.tab(prt.data,std.data))
There are two steps to run the NMA/NMR model. The first step is to create a JAGS model using crosnma.model() which produces the JAGS code and the data. In the second step, the output of that function will be used in crosnma.run() to run the analysis through JAGS.
Naïve network meta-analysis
We start by indicating the names of the datasets on participant-level (prt.data) and aggregate data (std.data). Then, the name of the variables on each dataset needs to be given respectively in prt.data and std.data. Next, the reference treatment needs to be assigned (we set it to drug A). By choosing trt.effect=random, we are assigning a normal distribution to each relative treatment effect to allow the synthesis across studies, see the table in Section 2.1. Finally, the different designs; RCT and NRS are combined with the information taken at face-value; method.bias = 'naive'.
Optionally, we can specify a prior to the common heterogeneity of the treatment effect across studies. We indicate that distribution in the argument prior as tau.trt='dunif(0,3)', see below.
# jags model: code+data mod1 <- crosnma.model(prt.data=prt.data, std.data=std.data, trt=c('trt','trt'), study=c('study','study'), outcome=c('outcome','outcome'), n='n', design=c('design','design'), trt.effect='random', reference='A', method.bias = 'naive', #---------- assign a prior ---------- prior=list(tau.trt='dunif(0,3)') )
Next, we fit the NMA model using crosnma.run()which requires us to set the number of adaptations, iterations, thinning and chains.
# run jags jagsfit1 <- crosnma.run(model=mod1, n.adapt = 50, n.iter=500, n.burnin = 200, thin=1, n.chains=2)
We summarize the estimated parameters in the following table.
knitr::kable(summary(jagsfit1,expo=F))
The estimated OR of B vs A can be obtained as exp(d.B) and similarly for exp(d.C) and exp(d.D) are the ORs of C and D relative to A, respectively. The value of tau refers to the estimates of the heterogeneity standard deviation in the relative treatment effects across studies.
We need also to assess the convergence of the MCMC chains either by visually inspect the trace plot or checking the Gelman and Rubin statistic, Rhat (it should be approximately 1) in the table above.
coda::traceplot(jagsfit1$samples)
Naïve network meta-regression
In this part, we run NMR model by adding age as a covariate from both datasets. We set a list of elements covariate=list(c('age'),c('age')) representing the names of the covariate in prt.data and std.data, respectively.
# jags model: code+data mod2 <- crosnma.model(prt.data=prt.data, std.data=std.data, trt=c('trt','trt'), study=c('study','study'), outcome=c('outcome','outcome'), n='n', design=c('design','design'), reference='A', trt.effect='random', #---------- meta-regression ---------- covariate = list(c('age'),c('age')), split.regcoef = F, #---------- bias adjustment ---------- method.bias='naive' )
The MCMC is run under the same set up as in the network meta-analysis.
# run jags jagsfit2 <- crosnma.run(model=mod2, n.adapt = 50, n.iter=500, n.burnin = 200, thin=1, n.chains=2)
and the output table is presented below
knitr::kable(summary(jagsfit2,expo=FALSE))
Now, we additionally estimate b_1 which indicates the mean effect of age and tau.b_1 which refers to the heterogeneity standard deviation in the effect of age across studies. Here, we obtain a single estimate because we choose to not split the within- and between-study age coefficients $(\beta^W_{1,jk} = \beta^B_{1,jk}=\beta_{1,jk})$ to improve the convergence of MCMC.
Again, we check convergence with trace plots
coda::traceplot(jagsfit2$samples)
Using non-randomized studies as a prior in network meta-regression
To run NMA with a prior from NRS, two additional arguments are needed: we indicate using NRS as a prior by setting method.bias='prior'. That means that the model runs internally NMA with only NRS data which are then used to construct informative priors. This requires defining MCMC settings (the number of adaptations, iterations, burn-ins, thinning and chains) in the argument run.nrs.
In this method, the prior for the basic parameters is set to a normal distribution. For basic parameters not examined in the NRS, the code sets a minimally informative prior d~dnorm(0, 1e-4). To account for possible bias, the means of the distribution can be shifted by mean.shift to reflect the potential bias in NRS and/or the variance can be inflated by var.infl to control the influence of NRS on the final estimation. Both should be provided in run.nrs.
# jags model: code+data mod3 <- crosnma.model(prt.data=prt.data, std.data=std.data, trt=c('trt','trt'), study=c('study','study'), outcome=c('outcome','outcome'), n='n', design=c('design','design'), reference='A', trt.effect='random', #---------- meta-regression ---------- covariate = list(c('age'),c('age')), split.regcoef = F, #---------- bias adjustment ---------- method.bias='prior', run.nrs=list(var.infl=0.6, n.adapt = 500, n.iter=10000, n.burnin = 4000, thin=1, n.chains=2))
# run jags jagsfit3 <- crosnma.run(model=mod3, n.adapt = 50, n.iter=500, n.burnin = 200, thin=1, n.chains=2)
knitr::kable(summary(jagsfit3,expo=F))
coda::traceplot(jagsfit3$samples)
Bias-adjusted model 1
In this part, the overall relative treatment effects are estimated from both NRS and RCT with adjustment to study-specific bias.
To fit the model, we set method.bias='adjust1' and we need to provide the name of the bias variable bias=c('bias','bias') in prt.data and std.data, respectively. By default, the effect of bias is assumed to be additive bias.type='add' and equal across studies bias.effect='common'. We also use the year of study publication to estimate the study-probability of bias, bias.covariate = c('year','year').
# jags model: code+data mod4 <- crosnma.model(prt.data=prt.data, std.data=std.data, trt=c('trt','trt'), study=c('study','study'), outcome=c('outcome','outcome'), n='n', design=c('design','design'), reference='A', trt.effect='random', #---------- bias adjustment ---------- method.bias='adjust1', bias=c('bias','bias'), bias.type='add', bias.effect='common', bias.covariate = c('year','year') )
# run jags jagsfit4 <- crosnma.run(model=mod4, n.adapt = 50, n.iter=500, n.burnin = 200, thin=1, n.chains=2)
The results are presented below
knitr::kable(summary(jagsfit4,expo=F))
The parameter g refers to the mean bias effect, common for all studies.
The trace plots are shown below
coda::traceplot(jagsfit4$samples)
Bias-adjusted model 2
The arguments for method.bias='adjust2' are similar to the ones used before in method.bias='adjust1'.
# jags model: code+data mod5 <- crosnma.model(prt.data=prt.data, std.data=std.data, trt=c('trt','trt'), study=c('study','study'), outcome=c('outcome','outcome'), n='n', design=c('design','design'), reference='A', trt.effect='random', #---------- bias adjustment ---------- method.bias='adjust2', bias=c('bias','bias'), bias.type='add', bias.effect='common', bias.covariate = c('year','year') )
# run jags jagsfit5 <- crosnma.run(model=mod5, n.adapt = 50, n.iter=500, n.burnin = 200, thin=1, n.chains=2)
knitr::kable(summary(jagsfit5,expo=F))
coda::traceplot(jagsfit5$samples)
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.