Past exam of the mathematics course of the University of Cambridge / 2018 / iii / Paper 218 (Statistical Learning in Practice)

The violated assumption is independence of random variables across observations: repeated measurements on the same rat form clustered data. A persistent rat-specific weight difference induces positive within-rat covariance, even after adjustment for time and diet. Treating all 176 measurements as independent is pseudoreplication.

For the claimed lack of identifiability in the second fit, Rat must be a categorical variable, as intended by the rat-effect interpretation. Each rat receives only one diet, so this is confounding of nested fixed factors. If $z_r$ is the rat-indicator column, the diet-$d$ indicator is ${\sum }_{r:D_r=d}{z}_r$. Thus a diet effect can be increased by $c_d$ and every corresponding rat effect decreased by $c_d$, leaving all fitted values unchanged. With an intercept, two diet columns, fifteen rat contrasts and time, there are 19 columns but only 17 independent columns: the intercept and rat contrasts already span the diet columns, while time varies within rats. This gives rank-deficient ordinary least squares.

The first model ignores within-rat dependence; the second cannot separate unrestricted fixed diet and rat effects. If Rat were instead encoded as a single numeric predictor, the asserted rank deficiency would not follow merely from nesting.

Write $r=1,\dots ,16$ for rats and $j$ for repeated observations. With diet 1 as the reference level in a regression factor, the random-intercept linear mixed model iswith the two families of normal distributions independent and the covariates treated as fixed. The reported estimates areThe corresponding estimated standard deviations are $36.577$ and $8.176$ grams. Individual realized random intercepts are not parameters in this marginal model, and their predictions are not supplied by this summary. For distinct measurements on the same rat, $\mathrm{Cov}(Y_{rj},Y_{rk})={\tau }^2$, while the marginal variance is ${\tau }^2+{\sigma }^2$; different rats are independent.

The estimated mean weight rises by $0.58568$ grams per day, adjusting for diet and the rat's random intercept. This common slope is imposed across all diets: there is no time-by-diet interaction term, so these estimates do not compare diet-specific growth slopes. The intercept refers to time zero for diet 1 and a zero rat effect.

Let $Z_{ir}=1_{\{\text{observation }i\text{ belongs to rat }r\}}$. After integrating out the random intercepts,The fixed-effect design matrix has four independent columns, so the given $A$ has dimensions $176\times 172$, with $A^{T}A=I_{172}$ and $A^{T}X=0$. Consequently the error contrasts $U=A^{T}Y$ have a multivariate normal distribution with mean zero and covariance matrix $W=A^{T}VA$.

The restricted likelihood from orthogonal error contrasts is their density. Equivalently, the estimates maximizeThis is restricted maximum likelihood, not maximization over individual rat effects. The unknown $\beta$ disappears because $A^{T}X=0$. Replacing $A$ by $AQ$ for an orthogonal $Q$ rotates $U$ and conjugates $W$, preserving the determinant and quadratic form; hence the objective is independent of the chosen orthonormal basis. The software's integrated REML log-likelihood can differ by an additive constant depending only on $X$, which has no effect on these variance estimates.

There are two defects. The mixed fit supplies a restricted maximum likelihood, whereas the ordinary fit supplies an ordinary log-likelihood, so their difference is not a likelihood-ratio test statistic. Also, $H_0:{\tau }^2=0$ lies on the boundary of ${\tau }^2\ge 0$; the usual Wilks theorem does not give a ${\chi }_1^2$ null law.

A valid approach first refits the mixed model by ordinary maximum likelihood estimation, setting REML to false, and fits the same fixed effects under $H_0$. For $V={\tau }^{2}Z{Z}^T+{\sigma }^{2}I$, maximizeover $\beta ,{\sigma }^2>0,{\tau }^2\ge 0$, and separately over $\beta ,{\sigma }^2>0$ with ${\tau }^2=0$. Set $T=2({\hat{\ell }}_1-{\hat{\ell }}_0)\ge 0$.

