---
title: "Bayesian Bivariate Latent CLass Meta-analysis Model"
author: "Ian Schiller and Nandini Dendukuri"
date: "Last compiled on `r format(Sys.time(), '%d %B, %Y, %Hh%M')`"
output:
  html_document:
    always_allow_html: true
    code_folding: hide
    code_download: true
    theme: united
    toc: yes
    toc_float: yes
    collapsed: false
    smooth_scroll: true
    toc_depth: 3
    number_sections: no
  pdf_document:
    toc: yes
  word_document:
    toc: yes
---

```{r, eval=T, echo=F,  results='hide', message=FALSE, warning=FALSE }
name_of_data = "Saliva_COVID" # Provide the name of your dataset excluding the ".txt" extension
digits=4                      # To control the number of digits to show in the output results
num.burn=100                  # Number of iterations to drop during the "burn-in" phase
num.iter=100                  # Number of iterations to keep to form the posterior samples
num.chains=3                  # Number of independent chains to run
Cond.Dep=FALSE                # To model condition dependence, set "Cond.Dep" equals to "TRUE"

################################  PROVIDING INITIAL VALUES --- OPTIONAL ##################################
###                                                                                                    ###
# If your.inits below is set to FALSE, the rest of the code in this block will be ignored and jags will###
# generate initial values while the user does not need to change anything further below.  If your.inits###
# is set to TRUE, then initial values can be provided for the following 4 parameters :                 ###
#                                                                                                      ###
# Mean logit-transformed sensitivity of index test (referred below as mu_se)                           ###
# Mean logit-transformed specificity of index test (referred below as mu_sp)                           ###
# Mean logit-transformed sensitivity of reference test (referred below as mu2_se)                      ###
# Mean logit-transformed specificity of reference test (referred below as mu2_sp)                      ###
#                                                                                                      ###
your.inits=TRUE                                                                                        
#                                                                                                      ###
# If your.inits was set to TRUE, then you will be required to provide initial values in the form of    ###
# a vector of length num.chain for each of the 4 parameters described above.                           ###
#                                                                                                      ###
mu_se = c(0,0.5,1)  # A VALUE OF ZERO TRANSPOSES INTO A VALUE OF 0.5 ON THE PROBABILITY SCALE
mu_sp = c(1,2,3)  # A VALUE OF ZERO TRANSPOSES INTO A VALUE OF 0.5 ON THE PROBABILITY SCALE
mu2_se = c(0,0.5,1) # A VALUE OF ZERO TRANSPOSES INTO A VALUE OF 0.5 ON THE PROBABILITY SCALE
mu2_sp = c(3,3.5,4) # A VALUE OF ZERO TRANSPOSES INTO A VALUE OF 0.5 ON THE PROBABILITY SCALE
#                                                                                                      ###
# For reproducibility, a seed value can be provided                                                    ###
seed = sample(0:(10^3), 1)  #259
###                                                                                                    ###
##########################################################################################################






```

```{r, eval=T, echo=F, results='hide', message=FALSE, warning=FALSE    }
# install.packages("remotes")
#remotes::install_github("rlesur/klippy")
require(klippy)
require(knitr)
require(plyr)
require(MCMCvis)
require(mcmcplots)
require(rjags)
require(DTAplots)
require(DT)  # THIS LIBRARY ALLOWS THE NICE DATA DISPLAY WITH THE SEARCH BAR OPTION.
require(mvtnorm)  # FOR GENERATION OF MULTIVARIATE NORMAL NUMBERS IN Geninits FUNCTION
```

```{r klippy, echo=FALSE, include=TRUE}
klippy::klippy()
```

```{r setup, include=FALSE, echo=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning=FALSE, fig.width=12, fig.height=8)
```

```{r, eval=T, echo=F,  results='hide', message=FALSE, warning=FALSE }
options(max.print=100000)
DATA <- read.table(paste(name_of_data,".txt", sep=""), header=TRUE)

if(c("tp") %in% colnames(DATA) | c("fp") %in% colnames(DATA) | c("fn") %in% colnames(DATA) | c("tn") %in% colnames(DATA)){
  DATA = rename(DATA, c("tp" = "n11", "fp" = "n10", "fn" = "n01", "tn" = "n00"))
}

N = dim(DATA)[1] # Number of studies
n11 = DATA$n11 # Both tests positive
n10 = DATA$n10 # Index test positive, reference test negative
n01 = DATA$n01 # Index test negative, reference test positive
n00 = DATA$n00 # Both tests negative
cell<-cbind(n11, n10, n01, n00)
n <- rowSums(cell) #  Study sample sizes
dataList = list(cell=cell,N=N, n=n)

```

## **`r paste(name_of_data)` DATASET**

The `r paste(name_of_data)` database includes `r sum(n)` participants spread across `r dim(DATA)[1]` studies as follow. `r if(name_of_data=="Saliva_COVID"){"The Saliva_COVID data is from Butler-Laporte et al(2021)"}`

```{r, eval=T, echo=F, class.source = 'fold-show', message=FALSE }
datatable(DATA, extensions = 'AutoFill')#, options = list(autoFill = TRUE))
```

