Past exam of the mathematics course of the University of Cambridge / 2013 / iii / Paper 30 (Applied Statistics)

Use the normal linear modelHere $Y$ is the response random vector, $X$ is the known design matrix, $\beta \in {\mathbb{R}}^p$ contains the unknown regression coefficients, and ${\sigma }^2>0$ is the common error variance. Conditional on $X$, the errors have normal distributions and are independent random variables. Require $\mathrm{rank}X=p\le n$ for identifiability of $\beta$ and invertibility of $X^{T}X$. Usually $n>p$ is needed to estimate the error variance from the regression residuals. An intercept, when included, is represented by a column of ones in $X$.

Full column matrix rank makes $X^{T}X$ a positive-definite matrix, and its matrix inverse is symmetric. Thus the hat matrix satisfiesAlso $PX=X$ and the image of $P$ is contained in the column space of $X$. These identities show that $P$ is the orthogonal projection matrix onto the column space of $X$.

The ordinary least squares estimator is $\hat{\beta }=(X^{T}X{)}^{-1}{X}^{T}Y$. Consequently the fitted values and regression residuals areAn affine transformation of a multivariate normal distribution is again multivariate normal, possibly with a singular covariance matrix. For a random vector with covariance matrix $\Sigma$, its transformed covariance matrix is $A\Sigma A^T$. Since $PX=X$, $P^2=P=P^T$ and $(I-P{)}^2=I-P$, these results giveBoth multivariate normal distributions are supported on their respective projected subspaces. In particular, ${\hat{Y}}_i$ has variance ${\sigma }^2{P}_{ii}$ and $e_i$ has variance ${\sigma }^2(1-P_{ii})$; the regression residuals need not be mutually independent.

Because both vectors are linear transformations of the same multivariate normal response, $(\hat{Y},e)$ is jointly multivariate normal. Its cross-covariance matrix isZero cross-covariance implies independence for jointly multivariate normal vectors, including singular ones. Therefore the fitted values and the entire vector of regression residuals are independent. The fitted-residual orthogonality identity gives the zero covariance; the normal distribution assumption is what upgrades it to independence.

The simple linear regression assumes a straight conditional mean, $E(Y_i\mid x_i)={\beta }_0+{\beta }_1{x}_i$, so that errors are centred at zero throughout the predictor range. The displayed regression residuals are predominantly positive at both ends and negative in the middle. This is a residual curvature diagnostic: the fitted straight line misses a curved conditional mean. The main concern is the shape of the mean function. The plot alone does not establish failure of the normal distribution assumption or a particular error variance model.

A natural next fit is quadratic regression, which remains a normal linear model in its unknown regression coefficients:The new design matrix has rows $(1,x_i,x_i^2)$ and must have matrix rank three; three distinct predictor values suffice. The curved regression residual pattern suggests trying a positive quadratic term, but its sign and adequacy should be checked after fitting. Inspect the new regression residuals to see whether the systematic curvature has disappeared.

Write $m_i=\mu +{\alpha }_i$. Apart from a constant, the log-likelihood isFor each group,Thus maximum likelihood estimation sets ${\hat{m}}_i={\overline{Y}}_i$. The corner-point constraint ${\alpha }_1=0$ identifies $m_1=\mu$, givingThe baseline mean is the first group mean, rather than the grand mean, because of the chosen identifiability constraint.

Group means are independent random variables with normal distributionsThe resulting univariate sampling distributions areThe same variance calculation holds for every $i\ge 2$, since the difference involves two independent group means. Therefore the standard errors satisfyReplacing $\sigma$ by a common estimated error standard deviation preserves this ratio. Although the individual group means are independent, the contrasts ${\hat{\alpha }}_i$ share the baseline mean and have covariance ${\sigma }^2/J$ for distinct $i\ge 2$.

