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 ∼ Bernoulli(p_ijk)

NMR model for AD studies

r_jk ∼ Binomial(p_.jk, n_jk)

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 d̃_Ak^NRS and variance V_Ak^NRS (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 ∼ 𝒩(d̃_Ak^NRS, V_Ak^NRS). We can adjust for potential biases associated with NRS by either shifting the mean of the prior distribution with a bias term ζ or by dividing the prior variance with a common inflation factor w, 0 < w < 1 controls NRS contribution. The assigned priors become d_Ak ∼ 𝒩(d̃_Ak^NRS + ζ, V_Ak^NRS/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 γ_2, jbkR_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 γ_1, jbk^R_j. We can add either multiplicative bias effects, additive bias effects, or both (in this case, δ_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

NMR model for AD studies

The bias indicator R_j follows the following distribution

R_j ∼ Bernoulli(π_j) The bias probabilities π_j are study-specific and can be estimated in two different ways. They are either given informative beta priors (Beta(a₁, a₂)) that are set according to the risk of bias for each study. π_j ∼ Beta(a₁, a₂)

The hyperparameters a₁ and a₂ 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₁/a₂ . When a₁/a₂ equals 1 (or a₁ = a₂), 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₁/a₂ is closer to 1, the more the mean of probability of bias (expected value of π_j = a₁/(a₁ + a₂)) 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 π_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 γ_1, jbk ∼ 𝒩(g_1, bk, τ_1, γ²),γ_2, jbk ∼ 𝒩(g_2, bk, τ_2, γ²)) or common γ_1, jbk = g_1, bk and γ_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.

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

NMR model for AD studies

Then the bias-adjusted relative treatment effect (δ_jbk^bias = δ_jbk + γ_jbk) can be assumed exchangeable across studies δ_jbk^bias ∼ 𝒩(g_bk + d_Ak − d_Ab, τ²/q_j) or fixed as δ_jbk^bias = g_bk + d_Ak − d_Ab. In this parameterisation, instead of assigning prior to the between-study heterogeneity in bias effect τ_γ, we model the RoB weight q_j = τ²/(τ² + τ_γ²) 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 δ_jbk by a “bias-adjusted” relative treatment effect θ_jbk. Then θ_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

NMR model for AD studies

where the bias-adjusted relative treatment effect (θ_jk) are modeled via random-effects model with a mixture of two normal distributions.

θ_jbk ∼ (1 − π_j)𝒩(d_Ak − d_Ab, τ²) + π_j𝒩(d_Ak − d_Ab + γ_jbk, τ² + τ_γ²)

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

θ_jbk = d_Ak − d_Ab + π_jγ_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)