```{r, eval=T, echo=F, class.source = 'fold-show', message=FALSE }
# Initial values


GenInits  = function(Seed=seed, n=N, mu_se_index=0, mu_sp_index=0, mu_se_ref=0, mu_sp_ref=0, cond_dep=Cond.Dep){
  
  ### SENSITIVITY OF INDEX TEST
  mu.1 = as.numeric(mu_se_index[1])
  
  ### SPECIFICITY OF INDEX TEST
  mu.2 = as.numeric(mu_sp_index[1])
  
  ### SENSITIVITY OF REFERENCE TEST
  mu2.1 = as.numeric(mu_se_ref[1])
  
  ### SPECIFICITY OF REFERENCE TEST
  mu2.2 = as.numeric(mu_sp_ref[1])
  
  
  mu = c(mu.1,mu.2)
  mu2 = c(mu2.1,mu2.2)
  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);
  Predicted.l = c( rmvnorm(1, mean = mu, sigma = matrix(ncol=2,c(1/prec[1], rho*sqrt(1/prec[1])*sqrt(1/prec[2]), rho*sqrt(1/prec[1])*sqrt(1/prec[2]), 1/prec[2]))) );
  Predicted.l2 = c( rmvnorm(1, mean = mu2, sigma = matrix(ncol=2,c(1/prec2[1], rho2*sqrt(1/prec2[1])*sqrt(1/prec2[2]), rho2*sqrt(1/prec2[1])*sqrt(1/prec2[2]), 1/prec2[2]))) );
  l = rmvnorm(n, mean = mu, sigma = matrix(ncol=2,c(1/prec[1], rho*sqrt(1/prec[1])*sqrt(1/prec[2]), rho*sqrt(1/prec[1])*sqrt(1/prec[2]), 1/prec[2])));
  l2 = rmvnorm(n, mean = mu2, sigma = matrix(ncol=2,c(1/prec2[1], rho2*sqrt(1/prec2[1])*sqrt(1/prec2[2]), rho2*sqrt(1/prec2[1])*sqrt(1/prec2[2]), 1/prec2[2])));  
  prev = rbeta(n,1,1)          
  
  
  
  if(cond_dep == TRUE){
    
    #s.c = 1/(1+exp(-l))
    #se = s.c[,1]
    #sp = s.c[,2]
    #s2.c2 = 1/(1+exp(-l2))
    #se2 = s2.c2[,1]
    #sp2 = s2.c2[,2]
    #
    #covp = runif(n,(-(1-se)*(1-se2)), (max(se,se2)-se*se2))
    #covn = runif(n,(-(1-sp)*(1-sp2)), (max(sp,sp2)-sp*sp2))
    
    covp = covn = rep(0,n)
    
    list(
      covp = covp,
      covn = covn,
      mu = mu,
      mu2 = mu2,
      prec = prec,
      prec2 = prec2,
      prev = prev,
      rho = rho,
      rho2 = rho2,
      Predicted.l = Predicted.l,
      Predicted.l2 = Predicted.l2,
      l = l,
      l2 = l2,
      .RNG.name="base::Wichmann-Hill",
      .RNG.seed=Seed
    )
  }
  else{
    list(
      mu = mu,
      mu2 = mu2,
      prec = prec,
      prec2 = prec2,
      prev = prev,
      rho = rho,
      rho2 = rho2,
      Predicted.l = Predicted.l,
      Predicted.l2 = Predicted.l2,
      l = l,
      l2 = l2,
      .RNG.name="base::Wichmann-Hill",
      .RNG.seed=Seed
    ) 
  }
}


Listinits = vector('list',num.chains)
for(i in 1:num.chains){
    tmp = (GenInits(Seed=seed, n=N, mu_se_index=mu_se[i], mu_sp_index=mu_sp[i], mu_se_ref=mu2_se[i], mu_sp_ref=mu2_sp[i], cond_dep=Cond.Dep))
    Listinits[[i]] = tmp
}

```

`r if(Cond.Dep==TRUE){"## **BAYESIAN BIVARIATE LATENT CLASS META-ANALYSIS MODEL ASSUMING CONDITIONNAL DEPENDENCE**"}else{"## **BAYESIAN BIVARIATE LATENT CLASS META-ANALYSIS MODEL ASSUMING CONDITIONNAL INDEPENDENCE**"}`

