Alternative Estimation for Network Meta-analysis

This vignette describes the bootstrap, jackknife, and permutation structure for network_meta(). The examples are intentionally small and use model_type = "wls" so the code path is easy to inspect. The same wrapper structure calls network_meta() for multilevel and multivariate network models; only the fitting arguments change.

library(mars)

Example Data

The resampling functions expect the same contrast-level data layout used by network_meta(): one row per observed treatment contrast, with an effect defined as treatment_2 - treatment_1.

nma_dat <- data.frame(
  study = c("s1", "s1", "s2", "s3", "s4", "s4", "s5", "s6"),
  trt1 = c("A", "A", "A", "A", "B", "A", "B", "A"),
  trt2 = c("B", "C", "B", "C", "C", "B", "C", "C"),
  yi = c(0.20, 0.42, 0.12, 0.48, 0.26, 0.15, 0.31, 0.44),
  vi = c(0.04, 0.05, 0.03, 0.06, 0.05, 0.04, 0.05, 0.04),
  stringsAsFactors = FALSE
)

The baseline model is fit with network_meta():

fit <- network_meta(
  data = nma_dat,
  study_id = "study",
  treatment_1 = "trt1",
  treatment_2 = "trt2",
  effect = "yi",
  variance = "vi",
  reference = "A",
  model_type = "wls"
)

fit$effects[, c("treatment_1", "treatment_2", "direct", "indirect", "total")]
#>   treatment_1 treatment_2    direct     indirect     total
#> 1           A           B 0.1530000  0.001426105 0.1544261
#> 2           A           C 0.4443243 -0.001927168 0.4423972
#> 3           B           C 0.2850000  0.002971051 0.2879711

Model Estimand

For a contrast row $i$ comparing treatment $a_i$ with treatment $b_i$ , network_meta() interprets the observed effect as treatment $b_i$ minus treatment $a_i$ :

$y_i = d_{b_i} - d_{a_i} + \varepsilon_i, \qquad \varepsilon_i \sim N(0, v_i).$

With reference treatment $r$ , $d_r = 0$ . The fixed-effect design row $x_i$ therefore has a $+1$ for $b_i$ , a $-1$ for $a_i$ , and zeros elsewhere, excluding the reference column. In matrix form, the common-effect weighted least-squares estimator used in this vignette is:

$\hat d = (X^\top W X)^{-1} X^\top W y, \qquad W = V^{-1}.$

For independent contrast rows, $V = \operatorname{diag}(v_i)$ . When multi-arm covariance is requested, network_meta() builds block covariance matrices by study and uses the same estimator with the block-structured $V$ .

For multilevel network meta-analysis, the same contrast equation is augmented with random components:

$y_i = x_i^\top d + z_i^\top u + \varepsilon_i, \qquad u \sim N(0, \Sigma_u), \quad \varepsilon_i \sim N(0, V_i).$

For multivariate network meta-analysis, comparisons are treated as multivariate outcomes within study:

$y_s \sim N(X_s d, S_s + \Psi),$

where $S_s$ is the within-study sampling covariance and $\Psi$ is the between-study heterogeneity covariance. The resampling functions below do not refit these equations themselves; they construct resampled data and then call network_meta() so the same model-building path is reused.

Bootstrap Structure

bootstrap_network_meta() resamples studies, not individual rows. This matters for network data because a multi-arm study can contribute more than one contrast. Every sampled study contributes all of its contrast rows. If a study is sampled more than once, the copied study blocks receive new study IDs so they are treated as independent bootstrap blocks inside network_meta().

If $S = \{1,\ldots,m\}$ indexes studies, bootstrap replicate $b$ samples study IDs $S_b^\ast = (s_1^\ast,\ldots,s_m^\ast)$ with replacement. The model is refit on all contrast rows belonging to the sampled studies:

$\hat\theta_b^\ast = T(D_b^\ast),$

where $T(\cdot)$ is the complete network_meta() estimation pipeline and $\theta$ can be a basic effect, total contrast, direct estimate, indirect estimate, ranking quantity, fit statistic, or heterogeneity summary.

set.seed(1)
boot_one <- bootstrap_network_meta(
  data = nma_dat,
  study_id = "study",
  treatment_1 = "trt1",
  treatment_2 = "trt2",
  effect = "yi",
  variance = "vi",
  reference = "A",
  model_type = "wls"
)
boot_one$effects[, c("treatment_1", "treatment_2", "total")]
#>   treatment_1 treatment_2     total
#> 1           A           B 0.1618535
#> 2           A           C 0.4414566
#> 3           B           C 0.2796031

Replicated bootstraps use the same future.apply helper as the existing MARS bootstrap framework.

set.seed(2)
boot_reps <- replicate_network_bootstrap(
  number_bootstraps = 5,
  data = nma_dat,
  study_id = "study",
  treatment_1 = "trt1",
  treatment_2 = "trt2",
  effect = "yi",
  variance = "vi",
  reference = "A",
  model_type = "wls"
)

