Past exam of the mathematics course of the University of Cambridge / 2018 / iii / Paper 207 (Statistics in Medicine)

Write $Z(u)$ for the continuous-time Markov chain and ${\mathbb{E}}_r$ for expected value conditional on $Z(0)=r$. The occupation time of a continuous-time Markov chain in state $s$ is the integral of its indicator random variable. By Tonelli theorem,The interchange is valid because the integrand is nonnegative; indeed the random integral is bounded by $t$.

Let $N(u)$ count all jumps of the continuous-time Markov chain. Its Markov-chain entry count into state $s$ isThusThis is a Lebesgue-Stieltjes integral against a counting process, so each jump contributes one. The initial occupation of state $s$ is not counted as an entry. An ordinary time integral of an indicator random variable asserting a jump at exactly time $u$ would be zero, since individual jump times have zero Lebesgue measure. The requested integral must therefore use this counting-measure interpretation. Assume non-explosion and finite expected jump counts on bounded intervals, as holds for a finite-state continuous-time Markov chain.

For $i\ne s$, the Markov property and the definition of transition intensity giveSumming over the possible pre-jump states gives the density of the expected Markov-chain entry count:For a finite-state continuous-time Markov chain, the chance of two or more jumps in this short interval is $o(du)$. There is no $q_{ss}$ contribution: diagonal entries of the infinitesimal generator encode exit rates, not jumps into the same state. With repeated entries possible, this density is an expected event rate rather than the probability density function of a single normalized random time.

Integrate the expected Markov-chain entry count density and use the occupation time of a continuous-time Markov chain formula:Each term is a transition intensity multiplied by the expected exposure time to that transition. This counts every entry, including returns, rather than only the first visit. For countably many states the same nonnegative interchange follows from Tonelli theorem, provided the resulting expected count is finite. The sum excludes $i=s$.

Use an irreversible continuous-time multi-state model $1\to 2\to 3$, with state 3 absorbing. Let $a=q_{12}(50)$ and $b=q_{23}(50)$, and define log-linear age effects byFor a closed-form panel-observed multi-state likelihood, freeze the age covariate at the start of each observation interval. The transition intensities are then $(a,b)$ over the first interval and $(\lambda ,\mu )=(a{e}^{5{\beta }_{12}},b{e}^{5{\beta }_{23}})$ over the second. Conditional on the initial state, the likelihood contribution isIndeed, $p_{11}(5)=e^{-5a}$ and $p_{13}(1)=1-p_{11}(1)-p_{12}(1)$. In this irreversible continuous-time Markov chain, the intermediate state has probabilityThe second expression is the continuous limit of the first and handles equal rates. The observed two-state jump across an interval is compatible with adjacent instantaneous transitions: an unobserved visit to state 2 occurs between the examinations.

The four-parameter panel-observed multi-state likelihood requires the following additional assumptions.

If age is instead updated continuously, the model is time-inhomogeneous and the displayed homogeneous transition formula is not exact. For example, setThen the exact irreversible likelihood contribution isThus the simple closed form uses an explicit frozen-age approximation; adjacency alone and continuously evolving age do not determine that form.

Under the frozen-age exponential distribution convention, the mean waiting time is the reciprocal of the outgoing transition intensity. The beliefs therefore specifyBoth rates double over ten years. ConsequentlyThese calibrate the age-specific constant-rate means; they are not exact mean ages of future events when age-dependent rates continue increasing throughout the waiting time.

For a literal continuously ageing interpretation, a rate $q{e}^{\beta u}$ gives a Gompertz distribution for the waiting time. For $\beta >0$, its mean isHere $E_1$ is the exponential integral. Solving $m(a,{\beta }_{12})=10$ and $m(a{e}^{10{\beta }_{12}},{\beta }_{12})=5$, and the analogous pair with means 2 and 1, instead giveswith rates and slopes in inverse years. These are an alternative continuous-age calibration, not the frozen-age answer boxed above.

