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

The first subscript in ${\beta }_{t,\tau }$ is the calendar time $t$, and the second is the infection age $\tau$, the time elapsed since the source individual became infected. Thus ${\beta }_{t,\tau }$ is the rate at which an individual of infection age $\tau$ generates infections at time $t$. The corresponding discrete infectious disease renewal equation isup to a separately modelled term for imported infection. The upper limit may instead be a fixed maximal infectious age, with unavailable terms set to zero.

The instantaneous reproduction number at time $t$ isIt is the expected number of secondary infections that one infected individual would produce if the transmission conditions at time $t$ applied throughout that individual's infectious life. The case reproduction number, also called the effective reproduction number in this question, isthe expected number actually produced by a person infected at $t$ as calendar-time conditions subsequently change. The instantaneous quantity is easier to estimate in real time because it depends on current and past incidence; the case quantity depends on future conditions.

Assume the infectivity profile is separable:Then $g_{\tau }$ is the discretized generation-interval distribution, $R_t$ is the instantaneous reproduction number, and the infectious disease renewal equation becomesHence $R_t={\Delta }_t/{\Lambda }_t$ whenever the total infectiousness ${\Lambda }_t$ is positive.

If incidence grows exponentially, ${\Delta }_t=C{e}^{\rho t}$, substitution in the infectious disease renewal equation gives the discrete Euler-Lotka equationFor $X\sim \mathrm{Geometric}(p)$ on $\{0,1,\dots \}$, $g_i=\mathbb{P}(X=i-1)=p(1-p{)}^{i-1}$. The finite geometric series therefore yieldsand consequently

Treat the diagnosis counts as the observed incidence series and define their total infectiousness byA conditionally independent Poisson observation model for the renewal process isInitial infections before the observation window can be included in the definition of $s_t$.

Use the shape-rate convention $R_k\sim \mathrm{Gamma}({\alpha }_0,{\beta }_0)$. The factors involving $R_k$ in the posterior density areBy Poisson-gamma conjugacy,so its posterior mean is $(I_k+{\alpha }_0)/(s_k+{\beta }_0)$.

Use a three-state continuous-time multi-state model with transient infected state $I$ and absorbing recovered and dead states $R$ and $D$. If recovery and death have constant transition intensities $\gamma$ and $\mu$, its Q-matrix iswith state order $(I,R,D)$. This is also a competing risks model: recovery and death are the two mutually exclusive first events.

The two event clocks have independent exponential distributions with rates $\mu$ and $\gamma$. The competing exponential clocks identity gives

The minimum of the recovery and death clocks has an exponential distribution with rate $\mu +\gamma$. Its expected value is therefore

The transition probability matrix has entriesunder the time-homogeneous Markov property. For the finite-state model in part d, the transition semigroup of a continuous-time Markov chain is

Assume a negative test identifies the recovered state. Person 1 remains infected through day 7 and recovers during $(7,14]$, so their interval-censored contribution isPerson 2 remains infected through day 7 and makes an exact $I\to D$ transition at day 10. A transition at an exact time contributes a transition probability density, givingAssuming independent individuals, their joint likelihood contribution is

The model assumes a time-homogeneous Markov property: transition intensities depend only on the current state, not on calendar time, infection age, or earlier history. A Semi-Markov multi-state model could let recovery and death hazards depend on time since infection or entry into the current state. Weekly tests reveal recovery only by interval censoring, however, so the relevant entry and transition times are latent. Fitting the less restrictive model then requires integrating over unobserved paths, and the available data may contain too little information to identify a flexible duration-dependent hazard.

Add a hospital state $H$ and estimate the enlarged Q-matrix from the cohort. If infection starts in state $I$, the expected hospital occupancy generated by one infection is the Markov reward model integralEquivalently, if $Q_T$ is the subgenerator on all transient states, $m_H$ is the $(I,H)$ entry of the fundamental matrix of an absorbing continuous-time Markov chain $(-Q_T{)}^{-1}$. The number infected in the population has expectation $Np$, so linearity of expectation gives total expected hospital use $Np,m_H$ days. If only occupancy within a finite period counts, replace the upper integration limit by the remaining follow-up time for each infection and average over infection times.

Let ${\hat{p}}_k=X_k/n_k$, where $X_k\sim \mathrm{Binomial}(n_k,p_k)$ independently, and putThe Wald statistic is $Z=({\hat{p}}_1-{\hat{p}}_0)/\sqrt{\hat{V}}$. The central limit theorem and Slutsky theorem give $Z\dot{\sim }N(0,1)$ under $H_0$. When $p_1-p_0={\delta }^{\ast }>0$,

