## Initial setup
beetles <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/701/data/beetles.txt")
with(beetles, plot(Dose, y/n, pch=19, ylab="Fraction killed"))

## Model
Data <- c(beetles, N=nrow(beetles))
model <- function() {
  ## Likelihood
  for (i in 1:N) {
    y[i] ~ dbinom(theta[i], n[i])
    logit(theta[i]) <- beta[1] + beta[2]*Dose[i]
  }
  for (j in 1:2) {
    beta[j] ~ dnorm(0, 0.0001)
  }
}
fit <- jags(model=model, data=Data, param="beta", n.chains=1, n.iter=11000, n.burn=1000, n.thin=1, DIC=FALSE)
fit

## Posterior summary
attach.jags(fit)
post.sum <- function(x) c(mean(x), quantile(x, c(.05, .95)))
post.sum(beta[,2]*0.05)
post.sum(exp(beta[,2]*0.05))

## Plot
xx <- seq(min(beetles$Dose), max(beetles$Dose), len=99)
f <- function(beta, x) {
  eta <- beta[1] + beta[2]*x
  exp(eta)/(1+exp(eta))
}
Y <- apply(beta, 1, f, xx)
plot(0, 0, xlim=c(1.65, 1.9), ylim=c(0, 1), xlab="Dose", ylab="P(Kill)", type="n", las=1)
yhat <- apply(Y, 1, mean)
lower <- apply(Y, 1, quantile, prob=.05)
upper <- apply(Y, 1, quantile, prob=.95)
drawBand(xx, lower, upper)
lines(xx, yhat)
with(beetles, points(Dose, y/n, pch=19))

## Lethal doses
f <- function(beta, prob) {
  eta <- log(prob/(1-prob)) 
  (eta - beta[1])/beta[2]
}
LD <- apply(beta, 1, f, 0.5)
post.sum(LD)
LD <- apply(beta, 1, f, 0.99)
post.sum(LD)