## OVERVIEW

This is a third article in our series of blog articles illustrating Rjags programs for diagnostic meta-analysis models. We previously looked at the Bayesian bivariate model and the Bayesian HSROC model. In this article we describe how to carry out Bayesian inference for the Bivariate latent class meta-analysis model (Xie, Sinclair and Dendukuri(2017)) when a perfect reference test is not available.

This blog will guide you through each step of a program in the R language (R Core Team(2020)). The program uses the rjags package (an interface to the JAGS library (Plummer(2003))) in order to draw samples from the posterior distribution of the parameters using Markov Chain Monte Carlo (MCMC) methods. The program 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

To illustrate the method, we apply the script to a data from a systematic review of studies evaluating the accuracy of GeneXpertTM (Xpert) test for tuberculosis (TB) meningitis (index test) and culture (reference test) for 29 studies (Kohli et al(2018)). Each study provides a two-by-two table of the form

reference_test_positive Reference_test_negative
Index test positive t11 t10
Index test negative t01 t00

where:

• t11 is the number of individuals who tested positive on both Xpert and culture tests
• t10 is the number of individuals who tested positive on Xpert and negative on culture test
• t01 is the number of individuals who tested negative on Xpert and positive on culture test
• t00 is the number of individuals who tested negative on both Xpert and culture tests
t11 <- c(7, 6, 2, 5, 0, 1, 4, 3, 3, 1, 7, 35, 103, 1, 2, 22, 6, 25, 27, 15, 5, 2, 2, 11, 2, 0, 8, 3, 1)
t10 <- c(5, 4, 0, 0, 0, 0, 6, 0, 3, 0, 0, 5, 6, 1, 0, 1, 1, 3, 11, 3, 9, 3, 0, 2, 4, 0, 2, 0, 2)
t01 <- c(5, 4, 0, 1, 2, 0, 8, 1, 1, 0, 0, 1, 18, 2, 0, 5, 0, 20, 25, 7, 10, 1, 1, 2, 2, 3, 5, 0, 0)
t00 <- c(63, 115, 2, 44, 132, 4, 83, 250, 67, 14, 150, 119, 252, 107, 8, 31, 148, 687, 204, 205, 115, 53, 4, 118, 22, 16, 186, 28, 43)
cell <- cbind(t11, t10, t01, t00)
n <- length(t11)  #  Number of studies
n.study <- rowSums(cell) #  Study sample sizes
dataList <- list(t12=cell,n=n, n.study=n.study)

## THE MODEL

Before moving any further, let’s look at a bit of theory behind the latent class meta-analysis model. In a latent class model, the target condition is assumed unknown.

Within-study level: In each study, the entries of the 2x2 table defined in the previous section are assumed to follow a multinomial distribution with probabilities p11, p10, p01 and p00 expressed as follows :

$(t11, t10, t01, t00) \sim \text{Multinomial}((p11, p10, p01, p00), n.study) \text{, where} \\ p11 = prev*(se1*se2 + covp) + (1-prev)*((1-sp1)*(1-sp2) + covn) \\ p10 = prev*(se1*(1-se2) - covp) + (1-prev)*((1-sp1)*sp2 - covn) \\ p01 = prev*((1-se1)*se2 - covp) + (1-prev)*(sp1*(1-sp2) - covn) \\ p00 = prev*((1-se1)*(1-se2) + covp) + (1-prev)*(sp1*sp2 + covn) \\ -(1-s1)*(1-s2) < covp < max(s1,s2) – s1*s2 \\ -(1-c1)*(1-c2) < covn < max(c1,c2) – c1*c2$

where prev denotes the prevalence of the latent target condition, se1 and sp1 denote the sensitivity and specificity of the index test, se2 and sp2 denote the sensitivity and specificity of the reference test, covp is the covariance between the tests among those who have the target condition, covn is the covariance between the tests among those who do not have the target condition and n.study is the total sample size in all four cells. The covariance terms adjust for the possibility that both tests make the same false positive or false negative errors. This is referred to as conditional dependence.

