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

Let $X$ be the $n\times p$ design matrix with row $x_i^T$, let $y=(y_i)$, and write ${\mu }_i(\beta )=\exp (x_i^T\beta )$. Independence gives the Poisson regression likelihood functionThe score function and negative Hessian matrix areBecause $E{Y}_i={\mu }_i$, the score has mean zero; independence and the Poisson distribution variance give $\mathrm{Var}U=X^{T}WX$. Thus the observed and expected Fisher information matrices are the same here.

Assume $X$ has full column matrix rank $p$. For a nonzero vector $a$, $a^{T}Ia={\sum }_i{\mu }_i(x_i^{T}a{)}^2>0$. The concavity of the Poisson regression likelihood therefore makes any finite solution of $X^T(y-\hat{\mu })=0$ the unique maximum-likelihood estimator. If the design is rank deficient, parameters differing by a vector in $\ker X$ give the same distribution, so only identifiable combinations can be estimated. Full rank alone does not ensure a finite maximizer: for an intercept-only model with all observations zero, $\ell =-n{e}^{\beta }$ approaches its supremum as $\beta \to -\infty$. More generally, directions that decrease only fitted means attached to zero counts can produce boundary fits. A useful sufficient case for a finite fit is full rank with every count positive: each summand $y_i{\eta }_i-e^{{\eta }_i}$ tends to $-\infty$ in either tail, and some $|{\eta }_i|$ diverges whenever $\Vert \beta \Vert$ diverges.

For an interior fit, Newton's method and Fisher scoring give the same update:This is also iteratively reweighted least squares. Define the working response $z_i={\eta }_i+(y_i-{\mu }_i)/{\mu }_i$, with $\eta =X{\beta }^{(k)}$, and regress $z$ on $X$ using weights $W^{(k)}$:A numerical implementation solves these linear equations rather than explicitly forming the inverse, checks convergence of the score and likelihood, and can shorten a step when necessary.

Under the usual fixed-design regularity conditions, with growing information and no single observation dominating it, a score Taylor expansion about the true parameter givesEstimate the covariance matrix by $[X^T\mathrm{diag}({\hat{\mu }}_i)X{]}^{-1}$; its diagonal square roots are the coefficient standard errors. This justifies large-sample Wald confidence intervals and Wald tests for specified coefficients or contrasts. Finite-sample unbiasedness is not being asserted. A likelihood-ratio test compares nested fits through twice their maximized log-likelihood difference. The Poisson deviance against the saturated fit iswith $0\log 0=0$. Residual diagnostics and a check of overdispersion assess whether the assumed Poisson variance is adequate; the variance correction required when it is not is derived in the next part.

The working score function and its derivative areAt an interior working maximizer, $U(\hat{\beta })=0$. A Taylor expansion at ${\beta }_0$ yieldsfor an intermediate value ${\beta }^{\ast }$. When the last term is negligible relative to the linear one, this givesThis is the Poisson score linearization under proportional variance. Correct specification of the mean alone implies $EU({\beta }_0)=0$. The stipulated independence and variance relation giveTaking the mean and variance of the linear approximation therefore provesThese are first-order conclusions, not exact finite-sample identities. The mean and variance assumptions alone do not guarantee that the approximation is accurate: an interior consistent solution, increasing information, control of leverage, and a suitable central limit condition are needed for normal inference.

The quasi-likelihood score for variance $\varphi {\mu }_i$ is $U/\varphi$, so its zero is unchanged by the unknown common dispersion parameter. Thus the working Poisson estimating equation remains appropriate for the mean parameter, but its unadjusted Poisson variance estimate is wrong when $\varphi \ne 1$. A usual Pearson dispersion estimator for this one-parameter model isIn particular overdispersion multiplies coefficient standard errors by approximately $\sqrt{\varphi }$, even though it leaves the coefficient estimating equation unchanged.

Treat the counts as independent binomial distributions, $R_i\sim \mathrm{Bin}(n_i,{\pi }_i)$, conditional on the known numbers of operations. The null hypothesis is that the probabilities ${\pi }_i$ are all equal, with no assumed numerical value. Centre 4 has no operations in this period and hence contributes likelihood one, supplying no information about its mortality probability. The empty groups in a binomial homogeneity test principle leaves 11 contributing centres and 10 test degrees of freedom.

