Past exam of the mathematics course of the University of Cambridge / 2010 / iii / Paper 37 (Applied Statistics)

The least-squares estimator minimizes $(Y-Xb{)}^T(Y-Xb)$. Differentiation gives $X^{T}X\hat{\beta }=X^{T}Y$, and the full column rank of $X$ makes $X^{T}X$ positive definite. HenceThe fitted values are $\hat{Y}=X\hat{\beta }=HY$, where the hat matrix is $H=X(X^{T}X{)}^{-1}{X}^T$. The regression residuals are $\hat{\epsilon }=Y-\hat{Y}=(I-H)Y$. Since $H$ is the orthogonal projection matrix onto the column space of $X$, $(I-H)X=0$ andThis is a singular multivariate normal distribution supported on the $(n-p)$-dimensional orthogonal complement of the column space: its components are generally correlated. The equality $\hat{Y}=X\hat{\beta }$ puts the fitted vector in that column space, while $X^T\hat{\epsilon }=0$ expresses fitted-residual orthogonality. No residual component lies in an explanatory direction that could improve the least-squares fit.

For the calibration model, write $R_i$ for a reading and $a_i$ for the known amount. The first fit assumes $R_i=\alpha +\beta a_i+{\epsilon }_i$ with independent $N(0,{\sigma }^2)$ errors. Its residual-versus-fitted plot has systematic changes in group means and an increasing within-group spread. These suggest both a nonlinear mean relationship and heteroscedasticity, rather than the random scatter expected from a satisfactory normal linear model.

The Box–Cox transformation plot is a profile log-likelihood for the power parameter $\lambda$. Its maximum is near $0.94$; the horizontal line gives an approximate 95% likelihood-ratio confidence interval, which excludes the untransformed value $1$. Since affine rescaling of the response does not change the fitted mean space, fitting the raw power $R_i^{0.94}$ is equivalent to using $(R_i^{0.94}-1)/0.94$ with an intercept. The second model isIts residual spread is more even, but observation 17 has a conspicuously negative residual near $-9.95$. This motivates checking for an outlier or recording error. The third fit omits that observation; omission should be justified by that investigation, since choosing exclusions from residuals can affect inference. Some group pattern still remains, so the transformation alone does not certify a perfect model.

At an amount of $3$ nanograms, the third fit predicts the transformed meanSimply applying the inverse power gives ${\hat{m}}^{1/0.94}=118.5187$. The requested expected reading needs bias correction after an inverse transformation, since $\mathbb{E}[g(Z)]\ne g(\mathbb{E}[Z])$ for nonlinear $g$. With $q=1/0.94$ and $s=2.832$, a second-order Taylor expansion givesHere the correction uses the residual variance, rather than the variance of the estimated fitted mean. The power-transformed Gaussian model is an approximation for positive readings; the fitted mean is over 31 residual standard deviations above zero, so the negative tail has negligible effect on this local approximation. The ordinary inverse-transformed prediction is therefore essentially $118.52$ at the displayed precision.

The data do not cross every age with every gender: only age 21 has both genders, while age 30 has only women and age 40 only men. Thus a gender effect separate from age can be estimated under an additive model, but age-by-gender interactions cannot all be estimated or checked. There is only one observation in each recorded combination of age, gender, policy, and points, so there is no within-cell replication for separating pure error from lack of fit. No sampling method or additional risk characteristics are supplied, limiting comparisons between policyholder categories.

There is also a source discrepancy: the printed table gives $268$ for one comprehensive premium, whereas the R data vector gives $368$ at the corresponding position. The supplied regression output is treated as referring to the R data. This should be resolved against the original data before substantive use.

Let $Y_{agpq}$ be the premium at age $a\in \{21,30,40\}$, gender $g\in \{F,M\}$, policy $p\in \{3rd,\mathrm{comp}\}$, and points $q\in \{0,3,6,9\}$, for the combinations actually observed. The additive normal linear model isIts corner-point constraints are ${\alpha }_{21}={\gamma }_F={\delta }_{3rd}={\eta }_0=0$. Thus $\mu$ is the reference-category mean; the other coefficients are additive contrasts. The model assumes a common variance and excludes all interactions. It has $1+2+1+1+3=8$ estimable mean parameters and $32-8=24$ residual degrees of freedom. The variance estimate is

The reduced model combines the 0-, 3-, and 6-point categories. Relative to the full model, the null hypothesis is ${\eta }_3={\eta }_6=0$, leaving the 9-point contrast unrestricted. These are two linear restrictions on an otherwise unchanged normal linear model.

