Past exam of the mathematics course of the University of Cambridge / 2009 / iii / Paper 36 / 3

At several doses the three replicate plates differ much more than the square-root-of-mean scale expected from a Poisson distribution. For example the group at dose $100$ has sample mean about $42.67$ and sample variance about $274.33$, rather than variance near $42.67$. Its counts span a range of $33$, compared with a Poisson standard deviation about $6.53$. The dose-$33$ and dose-$1000$ groups likewise have variances above their means. Thus extra plate-to-plate variation is plausible. With only three replicates per dose this is evidence suggesting overdispersion, rather than a precise determination of it; not every group has excess sample variance.

Let ${\eta }_i=\alpha +\beta \log (x_i+10)+\gamma x_i$. The Poisson-lognormal random-effect model assigns each plate a random intensity ${\mu }_{ij}=e^{{\eta }_i+{\lambda }_{ij}}$. Conditional on that intensity its Poisson distribution has mean and variance both ${\mu }_{ij}$, but marginalizing over the normal distribution of the random effect givesThe law of total variance yieldsThis exceeds the mean whenever $\tau >0$. The independent plate effects supply additional heterogeneity without inducing dependence between plates conditional on the global parameters. Also $e^{{\eta }_i}$ is the median random intensity, not its marginal mean; the factor $e^{{\tau }^2/2}$ matters for population predictions.

Centre and scale both predictors, particularly the raw dose whose scale is much larger than that of $\log (x_i+10)$. For example use $z_{1i}=[\log (x_i+10)-c_1]/s_1$, $z_{2i}=(x_i-c_2)/s_2$ and the linear predictor ${\eta }_i=a_0+b_1{z}_{1i}+b_2{z}_{2i}$. This reduces the intercept-slope correlations and numerical scale disparities; further orthogonalizing correlated columns or updating coefficients jointly can improve Markov chain Monte Carlo mixing. The exact old prior can be retained through the corresponding coefficient transformation; assigning new independent priors to transformed coefficients would instead change the prior.

For weakly informed plate effects or small $\tau$, use the non-centred Gaussian random-effect parameterisation ${\lambda }_{ij}=\tau z_{ij}$ with $z_{ij}\sim N(0,1)$. This often reduces dependence between the scale and local effects; it is not guaranteed to outperform the centred version when plate effects are strongly informed. Sensibly calibrated proper priors also avoid the enormous initial log intensities permitted by the stated coefficient ranges. Use dispersed chains and Markov chain Monte Carlo convergence diagnostics to assess stationarity, agreement between chains and effective sample size of a Markov chain. Reparameterization can improve mixing; simply retaining more iterations does not remove poor mixing.

The regression coefficients have independent proper uniform distributions on very wide finite intervals. These are flat in the coefficients, not flat in fitted rates, and can be highly informative on the rate scale. The scale $\tau$ is uniform on $(0,100)$, and WinBUGS receives the corresponding precision $1/{\tau }^2$ for each normal distribution random effect. Writing $v={\tau }^2$, the induced prior density isIt is proper and integrable at zero.

The usual scale-invariant variance prior $p(v)\propto 1/v$ would instead cause an improper posterior from a log-uniform random-effect scale prior. After integrating out the plate effects, the likelihood function tends as $v\downarrow 0$ to the ordinary positive Poisson likelihood. On any compact interior set of coefficient values this limit is positive and bounded away from zero. Hence the posterior normalizing integral contains a divergent factor ${\int }_0^{\epsilon }dv/v$. There is no posterior distribution under that prior, even if formal full conditional distributions appear usable. Thus a proper prior with integrable mass near zero is needed here. The reciprocal-scale rule is the standard pure-scale Jeffreys prior; it should not be confused with the model-specific information prior for an additive variance component.

Delete the plate random effects from the linear predictor and remove their stochastic nodes and the unused scale/precision nodes. The observation layer then remains a Poisson distribution withFor example the essential WinBUGS observation block can be written asRetain the coefficient priors. This is the $\tau =0$ sampling model implemented directly; attempting to keep 1/(tau*tau) while fixing tau to zero would divide by zero.

Under the common observation-level Bayesian deviance convention, DIC is $\overline{D}+p_D$, with the effective parameter count in DIC defined by $p_D=\overline{D}-D(\overline{\psi })$. The no-random-effect model's effective count is near three, corresponding to its three fitted regression coefficients. With plate random effects, the effective count rises to $13.6$: the local effects add flexibility, but their partial pooling means they do not count as eighteen fully independent unrestricted fitted parameters. The effective count need not be an integer or equal the raw parameter count; the scale parameter controls that shrinkage.

Adding the plate effects reduces the average Bayesian deviance by $28.6$ at an effective-complexity cost of $10.7$. Therefore its DIC is lower byThe reported criterion favours the overdispersed model. This is a fit-complexity comparison, not a posterior probability or Bayes factor. In a hierarchical Bayesian model, the interpretation depends on whether deviance is conditional on local effects or has marginalized them out; these figures use the supplied observation-level convention.

For each plate use its conditional fitted intensity, including the fitted random effect. At a posterior draw $s$ form the Pearson residualA simpler plot uses ${\hat{\mu }}_{ij}$ in this expression, but posterior draws also display uncertainty in the fitted values. Look for residuals centred around zero, unexpectedly large tails or isolated outliers, and patterns with dose, fitted intensity or plate grouping. For known conditional parameters the Poisson distribution gives mean zero and variance one for this residual; approximate normality is reasonable only at adequate count sizes. Fitting the same data and estimating latent effects changes its calibration, so compare the diagnostics with residuals from conditional posterior predictive checks, rather than treating them as independent exact standard normals.

This checks the plate-level Poisson sampling layer given the effects and can reveal remaining unmodelled variation. It does not by itself establish the population dose-response curve: fitted plate effects can absorb departures from that curve.

Use conditional versus marginal posterior predictive checks to replicate new plates, rather than recycling the fitted plate effects. For each joint posterior draw of $(\alpha ,\beta ,\gamma ,\tau )$, draw fresh independent ${\lambda }_{ij}^{\mathrm{rep}}\sim N(0,{\tau }^2)$ and thenUse the same doses and three-plate design, giving the posterior predictive distribution for repeated experiments. Compare observed dose means and their trend, the high-dose downturn, within-dose spreads, and outliers with the replicated curves and intervals. For example compare a discrepancywith $T(y^{\mathrm{rep}};\psi )$ draw by draw. The resulting posterior predictive tail fraction or graphical rank assesses whether the observed trend is unusual under this model. Fresh plate effects are needed to assess the population dose-response assumption. Reusing the fitted effects would mainly check the conditional sampling layer and could hide systematic departures from the proposed curve.