Bayesian inference for the Bivariate model for diagnostic meta-analysis

Contents

OVERVIEW

MOTIVATING EXAMPLE AND DATA

THE MODEL

INITIAL VALUES

COMPILING AND RUNNING THE MODEL

MONITORING CONVERGENCE

SUMMARY STATISTICS

SROC plot

SCRIPT

REFERENCES

AUTHORS: Ian Schiller and Nandini Dendukuri

OVERVIEW

We are often requested for the programs we have written for estimating diagnostic meta-analysis models. We felt a blog article would be a good way to share them.  In this article we describe how to carry out Bayesian inference for the Bivariate meta-analysis model1,2 when a perfect reference test is available.

Below we explain each step of a program in the  R language3 (Link to whole program).The program uses the rjags package (an interface to the JAGS4 library) in order to draw samples from the posterior distribution of the parameters using Markov Chain Monte Carlo (MCMC) methods. It relies on the mcmcplots package to assess convergence of the meta-analysis model.  It returns summary statistics on the mean sensitivity and specificity across studies and between-study heterogeneity. A summary plot in ROC space can be created with the help of the DTAplots package.

MOTIVATING EXAMPLE AND DATA

We illustrate the script by applying it to a data from a systematic review of studies evaluating the accuracy of Anti-cyclic citrullinated peptide (Anti-CCP) test for rheumatoid arthritis (index test) compared to the 1987 revised American college or Rheumatology (ACR) criteria or clinical diagnosis (reference test) for 37 studies5.  Each study provides a two-by-two cell of the form

Perfect reference test
Test under evaluation
TP FP
FN TN

where,

TP = Number of individuals who tested positive on both tests.

FP = Number of individuals who tested positive on the index test and negative on the reference test

FN = Number of individuals who tested negative on the index test and positive on the reference test

TN = Number of individuals who tested negative on both tests.

The data can be read into R as follows:

TP=c(115, 110, 40, 23, 236, 74, 89, 90, 31, 69, 25, 43, 70, 167, 26, 110, 26, 64, 71, 68, 38, 42, 149, 147, 47, 24, 40, 171, 72, 7, 481, 190, 82, 865, 139, 69, 90)

FP=c(17, 24, 5, 0, 20, 11, 4, 2, 0, 8, 2, 1, 5, 8, 8, 3, 1, 14, 2, 14, 3, 2, 7, 10, 7, 3, 11, 26, 14, 2, 23, 12, 13, 79, 7, 5, 7)

FN=c(16, 86, 58, 7, 88, 8, 29, 50, 22, 18, 10, 63, 17, 98, 15, 148, 20, 15, 58, 35, 0, 60, 109, 35, 20, 18, 46, 60, 77, 9, 68, 105, 71, 252, 101, 107, 101)

TN=c(73, 215, 227, 39, 231, 130, 142, 129, 75, 38, 40, 120, 228, 88, 15, 118, 56, 293, 66, 132, 73, 96, 114, 106, 375, 79, 146, 443, 298, 51, 185, 408, 301, 2218, 464, 133, 313)

pos= TP + FN
neg= TN + FP
n <- length(TP) #  Number of studies
dataList = list(TP=TP,TN=TN,n=n,pos=pos,neg=neg)

THE MODEL

First, a bit about the theory behind the hierarchical meta-analysis model.

Within-study level: In each study, TP is assumed to follow a Binomial distribution with probability equal to the sensitivity and TN is assumed to follow a Binomial distribution with probability equal to the specificity.

Between-study level: The logit transformed sensitivity \mu_{1i} and the logit transformed specificity \mu_{2i} in each study jointly follow a bivariate normal distribution:

\left(\begin{matrix}\mu_{1i}\\\mu_{2i}\\\end{matrix}\right)\sim N\left(\left(\begin{matrix}\mu_1\\\mu_2\\\end{matrix}\right),\ \sum_{12}\right)  with \sum_{12}=\left(\begin{matrix}\tau_1^2&\rho\tau_1\tau_2\\\rho\tau_1\tau_2&\tau_2^2\\\end{matrix}\right)

