dobson:83 analyses binary dose-response data published by bliss:35, in which the numbers of beetles killed after 5 hour exposure to carbon disulphide at N=8 different concentrations are recorded:
We assume that the observed number of deaths 
at each concentration 
is binomial with sample size 
and true rate 
. Plausible models for include the logistic, probit and extreme value (complimentary log-log) models, as follows
The corresponding graph is shown in Figure 17.
Figure 17:
Graphical model for beetles example
The fit of each model may be assessed by calculating the deviance D as follows (see also the seeds example). The log-likelihood for an observation arising from a binomial model with denominator and success probability is
The saturated log likelihood for the binomial model is
The deviance is calculated by computing a node 
llike.sat 
llike 
) within BUGS. This will yield a posterior distribution for D, the minimum value of which corresponds to the classical deviance obtained using maximum likelihood estimation. This may be compared with a 
distribution to assess model fit.
model beetles;
const
N = 8; # number of doses
var
r[N],p[N],x[N],n[N],alpha,alpha.star,beta,r.hat[N],llike[N],llike.sat[N],D;
data x, n, r in "beetles.dat";
inits in "beetles.in";
{
for (i in 1:N) {
r[i] ~ dbin(p[i], n[i]);
logit(p[i]) <- alpha.star + beta*(x[i]-mean(x[]));
# alternative link functions:
# probit(p[i]) <- alpha.star + beta*(x[i]-mean(x[]));
# cloglog(p[i]) <- alpha.star + beta*(x[i]-mean(x[]));
# log likelihood for sample i & saturated log-likelihood:
llike[i] <- r[i]*log(p[i]) + (n[i]-r[i])*log(1-p[i]);
llike.sat[i] <- r[i]*log(r[i]/n[i]) + (n[i]-r[i])*log(1-r[i]/n[i]);
r.hat[i] <- p[i]*n[i]; # fitted values
}
alpha.star ~ dnorm(0.0, 1.0E-3);
beta ~ dnorm(0.0, 1.0E-3);
alpha <- alpha.star - beta*mean(x[]);
D <- 2 * (sum(llike.sat[]) - sum(llike[]));
}
Note that we have standardized each dose about the mean: this gives approximately uncorrelated regression coefficients, and greatly improves convergence. Figure 18 shows the sample traces (ploted using CODA) for alpha and beta after a 10000 iteration BUGS run. The first run (chain:beetles1) used the centered parameterization, whilst the second run (chain:beetles2) used the uncentered parameterization. The chains from the latter run have still not converged, and exhibit very high autocorrelations (see Figure 19), whilst those from the former run converge almost immediately, and exhibit rapid mixing rather than a slow, `snaking' trace; this is reflected by the much lower autocorrelation.
Figure 18:
Traces for the models with centered covariates (beetles1) and uncentered covariates (beetles2)
Figure 19:
Plot of the within-chain autocorrelations for the models with centered covariates (beetles1) and uncentered covariates (beetles2)
Analysis
1000 iterations took 2 seconds after a 500 iteration burn-in. The BUGS posterior means and standard errors for the regression coefficients and fitted values, plus the minimum deviance for each model are given below, as are Dobson's maximum likelihood estimates (MLE). Note that the equivalent of MLE's may be obtained from the BUGS output by taking the values of alpha, beta and r.hat sampled at the iteration yielding the minimum value for the deviance. These are also given in the table below.
Comparison of the minimum deviances indicates that the extreme value model fits the data considerably better than do the logistic or probit models. This appears to be due to a smaller discrepancy between observed ( ) and fitted ( 
) values at the lower concentrations for the extreme value model.