Examples: bias adjustment thresholds for network meta-analysis

David Phillippo, University of Bristol

Introduction

This vignette provides worked examples for the nmathresh package, recreating exactly the analyses in the paper by Phillippo et al. (2018).

Thrombolytics

Caldwell et al. (2005) present a NMA of six thrombolytic treatments, with data taken from two systematic reviews (Boland et al. 2003, Keeley et al. 2003).
Thrombolytic treatment network, from Phillippo et al. (2018).

Thrombolytic treatment network, from Phillippo et al. (2018).

The results of the NMA are available in the nmathresh package, as Thrombo.post.summary (for the posterior summaries) and Thrombo.post.cov (for the posterior covariance matrix). The posterior summaries were generated using the coda package from saved WinBUGS chains and are stored as summary.mcmc objects, but the coda package is not required for our analysis.

library(nmathresh)

# library(coda)   # Not required - but prints the summary in a nicer format
# Thrombo.post.summary

Study level threshold analysis

To run a study level threshold analysis, we require the study data. This is available in nmathresh as a tab-delimited text file, and is read in like so:

dat.raww <- read.delim(system.file("extdata", "Thrombo_data.txt", package = "nmathresh"))

# print first few rows
head(dat.raww)
##    r.1   n.1  r.2   n.2 r.3   n.3 t.1 t.2 t.3 n_arms studyID
## 1 1472 20251  652 10396 723 10374   1   3   4      3       1
## 2    9   130    6   123  NA    NA   1   2  NA      2       2
## 3    5    63    2    59  NA    NA   1   2  NA      2       3
## 4    3    65    3    64  NA    NA   1   2  NA      2       4
## 5  887 10396  929 10372  NA    NA   1   2  NA      2       5
## 6 1455 13780 1418 13746  NA    NA   1   2  NA      2       6
n <- nrow(dat.raww)   # number of studies

Thresholds will be derived on the log odds ratio scale, so we derive log odds ratios against arm 1 as a reference.

# Log OR for two-arm trials, arm 2 vs. 1
dat.raww$lor.1 <- with(dat.raww, log(r.2 * (n.1 - r.1) / ((n.2 - r.2) * r.1)) )

dat.raww$k.1 <- dat.raww$t.2   # Arm 2 treatment
dat.raww$b.1 <- dat.raww$t.1   # Reference treatment


# Log OR for three-arm trials, arm 3 vs. 1
dat.raww$lor.2 <- with(dat.raww, log(r.3 * (n.1 - r.1) / ((n.3 - r.3) * r.1)) )

dat.raww$k.2 <- dat.raww$t.3   # Arm 3 treatment (NA if only 2 arms)
dat.raww$b.2 <- ifelse(is.na(dat.raww$k.2), NA, dat.raww$t.1)   # Reference treatment

The likelihood covariance matrix is then constructed in a block diagonal manner, with the aid of Matrix::bdiag.

# LOR variances and covariances, likelihood covariance matrix V
V.diag <- as.list(rep(NA, n))
attach(dat.raww)
for (i in 1:n){
  if (dat.raww$n_arms[i] == 2){
    V.diag[[i]] <- 1/r.1[i] + 1/r.2[i] + 1/(n.1[i]-r.1[i]) + 1/(n.2[i]-r.2[i])
  }
  else if (dat.raww$n_arms[i] == 3){
    v1 <- 1/r.1[i] + 1/r.2[i] + 1/(n.1[i]-r.1[i]) + 1/(n.2[i]-r.2[i])
    v2 <- 1/r.1[i] + 1/r.3[i] + 1/(n.1[i]-r.1[i]) + 1/(n.3[i]-r.3[i])
    # Covariance term
    c1 <- 1/r.1[i] + 1/(n.1[i] - r.1[i])
    V.diag[[i]] <- matrix(c(v1, c1, c1, v2), nrow = 2)
  }
}
detach(dat.raww)

library(Matrix)
V <- bdiag(V.diag)

The raw data was imported in wide format, with one row per study. It is much easier to work with the data in long format, with one row per data point (contrast).

# Reshape the data to have one row per contrast per study
dat.rawl <- reshape(dat.raww, varying = c("lor.1", "b.1", "k.1", "lor.2", "b.2", "k.2"), 
                    timevar = "c", idvar = "studyID", direction = "long")

# Sort data by study and contrast, removing NA rows
dat.rawl <- dat.rawl[order(dat.rawl$studyID, dat.rawl$c, dat.rawl$b, na.last = NA), ]

N <- nrow(dat.rawl)   # number of data points
K <- length(unique(c(dat.rawl$b, dat.rawl$k)))   # Number of treatments

Construct the design matrix.

# Construct the design matrix, with a row for each contrast and K-1 columns (parameters)
X <- matrix(0, nrow = N, ncol = K-1)

for (i in 1:N){
  X[i, dat.rawl$k[i]-1] <- 1
  if (dat.rawl$b[i] != 1){
    X[i, dat.rawl$b[i]-1] <- -1
  }
}