Put $v_k=p_k(1-p_k)$ and assume $0<p_k<1$. For fixed $n$ and a fixed nonzero difference $\Delta =p_0-p_1$, the approximate statistical power of the Wald test increases as its noncentrality $|\Delta |/\sqrt{V}$ increases, whereDifferentiating with respect to the continuous allocation givesHence the unique minimum and its associated minimum variance areThis is Neyman allocation. Integer allocations use a feasible adjacent integer to the continuous optimum, comparing the resulting variances; $n_0,n_1\ge 1$ must still hold. It is a large-sample normal approximation argument, not an exact finite-sample power theorem for arbitrary response-adaptive randomization. Under the null $p_0=p_1$, the approximate power is the significance level for any allocation, although the same allocation still minimizes the stated variance.

The difference of the Bernoulli distribution variances factors asThe first factor is positive, whereas $p_0>1-p_1$ makes the second negative. Thus $v_0<v_1$ and $R^{\ast }=\sqrt{v_0/v_1}<1$. Therefore Neyman allocation assigns more patients to treatment 1, the inferior treatment. Minimizing the variance of a treatment comparison and maximizing the number of desirable responses are different objectives.

For the plug-in Neyman allocation, the estimated Bernoulli distribution variances are $6/25$ and $3/16$, soThis uses the target ratio directly as the allocation odds; a procedure that actively corrects previous allocation imbalances would need an additional rule.

For the randomized play-the-winner rule RPW$(1,1)$, return the drawn ball and add one ball of the same treatment after a success or of the opposite treatment after a failure. Treatment 0 has three successes and two failures, and treatment 1 has one success and three failures. The urn therefore has $1+3+3=7$ treatment-0 balls and $1+1+2=4$ treatment-1 balls. HenceFor dynamic programming with only the final patient left, there is no future value from learning. Independent uniform prior distributions and Beta-binomial conjugacy give posteriors $p_0\mid \mathcal{D}\sim \mathrm{Beta}(4,3)$ and $p_1\mid \mathcal{D}\sim \mathrm{Beta}(2,4)$. Their posterior predictive probabilities of success are $4/7$ and $2/6=1/3$. The optimal terminal action therefore givesThis is optimal for expected successes with no imposed lower bound on the randomization probabilities.

A second opportunity to reject raises the trial's overall Type I error above the nominal 10% in general. If $A$ and $B$ denote rejection at the first and second analyses, the overall error isThe cumulative-data test statistics are correlated, so neither 20% nor $1-0.9^2=19\%$ is the general answer. The latter requires independent tests, which is not the situation here.

The rejection thresholds and the stopping/continuation rule must be calibrated jointly. A prespecified group sequential design or alpha-spending function can allocate a total error budget of 10% across the analyses, using their joint null distribution. A conservative alternative planned in advance is two tests each at 5%, controlled by the Bonferroni inequality. An unplanned change to a design with a remaining error budget can be justified by the conditional error principle. Here, however, an original single-look test at 10% has no remaining conditional rejection probability after its final nonrejection. Keeping its entire original rejection region and adding a new positive-probability rejection opportunity cannot in general retain a 10% overall error rate. The two-look rule therefore needs prospective joint calibration; an already observed result cannot retrospectively create an unused error budget.

Let $Y(1)$ and $Y(0)$ be the potential outcomes for the same individual under well-defined interventions setting the binary exposure to 1 and 0. The average treatment effect, also called the average causal effect on the additive scale, isThe average is over the target population. This compares hypothetical outcomes under two interventions, rather than the observed difference $\mathbb{E}[Y\mid X=1]-\mathbb{E}[Y\mid X=0]$, which can be distorted by confounding. For binary $Y$, the displayed contrast is a causal risk difference. Consistency in causal inference and a well-defined intervention, including an appropriate no-interference in causal inference assumption, give these potential outcomes their meaning.

An instrumental variable $G$ supplies variation in an exposure $B$ that can be used for causal inference about its effect on $Y$. Its three core requirements are:

Relevance, independence, and exclusion are the core instrumental-variable assumptions. They do not by themselves identify an arbitrary population-average effect. For a binary instrument and treatment, instrumental-variable monotonicity and the usual consistency in causal inference assumptions can identify a local average treatment effect; stronger effect-homogeneity or structural-model assumptions are needed for other effect interpretations.

Under the displayed causal directed acyclic graphs, with nonzero effects along the relevant arrows and no exact cancellation of the effects on $B$, the valid instrumental variables are scenarios (i), (ii), and (iv). The path arguments are:

An upstream intermediate variable need not invalidate an instrumental variable if its routes to the outcome all pass through the target exposure; an outcome intermediate reached without that exposure does invalidate it.

