Past exam of the mathematics course of the University of Cambridge / 2006 / iii / Paper 44 (Biostatistics)

Use a longitudinal study with repeated measurements of the same individuals. Recruit people across several starting ages and measure the same ear, with a standardized anatomical definition and instrument, at baseline and subsequent visits over several years. Repeat measurements at each visit to estimate measurement error; where practical, conceal previous measurements from the observer. Record sex, body size, birth cohort and observer, and document reasons for loss to follow-up.

A simple random-intercept linear mixed model for ear length $Y_{ij}$ at elapsed follow-up time $t_{ij}$ isThe individual intercept ${\alpha }_i$ absorbs persistent differences between people. Alternatively, differences $Y_{ij}-Y_{i0}$ remove that intercept directly. Estimate the within-person regression coefficient $\beta$ and its confidence interval, accounting for the correlation of repeated measurements. Compare it with the cross-sectional slope rather than treating the two slopes as automatically identical.

This longitudinal study directly checks whether an individual's ears tend to lengthen as time passes. Several entry-age cohorts help assess whether the change varies with starting age. Calendar-period changes and selective dropout can still affect interpretation, so repeated observation alone does not eliminate every possible source of confounding or selection bias.

The association is strongly statistically significant under the reported regression assumptions, but age alone has only moderate predictive ability. The reported confidence interval for the age regression coefficient excludes zero by a wide margin. If it is the usual two-sided 95% interval from a normal linear model with an intercept and one predictor, there are $206-2=204$ degrees of freedom. Taking $t_{0.975,204}\simeq 1.97$ givesThus the two-sided p-value against zero slope is far below $0.001$. This conclusion concerns a statistical association in the sampled population; it does not establish an individual growth rate. Clustering by practitioner, nonconstant variance or nonindependent observations would require an appropriate standard error rather than blind reliance on this calculation.

The scatter plot has substantial vertical spread at each age. A narrow slope confidence interval describes uncertainty about the average trend, whereas a prediction interval for a new person's ear length must also include the residual variance. The reported interval even permits an approximate numerical assessment: write $S_{xx}={\sum }_i(x_i-\bar{x}{)}^2$ and let $s_e^2=\mathrm{RSS}/(n-2)$. In simple linear regression, $t^2={\hat{\beta }}^2{S}_{xx}/s_e^2$, while the fitted and residual sums of squares are ${\hat{\beta }}^2{S}_{xx}$ and $(n-2)s_e^2$. Hence the coefficient of determination isUnder these assumptions, age explains only about 27% of the sample variation; roughly 73% remains unexplained. The calculation is approximate because the printed slope and interval are rounded. For a new individual aged $x_0$, the usual prediction interval is based onThe leading one represents individual residual variation. The provided summaries do not determine $s_e$ and $S_{xx}$ separately, so they do not supply a numerical prediction interval.

Record sex, height or another measure of body size, ancestry, family characteristics and practitioner or measurement method. Use these covariates in a multiple linear regression, allowing scientifically justified nonlinear age effects or interactions, and assess prediction error on held-out observations. A guessed overall sex proportion cannot replace individual sex measurements or identify sex-adjusted effects.

There is a source inconsistency: the PDF's numerical ear-length summary is ten times the scale of its plotted measurements and fitted line. Indeed, the fitted line at the reported mean age gives about $67.7$ mm, whereas the printed mean is about ten times larger. Missing decimal points in the summary are a plausible explanation. The significance and approximate $R^2$ above use the mutually consistent slope, interval and graph; they do not treat the inconsistent summary as valid data.

A cross-sectional age association need not equal within-person growth. The following alternatives each suggest a concrete extension of the study.

These explanations are hypotheses to be checked, not established facts about the participants. Neither a small p-value nor an approximately linear scatter plot excludes them.

