This vignette provides worked examples for the `nmathresh`

package, recreating exactly the analyses in the paper by 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
```

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
```

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:

```
## 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.

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.

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.