For finite-sample calibration, use the location-scale invariant simulation test for a Gaussian variance component. With the actual $X,Z$ fixed, simulate $B$ independent vectors $Y^{(b)}\sim N(0,I_{176})$, refit both models by ordinary maximum likelihood estimation to each, and calculate $T_b$ in exactly the same way. Under $H_0$, $Y=X\beta +\sigma z$. Translating by a vector in the column space of $X$ and multiplying by a positive scalar preserves the statistic: both maximized log-likelihoods acquire the same scale constant. Thus this simulation has the correct null law without knowing $\beta$ or $\sigma$. A conservative Monte Carlo p-value isIn the usual regular limit with increasing independent groups, a variance-component likelihood-ratio test at a boundary instead uses $\frac{1}{2}{\delta }_0+\frac{1}{2}{\chi }_1^2$: for $T>0$, its approximate tail probability is $\frac{1}{2}\mathrm{Pr}({\chi }_1^2\ge T)$, and at $T=0$ the nonrandomized p-value is one. This approximation is not an exact guarantee for sixteen rats. Alternatively one can compare consistently defined restricted likelihoods with the same $X$ and simulate their restricted-ratio null distribution, as in the RLRsim documentation.

In the terminology used for positive-definite kernels, positive definiteness means positive semidefiniteness of every finite kernel matrix: $k$ is symmetric and, for every $m$, points $x_1,\dots ,x_m$ and real $c_1,\dots ,c_m$,The linear kernel is symmetric, andThus it is a positive-semidefinite kernel, with the identity as its feature map. It need not be strictly positive definite: a zero vector or linearly dependent points can make this quadratic form zero for nonzero $c$. This distinction is essential when interpreting the terminology.

Fix $C>0$ and use the unnormalized sum convention for the penalty. The linear soft-margin support vector machine solvesEliminating the slack variables of a support vector machine gives $\frac{1}{2}\Vert w{\Vert }^2+C{\sum }_i\max \{0,1-y_{i}f(x_i)\}$, where $f(x)=w^{T}x+b$ and the maximum is the hinge loss. The classifier uses the sign of $f$ and its support-vector-machine decision boundary is $f=0$.

For clarity about dual support vectors and margin degeneracy, let ${\alpha }_i$ be the multipliers of the margin constraints. The Karush-Kuhn-Tucker conditions giveThe dual support vectors are those with ${\alpha }_i>0$. A coefficient strictly between $0$ and $C$ puts a point exactly on its margin boundary; a point with positive slack has ${\alpha }_i=C$. A zero coefficient implies zero slack and signed margin at least one. Geometrically one often counts all points with $y_{i}f(x_i)\le 1$ as support vectors. At a degeneracy, equality can occur with ${\alpha }_i=0$, so the two descriptions need not coincide. Either convention supports the bound in the following part when stated consistently.

Choose a candidate grid of $C>0$. For each candidate and each $i$, fit the soft-margin support vector machine to all observations except $i$, keeping exactly the same $C$ in the unnormalized sum objective. Let $f_{-i,C}$ be its fitted score. Use a fixed rule for assigning a label when the score is zero. The binary misclassification estimate from Leave-one-out cross-validation isChoose a minimizing candidate, break tuning ties by a predetermined rule, and refit on all $n$ observations. Any estimated preprocessing must also be fitted inside each training fold to avoid data leakage. If a normalized average-loss parameterization is used instead, translate its parameter so that the same $C$ is retained when the sample size decreases.

First use the dual convention for support vectors, and interpret the stipulated uniqueness as applying to each training fit. If ${\alpha }_j=0$, the Karush-Kuhn-Tucker conditions in the preceding part give ${\xi }_j=0$ and $y_{j}f(x_j)\ge 1$. Delete the $j$th primal variable and its zero dual coefficient. All remaining primal constraints, stationarity equations and complementary-slackness equations still hold with the same $w,b$. These conditions certify a global minimizer of the reduced convex optimization problem, since the objective is convex and the Slater condition holds by choosing positive large slacks. Uniqueness therefore gives exactly the same fitted score and support-vector-machine decision boundary. The omitted point is correctly classified.

