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

The paper uses $F$ for the survivor function, rather than for the cumulative distribution function. For a nonnegative survival time $T$ with an absolutely continuous probability distribution, defineWhere $F(t)>0$, the hazard function is the instantaneous event rate conditional on survival:The integrated hazard is $H(t)={\int }_0^{t}h(u)du$. Since $h=-F^{\prime }/F$, integration gives $H(t)=\log F(0)-\log F(t)$. With the usual time origin and no event at zero, $F(0)=1$, soConversely, differentiating $F=e^{-H}$ gives $F^{\prime }=-hF$, so these relationships determine each of the survivor function, hazard function and integrated hazard from either of the others. A positive limiting survivor function is possible if the integrated hazard remains finite; a proper finite event time instead has $F(t)\to 0$.

Assume independent individuals and independent censoring: conditional on any modelled covariates, the censoring mechanism gives no additional information about the latent event time and has no parameter shared with its survival distribution. An observed event at $x_i$ contributes the probability density function $f(x_i)$. A right-censored observation contributes the survivor function $F(x_i)$, since it establishes only $T_i>x_i$.

More formally, let $C_i$ be an independent censoring time with density $c_i$ and survivor function $G_i$. The observed pair $X_i=\min (T_i,C_i)$, $v_i=1_{\{T_i\le C_i\}}$ has a parameter-dependent factor $f(x_i)G_i(x_i)$ if $v_i=1$ and $F(x_i)c_i(x_i)$ if $v_i=0$. The factors involving only the censoring distribution can be omitted from the event-time likelihood. ThusUsing $f=hF$ also gives the useful survival likelihood formEvery individual therefore contributes survival exposure, whether or not an event is observed. These formulas do not justify ignoring informative censoring; that situation requires a model for the observation mechanism as well.

For an exponential distribution of rate $\theta >0$, the probability density function and survivor function are $f(t)=\theta e^{-\theta t}$ and $F(t)=e^{-\theta t}$. Setwhere $D$ is the observed event count and $E$ is total person-time at risk. The log-likelihood, up to a constant, becomesWhen $D>0$ and $E>0$, its score function is $D/\theta -E$, which vanishes atIts second derivative is negative, so this is the unique maximum-likelihood estimator. This is exponential-rate estimation from censored exposure: the estimate is observed events divided by all observed exposure, including exposure from right-censored individuals.

If $D=0$ and $E>0$, $\ell =-\theta E$ decreases strictly. There is no maximizer in the open parameter space $\theta >0$; its supremum occurs as $\theta \downarrow 0$, or the extended boundary estimate is $\hat{\theta }=0$ if a zero rate is permitted. The formula should not be interpreted as an interior regular fit in this case.

Differentiating the log-likelihood again givesThus the observed Fisher information for the rate is proportional to the total number of observed events. At a positive interior maximum-likelihood estimator, the large-sample variance estimate and standard error arePrecision is therefore driven by the number of observed events rather than by the number recruited alone. Heavy right censoring can leave a large study with little information about its hazard function. Under the regularity conditions for the Fisher information, the expected information is $\mathbb{E}[D]/{\theta }^2$.

A right-censored observation is not wholly uninformative: its exposure contributes to $E$, changes the maximum-likelihood estimator, and affects the expected number of events under a sampling design. Its contribution $-\theta x_i$ happens to have zero second derivative in this rate parameterization. For example, if no events occur, the likelihood still declines as $e^{-\theta E}$ and disfavors large rates, despite its zero curvature. The usual interior variance approximation is then invalid.

In a frailty model, an individual's unobserved susceptibility to an event is represented by a nonnegative frailty random variable multiplying that individual's hazard function. Conditional on its value, subjects with the same measured covariates can therefore have different event rates. The frailty is a persistent latent characteristic of the individual, not newly sampled noise at each time.

Here the conditional hazard function is $u{e}^{\beta z}{h}_0(t)$, with baseline hazard $h_0$. A multiplier $u>1$ raises susceptibility and $u<1$ lowers it. Normalization $\mathbb{E}U=1$ fixes the frailty scale: without it, multiplying every frailty by a constant and dividing the baseline hazard by that constant would give the same individual hazards. The variance of $U$ measures unobserved heterogeneity. An atom at $U=0$, if allowed, would give a permanently zero conditional hazard.

The hazard ratio between groups is $e^{\beta }$ for individuals compared at the same frailty. It need not equal the population hazard ratio, because survival selection changes the frailty distribution among those still alive or event-free.

