---
title: "Bayesian Latent Class Random Effects Model with Gaussian Quadrature"
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_RE" # 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
gauss.quad.pts = 19                    # Number of quadrature points

################################  PRIOR DISTRIBUTION ##################################
###                                                                                 ###
# 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 is used for the sensitivity and specificity of each 
# tests.  The parameters of the beta prior are provided following the approach found
# at https://nandinidendukuri.com/blogposts/2022-07-30-lca-remodel-prior/
#
# For the sensitivity of test 1 (index test) the beta prior is Beta(alpha.se1, beta.se1)
# where
  alpha.se1=21.96  
  beta.se1=5.49
# For the specificity of test 1 (index test) the beta prior is Beta(alpha.sp1, beta.sp1)
# where
  alpha.sp1=4.1  
  beta.sp1=1.76  
# For the sensitivity of test 2 (reference test) the beta prior is Beta(alpha.se2, beta.se2)
# where 
  alpha.se2=4.44  
  beta.se2=13.31
# For the specificity of test 2 (index test) the beta prior is Beta(alpha.sp2, beta.sp2)
# where
  alpha.sp2=71.25  
  beta.sp2=3.75
###                                                                                 ###
#######################################################################################     
    
################################  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 5 parameters :                                  ###
#                                                                                   ###
# Prevalence of the target condition (referred below as prev)                       ###
# Sensitivity of index test (referred below as se1)                                 ###
# Specificity of index test (referred below as sp1)                                 ###
# Sensitivity of reference test (referred below as se2)                             ###
# Specificity of reference test (referred below as sp2)                             ###
#                                                                                   ###
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 5 parameters   ###
# described above.                                                                  ###
#                                                                                   ###
prev=rep(0.9, num.chains)
se1=c(0.5, 0.6, 0.7)
sp1=rep(.9, num.chains)

se2=rep(.5, num.chains)
sp2=c(.9, 0.95, .99)

#                                                                                   ###
# 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)   # THIS PACKAGE CONTAINS THE MCMCsummary FUNCTION USED IN THIS SCRIPT
require(mcmcplots) # USED FOR THE CREATION OF THE CONVERGENCE PLOTS
require(rjags)     # PACKAGE TO RUN THE jags MODEL
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
require(statmod)   # FOR GAUSSIAN QUADRATURE 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$n00, DATA$n01, DATA$n10, DATA$n11)
N = sum(y)
gaussquad <- gauss.quad.prob(n=gauss.quad.pts, dist="normal", mu=0, sigma=1)
n.pattrn <- dim(DATA)[2]
n.gauss.quad = length(gaussquad$nodes)





