---
title: "Examples: bias adjustment thresholds for network meta-analysis"
author: "David Phillippo, University of Bristol"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{"Examples: bias adjustment thresholds for network meta-analysis"}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# 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).
```{r, echo=FALSE, out.width='60%', fig.cap="Thrombolytic treatment network, from Phillippo et al. (2018)."}
knitr::include_graphics("Thrombo_network.png")
```

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.
```{r}
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:
```{r}
dat.raww <- read.delim(system.file("extdata", "Thrombo_data.txt", package = "nmathresh"))

# print first few rows
head(dat.raww)

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.
```{r}
# 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`.
```{r}
# 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).
```{r}
# 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.
```{r}
# 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.
```{r}
# 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)
```
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:
```{r}
thresh
```

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`.)
```{r, fig.width=13, fig.height=5.5, out.width='100%', dpi=300}
# 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:
```{r, fig.width=6, fig.height=6, out.width='60%', dpi=300}
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
```{r, include=FALSE}
# Reset environment
rm(list = ls())
```
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.
```{r}
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)
```
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.
```{r}
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)
```

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.
```{r, fig.width = 12, fig.height = 3.5, out.width = '100%', dpi = 300}
# 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
```{r, include=FALSE}
# Reset environment
rm(list=ls())
```
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.
```{r, echo=FALSE, out.width='60%', fig.cap="Social anxiety treatment network, from Phillippo et al. (2018)."}
knitr::include_graphics("Social_Anxiety_network.png")
```

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.
```{r, eval=FALSE}
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:
```{r}
# Read study data
dat.raww <- read.delim(system.file("extdata", "SocAnx_data.txt", package = "nmathresh"))

# Print first few rows
head(dat.raww)

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.
```{r}
# 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`.
```{r}
# 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$).
```{r}
thresh <- nma_thresh(mean.dk = SocAnx.post.summary$statistics[ind.d, "Mean"], 
                     lhood = lik.cov, 
                     post = SocAnx.post.cov,
                     nmatype = "random", 
                     opt.max = FALSE)
```

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).
```{r, fig.width=14, fig.height=9.5, out.width='100%', dpi=300}
# 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
```{r, include = FALSE}
# Reset environment
rm(list=ls())
```
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.
```{r}
trt.dat <- read.delim(system.file("extdata", "SocAnx_data.txt", package = "nmathresh"))[, 1:6]

head(trt.dat)   # Print first few rows

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.
```{r}
# 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).
```{r}
# 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)
```

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"`.
```{r}
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)
```

The results are presented on a forest plot using `thresh_forest`. Again, we only display contrasts with thresholds smaller than 2 SMD for brevity.
```{r, fig.width=13, fig.height=11, out.width='100%', dpi=300}
# 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$.
```{r}
# 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.
```{r}
# 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.
```{r, fig.height=6, fig.width=8, out.width='80%', dpi=300}
## 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:
```{r, fig.width=8, fig.height=5.5, out.width='80%', dpi=300}
# 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`.
```{r, fig.width=8, fig.height=8, out.width='80%', dpi=300}
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](http://dx.doi.org/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.
```{r, eval=FALSE}
# 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.
```{r, eval=FALSE}
# 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)]))
```