Under the common-probability model the maximum-likelihood estimator isUse the two-column contingency table of deaths and survivors, deleting its all-zero row. The Pearson chi-squared test of homogeneity has statisticUnder the null, the large-sample reference distribution is ${\chi }_{10}^2$. The expected counts are adequate for this approximation here, including at the smallest nonempty centre. Calculating givesAn alternative is the likelihood-ratio test against separate centre probabilities, with binomial devianceThe unrestricted probabilities are $r_i/n_i$; subtracting one fitted common probability leaves the same 10 degrees of freedom. Here $D=39.5033$ and $p\approx 2.07\times 10^{-5}$. Thus there is strong evidence against a common mortality probability. This omnibus test does not identify Centre 1 as the sole source of the heterogeneity. If its comparison with the remaining centres is the primary scientific question, specify a separate planned one-versus-rest comparison; inspection of many pairwise comparisons calls for consideration of multiple testing.

Use the logit link $g(\pi )=\log \{\pi /(1-\pi )\}$ and independent counts $R_{ij}\sim \mathrm{Bin}(n_{ij},{\pi }_{ij})$. This gives a grouped-binomial logistic regression with 24 cells and 13 free coefficients: an intercept, 11 centre contrasts, and one period contrast. Reference treatment coding makes the intercept the log odds for Centre 1 in the earlier of the two included periods. The two factors are assumed to have been coded as specified in the question.

The binomial log-likelihood, apart from known constants, isMaximize by Fisher scoring or iteratively reweighted least squares. Using response proportions $r_{ij}/n_{ij}$ with binomial trial weights $n_{ij}$ is equivalent to using the count response $(r_{ij},n_{ij}-r_{ij})$; omitting those weights would give the wrong likelihood and standard errors. The coefficient covariance matrix is approximately the inverse of $X^T\mathrm{diag}\{n_{ij}{\hat{\pi }}_{ij}(1-{\hat{\pi }}_{ij})\}X$.

The fitted reference log odds are $-0.809585$, and the period change is $-0.426758$, givingThese are fitted probabilities under a common period effect, so they need not equal the separate empirical proportions at Centre 1. The estimated period odds ratio isThus the later-period death odds are estimated to be about 35% lower at each centre, with strong evidence of a period effect: the normal Wald statistic is $-5.415$. The displayed “t value” is a coefficient divided by its standard error; with fixed binomial dispersion its usual large-sample calibration is standard normal, not a Student distribution with an estimated Gaussian-error variance. The odds ratio is not a mortality probability ratio.

Each centre coefficient is a log odds ratio relative to Centre 1, adjusted for period. All its point estimates are negative, so all other centres have lower fitted mortality probabilities in either period. For example, Centre 11's odds ratio is $e^{-1.079474}=0.3398$, or about $2.94$ times greater odds at Centre 1 than at Centre 11. Individual unadjusted two-sided Wald tests find significant differences at the 5% level for centres 2, 4, 6, 7, 8, 9, 11 and 12; the coefficients for 3, 5 and 10 do not individually reach that threshold. These separate tests are correlated comparisons against a common reference and should not be presented as 11 independent discoveries without considering multiple testing. They adjust only for time period, not for differing patient case mix.

For fit assessment, the residual binomial deviance is $15.73891$ on $24-13=11$ statistical degrees of freedom. Its ${\chi }_{11}^2$ upper tail is approximately $0.151$, so there is no appreciable evidence of lack of fit to these binomial cell counts. In particular, comparing with a saturated centre-by-period model gives this same deviance and 11 additional parameters; the data do not strongly require a centre-period interaction. The moderate deviance residuals are consistent with that assessment, although a residual plot by centre and period would be more informative than the five-number summary alone.