```{r, eval=T, echo=T, message=FALSE }
# Bayesian bivariate model
if(Cond.Dep == TRUE){
  
modelString =

"model {
#=================================================================
############ LIKELIHOOD	############
#=================================================================
for(i in 1:N) {	
# The 4 cells in the two-by-two table of the new test vs. the imperfect reference standard in each study in the meta-analysis	
cell[i,1:4] ~ dmulti(prob[i,1:4],n[i])
# Multinomial probabilities of the 4 cells of the two-by-two table expressed in terms of the disease prevalence (prev), the sensitivity (se1, se2) and specificity (sp1, sp2) of the two tests and the covariance between them among disease positive (covs12) and disease negative (covc12) groups
se[i] <- p[1,i]
sp[i] <- 1-p[2,i]		
prob[i,1] <- prev[i]*( p[1,i]     * se2[i]     + covp[i] ) + (1-prev[i])*( p[2,i]     *(1-sp2[i]) + covn[i] )
prob[i,2] <- prev[i]*( p[1,i]     * (1-se2[i]) - covp[i] ) + (1-prev[i])*( p[2,i]     *sp2[i]     - covn[i] )
prob[i,3] <- prev[i]*( (1-p[1,i]) * se2[i]     - covp[i] ) + (1-prev[i])*( (1-p[2,i]) *(1-sp2[i]) - covn[i] )
prob[i,4] <- prev[i]*( (1-p[1,i]) * (1-se2[i]) + covp[i] ) + (1-prev[i])*( (1-p[2,i]) *sp2[i]     + covn[i] )


#=================================================================
############ CONDITIONAL DEPENDENCE	############
#=================================================================
#=======================================
# upper limits of covariance parameters
#=======================================
us[i]<-min(se[i],se2[i])-(se[i]*se2[i]);
uc[i]<-min(sp[i],sp2[i])-(sp[i]*sp2[i]);
ls[i]<- -(1-se[i])*(1-se2[i])
lc[i]<- -(1-sp[i])*(1-sp2[i])
#==============================================================
# prior distribution of transformed covariances on (0,1) range
#==============================================================
covp[i]~dunif(ls[i],us[i]);
covn[i]~dunif(lc[i],uc[i]);


#=================================================================
############ HIERARCHICAL PRIORS ############	
#=================================================================
#=======================================
# Hierarchichal prior for xpert
#=======================================
logit(p[1, i]) <- l[i,1]
logit(p[2, i]) <- -l[i,2]
l[i,1:2] ~ dmnorm(mu[], T[,])
#=======================================
# Prior distribution on prevalence
#=======================================
prev[i] ~ dbeta(1,1)
}
# ==================================	
# Hierarchical prior for CULTURE
# ==================================	
for(j in 1:N) {
logit(se2[j]) <- l2[j,1]
logit(sp2[j]) <- l2[j,2]
l2[j,1:2] ~ dmnorm(mu2[1:2], T2[1:2,1:2])			
}


#=================================================================
############  HYPER PRIOR DISTRIBUTIONS    ################
#=================================================================
### 			INDEX TEST
#=================================================================
### prior for the logit transformed sensitivity (mu[1]) and specificity (mu[2])
mu[1] ~ dnorm(0,0.25)
mu[2] ~ dnorm(0,0.25) 
T[1:2,1:2]<-inverse(TAU[1:2,1:2])
#### BETWEEN-STUDY VARIANCE-COVARIANCE MATRIX
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]	
#### prec = between-study precision in the logit(sensitivity) and logit(specificity)
prec[1] ~ dgamma(2,0.5)
prec[2] ~ dgamma(2,0.5)
rho ~ dunif(-1,1)
#=================================================================
### 			REFERENCE TEST
#=================================================================
### prior for the logit transformed sensitivity (mu2[1]) and specificity (mu2[2])
mu2[1] ~ dnorm(0,0.25)
mu2[2] ~ dnorm(0,0.25)
T2[1:2,1:2]<-inverse(TAU2[1:2,1:2])
#### BETWEEN-STUDY VARIANCE-COVARIANCE MATRIX
TAU2[1,1] <- tau2[1]*tau2[1]
TAU2[2,2] <- tau2[2]*tau2[2]
TAU2[1,2] <- rho2*tau2[1]*tau2[2]	
TAU2[2,1] <- rho2*tau2[1]*tau2[2]
#### prec = between-study precision in the logit(sensitivity) and logit(specificity)
prec2[1] ~ dgamma(2,0.5)
prec2[2] ~ dgamma(2,0.5)
rho2 ~ dunif(-1,1)


#=================================================================
############ OTHER PARAMETERS OF INTEREST ############
#=================================================================
### 			INDEX TEST
#=================================================================
#### SUMMARY SENSITIVITY AND SPECIFICITY OF INDEX TEST
Summary_Se<-1/(1+exp(-mu[1]))
Summary_Sp<-1/(1+exp(-mu[2]))
#### BETWEEN_STUDY VARIANCE IN THE LOGIT(SENSITIVITY) AND LOGIT(SPECIFICITY)
tau.sq[1] <- pow(tau[1], 2)
tau.sq[2] <- pow(tau[2], 2)
#### BETWEEN_STUDY STANDARD DEVIATION IN THE LOGIT(SENSITIVITY) AND LOGIT(SPECIFICITY)
tau[1]<-pow(prec[1],-0.5)
tau[2]<-pow(prec[2],-0.5)
#### PREDICTED SENSITIVITY AND SPECIFICITY OF XPERT IN A FUTURE STUDY
Predicted.l[1:2] ~ dmnorm(mu[],T[,])
Predicted_Se <- 1/(1+exp(-Predicted.l[1]))
Predicted_Sp <- 1/(1+exp(-Predicted.l[2]))
#=================================================================
### 			REFERENCE TEST
#=================================================================
#### SUMMARY SENSITIVITY AND SPECIFICITY OF REFERENCE TEST
Summary_Se2<-1/(1+exp(-mu2[1]))
Summary_Sp2<-1/(1+exp(-mu2[2]))
#### BETWEEN_STUDY STANDARD DEVIATION IN THE LOGIT(SENSITIVITY) AND LOGIT(SPECIFICITY)
tau2[1] <-pow(prec2[1],-0.5)
tau2[2] <-pow(prec2[2],-0.5)
#### BETWEEN_STUDY VARIANCE IN THE LOGIT(SENSITIVITY) AND LOGIT(SPECIFICITY)
tau.sq2[1] <- pow(tau2[1], 2)
tau.sq2[2] <- pow(tau2[2], 2)
#### PREDICTED SENSITIVITY AND SPECIFICITY OF CULTURE IN A FUTURE STUDY
Predicted.l2[1:2] ~ dmnorm(mu2[],T2[,])
Predicted_Se2 <- 1/(1+exp(-Predicted.l2[1]))
Predicted_Sp2 <- 1/(1+exp(-Predicted.l2[2]))	
}"

}else{
  
  modelString =

"model {
#=================================================================
############ LIKELIHOOD	############
#=================================================================
for(i in 1:N) {	
# The 4 cells in the two-by-two table of the new test vs. the imperfect reference standard in each study in the meta-analysis	
cell[i,1:4] ~ dmulti(prob[i,1:4],n[i])
# Multinomial probabilities of the 4 cells of the two-by-two table expressed in terms of the disease prevalence (prev), the sensitivity (se1, se2) and specificity (sp1, sp2) of the two tests and the covariance between them among disease positive (covs12) and disease negative (covc12) groups
se[i] <- p[1,i]
sp[i] <- 1-p[2,i]		
prob[i,1] <- prev[i]*( p[1,i]     * se2[i]     + covp[i] ) + (1-prev[i])*( p[2,i]     *(1-sp2[i]) + covn[i] )
prob[i,2] <- prev[i]*( p[1,i]     * (1-se2[i]) - covp[i] ) + (1-prev[i])*( p[2,i]     *sp2[i]     - covn[i] )
prob[i,3] <- prev[i]*( (1-p[1,i]) * se2[i]     - covp[i] ) + (1-prev[i])*( (1-p[2,i]) *(1-sp2[i]) - covn[i] )
prob[i,4] <- prev[i]*( (1-p[1,i]) * (1-se2[i]) + covp[i] ) + (1-prev[i])*( (1-p[2,i]) *sp2[i]     + covn[i] )


#=================================================================
############ CONDITIONAL DEPENDENCE	############
#=================================================================
#=======================================
# upper limits of covariance parameters
#=======================================
us[i]<-min(se[i],se2[i])-(se[i]*se2[i]);
uc[i]<-min(sp[i],sp2[i])-(sp[i]*sp2[i]);
ls[i]<- -(1-se[i])*(1-se2[i])
lc[i]<- -(1-sp[i])*(1-sp2[i])
#==============================================================
# prior distribution of transformed covariances on (0,1) range
#==============================================================
covp[i]<-0 #~dunif(ls[i],us[i]);
covn[i]<-0 #~dunif(lc[i],uc[i]);


#=================================================================
############ HIERARCHICAL PRIORS ############	
#=================================================================
#=======================================
# Hierarchichal prior for xpert
#=======================================
logit(p[1, i]) <- l[i,1]
logit(p[2, i]) <- -l[i,2]
l[i,1:2] ~ dmnorm(mu[], T[,])
#=======================================
# Prior distribution on prevalence
#=======================================
prev[i] ~ dbeta(1,1)
}
# ==================================	
# Hierarchical prior for CULTURE
# ==================================	
for(j in 1:N) {
logit(se2[j]) <- l2[j,1]
logit(sp2[j]) <- l2[j,2]
l2[j,1:2] ~ dmnorm(mu2[1:2], T2[1:2,1:2])			
}


#=================================================================
############  HYPER PRIOR DISTRIBUTIONS    ################
#=================================================================
### 			INDEX TEST
#=================================================================
### prior for the logit transformed sensitivity (mu[1]) and specificity (mu[2])
mu[1] ~ dnorm(0,0.25)
mu[2] ~ dnorm(0,0.25) 
T[1:2,1:2]<-inverse(TAU[1:2,1:2])
#### BETWEEN-STUDY VARIANCE-COVARIANCE MATRIX
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]	
#### prec = between-study precision in the logit(sensitivity) and logit(specificity)
prec[1] ~ dgamma(2,0.5)
prec[2] ~ dgamma(2,0.5)
rho ~ dunif(-1,1)
#=================================================================
### 			REFERENCE TEST
#=================================================================
### prior for the logit transformed sensitivity (mu2[1]) and specificity (mu2[2])
mu2[1] ~ dnorm(0,0.25)
mu2[2] ~ dnorm(0,0.25)
T2[1:2,1:2]<-inverse(TAU2[1:2,1:2])
#### BETWEEN-STUDY VARIANCE-COVARIANCE MATRIX
TAU2[1,1] <- tau2[1]*tau2[1]
TAU2[2,2] <- tau2[2]*tau2[2]
TAU2[1,2] <- rho2*tau2[1]*tau2[2]	
TAU2[2,1] <- rho2*tau2[1]*tau2[2]
#### prec = between-study precision in the logit(sensitivity) and logit(specificity)
prec2[1] ~ dgamma(2,0.5)
prec2[2] ~ dgamma(2,0.5)
rho2 ~ dunif(-1,1)


#=================================================================
############ OTHER PARAMETERS OF INTEREST ############
#=================================================================
### 			INDEX TEST
#=================================================================
#### SUMMARY SENSITIVITY AND SPECIFICITY OF INDEX TEST
Summary_Se<-1/(1+exp(-mu[1]))
Summary_Sp<-1/(1+exp(-mu[2]))
#### BETWEEN_STUDY VARIANCE IN THE LOGIT(SENSITIVITY) AND LOGIT(SPECIFICITY)
tau.sq[1] <- pow(tau[1], 2)
tau.sq[2] <- pow(tau[2], 2)
#### BETWEEN_STUDY STANDARD DEVIATION IN THE LOGIT(SENSITIVITY) AND LOGIT(SPECIFICITY)
tau[1]<-pow(prec[1],-0.5)
tau[2]<-pow(prec[2],-0.5)
#### PREDICTED SENSITIVITY AND SPECIFICITY OF XPERT IN A FUTURE STUDY
Predicted.l[1:2] ~ dmnorm(mu[],T[,])
Predicted_Se <- 1/(1+exp(-Predicted.l[1]))
Predicted_Sp <- 1/(1+exp(-Predicted.l[2]))
#=================================================================
### 			REFERENCE TEST
#=================================================================
#### SUMMARY SENSITIVITY AND SPECIFICITY OF REFERENCE TEST
Summary_Se2<-1/(1+exp(-mu2[1]))
Summary_Sp2<-1/(1+exp(-mu2[2]))
#### BETWEEN_STUDY STANDARD DEVIATION IN THE LOGIT(SENSITIVITY) AND LOGIT(SPECIFICITY)
tau2[1] <-pow(prec2[1],-0.5)
tau2[2] <-pow(prec2[2],-0.5)
#### BETWEEN_STUDY VARIANCE IN THE LOGIT(SENSITIVITY) AND LOGIT(SPECIFICITY)
tau.sq2[1] <- pow(tau2[1], 2)
tau.sq2[2] <- pow(tau2[2], 2)
#### PREDICTED SENSITIVITY AND SPECIFICITY OF CULTURE IN A FUTURE STUDY
Predicted.l2[1:2] ~ dmnorm(mu2[],T2[,])
Predicted_Se2 <- 1/(1+exp(-Predicted.l2[1]))
Predicted_Sp2 <- 1/(1+exp(-Predicted.l2[2]))	
}"
  
}
```

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE, results='hide' }
writeLines(modelString,con="model.txt")

  # Compile the model