For nested normal models, the reduction in residual sum of squares divided by ${\sigma }^2$ has a chi-squared distribution with degrees of freedom equal to the number of restrictions. It is independent of the full model's residual variance estimate. Hence the F-test statistic isIt is below the supplied 5% critical value $3.402826$, so do not reject equal effects for 0, 3, and 6 points. There is insufficient evidence that separating those categories improves this additive model; this does not prove that their true premiums are identical.

Let $h(q)=1$ for $q=0,3,6$ and $h(9)=2$. The new normal linear model allows an age-by-policy interaction:The corner-point constraints are ${\alpha }_{21}={\delta }_{3rd}={\gamma }_F={\xi }_1=0$ and ${\kappa }_{21,p}={\kappa }_{a,3rd}=0$. Compared with the six-parameter reduced additive model, there are two additional interaction coefficients. Its F-test givesThus the interaction model significantly improves the reduced additive model. Estimated premiums decrease with age. Men have an estimated premium $94.375$ higher than women at fixed other characteristics, and 9 points adds $41.708$ relative to 0, 3, or 6. Comprehensive cover costs more than third-party cover, but its increment depends on age: $175.375$ at 21, $148.500$ at 30, and $79.500$ at 40. These are model-based comparisons, including extrapolations into the unobserved age-gender combinations.

For a 40-year-old woman with comprehensive cover and 6 points, the gender and points contrasts vanish. The fitted premium is

The Bernoulli distribution probability mass function isWriting ${\theta }_i=\log [p_i/(1-p_i)]$ gives $p_i=e^{{\theta }_i}/(1+e^{{\theta }_i})$ and $\log (1-p_i)=-\log (1+e^{{\theta }_i})$. Therefore the exponential family representation hasThe logit link is the canonical link.

Independence gives the log-likelihoodDifferentiating yields the score equationsWhen a finite maximum-likelihood estimator exists, multiplying these equations by ${\hat{\beta }}^T$ givesThe qualification about a finite estimator matters under complete separation in logistic regression, where coefficient estimates can escape to infinity.

For separate binary observations, the saturated model fits each outcome exactly and has maximized log-likelihood zero. With ${\hat{p}}_i=p_i(\hat{\beta })$, the residual deviance is thereforeUse the identity in part (b) to replace $y_i$ by ${\hat{p}}_i$ in its first term:Thus the deviance is twice the sum of the fitted binary information entropies. A deviance goodness-of-fit test against ${\chi }_{n-p}^2$ is unreliable for individual ungrouped binary observations: there is only one trial per saturated-model parameter, so the usual large-cell-count approximation fails. Deviance differences between fixed-dimensional nested logistic regression models can still support a likelihood-ratio test, subject to its usual regularity assumptions.

The three tree models are respectively $\mathrm{logit}{p}_i={\beta }_0$, $\mathrm{logit}{p}_i={\beta }_0+{\beta }_1{\log }_2{T}_i$, and $\mathrm{logit}{p}_i={\beta }_0+{\beta }_1{\log }_2{T}_i+{\beta }_2{S}_i$, with independent Bernoulli random variables conditional on the covariates. Testing the additional severity effect givesUnder $H_0:{\beta }_2=0$, its reference distribution is asymptotically ${\chi }_1^2$. The $p$-value is approximately $1.21\times 10^{-21}$, giving overwhelming evidence that severity improves the fitted model. Holding severity fixed, doubling $T$ increases ${\log }_{2}T$ by one and multiplies the odds byThis is an odds ratio, not a multiplication of the event probability itself.

The line compares the two pooled accident rates, dividing the total number of accidents by total observation years in each period. It gives the unadjusted incidence rate ratiocorresponding to an estimated reduction of about $50.3\%$. Pooling does not adjust for differences between locations or for the unequal location composition of observation time before and after installation.

Let $Y_{it}$ be the accident count at location $i$ in period $t\in \{0,1\}$, and $e_{it}$ its exposure in years, where $t=1$ denotes after installation. The Poisson exposure model isThe offset $\log e_{it}$ has its coefficient fixed at one. The model assumes constant rates within each period, a common before rate across locations, a common treatment rate ratio, independent counts, and Poisson variance equal to the mean. The incidence rate ratio is $r=e^{\beta }$, soAn approximate 95% Wald confidence interval transforms the interval for $\beta$:

The second Poisson exposure model allows each location its own baseline accident rate:The last equality is the corner-point constraint. The common after-to-before incidence rate ratio remains $e^{\beta }$. This model is appropriate because the sites can have different background risks; treating those differences as unexplained Poisson variation made the first model fit poorly, with deviance $50.863$ on 14 residual degrees of freedom.