Condition on a frailty random variable value $U=u$. Integrating the conditional hazard function gives $u{e}^{\beta z}{H}_0(t)$, so the conditional survivor function isAverage this conditional probability over the common frailty density $g$:By definition the Laplace transform of $g$ is $\tilde{g}(s)={\int }_0^{\infty }{e}^{-su}g(u)du$. HenceThis averages survivor functions, not hazards. Differentiating its logarithm, when justified, givesThus the population hazard function uses the mean frailty among current survivors, explaining why a conditional proportional hazards model may cease to be proportional after frailty is integrated out.

For ${\sigma }^2>0$, a gamma distribution with unit mean and variance ${\sigma }^2$ has shape and rate $\psi =1/{\sigma }^2$. Substituting its Laplace transform into the preceding result givesIts negative logarithmic derivative is the population hazard function:Consequently, at times with $h_0(t)>0$, the gamma frailty hazard ratio isAt a time when both hazards vanish this formula is a continuous extension; the literal ratio $0/0$ is undefined.

Write $a=e^{\beta }>0$ and $x={\sigma }^2{H}_0(t)\ge 0$. Two useful identities areFor small cumulative exposure $x$, $r=a\{1+(1-a)x+O(x^2)\}$, so it starts at the conditional hazard ratio $e^{\beta }$. If $\beta >0$, it decreases towards one; if $\beta <0$, it increases towards one; and if $\beta =0$, it equals one throughout. It stays between $e^{\beta }$ and one and never crosses one.

For fixed ${\sigma }^2>0$, the late-time limit is one provided $H_0(t)\to \infty$. If instead the integrated hazard tends to a finite value $K$, the limiting hazard ratio is $a(1+{\sigma }^{2}K)/(1+a{\sigma }^{2}K)$, which need not be one. Elapsed time alone does not imply unbounded exposure.

For small ${\sigma }^2$ at a fixed time, $r\to e^{\beta }$. At ${\sigma }^2=0$ the frailty degenerates to $U=1$, giving an ordinary proportional hazards model. For large ${\sigma }^2$ at a fixed time with $H_0(t)>0$, $r\to 1$; at time zero, however, it remains $e^{\beta }$. Thus that variance limit is not uniform near zero. Larger heterogeneity accelerates the movement of the population hazard ratio towards one. The mechanism is survival selection: frailer individuals leave the population sooner, particularly in the higher-hazard group, while each individual's conditional multiplier remains fixed.

Use $F(t)=\mathrm{Pr}(T>t)$, a nonincreasing right-continuous survivor function. The empirical likelihood treats an exact event as a probability mass rather than maximizing an unrestricted probability density function at isolated data points. Under a noninformative observation mechanism, the event-time likelihood contributions areThe second, third and fourth rows account for right censoring, left censoring and interval censoring, respectively. A censoring mechanism with unknown dependence on $T$ cannot simply be dropped; these are the event-time factors under independent censoring or another ignorable mechanism.

There is a useful finite-dimensional reduction before optimization. Let $A_i$ be the set in which observation $i$ tells us its event time lies. Partition the time axis using all observed endpoints, keeping endpoint singletons when required. Within each cell, membership in every $A_i$ is constant. Assign cell masses $w_j\ge 0$, with ${\sum }_j{w}_j=1$, including a tail cell when necessary, and let $A_{ij}$ be the cell-membership indicator. ThenMass locations within a cell are not identifiable and need not be separate parameters. Repeated observation sets can be grouped and their factors raised to their counts. Cells with the same membership column can also be combined. If one column is componentwise dominated by another, transferring its mass to the dominating column weakly increases every factor; hence an unconstrained maximizing distribution can omit the dominated cell. This proves likelihood support reduction under censoring directly. The remaining log-likelihood is a concave function of the probability masses on a simplex, often far smaller than the original unrestricted distribution problem. Added constraints can restrict which transfers remain permissible.

For right censoring alone, take a grid of all distinct observation times $s_j$ and initially allow a discrete hazardLet $d_j$ be the number of events at $s_j$ and $n_j=\#\{i:x_i\ge s_j\}$ the risk set immediately before that time. For tied recorded times, use the convention that events precede censoring, so those censored at $s_j$ remain in this risk set. The mass at $s_j$ is $q_j{\prod }_{k<j}(1-q_k)$, and the probability of surviving past $s_j$ is ${\prod }_{k\le j}(1-q_k)$. Substitution into the exact-event and right-censoring factors yieldsThe exponent $n_j-d_j$ counts all observed individuals whose information requires survival beyond that jump. In particular, censoring observations affect the risk set even though they give no event factor.