For the geometric convention, a non-support point has strict margin $y_{j}f(x_j)>1$. One can obtain the same conclusion using just uniqueness of the full fit: if the reduced problem had another optimizer, the convex segment from the original optimizer toward it would also consist of reduced optimizers. A sufficiently short segment retains strict margin at $j$, hence has zero omitted hinge loss and would give another full optimizer, a contradiction.

Only a support vector can therefore cause a held-out error. Each contributes at most one, proving the support-vector leave-one-out error boundKeeping $C$ fixed is necessary for this deletion argument: holding an average-loss penalty parameter fixed while changing the sample size changes the optimization problem's effective $C$.

Conditional on the covariates, the binary responses are independent with Bernoulli distributions:Thus $p_i=(1+e^{-{\eta }_i}{)}^{-1}$, the logistic regression model. The maximum-likelihood estimates areUsing the reported age standard error, an approximate 95% Wald confidence interval isThis interval concerns the age coefficient in log odds per year, with income and loan amount held fixed.

The sigmoid function is strictly increasing. Therefore the risky region is $\hat{\eta }(x)\ge \log \{p/(1-p)\}$. Its decision boundary is the hyperplanePoints on or above this score threshold are assigned to the risky class. For $p=1/2$, the right side is zero.

The feedforward neural network has three input coordinates, two hidden sigmoid function units and two output softmax function units. Every input is connected to each hidden unit, and each hidden unit to each output unit; each hidden and output unit also has a bias. A layer diagram isAlgebraically, with $x_i=({\mathrm{income}}_i,{\mathrm{age}}_i,{\mathrm{amount}}_i{)}^T$,with responses conditionally independent. The fourteen displayed weight and bias parameters comprise six input weights, two hidden biases, four output weights and two output biases. The one-hot encoding supplies targets $t_{i0}=1-Y_i$, $t_{i1}=Y_i$, and the categorical cross-entropy loss is the negative of $\ell ={\sum }_i{\sum }_c{t}_{ic}\log p_{ic}$. There is softmax non-identifiability under a common shift of both output logits, but the predicted probabilities are well defined.

For one observation, let ${\ell }_i={\sum }_{c=0}^1{t}_{ic}\log p_{ic}$ and define ${\delta }_{ic}=t_{ic}-p_{ic}$. Differentiating the softmax function and its log-likelihood gives $\partial {\ell }_i/\partial o_{ic}={\delta }_{ic}$. The sigmoid-softmax network gradients areThese equations are backpropagation: use the chain rule through the output affine map, then ${\sigma }^{\prime }(u)=\sigma (u)(1-\sigma (u))$ through the hidden units, then the input affine map.

For a learning rate $\eta >0$, select one observation and update every parameter by $\theta \leftarrow \theta +\eta \nabla {\ell }_i(\theta )$, computing all entries of the gradient from the old parameters before applying the update. A uniformly selected observation gives an unbiased gradient of the average log-likelihood; multiplying by 3000 would instead estimate the gradient of the sum. Equivalently perform stochastic gradient descent on the negative log-likelihood. Batch size one makes one update per observation, and five epochs make five passes through the data, typically with a fresh random ordering on each pass.

Let $T_i=1$ denote Progabide and $T_i=0$ placebo; write $a_i$ for age and $b_i$ for baseline count. The fitted Poisson regression specifies independentIts log-likelihood isThe coefficient table in the original PDF, whose TeX transcription is truncated, gives the maximum-likelihood estimatesFor the treatment interaction contrast in a Poisson regression, the baseline slope of the log mean is $0.0186824$ on placebo and $0.0195506$ on Progabide. At fixed age and baseline $b$, the fitted drug-to-placebo mean ratio is $\exp (-0.3592068+0.0008682b)$. Increasing baseline by one multiplies that ratio by $e^{0.0008682}\approx 1.0008686$. Thus the interaction term describes how the treatment ratio changes with baseline, rather than a uniform additive change in seizure counts. The displayed call makes clear that the source's summary(model1) line is an object-name typo for the epilepsy fit.