For equal arm size $n$, write $v^{\ast }=p_1(1-p_1)+p_0(1-p_0)$ at the clinically relevant alternative. A one-sided level-$\alpha$ Wald test rejects when $Z>z_{1-\alpha }$, and its approximate statistical power isEquating this to $1-\beta$ gives the per-arm sample sizerounded up. If the design uses a null-based critical standard error $v_0$ but an alternative standard error $v^{\ast }$, the corresponding more general formula is

A group sequential design can reduce the expected sample size by stopping at an interim analysis once efficacy or futility is sufficiently clear. It does not generally reduce the prespecified maximum sample size: repeated opportunities to reject inflate the Type I error, so valid sequential stopping boundaries usually require a modest increase in maximum information relative to a fixed-sample design with the same power. The benefit is a smaller expected sample size under alternatives that often cross an early boundary, and sometimes under the null through early futility stopping.

Letwhere $m_j$ is the cumulative sample size per arm. The canonical joint distribution for group sequential test statistics isThus $\mathrm{Corr}(Z_1,Z_2)=\sqrt{m_1/m_2}$. This correlation arises because the second statistic reuses all first-stage observations.

Under $H_0$, $(Z_1,Z_2)$ is a bivariate standard normal distribution with correlation $\rho =\sqrt{m_1/m_2}$. Rejection occurs either at stage 1 through $Z_1\ge u_1$, or at stage 2 through $Z_2\ge u_2$ after continuation $l_1<Z_1<u_1$. Hence the Type I error isThe final lack-of-benefit boundary $l_2$ affects acceptance, but not the probability of crossing an efficacy boundary.

Response-adaptive randomization can assign a larger proportion of later participants to the treatment currently estimated to be better, improving outcomes for participants within the trial. Its allocation probabilities depend on earlier outcomes, which complicates statistical inference; delayed responses and calendar-time trends can also make adaptation ineffective or biased, and an allocation aimed at patient benefit need not maximize power.

Under Neyman allocation, sample sizes are proportional to the arm standard deviations. HereFor total size $n_{\max }$, the minimized asymptotic variance isEqual allocation givesThe Neyman allocation therefore reduces the large-sample variance by $0.04/n_{\max }$, about $5.9\%$ of the equal-allocation variance.

Let $K\in {\mathbb{R}}^{n\times n}$ have entries $K_{ij}=\kappa (x_i,x_j)$, let $k_{\ast }\in {\mathbb{R}}^n$ have entries $(k_{\ast }{)}_i=\kappa (x_i,x_{\ast })$, and let $k_{\ast \ast }=\kappa (x_{\ast },x_{\ast })$. The Gaussian process prior and independent Gaussian noise implyApplying the conditional multivariate normal distribution gives the Gaussian process regression posteriorwhere

Set $x_{\ast }=0.5$ in the Gaussian process regression posterior and compute $m_{\ast }$ and $v_{\ast }$ as in part a. Standardization of a normal random variable then giveswhere $\Phi$ is the standard normal cumulative distribution function.

Let $N_m=|\{i:c_i=m\}|$. The label likelihood for the mixture weights is proportional to ${\prod }_{m=1}^3{\pi }_m^{N_m}$. Multiplication by the $\mathrm{Dirichlet}({\alpha }_1,{\alpha }_2,{\alpha }_3)$ density and Dirichlet-multinomial conjugacy givesConditional on the labels, the observations and component means provide no further information about $\pi$.

For component $m$, stack its $N_m$ assigned observations as $y_1^{(m)},\dots ,y_{N_m}^{(m)}\in {\mathbb{R}}^n$. The prior is $f_m\sim N_n(0,K)$ and each assigned vector is conditionally $N_n(f_m,{\sigma }^2{I}_n)$. Normal-normal conjugacy giveswhereIf $N_m=0$, this reduces to the prior $N_n(0,K)$.

Bayes theorem turns the categorical distribution prior probabilities and the component multivariate normal densities intoThese probabilities define the label update in the Gibbs sampler.

A Dirichlet process mixture model avoids fixing the number of occupied functions. Let $G_0=N_n(0,K)$ be the finite-dimensional Gaussian process law on the common input grid and specifyA draw from a Dirichlet process is almost surely discrete, so several $f_i$ coincide and thereby form clusters. The number of occupied clusters is random and can grow with the data.

The risk set at event time $a_j$ contains individuals still under observation and event-free immediately before $a_j$. For group $k$ its size isThe use of $\ge$ keeps the individual who experiences the event at $a_j$ in the risk set just before that event.

Under the null hypothesis of equal event-time distributions, every member of the combined risk set has the same instantaneous chance of being the next event. Conditional on one event at $a_j$ and on the two risk-set sizes,so