if(your.inits==FALSE){
  jagsModel = jags.model("model.txt",data=dataList,n.chains=num.chains,n.adapt=0)
}else{
  jagsModel = jags.model("model.txt",data=dataList,n.chains=num.chains,n.adapt=0, inits=Listinits)
}

#jagsModel$state(internal=FALSE)

# Burn-in iterations 
update(jagsModel,n.iter=num.burn)
# Parameters to be monitored
if(Cond.Dep == TRUE){
  
parameters = c( "se","sp", "mu", "tau.sq", "tau","rho", "Summary_Se" , "Summary_Sp",   "Predicted_Se", "Predicted_Sp", "prec",
                "se2","sp2", "mu2", "tau.sq2", "tau2","rho2", "Summary_Se2" , "Summary_Sp2",   "Predicted_Se2", "Predicted_Sp2", "prec2", 
                "prev", "covp", "covn")
}else{
  parameters = c( "se","sp", "mu", "tau.sq", "tau","rho", "Summary_Se" , "Summary_Sp",   "Predicted_Se", "Predicted_Sp", "prec",
                "se2","sp2", "mu2", "tau.sq2", "tau2","rho2", "Summary_Se2" , "Summary_Sp2",   "Predicted_Se2", "Predicted_Sp2", "prec2", 
                "prev")
}

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

