Unadjusted network meta-regression (NMR)

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.

NMR model for IPD studies

\[ y_{ijk} \sim Bernoulli(p_{ijk}) \] \[\begin{equation} logit(p_{ijk}) = \begin{cases} u_{jb} +\beta_{0j} x_{ijk} & \text{if $k=b$}\\ u_{jb} +\delta_{jbk} + \beta_{0j}x_{ijk}+\beta^w_{1,jbk}x_{ijk} + (\beta^B_{1,jbk}-\beta^w_{1,jbk}) \bar{x}_{j} & \text{if $k\ne b$} \end{cases} \end{equation}\]

NMR model for AD studies

\[ r_{jk} \sim Binomial(p_{.jk},n_{jk}) \] \[\begin{equation} logit(p_{.jk}) = \begin{cases} u_{jb} & \text{if $k=b$}\\ u_{jb} +\delta_{jbk} +\beta^B_{1,jbk} \bar{x}_{j} & \text{if $k\ne b$} \end{cases} \end{equation}\]

Using non-randomized studies (NRS) to construct priors for the treatment effects

First, the (network) meta-regression with only NRS data estimates the relative treatment effects with posterior distribution of mean \(\tilde{d}^{NRS}_{Ak}\) and variance \(V^{NRS}_{Ak}\) (use run.nrs in crossnma.model() to control this process). The posteriors of NRS results are then used as priors for the corresponding basic parameters in the RCT model, \(d_{Ak} \sim \mathcal{N}(\tilde{d}^{NRS}_{Ak},V^{NRS}_{Ak})\). We can adjust for potential biases associated with NRS by either shifting the mean of the prior distribution with a bias term \(\zeta\) or by dividing the prior variance with a common inflation factor \(w, 0<w<1\) controls NRS contribution. The assigned priors become \(d_{Ak} \sim \mathcal{N}(\tilde{d}^{NRS}_{Ak}+\zeta,V^{NRS}_{Ak}/w)\).

Bias-adjusted model 1

We incorporate judgments about study risk of bias (RoB) in bias-adjusted model 1 and model 2. Each judgment about the risk of bias in a study is summarized by the index \(R_j\) which takes binary values 0 (no bias) or 1 (bias). In bias-adjusted model 1, we extend the method introduced by Dias et al. (2010) by adding a treatment-specific bias term \(\gamma_{2,jbk} R_j\) to the relative treatment effect on both the AD and IPD parts of the model. A multiplicative model can also be employed, where treatment effects are multiplied by \(\gamma_{1,jbk}^{R_j}\). We can add either multiplicative bias effects, additive bias effects, or both (in this case, \(\delta_{jbk}\) should be dropped from the additive part). The models in previous section are extended to adjust for bias as follows.

NMR model for IPD studies

\[\begin{equation} logit(p_{ijk}) = \begin{cases} u_{jb} +\beta_{0j} x_{ijk} & \text{if $k=b$}\\ u_{jb} +\overbrace{\delta_{jbk} \gamma_{1,jbk}^{R_j}}^{\text{multiplicative}}+\overbrace{\delta_{jbk}+\gamma_{2,jbk} R_j}^{\text{additive}}+ \beta_{0j}x_{ijk}+\beta^w_{1,jbk} x_{ijk}+ (\beta^B_{1,jbk}-\beta^w_{1,jbk}) \bar{x}_{j} & \text{if $k\ne b$} \end{cases} \end{equation}\]

NMR model for AD studies

\[\begin{equation} logit(p_{.jk}) = \begin{cases} u_{jb} & \text{if $k=b$}\\ u_{jb} +\overbrace{\delta_{jbk} \gamma_{1,jbk}^{R_j}}^{\text{multiplicative}}+\overbrace{\delta_{jbk}+\gamma_{2,jbk} R_j}^{\text{additive}}+ \beta^B_{1,jbk} \bar{x}_{j} & \text{if $k\ne b$} \end{cases} \end{equation}\]

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_1,a_2)}\)) that are set according to the risk of bias for each study. \[ \pi_j \sim Beta(a_1, a_2) \]