Assume the incubation times are independent observations from the same upper-truncated distribution. Conditional on lying below $M$, the gamma distribution density is $g(t;a,s)/G(M;a,s)$ on its support. With $t_{(n)}={\max }_i{t}_i$, the likelihood function isIn particular, the normalizing factor $G(M;a,s{)}^{-n}$ cannot be omitted. On the admissible region the log-likelihood isThe profile likelihood eliminates the nuisance parameters $a,s$ by fitting them separately at each candidate endpoint:For fixed $a,s$ and $M>t_{(n)}$,Equivalently, for $t_{(n)}<M_1<M_2$, every nuisance-parameter fit has $L(a,s,M_1)>L(a,s,M_2)$. Taking suprema shows that $L_p$ is nonincreasing. If the nuisance maximum at $M_2$ is attained and finite, evaluating its maximizing parameters at $M_1$ shows strict decrease there as well. Values below the largest observation are inadmissible. Thus upper-endpoint estimation in a truncated distribution gives the usual boundary estimateThere is a support-convention qualification: the literal strict support $t<M$ excludes $M=t_{(n)}$, so it gives a supremum approached as $M\downarrow t_{(n)}$, rather than an attained maximum-likelihood estimator. Using the equivalent closed-endpoint density convention for a continuous distribution yields the usual boundary estimate above. The argument presumes a meaningful finite nuisance fit; degenerate samples can instead have an unbounded density likelihood. Since $M$ changes the support, regular interior-parameter asymptotic formulas are not automatically applicable.

For the log-transformed location-scale model, let $F_W$ denote the distribution function associated with $f$. The change-of-variables formula for a probability density gives the untruncated density and distribution functionConsequently its upper-truncated likelihood is ${\prod }_{i}p(t_i)/F(M{)}^n$ on the observed support, and the profile likelihood maximizes this over $\alpha ,\sigma ,q$. Restricting to the corresponding subfamilies gives the gamma distribution and lognormal distribution profiles. Write the plotted likelihood ratio asIt measures how well each endpoint fits relative to that model's own best endpoint, after nuisance refitting. It is not a probability density function for $M$ and need not integrate to one. Endpoints whose ratios are near one are relatively compatible with the observations; a rapidly falling profile indicates more information favoring endpoints near the maximum.

All three profiles favor the smallest admissible endpoint, approximately 14 days; their information about larger endpoints differs substantially. Reading the actual PDF legend and curves, the dotted gamma distribution profile falls initially but is nearly flat around $0.60$ from about 18 days onward. The solid generalized log-gamma profile declines much further, to roughly $0.20$ at 19 days and $0.18$ at 21 days. The dashed lognormal distribution profile falls most strongly, to roughly $0.14$ at 18 days and $0.08$ at 21 days. These are graphical readings, not recomputed fits.

Thus the gamma distribution fit gives weak discrimination between moderate and very large endpoints: much of its untruncated mass can already be below these endpoints, making further changes to the truncation have little effect. Over the displayed range, the lognormal distribution fit gives stronger relative evidence against large endpoints, and the more flexible generalized log-gamma fit is intermediate. The solid curve remains close to the dashed curve near the maximum but has appreciably more relative support farther out. The plot shows slowing declines rather than a sharply determined finite upper endpoint; it cannot establish a biological maximum.

Normalization is essential to this comparison. Although the generalized family includes the two special families and therefore has at least as large a raw maximized likelihood at any fixed $M$, its normalization constant can also be larger. Its normalized profile need not lie above the normalized profiles of its submodels. Their individual maxima have all been replaced by one, so this plot alone cannot compare overall goodness of fit or select a family while accounting for extra fitted parameters.

A set $\{M:R(M)\ge c\}$ is a relative-likelihood support set once $c$ is chosen. Turning it into a calibrated confidence interval requires an appropriate sampling argument, simulation or other justified method. In particular, a usual Wilks theorem chi-squared cutoff cannot simply be assumed for this changing-support endpoint problem. The original incubation observations are not supplied here, so numerical refitting or calibrated limits cannot be recovered from the graph alone.