# To display if conditional dependence was selected
covp_print <- '$covp \\sim U \\left(-(1-s_1)\\times(1-s_2), max(s_1,s_2)-s_1\\times s_2 \\right)$'
covn_print <- '$covn \\sim U \\left(-(1-c_1)\\times(1-c_2), max(c_1,c_2)-c_1\\times c_2 \\right)$'
```

The posterior samples were obtained by running a latent-class meta-analysis model (Xie, Sinclair and Dendukuri(2017)) with `r num.chains` independent chains each having their own starting values to assess convergence of the MCMC algorithm. For each chain, after discarding the first `r num.burn` iterations (burn-in), we kept the next `r options(scipen = 1000000) ; num.iter` iterations to form a posterior samples of `r options(scipen = 1000000) ; 3*num.iter` iterations.

\*\* Prior Distributions\*\*

A hierarchical prior distribution structure is specified for the sensitivity and specificity of the index test. The prior distributions are as followed

-   $mu[1] \sim N(\mu=0,\sigma^2=4)$
-   $mu[2] \sim N(\mu=0,\sigma^2=4)$
-   $prec[1] \sim \Gamma(\alpha=2,s=\frac{1}{2})$ where $\alpha$ and $s$ are the shape and scale parameters respectively
-   $prec[2] \sim \Gamma(\alpha=2,s=\frac{1}{2})$ where $\alpha$ and $s$ are the shape and scale parameters respectively
-   $rho \sim U(-1,1)$

A similar hierarchical prior distribution is specified over the sensitivity and specificity of the reference standard. The prior distributions are as followed

-   $mu2[1] \sim N(\mu=0,\sigma^2=4)$
-   $mu2[2] \sim N(\mu=0,\sigma^2=4)$
-   $prec2[1] \sim \Gamma(\alpha=2,s=\frac{1}{2})$ where $\alpha$ and $s$ are the shape and scale parameters respectively
-   $prec2[2] \sim \Gamma(\alpha=2,s=\frac{1}{2})$ where $\alpha$ and $s$ are the shape and scale parameters respectively
-   $rho2 \sim U(-1,1)$

Additional prior distributions must be specified for the prevalence in each study `r if(Cond.Dep == TRUE){"and for the covariance parameters"}` as followed

-   $prev \sim \beta(1,1)$, for $i \in 1, ..., n$ where $n$ is the number of studies in the meta-analysis
-   `r if (Cond.Dep == TRUE) {covp_print}`
-   `r if (Cond.Dep == TRUE) {covn_print}`

\*\* Seed value for reproducibility : \*\* `r seed`

## **POSTERIOR RESULTS**

Posterior `mean`, standard deviation (`sd`) as well as posterior median (`50%`) and the 95% credible interval (`2.5%` and `97.5%`) statistics are presented below. Convergence statistics are also provided.

-   `Rhat` is the Gelman-Rubin statistic (Gelman and Ruben(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.

-   `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 resulting in poor posterior estimates.