The supplementary excerpt uses an adjusted observational study and a two-sample Mendelian randomization analysis.

The observational study examines co-occurrence of hypercalcaemia and migraine diagnoses in electronic health records, adjusting for age, sex, and ancestry. Its adjusted odds ratio is 1.58, with $P=4.75\times 10^{-13}$.

The two-sample Mendelian randomization analysis combines the associations of eight genetic variants with serum calcium and with migraine, using a weighted genetic risk score. It estimates an odds ratio of 1.80 per genetically predicted 1 mg/dL higher serum calcium, with 95% confidence interval $(1.31,2.46)$. The weighted-median Mendelian randomization estimator and MR-Egger regression provide sensitivity analyses.

The measured and instrumented exposure is serum calcium, rather than dietary calcium intake itself. Calcium intake and serum calcium are not interchangeable interventions, and an odds ratio of 1.80 is an odds contrast rather than necessarily an 80% increase in probability. The shortened excerpt supplies two approaches; a further analysis in the full research paper is outside this excerpt.

For the electronic-record observational study, three strengths are:

Three weaknesses are:

For two-sample Mendelian randomization, three strengths are:

Three weaknesses are:

Thus the agreement supports the hypothesis, but does not establish the effect of a dietary intervention.

Two useful ways to strengthen the causal inference are:

These are complementary: the first changes the proposed exposure experimentally, and the second probes the validity and reproducibility of its genetic evidence.

The survivor function of a nonnegative event time isHere the paper uses $F$ for the survival function, rather than for the cumulative distribution function; the latter is $1-F(t)$ for a continuous distribution. The survivor function is nonincreasing, starts at 1 when there is no atom at zero, and tends to zero for an almost surely finite event time.

With no censoring or new entry within the interval, all $r$ individuals can be classified at its end: $d$ have experienced the event and $r-d$ remain event-free. Thus the empirical conditional survival probability isThis assumes $r>0$ and that the observed individuals represent the group whose conditional survival probability is being estimated.

An individual is right-censored at $t$ when follow-up establishes that the event has not occurred through $t$, but does not reveal the subsequent event time. The observed information is $T>t$ for the continuous-time convention.

Knowing only the total number censored inside an interval does not specify the risk sets at the event times. For example, start with two individuals, one event and one censoring. If censoring precedes the event, the Kaplan–Meier estimator contributes $1-1/1=0$; if the event precedes censoring, it contributes $1-1/2=1/2$. The totals $r=2,d=1,c=1$ are identical.

Therefore the timing of right censoring relative to the events is needed. Independent censoring is also required for these changing risk sets to yield a valid estimate of the underlying survivor function.

The estimate remains $(r-d)/r$. Individuals right-censored at $t_R$ are observed throughout the preceding interval, so they contribute to each risk set until the end. They must not be removed from the denominator at $t_L$.

If events and right censoring share the recorded time $t_R$, use the usual Kaplan–Meier estimator convention: count events first, using a risk set that still includes those censored at that time, and then remove the censored observations. With no earlier censoring, the product of event factors over the interval telescopes to $(r-d)/r$.

A potential event time is an allowed support point for the discrete event-time distribution used to derive the Kaplan–Meier estimator. A positive probability mass may be assigned there before fitting, even if no event was observed at that time.

This is a device for nonparametric maximum likelihood estimation from exact event times and right censoring; it does not assert that the underlying biological event time is genuinely discrete. Its estimated mass can be zero at an unobserved potential time.

The set must contain every observed event time. Otherwise a discrete survival likelihood would assign probability zero to at least one exact event observation.

The set can include additional times without observed events. Include time zero when events at zero are present, and allow remaining tail probability beyond the fitted range. Estimation of a conditional event mass requires a nonempty risk set; once follow-up ends, later probability masses are not identified by the data. The usual Kaplan–Meier estimator retains the final estimated survival probability as its tail convention.

Let $d_j$ be the number of events at potential time $a_j$ and let $r_j$ be the size of the risk set immediately before that time, including observations censored at $a_j$. Define the discrete conditional event probabilityUnder independent censoring, and ignoring the separate censoring distribution, the survival likelihood factorizes asThis follows by multiplying a survival factor for each individual at risk who does not fail at $a_j$, and an event factor for each failure. The log factor $d_j\log q_j+(r_j-d_j)\log (1-q_j)$ is maximized at ${\hat{q}}_j=d_j/r_j$, with the corresponding boundary values when $d_j=0$ or $d_j=r_j$.