The hyperparameters \(a_1\) and \(a_2\) should be chosen in a way that reflects the risk of bias for each study. The degree of skewness in beta distribution can be controlled by the ratio \(a_1/a_2\) . When \(a_1/a_2\) equals 1 (or \(a_1=a_2\)), there is no skewness in the beta distribution (the distribution is reduced to a uniform distribution), which is appropriate for studies with unclear risk of bias. When the ratio \(a_1/a_2\) is closer to 1, the more the mean of probability of bias (expected value of \(\pi_j=a_1/(a_1+a_2))\) gets closer to 1 and the study acquires ‘major’ bias adjustment. The default beta priors are as follows: high bias RCT prior.pi.high.rct='dbeta(10, 1)', low bias RCT prior.pi.low.rct = 'dbeta(1, 10)', high bias NRS prior.pi.high.nrs = 'dbeta(30, 1)' and low bias NRS prior.pi.low.nrs = 'dbeta(1, 30)'. Alternatively, we can use the study characteristics \(z_j\) to estimate \(\pi_j\) through a logistic transformation (internally coded).

We combine the multiplicative and the additive treatment-specific bias effects across studies by assuming they are exchangeable \(\gamma_{1,jbk}\sim \mathcal{N}(g_{1,bk},\tau_{1,\gamma}^2 )\),\(\gamma_{2,jbk}\sim \mathcal{N}(g_{2,bk},\tau_{2,\gamma}^2 )\)) or common \(\gamma_{1,jbk}=g_{1,bk}\) and \(\gamma_{2,jbk}=g_{2,bk}\). Dias et al. (2010) proposed to model the mean bias effect \((g_{1,bk}, g_{2,bk})\) based on the treatments being compared.

\[\begin{equation} g_{m,bk} = \begin{cases} g_m & \text{if $b$ is inactive treatment}\\ 0 \text{ or } (-1)^{dir_{bk}} g_m^{act} & \text{if $b$ and $k$ are active treatments} \end{cases} \end{equation}\] where \(m={1,2}\). This approach assumes a common mean bias for studies that compare active treatments with an inactive treatment (placebo, standard or no treatment). For active vs active comparisons, we could assume either a zero mean bias effect or a common bias effect \(g_m^{act}\). The direction of bias \(dir_{bk}\) in studies that compare active treatments with each other should be defined in the data. That is set to be either 0, meaning that bias favors \(b\) over \(k\), or 1 , meaning that \(k\) is favored to \(b\). In crossnma.model(), the bias direction is specified by providing the unfavoured treatment for each study, unfav. To select which mean bias effect should be applied, the user can provide the bias.group column as data. Its values can be 0 (no bias adjustment), 1 (to assign for the comparison mean bias effect \(g_m\)) or 2 (to set bias \(g_m^{act}\)).

Another parameterisation of the logistic model with additive bias effect is

NMR model for IPD studies

\[\begin{equation} logit(p_{ijk}) = \begin{cases} u_{jb} +\beta_{0j} x_{ijk} & \text{if $k=b$}\\ u_{jb} +(1-R_j)\delta_{jbk}+\delta_{jbk}^{bias}R_j+ \beta_{0j}x_{ijk}+\beta^w_{1,jbk} x_{ijk}+ (\beta^B_{1,jbk}-\beta^w_{1,jbk}) \bar{x}_{j} & \text{if $k\ne b$} \end{cases} \end{equation}\]

NMR model for AD studies

\[\begin{equation} logit(p_{.jk}) = \begin{cases} u_{jb} & \text{if $k=b$}\\ u_{jb} +(1-R_j)\delta_{jbk}+\delta_{jbk}^{bias}R_j+ \beta^B_{1,jbk} \bar{x}_{j} & \text{if $k\ne b$} \end{cases} \end{equation}\]

Then the bias-adjusted relative treatment effect (\(\delta_{jbk}^{bias}=\delta_{jbk}+\gamma_{jbk}\)) can be assumed exchangeable across studies \(\delta_{jbk}^{bias} \sim \mathcal{N}(g_{bk}+d_{Ak}-d_{Ab}, \tau^2 /q_j)\) or fixed as \(\delta_{jbk}^{bias}=g_{bk} + d_{Ak}-d_{Ab}\). In this parameterisation, instead of assigning prior to the between-study heterogeneity in bias effect \(\tau_{\gamma}\), we model the RoB weight \(q_j=\tau^2/(\tau^2+\tau_{\gamma}^2\)) for each study. This quantity \(0<q_j<1\) quantifies the proportion of the between-study heterogeneity that is not explained by accounting for risk of bias. The values of \(v\) determine the extent studies at high risk of bias will be down-weighted on average. Setting \(v=1\) gives \(E(q_j )=v/(v+1)=0.5\), which means that high risk of bias studies will be penalized by 50% on average. In crossnma.model(), the user can assign the average down-weight \(E(q_j )\) to the argument down.wgt.

Bias-adjusted model 2