We are now ready to perform a threshold analysis at the study level. This is made easy using the nma_thresh function, which takes the posterior means and covariance matrix of the treatment effect parameters (\(d_k\)), the likelihood covariance matrix, and the design matrix. We specify nmatype = "fixed" to derive thresholds for the FE model, and opt.max = FALSE since the optimal treatment is the one which minimises the log odds.

# Now we can perform thresholding at the study level
thresh <- nma_thresh(mean.dk = Thrombo.post.summary$statistics[1:(K-1), "Mean"], 
                     lhood = V, 
                     post = Thrombo.post.cov, 
                     nmatype = "fixed",
                     X = X,
                     opt.max = FALSE)
## Likelihood for N = 15 data points.
## Number of treatments is K = 6.
## Current optimal treatment is k* = 3.

The nma_thresh function prints some basic details, which can be used to verify that the input was as expected (the number of data points and treatments, and the base-case optimal treatment). These are also shown when the threshold object is printed:

thresh
## A thresh object. For help, see ?'thresh-class'.
## 
## Base-case optimal treatment is k* = 3.
## 
##             lo lo.newkstar          hi hi.newkstar
## 1   -1.0318603           6   0.1197831           4
## 2   -0.1071759           4   5.2771792           5
## 3  -44.6044047           2 162.1861914           6
## 4 -111.2554454           2 404.5362135           6
## 5 -105.8657622           2 384.9387723           6
## 6   -0.3656429           2   1.3295151           6
## ... 9 further rows omitted ...

Finally, we will use the function thresh_forest to display the thresholds on a forest plot. We sort the rows of the plot to display those with smallest thresholds first; this is achieved using the orderby option. (By default, thresh_forest calculates recommended output dimensions, which are useful for saving to PDF. This isn’t necessary here, so we set calcdim = FALSE.)

# Display using a forest plot, along with 95% confidence intervals for LORs
# Create row labels
dat.rawl$lab <- rep(NA, nrow(dat.rawl))
for (i in 1:nrow(dat.rawl)) {
  dat.rawl$lab[i] <- paste0(dat.rawl$studyID[i], " (", dat.rawl$k[i], " vs. ", dat.rawl$b[i], ")")
}

# Calculate 95% CIs
dat.rawl$CI2.5 <- dat.rawl$lor + qnorm(.025)*sqrt(diag(V))
dat.rawl$CI97.5 <- dat.rawl$lor + qnorm(.975)*sqrt(diag(V))

# Calculate the proportion of CI covered by invariant interval, for sorting.
# Coverage <1 means that the CI extends beyond the bias invariant threshold, and 
# the threshold is below the level of statistical uncertainty.
dat.rawl$coverage <- apply(cbind(thresh$thresholds$lo / (dat.rawl$CI2.5 - dat.rawl$lor), 
                             thresh$thresholds$hi / (dat.rawl$CI97.5 - dat.rawl$lor)), 
                         1, min, na.rm = TRUE)


# Plot
thresh_forest(thresh, 
              y = lor, CI.lo = CI2.5, CI.hi = CI97.5, 
              label = lab, orderby = coverage, data = dat.rawl,
              CI.title = "95% Confidence Interval", y.title = "Log OR", 
              label.title = "Study (Contrast)", xlab = "Log Odds Ratio", 
              xlim = c(-3, 2), refline = 0, digits = 2,
              calcdim = FALSE)

For a more detailed explanation and interpretation of these forest plots see Phillippo et al. (2018). For example, the estimated log odds ratio of treatment 5 vs. treatment 1 in study 11 is -0.06. A negative bias adjustment of -0.14 would move the estimate to the lower bound of the decision-invariant bias adjustment interval (-0.20), at which point treatment 5 would replace treatment 3 as optimal. Conversely, a positive bias adjustment of 0.87 would move the estimate to the upper bound of the decision-invariant bias adjustment interval, at which point treatment 2 would replace treatment 3 as optimal. The lower portion of the invariant interval is shaded red (and the row label is bold) as the lower threshold lies within the 95% confidence interval of the study 11 estimate; the decision is sensitive to the level of imprecision in this estimate. The threshold analysis highlights that the decision is sensitive to bias adjustments in several studies (particularly 34, 35, and 11), but also shows a robustness to bias adjustments in many other studies with wide invariant intervals.

Thresholds for multiple biases

Continuing at study level, we can consider thresholds for multiple bias adjustments simultaneously. Here, we shall derive thresholds for simultaneous bias adjustments to the two contrasts of the 3-arm study 1, which are data points 1 and 2 in the dataset. The threshold lines can be derived directly from the set of values \(u_{ak^*,m}\), which is provided in matrix form by the output of nma_thresh in $Ukstar. It is straightforward to calculate the gradient - thresh$Ukstar[,2] / thresh$Ukstar[,1] and intercept thresh$Ukstar[,2] of every threshold line, however it is rather tedious to construct the plot by hand. The function thresh_2d does all the hard work for us, taking the threshold object thresh and the row numbers of the datapoints to consider as its main arguments:

thresh_2d(thresh, 1, 2,
          xlab = "Adjustment in Study 1 LOR: 3 vs. 1",
          ylab = "Adjustment in Study 1 LOR: 4 vs. 1",
          xlim = c(-1.5, 0.5), ylim = c(-2, 14),
          ybreaks = seq(-2, 14, 2))

We can see that, rather than requiring a single very large positive bias adjustment of 5.277 in the log OR of treatment 4 vs. 1 was needed to to change the optimal treatment to 5, two much smaller bias adjustments of 0.144 to the 3 vs. 1 log OR and 0.021 to the 4 vs. 1 log OR are also capable of crossing the invariant threshold to make treatment 5 optimal.

Contrast level threshold analysis

We can also perform a threshold analysis at the contrast level. This does not require the original data, only the joint posterior distribution of the treatment effect parameters.

First, we must reconstruct the likelihood covariance matrix that would have produced the posterior distribution in a FE 1-stage Bayesian NMA, where there was one data point for each contrast compared in one or more studies. To do this, we specify the design matrix of the contrasts with data (see the network diagram), and then use the function recon_vcov to reconstruct the likelihood covariance matrix.

K <- 6   # Number of treatments

# Contrast design matrix is
X <- matrix(ncol = K-1, byrow = TRUE,
            c(1, 0, 0, 0, 0,
              0, 1, 0, 0, 0,
              0, 0, 1, 0, 0,
              0, 0, 0, 1, 0,
              0, -1, 1, 0, 0,
              0, -1, 0, 1, 0,
              0, -1, 0, 0, 1))

# Reconstruct using NNLS
lik.cov <- recon_vcov(Thrombo.post.cov, prior.prec = .0001, X = X)
## Likelihood precisions found using NNLS.
## Residual Sum of Squares: 26.9748
##              --------------------
##               RSS fixed: 0.295837
##              RSS fitted: 26.679
##              --------------------
## Kullback-Leibler Divergence of fitted from 'true' posterior: 6.75957000427023e-05

The KL divergence is very small, which means that the reconstructed likelihood covariance matrix results in a posterior which is very close to that coming from the original NMA.

The function nma_thresh then calculates the thresholds.

thresh <- nma_thresh(mean.dk = Thrombo.post.summary$statistics[1:(K-1), "Mean"], 
                     lhood = lik.cov, 
                     post = Thrombo.post.cov, 
                     nmatype = "fixed", 
                     X = X, 
                     opt.max = FALSE)
## Likelihood for N = 7 data points.
## Number of treatments is K = 6.
## Current optimal treatment is k* = 3.

We now construct the forest plot using thresh_forest. The function d_ab2i is useful here, which quickly converts between the two ways of referencing contrasts: from \(d_{ab}\) style, used when writing or presenting contrasts, to d[i] style, used when storing and referencing vectors.

# Get treatment codes for the contrasts with data
d.a <- d.b <- vector(length = nrow(X))
for (i in 1:nrow(X)){
  d.a[i] <- ifelse(any(X[i, ] == -1), which(X[i, ] == -1), 0) + 1
  d.b[i] <- ifelse(any(X[i, ] == 1), which(X[i, ] == 1), 0) + 1
}

# Transform from d_ab style contrast references into d[i] style from the full set of contrasts
# for easy indexing in R
d.i <- d_ab2i(d.a, d.b, K = 6)

# Create plot data
plotdat <- data.frame(lab = paste0(d.b, " vs. ", d.a),
                      contr.mean = Thrombo.post.summary$statistics[d.i, "Mean"],
                      CI2.5 = Thrombo.post.summary$quantiles[d.i, "2.5%"],
                      CI97.5 = Thrombo.post.summary$quantiles[d.i, "97.5%"])

# Plot
thresh_forest(thresh, contr.mean, CI2.5, CI97.5, label = lab, data = plotdat,
              label.title = "Contrast", xlab = "Log Odds Ratio", CI.title = "95% Credible Interval",
              xlim = c(-.3, .3), refline = 0, digits = 2, calcdim = FALSE)

The contrast-level threshold analysis gives very similar results to the study-level threshold analysis - largely because the “combined evidence” on most contrasts is only a single study anyway. The threshold analysis gives very small thresholds for the combined evidence on treatment contrasts 5 vs. 3 and 6 vs. 3, which is symptomatic of the lack of evidence for significant differences between these treatments.

Social Anxiety

Forty-one treatments for social anxiety were compared in a network-meta analysis by Mayo-Wilson et al. (2014, also National Collaborating Centre for Mental Health 2013). A random effects model was used, and treatments were modelled within classes.
Social anxiety treatment network, from Phillippo et al. (2018).

Social anxiety treatment network, from Phillippo et al. (2018).

The results of the NMA are available in the nmathresh package, as SocAnx.post.summary (for the posterior summaries) and SocAnx.post.cov (for the posterior covariance matrix). The posterior summaries were generated using the coda package from saved WinBUGS chains and are stored as summary.mcmc objects, but the coda package is not required for our analysis.