The survival function is the product of successive conditional survival probabilities. ThereforeThis is the Kaplan–Meier estimator. A product over no event times is 1. If everyone remaining experiences an event, the product becomes zero and stays zero; subsequent empty risk sets supply no additional factor.

At an added potential event time with no event and a nonempty risk set, $d_j=0$, so maximum likelihood estimation gives ${\hat{q}}_j=0$. Its factor in the Kaplan–Meier estimator is therefore $1-d_j/r_j=1$.

Adding or deleting such times also leaves the risk sets at actual event times unchanged: those depend on the observed event and censoring histories. Consequently the fitted survivor function depends only on observed event times and their risk sets. Beyond the last observed follow-up, the tail is not identified; the standard flat continuation is a convention, rather than an additional consequence of likelihood maximization.

Choose the potential right-censoring times to include every observed right censoring time up to $t$, and also the endpoints 0 and $t$. Then no individual is censored strictly inside $(c_{k-1},c_k)$. Use an explicit endpoint convention: let $r_k$ count individuals still event-free after any events at $c_k$, but before censoring at $c_k$. If $d_k$ events occur in $(c_{k-1},c_k]$, the number at the beginning of that interval, after censoring at its left endpoint, is $r_k+d_k$.

The no-censoring interval estimate from the preceding parts is therefore $r_k/(r_k+d_k)$. The possible event mass at zero contributes $r_0/(r_0+d_0)$. Multiplying yieldsFor a continuous event-time distribution, normally $d_0=0$. Within interval $k$, put $b_k=r_k+d_k$ and let its successive event counts be $e_1,\dots ,e_m$. Its Kaplan–Meier estimator factors telescope:The same calculation at zero provesThe paper's phrase “at risk at $c_k$” does not specify its side of an event at that endpoint. If it denotes the conventional pre-event risk set size $r_k^{-}$ instead, and $e_k$ events occur exactly at $c_k$, substitute $r_k=r_k^{-}-e_k$. The interval factor is then $(r_k^{-}-e_k)/(r_k^{-}+d_k-e_k)$; at zero it is $(r_0^{-}-d_0)/r_0^{-}$. This keeps the events-before-censoring convention consistent. After an empty risk set, stop the product; if the fitted survival has reached zero, retain zero rather than creating a $0/0$ factor.

Define the event counting process and at-risk process byUnder a common hazard function and independent censoring, the conditional event rate is $Y(u)h(u)$. Equivalently, the counting-process compensator satisfies $\mathbb{E}[dN(u)\mid {\mathcal{F}}_{u-}]=Y(u)dH(u)$. Estimate a cumulative hazard increment by the observed event increment divided by the current risk set:With no ties, $dN$ has mass $v_i$ at $x_i$, and $Y(x_i)={\sum }_{j:x_j\ge x_i}1$. Thus the Nelson–Aalen estimator isOnly event times contribute, so no empty risk set is divided into an event count. Right-censored observations contribute time to the risk sets, but no event increment.

Let $R_i=\{j:x_j\ge x_i\}$ be the risk set at an observed event time $x_i$. In the Cox proportional-hazards model, conditional on one member of this set failing then and on the preceding history, the probability that it is individual $i$ isThe baseline hazard cancels. Multiplying these conditional contributions gives the Cox partial likelihoodThe no-ties assumption avoids the need for a tied-event approximation. This construction assumes independent censoring conditional on the modelled covariates. The printed phrase calling $z$ a scalar parameter is read as a scalar observed covariate $z_i$; $\beta$ is the unknown regression parameter.

Taking logarithms of the Cox partial likelihood givesDifferentiate each term, using the chain rule for the logarithm:This score function is the sum over failures of the observed covariate minus its current risk-weighted mean. Define that mean as ${\overline{z}}_{\beta }(x_i)$; observations with $v_i=0$ have no direct event term, although they occur in the relevant risk sets.