Another way to incorporate the RoB of the study is by replacing \(\delta_{jbk}\) by a “bias-adjusted” relative treatment effect \(\theta_{jbk}\). Then \(\theta_{jbk}\) is modeled with a bimodal normal distribution as described in Section 2.5. For more details see Verde (2020).

NMR model for IPD studies

\[\begin{equation} logit(p_{ijk}) = \begin{cases} u_{jb} +\beta_{0j} x_{ijk} & \text{if $k=b$}\\ u_{jb} +\theta_{jbk} + \beta_{0j} x_{ijk}+\beta^w_{1,jbk} x_{ijk}+ (\beta^B_{1,jbk}-\beta^w_{1,jbk}) \bar{x}_{j} & \text{if $k\ne b$} \end{cases} \end{equation}\]

NMR model for AD studies

\[\begin{equation} logit(p_{jk}) = \begin{cases} u_{jb} & \text{if $k=b$}\\ u_{jb} +\theta_{jbk} +\beta^B_{1,jbk} \bar{x}_{j} & \text{if $k\ne b$} \end{cases} \end{equation}\]

where the bias-adjusted relative treatment effect (\(\theta_{jk}\)) are modeled via random-effects model with a mixture of two normal distributions.

\[ \theta_{jbk} \sim (1-\pi_j) \mathcal{N}(d_{Ak}-d_{Ab}, \tau^2) + \pi_j \mathcal{N}(d_{Ak}-d_{Ab}+\gamma_{jbk}, \tau^2+\tau_\gamma^2) \]

Alternatively, we can summarize these relative effects assuming a common-effect model

\[ \theta_{jbk}= d_{Ak}-d_{Ab}+\pi_j \gamma_{jbk} \]

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.

Description of the data

The data we use are fictitious but resemble real RCTs with IPD and aggregate data included in Tramacere and Filippini (2015). The studies provide either individual participant data ipddata (3 RCTs and 1 cohort study) or aggregate data stddata (2 RCTs). In total, four drugs are compared which are anonymized.

The ipddata contains 1944 participants / rows. We display the first few rows of the data set:

dim(ipddata)
#> [1] 1944   10
head(ipddata)
#>   id relapse treat design age sex rob unfavored bias.group year
#> 1  1       0     D    rct  22   1 low         1          1 2002
#> 2  1       0     D    rct  31   1 low         1          1 2002
#> 3  1       0     D    rct  34   1 low         1          1 2002
#> 4  1       0     D    rct  38   0 low         1          1 2002
#> 5  1       0     D    rct  46   0 low         1          1 2002
#> 6  1       0     D    rct  45   0 low         1          1 2002

For each participant, we have information for the outcome relapse (0 = no, 1 = yes), the treatment label treat, 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 id, the design of the study (needs to be either "rct" or "nrs"), the risk of bias rob on each study (can be set as low, high or unclear), the year of publication, the bias.group for the study comparison and the study unfavoured treatment unfavored.

The aggregate meta-analysis data must be in long arm-based format with the exact same variable names and an additional variable with the sample sizes:

stddata
#>   id   n relapse treat design  age sex  rob unfavored bias.group year
#> 1  1  25      19     A    rct 34.3 0.2 high         0          1 2010
#> 2  1  25      11     C    rct 34.3 0.3 high         1          1 2010
#> 3  2 126      97     A    rct 30.0 0.4 high         0          1 2015
#> 4  2 125      89     C    rct 30.0 0.5 high         1          1 2015

Analysis

There are two steps to run the NMA/NMR model. The first step is to create a JAGS model using crossnma.model() which produces the JAGS code and the data. In the second step, the output of that function will be used in crossnma() to run the analysis through JAGS.

# JAGS model: code + data
mod1 <- crossnma.model(treat, id, relapse, n, design,
  prt.data = ipddata, std.data = stddata,
  #---------- bias adjustment ----------
  method.bias = "naive",
  #---------- assign a prior ----------
  prior.tau.trt = "dunif(0, 3)",
  #---------- SUCRA ----------
  sucra = TRUE, small.values = "desirable"
  )
#> Both designs are combined naively without acknowledging design differences.

netgraph(mod1, cex.points = n.trts, adj = 0.5, plastic = FALSE,
 number = TRUE, pos.number.of.studies = c(0.5, 0.4, 0.5, 0.5, 0.6, 0.5))