Posterior estimates are provided for the following parameters :

**Parameters of index test**

-   Predicted sensitivity of index test in a future sutdy, noted as Predicted_Se
-   Predicted specificity in a future study, noted as Predicted_Sp
-   Summary sensitivity across all studies, noted as Summary_Se
-   Summary specificity across all studies, noted as Summary_Sp
-   Mean logit-transformed sensitivity, noted as mu\[1\]
-   Mean logit-transformed sensitivity, noted as mu\[2\]
-   Correlation between the mean logit-transformed sensitivity and the mean logit-transformed specificity, noted as rho
-   Sensitivity in individual studies, noted as se\[i\], where i is the study identifier
-   Specificity in individual studies, noted as sp\[i\], where i is the study identifier
-   Between-study standard deviation in the logit-transformed sensitivity, noted as tau\[1\]
-   Between-study standard deviation in the logit-transformed specificity, noted as tau\[2\]
-   Between-study variance in the logit-transformed sensitivity, noted as tau.sq\[1\]
-   Between-study variance in the logit-transformed specificity, noted as tau.sq\[2\]

**Parameters of reference test**

-   Predicted sensitivity in a future sutdy, noted as Predicted_Se2
-   Predicted specificity in a future study, noted as Predicted_Sp2
-   Summary sensitivity across all studies, noted as Summary_Se2
-   Summary specificity across all studies, noted as Summary_Sp2
-   Mean logit-transformed sensitivity, noted as mu2\[1\]
-   Mean logit-transformed sensitivity, noted as mu2\[2\]
-   Correlation between the mean logit-transformed sensitivity and the mean logit-transformed specificity, noted as rho2
-   Sensitivity in individual studies, noted as se2\[i\], where i is the study identifier
-   Specificity in individual studies, noted as sp2\[i\], where i is the study identifier
-   Between-study standard deviation in the logit-transformed sensitivity, noted as tau2\[1\]
-   Between-study standard deviation in the logit-transformed specificity, noted as tau2\[2\]
-   Between-study variance in the logit-transformed sensitivity, noted as tau.sq2\[1\]
-   Between-study variance in the logit-transformed specificity, noted as tau.sq2\[2\]

**Other parameters**

-   Prevalence in each study, noted as prev `r if(Cond.Dep==TRUE){"* Covariance measuring conditional dependence among disease positive, noted covp"}` `r if(Cond.Dep==TRUE){"* Covariance measuring conditional dependence among disease negative, noted covn"}`

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
res = MCMCsummary(output, digits=digits)
datatable(res, extensions = 'AutoFill')
```

### POSTERIOR ESTIMATES NEEDED TO CREATE SROC PLOT IN RevMan

The following posterior estimates are required in order to create the SROC plot in RevMan :

-   E(logitSe) : Mean logit-transformed sensitivity, mu\[1\] = `r round(mean(as.matrix(output[, "mu[1]"])), digits=digits)`
-   E(logitSp) : Mean logit-transformed sensitivity, mu\[2\] = `r round(mean(as.matrix(output[, "mu[2]"])), digits=digits)`
-   Var(logitSe) : Between-study variance in the logit-transformed sensitivity, tau.sq\[1\] = `r round(mean(as.matrix(output[, "tau.sq[1]"])), digits=digits)`
-   Var(logitSp) : Between-study variance in the logit-transformed specificity, tau.sq\[2\] = `r round(mean(as.matrix(output[, "tau.sq[2]"])), digits=digits)`
-   Corr(logits) : Correlation between the mean logit-transformed sensitivity and the mean logit-transformed specificity, rho = `r round(mean(as.matrix(output[, "rho"])), digits=digits)`
-   Cov(Es) : Covariance between the posterior samples of the mean logit-transformed sensitivity and mean logit-transformed specificity, `r round(cov(as.matrix(output[, "mu[1]"]), as.matrix(output[, "mu[2]"])), digits=digits)`
-   SE(E(logitSe)) : Estimated standard error of mean logit-transformed sensitivity = `r round(sd(as.matrix(output[, "mu[1]"])), digits=digits)`
-   SE(E(logitSp)) : Estimated standard error of mean logit-transformed specificity = `r round(sd(as.matrix(output[, "mu[2]"])), digits=digits)`

## **SUMMARY ROC PLOT**

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
SROC_rjags(X=output, model="Bivariate",n=N, study_col1="blue", study_col2=rgb(0, 0, 1, 0.15), dataset=cell, ref_std=FALSE,SROC_curve = F)

```

## **CONVERGENCE DIAGNOSTIC TOOLS**

Visual inspection of convergence for key parameters can be studied below. For a given parameter, panel (a) shows posterior density plot; (b) the running posterior mean value; and (c) the history plot. Each chain is identified by a different color. Similar results from all`` `r num.chains ` `` chains would suggest the algorithm has converged.

### **INDEX TEST**