Within the risk set $R_i$, the Cox proportional-hazards model makes the conditional failure probabilities proportional to $e^{\beta z_j}$. ThusThe score function sums observed minus expected covariates at event times. At a finite interior maximum of the Cox partial likelihood,The balance holds in aggregate, not separately at every event. The log partial likelihood is concave becauseWhen there is informative covariate variation in the risk sets, this gives a unique finite solution if it exists. A completely separated event pattern can instead have its maximum only as $\beta \to \pm \infty$; the question's finite $\hat{\beta }$ implicitly excludes that case.

In the Cox proportional-hazards model, replace the number at risk in the Nelson–Aalen estimator by the sum of relative hazards:The total event counting process has conditional rate $W_{\beta }(u)h_0(u)$, so the estimated baseline hazard increment is $dN(u)/W_{\hat{\beta }}(u)$. Integrating gives the Breslow estimatorWhen $\hat{\beta }=0$, all relative hazards are 1 and this reduces to the Nelson–Aalen estimator. The denominator is positive at each observed event time.

In the martingale residual$v_i$ is the observed event count up to the end of the individual's follow-up: 1 for a recorded event and 0 for right censoring. The second term is the fitted counting-process compensator, obtained by integrating that individual's fitted hazard function while they remain at risk:Thus the residual is observed events minus the fitted cumulative event intensity over follow-up. The fitted intensity is not an event probability and can exceed 1. Positive values indicate more events than predicted over that observed exposure; negative values indicate fewer. A right-censored individual has a nonpositive residual. With true model parameters and independent censoring, the underlying observed-minus-compensator process is a martingale; fitted terminal residuals are not independent identically distributed errors.

Fit the Cox proportional-hazards model with the other explanatory variables, initially omitting the variable $z$ whose form is being assessed, and plot the resulting martingale residuals against $z$. Add a smooth trend, for example using a cubic spline or a local smoother. A systematic pattern indicates that $z$ affects the hazard function; the shape helps suggest a linear term, a transformation, or a nonlinear spline term in the log hazard.

Refit with the proposed form and inspect martingale residuals against $z$ again. Under an adequate functional form, their conditional mean should have no systematic trend, allowing for model fitting and the other covariates. Compare plausible alternative forms using the partial likelihood and suitable model checks.

Use the smoothed pattern to assess the covariate's functional form, rather than treating these residuals as normal regression errors. Martingale residuals are asymmetric, bounded above by 1, and can be arbitrarily negative; the raw residual plot can therefore be misleading without smoothing. This diagnostic concerns the covariate form, and does not by itself test the proportional hazards assumption.

Insert the Breslow estimator into the fitted cumulative intensities and interchange the finite sums:Every observed event contributes exactly one after summing its fitted contribution over its risk set. ThereforeThis identity is algebraic and holds for any finite coefficient used consistently in the relative hazards and the Breslow estimator; it does not require the score equation to be solved exactly. It can fail if the baseline hazard is estimated or smoothed by a different procedure.

Condition on the event time $T_A=u$. Independence makes the chance that $T_B$ is still later equal to its survivor function $F_B(u)$. HenceThe continuous distributions give zero probability of a tie. The complementary probability is ${\pi }_B={\int }_0^{\infty }{F}_A(u)f_B(u)du=1-{\pi }_A$.

For the survivor functions, $F_A^{\prime }(u)=-f_A(u)$ and $F_B^{\prime }(u)=-f_B(u)$. The product rule givesIntegrating from 0 to $t$, using $F_A(0)=F_B(0)=1$, yieldsThe three terms partition the possible histories: A occurs first by time $t$, B occurs first by time $t$, or neither has occurred by $t$. The first two are the cumulative incidence functions for the two competing risks, while the third is the survival probability of their minimum. The initial value assumes nonnegative event times with no atom at zero.

Let ${\pi }_A=\mathbb{P}(T_A<T_B)$, and assume $0<{\pi }_A<1$ so that both conditional distributions are defined. By independence, the joint density of $T_A=u$ and $T_A<T_B$ is $f_A(u)F_B(u)$, whereas the joint density of $T_A=u$ and $T_B<T_A$ is $f_A(u)[1-F_B(u)]$. Normalizing givesBoth are densities of the latent event time $T_A$. The second is not the density of the observed first event $T_B$ in cases where B wins. Their integrals are 1 by the definitions of ${\pi }_A$ and its complement.