Past exam of the mathematics course of the University of Cambridge / 2002 / iii / Paper 36 (Applied Statistics)

Write $1$ for the all-ones vector and setThe centering assumption makes $P_0{P}_X=0$. Thus the full normal linear model has fitted-value orthogonal projection $P_0+P_X$, and its residual sum of squares is $R=Y^{\mathsf{T}}(I-P_0-P_X)Y$. Assume $\nu >0$, as is needed to estimate the unknown error variance from residuals.

For the null that all slopes vanish, the reduced fitted value is $\bar{Y}1$, with residual sum of squares $R_0=Y^{\mathsf{T}}(I-P_0)Y$. The F-test isReject for an upper-tail value exceeding the $1-\alpha$ quantile of this F-distribution. The required theorem is Cochran's theorem: for an isotropic normal distribution, projections onto orthogonal subspaces are independent, and the squared norm on a rank-$r$ subspace divided by ${\sigma }^2$ has a chi-squared distribution with $r$ degrees of freedom. Here the null makes the projected mean in $\mathrm{col}(X)$ zero, so $(R_0-R)/{\sigma }^2\sim {\chi }_p^2$, independently of $R/{\sigma }^2\sim {\chi }_{\nu }^2$.

Let $p_1,p_2$ be the numbers of columns in the two predictor blocks. To test the first block while retaining the second as nuisance, formThe partial F-test for nested linear models isAgain reject in the upper tail. The difference of the nested fitted-value projections has rank $p_1$, is orthogonal to the full residual projection, and annihilates the null mean, which gives the same independent chi-squared argument. If a block has no columns, the corresponding test is vacuous.

The phrase orthogonal statistical parameters refers to zero cross-block Fisher information, not to the numerical parameter vectors being perpendicular. In this model the cross-block information is $X_1^{\mathsf{T}}{X}_2/{\sigma }^2$, so the orthogonal coefficient blocks in a centered normal linear model condition isThen the two fitted subspaces are orthogonal, the block least-squares estimators have zero covariance and are independent because they are jointly Gaussian, and each block's extra sum of squares is unchanged by including the other block. Centering also makes both blocks orthogonal to the intercept.

Both fits are the same additive two-way analysis of variance model: a common intercept, five country effects represented by four contrasts, and eight item effects represented by seven contrasts. The three missing cells are omitted, leaving $37$ observations and model rank $1+4+7=12$, hence $25$ residual degrees of freedom. The residual sum of squares is $659.44$ in either ordering, giving

These are sequential sums of squares, so each row measures the improvement when that factor is added after the preceding factors. The country-first sum $1115.56$ describes its unadjusted contribution before accounting for item; its displayed $F=10.573$ and $p=3.73\times 10^{-5}$ are not the country-adjusted-for-item test. When country is added after item, its extra sum is $1616.74$, givingThus the additive normal linear model gives strong evidence of country differences after accounting for the different products. Conversely, item added after country has extra sum $16910.20$, soThe printed zero for this probability is numerical formatting, not a probability literally equal to zero. Product differences are much larger than residual variation, as expected from their different scales.

The explanation for the order dependence is that missing cells can destroy factor orthogonality in two-way ANOVA. A complete equally replicated crossed design would make centered country and item indicator columns orthogonal. Here the missing observations change the item mix across countries and the country mix across items, so the two factor subspaces are not orthogonal. The two sequential partitions still have the same total explained sum: $1115.56+16910.20=16409.02+1616.74=18025.76$ to the displayed precision. For inference about either factor conditional on the other, compare the full model with the reduced model that omits that factor, rather than interpreting its first-entry sum as adjusted evidence.

The omnibus country test does not identify the particular country contrasts. As an extension, fitting the supplied price table gives the item-adjusted UK-minus-US contrast about $21.33$ pounds, with standard error $2.824$ and an ordinary Gaussian-model $95\%$ confidence interval approximately $[15.51,27.14]$ pounds. This is a derived model contrast, with Student t-distribution reference on $25$ degrees of freedom; simultaneous exploration of many country pairs requires control for multiple hypothesis testing.

