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

The transformation $t=(x-\overline{x})/s_x$ centers day at zero and measures it in standard-deviation units. It reduces numerical scale differences and makes the intercept-like parameter refer to the middle of the observed time range. In nonlinear regression this can reduce posterior parameter dependence and improve Markov chain Monte Carlo mixing. It also makes priors on the growth coefficient easier to interpret; it does not change the curve family.

Write $A_i={\varphi }_{i1}$, $b_i={\varphi }_{i2}>0$ and $k_i=-{\varphi }_{i3}={\theta }_{i3}$. ThenFor $k_i>0$, $A_i$ is the limiting circumference and $1/(1+b_i)$ is the fraction of that asymptote attained at the centered day $x=\overline{x}$. Under the code's transformations,The latter expression is the logistic function. The curve also satisfies the logistic growth equation $d{m}_i/dx=(k_i/s_x)m_i(1-m_i/A_i)$. Its inflection point is at half the asymptote. The original prior permits $k_i\le 0$, in which case it does not describe increasing growth.

The priors are proper but difficult to justify as practically weak information. Uniform bounds $[-100,100]$ on log asymptotes cover fantastically small and large circumferences. The independent prior on the transformed location makes the centered fraction spend much of its mass near zero or one, while the slope prior permits decreasing growth. Independence also ignores likely relationships among curve parameters.

The $\mathrm{Gamma}(0.001,0.001)$ prior on error precision is proper, but its near-zero shape creates a highly skewed prior with consequential mass near very small precision. It is not automatically harmless. Use biologically scaled priors for asymptotes and timing, a positive growth-rate prior, and a defensible prior on the error standard deviation such as a suitably scaled half-normal distribution. Check simulated trajectories with prior predictive checks and assess prior sensitivity.

Model B replaces independent tree parameters by multivariate normal random effects. It permits correlated asymptote, timing and slope variation and supplies partial pooling, a useful improvement if the five trees can be regarded as exchangeable.

The Wishart distribution prior is on the precision matrix $\Omega$, not the covariance matrix. In the WinBUGS convention, dimension three and degrees of freedom three give a proper but very dispersed prior. The induced inverse-Wishart distribution has no finite mean at those degrees of freedom. Thus this is not a neutral choice of covariance prior, and estimating six covariance components from only five trees is weakly informed. The broad mean priors and possible negative growth slopes also remain. Check prior-predictive curves and sensitivity, and consider separate, interpretable scale and correlation priors.

The resulting posterior means of the second and third parameters are quite similar across trees, motivating a simpler shared-shape model; their uncertainty means this similarity is evidence to explore, rather than a proof of exact equality.

A common second parameter means the trees reach the same fraction of their asymptotic sizes at the centered day; a common third parameter means they have the same standardized growth rate. The trees may differ mainly in final size while sharing a growth shape. Similar environment and the overlapping posterior estimates make this biologically plausible. Equality is a stronger modelling assumption than partial pooling, and should be checked predictively.

Keep a tree-specific log asymptote $a_i$ with a normal population distribution, but use one shared location $b$ and one shared positive growth coefficient $k$. For example, with fixed prior constants chosen on meaningful scales:

mean.a ~ dnorm(m.a, precision.a.prior)
sd.a ~ dnorm(0, precision.sd.a) T(0,)
precision.a <- 1 / pow(sd.a, 2)
b ~ dnorm(m.b, precision.b.prior)
k ~ dnorm(m.k, precision.k.prior) T(0,)
sigma ~ dnorm(0, precision.sigma.prior) T(0,)
precision.y <- 1 / pow(sigma, 2)
for (i in 1:5) {
    a[i] ~ dnorm(mean.a, precision.a)
    for (j in 1:7) {
        fitted[i,j] <- exp(a[i]) / (1 + exp(-b - k*t[j]))
        Y[i,j] ~ dnorm(fitted[i,j], precision.y)
    }
}

Here the transformed days t[j] are supplied as data. This is Model C: only the log asymptote is a random effect, while location and growth rate are shared. The prior constants are fixed inputs, and the truncations give positive standard deviations and an increasing growth curve. One can retain the original error-precision prior for a comparison matching the supplied fits; the code above illustrates a more interpretable positive-scale alternative.

Smaller mean Bayesian deviance indicates better average likelihood fit. Model B has the smallest $\overline{D}$, improving on A by $3.9$ and on C by $3.1$. Model C has the smallest effective parameter count in DIC, reflecting its shared timing and growth coefficients. Model B's hierarchical covariance parameters and residual flexibility contribute to a larger effective count than C.

Adding fit and complexity gives the reported deviance information criterion: A has $260.1$, B $258.1$, and C $255.6$. Thus DIC favors Model C, whose simpler shape more than compensates for its slightly worse mean fit than B. The differences are modest, and the criteria must use the same observation likelihood and deviance convention. These numbers are not posterior model probabilities; predictive assessment remains useful.

Model A has fifteen separate curve parameters plus one error-scale parameter. With such broad priors, one might have expected more of those sixteen parameters to be effectively estimated, rather than $p_D=12$, and might also have expected its greater unpooled flexibility to give better fit than the displayed $\overline{D}=248.1$. Model B instead has a smaller mean deviance despite sharing information between trees.

This is a diagnostic concern, not a mathematical impossibility. Effective parameter count in DIC depends on posterior curvature and parameterization, while mean posterior deviance is not minimized deviance. A weakly identified nonlinear asymptote–slope trade-off, skewness, multiple posterior regions or slow Markov chain Monte Carlo mixing could explain the unexpected comparison. For A, the implied deviance at the posterior mean is $\overline{D}-p_D=236.1$; for B it is $230.3$. Inspect traces, effective sample sizes, independent chains and fitted curves, and extend sampling where needed. Fifty thousand iterations alone do not establish convergence.

Assess the observation-error assumption with standardized residual plots and quantile-quantile plots, preferably comparing the same discrepancies with posterior predictive replicated data. Inspect tail behavior, asymmetry, variance versus fitted size and serial patterns within a tree; fitted residuals are dependent and do not behave exactly like independent errors.

This does not test normality of the tree random effects. For that, simulate new tree-specific log asymptotes from the population hierarchy and compare between-tree discrepancies or use leave-one-tree-out predictions. There are only five trees, so distributional checks on their population are intrinsically weak. Conditioning on the five fitted effects can hide misspecification of their normal population law.