Each factor can now be maximized separately. For $0<d_j<n_j$, its derivative is $d_j/q_j-(n_j-d_j)/(1-q_j)$, giving ${\hat{q}}_j=d_j/n_j$; the second derivative is negative. The same result holds at the boundaries: ${\hat{q}}_j=0$ if $d_j=0$, and ${\hat{q}}_j=1$ if $d_j=n_j$. The Kaplan–Meier estimator is thereforeOnly event times contribute nontrivial factors. Residual mass after the last observation may be placed anywhere later without changing the empirical likelihood. A terminal plateau is therefore not evidence of permanent survival; the tail is unobserved.

For a constrained Kaplan–Meier estimator with $F(t)=p$, first split the support cells at $t$. In the probability-mass representation the exact optimization problem isThe last sum excludes a singleton mass at $t$, because $F(t)=\mathrm{Pr}(T>t)$. This is a concave maximization with linear equality constraints, so it directly gives a constrained nonparametric maximum-likelihood estimator. It also handles a time lying between observed events or beyond all follow-up.

An equivalent hazard formulation adds $t$ to the grid, permits zero-event jumps when needed, and imposesFor an interior solution using only positive-event jumps before $t$, add $\lambda \{{\sum }_{s_j\le t}\log (1-q_j)-\log p\}$ to the log-likelihood. Differentiating givesThe multiplier is found by solving the scalar equationwith denominators and fitted hazards in their admissible ranges. On an admissible event-only interior branch, this product is increasing in $\lambda$, making a one-dimensional root search straightforward.

A constraint can require mass at a time with no observed event. The event-only multiplier equation must not be used blindly. To include boundary cases, put $u_j=-\log (1-q_j)\ge 0$. Before $t$, the constraint becomes ${\sum }_j{u}_j=-\log p$, and each contribution isFor $d_j>0$ its multiplier equation again gives $q_j=d_j/(n_j+\lambda )$. For $d_j=0$, its derivative is $-n_j$, so the maximum has $u_j=0$ when $\lambda >-n_j$, and may have $u_j>0$ when $\lambda =-n_j$. Such a jump supplies mass needed by the constraint; an unobserved tail with $n_j=0$ is included by the same rule. For example, if every observation is right-censored after $t$ and there are no events, the ordinary estimate has $F(t)=1$. To impose $p<1$, put mass $1-p$ at $t$ and the remaining $p$ after all censoring times. The maximized likelihood is $p^n$, although there are no event-time factors at all. The mass-constrained optimization above automatically accounts for this possibility and for any remaining tail nonuniqueness.

A court-based randomized controlled trial would assign eligible offenders between the two sentencing policies using randomization, ideally stratified by court and important prespecified risk factors. With allocation concealment, the next assignment cannot be selected to suit an offender's prognosis. Baseline characteristics, both recorded and unrecorded, are balanced in distribution rather than confounded with judges' sentencing preferences. This makes a comparison of outcomes interpretable as a causal effect of assignment, subject to complete and comparable follow-up, whereas separate historical cohorts can differ in severity, social circumstances or calendar period.

The commission should compare both specified outcomes over the same 360-day period from sentencing, using an intention-to-treat analysis. An offender initially assigned community service remains in that group even if the order is later revoked and prison follows. Revocation is part of the policy's consequences; restricting the comparison to completers would select on a post-assignment outcome and lose the protection of randomization.

Within-court individual randomization also provides an objective allocation procedure and valid uncertainty estimates. It does not guarantee identical characteristics in a particular finite sample or eliminate loss to follow-up, unequal outcome ascertainment or limited generalizability. Eligibility, genuine uncertainty about the two policies, and a feasible authorized allocation mechanism would need agreement before the trial. If whole courts rather than individuals were assigned, that would be a cluster-randomized trial and its correlation structure would require a different sample size calculation.

Under the no-death assumption, use the mixture of the two community-service paths. Of 10,000 initially assigned offenders, 4,000 breach and leave their subsequent 60-day prison term at day 120. Their remaining 240 days of follow-up match the supplied reconviction window. The 6,000 completers finish community service at day 180 and have the supplied 180-day post-completion follow-up.

The expected value of the number reconvicted is consequentlyThus the initial-policy 360-day reconviction probability is $0.4(0.8)+0.6(0.4)=0.56$. This is a mixture calculation for the initially sentenced cohort, not the 40% rate of its selected completers alone. It uses the supplied windows and assumes they capture the reconvictions being counted, rather than adding unspecified reconvictions during the earlier order or custody periods.