The reduction from null to fitted deviance isIt tests the 12 centre-and-period coefficients jointly against an intercept-only model, with reference ${\chi }_{12}^2$ and $p\approx 1.24\times 10^{-16}$. It is not a centre-only test adjusted for period. For that test, fit the period-only model and compare its deviance with the full fit: the pooled binomial probabilities in each period give deviance $90.34578$, soThus centre heterogeneity persists after allowing for period. Dispersion one is the stipulated binomial variance, and the five scoring iterations describe numerical convergence rather than statistical evidence. The fitted associations alone do not establish a causal difference in surgical quality.

The first output uses the equal-variance, pooled Student's t-test, as shown by its seven statistical degrees of freedom. Its exact model assumes two independent samples from normal distributions, with means ${\mu }_X,{\mu }_Y$ and a common unknown variance ${\sigma }^2$. The hypotheses are $H_0:{\mu }_X={\mu }_Y$ against $H_1:{\mu }_X<{\mu }_Y$. The sample means are $\bar{x}=4.35$ and $\bar{y}=9.88$, and the unbiased within-sample sample variances areThe pooled sample variance and its resulting standard error areHenceTo justify its reference distribution, under the null $(\bar{X}-\bar{Y})/(\sigma \sqrt{1/4+1/5})$ is standard normal, while $(3s_X^2+4s_Y^2)/{\sigma }^2$ is ${\chi }_7^2$, independently of the two sample means. Dividing the first variable by the square root of the second divided by seven gives the exact pooled two-sample t statistic, with Student t-distribution $t_7$. Its one-sided p-value is the lower tailThe associated one-sided 95% confidence interval for ${\mu }_X-{\mu }_Y$ isThe printed missing lower bound means that this is a one-sided upper confidence limit, not that a finite lower endpoint should be estimated. With such small samples the normality and common-variance assumptions matter; a Welch analysis would use a different variance estimate and generally different degrees of freedom.

The second output is the Wilcoxon rank sum test for two independent continuous samples. Under a common-distribution null hypothesis, labels are exchangeable conditional on the pooled values, so the four ranks assigned to the first sample form a uniformly chosen four-element subset of $\{1,\dots ,9\}$. The first sample receives ranks $2,1,3,5$, giving raw rank sumEquivalently the Mann–Whitney U test statistic is $U_X=W_X-4\cdot 5/2=1$, the number of cross-sample pairs with an $X$ observation larger than a $Y$ observation. This output reports the raw rank sum, rather than the shifted statistic used by some software.

A first sample tending to be smaller has a small rank sum. Under the null, all $\left(\genfrac{}{}{0ex}{}{9}{4}\right)=126$ rank allocations have equal probability. Only $\{1,2,3,4\}$ and $\{1,2,3,5\}$ have rank sum at most 11. By the exact rank-sum tail one above the minimum, the exact lower-tail p-value isNo large-sample normal approximation or continuity correction is involved; the absence of ties makes this simple enumeration valid. For a common-shape location-shift family, the test's alternative is a negative location shift of $X$ relative to $Y$. Without that model, the null is equality of distributions and the rank comparison is not simply a test of equality of means or medians.

Both tests reject at the 5% level in the direction of smaller $X$ observations; only the pooled t-test rejects at 1%. The t-test uses the numerical spacing and targets a mean difference under its normal model. The rank-sum test uses ordering and an exact exchangeability calculation, avoiding the normality and common-normal-variance assumptions. Their different p-values therefore reflect both different statistics and different modelling assumptions, not a numerical contradiction.

The commands read the 52 cell values into a vector and build the matching two-factor layout. In expand.grid(country, goods) the first factor varies fastest, so the four countries are listed successively within each goods category. This matches row-by-row reading of the displayed table. Converting Goods to a factor is essential: its labels identify 13 categories and are not a quantitative trend. Country supplies the other categorical factor. The summaries suppress coefficient correlations, not coefficient standard errors.

The first fit is the additive two-factor normal linear modelThere are $1+12+3=16$ free coefficients and $52-16=36$ residual statistical degrees of freedom. This is two-way analysis of variance without an interaction term. For this complete balanced layout, ordinary least squares givesThese formulas follow from the row- and column-sum normal equations: fitted row and column totals equal their observed totals. The balanced-design coefficient standard errors are $s\sqrt{2/4}$ for a goods contrast, $s\sqrt{2/13}$ for a country contrast, and $s\sqrt{1/4+1/13-1/52}$ for the reference-cell intercept. These yield the common contrast standard errors and the different intercept standard error in each summary. In particular the intercept is the fitted Software–UK value, not the grand mean or the observed Software–UK entry.