The null model has only an intercept:For $\bar{Y}>0$, maximizing its log-likelihood gives $\hat{\mu }=\bar{Y}$ and ${\hat{\beta }}_0=\log \bar{Y}$. The saturated Poisson regression gives each observation its own fitted mean ${\hat{\mu }}_i=Y_i$, allowing zero as a limiting fitted mean. By the definition of Poisson deviance,Since ${\sum }_i(Y_i-\bar{Y})=0$,Use the continuous convention $0\log (0/\bar{Y})=0$. If all counts are zero, both the saturated and null likelihood maxima are attained as limits at zero mean and the deviance is zero.

Under a correctly specified Poisson regression with adequate expected counts, the residual Poisson deviance is approximately chi-squared distribution with $n-p=107-5=102$ residual degrees of freedom. The observed $D=319.86$ givesThus the ordinary Poisson model gives a very poor fit under this goodness-of-fit approximation. This is consistent with overdispersion, although a wrong mean model or dependence can also cause lack of fit; quasi-Poisson variance inflation alone does not establish that the mean model is adequate.

For Quasi-Poisson regression, $\mathrm{Var}(Y_i)=\varphi {\mu }_i$. The Pearson dispersion estimator isIt uses squared Pearson residuals, not the residual deviance divided by 102. The fitted coefficients are unchanged, and their estimated covariance matrix is multiplied by $\hat{\varphi }$. With $\hat{\varphi }=3.24$, standard errors increase by $1.8$, so the main treatment coefficient hasThe main Progabide coefficient is not significant at 5% after this dispersion adjustment. R uses a Student's t-distribution when dispersion is estimated, as described in its GLM summary documentation; the large-sample normal approximation also gives the same decision.

Because an interaction term is present, this tests the treatment contrast at baseline zero. An overall absence of treatment influence requires the joint null ${\beta }_T={\beta }_{TB}=0$. The printed marginal standard errors do not give their covariance, so they are insufficient to calculate that joint Wald test, or the treatment contrast test at an arbitrary nonzero baseline.

Use the stationary autoregressive process of order oneMultiplication by $X_{t-1}-\mu$ and expectation, using independence of the innovation from the past, gives the autocovariance equation $\gamma (1)=\varphi \gamma (0)$. With $\bar{X}=n^{-1}{\sum }_t{X}_t$, define the sample autocovariances using the same denominator:Substitution into that equation gives the Yule–Walker estimator for an autoregressive process of order oneThis is exactly the lag-one sample autocorrelation function under the divisor-$n$ convention, provided the sample has positive variance. Using different divisors for different lags would change the estimator.

The parameters are $\mu \in \mathbb{R}$, $|\varphi |<1$ and innovation variance ${\sigma }^2>0$. Stationarity gives $X_1\sim N(\mu ,{\sigma }^2/(1-{\varphi }^2))$, and the Markov property gives $X_t\mid X_{t-1}\sim N(\mu +\varphi (X_{t-1}-\mu ),{\sigma }^2)$. Multiplying the initial density and the transition densities yields the stationary Gaussian AR1 likelihood, with $n=31$:The initial-observation factor depends on all three parameters and must be retained for the full stationary likelihood function. Omitting it gives a conditional likelihood given $X_1$, a different objective.

By the Markov property, all factors in the stationary likelihood function that involve $X_2$ are the transitions from $X_1$ to $X_2$ and from $X_2$ to $X_3$. Hence conditioning on later observations adds no information once the two neighbors are known. Writing $u=X_2-\mu$, $u_1=x_1-\mu$ and $u_3=x_3-\mu$, its conditional density is proportional toThe quadratic coefficient of $u$ is $1+{\varphi }^2$ and the linear coefficient is $-2\varphi (u_1+u_3)$. Completing the square therefore gives the Gaussian AR1 bridgeHere ${\sigma }^2$ is the innovation variance, not the marginal stationary variance.

