2.1 AB-NMA model

Before describing the tipping point analysis method, we first present the model specification for the AB-NMA. In an NMA reviewing N studies and K treatments, let i denote the study, $i = 1, 2, \dots , N$ ; let k denote the treatment ( $k = 1, 2, \dots , K$ ), and $T_i$ is the subset of treatments investigated in study i. In study i for treatment k ( $k \in T_i$ ), let $y_{ik}$ denote the number of events for a binary outcome, $n_{ik}$ denote the total number of subjects, and $p_{ik}$ denote the underlying absolute risk. The AB-NMA model with a logit link can be specified as

$$ \begin{align*}\begin{aligned} &y_{ik} \sim Bin(n_{ik}, p_{ik}) \text{ for } k \in T_i;\\ &\text{logit}(p_{ik})=\mu_k+v_{ik};\\ &(v_{i1}, v_{i2}, \dots, v_{iK})^T \sim MVN(\mathbf{0}, \boldsymbol{\Sigma}),\\ \end{aligned}\end{align*} $$

where $\mu _k$ denotes the fixed effect of the treatment k, $v_{ik}$ denotes the random effect of the treatment k in study i, and $\boldsymbol {\Sigma }$ denotes the variance–covariance matrix of the vector of random effects $(v_{i1}, v_{i2}, \dots , v_{iK})^T$ . The population-averaged absolute risk of treatment k can be approximated as

$$ \begin{align*} E\left(p_{ik} \mid \mu_k, \sigma_k \right) & \approx \operatorname{expit}\left\{\mu_k / \sqrt{1 + \frac{256}{75 \pi^2} \sigma_{k}^2}\right\}, \end{align*} $$

where $\sigma _{k}$ denotes the standard deviation of the random effect of the treatment k, and the relative risk (RR), risk difference (RD), and odds ratio (OR) can be calculated based on the absolute risk using consistency equations.Reference Wang, Lin, Hodges, MacLehose and Chu ^{21 ,} Reference Zeger, Liang and Albert ³¹

As recommended by Wang et al.,Reference Wang, Lin, Hodges, MacLehose and Chu ²¹ we used the separation strategy and the variance shrinkage method to assign weakly informative priors to $\boldsymbol {\Sigma }$ . It can be decomposed as $\boldsymbol {\Sigma } = \boldsymbol {\Delta } \mathbf {R} \boldsymbol {\Delta }$ , where $\boldsymbol {\Delta } = \operatorname {diag}(\sigma _1, \dots , \sigma _K)$ , and $\mathbf {R}$ is the correlation matrix. By assuming the exchangeable correlation structure, the correlation matrix is a K by K matrix as follows:

$$ \begin{align*}\mathbf{R} = \begin{bmatrix} 1 & \rho & \rho & \cdots & \rho \\ \rho & 1 & \rho & \cdots & \rho \\ \rho & \rho & 1 & \cdots & \rho \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \rho & \rho & \rho & \cdots & 1 \end{bmatrix}, \end{align*} $$

where $\rho $ denotes the correlation.

Then, a uniform prior is used for the standard deviationReference Gelman ³² :

$$ \begin{align*}\sigma_k \sim U(0.0001, 10), \text{ for } k=1, 2, \dots, K.\end{align*} $$

Alternatively, a hierarchical half-Cauchy (HHC) prior for the standard deviation is recommended by Wang et al.Reference Wang, Lin, Hodges, MacLehose and Chu ²¹ :

$$ \begin{align*}\begin{aligned} &\sigma_k \sim HC(a), \text{ for } k=1, 2, \dots, K; \\ &a \sim U(0, 5).\\ \end{aligned}\end{align*} $$

A uniform prior is commonly used for the correlation parameter $\rho $ ,

$$ \begin{align*}\rho \sim U\left(-\frac{1}{K-1},1\right),\end{align*} $$

where the lower bound is chosen to guarantee that the correlation matrix $\mathbf {R}$ is positive definite. For the fixed effect, $\mu _k$ , we assign a non-informative prior:

$$ \begin{align*}\mu_k \sim N(0, 100^2).\end{align*} $$

2.2 Method for tipping point analysis

Our tipping point analysis focuses on the correlation parameter $\rho $ . A tipping point is defined as the value for $\rho $ where some substantive conclusion changes, e.g., the interval conclusion flips or the change in magnitude of the relative treatment effect exceeds some meaningful threshold. The proposed procedure for searching for a tipping point is as follows:

1. 1. Conduct the NMA with the model specified in Section 2.1 to estimate:

   • • Population-averaged point estimate and 95% credible interval of the treatment effect, e.g., RR or RD, for all pairs of treatments;

   • • The posterior median, 95% credible interval, and a series of other posterior percentiles for the correlation parameter $\rho $ .

   The number of percentiles for $\rho $ depends on how granular one wants to be when checking for a tipping point. We recommend using the following 11 percentiles of the posterior samples: 1%, 2.5%, 5%, 10%, 25%, 50%, 75%, 90%, 95%, 97.5%, and 99%. These were chosen to include the important posterior percentiles that are commonly reported in the summary of a Bayesian posterior distribution: 2.5%, 25%, 50%, 75%, and 97.5%,Reference Carlin and Louis ³³ supplemented by additional points to provide more granularity. While a tighter grid of percentiles may provide a more accurate estimate of the tipping point, it is more computationally intensive.

2. 2. Repeatedly conduct the NMA with the model described in Section 2.1, but instead of assigning a prior distribution to $\rho $ , set $\rho $ to each value of the estimated percentiles in Step 1. Thus, 11 NMA models with different assumed values for $\rho $ values are fitted separately, and the mean and 95% credible intervals of the treatment effects for all pairs of treatments are obtained in each model. Draw the interval conclusion of the treatment effect based on whether the 95% credible interval includes the null value, e.g., the null value is 1 for RR and 0 for RD.

3. 3. Search for a tipping point at the 11 data points based on either of the following criteria:

   • • Interval conclusion change: An interval conclusion change occurs when the interval conclusion of a treatment effect comparing a pair of treatments based on the result of one or more NMAs in Step 2 is opposite to the interval conclusion for this pair of treatments based on the result of the NMA in Step 1. For example, an interval conclusion can shift from strong evidence supporting the presence of a treatment effect (i.e., 95% credible interval does not include the null value) to weak evidence (i.e., 95% credible interval includes the null value), or vice versa. The 95% credible interval has been widely used for such conclusion because a posterior probability threshold of 0.975 is commonly accepted in determining the strength of evidence. However, alternative thresholds other than 0.975 can be carefully selected depending on the clinical context and specific study objectives.

   • • Magnitude change: A magnitude change happens when the percent change in the magnitude of the mean treatment effect comparing a pair of treatments in Step 2 relative to the result of the NMA in Step 1 exceeds some meaningful threshold. The choice of thresholds should ideally be informed by clinical expertise and tailored to the specific disease context and primary outcomes. To demonstrate the methodology, in the selected NMA datasets, we considered thresholds of $\pm $ 15% and $\pm $ 30%, corresponding to low and high thresholds.

A schematic summary of the above steps is presented in Figure 1. In Step 2, the models can be executed in parallel for a comprehensive search. Alternatively, to reduce the computational demands, we may also adopt a bisection approach to pinpoint the location of the tipping point by iteratively narrowing a range of search. Starting with a broad range, for example, from 1 to 99 percentiles, we only compare two models fitted with the $\rho $ values at the two endpoints of the range. If we observe either an interval conclusion change or a magnitude change, the search continues and the range is narrowed down in the next iteration; otherwise, the searching process ends, concluding no tipping point.

Figure 1 Steps of searching for tipping points in correlation parameters that alter conclusions about relative treatment effects in an NMA.

All analyses were conducted on JAGS version 4.3.0-gcc7.2.0 through the R package “rjags” in R version 4.2.2. The JAGS code for Bayesian models is provided in Appendix 1 of the Supplementary Material.

3.1 Dataset screening and selection

We extracted 453 NMA datasets published between 1999 and 2015Reference Petropoulou, Nikolakopoulou and Veroniki ³⁴ from the R package “nmadb,”Reference Papakonstantinou ³⁵ which provides the Application Programming Interface (API) for an NMA database. With computational considerations and with the purpose to include various types of networks with a varied sparsity to demonstrate our method on networks that reflect more typical real-world scenarios, we evaluated NMA datasets with the following inclusion criteria: (1) a binary outcome, (2) $\le $ 5 treatments, and (3) $\ge $ 30 studies. The detailed workflow for this screening process is illustrated in Figure 2. Using these criteria, we identified 14 candidate NMA datasets that collectively comprised 790 studies with 355,923 total participants.

Figure 2 A diagram of the dataset screening process.

Note: 453 datasets were extracted from the “nmadb” package in R. All selections were based on the record in the package.