Index chocolate by $a\in \{A,B,C,D\}$, day by $d$ in the three observed categories, and replicate by $r\in \{1,2\}$. The additive two-factor normal linear model isHere $Y_{adr}$ is the board count, $\mu$ is the mean for the reference chocolate and reference day, and ${\alpha }_a,{\gamma }_d$ are chocolate and day fixed effects. With corner-point constraints, set ${\alpha }_A=0$ and ${\gamma }_{d_0}=0$, where $d_0$ is the first level in the day factor. The printed coefficient-free output does not determine that factor ordering; the model is unchanged by a different reference category. There is no chocolate–day interaction term in this fit. The six free mean statistical parameters consist of one baseline, three chocolate contrasts and two day contrasts; the common error variance supplies a further statistical parameter.

Test $H_0:{\gamma }_d=0$ for both nonreference days, against at least one nonzero day fixed effect, conditional on chocolate. There are two added regression coefficients and $24-6=18$ residual statistical degrees of freedom. The nested-model F-test statistic isUnder $H_0$ and the normal linear model assumptions, $F\sim F_{2,18}$. Its 5% upper critical value is $3.555$, so do not reject the absence of day effects; the p-value is approximately $0.456$. Thus the missing row has two statistical degrees of freedom, mean square $4.875$, and the F-test value above. These observations do not provide evidence that including day improves the chocolate-adjusted mean model. They do not prove that every possible day effect or chocolate–day interaction term is absent.

The analysis of variance provides evidence of a chocolate effect both with day included ($F_{3,18}=3.9665$, p-value $0.02473$) and with day omitted ($F_{3,20}=4.039$, p-value $0.02133$). Because every chocolate–day cell has the same replication, balanced factorial orthogonality separates the two main effects; the chocolate row is meaningful despite being entered first.

The chocolate-only fitted group means are approximately $14.17,11.67,9.33,11.33$ boards for A, B, C, D respectively. Thus A has the largest fitted lecturing speed and C the smallest. Relative to A, the fitted differences are $-2.50,-4.83,-2.83$ boards. The printed individual Student t-tests give strong evidence for the A–C contrast (p-value $0.00245$); B and D versus A have p-values $0.08835$ and $0.05582$, respectively. Those latter contrasts are not significant at 5%, and the output does not test all other pairwise comparisons. Simultaneous claims would require accounting for multiple hypothesis testing.

There is little evidence of a day effect after adjusting for chocolate. A chocolate-only normal linear model is therefore a reasonable simpler summary, with residual standard deviation about $2.42$ boards and explained variation $R^2\simeq 0.377$. This is an association under the additive statistical model; a causal claim would additionally require an appropriate assignment of chocolate and checks of the regression residuals and possible interaction terms.

An exponential dispersion family has density or mass function, relative to a fixed measure,Here $\theta$ is the natural parameter, $b$ is the cumulant function, and $\varphi$ is the dispersion parameter. Assume the natural parameter lies in the interior of its domain and derivatives can pass through the normalizing integral. Differentiating normalization once and twice givesThe variance function is $V(\mu )=b^{\prime \prime }((b^{\prime }{)}^{-1}(\mu ))$, where the mean-to-natural-parameter inverse exists. The dispersion parameter may be fixed, as in a Poisson distribution, rather than estimated.

For $y=0,1,2,\dots$ and $\lambda >0$, the Poisson distribution has massMatch this to the exponential dispersion family withThen $\mu =b^{\prime }(\theta )=e^{\theta }=\lambda$ and $V(\mu )=b^{\prime \prime }(\theta )=\mu$. The canonical link function expresses the natural parameter in terms of the mean, so the Poisson canonical link is $g(\mu )=\log \mu$. Its conditional variance equals its conditional mean.

The coefficient labelled yr2 identifies year as a factor. With $z_i=0$ in the first year and $z_i=1$ in the second, the Poisson regression assumes independent daily counts conditional on year,The unknown mean statistical parameters are the first-year log daily rate ${\beta }_0$ and the second-year log rate ratio ${\beta }_1$; the dispersion parameter is fixed at one. Approximate 95% Wald confidence intervals areExponentiating yields the first-year fitted daily mean $5.978$, with confidence interval $(5.732,6.234)$, andThe second-year fitted daily mean is $e^{1.78810-0.09816}=5.419$. Thus this Poisson regression estimates a 9.35% fall, with an approximate interval for the percentage fall from 3.66% to 14.70%. These confidence intervals rely on the Poisson distribution and independence assumptions.