Thus the first intercept predicts $149.8269$ million euros for Software in the UK. The Books coefficient adds $254.7$ million euros to the fitted value in any country, relative to Software; every other goods coefficient has this same arithmetic-difference interpretation. The country coefficients subtract $193.7$, $57.0308$ and $97.3769$ million euros from the UK prediction for France, Germany and the rest of Europe respectively, regardless of goods. The model therefore assumes common absolute country differences across goods categories. For Software–France its prediction iswhich is impossible for positive spending and indicates a limitation of this additive original-scale specification.

The residual standard error is $80.98$ million euros, obtained from $s^2=\mathrm{RSS}/36$; direct calculation gives $\mathrm{RSS}=236088.7435$. Each coefficient's printed statistic is $\hat{\theta }/\mathrm{se}(\hat{\theta })$, with a $t_{36}$ reference distribution under this normal linear model. The two-sided p-values compare each goods category with Software, or each other country with the UK, after adjustment for the other factor. For example France and the rest of Europe differ significantly from the UK at 5%, but Germany does not. The raw-scale significant goods differences are Books, Clothing, Electronic Goods, Groceries and Leisure Travel; these are individual reference-level comparisons, not an overall test of the goods factor.

The coefficient of determination is $R^2=1-\mathrm{RSS}/\mathrm{TSS}=0.7743$. The omnibus F-test compares all 15 nonintercept coefficients with zero:with reference $F_{15,36}$ and $p=1.2544\times 10^{-7}$. This says at least one factor contrast is nonzero under the assumed error model. There is one observation per cell, so an unrestricted goods-country interaction would saturate the model, leaving zero residual degrees of freedom. Without replication or a specified interaction structure, the residual variation cannot be separated into interaction and independent error.

The second fit applies a logarithmic transformation and fits the same factor design to $Z_{gc}=\log Y_{gc}$. The Gaussian log-response model isIts intercept exponentiates to a fitted Software–UK conditional median of $e^{4.5267}\approx 92.45$ million euros. Goods coefficients now give multiplicative ratios to Software, and country coefficients give multiplicative ratios to the UK. The estimated country ratios areFor example the Books ratio is $e^{1.6415}\approx 5.16$ and the Sports Equipment ratio is $e^{-1.8364}\approx 0.1594$, consistently across countries. The model assumes common relative, rather than common absolute, country effects. It also guarantees positive retransformed fitted values. Under normal log errors, the original-scale conditional mean is $e^{\eta +{\sigma }^2/2}$, whereas $e^{\eta }$ is the conditional median and geometric mean. Correcting the retransformation bias here multiplies the naive predictions by approximately $e^{0.127406/2}=1.0658$.

The log-scale residual standard error is $0.3569$, with the same 36 residual degrees of freedom, corresponding to a typical multiplicative error multiplier $e^{0.3569}\approx 1.429$. France and the rest remain significant reference-country contrasts; Germany remains nonsignificant at 5%. The log-scale significant reference-goods contrasts are Books, Music, Videos/DVDs, Clothing, Sports Equipment, Electronic Goods, Groceries and Leisure Travel. Event Tickets, Toys, Video Games and Housewares do not individually differ significantly from Software at that level. These p-values remain subject to the usual multiple testing qualification.

The log-scale $R^2=0.9381$ is computed using variation in the logarithms, so it cannot by itself establish better fit on the original response scale. The change is plausible because spending varies widely in magnitude and relative discrepancies are more meaningful than constant absolute discrepancies, but it should be assessed with regression diagnostics, especially a residual-versus-fitted plot and a normal Q-Q plot. The minimum residual near $-1.01$ merits inspection; the command summaries alone do not verify independent normal errors. Moreover these responses are forecasts, so the inferential p-values have a literal sampling interpretation only if an appropriate stochastic model for their uncertainty is justified.