Quarantine restricts contacts of exposed people whose infection status is uncertain and monitors them for illness, reducing opportunities for transmission before infection is recognized. Recognized cases are managed by isolation rather than treated as merely uncertain exposures. Duration must be measured from the last possible infection-producing exposure, or otherwise account for uncertainty in exposure time. The mathematical calculations below concern the incubation period among infected people who develop symptoms, not a complete infectiousness model.

With a credible hard endpoint $M$, the truncated distribution hasA duration $q\ge M$ then covers all symptom-onset times under that model, provided the exposure time is correctly specified. If a nonzero late-onset probability $\epsilon$ is accepted, solveThe calculation is appropriate when external biological evidence genuinely supports a hard maximum, or when a deliberately bounded model is justified for the population being considered. Truncation caused by a study's ascertainment window instead requires an observation correction: it does not prove that the disease has a hard biological endpoint. Nor does $\hat{M}=t_{(n)}$ guarantee that future cases will develop symptoms by the current sample maximum. Endpoint and nuisance-parameter uncertainty need to enter a conservative duration calculation.

For an untruncated gamma distribution, lognormal distribution or other model with an unbounded upper tail, no finite duration eliminates the modeled late-onset probability. Choose a tail tolerance and solveThis is appropriate when a hard endpoint is unsupported and rare long incubation times remain possible. One can estimate a high quantile, assess its uncertainty and balance the modeled residual probability against the costs of longer restrictions. Different tail families can give very different high quantiles even when they fit the central observations similarly, which makes the model dependence seen in part b operationally relevant.

If an individual has already remained symptom-free for $a$ time units after a known exposure, a conditional additional duration $b$ instead uses $S(a+b)/S(a)\le \epsilon$, where $S$ is the relevant survivor function and $S(a)>0$. This conditions on being an infected eventual symptomatic case. Incubation data alone do not describe asymptomatic infection, the timing of infectiousness or compliance with restrictions; symptom surveillance, testing where informative and a transmission model may therefore be needed to translate late onset into an actual release-risk calculation. Use a hard bound only when it is supported; otherwise quantify and allow for the upper-tail uncertainty.

A frailty model represents unobserved individual heterogeneity by a random effect in the hazard function. In a proportional frailty model, an individual with frailty random variable $U=u$ has hazard $u{h}_0(t)$. The normalization $\mathbb{E}U=1$ identifies the scale of the baseline hazard: without it, multiplying all frailties by a constant and dividing $h_0$ by that constant would leave all individual hazards unchanged.

Let $H_0(t)={\int }_0^t{h}_0(v)dv$. Conditional survival follows by integrating the hazard function:Taking the expectation over the frailty law yieldsThis is the Laplace transform of the frailty density evaluated at the cumulative hazard. More generally, the same calculation integrates against a probability measure, allowing atoms. It averages conditional survival rather than substituting $\mathbb{E}U$ into its exponent. Indeed, the frailty distribution among survivors is proportional to $e^{-u{H}_0(t)}g(u)$, giving population hazardSurvival progressively selects smaller frailties, so the population hazard need not be the baseline or even retain the individual's proportional structure.

For the cure model, first take $0\le \pi <1$. The necessary mean-one frailty law and baseline hazard areThen $\mathbb{E}U=\pi \cdot 0+(1-\pi )/(1-\pi )=1$. A cured person's conditional hazard is zero, while a susceptible person's conditional hazard is $h_0/(1-\pi )=h_{\ast }$. Taking $U$ to be zero or one while leaving $h_0=h_{\ast }$ would instead have mean frailty $1-\pi$ and violate the required normalization.

