model {
for (i in 1:n) {
lbw[i] ~ dbern(p[i])
logit(p[i]) <- beta[1] + beta[2]*THM[i] +
beta[3]*male[i] + beta[4]*dep[i] +
beta[5]*smk[i] + beta[6]*eth[i]
}
for (k in 1:6) {
beta[k] ~ dnorm(0, 0.0001)
}
for (k in 2:6) {
OR[k] <- exp(beta[k])
}
# multivariate probit covariate imputation model
for (i in 1:n) {
Z[i,1:2] ~ dmnorm(mu[i,1:2],
Omega[1:2,1:2])I(lo[i,1:2],up[i,1:2])
mu[i,1] <- delta[1,1] + delta[2,1]*THM[i] +
delta[3,1]*male[i] + delta[4,1]*dep[i] +
delta[5,1]*area.smk[i] +
delta[6,1]*area.eth[i]
mu[i,2] <- delta[1,2] + delta[2,2]*THM[i] +
delta[3,2]*male[i] + delta[4,2]*dep[i] +
delta[5,2]*area.smk[i] +
delta[6,2]*area.eth[i]
}
for (i in 1:Nmis) { # Data file is ordered so subjects
# 1,...,Nmis have missing values
smk[i] <- step(Z[i,1]) # thresholded value of Z[i,1]
eth[i] <- step(Z[i,2]) # thresholded value of Z[i,2]
}
Sigma[1,1] <- 1
Sigma[2,2] <- 1
Sigma[1,2] <- corr
Sigma[2,1] <- corr
corr ~ dunif(-1, 1)
Omega[1:2, 1:2] <- inverse(Sigma[,])
for (k in 1:6) {
delta[k,1] ~ dnorm(0, 0.0001)
delta[k,2] ~ dnorm(0, 0.0001)
}

dummy <- smk.ful[1] + eth.ful[1]

}

Inits:
list(beta = c(1.5, 0.18, -0.27, 0.29, 0.69, 1.38), corr=0.5,
delta=structure(.Data = c(0,0,0,0,0,0, 0,0,0,0,0,0), .Dim=c(6,2)))

Data...

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   OR[2]   1.178   0.121   0.002129   0.956   1.172   1.425   1001   20000
   OR[3]   0.8095   0.08177   0.001518   0.6607   0.8052   0.9821   1001   20000
   OR[4]   1.007   0.101   0.001919   0.8243   1.001   1.221   1001   20000
   OR[5]   2.915   0.5202   0.01434   2.033   2.867   4.057   1001   20000
   OR[6]   4.116   0.7181   0.01833   2.897   4.053   5.716   1001   20000
   corr   -0.3039   0.05583   0.002219   -0.4096   -0.3046   -0.1904   1001   20000