---
title: "Bayesian Latent Class Fixed Effects 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 = "Strongyloides_data" # 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=5000                       # Number of iterations to drop during the "burn-in" phase
num.iter=5000                       # Number of iterations to keep to form the posterior samples
num.chains=3                        # Number of independent chains to run


############################################## PRIOR DISTRIBUTIONS ###############################################################
###                                                                                                                            ###
# A Beta prior distribution with parameter "alpha.prev" and "beta.prev" is used for 
# the prevalence parameter.  Assigning "alpha.prev = beta.prev = 1" is equivalent to
# uniform(0,1) weak prior
 alpha.prev = 1
 beta.prev = 1
###                                                                                                                            ###
# A Beta prior distribution Beta(alpha, beta) is used for the sensitivity and specificity 
# of each tests.  The Beta distribution parameters alpha and beta for each tests are given
# as followed:
 
  alpha.se2=4.44 # alpha parameter of the Beta prior distribution for the sensitivity of the reference test
  beta.se2=13.31 # beta  parameter of the Beta prior distribution for the sensitivity of the reference test
  alpha.sp2=71.25 # alpha parameter of the Beta prior distribution for the specificity of the reference test
  beta.sp2=3.75 # beta parameter of the Beta prior distribution for the specificity of the reference test
  alpha.se1=21.96 # alpha parameter of the Beta prior distribution for the sensitivity of the index test
  beta.se1=5.49 # beta parameter of the Beta prior distribution for the sensitivity of the index test
  alpha.sp1=4.1 # alpha parameter of the Beta prior distribution for the specificity of the index test
  beta.sp1=1.76 # alpha parameter of the Beta prior distribution for the specificity of the index test
###                                                                                                                            ###  
  alpha.covs12=3 # alpha parameter of the Beta prior distribution for covariance between both tests among the diseased subjects
  beta.covs12=3 # beta parameter of the Beta prior distribution for covariance between both tests among the diseased subjects
  alpha.covc12=3 # alpha parameter of the Beta prior distribution for covariance between both tests among the nondiseased subjects
  beta.covc12=3 # beta parameter of the Beta prior distribution for covariance between both tests among the nondiseased subjects
 
###                                                                                                                            ###
##################################################################################################################################     
    
################################  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 7 parameters :                 ###
#                                                                                                      ###
# Disease prevalence denoted by prev                                                                   ###
# Sensitivity of the test under evaluation denoted by se1                                              ###
# Specificity of the test under evaluation denoted by sp1                                              ###
# Sensitivity of the reference test denoted by se2                                                     ###
# Specificity of the reference test denoted by sp2                                                     ###
# Covariance between both tests among diseased subjects denoted by u.covs12                            ###
# Covariance between both tests among nondiseased subjects denoted by u.covc12                         ###
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 7 parameters described above.                           ###
#                                                                                                      ###
prev=rep(0.9, num.chains)
se1=rep(0.8, num.chains)
sp1=c(0.9, 0.85, 0.95)
se2=c(0.5, 0.3, 0.7)
sp2=rep(.95, num.chains)
u.covs12=rep(0.5, num.chains)
u.covc12=rep(0.5, num.chains)
#                                                                                                      ###
# 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")
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.
library(mvtnorm)  # FOR GENERATION OF MULTIVARIATE NORMAL NUMBERS IN Geninits FUNCTION
```



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



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

y <- c(DATA$n11, DATA$n10, DATA$n01, DATA$n00)

N = sum(y)


dataList = list(y=y,N=N,alpha.se2=alpha.se2,beta.se2=beta.se2,alpha.sp2=alpha.sp2,beta.sp2=beta.sp2,alpha.se1=alpha.se1,beta.se1=beta.se1,alpha.sp1=alpha.sp1,beta.sp1=beta.sp1,alpha.covs12=alpha.covs12,beta.covs12=beta.covs12,alpha.covc12=alpha.covc12,beta.covc12=beta.covc12,alpha.prev=alpha.prev,beta.prev=beta.prev)

```

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


The `r paste(name_of_data)` database `r if(name_of_data=="Strongyloides_data"){"is a toy example provided to run this script.  It's about a study conducted to estimate the prevalence of Strongyloides infection among a group of Cambodian refugees to Canada (Joseph, Gyorkos, and Coupal (1995)).  It"}` includes `r sum(y)` participants with diagnostic results on 2 tests. `r if(name_of_data=="Strongyloides_data"){"From a notation point of view, we suppose here that the Stool examination is the reference test and the Serology test is the index test.  The data is"}` as followed 

* n11 cell = Number of patients positive on both tests
* n10 cell = Number of patients positivie on first test (index test) and negative on second test (reference test)
* n01 cell = Number of patients  negative on first test (index test) and positive on second test (reference test)
* n00 cell = Number of patients negative on both tests




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



{if(your.inits==FALSE){          
      prior.list = c(alpha.se2,beta.se2,alpha.sp2,beta.sp2,alpha.se1,beta.se1,alpha.sp1,beta.sp1,alpha.covs12,beta.covs12,alpha.covc12,beta.covc12,alpha.prev,beta.prev)
}
else{prior.list=NULL}
}