Adding seven location effects to the first fitted model reduces the deviance by $50.863-16.275=34.588$, which is large relative to ${\chi }_7^2$. The supplied sequential analysis of deviance reports a slightly different location contribution because it enters locations before the period indicator. The full model is much better but its deviance $16.275$ on 7 degrees of freedom still gives an approximate goodness-of-fit $p$-value $0.0227$. That suggests residual lack of fit or overdispersion, although sparse accident counts make the asymptotic deviance goodness-of-fit test approximate.

After adjusting for location, adding the period effect reduces deviance by $9.750$ on one degree of freedom, giving the supplied likelihood-ratio test $p\simeq 0.002$. The adjusted incidence rate ratio isUnder this Poisson exposure model, the accident rate after installation is estimated to be about $54.2\%$ lower, with a statistically significant reduction. The remaining lack of fit merits residual checks and possibly a Quasi-Poisson model or a model allowing location-specific changes. This before-and-after observational study establishes an association under the fitted model; changes in traffic volume, broader trends, or regression to the mean at sites chosen for high prior counts could also contribute, so the estimated reduction is not by itself a causal treatment effect.

Let $I\in \{1,2\}$ denote eventual death or recovery, respectively; $T$ is the time from admission to that terminal event, and $\theta =\mathrm{Pr}(I=1)$. A patient still hospitalized at collection has right censoring: only $T>c$ is known, where $c$ is the elapsed time since admission. Both eventual outcomes are then unobserved. A recorded death or discharge is an observed terminal event in this competing risks model; discharge is not censoring of an otherwise observable later in-hospital death. Assume the observation deadline, conditional on relevant admission information, is noninformative for $(T,I)$.

Let $f_j(t)=f(t\mid I=j)$ and $S_j(t)=1-F(t\mid I=j)$. A recorded death observes both $T=t$ and $I=1$, so its likelihood contribution, omitting factors from the independent observation mechanism, isIts integration over $t$ gives the eventual death probability $\theta$.

A recorded discharge observes $T=t$ and $I=2$. Its likelihood contribution isThe mixing probability distinguishes the joint outcome-and-time density from the conditional recovery-time density.

For a patient still in hospital, sum over the unobserved eventual outcome. Its likelihood contribution is the mixture survival functionThe probability of eventual death among those still hospitalized is generally different from $\theta$, because the two conditional time distributions can have different tails.

This is an EM algorithm for a censored lognormal competing-risks mixture. Write $Y_i=\log T_i$ and suppose $Y_i\mid I_i=j\sim N({\mu }_j,{\sigma }_j^2)$, with mixing weights ${\pi }_1=\theta$, ${\pi }_2=1-\theta$. The missing data are $I_i$ and $Y_i$ for censored patients; neither is missing for recorded deaths or recoveries. The complete-data log-likelihood, apart from terms independent of the parameters, isThe lognormal Jacobian $-\log T_i$ is also parameter-independent and therefore does not affect the maximization.

At iteration $k$, set $z_{ij}=(\log c_i-{\mu }_j^{(k)})/{\sigma }_j^{(k)}$ for each censored patient and let $\overline{\Phi }(z)=1-\Phi (z)$. Bayes' rule gives the E-step class weightsGiven that class, $Y_i$ has a truncated normal distribution above $\log c_i$. With the Inverse Mills ratio $\psi (z)=\phi (z)/\overline{\Phi }(z)$, compute both moments needed for the E-step:Thus $\mathbb{E}[1_{\{I_i=j\}}\mid \text{data}]=w_{ij}$, $\mathbb{E}[1_{\{I_i=j\}}{Y}_i\mid \text{data}]=w_{ij}{m}_{ij}$, and the corresponding second moment is $w_{ij}{s}_{ij}$. For an observed outcome $j_i$ at time $t_i$, use $w_{ij}=1_{\{j=j_i\}}$, $m_{ij}=\log t_i$, and $s_{ij}=(\log t_i{)}^2$. These specify the full conditional expectation of ${\ell }_c$.

Put $W_j={\sum }_i{w}_{ij}$ and let $n$ be the total number of patients. Maximizing that expected complete-data log-likelihood gives the M-stepAll moments on the right are evaluated at the old parameters. Repeat the E-step and M-step until the observed log-likelihood stabilizes, using positive initial variances and mixing weights. The expectation-maximization algorithm increases that likelihood at each exact iteration, but need not find its global maximum; multiple starts help distinguish local optima. The resulting estimates give the case fatality probability $\hat{\theta }$ andThese are conditional event-time distributions, distinct from the cumulative incidence functions $\hat{\theta }\hat{F}(t\mid I=1)$ and $(1-\hat{\theta })\hat{F}(t\mid I=2)$.