Check residuals against fitted values, product and country, together with a normal Q-Q plot and influential observations. The homoscedasticity assumption on absolute prices may be doubtful because expensive products can have larger absolute variation. A model for logarithmic prices gives multiplicative country effects and can make the scale assumption more plausible. Country-by-item interaction terms would allow discounts to differ by product, but there is at most one observed value per cell: an unrestricted interaction saturates the $37$ observed means and leaves no pure-error degrees of freedom, so it cannot be tested against independent measurement noise from these data alone. Additional products or genuine cell replication would support that extension. The selected products and non-random missing entries also limit claims about a whole country's consumer prices; the conditional inference is for the stated additive model and observed price collection.

For logistic regression, set ${\eta }_i=x_i^{\mathsf{T}}\beta$ and $p_i=(1+e^{-{\eta }_i}{)}^{-1}$. Its Bernoulli log-likelihood, score function and Fisher information arewhere $W=\mathrm{diag}\{p_i(1-p_i)\}$. Solve $U(\beta )=0$ by Fisher scoring or iteratively reweighted least squares. One scoring step isEquivalently regress the working response $z_i={\eta }_i+(y_i-p_i)/[p_i(1-p_i)]$ on $X$ using weights $p_i(1-p_i)$. Iterate until the likelihood and coefficients stabilize, with rank and convergence checks. Separation in logistic regression can prevent a finite maximum: if a predictor direction strictly separates the zeros and ones, moving indefinitely along it improves the likelihood rather than producing a finite root.

Under the unrestricted model, each individual probability is estimated independently as ${\hat{p}}_i=y_i$. Every observed binary outcome then has probability one, so, using boundary limits,The resulting residual deviance is $D=-2\ell (\hat{\beta })$. It is a legitimate likelihood contrast and can be used to compare fitted models. However, individual Bernoulli deviance need not have a chi-squared calibration: the saturated model has one boundary-valued parameter for each single observation, and its dimension grows with $n$. The fixed-dimensional regularity behind Wilks theorem does not justify a ${\chi }_{n-p}^2$ reference here.

For a concrete demonstration, let the true model contain only an intercept and have success probability $1/2$. Its fitted probability is $\bar{y}$, andwhereas a ${\chi }_{n-1}^2$ variable divided by $n$ converges to one. Thus treating the raw Bernoulli residual deviance as an ordinary absolute goodness-of-fit statistic can reject a correct model systematically. In contrast, deviance differences between fixed-dimensional nested logistic models have their usual asymptotic chi-squared reference under the relevant regularity assumptions. For absolute fit use meaningful replicated covariate groups and grouped-binomial logistic regression, or assess residual patterns and prediction on held-out data; a model-specific parametric bootstrap can calibrate a chosen statistic.

The family-history table contains $462$ men, of whom $160$ have the response present, giving an overall observed fraction $160/462=0.3463$. The unadjusted fractions are $96/192=0.5$ with a positive family history and $64/270=0.2370$ without one. Their crude odds ratio is $(96/96)/(64/206)=3.219$; this is not the covariate-adjusted effect from the fitted logistic regression.

The full model has an intercept and nine slopes, hence $462-10=452$ residual degrees of freedom. The binomial distribution dispersion is fixed at one. Although the software labels coefficient-to-standard-error ratios as t values, here they are asymptotic normal Wald statistics, not statistics with an estimated Gaussian residual variance and an exact Student reference. The full conditional coefficients provide evidence for tobacco, cholesterol, family history, behaviour score and age: their absolute Wald ratios are approximately $2.99,2.92,4.06,3.22,3.73$. The other four ratios are smaller in magnitude, especially alcohol at about $0.027$. These are conditional associations; a weak adjusted coefficient does not by itself establish absence of an underlying association, and multicollinearity can affect the estimates.