where  \tau_1^2 and \tau_2^2 are the between-study variance parameters and  \rho  is the correlation between \mu_{1i}  and \mu_{2i}.    

To carry out Bayesian estimation, we need to specify prior distribution over the unknown parameters, e.g . \mu_1, \mu_2, \tau_1^2, \tau_2^2, and \rho. Instead of providing prior distribution over the between-study variance parameters \tau_1^2 and \tau_2^2 we instead placed prior distributions over the between-study precision parameters, \frac{1}{\tau_1^2} and \frac{1}{\tau_2^2}, respectively.

The model is written in JAGS syntax and stored in a character object we call modelString which is later saved in the current working directory via the R function writeLines as follows.

modelString =
  
  "model {
#=== LIKELIHOOD ===#
for(i in 1:n) {
TP[i] ~ dbin(se[i],pos[i])
TN[i] ~ dbin(sp[i],neg[i])

# === PRIOR FOR INDIVIDUAL LOGIT SENSITIVITY AND SPECIFICITY === #
logit(se[i]) <- l[i,1]
logit(sp[i]) <- l[i,2]
l[i,1:2] ~ dmnorm(mu[], T[,])
}

#=== HYPER PRIOR DISTRIBUTIONS OVER MEAN OF LOGIT SENSITIVITY AND SPECIFICITY === #
mu[1] ~ dnorm(0,0.3)
mu[2] ~ dnorm(0,0.3)
# Between-study variance-covariance matrix

T[1:2,1:2]<-inverse(TAU[1:2,1:2])
TAU[1,1] <- tau[1]*tau[1]
TAU[2,2] <- tau[2]*tau[2]
TAU[1,2] <- rho*tau[1]*tau[2]
TAU[2,1] <- rho*tau[1]*tau[2]

#=== HYPER PRIOR DISTRIBUTIONS FOR PRECISION OF LOGIT SENSITIVITY ===#
#=== AND LOGIT SPECIFICITY, AND CORRELATION BETWEEN THEM === #
prec[1] ~ dgamma(2,0.5)
prec[2] ~ dgamma(2,0.5)
rho ~ dunif(-1,1)

# ===  PARAMETERS OF INTEREST === #
# BETWEEN-STUDY VARIANCE (tau.sq) AND STANDARD DEVIATION (tau) OF LOGIT SENSITIVITY AND SPECIFICITY
tau.sq[1]<-pow(prec[1],-1)
tau.sq[2]<-pow(prec[2],-1)
tau[1]<-pow(prec[1],-0.5)
tau[2]<-pow(prec[2],-0.5)
# MEAN SENSITIVITY AND SPECIFICITY
Mean_Se <- 1/(1+exp(-mu[1]))
Mean_Sp <- 1/(1+exp(-mu[2]))

# PREDICTED SENSITIVITY AND SPECIFICITY IN A NEW STUDY
l.predicted[1:2] ~ dmnorm(mu[],T[,])
Predicted_Se <- 1/(1+exp(-l.predicted[1]))
Predicted_Sp <- 1/(1+exp(-l.predicted[2]))
}"

writeLines(modelString,con="model.txt")

INITIAL VALUES

Running multiple chains, each starting at different initial values, serves to assess convergence of the MCMC algorithm by examining if all chains reach the same solution.   You can use randomly generated initial values provided by rjags via the GenInits() function. In the illustration below we provide a bit more detail on how to provide initial values chosen by the user. 

Providing relevant initial values can speed up convergence of the MCMC algorithm.  Here we create a function that allows us to provide our own initial values for the logit-transformed sensitivity (mu[1]) and logit-transformed specificity (mu[2]) parameters.  By setting the chain argument in GenInits() equal to “myinits” we can provide specific initial values to mu[1] and mu[2] through the init_mu argument  or we can use the function to randomly generate initial values from distributions of our choice for mu[1] and mu[2].  In this example, we have allowed other parameters, such as between-study precisions and correlation between logit-transformed sensitivity and logit-transformed specificity, to be randomly generated by GenInits().  In our experience, choosing initial values for these parameters is not easy.  For reproducibility, a Seed argument is provided.