The log-scale analysis of variance adds Goods and then Country. In this balanced design their effect spaces are orthogonal after removing the intercept, so the sequential sums of squares equal the corresponding adjusted factor sums and are unchanged by reversing the two factors:The residual sum of squares is $4.58662$, giving $M{S}_E=4.58662/36=0.127406$. ThusTheir tiny upper-tail p-values show strong overall goods and country effects under the log-error model; a globally significant country factor is compatible with the nonsignificant individual Germany contrast. The omnibus fit statistic is $36.3433$ on $(15,36)$ degrees of freedom, whose p-value is approximately $3.14\times 10^{-17}$. The printed zero reflects numerical formatting, not an exactly zero probability.

If two cell entries were missing, retain their correct factor labels and mark those responses unavailable, rather than shifting the remaining scan values into different cells. Fit the additive models to the 50 observed cells, assuming the missingness mechanism permits those observations to be analyzed without selection bias. All factor levels remain represented, and the row-column incidence graph remains connected: any two goods rows still share an observed country, since only two edges of the complete layout have been removed, and every country still has many observed goods. By identifiability of an incomplete two-factor additive design, the design therefore still has rank $13+4-1=16$. HenceThe coefficients and their covariance must be recomputed by the observed-cell ordinary least squares equations. The balanced row-mean-plus-column-mean formulas and common within-factor standard errors no longer apply. Goods and Country effects are generally no longer orthogonal, so sequential ANOVA sums can depend on order; factor tests adjusted for the other factor should use nested model comparisons. Missing cells do not restore replication or make an unrestricted interaction test possible: that larger model still fits every observed cell exactly. If missingness depends on the unavailable spending values, additional modelling of that mechanism would be needed.

Use the time-homogeneous continuous-time multi-state model displayed in part ii. Write $Z(t)$ for the state at time $t$, $Q=(q_{rs})$ for its generator matrix, andFor $r\ne s$, $q_{rs}={\lim }_{h\downarrow 0}{p}_{rs}(h)/h$ is the transition intensity, and $q_{rr}=-{\sum }_{s\ne r}{q}_{rs}$. Condition on the observed initial state, and treat the clinic observation process and administrative end of follow-up as noninformative. A panel-observed interval contributes its transition probability, not the density of a jump at the examination date; the transitions inside that interval need not have been observed.

For the first patient, the two clinic intervals contribute $p_{11}(2)p_{13}(2)$. The subsequent death is observed at elapsed time $u=0.2$ after the last clinic state, so its contribution is the conditional death-time densityThe equality uses the permitted directions: starting from state 3, only living states 3 and 4 can be occupied before death. To obtain this density directly, survival in a living state $s$ until time $u$, followed by death in $(u,u+du)$, has probability $p_{3s}(u)q_{s5}du+o(du)$; add over the unobserved state $s$. Hence the mixed panel and exact-death likelihood contribution isIt includes the possibility of an unrecorded transition $3\to 4\to 5$ between the last clinic visit and death. Keeping only $p_{33}(0.2)q_{35}$ would exclude that allowed history, while $p_{35}(0.2)$ would be the probability of death by that time rather than the density of death at that time. With dates recorded in whole days, an exact interval-likelihood treatment integrates this density over the relevant day; the density formulation is its fine-time approximation.

For the second patient, the three successive intervals with state 1 at both ends, then the interval ending in state 2, and the last interval remaining in state 2 give the panel-observed multi-state likelihoodThere is no additional death contribution after the final administrative observation. Each likelihood contribution already accounts for being alive in the observed state at that visit. For independent patients their joint conditional likelihood contribution is $L_1{L}_2$.

For an explicit form, put ${\lambda }_r=-q_{rr}$. Irreversibility givesandThese formulas follow by integrating the density of the one possible intermediate jump time; at equal exit rates the limiting expression is $q_{rs}u{e}^{-{\lambda }_{r}u}$. ThusUsing the displayed fitted rates as numerical values gives $g_{35}(0.2)\approx 0.289314$ per year, $L_1\approx 0.017596$ per year and $L_2\approx 0.016048$, providing a check of the probability-versus-density distinction.