The null deviance is $596.1084$ on $461$ degrees of freedom and the full residual deviance is $472.14$ on $452$. Their difference $123.9684$ has a nominal nine-degree-of-freedom likelihood-ratio test reference for the nine slopes jointly. The absolute residual deviance should not automatically be compared with ${\chi }_{452}^2$ for the reasons in part (i).

The shown stepwise selection by the Akaike information criterion follows a backward path. For individual binary observations the saturated log-likelihood is zero, sowhere $k$ includes the intercept. Removing a one-parameter term improves AIC exactly when the increased deviance is less than two. The successive deletions are alcohol, adiposity, blood pressure and obesity. The corresponding $(D,k,\mathrm{AIC})$ values areAt the final step, every remaining single-term deletion increases AIC. This is a local stopping condition of the search, not a proof that it found the best subset among all possible models. Obesity is retained for several steps because its successive deletion costs differ after the other terms are removed; the procedure does not simply sort the initial Wald ratios.

For the selected fit, write $H=1$ for positive family history and $H=0$ otherwise. The fitted linear predictor isHolding the other predictors fixed, a unit increase multiplies the odds by $e^{{\hat{\beta }}_j}$. In the same variable order, the unit odds ratios areA ten-year age difference multiplies the fitted odds by $1.656$, and ten behaviour-score points by $1.449$. Family history has adjusted odds ratio about $2.48$, with ordinary fitted-model $95\%$ confidence interval $\exp (0.908171\pm 1.96\times 0.225603)\simeq [1.59,3.86]$. The intercept is the log odds at zero numerical covariates and absent family history; that extrapolated baseline is not a typical man's risk. Odds ratios are not risk ratios.

The selected model's residual degrees of freedom are $456$, consistent with its six coefficients. All five retained slopes have nominal normal-Wald $p$ values below $0.004$ in that model, but these tests and confidence intervals do not incorporate the prior model selection. Its improvement over the intercept-only fit is $596.1084-475.6856=120.4228$ on five nominal degrees of freedom; it loses only $3.5456$ deviance units relative to the nine-slope fit while using four fewer parameters. Selection optimism remains relevant to predictive performance. Check possible nonlinear continuous-predictor effects, influential observations and clinically motivated interactions; assess out-of-sample discrimination and agreement of predicted probabilities with observed frequencies by cross-validation or external data. The printed output alone does not supply those diagnostics, and the selected high-risk male population limits extrapolation to other populations. Nothing here establishes causal effects of these observational predictors.

Use $e_{ij}$ for the expected accident count, not for the accident rate. With a design row $z_{ij}$ consisting of an intercept, seven non-reference site indicators and the after-treatment indicator,Here $\log p_{ij}$ is a known generalized linear model offset, with coefficient fixed at one. Differentiating the Poisson regression log-likelihood gives the score function equationsThese nine equations estimate the intercept, seven site contrasts and one treatment contrast, with the stated reference-level constraints.

The Fisher information is $Z^{\mathsf{T}}WZ$, where $W=\mathrm{diag}(e_{ij})$. A Fisher scoring update adds $(Z^{\mathsf{T}}WZ{)}^{-1}{Z}^{\mathsf{T}}(y-e)$ to the current coefficient vector. In iteratively reweighted least squares, the working response is $z_{ij}^{\ast }=\log e_{ij}+(y_{ij}-e_{ij})/e_{ij}$; regress $z_{ij}^{\ast }-\log p_{ij}$ on the design with weights $e_{ij}$. Thus the reported four iterations solve the likelihood equations while accounting for unequal exposure periods.

The common adjusted rate ratio is obtained by exponentiating the after coefficient:This estimates a $54.2\%$ lower accident rate after the intervention, conditional on the additive site model. The standard error $0.2751810$ gives normal Wald statistic $-2.8369$, two-sided $p\simeq 0.0046$, and an approximate $95\%$ confidence interval for the rate ratioThe label t value is again a normal-reference Wald ratio under fixed Poisson dispersion, rather than an exact Student statistic.