library(nmathresh)
library(Matrix)

# library(coda)   # Not required - but prints the summary in a nicer format
# SocAnx.post.summary

Study level threshold analysis

For a study level analysis, we require the original study data. This is available in the nmathresh package, and is read in like so:

# Read study data
dat.raww <- read.delim(system.file("extdata", "SocAnx_data.txt", package = "nmathresh"))

# Print first few rows
head(dat.raww)
##   n.arms t.1 t.2 t.3 t.4 t.5       y.2    Var.2 y.3 Var.3 y.4 Var.4 y.5 Var.5
## 1      2   1   6  NA  NA  NA -0.916255 0.101295  NA    NA  NA    NA  NA    NA
## 2      2   1   7  NA  NA  NA -0.956381 0.151841  NA    NA  NA    NA  NA    NA
## 3      2   1   8  NA  NA  NA -0.768207 0.015674  NA    NA  NA    NA  NA    NA
## 4      2   1   8  NA  NA  NA -0.748858 0.063047  NA    NA  NA    NA  NA    NA
## 5      2   1   8  NA  NA  NA -1.248431 0.059678  NA    NA  NA    NA  NA    NA
## 6      2   1   8  NA  NA  NA -0.944181 0.035749  NA    NA  NA    NA  NA    NA
##    V
## 1 NA
## 2 NA
## 3 NA
## 4 NA
## 5 NA
## 6 NA
n <- nrow(dat.raww)   # Number of studies

# Turn wide study data into long with one row for each arm
dat.rawl <- reshape(dat.raww, varying=c("t.2","y.2","Var.2","t.3","y.3","Var.3",
                                        "t.4","y.4","Var.4","t.5","y.5","Var.5"),
                    timevar="arm", idvar="studyID", direction="long")

# Sort data by study and contrast, removing NA rows
dat.rawl <- dat.rawl[order(dat.rawl$studyID, dat.rawl$arm, dat.rawl$y, na.last=NA),]

K <- length(unique(c(dat.rawl$t.1, dat.rawl$t)))    # Number of treatments
N <- nrow(dat.rawl)    # Number of data points

As the posterior summaries contain several variables, we pick out the indices of those which we need for later.

# Get indices of d, sd, diff in the CODA data
vnames <- sub("(.*)\\[.*","\\1", rownames(SocAnx.post.summary$statistics))
ind.d <- which(vnames=="d")
ind.sd <- which(vnames=="sd")
ind.diff <- which(vnames=="diff")
ind.delta <- which(vnames=="delta")

We then construct the likelihood covariance matrix as a block diagonal matrix, with the aid of Matrix::bdiag.

# Create likelihood covariance matrix
V.diag <- as.list(rep(NA,n))
attach(dat.raww)
for (i in 1:n){
  if (n.arms[i] == 2){
    V.diag[[i]] <- Var.2[i]
  }
  else {
    V.diag[[i]] <- matrix(V[i], nrow=n.arms[i]-1, ncol=n.arms[i]-1)
    tempVar <- c(Var.2[i], Var.3[i], Var.4[i], Var.5[i])
    diag(V.diag[[i]]) <- tempVar[!is.na(tempVar)]
  }
}
detach(dat.raww)

lik.cov <- bdiag(V.diag)

Once this is done, it is a simple matter of using nma_thresh to derive thresholds. We specify nmatype = "random", as a random effects NMA was performed, and opt.max = FALSE as the optimal treatment minimises the SMD of symptoms of social anxiety. The posterior covariance matrix specified in post is the covariance matrix of the basic treatment parameters (\(d_k\)) and the random effects parameters (\(\delta_i\)).

thresh <- nma_thresh(mean.dk = SocAnx.post.summary$statistics[ind.d, "Mean"], 
                     lhood = lik.cov, 
                     post = SocAnx.post.cov,
                     nmatype = "random", 
                     opt.max = FALSE)
## Likelihood for N = 146 data points.
## Number of treatments is K = 41.
## Number of delta parameters is n.delta = 146.
## Current optimal treatment is k* = 41.

A forest plot is then constructed using thresh_forest. The full forest plot is very large, with 146 rows. To save space, we only display contrasts with thresholds smaller than 2 SMD. We achieve this easily by using the orderby option, which here specifies a call to the order function – ordering the rows by the variable ord (here derived as the coverage ratio of the invariant interval to the CI, so smallest thresholds first) and removing any NA rows with na.last = NA (here where thresholds are larger than 2 SMD).

# 95% CIs
dat.rawl$CI2.5 <- with(dat.rawl, y + qnorm(0.025)*sqrt(Var))
dat.rawl$CI97.5 <- with(dat.rawl, y + qnorm(0.975)*sqrt(Var))

# Study labels
dat.rawl$lab <- with(dat.rawl, paste0(studyID," (",t," vs ",t.1,")"))