In the illustration below, we set mu[1] and mu[2] equal to 0 (which represents an initial value of 0.5 on the probability scale) for the first chain.  For chains number 2 and 3, we simply set the argument chain to 2 and 3, respectively.  This ensures both chain 2 and chain 3 are randomly initialized.

# Initial values
GenInits = function(chain, init_mu=c(0,0), Seed=seed){    
  if(chain == "myinits"){
    mu = c(as.numeric(init_mu[1]),as.numeric(init_mu[2]))
  }
  else{
    mu=c(rnorm(1,0,1/sqrt(0.3)),rnorm(1,0,1/sqrt(0.3)))
  }
  list(
    mu=mu,
    prec = c(rgamma(1,2,0.5),rgamma(1,2,0.5)),
    rho = runif(1,-1,1),
    .RNG.name="base::Wichmann-Hill",
    .RNG.seed=Seed
  )
}
Listinits = list(
  GenInits(chain="myinits", init_mu=c(0,0)),
  GenInits(chain=2),
  GenInits(chain=3)
)

COMPILING AND RUNNING THE MODEL

The jags.model function compiles the model and checks for possible syntax errors.  Argument n.chains=3 creates 3 MCMC chains with different starting values given by the Listinits object defined above.  

library(rjags)
# Compile the model
jagsModel = jags.model("model.txt",data=dataList,n.chains=3,inits=Listinits)

The update function serves to carry out the “burn-in” step.  This corresponds to discarding initial iterations till the MCMC algorithm has reached its equilibrium distribution.  Here we chose to discard the first 5,000 iterations by providing the argument n.iter=5000.  This can be helpful to mitigate “bad” initial values. 

# Burn-in iterations
update(jagsModel,n.iter=5000)

The coda.samples function draws samples from the posterior distribution of the unknown parameters that were defined in the parameters object.  Here we are drawing n.iter=15000 iterations from the posterior distributions for each chain.

# Parameters to be monitored
parameters = c("se","sp",  "mu", "tau.sq", "tau", "rho",  "Mean_Se" , "Mean_Sp", "Predicted_Se", "Predicted_Sp")

# Posterior samples
output = coda.samples(jagsModel,variable.names=parameters,n.iter=15000)

MONITORING CONVERGENCE

Convergence of the MCMC algorithm can now be examined.   The code below will create a 3-panel plot consisting of the density plot, the running mean plot and a history plot for parameters that vary across studies, e.g. sensitivity and specificity in each study.  

library(mcmcplots)
# Plots to check convergence for parameters varying across studies:
# Index j refers to the studies
# Index n refers to the number of studies
# Index i follows the ordering of the parameters vector defined in the previous box.  If left unchanged:
# i = 1: sensitivity in study j (se)
# i = 2: specificity in study j (sp)

for(j in 1:n) {
for(i in c(1,2)) {
  #  tiff(paste(parameters[i],"[",j,"] Converged",".tiff",sep=""),width = 23, height = 23, units = "cm", res=200)
  par(oma=c(0,0,3,0))
  layout(matrix(c(1,2,3,3), 2, 2, byrow = TRUE))
  denplot(output, parms=c(paste(parameters[i],"[",j,"]",sep="")), auto.layout=FALSE, main="(a)", xlab=paste(parameters[i],"",sep=""), ylab="Density")
  rmeanplot(output, parms=c(paste(parameters[i],"[",j,"]",sep="")), auto.layout=FALSE, main="(b)")
  title(xlab="Iteration", ylab="Running mean")
  traplot(output, parms=c(paste(parameters[i],"[",j,"]",sep="")), auto.layout=FALSE, main="(c)")
  title(xlab="Iteration", ylab=paste(parameters[i],"[",j,"]",sep=""))
  mtext(paste("Diagnostics for ", parameters[i],"[",j,"]","",sep=""), side=3, line=1, outer=TRUE, cex=2)
  #  dev.off()
}
}