Between-study level: For the index test, the logit transformed sensitivity ($$\mu_{1i}$$) and the logit transformed specificity ($$\mu_{1i}$$) in each study jointly follow a bivariate normal distribution:

$\left(\begin{array}{rr} \mu_{1i} \\ \mu_{2i} \end{array}\right) \sim N \left( \left(\begin{array}{rr} \mu_{1i} \\ \mu_{2i} \end{array}\right), \Sigma_{12} \right), \text{ where } \Sigma_{12} = \left(\begin{array}{cc} \tau^{2}_{1} & \rho \tau_1 \tau_2 \\ \rho \tau_1 \tau_2 & \tau^{2}_{2} \end{array}\right)$ where $$\tau^{2}_{1}$$ and $$\tau^{2}_{2}$$ are the between-study variance parameters and $$\rho$$ is the correlation between $$\mu_{1i}$$ and $$\mu_{2i}$$. A similar model is used for the reference test.

To carry out Bayesian estimation, we need to specify prior distribution over the unknown parameters. A hierarchical prior distribution structure is specified for the sensitivity and specificity of the index test to account for both between- and within-study variability and the correlation between the sensitivity and specificity across studies. A similar hierarchical prior distribution is specified over the sensitivity and specificity of the reference standard. Additionally, prior distributions must be specified for the prevalence in each study and for the covariance parameters.

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

# Bayesian bivariate latent class model
modelString =
"model {

#=== LIKELIHOOD ===#
for(i in 1:n) {

t12[i,1:4] ~ dmulti(p12[i,1:4],n.study[i])

p12[i,1] <- prev[i]*( se[i]* se2[i]     +covp[i]) + (1-prev[i])*( (1-sp[i])*(1-sp2[i]) +covn[i])
p12[i,2] <- prev[i]*( se[i]* (1-se2[i]) -covp[i]) + (1-prev[i])*( (1-sp[i])*sp2[i]     -covn[i])
p12[i,3] <- prev[i]*( (1-se[i])*se2[i]     -covp[i]) + (1-prev[i])*( sp[i]*(1-sp2[i]) -covn[i])
p12[i,4] <- prev[i]*( (1-se[i])*(1-se2[i]) +covp[i]) + (1-prev[i])*( sp[i]*sp2[i] +covn[i])

# Hierarchical prior distribution for XPERT test
logit(se[i]) <- l[i,1]
logit(sp[i]) <- l[i,2]
l[i,1:2] ~ dmnorm(mu[1:2], T[1:2,1:2])

# Hierarchical prior distribution for CUTURE test
logit(se2[i]) <- l2[i,1]
logit(sp2[i]) <- l2[i,2]
l2[i,1:2] ~ dmnorm(mu2[1:2],T2[1:2,1:2])

# Prior distribution for prevalence
prev[i] ~ dbeta(1,1)

#=== CONDITIONAL DEPENDENCE STRUCTURE ===#
# lower and upper limits of covariance parameters
maxs12[i]<-min(se[i],se2[i])-(se[i]*se2[i])
maxc12[i]<-min(sp[i],sp2[i])-(sp[i]*sp2[i])
mins12[i]<- (1-se[i])*(se2[i]-1)
minc12[i]<- (1-sp[i])*(sp2[i]-1)

# covariance parameters
covp[i] ~ dunif(mins12[i], maxs12[i]);
covn[i] ~ dunif(minc12[i], maxc12[i]);
}

####    HYPER PRIOR DISTRIBUTIONS FOR XPERT ACCURACY    ####
# prior for the logit transformed sensitivity (mu) and logit transformed  specificity (mu) of XPERT test
mu ~ dnorm(0,0.3)
mu ~ dnorm(0,0.3)

