model
{

### inputs are:
### main study data wtot1, n1, wtot2, n2
### validation data xwtot1, m1, xwtot2, m2
### hyperparameters muSN, varSN, rhoSN, muSP, varSP, rhoSP, 
### mur, varr, rhor

# main study
wtot1 ~ dbin(p1, n1)
wtot2 ~ dbin(p2, n2)
p1 <- r[1]*SN[1] + (1-r[1])*(1-SP[1])
p2 <- r[2]*SN[2] + (1-r[2])*(1-SP[2])

# validation study
# note order (x,w)=(0,0),(0,1),(1,0),(1,1)
xwtot1[1:4] ~ dmulti(q1[1:4], m1)
xwtot2[1:4] ~ dmulti(q2[1:4], m2)
q1[1] <- (1-r[1])*SP[1]
q1[2] <- (1-r[1])*(1-SP[1])
q1[3] <- r[1]*(1-SN[1])
q1[4] <- r[1]*SN[1]
q2[1] <- (1-r[2])*SP[2]
q2[2] <- (1-r[2])*(1-SP[2])
q2[3] <- r[2]*(1-SN[2])
q2[4] <- r[2]*SN[2]

### work on logit scale
logit(SN[1]) <- lgtSN[1]; logit(SN[2]) <- lgtSN[2]
logit(SP[1]) <- lgtSP[1]; logit(SP[2]) <- lgtSP[2]
logit(r[1]) <- lgtr[1]; logit(r[2]) <- lgtr[2]

### bivariate normal priors
lgtSN[1:2] ~ dmnorm(muvecSN[1:2], prcSN[1:2,1:2])
lgtSP[1:2] ~ dmnorm(muvecSP[1:2], prcSP[1:2,1:2])
lgtr[1:2] ~ dmnorm(muvecr[1:2], prcr[1:2,1:2])

### both components have same mean, variance
### precision matrix determined by var and correlation
muvecSN[1] <- muSN; muvecSN[2] <- muSN
muvecSP[1] <- muSP; muvecSP[2] <- muSP
muvecr[1] <- mur; muvecr[2] <- mur

prcSN[1,1] <- 1/(varSN*(1-pow(rhoSN,2)))
prcSN[2,2] <- 1/(varSN*(1-pow(rhoSN,2)))
prcSN[1,2] <- (-(rhoSN))/(varSN*(1-pow(rhoSN,2)))
prcSN[2,1] <- (-(rhoSN))/(varSN*(1-pow(rhoSN,2)))

prcSP[1,1] <- 1/(varSP*(1-pow(rhoSP,2)))
prcSP[2,2] <- 1/(varSP*(1-pow(rhoSP,2)))
prcSP[1,2] <- (-(rhoSP))/(varSP*(1-pow(rhoSP,2)))
prcSP[2,1] <- (-(rhoSP))/(varSP*(1-pow(rhoSP,2)))

prcr[1,1] <- 1/(varr*(1-pow(rhor,2)))
prcr[2,2] <- 1/(varr*(1-pow(rhor,2)))
prcr[1,2] <- (-(rhor))/(varr*(1-pow(rhor,2)))
prcr[2,1] <- (-(rhor))/(varr*(1-pow(rhor,2)))

}