#### PREDICTED SENSITIVITY AND SPECIFICITY IN A FUTURE STUDY

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(i in 9:10) {
  # tiff(paste(parameters[i],".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 SENSITIVITY AND SPECIFICITY ACROSS ALL STUDIES

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(i in 7:8) {
  # tiff(paste(parameters[i],".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()
}

```

#### MEAN LOGIT-TRANSFORMED SENSITIVITY AND SPECIFICITY

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(j in 1:2) {
  for(i in c(3)) {
    #  tiff(paste(parameters[i],"[",j,"].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()
  }
}

```

#### CORRELATION BETWEEN MEAN LOGIT-TRANSFORMED SENSITIVITY AND MEAN LOGIT-TRANSFORMED SPECIFICITY

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(i in 6) {
  # tiff(paste(parameters[i],".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()
}

```

#### SENSITIVITY IN INDIVIDUAL STUDIES

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
# Plots to check convergence for parameters shared across studies:

for(i in c(1)) { #POSITION OF THE PARAMETER IN THE "parameters" OBJECT
  for(k in 1:N){
    #  tiff(paste(parameters[i],"[",k,"].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],"[",k,"]",sep="")), auto.layout=FALSE, main="(a)", xlab=paste(parameters[i],"[",k,"]",sep=""), ylab="Density")
    rmeanplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(b)")
    title(xlab="Iteration", ylab="Running mean")
    traplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(c)")
    title(xlab="Iteration", ylab=paste(parameters[i],"[",k,"]",sep=""))
    mtext(paste("Diagnostics for ", parameters[i],"[",k,"]",sep=""), side=3, line=1, outer=TRUE, cex=2)
    #  dev.off()
  }
}
```

#### SPECIFICITY IN INDIVIDUAL STUDIES

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
# Plots to check convergence for parameters shared across studies:

for(i in c(2)) { #POSITION OF THE PARAMETER IN THE "parameters" OBJECT
  for(k in 1:N){
    #  tiff(paste(parameters[i],"[",k,"].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],"[",k,"]",sep="")), auto.layout=FALSE, main="(a)", xlab=paste(parameters[i],"[",k,"]",sep=""), ylab="Density")
    rmeanplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(b)")
    title(xlab="Iteration", ylab="Running mean")
    traplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(c)")
    title(xlab="Iteration", ylab=paste(parameters[i],"[",k,"]",sep=""))
    mtext(paste("Diagnostics for ", parameters[i],"[",k,"]",sep=""), side=3, line=1, outer=TRUE, cex=2)
    #  dev.off()
  }
}
```

#### BETWEEN-STUDY VARIANCE IN THE LOGIT-TRANSFORMED SENSITIVITY AND SPECIFICITY

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(j in 1:2) {
  for(i in c(4)) {
    #  tiff(paste(parameters[i],"[",j,"].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()
  }
}

```

#### BETWEEN-STUDY STANDARD DEVIATION IN THE LOGIT-TRANSFORMED SENSITIVITY AND SPECIFICITY

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(j in 1:2) {
  for(i in c(5)) {
    #  tiff(paste(parameters[i],"[",j,"].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()
  }
}

```

#### BETWEEN-STUDY PRECISION IN THE LOGIT-TRANSFORMED SENSITIVITY AND SPECIFICITY

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(j in 1:2) {
  for(i in c(11)) {
    #  tiff(paste(parameters[i],"[",j,"].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()
  }
}

```

### **REFERENCE TEST**

#### PREDICTED SENSITIVITY AND SPECIFICITY IN A FUTURE STUDY

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(i in 20:21) {
  # tiff(paste(parameters[i],".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 SENSITIVITY AND SPECIFICITY ACROSS ALL STUDIES

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(i in 18:19) {
  # tiff(paste(parameters[i],".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()
}

```

#### MEAN LOGIT-TRANSFORMED SENSITIVITY AND SPECIFICITY

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(j in 1:2) {
  for(i in c(14)) {
    #  tiff(paste(parameters[i],"[",j,"].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()
  }
}

```

#### CORRELATION BETWEEN MEAN LOGIT-TRANSFORMED SENSITIVITY AND MEAN LOGIT-TRANSFORMED SPECIFICITY

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(i in 17) {
  # tiff(paste(parameters[i],".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()
}

```

#### SENSITIVITY IN INDIVIDUAL STUDIES

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
# Plots to check convergence for parameters shared across studies:

for(i in c(12)) { #POSITION OF THE PARAMETER IN THE "parameters" OBJECT
  for(k in 1:N){
    #  tiff(paste(parameters[i],"[",k,"].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],"[",k,"]",sep="")), auto.layout=FALSE, main="(a)", xlab=paste(parameters[i],"[",k,"]",sep=""), ylab="Density")
    rmeanplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(b)")
    title(xlab="Iteration", ylab="Running mean")
    traplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(c)")
    title(xlab="Iteration", ylab=paste(parameters[i],"[",k,"]",sep=""))
    mtext(paste("Diagnostics for ", parameters[i],"[",k,"]",sep=""), side=3, line=1, outer=TRUE, cex=2)
    #  dev.off()
  }
}
```

#### SPECIFICITY IN INDIVIDUAL STUDIES

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
# Plots to check convergence for parameters shared across studies:

for(i in c(13)) { #POSITION OF THE PARAMETER IN THE "parameters" OBJECT
  for(k in 1:N){
    #  tiff(paste(parameters[i],"[",k,"].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],"[",k,"]",sep="")), auto.layout=FALSE, main="(a)", xlab=paste(parameters[i],"[",k,"]",sep=""), ylab="Density")
    rmeanplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(b)")
    title(xlab="Iteration", ylab="Running mean")
    traplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(c)")
    title(xlab="Iteration", ylab=paste(parameters[i],"[",k,"]",sep=""))
    mtext(paste("Diagnostics for ", parameters[i],"[",k,"]",sep=""), side=3, line=1, outer=TRUE, cex=2)
    #  dev.off()
  }
}
```

#### BETWEEN-STUDY VARIANCE IN THE LOGIT-TRANSFORMED SENSITIVITY AND SPECIFICITY

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(j in 1:2) {
  for(i in c(15)) {
    #  tiff(paste(parameters[i],"[",j,"].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()
  }
}

```

#### BETWEEN-STUDY STANDARD DEVIATION IN THE LOGIT-TRANSFORMED SENSITIVITY AND SPECIFICITY

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(j in 1:2) {
  for(i in c(16)) {
    #  tiff(paste(parameters[i],"[",j,"].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()
  }
}

```

#### BETWEEN-STUDY PRECISION IN THE LOGIT-TRANSFORMED SENSITIVITY AND SPECIFICITY

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(j in 1:2) {
  for(i in c(22)) {
    #  tiff(paste(parameters[i],"[",j,"].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()
  }
}

```

### **PREVALENCE IN EACH STUDY**

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
# Plots to check convergence for parameters shared across studies:

for(i in c(23)) { #POSITION OF THE PARAMETER IN THE "parameters" OBJECT
  for(k in 1:N){
    #  tiff(paste(parameters[i],"[",k,"].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],"[",k,"]",sep="")), auto.layout=FALSE, main="(a)", xlab=paste(parameters[i],"[",k,"]",sep=""), ylab="Density")
    rmeanplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(b)")
    title(xlab="Iteration", ylab="Running mean")
    traplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(c)")
    title(xlab="Iteration", ylab=paste(parameters[i],"[",k,"]",sep=""))
    mtext(paste("Diagnostics for ", parameters[i],"[",k,"]",sep=""), side=3, line=1, outer=TRUE, cex=2)
    #  dev.off()
  }
}

```

`r if(Cond.Dep == TRUE){"### **COVARIANCE AMONG DISEASE POSITIVE**"}`

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
# Plots to check convergence for parameters shared across studies:

if(Cond.Dep == TRUE){
  
for(i in c(24)) { #POSITION OF THE PARAMETER IN THE "parameters" OBJECT
  for(k in 1:N){
    #  tiff(paste(parameters[i],"[",k,"].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],"[",k,"]",sep="")), auto.layout=FALSE, main="(a)", xlab=paste(parameters[i],"[",k,"]",sep=""), ylab="Density")
    rmeanplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(b)")
    title(xlab="Iteration", ylab="Running mean")
    traplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(c)")
    title(xlab="Iteration", ylab=paste(parameters[i],"[",k,"]",sep=""))
    mtext(paste("Diagnostics for ", parameters[i],"[",k,"]",sep=""), side=3, line=1, outer=TRUE, cex=2)
    #  dev.off()
  }
}
  
}
```

`r if(Cond.Dep == TRUE){"### **COVARIANCE AMONG DISEASE NEGATIVE**"}`

```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
# Plots to check convergence for parameters shared across studies:

if(Cond.Dep == TRUE){
  
for(i in c(25)) { #POSITION OF THE PARAMETER IN THE "parameters" OBJECT
  for(k in 1:N){
    #  tiff(paste(parameters[i],"[",k,"].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],"[",k,"]",sep="")), auto.layout=FALSE, main="(a)", xlab=paste(parameters[i],"[",k,"]",sep=""), ylab="Density")
    rmeanplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(b)")
    title(xlab="Iteration", ylab="Running mean")
    traplot(output, parms=c(paste(parameters[i],"[",k,"]",sep="")), auto.layout=FALSE, main="(c)")
    title(xlab="Iteration", ylab=paste(parameters[i],"[",k,"]",sep=""))
    mtext(paste("Diagnostics for ", parameters[i],"[",k,"]",sep=""), side=3, line=1, outer=TRUE, cex=2)
    #  dev.off()
  }
}
  
}
```

## **REFERENCES**

-   Brooks, S. P., and Gelman, A. 1998. "General methods for monitoring convergence of iterative simulations." Journal of Computational and Graphical Statistics, no. 7: 434--55.

`r if(name_of_data=="Saliva_COVID"){"* Butler-Laporte, G, Lawandi, A, Schiller, I, Yao, M, Dendukuri, N, McDOnald, E. G. and Lee, T. C. 2021. “Comparison of Saliva and Nasopharyngeal Swab Nucleic Acid Amplification Testing for Detection of SARS-CoV-2; A Systematic Review and Meta-analysis.” JAMA Intern Med."}`

-   Gelman, A, Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A and Rbin, D. B. 2013. "Bayesian Data Analysis." Chapman & Hall/CRC Press, London, Third Edition.

-   Gelman, A, and Rubin, D. B. 1992. "Inference from iterative simulation using multiple sequences." Statistical Sciences, no. 7: 457--72.

-   Xie, X, A Sinclair, and N Dendukuri. 2017. "Evaluating the accuracy and economic value of a new test in the absence of a perfect reference test." Res Synth Methods, no. 8(3): 321--32.