# Between-study variance-covariance matrix of XPERT test
T[1:2,1:2]<-inverse(TAU[1:2,1:2])
TAU[1,1] <- tau*tau
TAU[2,2] <- tau*tau
TAU[1,2] <- rho*tau*tau
TAU[2,1] <- rho*tau*tau

# prior for the between-study precision of the logit transformed sensitivity (prec) and logit transformed specificity (prec) of XPERT test
prec ~ dgamma(2,0.5)
prec ~ dgamma(2,0.5)
rho ~ dunif(-1,1)

####    HYPER PRIOR DISTRIBUTIONS FOR CULTURE ACCURACY    ####
# prior for the logit transformed sensitivity (mu) and logit transformed  specificity (mu) of CULTURE test
mu2 ~ dnorm(0,0.3)
mu2 ~ dnorm(0,0.3)

# Between-study variance-covariance matrix of CULTURE test
T2[1:2,1:2] <- inverse(TAU2[1:2,1:2])
TAU2[1,1] <- tau2*tau2
TAU2[2,2] <- tau2*tau2
TAU2[1,2] <- rho2*tau2*tau2
TAU2[2,1] <- rho2*tau2*tau2

# prior for the between-study precision of the logit transformed sensitivity (prec2) and logit transformed specificity (prec2) of CULTURE test
prec2 ~ dgamma(2,0.5)
prec2 ~ dgamma(2,0.5)
rho2 ~ dunif(-1,1)

####  OTHER PARAMETERS OF INTEREST

# Between-study standard deviation of the logit transformed sensitivity (tau) and specificity (tau) of XPERT test
tau<-pow(prec,-0.5)
tau<-pow(prec,-0.5)
tau.sq<-pow(prec,-1)
tau.sq<-pow(prec,-1)

# Between-study standard deviation of the logit transformed sensitivity (tau2) and specificity (tau2) of CULTURE test
tau2 <- pow(prec2,-0.5)
tau2 <- pow(prec2,-0.5)

# Mean sensitivity and mean specificity of XPERT
Mean_Se <-1/(1+exp(-mu))
Mean_Sp <-1/(1+exp(-mu))

# Predicted sensitivity and predicted specificity of XPERT in a future study
l.predict[1:2] ~ dmnorm(mu[],T[,])
Predicted_Se <- 1/(1+exp(-l.predict))
Predicted_Sp <- 1/(1+exp(-l.predict))

# Mean sensitivity and mean specificity of CULTURE
Mean_Se2<-1/(1+exp(-mu2))
Mean_Sp2<-1/(1+exp(-mu2))

# Predicted sensitivity and predicted specificity of CULTURE in a future study
l.predict2[1:2] ~ dmnorm(mu2[],T2[,])
Predicted.Se2 <- 1/(1+exp(-l.predict2))
Predicted.Sp2 <- 1/(1+exp(-l.predict2))
}"

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

## INTIAL 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. There are different ways of providing initial values. In the illustration below we give details on how to provide initial values chosen by the user by using a homemade function we are calling GenInits().

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) and logit-transformed specificity (mu) of the test under evaluation, as well as the logit-transformed sensitivity (mu2) and logit-transformed specificity (mu2) of the reference test. By setting the chain argument in GenInits() equal to myinits we can provide specific initial values to mu, mu and mu2, mu2 through the init_mu and init_mu2 arguments respectively or we can use the function to randomly generate initial values from distributions of our choice for mu, mu, mu2 and mu2. In this example, we have allowed other parameters, such as between-study precisions (prec and prec for test under evaluation and prec2 and prec2 for reference test), study prevalence of the latent target condition (prev[i], for i=1, …, n) and correlation between logit-transformed sensitivity and logit-transformed specificity (rho for test under evaluation rho2 for reference test), to be randomly generated by GenInits(). In our experience, choosing initial values for these parameters is not easy. Initial values for the covariance between the tests (covp among those who have the target condition and covn among those who do not have the target condition) are initialised to be equal to zero. For reproducibility, a Seed argument is provided.