GenInits  = function(Seed=seed, prev=prev, se1=se1, sp1=sp1, se2=se2, sp2=sp2, u.covs12=u.covs12, u.covc12=u.covc12, inits=your.inits, prior=prior.list){
  
  
  if(inits == TRUE){
    se = c(se1,se2)
    sp = c(sp1,sp2)
    prev = prev
    u.covs12=u.covs12
    u.covc12=u.covc12
    list(
      se = se,
      sp = sp,
      prev = prev,
      u.covs12=u.covs12,
      u.covc12=u.covc12
    )
  }
  else{
    # INITIAL VALUES RANDOMLY SELECTED BASED ON THE PRIOR DISTRIBUTION OF EACH PARAMETERS
    
    se2 <- rbeta(1,shape1=prior[1],shape2=prior[2])
    sp2 <- rbeta(1,shape1=prior[3],shape2=prior[4])
    se1 <- rbeta(1,shape1=prior[5],shape2=prior[6])
    sp1 <- rbeta(1,shape1=prior[7],shape2=prior[8])
    u.covs12 <- rbeta(1,shape1=prior[9],shape2=prior[10])
    u.covc12 <- rbeta(1,shape1=prior[11],shape2=prior[12])
    se = c(se1,se2)
    sp = c(sp1,sp2)
    #alpha.prev = prior[13]
    #beta.prev = prior[14]
    prev <- rbeta(1,shape1=prior[13],shape2=prior[14])
    prev = prev
    list(
      se = se,
      sp = sp,
      prev = prev,
      u.covs12=u.covs12,
      u.covc12=u.covc12,
      .RNG.name="base::Wichmann-Hill",
      .RNG.seed=Seed
    )
  }
}


Listinits = vector('list',num.chains)
for(i in 1:num.chains){
    tmp = (GenInits(Seed=seed, prev=prev[i], se1=se1[i], sp1=sp1[i], se2=se2[i], sp2=sp2[i], u.covs12=u.covs12[i], u.covc12=u.covc12[i], inits=your.inits, prior=prior.list))
    Listinits[[i]] = tmp
}

```


## **BAYESIAN LATENT CLASS RANDOM EFFECTS MODEL**

```{r, eval=T, echo=T, message=FALSE }
# Bayesian bivariate model

modelString =

"model {

  # LIKELIHOOD
	y[1:4]~dmulti(p12[1:4],N)
	
	
	# probabilities of observing different cross-classifications of two dichotomous diagnostic tests
	p12[1]<- prev*(se[1]*se[2]+covs12)+(1-prev)*((1-sp[1])*(1-sp[2])+covc12)
  p12[2]<- prev*(se[1]*(1-se[2])-covs12)+(1-prev)*((1-sp[1])*sp[2]-covc12)
  p12[3]<- prev*((1-se[1])*se[2]-covs12)+(1-prev)*(sp[1]*(1-sp[2])-covc12)
  p12[4]<- prev*((1-se[1])*(1-se[2])+covs12)+(1-prev)*(sp[1]*sp[2]+covc12)
	

  #=======================================
  # upper limits of covariance parameters
  #=======================================
  us<- min(se[1],se[2])-(se[1]*se[2])
  uc<- min(sp[1],sp[2])-(sp[1]*sp[2])
  
  #==============================================
  # adjustment of range of covariance parameters
  #==============================================
  covs12<- u.covs12*us
  covc12<- u.covc12*uc

  #==================================================
  # Prior distributions of prevalence, sensitivities
  # and specificities
  #==================================================
  prev~dbeta(alpha.prev,beta.prev)
  se[1]~dbeta(alpha.se1,beta.se1)
  sp[1]~dbeta(alpha.sp1,beta.sp1)
  se[2]~dbeta(alpha.se2,beta.se2)
  sp[2]~dbeta(alpha.sp2,beta.sp2)
  
  #==============================================================
  # prior distribution of transformed covariances on (0,1) range
  #==============================================================
  u.covs12~ dbeta(alpha.covs12,beta.covs12)
  u.covc12~ dbeta(alpha.covc12,beta.covc12)

}"


  

  