Assume equal allocation of $n$ independent offenders to each arm and a two-sided statistical test of significance level $\alpha$. Write $p_0=0.60$ for P30 and $p_1=0.56$ for initial community-service assignment, with $\Delta =p_0-p_1=0.04$ and $\bar{p}=(p_0+p_1)/2=0.58$. The two event counts have independent binomial distributions. Under the stated alternative, the difference of sample proportions is approximately a normal distribution with mean $\Delta$ and varianceThe pooled two-proportion test has approximate upper rejection threshold $z_{1-\alpha /2}\sqrt{2\bar{p}(1-\bar{p})/n}$ at the planning rates, where $z_q$ is the standard normal quantile. To obtain statistical power $0.80$ against a positive difference, place the alternative mean $z_{0.80}$ alternative-standard-errors above that threshold:Solving yields the usual sample size for comparing two proportions:per arm, rounded upwards, with total sample size $2n$. This targets the dominant rejection tail under the planning alternative; an exact power calculation can include the very small opposite-tail probability or use the binomial distributions directly.

For illustration, at $\alpha =0.05$, $z_{0.975}\simeq 1.960$ and $z_{0.80}\simeq 0.842$, giving $n\simeq 2388.8$. ThusFor a prespecified one-sided test, replace $z_{1-\alpha /2}$ by $z_{1-\alpha }$. The calculation assumes complete endpoint ascertainment, independent individual allocation and that the two initial-policy rates are appropriate; losses or clustered allocation would require corresponding design adjustments.

Interpret the quoted daily death rates as piecewise constant hazard functions, and use the custody schedules described in the problem. This means community service itself has the ordinary community hazard, a completer has no prison exposure, and a breacher is imprisoned from day 60 to day 120. Subsequent imprisonment after another conviction is not specified, so it is not added to the exposure schedule. The abbreviated “CS180” denotes the same initial community-service policy.

For a completer, all 360 days have the baseline rate, giving integrated hazardFor a breacher, the four periods are 60 community days, 60 prison days, 30 high-risk days immediately after release, and 210 further community days. ThusThe piecewise-exponential survival model gives death probabilities $1-e^{-H_C}$ and $1-e^{-H_B}$. Average these over the stated 60%/40% mixture of paths. Equivalently, everyone first survives the common initial 60-day community interval, and survivors then take the two paths in those proportions. The expected value among 10,000 initial assignments isThis accounts for removal of people after death. Since the daily hazards and endpoint death probabilities are small, the usual person-time approximation replaces $1-e^{-H}$ by $H$ and givesBoth calculations use the same specified schedules. The latter treats all planned exposure as if contributed and is slightly larger. Additional custody histories would change the exposure-weighted integrated hazard; the provided reconviction proportions alone do not determine those histories.

Mortality requires its own sample size calculation, rather than using the much commoner reconviction endpoint. Under the same specified exposure assumptions, the P30 schedule has 30 prison days, 30 high-risk post-release days, and 300 ordinary days. Its integrated hazard isHence its 360-day death probability is $p_P=1-e^{-0.0145}\simeq 0.0143954$. For initial community-service assignment the preceding mixture gives $p_C\simeq 0.0127180$. The absolute risk difference is therefore only about $0.0016774$, or 16.8 fewer deaths per 10,000 initial assignments. The person-time approximation gives 145 versus 128 deaths per 10,000, a difference of 17.

As an illustrative design target, take two-sided significance level $0.05$, statistical power $0.80$ and equal individual allocation. Apply the sample size for comparing two proportions formula with these mortality probabilities:It gives about 74,600 participants per arm, or 149,200 altogether; using the simpler exposure approximations gives about 73,130 per arm. At 30,000 sentences per year, even if every offender were eligible and enrolled, this requires roughly five years of recruitment and another year to complete the final cohort's follow-up.

For scale, one year's equally allocated cohort has 15,000 per arm, with roughly 216 expected deaths after P30 and 191 after community-service assignment. The standard error of the risk difference is about $0.001335$, so the expected standardized separation is only about $1.26$, below the usual two-sided 5% critical value. A one-year trial would consequently have poor statistical power for this mortality difference.

Thus mortality discrimination is much more demanding than reconviction discrimination. National recruitment and reliable linkage to death records could make a multiyear study feasible, but restricted eligibility, incomplete follow-up, or lower recruitment would lengthen it. A trial sized only for the 4-percentage-point reconviction difference would not provide a precise mortality comparison. The numerical conclusion depends on the supplied daily hazards, custody schedules, and the chosen power and significance targets; the prompt does not specify a mortality design target of its own.