The intercept gives the before rate at reference site one: $e^{0.2707792}\simeq 1.311$ accidents per year. The seven site coefficients are log baseline-rate ratios relative to that site. In site order two through eight, their exponentials are about $0.615,2.767,1.711,0.769,1.797,0.615,1.221$. The third site has the largest estimated rate and a nominal Wald ratio about $3.118$ against the reference. These are adjusted baseline comparisons, not differences in absolute accident counts. The overall unadjusted before/after rate ratio is $(15/18)/(114/68)\simeq 0.497$, and it differs from $0.458$ because site risks and their relative exposure weights differ.

An independent way to solve the fit is the profile likelihood for a common Poisson rate ratio with unequal exposures. Put $t_i=y_{i1}+y_{i2}$ and ${\lambda }_i=\exp (\mu +{\alpha }_i)$. For fixed $r$, the site score gives ${\hat{\lambda }}_i=t_i/(p_{i1}+r{p}_{i2})$. The remaining equation isIts positive root reproduces the reported treatment estimate. Comparing the same site model with $r=1$ to the fitted model gives a deviance improvement about $9.750$ on one degree of freedom, with nominal likelihood-ratio $p\simeq 0.0018$.

There are $16$ count observations and $9$ fitted mean parameters, leaving $7$ residual degrees of freedom. The null deviance $132.9485$ on $15$ degrees of freedom is for the intercept-only rate model, still including the exposure offset. The residual Poisson deviance is $16.27524$ on $7$ degrees of freedom. A chi-squared comparison gives approximately $p=0.023$, suggesting potential lack of fit of a single common treatment ratio; the extreme deviance residuals near $-2.03$ and $2.14$ also deserve inspection. Several after-cell expected counts are below one or near two, so this absolute-fit calibration is only approximate. Site-specific treatment effects, overdispersion or temporal dependence are possible explanations to investigate, with a parametric bootstrap available to calibrate the sparse-count fit statistic. A before/after association alone does not isolate a causal effect from secular changes or regression to the mean at selected sites.

The Poisson regression margin-matching score equations from the site indicator for each $i\ge 2$ directly giveSubtracting these seven row balances from the intercept equation ${\sum }_{ij}{e}_{ij}={\sum }_{ij}{y}_{ij}$ gives the missing reference row balance. ThereforeThe after-indicator score similarly gives ${\sum }_i{e}_{i2}={\sum }_i{y}_{i2}$. Subtracting it from the intercept balance yields the before-column balance, soThus the fitted before and after count totals are $114$ and $15$, and each fitted site total equals its observed total. These identities concern expected counts $e_{ij}=p_{ij}{\hat{\mu }}_{ij}$; unequal exposures do not imply corresponding unweighted sums of fitted rates equal sums of observed rates. They follow from the unpenalized canonical Poisson likelihood and the available intercept/factor columns, not from balanced exposure periods.

The first approach specifies marginal means, variances and an exchangeable correlation for each student's repeated counts. It is a moment model suitable for a generalized estimating equation; equality of the mean and variance alone does not specify a full Poisson joint distribution. Different students supply independent sampling clusters, while the three counts within a student are dependent.

The second approach is a Poisson generalized linear mixed model with a positive shared random multiplier $U_i=e^{b_i}$. Conditional on it, the three counts are independent Poisson variables; integrating over it gives a fully specified dependent joint law. The gamma random-intercept Poisson model uses Gamma shape ${\tau }^2/\theta$ and rate $\tau /\theta$, so $\mathbb{E}{U}_i=\tau$ and $\mathrm{Var}(U_i)=\theta$. Treat the covariates as fixed, with the same multiplier distribution across covariate/treatment groups, as required by this model.