```


The posterior samples were obtained by running a Bayesian latent-class fixed effects model (Dendukuri and Joseph, 2001) 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**

Prior distribution for the prevalence parameter follows a Beta(alpha.prev,beta.prev) distribution where

* alpha.prev = `r alpha.prev`
* beta.prev = `r beta.prev`

Prior distribution for the sensitivity of the test under evaluation follows a Beta(alpha.se1,beta.se1) distribution where

* alpha.se1 = `r alpha.se1`
* beta.se1 = `r beta.se1`

Prior distribution for the specificity of the test under evaluation follows a Beta(alpha.sp1,beta.sp1) distribution where

* alpha.sp1 = `r alpha.sp1`
* beta.sp1 = `r beta.sp1`

Prior distribution for the sensitivity of the reference test follows a Beta(alpha.se2,beta.se2) distribution where

* alpha.se2 = `r alpha.se2`
* beta.se2 = `r beta.se2`

Prior distribution for the specificity of the reference test follows a Beta(alpha.sp2,beta.sp2) distribution where

* alpha.sp2 = `r alpha.sp2`
* beta.sp2 = `r beta.sp2`

Prior distribution for the covariance between the two tests among diseased subjects follows a Beta(alpha.covs12,beta.covs12) distribution where

* alpha.covs12 = `r alpha.covs12`
* beta.covs12 = `r beta.covs12`

Prior distribution for the covariance between the two tests among nondiseased subjects follows a Beta(alpha.covc12,beta.covc12) distribution where

* alpha.covc12 = `r alpha.covc12`
* beta.covc12 = `r beta.covc12`


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


## **INITIAL VALUES**


```{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, inits=Listinits)
}else{
  jagsModel = jags.model("model.txt",data=dataList,n.chains=num.chains,n.adapt=0, inits=Listinits)
}
```


`r ifelse(your.inits==FALSE, "Initial values were generated based on the prior distributions and were as followed" , "Initial values were provided by the user and were as followed" )`


```{r, eval=T, echo=T,  message=FALSE }

for (i in 1:num.chains) { 
  print("############################")
  print(paste("#INITIAL VALUES FOR CHAIN ",i,"#",sep=""))
  print("############################")
  se_parm = matrix(unlist(jagsModel$state()[[i]][2]),byrow=FALSE,ncol=2)
  colnames(se_parm) <- c("se[1]", "se[2]")
  rownames(se_parm) <- c("sensitivity")
  #datatable(aparm, extensions = 'AutoFill', options = list(autoFill = TRUE)
  print(" SENSITIVITY PARAMETERS")
  print(se_parm)
  print("****************************")
  sp_parm = matrix(unlist(jagsModel$state()[[i]][3]),byrow=FALSE,ncol=2)
  colnames(sp_parm) <- c("sp[1]", "sp[2]")
  rownames(sp_parm) <- c("sêcificity")
  print(" SPECIFICITY PARAMETERS")
  print(sp_parm)
  print("****************************")
  prevparm = unlist(jagsModel$state()[[i]][1])
  print(" PREVALENCE")
  print(prevparm)
  print("****************************")
  covsparm = unlist(jagsModel$state()[[i]][5])
  print(" COVARIANCE AMONG DISEASED SUBJECTS")
  print(covsparm)
  print("****************************")
  covcparm = unlist(jagsModel$state()[[i]][4])
  print(" COVARIANCE AMONG NONDISEASED SUBJECTS")
  print(covcparm)
  #print("******************************")
}
```

## **POSTERIOR RESULTS** 


```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE, results='hide' }
#jagsModel$state(internal=FALSE)

# Burn-in iterations 
update(jagsModel,n.iter=num.burn)
# Parameters to be monitored
parameters = c( "se","sp", "covs12", "covc12", "prev")

# Posterior samples
posterior_results = coda.samples(jagsModel,variable.names=parameters,n.iter=num.iter)
output = 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** `r if(name_of_data=="Strongyloides_data"){"Serology test"}`

* Sensitivity, denoted as se[1]
* Specificity, denoted as sp[1]

**Parameters of reference test** `r if(name_of_data=="Strongyloides_data"){"Stool examination"}`

* Sensitivity, denoted as se[2]
* Specificity, denoted as sp[2]

**Disease prevalence**

* Denoted as prev

**Covariance between the two tests among diseased subjects**

* Denoted as covs12

**Covariance between the two tests among nondiseased subjects**

* Denoted as covc12

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



## **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** `r if(name_of_data=="Strongyloides_data"){"SEROLOGY TEST"}`

#### SENSITIVITY AND SPECIFICITY
```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(k in 1) {
  for(i in 1:2) {
    #  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],"[",k,"]",sep="")), auto.layout=FALSE, main="(a)", xlab=paste(parameters[i],"",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()
  }
}


```




### **REFERENCE TEST** `r if(name_of_data=="Strongyloides_data"){"STOOL EXAMINATION"}`

#### SENSITIVITY AND SPECIFICITY
```{r, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
for(k in 2) {
  for(i in 1:2) {
    #  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],"[",k,"]",sep="")), auto.layout=FALSE, main="(a)", xlab=paste(parameters[i],"",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()
  }
}

```



### **PREVALENCE**

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

```

### **COVARIANCE BETWEEN THE TWO TESTS AMONG DISEASED SUBJECTS**

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

```

### **COVARIANCE BETWEEN THE TWO TESTS AMONG NONDISEASED SUBJECTS**

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

```

## **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=="Strongyloides_data"){"* Joseph, L, Gyorkos, T, and Coupal, L. 1995 “Bayesian estimation of disease prevalence and the parameters of diagnostic tests in the absence of a gold standard” American Journal of Epidemiology 141, 263-272"}`

* 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.

* Dendukuri, N, Joseph, L. "Bayesian approaches to modeling the conditional dependence between multiple diagnostic tests." Biometrics. (2001) 57:158–67.