# Forest plot - all contrasts, very large
# thresh_forest(thresh, y, CI2.5, CI97.5, label = lab, data = dat.rawl,
#               label.title = "Study (Contrast)", xlab = "Standardised Mean Difference", 
#               xlim = c(-4, 3), y.title = "SMD", refline = 0, digits = 2)

# Forest plot - only contrasts with thresholds <2 SMD
cutoff <- 2
absmin <- function(x) min(abs(x))
dat.rawl$coverage <- apply(thresh$thresholds[, c("lo", "hi")] / 
                             (dat.rawl[, c("CI2.5", "CI97.5")] - dat.rawl$y), 
                           1, absmin)
dat.rawl$ord <- ifelse(thresh$thresholds$lo > -cutoff | thresh$thresholds$hi < cutoff, 
                       dat.rawl$coverage, NA)

thresh_forest(thresh, y, CI2.5, CI97.5, label = lab, 
              orderby = list(ord, na.last = NA), data = dat.rawl,
              label.title = "Study (Contrast)", xlab = "Standardised Mean Difference", 
              xlim = c(-4, 3), y.title = "SMD", refline = 0, digits = 2, calcdim = FALSE)

The treatment decision is robust to bias adjustments in the vast majority of the study data. However, the decision shows sensitivity to the level of imprecision in two study data points (96, comparing treatments 41 vs. 2, and 81, comparing 36 vs. 2). The overall NMA results showed no evidence of a significant difference between treatments 41 and 36 (mean difference -0.12 (-0.59, 0.35)), which is borne out in this threshold analysis. Note that bias adjustments in general may be plausible beyond the range of the confidence interval; the results should therefore be interpreted in light of the magnitude and direction of possible bias.

Contrast level threshold analysis

We can also perform a threshold analysis on the social anxiety NMA at the contrast level. We do not need the original data for a contrast level analysis, we can proceed using only the joint posterior distribution of the treatment effect parameters.

First, we reconstruct the likelihood covariance matrix. Since there are many studies, we won’t construct the contrast design matrix by hand (though this is perfectly possible); instead, we use the treatment details from the study data.

trt.dat <- read.delim(system.file("extdata", "SocAnx_data.txt", package = "nmathresh"))[, 1:6]

head(trt.dat)   # Print first few rows
##   n.arms t.1 t.2 t.3 t.4 t.5
## 1      2   1   6  NA  NA  NA
## 2      2   1   7  NA  NA  NA
## 3      2   1   8  NA  NA  NA
## 4      2   1   8  NA  NA  NA
## 5      2   1   8  NA  NA  NA
## 6      2   1   8  NA  NA  NA
K <- with(trt.dat, length(unique(c(t.1, t.2, t.3, t.4, t.5)))) -1  # Number of treatments , -1 for NA

# work out which contrasts have data
contr.ab <- data.frame(a = c(), b = c())

for (i in 1:nrow(trt.dat)) {
   rowi <- trt.dat[i, which(!is.na(trt.dat[i, 2:6]))+1] # non NA elements of ith row
   
   # get contrast from all combinations of treatments
   trtcomb <- combn(rowi, 2, function(x) sapply(x, as.numeric))
   
   a <- apply(trtcomb, 2, min)
   b <- apply(trtcomb, 2, max)
   
   # remove contrasts of treatments against themselves
   iseq <- a == b
   a <- a[!iseq]
   b <- b[!iseq]
   
   if (!all(iseq)) contr.ab <- rbind(contr.ab, cbind(a, b))
}

contr.ab <- unique(contr.ab[order(contr.ab$a, contr.ab$b), ])


# Contrast design matrix
X <- matrix(0, nrow = nrow(contr.ab), ncol = K-1)
for (i in 1:nrow(X)) {
  if (contr.ab[i, "a"] > 1) X[i, contr.ab[i, "a"] - 1]  <- -1
  if (contr.ab[i, "b"] > 1)   X[i, contr.ab[i, "b"] - 1]    <- 1
}

As the posterior summaries contain several variables, we pick out the indices of those which we need.

# Get indices of d, sd, diff in the CODA data
vnames <- sub("(.*)\\[.*","\\1", rownames(SocAnx.post.summary$statistics))
ind.d <- which(vnames=="d")
ind.sd <- which(vnames=="sd")
ind.diff <- which(vnames=="diff")
ind.delta <- which(vnames=="delta")

The likelihood covariance matrix is then reconstructed using recon_vcov. The original NMA models treatments within classes, with informative gamma priors on the class precisions. We cannot incorporate this exactly into the estimation, but we make an approximation by setting the prior precision to the mean of the prior gamma distribution for all but treatment 3 (which had prior precision 0.0001).

# Class model is used, so use the prior mean precision from the gamma distribution
prior.prec <- rep(3.9/0.35, 40)
# Other than for treatment 3
prior.prec[2] <- 0.0001

