# Bayesian model for individual cow level fecal samples testing for Salmonella species in dairy cow 
# Two dependent tests, six populations.
# Estimate Se and Sp of fecal culture (CUL) and PCR screening followed by culture (PCR-CUL)

model{
  y1[1:Q, 1:Q] ~ dmulti(p1[1:Q, 1:Q], n1)
  y2[1:Q, 1:Q] ~ dmulti(p2[1:Q, 1:Q], n2)
  y3[1:Q, 1:Q] ~ dmulti(p3[1:Q, 1:Q], n3)
  y4[1:Q, 1:Q] ~ dmulti(p4[1:Q, 1:Q], n4)  
  y5[1:Q, 1:Q] ~ dmulti(p5[1:Q, 1:Q], n5)
  y6[1:Q, 1:Q] ~ dmulti(p6[1:Q, 1:Q], n6)
  
  p1[1,1] <- pi1*eta11 + (1-pi1)*theta11
  p1[1,2] <- pi1*eta12 + (1-pi1)*theta12
  p1[2,1] <- pi1*eta21 + (1-pi1)*theta21
  p1[2,2] <- pi1*eta22 + (1-pi1)*theta22
  
  p2[1,1] <- pi2*eta11 + (1-pi2)*theta11
  p2[1,2] <- pi2*eta12 + (1-pi2)*theta12
  p2[2,1] <- pi2*eta21 + (1-pi2)*theta21
  p2[2,2] <- pi2*eta22 + (1-pi2)*theta22
  
  p3[1,1] <- pi3*eta11 + (1-pi3)*theta11
  p3[1,2] <- pi3*eta12 + (1-pi3)*theta12
  p3[2,1] <- pi3*eta21 + (1-pi3)*theta21
  p3[2,2] <- pi3*eta22 + (1-pi3)*theta22
  
  p4[1,1] <- pi4*eta11 + (1-pi4)*theta11
  p4[1,2] <- pi4*eta12 + (1-pi4)*theta12
  p4[2,1] <- pi4*eta21 + (1-pi4)*theta21
  p4[2,2] <- pi4*eta22 + (1-pi4)*theta22
  
  p5[1,1] <- pi5*eta11 + (1-pi5)*theta11
  p5[1,2] <- pi5*eta12 + (1-pi5)*theta12
  p5[2,1] <- pi5*eta21 + (1-pi5)*theta21
  p5[2,2] <- pi5*eta22 + (1-pi5)*theta22
  
  p6[1,1] <- pi6*eta11 + (1-pi6)*theta11
  p6[1,2] <- pi6*eta12 + (1-pi6)*theta12
  p6[2,1] <- pi6*eta21 + (1-pi6)*theta21
  p6[2,2] <- pi6*eta22 + (1-pi6)*theta22
  
  eta11 <- lambdaD*eta1
  eta12 <- eta1 - eta11
  eta21 <- gammaD*(1-eta1)
  eta22 <- 1 - eta11 - eta12 - eta21
  theta11 <- 1 - theta12 - theta21 - theta22
  theta12 <- gammaDc*(1-theta1)
  theta21 <- theta1 - theta22
  theta22 <- lambdaDc*theta1
  eta2 <- eta11 + eta21           #Sensitivity test 2
  theta2 <- theta22 + theta12     #Specificity test 2
  rhoD <- (eta11 - eta1*eta2) / sqrt(eta1*(1-eta1)*eta2*(1-eta2))
  rhoDc <- (theta22 - theta1*theta2) / sqrt(theta1*(1-theta1)*theta2*(1-theta2))
 
 # Priors
  pi1 ~ dbeta(1, 1)
  pi2 ~ dbeta(1, 1)
  pi3 ~ dbeta(1, 1)
  pi4 ~ dbeta(1, 1)
  pi5 ~ dbeta(1, 1)
  pi6 ~ dbeta(1, 1)  
  eta1 ~ dbeta(1, 1)             #Sensitivity test 1
  theta1 ~ dbeta(212.12, 3.13)   #Specificity test 1, Mode=0.99, 95% sure Spcul > 0.97
  
 # Priors for conditional Se2 and Sp2
 #PCR sensitivity estimates
  lambdaD ~ dbeta(1, 1)   ### dbeta(6.25, 1.92), Mode=0.85, 95%tile=0.50
  gammaD ~ dbeta(1, 1)   ### dbeta(6.25, 1.92)

 #PCR specificity estimates  
 lambdaDc ~ dbeta(1, 1)  ###  dbeta(4.77, 1.19), Mode=0.95, 95%tile=0.50
 gammaDc ~ dbeta(1, 1)  ###  dbeta(4.77, 1.19)
  
# Renaming variables 
  Secul <- eta1
  Spcul <- theta1
  Sepcr <- eta2
  Sppcr <- theta2
}

# Data
list(n1=40, n2=40, n3=40, n4=40, n5=40, n6=40, Q=2,
y1=structure(.Data=c(0,1,0,39),.Dim=c(2,2)), y2=structure(.Data=c(0,1,0,39),.Dim=c(2,2)), 
y3=structure(.Data=c(3,1,0,36),.Dim=c(2,2)), y4=structure(.Data=c(21,1,0,18),.Dim=c(2,2)), 
y5=structure(.Data=c(2,1,0,37),.Dim=c(2,2)), y6=structure(.Data=c(28,1,0,11),.Dim=c(2,2)))

# Data (EXPANDED)
list(n1=140, n2=70, n3=110, n4=190, n5=60, n6=100, Q=2,
y1=structure(.Data=c(0,4,0,136),.Dim=c(2,2)), y2=structure(.Data=c(0,2,0,68),.Dim=c(2,2)), 
y3=structure(.Data=c(11,4,0,95),.Dim=c(2,2)), y4=structure(.Data=c(104,6,0,80),.Dim=c(2,2)), 
y5=structure(.Data=c(3,2,0,55),.Dim=c(2,2)), y6=structure(.Data=c(65,3,0,32),.Dim=c(2,2)))

#Initial values
list(pi1=0.15, pi2=0.30, pi3=0.13, pi4=0.65, pi5=0.23, pi6=0.73, eta1=0.55, theta1=0.98, lambdaD=0.50, lambdaDc=0.50, gammaD=0.50, gammaDc=0.50)
list(pi1=0.05, pi2=0.03, pi3=0.30, pi4=0.10, pi5=0.03, pi6=0.30, eta1=0.75, theta1=0.80, lambdaD=0.30, lambdaDc=0.30, gammaD=0.70, gammaDc=0.70)
list(pi1=0.20, pi2=0.13, pi3=0.05, pi4=0.25, pi5=0.50, pi6=0.15, eta1=0.65, theta1=0.75, lambdaD=0.70, lambdaDc=0.70, gammaD=0.30, gammaDc=0.30)