boot_stats <- summary(boot_reps)
boot_stats$effects_summary[boot_stats$effects_summary$effect_type == "total", ]
#>   effect_type treatment_1 treatment_2 point_estimate  ci_lower  ci_upper
#> 7       total           A           B      0.1473134 0.1325665 0.1792851
#> 8       total           A           C      0.4423972 0.4303405 0.4539479
#> 9       total           B           C      0.2967089 0.2674034 0.3064385

The full summary table retains named scalar estimands for programmatic use, while effects_summary reshapes the direct, indirect, and total rows into the same treatment-pair structure used by network_meta() output:

basic_effects.* are treatment effects against the network reference.
direct.*, indirect.*, and total.* are comparison-level estimates.
ranking.* terms summarize treatment ranking outputs when available.
fit_stats.* terms summarize model-fit quantities.

Jackknife Structure

jackknife_network_meta() removes one study at a time and then calls network_meta() on the remaining network. Some omissions can disconnect a network or make a model rank deficient. The on_error argument controls whether those failed fits stop immediately, warn and continue, or are kept silently as failure objects.

For omitted study $s$ , the jackknife replicate is:

$\hat\theta_{(-s)} = T(D \setminus D_s).$

The jackknife standard error reported by summary() is:

$\widehat{SE}_{jack}(\theta) = \sqrt{\frac{m - 1}{m} \sum_{s = 1}^m \left(\hat\theta_{(-s)} - \bar\theta_{(\cdot)}\right)^2},$

where $m$ is the number of leave-one-study-out replicates and $\bar\theta_{(\cdot)}$ is their mean.

jack_reps <- replicate_network_jackknife(
  data = nma_dat,
  study_id = "study",
  treatment_1 = "trt1",
  treatment_2 = "trt2",
  effect = "yi",
  variance = "vi",
  reference = "A",
  model_type = "wls"
)
names(jack_reps)
#> [1] "s1" "s2" "s3" "s4" "s5" "s6"

jack_stats <- summary(jack_reps)
jack_stats$effects_summary[jack_stats$effects_summary$effect_type == "total", ]
#>   effect_type treatment_1 treatment_2 point_estimate  ci_lower  ci_upper
#> 7       total           A           B      0.1531677 0.1436727 0.1683695
#> 8       total           A           C      0.4427489 0.4342128 0.4484432
#> 9       total           B           C      0.2853494 0.2787765 0.2992431
#>   jackknife_se
#> 7   0.01854617
#> 8   0.01183766
#> 9   0.01866164

For jackknife summaries, jackknife_se is computed from the leave-one-study-out replicate distribution. The interval columns are still empirical quantiles of the replicate estimates, matching the package’s existing alternative-estimation summary style.

Permutation Structure

permutation_network_meta() implements an effect-value permutation. The study IDs, treatment labels, network geometry, and sampling variances stay fixed; only the effect-size column is shuffled across contrast rows.

This is intentionally different from a row-order permutation. Reordering rows would usually produce the same network meta-analysis fit, because network_meta() is order-invariant. Shuffling effect values creates a null-like replicate in which the observed network structure is preserved but the association between comparisons and observed effects is broken.

For permutation replicate $p$ , let $\pi_p$ be a random permutation of the row indices. The permuted data replace $y_i$ with $y_{\pi_p(i)}$ , while $(a_i, b_i, v_i, s_i)$ remain fixed:

$D_p^\pi = \{(s_i, a_i, b_i, y_{\pi_p(i)}, v_i): i = 1,\ldots,n\}, \qquad \hat\theta_p^\pi = T(D_p^\pi).$

set.seed(3)
perm_reps <- replicate_network_permutation(
  number_permutations = 5,
  data = nma_dat,
  study_id = "study",
  treatment_1 = "trt1",
  treatment_2 = "trt2",
  effect = "yi",
  variance = "vi",
  reference = "A",
  model_type = "wls"
)

perm_stats <- summary(perm_reps)
perm_stats$effects_summary[perm_stats$effects_summary$effect_type == "total", ]
#>   effect_type treatment_1 treatment_2 point_estimate  ci_lower  ci_upper
#> 7       total           A           B      0.2282021 0.1339259 0.2566145
#> 8       total           A           C      0.3920264 0.3533489 0.4228029
#> 9       total           B           C      0.1763357 0.1224416 0.2854617

The permutation output should be interpreted as a resampling diagnostic for the implemented null structure, not as a full design-consistency test. It answers: what kinds of network estimates appear when the same comparison graph and variance structure are retained, but the observed effects are exchangeable across rows?

Unified Dispatcher

The convenience dispatcher mirrors mars_alt_estimation() and routes to the network-specific replicate functions.

set.seed(4)
alt <- network_alt_estimation(
  type = "bootstrap",
  number_bootstraps = 3,
  data = nma_dat,
  study_id = "study",
  treatment_1 = "trt1",
  treatment_2 = "trt2",
  effect = "yi",
  variance = "vi",
  reference = "A",
  model_type = "wls"
)

length(alt)
#> [1] 3
summary(alt)$number_success
#> [1] 3

The fitted replicate objects are ordinary nma_mars objects. This keeps the implementation simple: resampling changes the data, and network_meta() remains responsible for constructing the design matrix, fitting the model, computing direct/indirect/total effects, rankings, heterogeneity summaries, and diagnostics.