lik.cov <- recon_vcov(SocAnx.post.cov[1:(K-1), 1:(K-1)], prior.vcov = diag(1/prior.prec), X = X)
## Likelihood precisions found using NNLS.
## Residual Sum of Squares: 735.679
##              --------------------
##               RSS fixed: 0.253083
##              RSS fitted: 735.425
##              --------------------
## Warning in recon_vcov(SocAnx.post.cov[1:(K - 1), 1:(K - 1)], prior.vcov =
## diag(1/prior.prec), : Returned some infinite variances. These will be not be
## included in KL calculation.
## Kullback-Leibler Divergence of fitted from 'true' posterior: 1.54967637538608

The KL divergence is 1.55, which indicates that the reconstructed likelihood covariance matrix results in a posterior that is reasonably close to the true posterior (values less than 1 indicate negligible differences, values greater than 3 indicate considerable differences). The small differences arise because the reconstructed likelihood is restricted to having independent data points (one for each contrast), whereas in reality there are multi-arm trials and a hierarchical class model, which cannot be fully characterised by the independent data points.

Now that we have the reconstructed likelihood covariance matrix, we derive the contrast level thresholds using nma_thresh with nmatype = "fixed".

thresh <- nma_thresh(mean.dk = SocAnx.post.summary$statistics[ind.d, "Mean"], 
                     lhood = lik.cov,
                     post = SocAnx.post.cov[1:(K-1), 1:(K-1)],
                     nmatype = "fixed", 
                     X = X, 
                     opt.max = FALSE)
## Likelihood for N = 84 data points.
## Number of treatments is K = 41.
## Current optimal treatment is k* = 41.

The results are presented on a forest plot using thresh_forest. Again, we only display contrasts with thresholds smaller than 2 SMD for brevity.

# Get indices of contrasts in likelihood
d.a <- d.b <- vector(length = nrow(X))
for (i in 1:nrow(X)) {
  d.a[i] <- ifelse(any(X[i, ] == -1), which(X[i, ] == -1), 0) + 1
  d.b[i] <- ifelse(any(X[i, ] == 1), which(X[i, ] == 1), 0) + 1
}

d.i <- d_ab2i(d.a, d.b, K = K)

# Contrast labels and credible intervals (from the posterior summary)
plotdat <- data.frame(
  contr = paste0(d.b, " vs. ", d.a),
  contr.mean = SocAnx.post.summary$statistics[ind.diff[d.i], "Mean"],
  CI2.5 = SocAnx.post.summary$quantiles[ind.diff[d.i], "2.5%"],
  CI97.5 = SocAnx.post.summary$quantiles[ind.diff[d.i], "97.5%"]
)

cutoff <- 2
absmin <- function(x) min(abs(x))
plotdat$ord <- ifelse(thresh$thresholds$lo > -cutoff | thresh$thresholds$hi < cutoff,
                      apply(thresh$thresholds[, c("lo", "hi")], 1, absmin), NA)

thresh_forest(thresh, contr.mean, CI2.5, CI97.5, 
              label = contr, orderby = list(ord, na.last = NA), data = plotdat,
              label.title = "Contrast", xlab = "SMD", xlim = c(-4, 3), 
              CI.title = "95% Credible Interval",
              refline = 0, digits = 2, calcdim = FALSE)

The treatment decision is robust to the level of imprecision in the combined data available on each contrast. A standardised mean difference of more than 0.8 may be considered large in the context of behavioural sciences (Cohen, 1988). All but five thresholds are larger than this, and for each of these the new optimal treatment at the threshold is 36. Rather than performing a long and laborious qualitative assessment of all 84 contrasts and 100 studies, attention can be focused on the smaller number of contrasts (for example the 5 studies with thresholds smaller than 0.8 SMD) where plausible adjustments to the data may cause a change in treatment recommendation. For a more detailed explanation and interpretation of these forest plots see Phillippo et al. (2018).

More complex analyses

More complex analyses are possible by manipulating the set of \(u_{ak^*,m}\) values, which are contained in the thresh object output by nma_thresh in the matrix Ukstar. For example, continuing with the contrast-level analysis of the social anxiety NMA, we could consider a common bias in all pharmacological treatments against inactive control, or all psychological treatments against inactive control.

Firstly, consider a common pharmacological treatment bias. The overall influence on the contrast between treatments \(a\) and \(k^*\) of a common adjustment to all drug data points is found simply by summing over the individual influences of each drug data point (since a single common adjustment is to be made to all efficacy estimates). Therefore the point where treatment \(a\) replaces \(k^*\) as optimal is found by summing over the individual threshold solutions \(u_{ak^*,m}\) for each drug data point \(m\).

# Drug treatments (+ combi therapies)
drugtrts <- c(9:23, 37:41)

# Which data points compare these to an inactive trt?
drugdats <- which(contr.ab$a %in% 1:3 & contr.ab$b %in% drugtrts)

# Get U solutions by summing the indivual solutions of drug data points
U.drugs <- 1 / (rowSums(1 / thresh$Ukstar[,drugdats]))