Apply EM for a missing observation in a Gaussian AR1 process, starting from ${\theta }^{old}=({\mu }^{old},{\varphi }^{old},({\sigma }^2{)}^{old})$ in the stationary parameter space. The Gaussian AR1 bridge gives the E-step quantitiesThus $\mathbb{E}[X_2\mid x_{obs}]=m$, $\mathbb{E}[X_2^2\mid x_{obs}]=m^2+v$, and $\mathbb{E}[X_1{X}_2\mid x_{obs}]=x_{1}m$, $\mathbb{E}[X_2{X}_3\mid x_{obs}]=m{x}_3$.

Let $x_2^{\ast }=m$ and $x_t^{\ast }=x_t$ otherwise. Keeping $m,v$ fixed throughout the M-step, defineThe last term comes from the missing observation's conditional variance in the two adjacent innovation squares. The expected complete-data log-likelihood isFor fixed $\varphi$, differentiating the quadratic in $\mu$ givesAlso the maximizing variance is ${\sigma }^2=R(\mu ,\varphi )/31$. The M-step is therefore a one-dimensional profile likelihood maximization:Repeat the E- and M-steps until the observed-data likelihood and parameters stabilize. Exact M-steps make that likelihood nondecreasing; different initializations help detect different local optima, and convergence alone is not a guarantee of a global maximum. At an interior local maximum, impute the missing temperature by its conditional mean $m$ evaluated at the fitted parameters, retaining $v$ as its conditional uncertainty. Replacing $X_2$ by $m$ and discarding $v$ would not implement the expectation-maximization algorithm.

For the residual-versus-fitted plot, compute fitted values ${\hat{y}}_i={\hat{\beta }}_0+{\hat{\beta }}_1{x}_{i1}+{\hat{\beta }}_2{x}_{i2}$ and regression residuals $e_i=y_i-{\hat{y}}_i$; the plotted point is $({\hat{y}}_i,e_i)$.

For the normal Q-Q plot, let $H=X(X^{T}X{)}^{-1}{X}^T$ be the hat matrix for the full three-column design matrix and let $s^2=\mathrm{RSS}/(n-3)$. Form the standardized regression residualsand order them as $r_{(1)}\le \cdots \le r_{(n)}$. R plots $r_{(j)}$ against ${\Phi }^{-1}(p_j)$, where $p_j$ are its plotting positions. For $n>10$, $p_j=(j-1/2)/n$; $\Phi$ is the standard normal distribution function. This standardization is specified in R's diagnostic-plot documentation.

The original PDF shows a residual cloud around zero without compelling systematic curvature or a clear funnel, and a Q-Q pattern approximately following the reference line. Observations 81 and 83 are conspicuous positive tail points, with a smaller lower-tail departure. The plots show no decisive violation of linearity, constant variance or normality, but the extreme residuals merit inspection. They do not test independence across observations and cannot establish that all model assumptions hold.

The two predictors have correlation coefficient $\rho =0.9813846$, so there is severe multicollinearity. Their centered slope Gram matrix isThus either partial slope has variance ${\sigma }^2/[n(1-{\rho }^2)]$. Relative to a univariate regression with the same error variance, the two-predictor variance inflation factor and standard-error factor areThe small eigenvalue $n(1-\rho )$ makes the contrast between the two slope effects particularly poorly estimated. The univariate regressions test marginal association with the response; the bivariate tests ask whether one predictor adds an effect while holding its nearly identical companion fixed. Their much larger standard errors therefore make the partial tests insignificant without contradicting the marginal significance. The two fits also estimate slightly different residual variances, explaining why the reported standard-error ratio is not exactly 5.21.

The command performs a nested-model F-test for $H_0:{\beta }_1=0$ against $H_1:{\beta }_1\ne 0$, retaining the intercept and $x_2$ coefficient in both models. Write $RS{S}_2$ and $RS{S}_3$ for the residual sum of squares in the reduced and full fits. Under $H_0$ and the independent normal-error assumptions,By Cochran's theorem, the numerator and denominator are independent scaled chi-squares with one and $n-3$ degrees of freedom. Equivalently, for this single-coefficient restriction, $F=t_{{\beta }_1}^2$, so its p-value equals that of the two-sided coefficient Student's t-test in model3.