Then we create the same 3-panel plot for the parameters that are shared across studies.

# Plots to check convergence for parameters shared across studies:
# Index j = 1 refers to sensitivity and j = 2 refers to specificity
# Index i follows the ordering of the parameters vector defined above.  If left unchanged :
# i = 3 : logit-transformed mean sensitivity of index test (j=1) and logit-transformed mean specificity of index test (j=2) (mu)
# i = 4 : between-study variance in logit-transformed sensitivity of index test (j=1) and between-study variance in logit-transformed specificity of index test (j=2) (tau.sq)
# i = 5 : between-study standard deviation in logit-transformed sensitivity of index test (j=1) and between-study standard deviation in logit-transformed specificity of index test (j=2) (tau)
for(j in 1:2){
for(i in c(3,4,5)) {
  #  tiff(paste(parameters[i],"[",j,"] Converged",".tiff",sep=""),width = 23, height = 23, units = "cm", res=200)
  par(oma=c(0,0,3,0))
  layout(matrix(c(1,2,3,3), 2, 2, byrow = TRUE))
  denplot(output, parms=c(paste(parameters[i],"[",j,"]",sep="")), auto.layout=FALSE, main="(a)", xlab=paste(parameters[i],"",sep=""), ylab="Density")
  rmeanplot(output, parms=c(paste(parameters[i],"[",j,"]",sep="")), auto.layout=FALSE, main="(b)")
  title(xlab="Iteration", ylab="Running mean")
  traplot(output, parms=c(paste(parameters[i],"[",j,"]",sep="")), auto.layout=FALSE, main="(c)")
  title(xlab="Iteration", ylab=paste(parameters[i],"[",j,"]",sep=""))
  mtext(paste("Diagnostics for ", parameters[i],"[",j,"]","",sep=""), side=3, line=1, outer=TRUE, cex=2)
  #  dev.off()
}
}

# Index j not needed for the following shared parameters
# i = 6 : correlation between logit-transformed sensitivity of index test and logit-transformed specificity of index test (rho)
# i = 7 : mean sensitivity of index test (Mean_Se)
# i = 8 : mean specificity of index test (Mean_Sp)
# i = 9: predicted sensitivity of index test in a future study (Predicted_Se)
# i = 10: predicted specificity of index test in a future study (Predicted_Sp)


for(i in 6:10) {
  #  tiff(paste(parameters[i]," Converged",".tiff",sep=""),width = 23, height = 23, units = "cm", res=200)
  par(oma=c(0,0,3,0))
  layout(matrix(c(1,2,3,3), 2, 2, byrow = TRUE))
  denplot(output, parms=c(paste(parameters[i], sep="")), auto.layout=FALSE, main="(a)", xlab=paste(parameters[i],"",sep=""), ylab="Density")
  rmeanplot(output, parms=c(paste(parameters[i], sep="")), auto.layout=FALSE, main="(b)")
  title(xlab="Iteration", ylab="Running mean")
  traplot(output, parms=c(paste(parameters[i], sep="")), auto.layout=FALSE, main="(c)")
  title(xlab="Iteration", ylab=paste(parameters[i], sep=""))
  mtext(paste("Diagnostics for ", parameters[i], sep=""), side=3, line=1, outer=TRUE, cex=2)
  #  dev.off()
}

For example, in the figure below, the posterior densities (a), running posterior mean values (b) and history plots (c) for the parameter Mean_Se (mean sensitivity of the Anti-CCP test) are displayed for each chain represented by a different color.  The very similar results from all three chains suggest the algorithm has converged to the same solution in each case.

SUMMARY STATISTICS

Once convergence of the MCMC algorithm is achieved, summary statistics of the parameters of interest may be extracted from their posterior distributions. 

# Summary statistics
summary(output)