Let $a_{ij}=\exp ({\beta }_0+{\beta }^{\mathsf{T}}{x}_{ij})$ and $m_{ij}=\mathbb{E}{Y}_{ij}$. The law of total expectation, law of total variance and law of total covariance giveThus each marginal is a negative binomial distribution by the Poisson-gamma mixture, rather than generally a Poisson variable. Put $\kappa =\theta /{\tau }^2$. Its variance is $m_{ij}+\kappa m_{ij}^2$, and the marginal correlation isIt is positive but usually depends on the two marginal means, so it need not be the common exchangeable correlation proposed in the first approach. The first method focuses on population mean effects with a chosen working correlation matrix; the second models latent student heterogeneity and separates conditional Poisson variation from additional marginal variation. Correctly specified likelihood inference for the second uses more distributional assumptions, while the first can retain valid mean inference despite a wrong working covariance by using a cluster-level sandwich variance.

For the marginal mean model, $e^{{\beta }_0}$ is the expected term count for a student with all recorded covariates at their reference or zero values. If the treatment indicator is one for the new therapy and zero for the reference treatment, then $e^{{\beta }_1}$ is the ratio of population mean counts under the two treatments, holding the remaining covariates fixed. Values below one mean a reduced expected count; $100(1-e^{{\beta }_1})\%$ is the corresponding percentage reduction.

For the random-intercept model, $e^{{\beta }_0}$ is the conditional reference count at $b_i=0$, or multiplier $U_i=1$; the conditional baseline for student $i$ is $U_i{e}^{{\beta }_0}$. It is not automatically the mean-population baseline, which is $\tau e^{{\beta }_0}$. The conditional treatment rate ratio for students with the same multiplier is $e^{{\beta }_1}$. Because the multiplier has the same distribution in the treatment groups and the mean uses a logarithmic link function, that is also the marginal treatment rate ratio. This is the property that marginal and conditional slopes agree for an independent log-link random intercept.

The parameter $\theta$ is the variance of the positive student multiplier, not a Poisson sampling variance or a treatment coefficient. Greater $\theta$ at fixed $\tau$ means greater heterogeneity in underlying propensity and larger between-term covariance. Its extra contribution to the count variance is $\theta a_{ij}^2$. If $\tau =1$ is imposed, $\theta$ is also the relative multiplier variance; otherwise the relative heterogeneity is $\theta /{\tau }^2$.

There is a necessary identifiability qualification: if the multiplier mean $\tau$ is also unrestricted, scale identifiability in a gamma random-intercept Poisson model shows that replacingleaves the response law unchanged. Consequently only ${\beta }_0+\log \tau$ and $\theta /{\tau }^2$, together with the slopes, are identified from these counts without a scale convention. Fixing the multiplier mean, usually to one, is needed to interpret the conditional intercept and heterogeneity scale separately.

Integrating the conditional mean in the gamma random-intercept Poisson model givesand henceThe marginal model's logarithmic mean family is therefore correctly specified, with a shifted intercept and identical slopes. Under the usual independent-cluster and estimating-equation regularity conditions, fitting that marginal mean consistently estimatesThus the treatment and other slopes estimate the intended marginal log rate ratios; indeed they also coincide with the conditional slopes under the stated common random-effect distribution. The intercept consistently estimates the marginal baseline, but not the conditional $b_i=0$ intercept unless $\tau =1$. If the first investigator interprets the intercept as a conditional latent-student baseline without recognizing this shift, that interpretation is incorrect.

The assumed Poisson marginal variance and constant correlation generally fail under this random-effect law. They need not invalidate mean-coefficient consistency for a suitable generalized estimating equation, because the expected residual vector is still zero at the true marginal mean. They do invalidate naive standard errors based only on that assumed covariance; robust inference requires the correction in part (iv). The simple unchanged-slope result relies on the random multiplier distribution not changing with covariates, rather than on a universal equivalence of marginal and conditional regression models.

