Past exam of the mathematics course of the University of Cambridge / 2016 / iii / Paper 206 / 6

Write $t_{ij}$ for stirring rate of observation $i$ in furnace $j$. The independent random-intercept and random-slope model isErrors and random effects are independent. The separate R terms (1 | furnace) and (0 + stir | furnace) impose independent random intercepts and random slopes; (1 + stir | furnace) would instead estimate their covariance as well. The estimates areThe corresponding estimated standard deviations are $6.9813$, $0.1147$, and $6.0249$; they are not additional parameters.

For a furnace with predictor vector $t_j$, the marginal formulation, obtained by integrating the Gaussian random effects, isindependently across furnaces. In particular,In stacked notation $Y=X\beta +Zb+\epsilon$, $b\sim N(0,D)$, giving the Gaussian linear mixed model marginal law $N(X\beta ,V)$ with $V=ZD{Z}^T+{\sigma }^{2}I$.

The residual-versus-fitted plot on page 17 has residuals of both signs over the fitted range, with no convincing smooth curvature or clear fan shape. Its quantile-quantile plot is roughly straight in the central region, with some tail departures and a few large negative residuals. There is no decisive visible violation, but investigate those observations. These plots mainly address conditional observation errors; they do not validate the distribution of furnace random effects or independence within furnaces. With only three furnaces, normality and the variance of the random effects are especially difficult to assess.

For the Gaussian linear mixed model, temporarily regard $\beta ,D,R$ as known, with $R={\sigma }^{2}I$ here and $V=ZD{Z}^T+R$. The joint multivariate normal distribution isBy the conditional multivariate normal distribution formula,A nonsingular Gaussian density has its mode at its mean. Thus the conditional mode of Gaussian random effects equals the conditional expectation. Replacing the unknown parameters by fitted values givesFor positive-definite $D$ and $R$, an equivalent derivation minimizes the negative conditional log density,Differentiation yieldsThe inverse-covariance penalty shrinks furnace intercepts and slopes toward the population effects. For furnace $j$, set $Z_j=(1,t_j)$ and $D_j=\mathrm{diag}({\tau }_0^2,{\tau }_1^2)$; the same formula provides its two-component conditional mode from that furnace's observations, conditional on the fitted common parameters. If a variance component is zero, the covariance form remains meaningful and its corresponding random effect is zero, whereas the formula with $D^{-1}$ needs a limiting interpretation. With known covariance parameters these are best linear unbiased predictors; plug-in variance estimates give their empirical version.

The restricted maximum likelihood procedure estimates covariance parameters from error contrasts that remove the unknown fixed effects. Let $X$ have rank $p$ and choose an $n\times (n-p)$ full-rank matrix $A$ with $A^{T}X=0$. Then $A^{T}Y\sim N(0,A^{T}V(\vartheta )A)$ has a likelihood involving the variance parameters $\vartheta$ but not $\beta$. Maximizing that likelihood is restricted maximum likelihood. An equivalent restricted likelihood from orthogonal error contrasts, up to terms independent of $\vartheta$, iswhereAfter estimating $\vartheta$, use generalized least squares for $\hat{\beta }=(X^T{\hat{V}}^{-1}X{)}^{-1}{X}^T{\hat{V}}^{-1}Y$ and the conditional modes from part (b). In an independent-error normal linear model, this covariance procedure reduces to estimating ${\sigma }^2$ by $\mathrm{RSS}/(n-p)$ instead of $\mathrm{RSS}/n$.

The two models being compared have different fixed-effect spaces, so their error contrasts, dimensions, and restricted likelihoods differ. Their REML values are not likelihoods of the same transformed data and must not be used directly for this fixed-effect likelihood-ratio test. Refit both with ordinary maximum likelihood, integrating out the random effects but retaining the original response vector, and maximize each model's marginal likelihood.

The null hypothesis is ${\beta }_1=0$, leaving the random slope variance free; the alternative permits a nonzero population stirring effect. The supplied test has $2({\ell }_1-{\ell }_0)=7.9588$, one extra fixed coefficient, and approximate ${\chi }_1^2$ p-value $0.004786$. Retain the fixed stirring-rate effect at 5%. A mean-zero random slope cannot substitute for an estimated nonzero population slope. The asymptotic calibration assumes regular nuisance covariance parameters and should be interpreted cautiously with only three furnaces; simulation can provide a better finite-sample calibration.

The comparison removes the random slope and tests $H_0:{\tau }_1^2=0$ against $H_1:{\tau }_1^2>0$, while retaining the same fixed intercept and slope. A variance cannot be negative, so zero is a boundary point of the parameter space. The regular Wilks theorem assumptions behind an ordinary ${\chi }_1^2$ likelihood-ratio test fail. Consequently the default chi-squared p-value from anova is not reliably calibrated for this variance-component likelihood-ratio test at a boundary. This is different from the fixed-slope test in part (c). The reference to strength_model1 means the earlier object strength_model; no separate model 1 object was defined.

A useful alternative is a parametric bootstrap. Fit the null random-intercept model, then simulate data on the actual predictor values and furnace grouping by generating null intercept effects and observation errors from its fitted Gaussian distribution. Refit null and alternative models to each simulated dataset using the same likelihood method and compute the same likelihood-ratio statistic. If $T_{\mathrm{obs}}$ is the observed statistic, estimate its upper-tail probability byThis incorporates the zero-variance boundary and the actual small group count. It is a fitted-null simulation approximation because nuisance parameters are estimated, not a universally exact test. Since both models have the same fixed-effect space, a similarly calibrated restricted-likelihood statistic is also possible. Under suitable asymptotics for a single identifiable added variance component the familiar limiting law is $\frac{1}{2}{\chi }_0^2+\frac{1}{2}{\chi }_1^2$, but using the bootstrap avoids relying on that approximation with only three groups. There is no supplied statistic or calibrated bootstrap result here from which to decide whether the slope variance is required.

The generalized additive mixed model fitswhere $u_j$ denotes the furnace random intercept and $f$ is a centered penalized cubic regression spline. Its wiggly components have a mixed-model representation; there is no explicit furnace-specific random slope in this formula.

The plotted smooth gives no evidence that a nonparametric fixed effect is needed. It is essentially a straight rising line, and s(stir,1) indicates effective degrees of freedom close to one. The wider bands at extreme stirring rates reflect uncertainty, not detected nonlinearity. A linear fixed stirring effect is adequate on this evidence.

If the three supplied Akaike information criterion values are on a common valid likelihood basis, choose strength_model3, the fixed linear slope plus random-intercept model: its AIC $206.6408$ is smallest, versus $208.3243$ for the additional random slope and $210.3425$ for the smooth. The differences are $1.6835$ and $3.7017$, respectively, so especially the random-slope comparison represents modest support rather than decisive evidence. Retaining the population slope is consistent with part (c); this choice removes only its random variation.

There is a likelihood-basis caveat in the literal printed commands. lmer fits by REML by default, whereas the Gaussian gamm call defaults to ML; direct AIC on these objects need not compare the same type of likelihood. Moreover the earlier ML ANOVA prints $210.22$ for strength_model, whereas this table prints $208.3243$, so these numbers should not silently be treated as one identical ML fit. To obtain a defensible common comparison, fit all candidates to the same observations by ML:The supplied ranking favours the simpler random-intercept model; the displayed original calls alone do not certify that ranking as a consistent ML comparison. This distinction is particularly relevant when the smooth changes the fixed-effect space. The Gaussian gamm ML default and mixed-model representation are documented at stat.ethz.ch/R-manual/R-devel/library/mgcv/html/gamm.html.