# Run JAGS
jagsfit1 <- crossnma(mod1, n.iter = 5000, n.burnin = 2000, thin = 1)
jagsfit1
#>           Mean    SD  2.5%   50% 97.5%  Rhat n.eff
#> exp(d.A)     .     .     .     .     .     .     .
#> exp(d.B) 0.470 0.194 0.322 0.468 0.692 1.009   865
#> exp(d.C) 0.629 0.211 0.416 0.629 0.944 1.015   582
#> exp(d.D) 0.341 0.265 0.202 0.341 0.571 1.004   671
#> tau      0.188 0.182 0.002 0.140 0.673 1.044   241
#> SUCRA.A  0.006 0.047 0.000 0.000 0.000 1.010  2547
#> SUCRA.B  0.680 0.153 0.333 0.667 1.000 1.000  1349
#> SUCRA.C  0.368 0.128 0.333 0.333 0.667 1.004  1382
#> SUCRA.D  0.947 0.142 0.667 1.000 1.000 1.000  1013
#> Mean and quantiles are exponentiated

par(mar = rep(2, 4), mfrow = c(2, 3))
plot(jagsfit1)

# JAGS model: code + data
mod2 <- crossnma.model(treat, id, relapse, n, design,
  prt.data = ipddata, std.data = stddata,
  #---------- bias adjustment ----------
  method.bias = "naive",
  #----------  meta-regression ----------
  cov1 = age,
  split.regcoef = FALSE
  )
#> Both designs are combined naively without acknowledging design differences.

# Run JAGS
jagsfit2 <- crossnma(mod2, n.iter = 5000, n.burnin = 2000, thin = 1)

league(jagsfit2, cov1.value = 38, digits = 2)
#>                                                             
#>                    A 0.41 (0.21 to 0.79) 0.63 (0.30 to 1.33)
#>  2.44 (1.27 to 4.87)                   B 1.51 (0.75 to 3.38)
#>  1.59 (0.75 to 3.32) 0.66 (0.30 to 1.33)                   C
#>  2.68 (0.87 to 6.66) 1.11 (0.34 to 2.79) 1.67 (0.56 to 4.57)
#>                     
#>  0.37 (0.15 to 1.14)
#>  0.90 (0.36 to 2.95)
#>  0.60 (0.22 to 1.79)
#>                    D

league(jagsfit2, cov1.value = 38, digits = 2, direction = "long")
#>  Treatment Comparator median 2.5% 97.5%
#>          B          A   0.41 0.21  0.79
#>          C          A   0.63 0.30  1.33
#>          D          A   0.37 0.15  1.14
#>          A          B   2.44 1.27  4.87
#>          C          B   1.51 0.75  3.38
#>          D          B   0.90 0.36  2.95
#>          A          C   1.59 0.75  3.32
#>          B          C   0.66 0.30  1.33
#>          D          C   0.60 0.22  1.79
#>          A          D   2.68 0.87  6.66
#>          B          D   1.11 0.34  2.79
#>          C          D   1.67 0.56  4.57

# JAGS model: code + data
mod3 <- crossnma.model(treat, id, relapse, n, design,
  prt.data = ipddata, std.data = stddata,
  reference = "D",
  #----------  meta-regression ----------
  cov1 = age,
  split.regcoef = FALSE,
  #---------- bias adjustment ----------
  method.bias = "prior",
  run.nrs.trt.effect= "common",
  run.nrs.var.infl = 0.6, run.nrs.mean.shift = 0,
  run.nrs.n.iter = 10000, run.nrs.n.burnin = 4000,
  run.nrs.thin = 1, run.nrs.n.chains = 2
  )

# Run JAGS
jagsfit3 <- crossnma(mod3, n.iter = 5000, n.burnin = 2000, thin = 1)

heatplot(jagsfit3, cov1.value = 38,
  size = 6, size.trt = 20, size.axis = 12)

# JAGS model: code + data
mod4 <- crossnma.model(treat, id, relapse, n, design,
  prt.data = ipddata, std.data = stddata,
  #---------- bias adjustment ----------
  method.bias = "adjust1",
  bias.type = "add",
  bias.effect = "common",
  bias = rob,
  unfav = unfavored,
  bias.group = bias.group,
  bias.covariate = year
)
#> Bias effect is assumed common across studies

# Run JAGS
jagsfit4 <- crossnma(mod4, n.iter = 5000, n.burnin = 2000, thin = 1)

# JAGS model: code + data
mod5 <- crossnma.model(treat, id, relapse, n, design,
  prt.data = ipddata, std.data = stddata,
  #---------- bias adjustment ----------
  method.bias = "adjust2",
  bias.type = "add",
  bias = rob,
  unfav = unfavored,
  bias.group = bias.group
)

# Run JAGS
jagsfit5 <- crossnma(mod5, n.iter = 5000, n.burnin = 2000, thin = 1)