The candidate NMA datasets are named after the first author and the year of publication. Table 1 describes each of these datasets in terms of the number of studies, number of treatments, primary outcome, total number of events, total number of subjects, and the median (minimum, maximum) event rate across the studies. These networks encompass a diverse range of therapeutic areas such as cardiovascular (Ribeiro2011), respiratory (Puhan2009, Mills2010, and Eisenberg2008), gastroenterology (Yang2014b and Tadrous2014), neurology and psychiatry (Furukawa2014), hematology and oncology (Wang2015), and infectious disease (Edwards2009b). The median event rates in these studies spanned from 0.00 to 0.88. This broad range indicated the inclusion of both rare and common primary outcomes within the selected networks.

Table 1 A summary of selected network meta-analyses

Figure 3 shows the network plots for the 14 selected NMA datasets. In these plots, the nodes symbolize the treatments and are labeled with uppercase letters, and the edges between nodes indicate direct treatment comparisons. The weight of each edge is determined by the number of studies with a direct comparison between the connected treatments. Additionally, the size of each node is determined by the number of studies containing the specific treatment, offering an intuitive visual representation of the treatment’s prevalence within the network. The selected networks exhibited a diverse range of network plot configurations and sparsity characteristics. For example, Gafter-Gvili2005 and Edwards2009b had high sparsity demonstrated by star-shaped structures, while Hutton2012 and Puhan2009 had low sparsity with well-connected networks. Besides these two types of extreme networks, others showed moderate sparsity with more complex networks with loops and some groups with sparse connections.

Figure 3 Network plots of selected NMA datasets.

Note: The 4th plot in the first row represents the network selected as the case study. The nodes with uppercase letters indicate the distinct treatments in the network, and the edges indicate the direct comparisons in an RCT. The weight of each edge is determined by the number of studies with a direct comparison between the connected treatments.

Figure 4 displays the forest plot of the estimated correlation between treatments across the 14 selected NMA datasets using AB-NMA models as specified in Section 2.1. These networks exhibited a broad range of posterior median correlations, from 0.266 to 0.996. Additionally, the 95% credible intervals of the correlation differed in length, as evidenced by the varying lengths of the blue lines. Such diversity in the datasets reflected the comprehensive range of selected studies and allowed us to demonstrate the applicability of the tipping point analysis across diverse real-world examples.

Figure 4 Forest plot of the posterior median correlation estimates based on AB-NMA model in selected NMA datasets.

3.2 Dataset of the case study

To illustrate the proposed tipping point analysis in more detail, we selected one network from the selected 14 networks, the Trikalinos2009Reference Trikalinos, Alsheikh-Ali, Tatsioni, Nallamothu and Kent ³⁶ NMA dataset, as a detailed case study. The NMA study by Trikalinos et al. from 2009 in the Lancet evaluated catheter-based treatments, specifically percutaneous coronary intervention (PCI), for non-acute Coronary Artery Disease (CAD). The study focused on medical therapy (including administration of antiplatelet agents, $\beta $ blockers, nitrates, calcium channel blockers, or aggressive lipid-lowering treatment)Reference Bucher, Hengstler, Schindler and Guyatt ³⁷ and 3 main types of PCI interventions: Percutaneous Transluminal Balloon Coronary Angioplasty (PTCA), Bare Metal Stents (BMS), and Drug-Eluting Stents (DES). There were 4 outcomes of interest: all-cause mortality, myocardial infarction (MI), coronary artery bypass grafting (CABG), target lesion or vessel revascularization (TLR/TVR), and any revascularization. In our case study, we focus on the primary outcome, all-cause mortality.

The Trikalinos2009 NMA included 61 RCTs and a total of 25,388 participants. Since two of the trials, which reported results for 2 separately randomized strata, were entered as 4 distinct entries in the meta-analyses, this NMA dataset comprised 63 unique trial entries and 26,748 subjects. Of the 6 possible pairwise treatment comparisons, this NMA dataset contained trials with direct comparisons for 4 of these; there was no trial with a direct comparison of DES and PTCA, or DES and medical therapy.

In their analysis, Trikalinos et al. used a CB model and assumed the relative effects were exchangeable. In particular, there were inconsistencies in the results based on the direct pairwise comparison and the results based on NMA, which included indirect evidence. The RR of death for PTCA vs. medical therapy evaluated in direct comparisons had an opposite direction compared to the indirect RR of death when evaluated in the NMA. The overall conclusion was that no evidence showed PCI treatments (including PTCA, BMS, and DES) were more effective than medical therapy in treating non-acute CAD.