Use a generalized estimating equation with any well-behaved positive-definite working correlation matrix. Let $\gamma$ collect the marginal intercept and slopes, let $m_i(\gamma )$ be the vector of three correct marginal means, and put $D_i=\partial m_i/\partial {\gamma }^{\mathsf{T}}$. Write the working covariance as $V_i=A_i^{1/2}{R}_i(\alpha )A_i^{1/2}$, where the proposed Poisson variance gives $A_i=\mathrm{diag}(m_i)$; working independence is also allowed. SolveAt the true mean, each cluster residual has expectation zero, so the estimating equation is unbiased regardless of whether $V_i$ equals the actual covariance. With independent students, sufficient covariate variation and regularity as $m\to \infty$, this gives consistent mean coefficients.

Replace the model-based covariance by the cluster-level sandwich covariance matrix. Evaluated at the fitted values, defineThenThe residual outer products estimate the actual within-student variability, including correlations and extra-Poisson variance absent from the working model. Use these standard errors for asymptotic Wald tests and confidence intervals; exponentiating the treatment coefficient and its interval gives the mean-count rate ratio and interval. Treating the $3m$ observations as independent when forming the sandwich meat would miss the within-student dependence. The effective asymptotic replication is the number of students, so a small number of clusters requires further finite-sample care rather than assuming that three measurements per student repair the large-sample approximation.

For independent Poisson counts, the log-likelihood at mean vector $\mu$ isThe saturated model maximizes each cell separately at ${\hat{\mu }}_i=y_i$, including the boundary value zero. If $e_i$ is the fitted mean under the log-linear regression, twice the saturated-minus-fitted log-likelihood is the full Poisson deviancewith $0\log (0/e_i)=0$ understood by continuity. The fitted intercept has score function $\partial \ell /\partial \mu ={\sum }_i(y_i-e_i)$, since differentiating the logarithmic mean with respect to that intercept gives one. Its maximum-likelihood equation therefore gives ${\sum }_i{e}_i={\sum }_i{y}_i$. This is why Poisson deviance simplifies when an intercept is fitted:Without the fitted unpenalized intercept, one cannot generally discard the two linear terms. The individual logarithmic terms in the simplified sum need not each be positive; nonnegativity belongs to the full deviance, or to the full cell contributions before cancellation.

For the approximation, write $y_i=e_i+r_i$ and expand $\log (1+r_i/e_i)$. When relative residuals are small,Consequently the quadratic Pearson approximation to the Poisson deviance isthe Pearson chi-squared statistic for unit Poisson dispersion. Adequate expected counts and small relative residuals are what make the quadratic approximation useful; it is not an exact identity and can be poor for sparse cells, especially a zero count with positive fitted mean. Any further chi-squared goodness-of-fit reference needs its own large-count regularity assumptions.

With the negative-binomial size $\theta >0$ held fixed, collect all mean-independent terms into $c(y_i,\theta )$. The cell log-likelihood isIts derivative isFor $y_i>0$ it is positive below $y_i$ and negative above $y_i$, so the saturated maximum is at ${\mu }_i=y_i$. For $y_i=0$, the maximum is the boundary limit ${\mu }_i\downarrow 0$. Evaluating the saturated-minus-fitted log-likelihood gives the negative binomial devianceAll gamma-function and factorial terms cancel because the same known $\theta$ is used in both models. The zero-count logarithmic term is interpreted as zero; its remaining cell contribution is $2\theta \log (1+e_i/\theta )\ge 0$. No intercept score cancellation is required for this formula, so it applies to the specified design whether or not its covariates include an intercept.

For completeness, the fitted coefficients obey the score equations of the fixed-size negative binomial generalized linear model,The logarithmic link is not the canonical link for this family, but it produces the requested positive means. As a check, when $\theta \to \infty$,so the formula tends to the full Poisson deviance, including the linear terms when an intercept-balance equation is unavailable. The size $\theta$ used here is the negative-binomial size parameter; it is distinct from the random-multiplier variance denoted by the same letter in Question 4.