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 model (Reitsma and Zwinderman (2005)).
Below we explain each step of a program in the R language (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. 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 studies (Nishimura et al. (2007)). Each study provides a two-by-two table of the form
| Disease_positive | Disease_negative | |
|---|---|---|
| Inex test positive | TP | FP |
| Index test negative | FN | TN |
where:
TPis the number of individuals who tested positive on both testsFPis the number of individuals who tested positive on the index test and negative on the reference testFNis the number of individuals who tested negative on the index test and positive on the reference testTNis the number of individuals who tested negative on both tests
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
Before moving any further, let’s look at a bit of 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:
\[ \begin{pmatrix} \mu_{1i} \\ \mu_{2i} \\ \end{pmatrix} \sim N \begin{pmatrix} \begin{pmatrix} \mu_{1i} \\ \mu_{2i} \\ \end{pmatrix}, \Sigma_{12} \end{pmatrix} \text{with } \Sigma_{12} = \begin{pmatrix} \tau_{1}^2 & \rho \tau_1 \tau_2\\ \rho \tau_1 \tau_2 & \tau_2^2 \\ \end{pmatrix} \] 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}\).
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.
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 (TAU)
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
Summary_Se <- 1/(1+exp(-mu[1]))
Summary_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")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[1]) and logit-transformed specificity (mu[2]) parameters. 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 a normal distribution (corresponding to the prior distribution found in the model) for mu[1] and mu[2]. In this example, we have allowed other parameters, such as between-study precisions (prec[1] and prec[2]) and correlation between logit-transformed sensitivity and logit-transformed specificity (rho) 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 init_mu to NULL. This ensures both chain 2 and chain 3 are randomly initialized.
# Initial values
GenInits = function(init_mu=NULL, Seed){
if(is.null(init_mu) == "FALSE"){
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
)
}
seed=123
set.seed(seed)
Listinits = list(
GenInits(init_mu=c(0,0), Seed=seed),
GenInits(Seed=seed),
GenInits(Seed=seed)
)COMPILING AND RUNNING THE MODEL
The jags.model function compiles the model 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)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 74
Unobserved stochastic nodes: 43
Total graph size: 375
Initializing modelThe 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 30,000 iterations by providing argument n.iter=30000. This can be helpful to mitigate “bad” initial values.
# Burn-in iterations
update(jagsModel,n.iter=30000)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=30000 iterations from the posterior distributions for each chain.
# Parameters to be monitored
parameters = c( "se","sp", "mu", "tau.sq", "tau", "rho", "Summary_Se", "Summary_Sp", "Predicted_Se", "Predicted_Sp")
# Posterior samples
output = coda.samples(jagsModel,variable.names=parameters,n.iter=30000)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 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 Summary_Se (summary sensitivity of the Anti-CCP test) 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.
SUMMARY STATISTICS
Once convergence of the MCMC algorithm is achieved, summary statistics of the parameters of interest may be extracted from their posterior distributions using the MCMCsummary function found in the MCMCvis package.
# Summary statistics
library(MCMCvis)
MCMCsummary(output, params=c("se","sp", "mu", "tau.sq", "tau", "rho", "Summary_Se", "Summary_Sp", "Predicted_Se", "Predicted_Sp" ), round=4)Below is the summary statistics output for the parameters shared across studies. The 50% quantile is the posterior median which is commonly used to estimate the paramater of interest. 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 in the index test.
The table also displays Rhat, which is the Gelman-Rubin statistic (Gelman and Rubin (1992), Brooks and Gelman (1998)). It is enabled when 2 or more chains are generated. It evaluates MCMC convergence by comparing within- and between-chain variability for each model parameter. Rhat tends to 1 as convergence is approached.
Also, n.eff is the effective sample size (Gelman et al. (2013)). Because the MCMC process causes the posterior draws to be correlated, the effective sample size is an estimate of the sample size required to achieve the same level of precision if that sample was a simple random sample. When draws are correlated, the effective sample size will generally be lower than the actual numbers of draws.
# Summary statistics
library(MCMCvis)
MCMCsummary(output, params=c("se","sp", "mu", "tau.sq", "tau", "rho", "Summary_Se", "Summary_Sp", "Predicted_Se", "Predicted_Sp" ), round=4) mean sd 2.5% 50% 97.5% Rhat n.eff
se[1] 0.8663 0.0284 0.8063 0.8681 0.9166 1 9908
se[2] 0.5702 0.0346 0.5015 0.5705 0.6370 1 12240
se[3] 0.4252 0.0471 0.3339 0.4246 0.5178 1 9935
se[4] 0.7251 0.0699 0.5800 0.7288 0.8500 1 9196
se[5] 0.7275 0.0242 0.6795 0.7279 0.7737 1 13112
se[6] 0.8716 0.0335 0.7976 0.8741 0.9292 1 10292
se[7] 0.7432 0.0384 0.6648 0.7447 0.8156 1 10609
se[8] 0.6396 0.0390 0.5626 0.6402 0.7140 1 10917
se[9] 0.5843 0.0617 0.4611 0.5856 0.7014 1 10471
se[10] 0.7874 0.0407 0.7032 0.7895 0.8609 1 10715
se[11] 0.7011 0.0672 0.5599 0.7051 0.8222 1 10998
se[12] 0.4195 0.0458 0.3307 0.4188 0.5104 1 10828
se[13] 0.7814 0.0415 0.6956 0.7837 0.8569 1 10072
se[14] 0.6328 0.0287 0.5752 0.6332 0.6880 1 11970
se[15] 0.6695 0.0671 0.5307 0.6718 0.7937 1 9918
se[16] 0.4324 0.0303 0.3738 0.4324 0.4920 1 12529
se[17] 0.5736 0.0673 0.4371 0.5750 0.7022 1 10034
se[18] 0.7897 0.0427 0.7003 0.7919 0.8674 1 10940
se[19] 0.5556 0.0424 0.4732 0.5558 0.6392 1 11345
se[20] 0.6675 0.0439 0.5770 0.6689 0.7501 1 11510
se[21] 0.9017 0.0404 0.8092 0.9072 0.9655 1 7980
se[22] 0.4278 0.0474 0.3358 0.4279 0.5234 1 10531
se[23] 0.5811 0.0300 0.5224 0.5811 0.6388 1 12687
se[24] 0.8014 0.0284 0.7431 0.8025 0.8540 1 11379
se[25] 0.6850 0.0529 0.5761 0.6861 0.7843 1 10226
se[26] 0.5835 0.0685 0.4445 0.5854 0.7135 1 10103
se[27] 0.4889 0.0509 0.3912 0.4886 0.5888 1 11162
se[28] 0.7370 0.0283 0.6796 0.7379 0.7903 1 11741
se[29] 0.4946 0.0396 0.4171 0.4942 0.5733 1 11263
se[30] 0.5113 0.0995 0.3158 0.5130 0.6998 1 10164
se[31] 0.8723 0.0140 0.8435 0.8726 0.8984 1 13353
se[32] 0.6432 0.0274 0.5889 0.6437 0.6958 1 13699
se[33] 0.5436 0.0388 0.4670 0.5440 0.6182 1 10696
se[34] 0.7730 0.0124 0.7478 0.7730 0.7968 1 12724
se[35] 0.5786 0.0311 0.5180 0.5786 0.6388 1 12113
se[36] 0.4044 0.0358 0.3339 0.4038 0.4751 1 10536
se[37] 0.4777 0.0355 0.4079 0.4776 0.5475 1 10758
sp[1] 0.8357 0.0367 0.7562 0.8382 0.8996 1 8768
sp[2] 0.9094 0.0174 0.8723 0.9104 0.9410 1 9207
sp[3] 0.9758 0.0086 0.9561 0.9769 0.9893 1 9148
sp[4] 0.9665 0.0181 0.9216 0.9699 0.9913 1 8822
sp[5] 0.9244 0.0156 0.8908 0.9255 0.9521 1 8540
sp[6] 0.9218 0.0205 0.8764 0.9239 0.9563 1 10709
sp[7] 0.9650 0.0127 0.9363 0.9666 0.9849 1 9104
sp[8] 0.9740 0.0110 0.9481 0.9757 0.9908 1 8953
sp[9] 0.9779 0.0116 0.9493 0.9801 0.9940 1 9365
sp[10] 0.8698 0.0424 0.7739 0.8741 0.9387 1 9435
sp[11] 0.9495 0.0236 0.8923 0.9533 0.9832 1 10478
sp[12] 0.9809 0.0092 0.9589 0.9826 0.9942 1 8723
sp[13] 0.9705 0.0098 0.9488 0.9716 0.9865 1 9111
sp[14] 0.9294 0.0222 0.8797 0.9317 0.9653 1 8752
sp[15] 0.8102 0.0692 0.6528 0.8189 0.9204 1 7711
sp[16] 0.9722 0.0117 0.9444 0.9741 0.9894 1 8360
sp[17] 0.9688 0.0155 0.9304 0.9715 0.9906 1 8562
sp[18] 0.9527 0.0113 0.9280 0.9536 0.9720 1 11506
sp[19] 0.9652 0.0158 0.9277 0.9678 0.9880 1 9558
sp[20] 0.9148 0.0211 0.8690 0.9164 0.9510 1 10378
sp[21] 0.9444 0.0222 0.8923 0.9473 0.9792 1 7957
sp[22] 0.9738 0.0120 0.9447 0.9757 0.9910 1 8417
sp[23] 0.9475 0.0172 0.9088 0.9496 0.9753 1 8989
sp[24] 0.9198 0.0225 0.8699 0.9218 0.9573 1 10243
sp[25] 0.9767 0.0070 0.9612 0.9774 0.9884 1 9889
sp[26] 0.9610 0.0163 0.9220 0.9633 0.9853 1 9813
sp[27] 0.9392 0.0171 0.9006 0.9408 0.9676 1 10156
sp[28] 0.9449 0.0101 0.9237 0.9456 0.9631 1 10063
sp[29] 0.9570 0.0107 0.9335 0.9579 0.9752 1 9414
sp[30] 0.9619 0.0186 0.9155 0.9651 0.9880 1 10264
sp[31] 0.8940 0.0200 0.8518 0.8953 0.9296 1 9049
sp[32] 0.9690 0.0077 0.9517 0.9697 0.9820 1 7375
sp[33] 0.9594 0.0101 0.9374 0.9602 0.9767 1 9584
sp[34] 0.9651 0.0038 0.9575 0.9652 0.9722 1 8876
sp[35] 0.9812 0.0057 0.9685 0.9816 0.9907 1 7950
sp[36] 0.9659 0.0127 0.9369 0.9677 0.9856 1 9346
sp[37] 0.9755 0.0077 0.9582 0.9765 0.9880 1 8950
mu[1] 0.6557 0.1283 0.4045 0.6554 0.9109 1 34599
mu[2] 3.0858 0.1433 2.8088 3.0836 3.3760 1 16020
tau.sq[1] 0.5525 0.1496 0.3277 0.5296 0.9036 1 17090
tau.sq[2] 0.5637 0.1844 0.2891 0.5356 1.0055 1 7149
tau[1] 0.7369 0.0972 0.5724 0.7278 0.9506 1 15999
tau[2] 0.7414 0.1186 0.5377 0.7319 1.0028 1 6737
rho -0.4212 0.1595 -0.6985 -0.4329 -0.0776 1 9856
Summary_Se 0.6577 0.0288 0.5998 0.6582 0.7132 1 34685
Summary_Sp 0.9559 0.0060 0.9431 0.9562 0.9669 1 16270
Predicted_Se 0.6419 0.1561 0.3014 0.6584 0.8958 1 89201
Predicted_Sp 0.9444 0.0440 0.8277 0.9562 0.9900 1 87044SROC 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. The argument ref_std=TRUEis required to allow the plot to reflect the assumption of a perfect reference test.
# 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
Below is a copy of the complete script to run the Bayesian Bivariate meta-analysis model in rjags.
GenInits = function(init_mu=NULL, Seed=seed){
if(is.null(init_mu) == “FALSE”){
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
)
}
seed=123
Listinits = list(
GenInits(init_mu=c(0,0)),
GenInits(),
GenInits()
)
####################################
#### DATA FOR MOTIVATING EXAMPLE ###
####################################
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 ###
#################
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 (TAU)
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
Summary_Se <- 1/(1+exp(-mu[1]))
Summary_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 ###
######################
GenInits = function(init_mu=NULL, Seed){
if(is.null(init_mu) == "FALSE"){
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
)
}
seed=123
set.seed(seed)
Listinits = list(
GenInits(init_mu=c(0,0), Seed=seed),
GenInits(Seed=seed),
GenInits(Seed=seed)
)
#######################################
### COMPILING AND RUNNING THE MODEL ###
#######################################
library(rjags)
jagsModel = jags.model("model.txt",data=dataList,n.chains=3,inits=Listinits)
# Burn-in iterations
update(jagsModel,n.iter=30000)
# Parameters to be monitored
parameters = c( "se","sp", "mu", "tau.sq", "tau", "rho", "Summary_Se", "Summary_Sp", "Predicted_Se", "Predicted_Sp")
# Posterior samples
output = coda.samples(jagsModel,variable.names=parameters,n.iter=30000)
##############################
### MONITORING CONVERGENCE ###
##############################
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()
}
}
# 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()
}
##########################
### SUMMARY STATISTICS ###
##########################
library(MCMCvis)
library(DT)
datatable(MCMCsummary(output, round=4), options = list(pageLength = 3) )
#################
### 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))REFERENCES
Citation
@online{schiller2021,
author = {Ian Schiller and Nandini Dendukuri},
title = {Bayesian Inference for the {Bivariate} Model for Diagnostic
Meta-Analysis},
date = {2021-04-13},
url = {https://www.nandinidendukuri.com/blogposts/2021-04-13-bayesian-bivariate-meta-analysis/},
langid = {en}
}