# latent class model with random effects for two observed tests
# applied to Strongyloides infection dataset

library(rjags)
library(mcmcplots)

modelString = 
"model
{

#==========================================================
# Likelihood
#==========================================================

for (i in 1:N) { # N=total sample size

	D1[i]~dbern(prev)  # D1=True (latent) disease status of ith subject, prev=prevalence
	D[i]<-D1[i]+1       

	t1[i]~dbern(p1[i,D[i]])  # t1=result of test 1
	t2[i]~dbern(p2[i,D[i]])  # t2=result of test 2
	
	r[i]~dnorm(0,1) # random effect
	p1[i,2]<-phi(a[1,2]+b[2]*r[i]) # sensitivity of test 1
	p1[i,1]<-phi(-a[1,1]-b[1]*r[i]) # 1-specificity of test 1
	p2[i,2]<-phi(a[2,2]+b[2]*r[i]) # sensitivity of test 2
	p2[i,1]<-phi(-a[2,1]-b[1]*r[i])# 1-specificity of test 2

}

#==================================================
# Prior distributions
#==================================================

prev~dbeta(1,1)

S[1]~dbeta(4.44,13.31)  # Informative prior for sensitivity of microscopy
C[1]~dbeta(71.25,3.75)  # Informative prior for specificty of microscopy
S[2]~dbeta(21.96,5.49)  # Informative prior for sensitivity of serology
C[2]~dbeta(4.1,1.76)    # Informative prior for specificty of serology

b[1]~dunif(0,3)
b[2]~dunif(0,3)

a[1,2]<-probit(S[1])*sqrt(1+b[2]*b[2])
a[2,2]<-probit(S[2])*sqrt(1+b[2]*b[2])
a[1,1]<-probit(C[1])*sqrt(1+b[1]*b[1])
a[2,1]<-probit(C[2])*sqrt(1+b[1]*b[1])

}"  
  
writeLines(modelString,con="model_Strongy.txt")  

#  Dataset
t1=rep(c(1,0),c(40,122))
t2=rep(c(1,0,1,0),c(38,2,87,35))
dataList=list(N=162,t1=t1,t2=t2)

jagsModel = jags.model("model_Strongy.txt",data=dataList,n.chains=3)  
update(jagsModel,n.iter=5000)

parameters = c("S","C","prev","a","b")  

output = coda.samples(jagsModel,variable.names=parameters,n.iter=10000)

traplot(output)
denplot(output)
summary(output)
gelman.diag(output)