Past exam of the mathematics course of the University of Cambridge / 2008 / iii / Paper 44 / 3

Write the new measurement as $Y_{i1}={\theta }_i+{\epsilon }_{i1}$, with independent ${\theta }_i\sim N(\mu ,{\tau }^2)$ and ${\epsilon }_{i1}\sim N(0,{\sigma }_i^2)$. Adding independent normal distributions gives the prior predictive distribution

Let $z_{0.999}={\Phi }^{-1}(0.999)$ be the standard normal distribution quantile. Standardizing the prediction in part (a) givesThis is a one-sided predictive limit; its upper-tail probability is $0.001$ under the specified model.

By normal-normal conjugacy, the posterior distribution isThe posterior mean is a precision-weighted average of the population mean and the observed personal measurement.

Integrate the next independent measurement error over the posterior distribution of ${\theta }_i$. The posterior predictive distribution isBoth uncertainty about the personal mean and variation in a new measurement contribute.

Taking the upper $0.999$ quantile of the posterior predictive distribution givesThe same predictive false-alarm probability is retained after updating with the first measurement.

After $k$ measurements the normal-normal conjugacy formulas areThe upper monitoring limit is $m_{ik}+z_{0.999}\sqrt{v_{ik}+{\sigma }_i^2}$. Its center moves toward the athlete's own mean and the population's influence diminishes. Its predictive width decreases toward $z_{0.999}{\sigma }_i$, since uncertainty about the mean vanishes but the next measurement's error remains. The absolute upper limit need not decrease monotonically: its center can move upward or downward with the observations. This is sequential normal predictive monitoring.

The hierarchical Bayesian model separates variability between individuals from variability between repeated measurements on an individual. In WinBUGS, dnorm takes a precision parameter, so inverse variance is the appropriate second argument.

Line (1) gives $\log {\sigma }_j^2\sim N(\varphi ,{\psi }^2)$, hence a positive log-normal distribution for each within-individual variance. Estimating the shared parameters supplies partial pooling of these variances without forcing them equal. Line (2) assigns a proper bounded flat prior to the between-individual standard deviation $\tau$, and line (3) does the same for the standard deviation $\psi$ on the log-variance scale. The broad bounded priors on $\mu$ and $\varphi$ are also proper.

This is a reasonable structure for heterogeneous repeated measurements, but flat does not mean uninformative on all transformed scales. The bounds must be wide enough for plausible values, while such enormous ranges can put substantial mass on implausible predictions. Use prior predictive checks and sensitivity analysis, and check the Markov chain Monte Carlo convergence diagnostics. Serial dependence or systematic trends would require an extension of the independent-measurement likelihood.

Add the new individual's latent mean and log variance to the same population hierarchy. Enter only that individual's first two measurements as observations; use the historical data and these two values to fit the model. The third measurement must be held out, or the predictive check reuses the value it is meant to test. For example, augment the unchanged population hyperpriors with this WinBUGS structure:

new.mean ~ dnorm(mu, invtau2)
new.logvar ~ dnorm(phi, invpsi2)
new.precision <- exp(-new.logvar)
for (k in 1:2) {
    new.Y[k] ~ dnorm(new.mean, new.precision)
}
future.Y ~ dnorm(new.mean, new.precision)
upper.indicator <- step(future.Y - held.out.third)

future.Y is unobserved. Its simulated values give the posterior predictive distribution, including uncertainty in the population parameters and personal variance. A central 99.9% prediction interval uses its $0.0005$ and $0.9995$ empirical quantiles. Flag the third value if it is outside that interval. Equivalently, estimate its predictive tail probability by the mean of upper.indicator and flag if that mean is below $0.0005$ or above $0.9995$. This is two-sided, unlike the earlier one-sided upper limit.

Conditionally on ${\lambda }_i$, the proposed precision gives ${\theta }_i\mid {\lambda }_i\sim N(\mu ,4{\tau }^2/{\lambda }_i)$. Thus, with independent $Z\sim N(0,1)$ and ${\lambda }_i\sim {\chi }_4^2$,By the Student t-distribution construction, $\boxed{({\theta }_i-\mu )/\tau \sim t_4}$. Here $\tau$ is a scale rather than the standard deviation: the prior variance of ${\theta }_i$ is $2{\tau }^2$.

The Student t-distribution has heavier tails than the normal distribution, so unusual personal means need not force the entire population's spread upward or be shrunk as strongly toward its center. The latent ${\lambda }_i$ allows a large local prior variance for an outlying individual. This normal scale mixture representation of Student t is a useful robust version of the hierarchical Bayesian model. It still assumes a symmetric, unimodal population; a genuinely multimodal population may need a different model.

Fit both hierarchies to the historical data and compare their posterior predictive checks. Use discrepancies sensitive to between-individual tails, such as the maximum standardized personal mean, tail counts and the spread of individual sample means. Simulate new individual means as well as repeated observations; conditioning only on fitted personal means would conceal a poor population model.

Leave whole individuals out for cross-validation, rather than just measurements, when the aim is prediction for a new individual. Compare held-out predictive log densities or interval calibration, with uncertainty in the comparison. Check both fits' Markov chain Monte Carlo convergence diagnostics and sensitivity to scale priors. Better tail prediction and adequate calibration would support the heavier-tailed model; the data are required to decide, so no empirical winner can be asserted here.