The Quasi-Poisson regression retains the same conditional mean but permitsThis is a mean–variance function specification through quasi-likelihood; it does not assign a full probability distribution to each count. For independent observations, the quasi-score equation is proportional to ${\sum }_i{x}_i(Y_i-{\mu }_i)=0$, so the mean statistical parameter estimates equal those from Poisson regression. The output estimates $\hat{\varphi }=3.515351$ through the Pearson dispersion estimator, and inflates the standard errors by approximately $\sqrt{\hat{\varphi }}=1.875$.

The large residual deviance relative to 728 residual statistical degrees of freedom also signals substantial overdispersion. Daily weather, traffic and other omitted conditions may produce greater count variation than a homogeneous Poisson distribution allows. The Quasi-Poisson regression accounts for that extra marginal variance. It still requires a correct conditional mean and an appropriate independence assumption; a common dispersion parameter alone does not repair serial correlation.

Under the more plausible Quasi-Poisson regression, the year Wald statistic is $-1.686$, with reported p-value $0.0922$. An approximate 95% confidence interval isThus the estimated fall is about 9.35%, but the data do not establish a reduction at the 5% level after allowing for overdispersion. The interval allows both a sizeable reduction and a small increase. The before–after comparison also lacks a contemporaneous randomized control, so causal inference about the campaign would require addressing other changes between years. The small Poisson regression p-value is not sufficient evidence of a campaign effect when its variance assumption is unsuitable.

Measurements from the same person share baseline strength and likely share an individual response to training. An ordinary linear model with independent errors treats all 600 readings as independent conditional on week; that fails to represent their within-person correlation. The regression residuals can remain associated even when the population mean is linear in time. A random intercept captures differences in baseline strength, and a random slope can additionally describe differences in progress.

For subject $i=1,\dots ,100$ and $t_j\in \{0,2,4,6,8,10\}$, the random-intercept linear mixed model isAll subject random effects and measurement errors are mutually independent. The fixed effects ${\beta }_0,{\beta }_1$ describe the population mean baseline and weekly gain. Conditional on $b_i$, a person's observations have independent normal distributions; after integrating out $b_i$, their covariance is ${\tau }^2$ at distinct weeks and their common variance is ${\tau }^2+{\sigma }^2$. Therefore their correlation is ${\tau }^2/({\tau }^2+{\sigma }^2)$.

This Gaussian linear mixed model permits correlated repeated readings while keeping different subjects independent. The fitted standard deviations are $\hat{\tau }=21.67593$ kg and $\hat{\sigma }=14.23496$ kg, giving within-person correlation about $0.699$. All persons still have the same latent weekly slope in this model.

Allow each subject to have a separate baseline and slope through a correlated random-intercept and random-slope model:The subject random effects are independent of all measurement errors, and subjects are independent. The within-person covariance between times $s,t$ is ${\tau }_0^2+(s+t)\rho {\tau }_0{\tau }_1+st{\tau }_1^2$, with an additional ${\sigma }^2$ when the same reading is used twice.

Prefer the model with a random slope. The maximum likelihood estimation refits improve twice the log-likelihood by $227.9604$ while adding two covariance statistical parameters. Both the Akaike information criterion ($5165.779$ to $4941.819$) and Bayesian information criterion ($5183.367$ to $4968.200$) strongly favour it. The printed nominal likelihood-ratio test is overwhelming. The usual reference chi-squared distribution with two statistical degrees of freedom is not an exact regular calibration, since zero slope variance is a boundary and the intercept–slope correlation is then unidentified; a design-specific parametric bootstrap could calibrate it. This qualification does not undermine the substantial descriptive improvement shown by both information criteria.