The quantityis the observed-minus-expected group-1 event count at time $a_j$. A positive value is local evidence that group 1 has the greater hazard function; a negative value points toward group 0.

Summing the observed-minus-expected contributions gives the unstandardized log-rank statisticUnder the null it is centered at zero. A two-sided Log-rank test compares its magnitude with the square root of its null variance.

Let $Y(t)$ be the at-risk process and $N(t)$ the counting process for observed events. Over a short interval, the multiplicative-intensity model giveswhere $h$ is the hazard function and $H$ the cumulative hazard function. Solving this relation for the infinitesimal hazard increment suggests $d\hat{H}(t)=dN(t)/Y(t)$. Summing over distinct event times gives the Nelson–Aalen estimatorwhere $d_j$ events occur among $r_j$ individuals at risk. Here there are no ties, so $d_j=1$.

Let $u_j^{(1)}=u_j$ and $u_j^{(0)}=1-u_j$. Applying the Nelson–Aalen estimator separately to group $k$ giveswith a zero contribution when the event at $a_{\ell }$ occurs in the other group.

The group-specific estimated cumulative hazard jumps only when that group experiences the event, so

Put $r_j=r_j^{(0)}+r_j^{(1)}$. The difference between the two Nelson–Aalen estimator increments isTherefore the log-rank weightsmake each summand of $T_W$ equal the corresponding summand of $T_0$, and hence $T_W=T_0$.

The variance of an estimated hazard increment is large when its group has few individuals in the risk set. The log-rank weights are near zero when either $r_j^{(0)}$ or $r_j^{(1)}$ is small and are largest when both groups retain substantial information. They therefore suppress noisy late-event comparisons and weight each observed-minus-expected event by its available information. Unit weights would instead give equal influence to unstable increments from depleted risk sets.

A proportional hazards family has hazard functions related bywhere the hazard ratio $c(z)$ is positive and independent of time. Equivalently, its cumulative hazard functions satisfy $H(t\mid z)=c(z)H_0(t)$ and its survivor functions satisfy $S(t\mid z)=S_0(t{)}^{c(z)}$.

Writing the two functions in the question as survivor functions, $S_2(t)=S_1(t{)}^{\lambda }$. Since $H_k(t)=-\log S_k(t)$,Differentiating at times where the hazard functions exist gives $h_2(t)=\lambda h_1(t)$. Their hazard ratio is therefore the constant $\lambda$, so they form a proportional hazards family.

The transformation is $T_k=a_k{U}_k^b$. Because $U_k$ has a unit-rate exponential distribution,Thus $T_k$ has a Weibull distribution, withConsequently $h_2(t)/h_1(t)=(a_1/a_2{)}^{1/b}$ is constant, proving proportional hazards.

Choose a parametric baseline hazard $h_0(t;\eta )$ and fitUnder independent right censoring, the full survival likelihood isEstimate $(\eta ,\beta )$ by maximum likelihood estimation and test $H_0:\beta =0$ with a likelihood-ratio test, Wald test, or score test. Equality of the two event-time distributions is exactly $\beta =0$ within this model.

A semiparametric proportional hazards model specifieswith finite-dimensional parameter $\beta$ but an unspecified baseline hazard $h_0$. A partial likelihood uses a component of the data likelihood that depends on $\beta$ while eliminating the nuisance function. In the Cox proportional-hazards model, conditioning on which member of each risk set experiences the event produces the Cox partial likelihood.

For each observed event $i$ with $v_i=1$, let $R_i=\{j:x_j\ge x_i\}$ be its risk set. The Cox partial likelihood isMaximize it to obtain $\hat{\beta }$, estimate its variance from the observed partial information, and test $H_0:\beta =0$ using a partial likelihood-ratio test, Wald test, or score test. A positive fitted coefficient means the group with $z=1$ has the larger hazard.

Letand let $C$ be the common partial likelihood contribution from the first $n-3$ observations. Since $t_{n-2}>x_{n-2}$, the three possible complete-data tail orderings and their partial likelihoods areTheir sum isWhen individual $n-2$ is right-censored at $x_{n-2}$, that individual leaves the risk set before the event at $x_{n-1}$, so the directly calculated Cox partial likelihood is also $C,a/(a+b)$. Summing over the unobserved compatible event orderings therefore reproduces the censored-data partial likelihood.

Empirical likelihood assigns unknown probability masses $p_j\ge 0$ to data-supported event times or intervals, imposes ${\sum }_j{p}_j=1$, and maximizes the product of each observation's probability. For event-time data, write $F(t)=\mathbb{P}(T\le t)$ and $S(t)=1-F(t)$ for the survivor function. Exact, right-censored, left-censored, interval-censored, and truncated observations contribute the probability of their respective compatible sets.