The frailty law is the measurewhere the $\delta$ terms are Dirac measures. This is not an ordinary Lebesgue density; the cure construction requires permitting a discrete frailty law. Its Laplace transform isWriting $H_{\ast }(t)={\int }_0^t{h}_{\ast }(v)dv$, we have $H_0=(1-\pi )H_{\ast }$, and thereforeIf $H_{\ast }(t)\to \infty$, the susceptible class eventually experiences the event and $S(t)\to \pi$. If $H_{\ast }(\infty )<\infty$, the limiting survival is larger than $\pi$ because some susceptible subjects also never experience the event. At $\pi =1$, everyone is cured: set $h_0\equiv 0$ and, for example, $U\equiv 1$ to retain mean one, giving $S\equiv 1$ without using the singular formula $1/(1-\pi )$.

The log-rank test compares event counts with their conditional expectations under equal group hazard functions, using the individuals currently in each risk set. Assume the usual independent censoring conditions and independent strata. At pooled event time $t_j$, let $n_{Aj},n_{Bj}$ be the group numbers at risk, $n_j=n_{Aj}+n_{Bj}$, and let $d_{Aj},d_{Bj}$ be the event counts with total $d_j$. Under the null hypothesis, conditioning on the risk counts and total events gives the hypergeometric distributionThus the null expected A count and its conditional variance areWhen there is only one person at risk, the event allocation is deterministic and its variance is zero. The unstandardized log-rank statistic is $U={\sum }_j(d_{Aj}-E_j)$, with variance estimate $V={\sum }_j{V}_j$. Under standard large-sample conditions, $U/\sqrt{V}$ is approximately standard normal under the null; the test detects systematic group differences in event rates. It is particularly effective for proportional-hazards alternatives, while crossing effects can cancel in the score.

For a stratified log-rank statistic, construct these risk sets and observed-minus-expected differences separately within each stratum $s$, then sum:This allows different baseline hazard functions between strata and avoids comparing an event with individuals in another stratum. The null is equal A and B hazards within each stratum, together with a censoring mechanism that does not invalidate the within-risk-set comparison.

In a paired stratum, both individuals initially form the risk set. A first event while both are still observed gives one nontrivial event allocation. After either individual leaves, every subsequent event is certain to belong to the only remaining group and has observed count equal to expected count. Therefore each pair supplies at most one informative table. A pair supplies none if the first exit is censoring before either observed event, if neither has an event during their joint follow-up, or if both fail simultaneously so that the tied event allocation is fixed. An event-censoring tie needs a specified convention; with events processed before same-time censoring, the other individual remains at risk for that event and the pair is informative. The untied continuous-time case avoids this ambiguity.

If A has the informative first event, the table at that time isThe phrase in the second column refers only to this event time; it does not mean that B can never fail. The expected A count is $1\cdot 1/2=1/2$, so this table contributes $1-1/2=1/2$ to $U$, with variance $1/4$. An informative B-first pair contributes $-1/2$. All other tables contribute zero. Therefore the requested unnormalized score isA positive score means more observed early events on A. There is no need to standardize it to answer the question.

To obtain an exact conditional law, put $m=d_A+d_B$. Under the within-pair equal-hazard null, conditional on an event while both are at risk, either subject supplies it with probability $1/2$. Equivalently, before the pair's first exit, the cause-specific event intensities for A and B are equal; integration over the common period of joint observation gives equal probabilities for the two informative outcomes. This argument requires noninformative observation of the failure process, not just equality of unadjusted marginal survival curves. Different independent pairs then have independent fair signs when informative. Conditioning on the set of informative pairs gives $m$ independent Bernoulli outcomes with probability $1/2$; conditioning only on their count mixes identical binomial laws and gives the same result:This proves the paired log-rank reduction to a sign test. The conditional mass function is $\left(\genfrac{}{}{0ex}{}{m}{k}\right)2^{-m}$, $0\le k\le m$, so an exact two-sided p-value sums the masses for $k$ with $|2k-m|\ge |d_A-d_B|$. If $m=0$, there is no information for a treatment comparison. Failure-time distances beyond their ordering play no role in this paired score.