{if(your.inits==FALSE){          
      dataList = list(y=y, N=N, nodes=gaussquad$nodes, weights=gaussquad$weights, n.pattrn=n.pattrn,
n.gauss.quad = length(gaussquad$nodes), alpha.prev=alpha.prev,beta.prev=beta.prev)
}
else{
  dataList = list(y=y, N=N, nodes=gaussquad$nodes, weights=gaussquad$weights, n.pattrn=n.pattrn,
n.gauss.quad = length(gaussquad$nodes),alpha.se1=alpha.se1,beta.se1=beta.se1,alpha.se2=alpha.se2,beta.se2=beta.se2,alpha.sp1=alpha.sp1,beta.sp1=beta.sp1,alpha.sp2=alpha.sp2,beta.sp2=beta.sp2,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_RE"){"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 dim(y)[1]` participants with diagnostic results on 2 tests. `r if(name_of_data=="Strongyloides_data_RE"){"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.se1,beta.se1,alpha.se2,beta.se2,alpha.sp1,beta.sp1,alpha.sp2,beta.sp2,alpha.prev,beta.prev)
}
else{prior.list=NULL}
}



GenInits  = function(Seed=seed, prev=prev, se1=se1, se2=se2, sp1=sp1, sp2=sp2, inits=your.inits, prior=prior.list){
  
  
  if(inits == TRUE){
    se <- c(se1, se2)
    sp <- c(sp1, sp2)
    b1 <- runif(1,0,3)
    b2 <- runif(1,0,3)
    b = c(b1,b2)
    prev = prev
    list(
      se = se,
      sp = sp,
      b.RE = b,
      prev = prev
    )
  }
  else{
    # INITIAL VALUES RANDOMLY SELECTED BASED ON THE PRIOR DISTRIBUTION OF EACH PARAMETERS

    se1 <- rbeta(1,prior[1], prior[1])
    sp1 <- rbeta(1,prior[5], prior[6])
    se2 <- rbeta(1,prior[3], prior[4])
    sp2 <- rbeta(1,prior[7], prior[8])
    se <- c(se1, se2)
    sp <- c(sp1, sp2)
    b1 <- runif(1,0,3)
    b2 <- runif(1,0,3)
    b = c(b1,b2)
    prev <- rbeta(1,shape1=prior[9],shape2=prior[10])
    prev = prev
    list(
      se = se,
      sp = sp,
      b.RE = b,
      prev = prev,
      .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], se2=se2[i], sp1=sp1[i], sp2=sp2[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 {

#===============================================================================
# Relation between observed data and conditional probability of a positive observation
#===============================================================================



  # LIKELIHOOD
	y[1:n.pattrn] ~ dmulti(q[1:n.pattrn], N)
	
	for(k in 1:n.gauss.quad){ 
  
    GQ[1,k] <- prev*((1-s1[k])*(1-s2[k])) + (1-prev)*c1[k]*c2[k]
    GQ[2,k] <- prev*((1-s1[k])*s2[k]) + (1-prev)*c1[k]*(1-c2[k])
    GQ[3,k] <- prev*(s1[k]*(1-s2[k])) + (1-prev)*(1-c1[k])*c2[k]
    GQ[4,k] <- prev*(s1[k]*s2[k]) + (1-prev)*(1-c1[k])*(1-c2[k])
    
    s1[k] <- phi(a[1,1] + b.RE[1]*nodes[k])
    s2[k] <- phi(a[2,1] + b.RE[2]*nodes[k])
    
    c1[k] <- phi(a[1,2] + b.RE[1]*nodes[k])
    c2[k] <- phi(a[2,2] + b.RE[2]*nodes[k])  
    
	}

  for(l in 1:n.pattrn) {		#PROBABILITY CELL LOOP
			q[l] <- inprod(GQ[l,1:n.gauss.quad], weights[1:n.gauss.quad])/sum(weights[1:n.gauss.quad])
	}
	
	
#==================================================
# Prior distributions of prevalence, sensitivities
# and specificities, as well as for the b-parameters
# associated with the random effects
#==================================================


prev~dbeta(alpha.prev,beta.prev)
se[1]~dbeta(21.96,5.49)
sp[1]~dbeta(4.1,1.76) 
se[2]~dbeta(4.44,13.31)
sp[2]~dbeta(71.25,3.75)
b.RE[1]~dunif(0,3)
b.RE[2]~dunif(0,3)


#=========================================================
# a-parameters used in the 
#=========================================================

a[1,1]<-probit(se[1])*sqrt(1+b.RE[1]*b.RE[1])
a[1,2]<-probit(sp[1])*sqrt(1+b.RE[1]*b.RE[1])
a[2,1]<-probit(se[2])*sqrt(1+b.RE[2]*b.RE[2])
a[2,2]<-probit(sp[2])*sqrt(1+b.RE[2]*b.RE[2])


}"


  

  

```

The posterior samples were obtained by running a Bayesian latent-class random effects model (Dendukuri and Joseph, 2001) using the Gauss-Hermite quadrature ( Qu, Tan and Kutner, 1996) 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`

Beta prior distributions are used for the sensitivity and specificity of both tests as followed (see [Specifying prior distribution for diagnostic test accuracy in latent class model with random effects](https://nandinidendukuri.com/blogposts/2022-07-30-lca-remodel-prior/))

For test 1 (index test `r if(name_of_data=="Strongyloides_data_RE"){"=Serology test"}`) 

* $se[1] \sim Beta(alpha.se1,beta.se1)$ 
* $sp[1] \sim Beta(alpha.sp1,beta.sp1)$

where

-   alpha.se1 = `r alpha.se1`
-   beta.se1 = `r beta.se1`
-   alpha.sp1 = `r alpha.sp1`
-   beta.sp1 = `r beta.sp1`

For test 2 (reference test `r if(name_of_data=="Strongyloides_data_RE"){"=Stool examination"}`) 

* $se[2] \sim Beta(alpha.se2,beta.se2)$ 
* $sp[2] \sim Beta(alpha.sp2,beta.sp2)$

where

-   alpha.se2 = `r alpha.se2`
-   beta.se2 = `r beta.se2`
-   alpha.sp2 = `r alpha.sp2`
-   beta.sp2 = `r beta.sp2`

Uniform prior distributions are used for the random effects parameters as followed

-   $b.RE[1] \sim U(0,3)$
-   $b.RE[2] \sim U(0,3)$

\*\* 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_GQ.txt")

  # Compile the model
if(your.inits==FALSE){
  jagsModel = jags.model("model_GQ.txt",data=dataList,n.chains=num.chains,n.adapt=0, inits=Listinits)
}else{
  jagsModel = jags.model("model_GQ.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]][3]),byrow=FALSE,ncol=2)
  colnames(se_parm) <- c("se[1]", "se[2]")
  rownames(se_parm) <- c("sensitivity")
  print(" SENSITIVITY PARAMETERS")
  print(se_parm)
  print("****************************")
  sp_parm = matrix(unlist(jagsModel$state()[[i]][4]),byrow=FALSE,ncol=2)
  colnames(sp_parm) <- c("sp[1]", "sp[2]")
  rownames(sp_parm) <- c("specificity")
  print(" SPECIFICITY PARAMETERS")
  print(sp_parm)
  print("****************************")
  bparm = as.matrix(unlist(jagsModel$state()[[i]][1]),ncol=2)
  rownames(bparm) <- c("b.RE[1]", "b.RE[2]")
  print(" b-PARAMETERS")
  print(bparm)
  print("****************************")
  prevparm = unlist(jagsModel$state()[[i]][2])
  print(" PREVALENCE")
  print(prevparm)
  #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", "a", "b.RE", "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_RE"){"Serology test"}`

-   Sensitivity, denoted as se[1], which is a function of parameters a[1,1] and b.RE[1]
-   Specificity, denoted as sp[1], which is a function of parameters a[1,2] and b.RE[1]

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

-   Sensitivity, denoted as se[2], which is a function of parameters a[2,1] and b.RE[2]
-   Specificity, denoted as sp[2], which is a function of parameters a[2,2] and b.RE[2]

**Disease prevalence**

-   Denoted as prev

```{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_RE"){"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()
  }
}


```

#### OTHER PARAMETERS

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

for(i in c(4)) { #POSITION OF THE PARAMETER IN THE "parameters" OBJECT
  for(k in 1){
    #  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, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
# Plots to check convergence for parameters shared across studies:

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

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

```

#### OTHER PARAMETERS

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

for(i in c(4)) { #POSITION OF THE PARAMETER IN THE "parameters" OBJECT
  for(k in 2){
    #  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, eval=T, echo=F, class.source = 'fold-hide', message=FALSE }
# Plots to check convergence for parameters shared across studies:

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

```

## **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_RE"){"* 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.

-   Qu, Y., Tan, M., and Kutner, M.H. 1996. "Random effects models in latent class analysis for evaluating accuracy of diagnostic tests." Biometrics, 52:797-810
