## Setup
Data <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/701/data/roots.txt")
Data$Photoperiod <- as.factor(Data$Photoperiod)
require(R2jags); invisible(runif(1))

## Exploratory
x <- 0:max(Data$Roots)
y <- sapply(x, function(a) sum(Data$Roots==a))
b <- barplot(y/sum(y), names=x, border="white", xlab="Roots", ylab="Percent", las=1)
lines(b, dpois(x, mean(Data$Roots)), col=pal(2)[2], lwd=3, type="o", pch=19)

## ZIP Model
n <- nrow(Data)
jData <- c(Data, n=n)
jData$Photoperiod <- as.numeric(jData$Photoperiod)
model <- function() {
  ## Likelihood
  for (i in 1:n) {
    Roots[i] ~ dpois(m[i])
    m[i] <- group[i] * mu[Photoperiod[i]]
    group[i] ~ dbern(pi[Photoperiod[i]])
  }
  for (j in 1:2) {
    mu[j] ~ dgamma(0.5, 0.0001)
    pi[j] ~ dunif(0, 1)
  }
}
init <- function() {
  list("group"=rep(1, n))
}
fit <- jags(model=model, data=jData, init=init, param=c("mu", "pi"), n.chains=1, n.iter=11000, n.burn=1000, n.thin=1, DIC=FALSE)
fit

## Posterior
attach.jags(fit)
dnplot(mu, labs=levels(Data$Photoperiod), xlab=expression(mu))
dnplot(pi, labs=levels(Data$Photoperiod), xlab=expression(pi))

post.sum <- function(x) c(mean(x), quantile(x, c(.05, .95)))
post.sum(mu[,1]-mu[,2])
post.sum(pi[,1]/pi[,2])

## Plot of distribution
x <- 0:max(Data$Roots)
y1 <- sapply(x, function(a) sum(Data$Roots[Data$Photoperiod=="8"]==a))
y2 <- sapply(x, function(a) sum(Data$Roots[Data$Photoperiod=="16"]==a))

par(mfrow=c(1,2), oma=c(0,0,2,0), mar=c(5,4,1,1))
Y1 <- t(outer(x, mu[,1], dpois))
Y <- Y1*pi[,1]
Y[,1] <- Y[,1]+1-pi[,1]
plot(x-0.1, y1/sum(y1), type="h", lwd=5, col=pal(2)[1], ylim=c(0, 0.25), xlab="Roots", ylab="Probability")
lines(x+0.1, apply(Y, 2, mean), type="h", lwd=5, col=pal(2)[2])
mtext("8 hr")

Y1 <- t(outer(x, mu[,2], dpois))
Y <- Y1*pi[,2]
Y[,1] <- Y[,1]+1-pi[,2]
plot(x-0.1, y2/sum(y2), type="h", lwd=5, col=pal(2)[1], ylim=c(0, 0.5), xlab="Roots", ylab="Probability")
lines(x+0.1, apply(Y, 2, mean), type="h", lwd=5, col=pal(2)[2])
mtext("16 hr")
toplegend(c("Observed", "Fit"), ncol=2, lwd=5, col=pal(2))