The preferred fit estimates population mean baseline strength ${\hat{\beta }}_0=61.09286$ kg and weekly gain ${\hat{\beta }}_1=2.83143$ kg/week. The between-person baseline standard deviation is ${\hat{\tau }}_0=15.785582$ kg, and the between-person slope standard deviation is ${\hat{\tau }}_1=2.738666$ kg/week. Their estimated correlation $\hat{\rho }=0.117$ is weakly positive, giving random-effect covariance about $5.058$ kg${}^2$/week; its uncertainty is not supplied. The measurement-error standard deviation is $\hat{\sigma }=9.923280$ kg. The printed $-0.038$ instead describes the correlation between the estimated fixed effects, not between the subject random effects. These statistical parameter estimates come from restricted maximum likelihood; the model comparison uses the separate maximum likelihood estimation fits.

The population mean fitted trajectory is $61.09286+2.83143t$ kg. HenceUsing the reported slope standard error, an approximate 95% confidence interval for the population mean gain is $28.3143\pm 1.96(2.984464)=(22.46,34.16)$ kg; a Student t confidence interval with 499 statistical degrees of freedom is almost identical.

Successful individuals can improve far more than the average. Their latent ten-week gains have fitted normal distributionA clear estimate for an upper-performing group uses a between-person slope quantile: the 95th percentile is $28.3143+1.645(27.38666)\simeq 73.4$ kg, and the 97.5th percentile is about $82.0$ kg. The upper 5% of fitted underlying gains begin around 73 kg. These are person-to-person performance quantiles, not confidence intervals for the population mean. Identifying the best observed trainee would require that person's data or fitted subject random effects; the aggregate output cannot identify a literal maximum. If performance means the observed difference of endpoint readings, add the measurement-error variance $2{\sigma }^2$ to the latent gain variance.

Observed heterogeneity is variation in outcome propensities explained by measured covariates, such as age or sex. In logistic regression, subjects with different recorded predictor values may have different success probabilities, even before allowing for their previous outcomes. This is heterogeneity visible through observed predictors.

Unobserved heterogeneity is persistent variation between subjects arising from unmeasured characteristics. A subject-specific latent variable or random intercept can represent it. Subjects with high latent success propensities tend to succeed repeatedly, so their observed outcomes can have positive serial correlation even when outcomes are conditionally independent given that propensity. Apparent persistence therefore need not imply true state dependence.

True state dependence, also called true contagion, means that a previous outcome changes the distribution of a subsequent outcome after controlling both measured covariates and persistent unobserved heterogeneity. For example, an earlier successful week may make later success more likely through habit formation. A lagged-outcome coefficient in an inadequate statistical model can also reflect omitted subject differences, so positive observed persistence alone does not distinguish true state dependence from unobserved heterogeneity.

Let $Y_{it}\in \{0,1\}$ indicate a smoking-free week, $S_i$ the recorded sex indicator, $A_i$ age in years and $T_i$ assigned treatment. The screening history supplies $Y_{i0}=0$; earlier lag initializations are also zero. A history-dependent logistic regression for the first model isHere $t=1,\dots ,10$, ${\mathcal{H}}_{i,t-1}$ is the observed past, and the five unknown regression coefficients have time-invariant values. Conditional on baseline covariates, this model has the first-order Markov property: only the immediately preceding outcome enters the current conditional probability. Distinct subjects have independent histories. There is no subject random effect or additional time trend in this fit.

Using the chain rule for probabilities, the individual conditional likelihood isThe lagged values in ${\eta }_{it}$ are the individual's actual preceding outcomes. This product is a sequential conditional likelihood, not an assertion of unconditional independence of the ten readings. The full conditional likelihood is ${\prod }_i{L}_i$; no extra Bernoulli distribution factor is attached to the fixed screening history.

The second logistic regression adds ${\gamma }_2{Y}_{i,t-2}$ while retaining the first lag and baseline covariates. The models are nested under $H_0:{\gamma }_2=0$. Their likelihood-ratio test statistic is the reduction in binomial deviance,Under the null and regular large-sample conditions for the correctly specified conditional likelihood, this has an approximate chi-squared distribution with one statistical degree of freedom. Since $9.4>3.841$, reject at 5%; the approximate p-value is $0.0022$. Prefer the two-lag model to the one-lag model. Its additional lag captures statistically useful information in the history.

