---
title: "Bayesian Bivariate 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
```

```{r, eval=T, echo=F, results='hide', message=FALSE, warning=FALSE    }
# install.packages("remotes")
#remotes::install_github("rlesur/klippy")
library(klippy)
library(knitr)
library(plyr)
library(MCMCvis)
library(mcmcplots)
library(rjags)
library(DTAplots)
library(DT)  # THIS LIBRARY ALLOWS THE NICE DATA DISPLAY WITH THE SEARCH BAR OPTION.
```

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

n = dim(DATA)[1] # Number of studies
TP = DATA$tp # True positive
FP = DATA$fp # False positive
FN = DATA$fn # False negative
TN = DATA$tn # True negative
pos=TP + FN
neg=TN + FP 
dataList = list(TP=TP,TN=TN,n=n,pos=pos,neg=neg)

```

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

The `r paste(name_of_data)` database includes `r sum(c(TP, FP, FN, TN))` 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))
```

## **BAYESIAN BIVARIATE META-ANALYSIS MODEL SPECIFICATION**

```{r, eval=T, echo=T, message=FALSE }
# Bayesian bivariate 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 POOLED LOGIT SENSITIVITY AND SPECIFICITY === #
mu[1] ~ dnorm(0,0.25)
mu[2] ~ dnorm(0,0.25) 
# Between-study variance-covariance matrix  

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

# ===  PARAMETERS OF INTEREST === #
# BETWEEN-STUDY STANDARD DEVIATION (tau) AND VARIANCE (tau.sq) OF LOGIT SENSITIVITY AND SPECIFICITY
tau[1]<-pow(prec[1],-0.5)
tau[2]<-pow(prec[2],-0.5)
tau.sq[1]<-pow(prec[1],-1)
tau.sq[2]<-pow(prec[2],-1)
# SUMMARY 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]))
}
"
```

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

# Compile the model 
jagsModel = jags.model("model.txt",data=dataList,n.chains=num.chains)

# Burn-in iterations 
update(jagsModel,n.iter=num.burn)
# 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=num.iter)

```

The posterior samples were obtained by running the bivariate model (Reitsma et al(2005)) 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.

This model uses prior distributions on the log-odds scale over the pooled sensitivity and specificity parameters (normal distribution with mean 0 and a variance of 4), their corresponding between-study precision on the log odds scale (gamma distribution with a shape parameter of 2 and a rate parameter of 0.5), and the correlation terms (uniform distribution whose upper and lower limits are -1 and 1, respectively). The precision is defined as the reciprocal of the variance.

## **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 :

-   Predicted sensitivity 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\]

```{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=cbind(TP,FP,FN,TN), ref_std=TRUE,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.

### 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()
  }
}

```

## **REFERENCES**

-   Brooks, S. P., and A Gelman. 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, J. B. Carlin, H. S. Stern, D. B. Dunson, A Vehtari, and D. B. Rbin. 2013. "Bayesian Data Analysis." *Chapman & Hall/CRC Press, London, Third Edition*.

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

-   Reitsma, J. B., A. S. Glas, A. W. Rutjes, R. J. Scholten, P. M. Bossuyt, and A. H. Zwinderman. 2005. "Bivariate analysis of sensitivity and specificity produces informative summary measures in diagnostic reviews." *Journal of CLinical Epidemiology*, no. 58: 892--990.
