Sensitivity of treatment recommendations to bias in network meta‐analysis

Summary Network meta‐analysis (NMA) pools evidence on multiple treatments to estimate relative treatment effects. Included studies are typically assessed for risk of bias; however, this provides no indication of the impact of potential bias on a decision based on the NMA. We propose methods to derive bias adjustment thresholds which measure the smallest changes to the data that result in a change of treatment decision. The methods use efficient matrix operations and can be applied to explore the consequences of bias in individual studies or aggregate treatment contrasts, in both fixed and random‐effects NMA models. Complex models with multiple types of data input are handled by using an approximation to the hypothetical aggregate likelihood. The methods are illustrated with a simple NMA of thrombolytic treatments and a more complex example comparing social anxiety interventions. An accompanying R package is provided.


Proof of Lemma 1a
Recall We have that * k ɶ is the new optimal treatment under the bias-adjusted data, so by definition

Proof of Lemma 1b
Let us choose a value of the bias adjustment β so that the posterior expectation of only one contrast changes sign, say ( ) Following on from the proof of part a) we see that { } B b = , a singleton set. Since * k B ∈ ɶ by definition, we therefore know that * k b = ɶ .

Proof
Using Bayes rule (Bayes 1763), the posterior distribution is  Then since and are symmetric, and noting proportionality with ,   1  exp  2  2 And completing the square, which we recognise as the multivariate Normal as required. See also Gelman et al. (2013, p. 71).

A.3 Derivation of bias adjustment thresholds for the basic FE model
We derive thresholds for the basic FE model described in section 2.3.2 by rearranging the expression for the posterior mean from equation (8). Under the bias adjusted data, this becomes Notice here that we can think of  Furthermore, from Lemma 1 (appendix A.1) and equation (A.7) we know that the new optimal treatment at the threshold is * b k = ɶ . Note that, although it is in theory possible for there to be two or more treatment effects that are both exactly maximal at the threshold value of β , this happens with probability zero.

A.4 Theorem: Posterior distribution for RE model
Given the RE model

Proof
Similarly to the proof for the posterior distribution of the FE model, the posterior distribution is which we recognise as the multivariate Normal as required. See also Gelman et al. (2013, p. 582).

A.5 Bias adjustment thresholds for the basic RE model
The joint posterior distribution of d and δ (conditional on 2 τ Σ ) is given in equation (9)(with proof in appendix A.4). We partition the joint covariance matrix as Under bias-adjusted data m = + y y β ɶ , the joint posterior mean of d and δ becomes where we make use of a subscript colon notation to indicate subvectors (i.e. for a general vector x of length l , we define : , , ) ( for 1 p q l ≤ < ≤ ). Then following the same arguments as the basic FE case (given in appendix A.3), the influence matrix is 1 * − = H B V and the thresholds are given by equations (5) and (6).

A.6 Bias adjustment thresholds for the extended FE model
The extended FE model with parameter vector where the posterior covariance matrix is ( ) Then the joint posterior expectation of the basic treatment effect parameters d under the bias- where we make use of a subscript colon notation to denote subvectors (as in appendix A.5) and similarly to subset rows of a matrix.
From here we continue exactly as in the basic case (appendix A.3) to derive bias adjustment thresholds using the threshold equations (5) and (6), where the influence matrix is now

A.7 Bias adjustment thresholds for the extended RE model
We extend the RE model to include additional parameters µ , which are given a Normal prior distribution, and have an associated design matrix M . We also allow for random effects terms to only be included for certain data points, using a design matrix L , for example in the case of noting that, as in the basic RE case, we can partition the posterior covariance matrix into blocks.
Following the same arguments as before (appendix A.2), we see that the thresholds are given by equations (5) and (6) where the influence matrix is now

A.8 Bias adjustment thresholds for RE models including class effects
Further extending the RE model of section 2.3.4 to include class effects, we have where z is the vector of class effect parameters and the matrix Z is a design matrix for the treatment clasfses. Each row of Z corresponds to a treatment, which is assigned a class by a 1 in the corresponding column and zeros elsewhere in the row. As before, we must work with the posterior distribution assuming that 2 τ Σ is known, fixed, and invariant to bias adjustments; furthermore, we now assume that the between-treatment covariance matrix d Σ is also known, fixed, and invariant to bias adjustments. Both of these assumptions can be tested using sensitivity analyses.
The posterior distribution is then noting that, as in the basic RE case, we can partition the posterior covariance matrix into blocks.
Proof of the posterior distribution follows closely that of the basic RE model in appendix A.4.
The joint posterior expectation of the basic treatment effect parameters d under the bias-adjusted data y ɶ is then Following the same arguments as before (appendix A.3), we see that the thresholds are given by equations (5) and (6) Equating this with the true joint posterior distribution ( ) N , η Σ reported by the original NMA, we see that In order to calculate bias adjustment thresholds using the results of section 2.3.2, we require the influence matrix The KL divergence is always non-negative and smaller values are desirable, indicating that the approximation reconstructs the posterior distribution well. Noting that the KL divergence is equivalent to the expected value of a log Bayes factor, we refer to Kass and Raftery (1995) for interpretation: for example, a KL divergence less than 1 is negligible, and values greater than 3 may be considered large. If evaluation of the KL divergence suggests a bad approximation, the contrast-level threshold analysis may give inaccurate results.
Once the hypothetical likelihood covariance matrix has been reconstructed, the thresholds are then evaluated as before using equations (5) and (6)  nmatype parameter (either "fixed" or "random" respectively). This function takes the posterior mean of the relative treatment effects, the likelihood and posterior covariance matrices, and the design matrix/matrices as inputs. The R function nma_thresh is also used to perform contrast-level analysis (section 2.4). In this scenario the hypothetical likelihood covariance matrix is constructed from the prior and posterior covariance matrices and the design matrix using the R function recon_vcov, either exactly or using non-negative least squares (NNLS) (Lawson and Hanson 1995) via the function nnls from the package nnls (Mullen and van Stokkum 2012). The R function thresh_forest takes the results from the NMA and threshold analysis and presents them graphically on a forest plot.

A.10.1 Mathematical derivation
Here we briefly describe mathematically the derivation of bias adjustment thresholds in a vectorised manner, allowing for highly efficient computation that does not rely on looping.
First, define a matrix U which contains the elements * , The NNLS problem is then