The cumulative-response logistic model uses $C_{it}={\sum }_{s<t}{Y}_{is}$ instead of the two individual lag predictors. It and the two-lag history-dependent logistic regression are nonnested: the entire accumulated history generally cannot be represented using only the two recent outcomes. Their binomial deviance difference therefore has no ordinary nested chi-squared distribution calibration.

Use the Akaike information criterion, which up to the same saturated-model constant equals $D+2k$ for these Bernoulli distribution conditional likelihoods. The two-lag fit has six coefficients, while the cumulative fit has five:Thus prefer the cumulative-history model, whose Akaike information criterion is lower by 35 despite its smaller number of statistical parameters. This is a model-selection comparison of conditional histories, not a nested likelihood-ratio test.

The preferred cumulative-response logistic model estimatesIts logit link describes conditional smoking-free probability given baseline predictors and prior successful weeks. Holding the other predictors fixed, males have $e^{-0.364092}=0.695$ times the female success odds; the reported p-value is $0.0152$, and the approximate 95% odds ratio confidence interval is $(0.518,0.932)$. Age has estimated odds ratio $e^{0.001224}=1.0012$ per additional year, with p-value $0.918$, giving little evidence for an age association in this fit.

The combined treatment has conditional success odds ratio $e^{0.327118}=1.387$ versus the reference treatment, with approximate 95% confidence interval $(1.035,1.859)$ and p-value $0.0285$. The fitted conditional odds are about 39% higher for the combined treatment. The trial randomization supports treatment comparisons, but conditioning on accumulated post-treatment outcomes means this coefficient is not directly the marginal total treatment effect.

Each previous successful week multiplies current success odds by $e^{0.399296}=1.491$, with approximate 95% confidence interval $(1.329,1.673)$ and very small p-value $1.06\times 10^{-11}$. This is strong fitted persistence. It can reflect true state dependence, unobserved heterogeneity, or an omitted calendar-time trend; this fit alone cannot distinguish them. The baseline intercept implies success probability ${\mathrm{logit}}^{-1}(-1.417630)\simeq 0.195$ for a reference-treatment female aged zero with no previous success. That age is outside the study's useful interpretation range, so the intercept chiefly anchors the regression. The Bernoulli distribution dispersion parameter is fixed at one, and the residual binomial deviance is 1122 on 995 statistical degrees of freedom. With individual binary outcomes, comparing that binomial deviance mechanically to a chi-squared distribution is not a reliable general goodness-of-fit test.

Take the event to mean three smoking-free weeks in succession, $Y_{i1}=Y_{i2}=Y_{i3}=1$. The initial cumulative count is zero. Along this path it equals $0,1,2$ in weeks one, two, three, respectively. In the cumulative-response logistic model, defineThe chain rule for probabilities then givesThe three conditional probabilities are approximately $0.19891,0.27015,0.35559$. If “stop during the first three weeks” instead means at least one smoking-free week by week three, its different event has probability $1-(1-{\pi }_0{)}^3\simeq 0.48590$: along the all-failure path the cumulative count stays zero. Stating the event resolves this wording ambiguity.

Let $q_{rs}$ be the transition intensity from state $r$ to $s$, and let $Q$ be the transition intensity matrix, with $q_{rr}=-{\sum }_{s\ne r}{q}_{rs}$. For a continuous-time multi-state model with the time-homogeneous Markov property, the transition probability matrix is $P(u)=e^{uQ}$, with entry $p_{rs}(u)$. Condition on each recorded initial state. Panel visits contribute transition probabilities; an exact entry into the absorbing death state contributes a statistical probability densityThis mixed panel and exact-death likelihood sums over the living state just before death. It accounts for survival until the event; replacing its final factor by $p_{r3}(u)$ would count deaths throughout the interval.

Using the actual visit times recovered from the PDF, the three individual likelihood contributions areHere $g_{23}$ means $g_{r3}$ with $r=2$, not a transition probability. Subjects 8 and 9 supply no event-density factor after their last panel observation. Assume independent subjects and noninformative examination and censoring times; conditional on their observation schedule, its distribution supplies no additional statistical parameter-dependent factor. The patients' covariates can be incorporated by using their own $Q_i$ in these same expressions.