Let $T_i=X_i+Y_i$. Since the two counts are independent Poisson distributions with means ${\lambda }_i$ and $\beta {\lambda }_i$, their total is a Poisson distribution with mean $(1+\beta ){\lambda }_i$. For a fixed observed total $t_i=x_i+y_i$, divide the joint probability mass function by that of the total:Define $p=\beta /(1+\beta )$, so $1-p=1/(1+\beta )$. ThenFor independent patients, the paired Poisson conditional likelihood isEvery nuisance rate ${\lambda }_i$ cancels, while each patient's own baseline count level is retained through the conditioning total. A patient with total zero would contribute the constant one.

For these data, ${\sum }_i{x}_i=161$, ${\sum }_i{y}_i=23$, and ${\sum }_i{t}_i=184$. Maximizing this conditional likelihood givesThese are the fit under the model without a cure component, not the estimates of the mixture model in the later parts.

Work with the stipulated approximate conditional likelihood. Let $Z_0=\{i:y_i=0\}$ and $Z_{+}=\{i:y_i>0\}$, and putTerms independent of the parameters, including the positive-count binomial coefficients, can be omitted. The log-likelihood isFor an interior maximum-likelihood estimator, set its two score functions to zero. Differentiating with respect to $\theta$ givesDifferentiating with respect to $p$ givesHere there are seven zero-count patients and five positive-count patients, with ${\sum }_{Z_{+}}{y}_i=23$ and ${\sum }_{Z_{+}}(t_i-y_i)=36$. Thus the positive-count contribution to the second equation is $23/p-36/(1-p)$.

An interpretable way to solve the equations uses a posterior cure weight for each observed patient:For a zero count, $1-q_i=(1-\theta )(1-p{)}^{t_i}/A_i$ is the posterior probability of belonging to the non-cured component under this approximate conditional model. The two score functions can then be rewritten asTherefore the conditional cure-weight equations for paired Poisson counts areThey can be solved by an expectation-maximization algorithm: calculate the weights at the current parameters, update $\theta$ to their mean and $p$ to the weighted event fraction, and iterate. A solution of the equations must be checked as a maximum; the equations alone do not guarantee a global maximum. In particular the boundary $\theta =0$ gives the no-cure estimate $p=1/8$ and should be compared. Boundaries $p=0$, $p=1$ or $\theta =1$ have zero likelihood for these positive counts with positive pre-treatment counts. Numerical maximization gives an interior fit about $(0.57549,0.38614)$, with larger log-likelihood than the no-cure boundary, consistent with the PDF's rounded estimates.

The word “approximate” is important. If cure is specified before conditioning in the original Poisson distributions, the conditional cure probability given $T_i=t_i$ would actually bewhich still depends on ${\lambda }_i$. This follows because a cured total has mean ${\lambda }_i$ and a non-cured total has mean $(1+\beta ){\lambda }_i$. The likelihood supplied in the question instead treats $\theta$ as the mixture probability after conditioning; the score and weight equations above maximize exactly that stipulated approximation, without claiming nuisance-rate elimination for the full cure model.

Use the original PDF's value $\hat{p}=0.386$; the TeX aid misreads its last digit. In the fitted conditional mixture, $\hat{\theta }=0.575$ means that about 57.5% are assigned to the structural-zero, or modelled cured, component. Its asymptotic 95% confidence interval of $(0.30,0.81)$ shows substantial uncertainty about that fraction. An observed post-treatment zero alone is not proof of cure, since the non-cured Poisson distribution can also produce zero.

The parameter $p$ is a conditional allocation probability, not the ratio of post-treatment to pre-treatment means. That ratio among non-cured patients is $\beta =p/(1-p)$. HenceThe map $p\mapsto p/(1-p)$ is strictly increasing, so transform the endpoints of the supplied 95% confidence interval:No change in the non-cured mean is $\beta =1$, equivalently $p=1/2$, which lies in this interval. Thus the point estimate suggests a reduction, but this interval does not firmly establish one among the non-cured patients at the corresponding two-sided 5% level. It is not correct to describe $p=0.386$ itself as a 61.4% reduction in their mean count.

Under a full pre-conditioning mixture with cure independent of the baseline rate, the population mean ratio would be $(1-\theta )\beta$, whose plug-in value is about $0.425(0.629)=0.267$. This combines disappearance of counts in one component with a change of rate in the other; it is distinct from the conditional estimate's direct interpretation. The two marginal confidence intervals cannot simply be multiplied to obtain a calibrated joint 95% interval for this derived quantity.

Only 12 patients were observed, so the asymptotic intervals, particularly for a mixture parameter near a boundary, should be interpreted cautiously. The data describe one-minute counts and do not establish permanent cure. Nor does an uncontrolled before-and-after comparison alone prove a causal effect of the drug: natural count variability, secular changes and regression to the mean remain possible. The structural-zero interpretation and the reduction among non-cured patients are conclusions within the specified statistical model.