# Which contrasts do the rows of Ukstar correspond to?
Ukstar.ab <- d_i2ab(1:(K*(K-1)/2), K)
Ukstar.ab <- Ukstar.ab[Ukstar.ab$a == thresh$kstar | Ukstar.ab$b == thresh$kstar, ]

We then derive threshold values and new optimal treatments at each threshold by picking out the minimum positive value and maximum negative value from the overall drug threshold solution vector U.drugs. The function get.int makes this simple, with some additional handling of infinite thresholds behind the scenes.

# Thresholds are then
thresh.drugs <- get.int(U.drugs, thresh$kstar, 1:K, Ukstar.ab)

Now we plot the invariant interval, along with the pharmacological treatment estimates.

## Function to plot the common invariant interval with the data
plotII <- function(thresh, contr.mean, CrI.lo, CrI.hi, rowlabs, xlim, xlab, ylab, ...){
  
  yaxs <- length(contr.mean):1
  
  # split plot in two
  layout(matrix(1:2,nrow=1), widths=c(.2,1))
  
  # plot row labels on left side
  gp <- par(mar=c(5,4,1,0))
  plot(rep(0,length(yaxs)), yaxs, pch=NA, ylim=c(.5,yaxs[1]+.5), ylab="",
       xlab="", yaxt="n", xaxt="n", bty="n")
  text(0, yaxs, labels=rowlabs,xpd=NA)
  title(ylab=ylab, line=2)
  
  # fake plot for setup of right side
  par(mar=c(5,1,1,2))
  plot(contr.mean, yaxs, pch=NA, yaxt="n", xlim=xlim,
       ylim=c(.5,yaxs[1]+.5), ylab="", xlab="",...)
  title(xlab=xlab, line=2)
  
  # reference line
  abline(v=0, lty=2, col="gray")
  
  # combined invariant region
  polygon(rep(c(contr.mean + thresh$lo, rev(contr.mean) + thresh$hi), each=2), 
          c(rep(yaxs,each=2) + c(.5,-.5), rep(rev(yaxs),each=2) + c(-.5,.5)),
          col=rgb(.7,.8,.9,.7),
          border=rgb(.6,.7,.8))
  
  # credible intervals
  segments(y0=yaxs, x0=CrI.lo, x1=CrI.hi, lend=1)
  
  # contrast means
  points(contr.mean, yaxs, pch=21, col="black", bg="white")
  
  # write new optimal treatments at either side
  text(xlim[1], round(length(yaxs)/2), 
       labels=as.expression(substitute(tilde(k)*"* = "*xx,list(xx=thresh$lo.newkstar))), 
       pos=4)
  text(xlim[2], round(length(yaxs)/2),
       labels=as.expression(substitute(tilde(k)*"* = "*xx,list(xx=thresh$hi.newkstar))),
       pos=2)
  
  # write invariant interval below plot
  with(thresh, title(xlab=paste0("Invariant interval about zero: ",lo.newkstar,
                                 " (",formatC(lo,digits=2,format="f"),", ",
                                 formatC(hi,digits=2,format="f"),") ",
                                 hi.newkstar), line=3))
  
  par(gp)
}


plotII(thresh.drugs, 
       SocAnx.post.summary$statistics[ind.diff[drugdats], "Mean"], 
       SocAnx.post.summary$quantiles[ind.diff[drugdats], "2.5%"],
       SocAnx.post.summary$quantiles[ind.diff[drugdats], "97.5%"],
       rowlabs = paste0(contr.ab[drugdats, "b"], " vs. ", contr.ab[drugdats, "a"]),
       xlim = c(-4, 1.5), ylab = "Drug vs. Inactive Contrasts", xlab = "SMD")

Similarly for a common psychological treatment bias:

# Psych treatments
psychtrts <- c(4:8, 24:36)

# Which data points compare these to an inactive trt?
psychdats <- which(contr.ab$a %in% 1:3 & contr.ab$b %in% psychtrts)

# Get U solutions by summing the influences of drug data points
U.psych <- 1 / (rowSums(1 / thresh$Ukstar[,psychdats]))

# Thresholds are then
thresh.psych <- get.int(U.psych, thresh$kstar, 1:K, Ukstar.ab)

plotII(thresh.psych, 
       SocAnx.post.summary$statistics[ind.diff[psychdats], "Mean"], 
       SocAnx.post.summary$quantiles[ind.diff[psychdats], "2.5%"],
       SocAnx.post.summary$quantiles[ind.diff[psychdats], "97.5%"],
       rowlabs=paste0(contr.ab[psychdats,"b"]," vs. ",contr.ab[psychdats,"a"]),
       xlim=c(-3,2), ylab="Psych vs. Inactive Contrasts", xlab="SMD")

For a more detailed explanation and interpretation of these forest plots see Phillippo et al. (2018). The magnitude of the thresholds for both common pharmacological and common psychological treatment biases are large (Cohen, 1988). Any such biases - if they exist - are likely to be much smaller than these thresholds.