For the progressive structure used in the subsequent output, put $a=q_{12}$, $b=q_{13}$, $c=q_{23}$ and $\lambda =a+b$. The progressive illness-death model permits no recovery, soWhen $\lambda =c$, the continuous limit is $p_{12}(u)=au{e}^{-cu}$. Thus the likelihood contributions simplify toIntermediate unobserved disease transitions remain integrated into each panel transition probability.

The fitted progressive illness-death model has three allowed arrows: $1\to 2$ with estimated transition intensity $0.1849$, $1\to 3$ with $0.01935$, and $2\to 3$ with $0.06143$, all in years${}^{-1}$. State 3 is an absorbing state, and state 2 has no return arrow.

Once in state 2, the holding time has an exponential distribution with rate $q_{23}$, so the mean holding time from a transition intensity matrix isApply a confidence interval for an inverse rate to the printed rate interval $(0.03552,0.1063)$: since inversion reverses order, the approximate 95% interval for the mean isThe time-homogeneous Markov property makes the future depend only on the current state; the exponential distribution also has the memoryless property. For a person currently in state 2,Equivalently, $p_{22}(2)=p_{22}(1{)}^2\simeq 0.9404163^2$. Thus the fitted two-year death probability is approximately 11.6%.

Let $z_i=(S_i,E_i,A_i{)}^T$ contain sex, the education indicator, and age at diagnosis. Write $\overline{z}$ for the sample means. The output's baseline transition intensities are evaluated at those means, so use the centred log-linear transition intensity modelThe full transition intensity matrix isThe fitted centred baselines are $(q_{12},q_{13},q_{23})(\overline{z})=(0.1821,0.0126,0.0450)$, and fitted slope vectors in sex–education–age order areThe two sex coefficients displayed as zero are fixed by the specified constraints; they are not estimated to be exactly zero. There are three free baseline transition intensities and seven free covariate slopes. The sample means are not printed, so uncentred intercepts at $z=0$ cannot be recovered numerically from this output.

Conditional on fixed covariates, subjects follow independent, correctly classified continuous-time multi-state models obeying the time-homogeneous Markov property, with $P_i(u)=e^{u{Q}_i}$. The progression is irreversible and death absorbing, with constant transition intensities during follow-up for each subject. In particular, the model uses fixed age at diagnosis rather than attained age. It assumes noninformative observation and censoring, and treats death times as exact through the mixed panel and exact-death likelihood. The Markov property rules out an additional effect of elapsed time in the current state after conditioning on it and the recorded predictors.

Under $H_0$, the seven free covariate coefficients all equal zero. The likelihood-ratio test statistic isThe larger fit has ten free statistical parameters, compared with three in the constant-rate fit, so the reference chi-squared distribution has seven statistical degrees of freedom. Its 95th percentile is $14.06714$. Therefore reject $H_0$ at 5%; the p-value is approximately $0.0043$, and the covariate model is preferred. The two constrained sex coefficients add no free statistical parameters.

For the dementia-to-death arrow, the hazard ratio for higher versus lower education is $e^{-1.490}=0.2254$, with approximate 95% confidence interval $(e^{-2.547},e^{-0.4327})=(0.0783,0.6488)$. Thus higher education is associated with a roughly 77.5% lower fitted death transition intensity, conditional on diagnosis age; the interval excludes one. This is an observational association and is not a direct multiplicative statement about death probability.

Each additional year of age at diagnosis multiplies the dementia-to-death transition intensity by $e^{0.07984}=1.0831$, with confidence interval $(e^{-0.007765},e^{0.1674})=(0.9923,1.1822)$. The estimate suggests an increase, but the interval includes one, so this individual age coefficient is not significant at 5%. Sex has hazard ratio one for this transition by the model's imposed constraint. The output supplies no estimated sex effect or test of that constraint for this arrow. The nonzero sex coefficient for $1\to 2$ must not be mistaken for a sex effect on $2\to 3$.