Using the printed $t$ value gives $F=0.800^2=0.640$ and $p\approx 0.43$. The printed estimate and standard error instead giveThis discrepancy is present in the PDF and is too large to be explained by rounding at the shown precision; the underlying ANOVA output is omitted, so an exact unique statistic cannot be recovered. The univariate correlation and $t$ value imply $n\approx 2+2.575^2(1-0.2517297^2)/0.2517297^2\approx 100.006$, consistent with $n=100$ after rounding. With denominator degrees of freedom 97, the two versions give p-values approximately $0.4257$ and $0.4298$, respectively. Either way the partial effect is not significant at 5%.

Let $\lambda$ denote the effective penalty in the Gaussian squared-error objective. With the loss normalization used in the glmnet documentation, alpha zero gives ridge regression:The intercept is fitted separately and unpenalized. The constant column supplied by model.matrix is not an additional penalized slope; glmnet excludes zero-variance predictor columns. Centering and the stated predictor norms make the two slope columns already standardized under the divisor-$n$ convention.

Let $X_s$ consist of those two columns, $z=y-\bar{y}1$, and $G=X_s^T{X}_s$. Minimizing over the intercept gives ${\hat{\beta }}_0=\bar{y}$, and the closed-form ridge regression estimator under this normalization isThus the right-hand limit in the plotted log-$\lambda$ coordinate is $\lambda \to \infty$, givingThe left-hand limit is $\lambda \to 0$, giving the unpenalized full-model slopesThe plotted curves are these two slopes. If the objective is instead written as unnormalized RSS plus a quadratic penalty, its penalty coefficient is $n\lambda$, not $\lambda$.

There is also a software response-scale convention to distinguish from this mathematical parameterization. The Gaussian fitting engine divides a supplied penalty by the response standard deviation $s_y=\sqrt{n^{-1}{\sum }_i(y_i-\bar{y}{)}^2}$ while fitting a standardized response, then rescales the coefficients. For pure ridge with these already standardized predictors, this corresponds to effective $\lambda ={\lambda }_{\mathrm{API}}/s_y$ in the displayed objective. Such a fixed positive rescaling of the horizontal parameter for the observed dataset does not change either asymptote. The fixed-penalty variance formulas below refer to $\lambda$ in the displayed mathematical objective.

Condition on the fixed design matrix and a prespecified $\lambda$. Put $M=G+n\lambda I_2$. The covariance and bias of a ridge regression estimator follow from its being a linear transformation of the response:Here centering the response introduces no correction in this expression because $X_s^{T}1=0$. For the equal-norm two-predictor design, either diagonal entry simplifies toHowever,The variances quantify sampling variability, but alone do not give valid confidence intervals for the unshrunk slopes. Centering a normal interval at the ridge estimate ignores its bias, which depends on the unknown slopes and need not be small relative to its standard error. Intervals require an appropriate bias correction, a bound on that bias, or another justified inferential procedure. A data-selected penalty introduces additional selection dependence not covered by the fixed-$\lambda$ formula. Similarly, holding a software penalty fixed while converting it using an estimated response scale does not hold the effective penalty fixed over repeated samples; the displayed covariance is not an unconditional variance formula for that nonlinear procedure.

One suitable criterion among the three unpenalized normal linear models is the Akaike information criterion. With $d$ estimated mean coefficients and the variance also estimated, minimizing AIC is equivalent to minimizingThe values of $d$ are $2,2,3$ for model1, model2 and model3. This compares fit and complexity rather than choosing solely by individual coefficient significance. The Bayesian information criterion or a common-fold predictive comparison would also be reasonable, depending on the selection objective.

For the ridge regression tuning parameter, use K-fold cross-validation with squared prediction loss:Use the same folds for every candidate, estimate any centering or scaling within each training fold, choose a minimizing $\lambda$, and refit on all observations. This estimates mean squared prediction error while avoiding data leakage. The printed summaries do not supply all residual sums of squares or validation losses, so they do not determine which model or penalty wins.