The permitted transitions form the following diagram. Each arrow is labelled by its fitted transition intensity in units of inverse years; the absence of an arrow is a structural restriction, not an estimated small effect.

State 1 can move directly to 2, 3 or 5; state 2 can move to 4 or 5; state 3 can move to 4 or 5; and state 4 can move only to 5. There are no recovery transitions, and state 5 is an absorbing state. The rates of leaving the four living states are the sums of their outgoing arrows:The diagonal generator entries are their negatives, so every row of $Q$ sums to zero. These diagonal entries are not arrows or negative event rates. For example the direct $1\to 5$ rate is $0.004$ per year, while death after reaching another state can occur at a much larger rate. The printed fit has no covariates on the intensities, so it estimates one pooled rate system and contains no treatment-effect comparison, despite the trial's treatment randomization. The reported $-2\ell$ describes the fitted likelihood level; by itself it is not a goodness-of-fit test without a reference or competing model.

For a homogeneous continuous-time Markov chain, the holding time in a living state $r$ before the next transition has survival probability $e^{-{\lambda }_{r}t}$, where ${\lambda }_r={\sum }_{s\ne r}{q}_{rs}=-q_{rr}$. Indeed, conditional on still occupying that state, the probability of leaving in the next $dt$ is ${\lambda }_{r}dt+o(dt)$, so its survival function satisfies $S^{\prime }(t)=-{\lambda }_{r}S(t)$ and $S(0)=1$. Integrating gives the mean holding time from a transition intensity matrix:This is the mean time until any allowed exit from that state, not necessarily the patient's remaining lifetime.

For confidence limits, use the supplied standard error $s_r$ of the corresponding diagonal generator estimate. Changing its sign does not change its standard error, so it is also the standard error of ${\hat{\lambda }}_r$. The delta method, with derivative $d(1/\lambda )/d\lambda =-1/{\lambda }^2$, givesThe resulting approximate 95% Wald confidence intervals areThe diagonal standard errors already incorporate the covariance among outgoing rates. They should not be replaced by the square root of a sum of squared off-diagonal standard errors, which would silently assume those rate estimates independent.

Alternatively construct a positive rate confidence interval and apply the confidence interval for an inverse rate transformation. Here every normal rate interval has a positive lower endpoint, givingThese approximate intervals are, respectively, $(10.863,12.830)$, $(4.497,7.161)$, $(2.397,3.117)$ and $(1.548,2.631)$ years. They are asymmetric because inversion is nonlinear; positive log-rate or likelihood-based intervals could also be inverted. The symmetric delta intervals and the inverted-rate intervals are two legitimate first-order approaches and need not have numerically identical endpoints.

The death state has zero exit rate, so its holding time is infinite. There is no finite mean-sojourn estimate or ordinary finite confidence interval for this structurally absorbing state. All intervals above concern the mean parameter, not a prediction interval for one patient's realized holding time.

The transition matrix over five years is $P(5)=e^{5Q}$. Its first row conditions on starting in state 1. The displayed approximate probabilities mean that, five years later, about 65% remain disease-free, 6% have local recurrence only, 10% have distant metastasis only, 2% have both, and 16% have died. They are probabilities of the states occupied at five years; they are not the probabilities of the next jump destinations or annual transition intensities. In particular the positive probability of state 4 does not require a direct $1\to 4$ arrow: the paths $1\to 2\to 4$ and $1\to 3\to 4$ contribute to it. The probability of death includes deaths after any allowed sequence of intermediate states.

The printed row sums to $0.99$ rather than one because the entries are approximate; an exact stochastic transition row sums to one. Exponentiating the displayed three-decimal generator gives approximately $(0.65377,0.06490,0.09775,0.02532,0.15826)$, illustrating that the edited, low-precision table should not be treated as an exact matrix exponential of its separately rounded rates.

State 5 is an absorbing state, so once entered it can never be left. Its generator row is zero, and the last row of every positive-time transition matrix is thereforeThis explains all four leading zeroes in the printed last row.