## Exploratory
puromycin <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/701/data/puromycin.txt")
xyplot(Rate~Conc|State, puromycin, layout=c(2,1), pch=19)

## Model
require(R2jags); invisible(runif(1))
Data <- c(puromycin, n=nrow(puromycin))
Data$State <- as.numeric(Data$State)
model <- function() {
  for (i in 1:n) {
    Rate[i] ~ dnorm(mu[i], tau)
    mu[i] <- (Vm[State[i]]*Conc[i])/(Km[State[i]]+Conc[i])
  }
  for (j in 1:2) {
    Vm[j] ~ dunif(0,500)
    Km[j] ~ dunif(0,2)
  }
  tau ~ dgamma(0.001, 0.001)
}
fit <- jags(model=model, data=Data, param=c("Km", "Vm"), n.chains=1, n.iter=11000, n.burn=1000, n.thin=1, DIC=FALSE)
fit

## Summarizing the posterior
post.sum <- function(x) c(mean(x), quantile(x, c(.05, .95)))
attach.jags(fit)
dnplot(Km, labs=levels(puromycin$State), xlab="K")
post.sum(Km[,1] - Km[,2])
dnplot(Vm, labs=levels(puromycin$State), xlab="Vm")
post.sum(Vm[,1] - Vm[,2])

E <- Vm/Km
dnplot(E, labs=levels(puromycin$State), xlab="Efficiency", ylab="")
post.sum(E[,1] - E[,2])

## Plotting the model fit
xx <- seq(min(puromycin$Conc), max(puromycin$Conc), len=99)
f <- function(x, s) (x[1]*s)/(x[2]+s)
par(mfrow=c(1,2))
for (j in 1:2) {
  A <- apply(cbind(Vm[,j], Km[,j]), 1, f, s=xx)
  plot(0, 0, xlim=c(0, 1.1), ylim=c(0, 250), xlab="Concentration", ylab="Rate", type="n", main=levels(puromycin$State)[j], las=1)
  yhat <- apply(A, 1, mean)
  lower <- apply(A, 1, quantile, prob=.05)
  upper <- apply(A, 1, quantile, prob=.95)
  drawBand(xx, lower, upper)
  lines(xx, yhat)
  with(subset(puromycin, State==levels(State)[j]), points(Conc, Rate, pch=19, cex=0.7))
}