Use the usual semiparametric Cox proportional-hazards model with observed covariates $z_i$:The covariates are implicit in the question's regression fit, even though the displayed data notation only lists exit times and event indicators. Assume entry at a common origin, time-constant covariates, distinct observed times and independent censoring. Write $R_i=\{j:x_j\ge x_i\}$ and $w_j(\beta )=e^{{\beta }^T{z}_j}$. Conditional on a single event at $x_i$ in its risk set, the chance that individual $i$ supplies it isThe baseline hazard cancels. Consequently estimate $\beta$ by maximizing Cox partial likelihoodFor example, the score and information of Cox partial likelihood follow by differentiation. Define $S_i^{(0)}={\sum }_{j\in R_i}{w}_j$, $S_i^{(1)}={\sum }_{j\in R_i}{w}_j{z}_j$ and $S_i^{(2)}={\sum }_{j\in R_i}{w}_j{z}_j{z}_j^T$. ThenWhen a finite identifiable maximum exists, iterate $\beta \leftarrow \beta +I(\beta {)}^{-1}U(\beta )$ until convergence. The observed Fisher information also gives the usual asymptotic coefficient uncertainty. Censored subjects enter the relevant risk sets, even though they supply no numerator event factor.

After fitting $\hat{\beta }$, estimate the cumulative hazard of the reference subject by the Breslow estimatorIts jump at an event is the observed event count divided by the total fitted hazard multiplier at risk. The estimate is an integrated baseline hazard, with hazard represented by a measure placing these jumps at observed event times. If a smooth ordinary hazard function is wanted, smooth these increments with a chosen kernel and bandwidth; differentiating the step function as an ordinary function would incorrectly give zero away from the jumps. The fitted individual cumulative hazard is ${\hat{H}}_i(t)=e^{{\hat{\beta }}^T{z}_i}{\hat{H}}_0(t)$ over that individual's follow-up. The invariances below concern this Cox partial likelihood fit and this Breslow estimator; they need not hold for a fully parametric proportional-hazards fit or an arbitrarily smoothed hazard estimator.

Since the observed times are ordered, the final risk set is $R_m=\{m\}$. Its Cox partial likelihood factor is thereforefor either value of the last event indicator. Every earlier risk set and earlier event factor is unaffected by changing $v_m$. Thus the entire partial likelihood as a function of $\beta$ is unchanged, not merely its value at one fitted coefficient. Its maximizing set, its score function and its observed Fisher information are identical. In particular, any consistently chosen finite fitted $\hat{\beta }$ is unchanged. The lone terminal event conveys no relative-hazard information because there is no other individual with whom to compare it.

Use the right-continuous Breslow estimator, so the increase from the value at $x_{m-1}$ to the value at $x_m$ means the interval $(x_{m-1},x_m]$. There is no event between these consecutive exit times, and at $x_m$ the risk set contains only subject $m$. Henceand the fitted individual cumulative hazard increases byIt is the subject's integrated hazard, rather than the reference baseline increment, that equals the event indicator.

For the final unheaded request, a fitted martingale residual at a subject's exit isThe revised observation leaves the last subject at risk until a later event time, still after every other subject has exited. Every earlier risk set is exactly the same as before, and the new final event again contributes the constant factor one. Thus part a gives unchanged coefficients, and every earlier Breslow estimator increment is unchanged. No additional events occur during the extended interval because all other subjects have already left observation.

For $i<m$, the event indicator and fitted cumulative hazard at $x_i$ remain unchanged, hence so does ${\hat{M}}_i$. For the terminal subject, let $C$ denote its old fitted integrated hazard at its censored exit. That old residual is $-C$. Part b shows that the new final event adds one to its fitted integrated hazard, giving new residualTherefore all fitted martingale residuals remain unchanged. Their invariance does not assert that the fitted survivor function at later times or the terminal event record is unchanged: the new late baseline hazard jump is real. This is the terminal-observation invariance of Cox martingale residuals, which relies on a unique final subject and unchanged covariate history.