We may also assess the impact of adjusting for common biases of this nature simultaneously. Again, this is done by manipulating the set of \(u_{ak^*,m}\) values to derive threshold lines (the boundaries where treatment \(a\) replaces treatment \(k^*\) as optimal). With adjustment for common pharmacological treatment bias on the \(x\)-axis and adjustment for common psychological treatment bias on the \(y\)-axis, the \(y\) intercept of each line is the sum of all \(u_{ak^*,m}\) values over the psychological data points (i.e. the one-dimensional threshold solutions). Similarly, the \(x\) intercept of each line is the sum of all \(u_{ak^*,m}\) values over the pharmacological data points.

Instead of deriving the threshold lines by hand using U.drugs and U.psych, we will assemble a bare-bones thresh object for input to thresh_2d, to take advantage of the automated plotting routines. The Ukstar matrix has a column for the threshold solutions for the combined drug data points U.drugs, and a column for threshold solutions for the combined psychological data points U.psych. thresh_2d also uses kstar, so we include that too. Internally, thresh_2d uses the Ukstar matrix to derive the threshold lines with gradient - U.psych / U.drugs and intercept U.psych.

thresh.drugpsych <- list(
  Ukstar = cbind(U.drugs, U.psych),
  kstar = thresh$kstar
)

thresh_2d(thresh.drugpsych, 1, 2,
          xlab = "SMD adjustment: Drug vs. Inactive",
          ylab = "SMD adjustment: Psych vs. Inactive",
          xlim = c(-6, 2), ylim = c(-6, 2),
          breaks = -6:2)

References

Boland A, Dundar Y, Bagust A, Haycox A, Hill R, Mujica Mota R, et al. Early thrombolysis for the treatment of acute myocardial infarction: a systematic review and economic evaluation. Health technology assessment 2003;7:1-136.

Caldwell DM, Ades AE, Higgins JPT. Simultaneous comparison of multiple treatments: combining direct and indirect evidence. Brit Med J 2005;331:897-900.

Cohen J. Statistical Power Analysis for the Behavioral-Sciences. Percept Motor Skill 1988.

Keeley EC, Boura JA, Grines CL. Primary angioplasty versus intravenous thrombolytic therapy for acute myocardial infarction: a quantitative review of 23 randomised trials. Lancet 2003;361:13-20.

Mayo-Wilson E, Dias S, Mavranezouli I, Kew K, Clark DM, Ades AE, et al. Psychological and pharmacological interventions for social anxiety disorder in adults: a systematic review and network meta-analysis. Lancet Psychiatry 2014;1:368-76.

National Collaborating Centre for Mental Health. Social Anxiety Disorder: Recognition, Assessment and Treatment. Leicester and London: The British Psychological Society and the Royal College of Psychiatrists; 2013.

Phillippo DM, Dias S, Ades AE, Didelez V and Welton NJ. Sensitivity of treatment recommendations to bias in network meta-analysis. Journal of the Royal Statistical Society: Series A (Statistics in Society) 2018;181:843-867. DOI: 10.1111/rssa.12341

Appendix: WinBUGS code and reading CODA output

Thrombolytics example

The WinBUGS code and data for the thrombolytics example are provided by Caldwell et al. (2005, see supplementary materials). After running the model until convergence in WinBUGS, samples from the posterior are saved into text files in the CODA format. We read these in to R with the coda package, and derive the posterior summaries included in the nmathresh package.

# Use coda package to read in the CODA files generated by WinBUGS
# The CODA files need only contain the dd parameter (contrasts).
library(coda)
dat.CODA <- mcmc.list(lapply(c("Coda1.txt",
                               "Coda2.txt",
                               "Coda3.txt"),
                      read.coda, index.file = "CodaIndex.txt"))


# Posterior summary table
Thrombo.post.summary <- summary(dat.CODA)

# Posterior covariance matrix of basic treatment effects d_k = d_1k
Thrombo.post.cov <- cov(as.matrix(dat.CODA[,1:5]))

Social anxiety example

The WinBUGS code and data for the thrombolytics example are provided by Mayo-Wilson et al. (2014, see supplementary materials). After running the model until convergence in WinBUGS, samples from the posterior are saved into text files in the CODA format. We read these in to R with the coda package, and derive the posterior summaries included in the nmathresh package.

# Use coda package to read in the CODA files generated by WinBUGS
# The CODA files need only contain the d, delta, and sd parameters.
library(coda)
dat.CODA <- mcmc.list(lapply(c("Coda1.txt",
                               "Coda2.txt"),
                             read.coda, index.file = "CodaIndex.txt"))

# Get indices of parameters
vnames <- sub("(.*)\\[.*","\\1", varnames(dat.CODA))
ind.d <- which(vnames=="d")
ind.diff <- which(vnames=="diff")
ind.delta <- which(vnames=="delta")
ind.sd <- which(vnames=="sd")

# Posterior summary table
SocAnx.post.summary <- summary(dat.CODA)

# Posterior covariance matrix of d and delta parameters
SocAnx.post.cov <- cov(as.matrix(dat.CODA[,c(ind.d, ind.delta)]))