## Setup
Data <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/701/data/birthweight.txt")
attach(Data)
n <- nrow(Data)

## Ignorable missingness model
jData <- list("lbw", "THM", "Male", "Dep", "Smk", "Eth", "n")
model <- function() {
  for (i in 1:n) {
    lbw[i] ~ dbern(p[i])
    logit(p[i]) <- b[1] + b[2]*THM[i] + b[3]*Male[i] + b[4]*Dep[i] + b[5]*Smk[i] + b[6]*Eth[i]
  }
  for (k in 1:6) {
    b[k] ~ dnorm(0, 0.0001)
  }

  # Imputation
  for (i in 1:n) {
    Smk[i] ~ dbern(theta[1])
    Eth[i] ~ dbern(theta[2])
  }
  for (j in 1:2) {
    theta[j] ~ dunif(0,1)
  }
  Smk14 <- Smk[14]
  Eth14 <- Eth[14]
  Smk1476 <- Smk[1476]
  Eth1476 <- Eth[1476]
}
require(R2jags); invisible(runif(1))
fit <- jags(model=model, init=init, param=c("b", "Smk14", "Eth14", "Smk1476", "Eth1476"), data=jData, n.iter=5000, n.thin=1)
attach.jags(fit)

## Convergence
max(gelman.diag(as.mcmc(fit))$psrf)
min(effectiveSize(as.mcmc(fit)))

## Posteriors
psm(exp(b)[,-1], labels=c("THM", "Male", "Dep", "Smk", "Eth"))
B <- psm(exp(b))[-1,]
rownames(B) <- c("THM", "Male", "Deprived", "Smoking", "Ethnicity")
CIplot(B, xlab="Odds ratio")
abline(v=1, col="gray")

## Imputation
data.frame(Data[c(14, 1476),], SmkImp=c(mean(Smk14), mean(Smk1476)), EthImp=c(mean(Eth14), mean(Eth1476)))
mean(Eth, na.rm=TRUE)
mean(Smk, na.rm=TRUE)

## Comparisons with simple logistic regression
eci <- function(fit, j) {c(coef(fit)[j], confint(fit, j))}
exp(eci(glm(lbw~THM, family=binomial), 2))
exp(eci(glm(lbw~THM+Male+Dep+Smk+Eth, family=binomial), 2))

## Simple imputation data
IData <- Data
IData$Smk[is.na(Smk)] <- rbinom(sum(is.na(Smk)), size=1, prob=mean(Smk, na.rm=TRUE))
IData$Eth[is.na(Eth)] <- rbinom(sum(is.na(Eth)), size=1, prob=mean(Eth, na.rm=TRUE))
ifit <- glm(lbw~THM+Male+Dep+Smk+Eth, family=binomial, data=IData)
exp(eci(ifit, 5))
exp(eci(ifit, 6))