Chapter 6 Exercises
Beetles
Solutions

1. Run the code from Example 6.5.1 to analyse the "beetles" data with
centered doses of carbon disulphide. Create two sets of initial values so that a two-chain analysis can be performed. When do the chains converge? And what is the degree of autocorrelation in the chains? Select "Correlations" from the "Inference" menu and create a scatter-plot of alpha vs beta to examine their posterior correlation [Hint: don't use too many samples as this plot requires a substantial amount of memory.]

model {
for (i in 1:8) {
y[i] ~ dbin(p[i], n[i])
logit(p[i]) <- alpha + beta*(x[i] - mean(x[]))
phat[i] <- y[i]/n[i]
yhat[i] <- n[i]*p[i]
}
alpha ~ dnorm(0, 0.0001)
beta ~ dnorm(0, 0.0001)
}

list(x = c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839),
n = c(59, 60, 62, 56, 63, 59, 62, 60),
y = c(6, 13, 18, 28, 52, 53, 61, 60))
   
list(alpha = 50, beta = -50)
list(alpha = -50, beta = 50)

With centred doses convergence occurs almost instantly:

[exercises-ch6-beetles-solutions0][exercises-ch6-beetles-solutions1]

[exercises-ch6-beetles-solutions2][exercises-ch6-beetles-solutions3]
Autocorrelation in the chains is negligible (note that the first bar shows perfect correlation at "lag-0", because this is the extent to which each sample is correlated with itself -- this will always occur):

[exercises-ch6-beetles-solutions4][exercises-ch6-beetles-solutions5]
There is a very small amount of positive posterior correlation between alpha and beta:

[exercises-ch6-beetles-solutions6]

2. Adapt the code to run the analysis with un-centred doses of carbon disulphide. What is the impact on convergence, autocorrelation, and posterior correlation between alpha and beta?

model {
for (i in 1:8) {
y[i] ~ dbin(p[i], n[i])
logit(p[i]) <- alpha + beta*x[i]
phat[i] <- y[i]/n[i]
yhat[i] <- n[i]*p[i]
}
alpha ~ dnorm(0, 0.0001)
beta ~ dnorm(0, 0.0001)
}

list(x = c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839),
n = c(59, 60, 62, 56, 63, 59, 62, 60),
y = c(6, 13, 18, 28, 52, 53, 61, 60))
   
list(alpha = 50, beta = -50)
list(alpha = -50, beta = 50)

Convergence is much slower with un-centred doses. The history and BGR plots both suggest that at least ~5000 iterations are required for convergence. It would be prudent to run the chains for considerably longer: (a) to make sure that convergence has occurred; and (b) to ensure accurate inferences -- the summary statistics for the latter 5000 iterations show MC errors of around 8% of the corresponding posterior standard deviations.

[exercises-ch6-beetles-solutions7][exercises-ch6-beetles-solutions8]

[exercises-ch6-beetles-solutions9][exercises-ch6-beetles-solutions10]
   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   alpha   -61.88   4.416   0.368   -69.62   -62.56   -53.58   5001   10000
   beta   34.93   2.483   0.2069   30.27   35.31   39.27   5001   10000
   
There is huge autocorrelation in the chains, even at lag-50 (the correlation between samples 50 iterations apart):

[exercises-ch6-beetles-solutions11][exercises-ch6-beetles-solutions12]

There is almost perfect (negative) posterior correlation between alpha and beta. This explains the above difficulties. Because the Gibbs sampler can only move perpendicular to the axes (in the plot below, say), conditional on any given value for alpha, there is a very limited range of plausible values for beta, and vice versa. Hence, exploration of the joint distribution is very slow as only a small distance can be covered in each Gibbs iteration. Centering covariates reduces the posterior correlation between intercept and gradient parameters, as illustrated in the previous alpha-beta scatter-plot.

[exercises-ch6-beetles-solutions13]


3. Try fitting probit versions of the centred model, using both the probit link and the phi function [note that these can be sensitive to initial values] . Is the model fit visually different compared to the logistic model?

model {
for (i in 1:8) {
y[i] ~ dbin(p[i], n[i])
# probit (p[i]) <- alpha + beta*(x[i] - mean(x[]))
p[i] <- phi(alpha + beta*(x[i] - mean(x[])))
phat[i] <- y[i]/n[i]
yhat[i] <- n[i]*p[i]
}
alpha ~ dnorm(0, 0.0001)
beta ~ dnorm(0, 0.0001)
}

list(x = c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839),
n = c(59, 60, 62, 56, 63, 59, 62, 60),
y = c(6, 13, 18, 28, 52, 53, 61, 60))
   
list(alpha = 5, beta = -5)
list(alpha = -5, beta = 5)

Both versions of the probit model struggle with the initial values previously specified. Using the less extreme values shown above, both models run without a problem. (Note that the phi formulation is generally slower but is typicaly more stable, leading to fewer numerical issues.) Visually, the model fit is almost identical with that from the logistic model. The parameter estimates are quite different but this is to be expected since logistic and probit scales are different.

[exercises-ch6-beetles-solutions14][exercises-ch6-beetles-solutions15]

Probit model:
   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   alpha   0.4473   0.07737   6.024E-4   0.2959   0.4467   0.6016   1   20000
   beta   19.83   1.49   0.01213   17.01   19.81   22.81   1   20000
Logistic model:
   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   alpha   0.7501   0.1382   0.001125   0.4824   0.7491   1.029   1   20000
   beta   34.59   2.929   0.02442   29.14   34.53   40.58   1   20000