In the illustration below, for the first chain, we set mu, mu and mu2 equal to 0 (which represents an initial value of 0.5 on the probability scale) while mu2=4, representing a value of 98.2% on the probability scale. For chains number 2, 3, 4 and 5, we simply set the argument chain to 2, 3, 4 and 5, respectively. This ensure the remaining 4 chains are randomly initialized.

# Initial values
GenInits = function(chain, init_mu=c(0,0), init_mu2=c(0,0), Seed=seed, n){
if(chain == "myinits"){
mu = c(as.numeric(init_mu),as.numeric(init_mu))
mu2 = c(as.numeric(init_mu2),as.numeric(init_mu2))
}
else{
mu=c(rnorm(1,0,2),rnorm(1,0,2))
mu2=c(rnorm(1,0,2),rnorm(1,0,2))
}
list(
mu=mu,
mu2=mu2,
covn = rep(0,n),
covp = rep(0,n),
prev = runif(n),
prec = c(rgamma(1,2,0.5),rgamma(1,2,0.5)),
prec2 = c(rgamma(1,2,0.5),rgamma(1,2,0.5)),
rho = runif(1,-1,1),
rho2 = runif(1,-1,1) ,
.RNG.name="base::Wichmann-Hill",
.RNG.seed=Seed
)
}
seed=24050; set.seed(seed)
Listinits = list(
GenInits(chain= "myinits", init_mu=c(0,0), init_mu2=c(0,4), n=n),
GenInits(chain=2, n=n),
GenInits(chain=3, n=n),
GenInits(chain=4, n=n),
GenInits(chain=5, n=n)
)

## COMPILING AND RUNNING THE MODEL

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

library(rjags)
## Loading required package: coda
## Linked to JAGS 4.3.0
## Loaded modules: basemod,bugs
# Compile the model
jagsModel = jags.model("model.txt",data=dataList,n.chains=5,inits=Listinits,n.adapt=0)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 29
##    Unobserved stochastic nodes: 157
##    Total graph size: 1727
##
## Initializing model

The update function serves as 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 5,000 iterations by providing argument n.iter=5000. This can be helpful to mitigate “bad” initial values.

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

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( "prev", "se","sp", "mu", "mu2",  "tau", "tau2", "rho", "rho2", "Mean_Se", "Mean_Se2", "Mean_Sp", "Mean_Sp2", "Predicted_Se", "Predicted_Sp")

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

## 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. prevalence, sensitivity and specificity in each study.

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

for(j in 1:n) {
for(i in c(1,2,3)) {
#  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 = 4 : logit-transformed mean sensitivity of index test (j=1) and logit-transformed mean specificity of index test (j=2) (mu)
# i = 5 : logit-transformed mean sensitivity of reference test (j=1) and logit-transformed mean specificity of reference test (j=2) (mu2)
# i = 6 : 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)
# i = 7 : between-study standard deviation in logit-transformed sensitivity of reference test (j=1) and between-study standard deviation in logit-transformed specificity of reference test (j=2) (tau2)
for(j in 1:2) {
for(i in c(4,5,6)) {
#  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 = 8 : correlation between logit-transformed sensitivity of index test and logit-transformed specificity of index test (rho)
# i = 9 : correlation between logit-transformed sensitivity of reference test and logit-transformed specificity of reference test (rho2)
# i = 10: mean sensitivity of index test (Mean_Se)
# i = 11: mean sensitivity of reference test (Mean_Se2)
# i = 12: mean specificity of index test (Mean_Sp)
# i = 13: mean specificity of reference test (Mean_Sp2)
# i = 14: predicted sensitivity of index test in a future study (Predicted_Se)
# i = 15: predicted specificity of index test in a future study (Predicted_Sp)

for(i in 8:15) {
#  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 (prevalence of the 11th study) are displayed for each chain represented by different color. The very similar results from all five chains suggest the algorithm has converged to the same solution in each case.