Below is the summary statistics output for the parameters shared across studies. The 2.5% and 97.5% quantiles for the Predicted_Se and Predicted_Sp can be used to define prediction intervals for the mean sensitivity and mean specificity, respectively. If the prediction interval is much wider than the credible interval it would be indicative of high between-study heterogeneity.

Iterations = 5001:20000
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 15000

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                Mean       SD  Naive SE Time-series SE
Mean_Se       0.6580 0.028561 1.346e-04      2.163e-04
Mean_Sp       0.9561 0.006013 2.835e-05      6.614e-05
Predicted_Se  0.6421 0.155076 7.310e-04      7.547e-04
Predicted_Sp  0.9446 0.043617 2.056e-04      2.076e-04
mu[1]         0.6570 0.127495 6.010e-04      9.671e-04
mu[2]         3.0892 0.143139 6.748e-04      1.592e-03
rho          -0.4202 0.159652 7.526e-04      2.149e-03
tau.sq[1]     0.5546 0.151128 8.725e-04      2.040e-03
tau.sq[2]     0.5699 0.186720 1.078e-03      3.853e-03

2. Quantiles for each variable:

                2.5%     25%     50%     75%    97.5%
Mean_Se       0.6009  0.6391  0.6583  0.6772  0.71320
Mean_Sp       0.9433  0.9523  0.9564  0.9602  0.96694
Predicted_Se  0.3042  0.5402  0.6583  0.7600  0.89310
Predicted_Sp  0.8284  0.9298  0.9563  0.9731  0.99002
mu[1]         0.4091  0.5712  0.6559  0.7410  0.91098
mu[2]         2.8110  2.9934  3.0881  3.1828  3.37582
rho          -0.6945 -0.5356 -0.4334 -0.3192 -0.07313
tau.sq[1]     0.3250  0.4470  0.5313  0.6335  0.90124
tau.sq[2]     0.2898  0.4327  0.5337  0.6604  0.99862

SROC plot

The output from rjags can be used to obtain a plot in the summary receiver operating characteristic (SROC) space with 95% credible region (red line) and 95% prediction region (black dotted line) around the summary estimates of sensitivity and specificity (solid black circle) using the SROC_rjags function found in the DTAplots package.  Individual studies can also be represented (blue circles) and are proportional to their sample size.

# SROC plot
library(DTAplots)
SROC_rjags(X=output, model="Bivariate",n=n, study_col1="blue", study_col2=rgb(0, 0, 1, 0.15), dataset=cbind(TP,FP,FN,TN), ref_std=TRUE, Sp.range=c(0,1))

SCRIPT

For a copy of the complete script, please follow the link below

Link to whole script

REFERENCES

  1. Reitsma, J.B., Glas, A. S., Rutjes, A. W., Scholten, R. J., Bossuyt, P. M., and Zwinderman, A. H. (2005).  Bivariate analysis of sensitivity and specificity produces informative summary measures in diagnostic reviews.  Journal of Clinical Epidemiology 58, 892-990.
  • Chu H, Chen S, Louis TA. Random effects models in a meta‐analysis of the accuracy of two diagnostic tests without a gold standard. Journal of the American Statistical Association 2009;104(486):512‐23
  • R Core Team (2020). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
  • Plummer, Martyn. (2003). JAGS: A Program for Analysis of Bayesian Graphical Models using Gibbs Sampling. 3rd International Workshop on Distributed Statistical Computing (DSC 2003); Vienna, Austria. 124.
  • Nishimura K, Sugiyama D, Kogata Y, Tsuji G, Nakazawa T, Kawano S, Saigo K, Morinobu A, Koshiba M, Kuntz KM, Kamae I, Kumagai S. Meta-analysis: diagnostic accuracy of anti-cyclic citrullinated peptide antibody and rheumatoid factor for rheumatoid arthritis. Ann Intern Med. 2007 Jun 5;146(11):797-808. doi: 10.7326/0003-4819-146-11-200706050-00008. PMID: 17548411.
  • Cochrane handbook chapter

Add Your Comment

Your email address will not be published. Required fields are marked *