Biostatistics

Contents

• Introduction and overview
• 0. Review: statistical inference for clinical data
• 1. Clinical research questions and study designs
• 2. Randomized controlled trials: design and sample size
• 3. Therapeutic evaluation and biases in clinical research
• 4. Censored data I: survival functions, Kaplan–Meier and the log-rank test
• 5. Censored data II: the Cox model and time-dependent covariates
• 6. Validating predictive models in health science
• 7. Introduction to causal inference
• 8. Complements: further designs and survival topics
• 9. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student trained in mathematics, computer science or physics who will analyse clinical data, build a risk score or work with physicians, and who needs the statistical language in which medicine decides what works.

Empirical first: give a drug with no effect at all to the patients of a cohort who happen to start it during follow-up, then compare “users” with “non-users” from the start of follow-up. The users can die at one third of the rate of the non-users (ch. 3, Proposition 3.10), only because they had to stay alive long enough to start the drug. Mantel and Byar (1974) used heart-transplant data to show how this waiting time biases a naive comparison in favour of transplantation.

Biostatistics is the statistics of clinical research. The course follows the two goals that the official description puts at its centre, prediction and causation, and the methods it singles out as not covered elsewhere in the master: the design of studies and trials, the analysis of censored (time-to-event) data, and causal inference.

It is the statistical companion of the other courses of the Santé track, such as Deep Learning for Medical Imaging and Representation Learning for Computer Vision and Medical Imaging, whose models end up evaluated in clinical studies. It shares the bootstrap and Monte Carlo tools of Computational Statistics and the causal graphs of Introduction to Probabilistic Graphical Models and Deep Generative Models.

What you will be able to do

• Translate a clinical question into an estimand, and choose a study design that can answer it.
• Compute and interpret risk differences, risk ratios, odds ratios, hazard ratios and numbers needed to treat, with their confidence intervals.
• Plan a randomized trial: randomization scheme, endpoints, sample size for binary, continuous and time-to-event outcomes, and interim analyses.
• Recognize selection, information, confounding and time-dependent biases, and quantify them in simple models.
• Estimate survival curves from censored data with the Kaplan–Meier and Nelson–Aalen estimators, and compare groups with the log-rank test.
• Fit and check a Cox proportional hazards model, use time-dependent covariates, and estimate cumulative incidences under competing risks.
• Validate a clinical prediction model: discrimination, calibration, Brier score, optimism correction and decision curves.
• Define causal effects with potential outcomes and estimate them from observational data by standardization, inverse probability weighting and doubly robust estimators.
• Read a causal diagram, choose an adjustment set, and emulate a target trial with observational data.
• Implement these methods from scratch in NumPy, and check them on simulated data where the truth is known.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com aims at introducing both the concepts and the methods used in clinical (or medical) research, as part of the health-science theme of the master. It emphasizes the principles and concepts that underlie the different goals of clinical research, prediction and causation, and it develops specific statistical methods to plan studies and to analyse their data.

The description focuses on notions and methods that are not covered by other courses of the master, such as design, survival analysis and causal inference, treated from both the theoretical and the applied point of view and illustrated on several practical examples. It announces 7 lectures of 3 hours, each comprising formal lectures and practical work, and a final 3-hour presentation of the group home projects.

Extensions beyond the syllabus: the review ch. 0, decision curves (section 6.6), doubly robust estimation and the E-value (ch. 7), group-sequential designs and competing risks (ch. 8), and the research directions of ch. 9.

Where this course fits in the programme

Official prerequisites: there are no formal prerequisites, but undergraduate courses in probability or statistics will help. The course is a second-semester Modelling course of 5 ECTS of the Santé track, taught by Raphaël Porcher (Université Paris Cité, as listed on master-mva.com).

• Refreshers — Probability and Statistics for ch. 0 and the asymptotics used throughout; Scientific Python, NumPy and PyTorch for the labs; Linear Algebra, Analysis and Optimization for the Newton iterations of the logistic and Cox models.
• First-semester courses that prepare it — Introduction to Statistical Learning (risk, overfitting, cross-validation); Computational Statistics (bootstrap, Monte Carlo); Introduction to Probabilistic Graphical Models and Deep Generative Models (directed graphs, d-separation); Medical Image Analysis Based on Generative, Geometric and Biophysical Models (the clinical side of the Santé track).
• Related second-semester courses — Deep Learning for Medical Imaging and Representation Learning for Computer Vision and Medical Imaging (models whose validation follows ch. 6); Functional Brain Imaging and Brain–Computer Interfaces; Bayesian Machine Learning; Responsible Machine Learning (fairness and causality); Stopping Times and Online Algorithms (the stopping problems behind interim analyses); Regulating AI (medical AI as a regulated product).
• Alongside and after the course — the Hospital Immersion: Paired Project with a Physician, where the questions of ch. 1 are met in practice, and the Research Internship.

Course outline

• Review: statistical inference for clinical data — Wald intervals and the delta method, effect measures of a two-by-two table, tests and power, logistic regression.
• Clinical research questions and study designs — describing, predicting and explaining; cohort and case-control studies; phases of drug development.
• Randomized controlled trials: design and sample size — why and how to randomize, concealment and blinding, endpoints, sample size, non-inferiority, multiplicity, intention to treat.
• Therapeutic evaluation and biases in clinical research — absolute and relative effects, selection and information bias, regression to the mean, confounding and Mantel–Haenszel, immortal time.
• Censored data I — survival function and hazard, censoring, likelihood, the Kaplan–Meier estimator and Greenwood’s formula, the log-rank test.
• Censored data II — the Cox model and its partial likelihood, checking proportional hazards, time-dependent covariates and landmark analysis.
• Validating predictive models in health science — discrimination, calibration, Brier score, optimism and internal validation, net benefit.
• Introduction to causal inference — potential outcomes, standardization, inverse probability weighting, the propensity score, doubly robust estimation, causal diagrams, target trial emulation, sensitivity analysis.
• Complements: further designs and survival topics — estimands, non-compliance, group-sequential designs, cluster trials, subgroups, restricted mean survival time, sample size of survival trials, competing risks, validation of survival models (self-study).
• Research directions and open problems — principal strata, machine learning for survival data, heterogeneous treatment effects, real-world evidence, prediction models with AI, platform trials (self-study).

Session plan

27 h in class = 9 sessions of 3 h.

The volume is the Université Paris-Saclay programme table (lectures 24 h, project 3 h). master-mva.com announces 7 lectures of 3 h, each combining formal lectures and practical work, followed by a 3 h presentation of the group home projects, which makes 24 h; the plan below uses the 3 extra lecture hours of the programme table for an eighth lecture session. The split of the content into sessions is a proposal.

Level 1 = foundations, 2 = core results and methods, 3 = advanced or research-level material; the level never decreases.

Self-study, outside the 27 h (about 9 h in total): ch. 8 Complements: further designs and survival topics (≈ 4 h, with Ex. 14, 15, 17 and Codes 8.6, 8.15) · ch. 9 Research directions and open problems (≈ 3 h) · the written-exam proposal of the Projects and assessment section, as a self-test before the final exam (≈ 2 h).

How to study this module

• Read ch. 0 even if inference feels familiar: the delta method and the score equations of the logistic model come back in every later chapter.
• Redo each worked example with a pencil before running its code; then change one number (the effect, the censoring, the sample size) and predict the new output before you run it again.
• Run every code block: they use NumPy, SciPy and scikit-learn only, with fixed seeds and simulated data, and each runs in a few seconds on a laptop (about ten for the Cox code of ch. 5). The survival and causal estimators are written from scratch so that every formula of the page is visible in the code.
• Read the sections of ch. 8 as their chapters are taught: 8.1–8.5 after session 3, 8.6–8.8 after session 6, 8.9 after session 7.
• After each session, solve the exercises of the Practice column; the ★★★ exercises are good starting points for the home project.
• Keep a running list of the biases met in ch. 3: ch. 5 and ch. 7 give the analyses that remove them, and the final exam is likely to ask for that link.

0. Review: statistical inference for clinical data

Clinical research reports almost every result as an effect measure with a confidence interval or a test. This chapter recalls the M1-level tools behind those numbers, on the objects used throughout the course: proportions, two-by-two tables and the logistic model.

0.1 Wald intervals and the delta method

For a proportion \hat p = X/n with X \sim \mathrm{Bin}(n, p), the central limit theorem gives the standard error \sqrt{\hat p (1 - \hat p)/n}. Ratios are better handled on the log scale, which needs the delta method.

The multivariate version replaces g'(\theta)^2\sigma^2 by \nabla g(\theta)^\top \Sigma\, \nabla g(\theta).

0.2 Two-by-two tables and effect measures

When the event is rare in both groups, 1 - p_a \approx 1 and the odds ratio is close to the risk ratio. With frequent events they differ, and the odds ratio is always further from 1 than the risk ratio.

The interval for the odds ratio is due to Woolf (1955). The intervals are computed on the log scale and then exponentiated, so they are not symmetric around the estimate.

0.3 Tests, errors and power

The two-sided Wald test of \theta = \theta_0 at level \alpha rejects exactly when \theta_0 lies outside the Wald interval of level 1 - \alpha. This duality is used in ch. 2 for non-inferiority trials.

0.4 Likelihood and logistic regression

The first component of U(\hat\beta) = 0 reads \sum_i (y_i - \hat p_i) = 0: the fitted risks average to the observed proportion. Chapter 6 uses this identity to explain why calibration can never be judged on the development data.

import numpy as np
from scipy import stats

# 2x2 table: events / patients in each arm (Example 0.5)
a, n1 = 30, 200     # treated: 30 events out of 200
c, n0 = 50, 200     # control: 50 events out of 200
p1, p0 = a / n1, c / n0
z = stats.norm.ppf(0.975)

rd = p1 - p0
se_rd = np.sqrt(p1 * (1 - p1) / n1 + p0 * (1 - p0) / n0)
rr = p1 / p0
se_logrr = np.sqrt(1 / a - 1 / n1 + 1 / c - 1 / n0)
odds = (a * (n0 - c)) / ((n1 - a) * c)
se_logor = np.sqrt(1 / a + 1 / (n1 - a) + 1 / c + 1 / (n0 - c))
print(f"RD {rd:.3f}  [{rd - z*se_rd:.3f}, {rd + z*se_rd:.3f}]")
print(f"RR {rr:.3f}  [{rr*np.exp(-z*se_logrr):.3f}, {rr*np.exp(z*se_logrr):.3f}]")
print(f"OR {odds:.3f}  [{odds*np.exp(-z*se_logor):.3f}, {odds*np.exp(z*se_logor):.3f}]")

# The same odds ratio from a logistic regression fitted by Newton-Raphson
x = np.r_[np.ones(n1), np.zeros(n0)]
y = np.r_[np.ones(a), np.zeros(n1 - a), np.ones(c), np.zeros(n0 - c)]
X = np.c_[np.ones_like(x), x]
beta = np.zeros(2)
for _ in range(25):
    p = 1 / (1 + np.exp(-X @ beta))
    score = X.T @ (y - p)
    info = X.T @ (X * (p * (1 - p))[:, None])
    beta += np.linalg.solve(info, score)
se = np.sqrt(np.diag(np.linalg.inv(info)))
print(f"logistic: OR {np.exp(beta[1]):.3f}, SE(log OR) {se[1]:.4f} vs Woolf {se_logor:.4f}")
# RD -0.100  [-0.178, -0.022]
# RR 0.600  [0.399, 0.902]
# OR 0.529  [0.320, 0.876]
# logistic: OR 0.529, SE(log OR) 0.2567 vs Woolf 0.2567

1. Clinical research questions and study designs

A study is designed to answer a question, and the same data can answer some questions and not others. This chapter separates the kinds of questions that clinical research asks, describes the observational and experimental designs that answer them, and proves the property that makes case-control studies work.

1.1 Describe, predict, explain

The official course description puts the second and third goals at its centre: ch. 6 is about prediction, ch. 7 about causation, and the designs of ch. 1–2 serve both.

1.2 Observational designs

The theorem goes back to Cornfield (1951). When the disease is rare, the odds ratio approximates the risk ratio (ch. 0), so a case-control study estimates a risk ratio without ever following a cohort.

import numpy as np
rng = np.random.default_rng(1)

# Source population: exposure E, disease D with risk 1 % if exposed, 1/3 % if not
N = 1_000_000
E = rng.random(N) < 0.10
D = rng.random(N) < np.where(E, 0.010, 0.010 / 3)

def odds_ratio(e, d):
    a, b = np.sum(e & d), np.sum(e & ~d)
    c, dd = np.sum(~e & d), np.sum(~e & ~d)
    return (a * dd) / (b * c)

rr = D[E].mean() / D[~E].mean()
print(f"cohort: RR {rr:.2f}, OR {odds_ratio(E, D):.2f}")

# Case-control study: all cases, 4 controls per case drawn among non-cases
cases = np.flatnonzero(D)
controls = rng.choice(np.flatnonzero(~D), size=4 * cases.size, replace=False)
idx = np.r_[cases, controls]
print(f"case-control ({cases.size} cases): OR {odds_ratio(E[idx], D[idx]):.2f}")
print(f"naive 'risk' among exposed in the case-control sample: {D[idx][E[idx]].mean():.2f}")
# cohort: RR 2.97, OR 2.99
# case-control (4093 cases): OR 2.92
# naive 'risk' among exposed in the case-control sample: 0.38

The theorem assumes that controls represent the source population of the cases and that exposure is measured alike in both groups; ch. 3 shows what happens otherwise.

1.3 Experimental designs and the phases of drug development

In an experimental study the investigator assigns the intervention. A new drug goes through four phases before and after it reaches the market.

Each design comes with a reporting guideline that lists what a paper must state: CONSORT for randomized trials (Schulz et al., 2010), STROBE for observational studies (von Elm et al., 2007), and TRIPOD for prediction models (ch. 6).

1.4 Framing the question

A clinical question is framed by the PICO scheme: Population, Intervention, Comparator, Outcome. Regulatory guidance now asks for a more precise target, the estimand, which also says how events such as stopping treatment or taking a rescue medication are handled; section 8.1 presents it.

2. Randomized controlled trials: design and sample size

The randomized controlled trial is the reference design to evaluate a treatment. This chapter explains why randomization works, how it is carried out, and how the number of patients is computed. It then covers non-inferiority, repeated looks at the data, and the choice of the analysis population.

2.1 Why randomize

This formalism goes back to Neyman (1923) and Rubin (1974); ch. 7 develops it for observational data.

The difference of the sample means is therefore an unbiased estimate of the average effect. The independence holds for unmeasured characteristics too, which no adjustment of an observational study can guarantee (ch. 3 and 7).

2.2 Randomization procedures, concealment and blinding

With n = 100 and k = 5, the probability of a split at least as unbalanced as 45/55 is 0.368, both exactly and by the formula. Block randomization removes this risk, which matters for small trials and for subgroups.

When blinding is impossible, as for surgery, the outcome can still be assessed by a blinded committee (the PROBE design: prospective randomized open blinded end-point).

2.3 Endpoints

2.4 Sample size for a superiority trial

Halving the expected effect multiplies the trial by four, and the formula is only as good as the assumed effect and variance.

import numpy as np
from scipy import stats
rng = np.random.default_rng(2)

def n_per_arm(p0, p1, alpha=0.05, power=0.90):
    z = stats.norm.ppf(1 - alpha / 2) + stats.norm.ppf(power)
    return int(np.ceil(z**2 * (p0 * (1 - p0) + p1 * (1 - p1)) / (p1 - p0) ** 2))

n = n_per_arm(0.30, 0.20)
print("n per arm:", n)

# Power of the two-sided z-test by simulation (10 000 trials)
def z_test(x1, x0, n):
    p1, p0 = x1 / n, x0 / n
    se = np.sqrt(p1 * (1 - p1) / n + p0 * (1 - p0) / n)
    return (p1 - p0) / se

reps = 10_000
z_alt = z_test(rng.binomial(n, 0.20, reps), rng.binomial(n, 0.30, reps), n)
z_null = z_test(rng.binomial(n, 0.30, reps), rng.binomial(n, 0.30, reps), n)
print(f"simulated power {np.mean(np.abs(z_alt) > 1.96):.3f}, "
      f"type I error {np.mean(np.abs(z_null) > 1.96):.3f}")

# Permuted blocks of size 4, stratified by centre (3 centres, 10 patients each)
def block_list(n_patients):
    seq = []
    while len(seq) < n_patients:
        seq += list(rng.permutation(["A", "A", "B", "B"]))
    return seq[:n_patients]
for centre in range(3):
    print("centre", centre, "".join(block_list(10)))
# n per arm: 389
# simulated power 0.899, type I error 0.050
# centre 0 BBAAAABBBA
# centre 1 BABAABABAB
# centre 2 BABABABABA

2.5 Non-inferiority and equivalence

A new treatment that is cheaper, safer or easier to take need not be better: it is enough that it is not much worse.

2.6 Multiplicity and interim analyses

Every additional test is another chance of a false positive.

Holm’s step-down procedure (Holm, 1979) controls the same error rate and is uniformly more powerful.

Repeated looks at accumulating data are a form of multiplicity: trials are monitored for safety and efficacy, and stopping at the first p < 0.05 inflates the type I error. Section 8.3 computes this inflation and the group-sequential boundaries that control it.

2.7 Analysis populations and non-compliance

Only the ITT comparison inherits Theorem 2.2. The other two compare groups selected after randomization, and adherence is itself prognostic. In the Coronary Drug Project (1980), patients who took at least 80 % of their placebo had a much lower five-year mortality than those who did not.

Under non-compliance the ITT effect is diluted towards zero (section 8.2): it estimates the effect of assigning a treatment, while a valid per-protocol comparison needs the causal methods of ch. 7.

Two common variants of the parallel-group design, cluster and crossover trials, are treated in section 8.4.

3. Therapeutic evaluation and biases in clinical research

An estimate can be precise and wrong. This chapter states what a treatment effect should be measured against, defines bias, reviews its main sources, and studies in detail the two named in the official syllabus: confounding bias and time-dependent bias.

3.1 Principles of therapeutic evaluation

A treatment is evaluated for efficacy, its effect under the controlled conditions of an explanatory trial, and for effectiveness, its effect in routine care, which a pragmatic trial measures with broad eligibility and usual practice (Schwartz and Lellouch, 1967).

A fair evaluation rests on the principles of ch. 2:

• a concurrent control group, since historical controls differ in case mix and care;
• randomized, concealed allocation and, when possible, blinded assessment;
• one pre-specified, clinically relevant primary endpoint;
• the analysis of all randomized patients (intention to treat, section 2.7).

Every decision then weighs benefit against harm, on the absolute scale.

Effects reported in subgroups are fragile; section 8.5 shows how to test them.

3.2 Bias and its main sources

Epidemiology classifies biases in three families:

• selection bias, from the way subjects enter or remain in the analysis;
• information bias, from errors in measuring exposure, outcome or covariates;
• confounding, from common causes of exposure and outcome (section 3.3).

Selection bias arises when hospital controls are used in a case-control study, when patients lost to follow-up have a different prognosis, or when only “healthy users” persist with a preventive drug.

The observed difference is shrunk towards zero: with \mathrm{Se} = 0.8 and \mathrm{Sp} = 0.9, it keeps 70 % of its size. Differential errors can bias in either direction, as when cases recall past exposures better than controls (recall bias).

A third source of error mimics a treatment effect when patients are selected on an extreme measurement.

Patients enrolled because their systolic pressure is 170 mmHg, in a population with mean 140 and test–retest correlation 0.7, will have an expected pressure of 161 at the next visit, without any treatment. Only a concurrent, randomized control group separates this drop from the effect of a drug.

3.3 Confounding

In therapeutic studies the most common form is confounding by indication: physicians treat the patients who need it most, who also have the worst prognosis.

The first remedy is to compare like with like within strata of the confounder, then to pool the stratum-specific comparisons.

The estimator remains consistent when the number of strata grows with fixed small strata, such as matched pairs, where the unconditional maximum likelihood estimator is inconsistent (Breslow, 1981). On the kidney-stone data, Code 3.12 gives stratum odds ratios of 2.08 and 1.23, a Mantel–Haenszel odds ratio of 1.45 in favour of A, and a crude odds ratio of 0.75.

Regression, standardization and propensity scores (ch. 7) also adjust for confounders, but only for measured ones: residual confounding remains the main weakness of observational studies.

3.4 Time-dependent bias: immortal time

Many exposures start during follow-up, such as a drug prescribed after discharge. Classifying subjects by their exposure status at the end of follow-up biases the comparison in favour of the exposure.

When the waiting time to treatment has the same mean as the survival time (\mu = \lambda), a useless drug appears to divide mortality by three. The never-treated rate is inflated too, because those subjects are selected for dying early.

The person-time split becomes a time-dependent covariate in the Cox model of ch. 5 (Mantel and Byar, 1974); landmark analyses and target trial emulation (ch. 7) are the other remedies. Lévesque et al. (2010) show the bias at work in a cohort of statin users.

import numpy as np
rng = np.random.default_rng(4)

# Kidney stones (Charig et al., 1986): successes and patients, open surgery (A) vs PCNL (B)
strata = {"small": (81, 87, 234, 270), "large": (192, 263, 55, 80)}
num = den = 0.0
for name, (sa, na, sb, nb) in strata.items():
    a, b, c, d = sa, na - sa, sb, nb - sb      # A success, A failure, B success, B failure
    n = na + nb
    print(f"{name}: A {sa/na:.2f}, B {sb/nb:.2f}, OR {a*d/(b*c):.2f}")
    num += a * d / n
    den += b * c / n
SA, NA = 81 + 192, 87 + 263
SB, NB = 234 + 55, 270 + 80
print(f"crude: A {SA/NA:.2f}, B {SB/NB:.2f}, OR {SA*(NB-SB)/((NA-SA)*SB):.2f}")
print(f"Mantel-Haenszel OR {num/den:.2f}")

# Immortal time bias: death T ~ Exp(lam), treatment start W ~ Exp(mu), no effect at all
lam, mu, n = 0.1, 0.1, 200_000
T = rng.exponential(1 / lam, n)
W = rng.exponential(1 / mu, n)
ever = W < T
naive = (ever.sum() / T[ever].sum()) / ((~ever).sum() / T[~ever].sum())
print(f"naive 'ever vs never' rate ratio {naive:.3f}, theory {lam/(2*lam+mu):.3f}")

# Correct person-time split: untreated time [0, min(W, T)), treated time [W, T)
pt_untreated = np.minimum(W, T).sum()
pt_treated = np.clip(T - W, 0, None).sum()
d_untreated, d_treated = (~ever).sum(), ever.sum()
print(f"person-time rate ratio {(d_treated/pt_treated)/(d_untreated/pt_untreated):.3f}")
# small: A 0.93, B 0.87, OR 2.08
# large: A 0.73, B 0.69, OR 1.23
# crude: A 0.78, B 0.83, OR 0.75
# Mantel-Haenszel OR 1.45
# naive 'ever vs never' rate ratio 0.334, theory 0.333
# person-time rate ratio 0.997

4. Censored data I: survival functions, Kaplan–Meier and the log-rank test

In many clinical studies the outcome is a time: time to death, to relapse, to discharge from hospital. For many patients follow-up ends before the event, and their time is only known to exceed the follow-up: it is censored. This chapter defines the objects that describe a time to event and estimates them without any parametric model.

4.1 Survival function and hazard

4.2 Censoring

Censoring comes from the end of the study (administrative censoring), from losses to follow-up and from withdrawals. Independence is plausible for the first and doubtful when patients leave because they feel worse.

A related phenomenon is left truncation: subjects enter observation at a time V > 0, for instance at the date of referral, and are seen only if T > V. Risk sets then include a subject only after its entry time, as in the (\text{start}, \text{stop}] format of ch. 5.

4.3 Likelihood under right censoring

For the exponential law, the log-likelihood is D\log\lambda - \lambda \sum_i x_i with D = \sum_i \delta_i the number of events, so that

\hat\lambda = \frac{D}{\sum_i x_i} = \frac{\text{number of events}}{\text{person-time}} .

The observed information is D/\hat\lambda^2, so \mathrm{se}(\log \hat\lambda) = 1/\sqrt{D}: the precision depends on the number of events, not on the number of patients. For the ten patients of Example 4.9, D = 5 over 78 months of follow-up gives \hat\lambda = 0.064 per month and an estimated median survival of \log 2/\hat\lambda = 10.8 months.

4.4 The Kaplan–Meier estimator

The estimator applies Proposition 4.2 in discrete time: each conditional probability of surviving an event time is estimated by the proportion of survivors among those at risk. Censored subjects count in the risk sets until they leave, and independent censoring makes those at risk representative. The curve drops only at event times; it is the nonparametric maximum likelihood estimator.

Intervals computed on the \log(-\log) scale stay within [0, 1] and have better coverage than the plain Wald interval.

4.5 Comparing two groups: the log-rank test

The test is due to Mantel (1966) and Peto and Peto (1972). It is the locally most powerful rank test against proportional-hazards alternatives (ch. 5) and loses power when the survival curves cross. Weighted versions put more weight on early differences (Gehan–Breslow) or, with suitable parameters, late ones (the Fleming–Harrington family), and the stratified log-rank test sums O - E and V over strata of a prognostic factor.

When hazards are not proportional, a single hazard ratio summarizes a difference poorly; the restricted mean survival time of section 8.6, the area under the survival curve up to a horizon, is then preferred.

import numpy as np
from scipy import stats
rng = np.random.default_rng(5)

def kaplan_meier(time, event):
    """Event times, KM survival and Greenwood standard error."""
    t = np.unique(time[event == 1])
    n = np.array([(time >= s).sum() for s in t])           # at risk just before s
    d = np.array([((time == s) & (event == 1)).sum() for s in t])
    S = np.cumprod(1 - d / n)
    with np.errstate(divide="ignore", invalid="ignore"):
        gw = np.cumsum(d / (n * (n - d)))
    return t, S, S * np.sqrt(gw)

def logrank(time, event, group):
    O = E = V = 0.0
    for s in np.unique(time[event == 1]):
        at_risk = time >= s
        n, n1 = at_risk.sum(), (at_risk & (group == 1)).sum()
        d = ((time == s) & (event == 1)).sum()
        d1 = ((time == s) & (event == 1) & (group == 1)).sum()
        O += d1
        E += d * n1 / n
        if n > 1:
            V += d * (n1 / n) * (1 - n1 / n) * (n - d) / (n - 1)
    chi2 = (O - E) ** 2 / V
    return O, E, chi2, stats.chi2.sf(chi2, 1)

# Example 4.9 by hand: months, 1 = death, 0 = censored
time = np.array([3, 4, 5, 5, 6, 8, 9, 11, 12, 15])
event = np.array([1, 0, 1, 1, 0, 1, 0, 1, 0, 0])
for s, S, se in zip(*kaplan_meier(time, event)):
    print(f"t={s:2d}  S={S:.3f}  se={se:.3f}")

# Simulated trial: exponential survival, hazard ratio 0.6, uniform censoring
n = 150
group = np.r_[np.zeros(n), np.ones(n)]
T = rng.exponential(1 / np.where(group == 1, 0.06, 0.10))
C = rng.uniform(0, 24, 2 * n)
X, delta = np.minimum(T, C), (T <= C).astype(int)
t, S, se = kaplan_meier(X[group == 0], delta[group == 0])
k = np.searchsorted(t, 12, side="right") - 1
print(f"control S(12) = {S[k]:.3f} +/- {1.96*se[k]:.3f}; truth {np.exp(-0.10*12):.3f}")
O, E, chi2, p = logrank(X, delta, group)
print(f"log-rank: O={O:.0f}, E={E:.1f}, chi2={chi2:.2f}, p={p:.1e}")
# t= 3  S=0.900  se=0.095
# t= 5  S=0.675  se=0.155
# t= 8  S=0.540  se=0.173
# t=11  S=0.360  se=0.187
# control S(12) = 0.268 +/- 0.089; truth 0.301
# log-rank: O=60, E=85.5, chi2=17.72, p=2.6e-05

5. Censored data II: the Cox model and time-dependent covariates

The log-rank test compares groups; clinical studies need to estimate an effect while accounting for several covariates, some of which change over time. The Cox model does this without specifying the baseline hazard. This chapter derives its partial likelihood, shows how to check its main assumption, and uses it to remove the immortal-time bias of ch. 3.

5.1 Proportional hazards

The ratio of the hazards of two subjects does not depend on time: this is the proportional hazards assumption. For a binary treatment indicator x_j, e^{\beta_j} is the hazard ratio of treated versus untreated patients with the same other covariates.

5.2 The partial likelihood

Let R(t) = \{j : X_j \ge t\} be the risk set at t, and assume first that event times are distinct.

With tied event times, Breslow’s approximation lets every tied event use the full risk set, which is what Code 5.7 does; Efron’s approximation is more accurate when ties are many.

With \hat\beta = 0 the Breslow estimator is the Nelson–Aalen estimator. Predicted survival curves are what a prognostic model reports, and section 8.9 shows how they are validated.

5.3 Checking proportional hazards

• Graphically, curves of \log(-\log \hat S_g(t)) against \log t should be parallel across the groups g.
• With residuals, the Schoenfeld residual of an event i is r_i = x_i - \bar x(\hat\beta, x_i). If the true coefficient varies with time, the residuals scaled by D\, I(\hat\beta)^{-1}, with D the number of events, have a mean close to \beta(t) - \hat\beta; Grambsch and Therneau (1994) turn this into a test.
• Remedies include stratifying the baseline hazard on the offending variable, letting the coefficient change over time through an interaction with a function of t, and reporting restricted mean survival times (section 8.6).

5.4 Time-dependent covariates

The same format handles left truncation, since a subject enters the risk sets only after its first start time.

A landmark analysis (Anderson, Cain and Gelber, 1983) is the simple alternative: choose a time s, classify subjects by their treatment status at s, and analyse survival from s among those still alive.

import numpy as np
rng = np.random.default_rng(7)

def cox(start, stop, event, X, iters=30):
    """Cox partial likelihood (Breslow ties), counting-process data (start, stop]."""
    X = np.atleast_2d(X.T).T
    beta = np.zeros(X.shape[1])
    times = np.unique(stop[event == 1])
    for _ in range(iters):
        U, I = np.zeros_like(beta), np.zeros((beta.size, beta.size))
        for t in times:
            R = (start < t) & (stop >= t)                    # risk set at t
            D = (stop == t) & (event == 1)                   # failures at t
            w = np.exp(X[R] @ beta)
            xbar = w @ X[R] / w.sum()
            cov = (X[R] * w[:, None]).T @ X[R] / w.sum() - np.outer(xbar, xbar)
            U += X[D].sum(0) - D.sum() * xbar
            I += D.sum() * cov
        step = np.linalg.solve(I, U)
        beta += step
        if np.abs(step).max() < 1e-9:
            break
    return beta, np.sqrt(np.diag(np.linalg.inv(I)))

# (a) Recovering known log hazard ratios: h(t|x) = 0.1 exp(0.7 x1 - 0.5 x2)
n = 1000
X = np.c_[rng.integers(0, 2, n), rng.standard_normal(n)]
T = rng.exponential(1 / (0.1 * np.exp(X @ np.array([0.7, -0.5]))))
C = rng.uniform(0, 15, n)
b, se = cox(np.zeros(n), np.minimum(T, C), (T <= C).astype(int), X)
print("beta", b.round(3), "se", se.round(3))

# (b) Immortal time (ch. 3): useless treatment started at W ~ Exp(0.1), death T ~ Exp(0.1)
n = 4000
T, W = rng.exponential(10, n), rng.exponential(10, n)
ever = (W < T).astype(float)
b_naive, _ = cox(np.zeros(n), T, np.ones(n, int), ever)
# split treated subjects into (0, W] untreated and (W, T] treated
s = np.r_[np.zeros(n), W[W < T]]
e = np.r_[np.minimum(W, T), T[W < T]]
d = np.r_[(T <= W).astype(int), np.ones(int(ever.sum()), int)]
z = np.r_[np.zeros(n), np.ones(int(ever.sum()))]
b_td, se_td = cox(s, e, d, z)
lo, hi = np.exp(b_td[0] + np.array([-1.96, 1.96]) * se_td[0])
print(f"HR naive {np.exp(b_naive[0]):.2f}; time-dependent {np.exp(b_td[0]):.2f} [{lo:.2f}, {hi:.2f}]")
# beta [ 0.714 -0.523] se [0.086 0.045]
# HR naive 0.28; time-dependent 0.93 [0.86, 1.01]

6. Validating predictive models in health science

A prediction model is useful only if its predictions are accurate for new patients. This chapter defines what “accurate” means for a risk model (discrimination, calibration, overall error and clinical usefulness), explains why performance measured on the development data is optimistic, and shows how to correct it.

6.1 Prediction models and their life cycle

A model goes through development, internal validation on the development population, external validation on patients from another time or place, and ideally an impact study. TRIPOD (Collins et al., 2015) and TRIPOD+AI (Collins et al., 2024) list what a report must contain. Performance has four aspects (Steyerberg et al., 2010), taken in turn below.

6.2 Discrimination

The AUC depends only on the ranks of the predictions, so it is unchanged by any increasing transformation of \hat p: a model can rank patients exactly as the true risk does, with the same AUC, and still predict odds that are all ten times too high.

For time-to-event outcomes, Harrell’s C (Harrell et al., 1982) is the proportion of usable pairs, those whose shorter time is an event, in which the subject who fails first has the higher predicted risk. It depends on censoring, which Uno et al. (2011) correct by inverse probability of censoring weights.

6.3 Calibration

Calibration can therefore only be assessed on new data or by resampling. A calibration plot shows the observed proportion against the predicted risk, by groups or with a smoother; Van Calster et al. (2019) call calibration “the Achilles heel of predictive analytics”.

6.4 Overall performance: the Brier score

The second term is irreducible: it depends on how predictable the outcome is. Predicting the prevalence \bar y for everyone scores \bar y(1 - \bar y), and the scaled Brier score is 1 - \frac{\mathrm{BS}}{\bar y\,(1 - \bar y)} .

6.5 Overfitting, optimism and internal validation

A single train–test split wastes data and gives unstable estimates; the bootstrap or repeated cross-validation is preferred (Steyerberg et al., 2001). The optimism grows with the number of candidate predictors per event. The rule of ten events per variable (Peduzzi et al., 1996) has given way to sample-size calculations that target, among other criteria, an expected calibration slope of at least 0.9 (Riley et al., 2020). Shrinkage, uniform (Van Houwelingen and Le Cessie, 1990) or penalized, reduces overfitting.

6.6 Clinical usefulness: net benefit

A model is used to decide whom to treat, and discrimination and calibration do not say whether these decisions improve.

In Code 6.14 the prevalence is 0.32. At t = 0.1 the model’s net benefit (0.232) is below that of treating everyone (0.240), so the model is useless there. At t = 0.2 it beats both alternatives (0.166 against 0.145 and 0).

Validation on patients from another setting (external validation), the updating of a model for a new setting, and the measures used for survival models are presented in section 8.9.

import numpy as np
from sklearn.metrics import roc_auc_score
rng = np.random.default_rng(9)

def fit_logistic(X, y, iters=25):
    """Maximum likelihood logistic regression with intercept (Newton-Raphson)."""
    Z = np.c_[np.ones(len(y)), X]
    b = np.zeros(Z.shape[1])
    for _ in range(iters):
        p = 1 / (1 + np.exp(-Z @ b))
        b += np.linalg.solve(Z.T @ (Z * (p * (1 - p))[:, None]), Z.T @ (y - p))
    return b

def predict(b, X):
    return 1 / (1 + np.exp(-(b[0] + X @ b[1:])))

beta = np.r_[0.8, -0.6, 0.5, 0.4, np.zeros(16)]          # 4 real predictors, 16 noise
def draw(n):
    X = rng.standard_normal((n, beta.size))
    return X, (rng.random(n) < 1 / (1 + np.exp(-(-1.0 + X @ beta)))).astype(int)

X, y = draw(200)                                          # development sample
Xv, yv = draw(20_000)                                     # large external sample
b = fit_logistic(X, y)
auc_app = roc_auc_score(y, predict(b, X))
print(f"events {y.sum()}, candidate predictors {X.shape[1]}")

# Bootstrap optimism correction (Efron 1983; Harrell et al. 1996)
optimism = []
for _ in range(200):
    i = rng.integers(0, len(y), len(y))
    bb = fit_logistic(X[i], y[i])
    optimism.append(roc_auc_score(y[i], predict(bb, X[i])) - roc_auc_score(y, predict(bb, X)))
p = predict(b, Xv)
slope = fit_logistic(np.log(p / (1 - p))[:, None], yv)[1]   # logistic recalibration
print(f"AUC apparent {auc_app:.3f}, optimism-corrected {auc_app - np.mean(optimism):.3f}, "
      f"external {roc_auc_score(yv, p):.3f}")
print(f"external: calibration slope {slope:.2f}, mean risk {p.mean():.3f} "
      f"vs observed {yv.mean():.3f}, Brier {np.mean((p - yv) ** 2):.3f}")

# Net benefit at threshold t (Vickers and Elkin, 2006)
for t in (0.1, 0.2, 0.4):
    tp, fp = np.mean((p >= t) & (yv == 1)), np.mean((p >= t) & (yv == 0))
    nb_all = yv.mean() - (1 - yv.mean()) * t / (1 - t)
    print(f"t={t}: net benefit model {tp - fp * t / (1 - t):.3f}, treat all {nb_all:.3f}")
# events 70, candidate predictors 20
# AUC apparent 0.857, optimism-corrected 0.790, external 0.724
# external: calibration slope 0.46, mean risk 0.341 vs observed 0.316, Brier 0.203
# t=0.1: net benefit model 0.232, treat all 0.240
# t=0.2: net benefit model 0.166, treat all 0.145
# t=0.4: net benefit model 0.064, treat all -0.140

7. Introduction to causal inference

Randomization answers causal questions by design (ch. 2). When a trial is impossible, the effect of a treatment must be estimated from observational data, where treated and untreated patients differ. This chapter states the assumptions that make this possible and the estimators that follow.

7.1 Causal estimands

Recall the potential outcomes Y(1), Y(0) of Definition 2.1 and the consistency relation Y = Y(A), which also assumes that one patient’s treatment does not affect another’s outcome (no interference).

These are marginal (population-average) effects. Conditional effects, given covariates, can differ from them even without confounding when the measure is an odds ratio or a hazard ratio, which are not collapsible (Exercise 16).

7.2 Identification by standardization

The formula, due to Robins (1986), standardizes the outcome to the covariate distribution of the whole population. In practice one fits an outcome model \hat m_a(l) and averages \hat m_1(L_i) - \hat m_0(L_i) over the sample (g-computation).

7.3 Inverse probability weighting

Weighting each patient by the inverse probability of the treatment received creates a pseudo-population in which treatment is independent of L. The Horvitz–Thompson (1952) form divides by n, the more stable Hájek form by the sum of the weights. Near-violations of positivity produce huge weights, which stabilization or truncation tame at the price of some bias.

7.4 The propensity score

A scalar score thus replaces a high-dimensional L for matching, stratification or weighting. The fitted model is judged by the balance it achieves, measured by standardized mean differences (below 0.1), not by its AUC (Austin, 2011).

7.5 Doubly robust estimation

The bias is a product of the two errors. This is why the estimator tolerates flexible, slowly converging machine-learning fits of both models, the starting point of double machine learning (ch. 9).

7.6 Causal diagrams

The criterion explains the checklist of Definition 3.5, and it explains two errors: adjusting for a mediator (the LDL cholesterol of Example 1.2) removes part of the effect, and adjusting for a collider creates an association.

7.7 Target trial emulation

Aligning eligibility, treatment assignment and the start of follow-up at time zero removes the immortal time of ch. 3 by construction (Hernán et al., 2016). When strategies differ only after time zero, such as “start within six months”, patients are cloned into every compatible strategy, censored when they deviate, and reweighted. The emulation still assumes no unmeasured confounding.

The E-value of VanderWeele and Ding (2017) measures how strong such confounding would have to be. For an observed risk ratio \mathrm{RR} > 1, E = \mathrm{RR} + \sqrt{\mathrm{RR}\,(\mathrm{RR} - 1)} is the minimum risk ratio that an unmeasured confounder would need with both treatment and outcome to explain the estimate away. For \mathrm{RR} = 2, E = 3.41; a protective ratio is first inverted.

import numpy as np
rng = np.random.default_rng(10)

def fit_logistic(Z, y, iters=25):
    b = np.zeros(Z.shape[1])
    for _ in range(iters):
        p = 1 / (1 + np.exp(-Z @ b))
        b += np.linalg.solve(Z.T @ (Z * (p * (1 - p))[:, None]), Z.T @ (y - p))
    return 1 / (1 + np.exp(-Z @ b))

# Confounders L1, L2; treatment depends on both; true average treatment effect = 2
n = 50_000
L1, L2 = rng.standard_normal(n), rng.standard_normal(n)
A = (rng.random(n) < 1 / (1 + np.exp(-(-1 + 1.2 * L1 - 0.6 * L2)))).astype(float)
Y = 2 * A + 2 * np.exp(L1) + L2 + rng.standard_normal(n)
one = np.ones(n)

def outcome_model(features):          # linear regression of Y on (1, A, features)
    F = np.c_[one, A, features]
    coef = np.linalg.lstsq(F, Y, rcond=None)[0]
    m1 = np.c_[one, one, features] @ coef
    m0 = np.c_[one, 0 * one, features] @ coef
    return m1, m0

def aipw(m1, m0, e):
    return np.mean(m1 - m0 + A * (Y - m1) / e - (1 - A) * (Y - m0) / (1 - e))

e_good = fit_logistic(np.c_[one, L1, L2], A)          # correct propensity model
e_bad = fit_logistic(np.c_[one, L2], A)               # omits L1
m1_bad, m0_bad = outcome_model(np.c_[L1, L2])         # linear in L1: misspecified
m1_good, m0_good = outcome_model(np.c_[np.exp(L1), L2])

w1, w0 = A / e_good, (1 - A) / (1 - e_good)
print(f"naive difference         {Y[A == 1].mean() - Y[A == 0].mean():.3f}")
print(f"g-formula, wrong model   {np.mean(m1_bad - m0_bad):.3f}")
print(f"g-formula, right model   {np.mean(m1_good - m0_good):.3f}")
print(f"IPW (Hajek)              {(w1 @ Y) / w1.sum() - (w0 @ Y) / w0.sum():.3f}")
print(f"AIPW, wrong outcome      {aipw(m1_bad, m0_bad, e_good):.3f}")
print(f"AIPW, wrong propensity   {aipw(m1_good, m0_good, e_bad):.3f}")
print(f"AIPW, both wrong         {aipw(m1_bad, m0_bad, e_bad):.3f}")
# naive difference         4.799
# g-formula, wrong model   2.312
# g-formula, right model   1.997
# IPW (Hajek)              2.001
# AIPW, wrong outcome      2.036
# AIPW, wrong propensity   1.996
# AIPW, both wrong         2.364

8. Complements: further designs and survival topics

This chapter is for self-study. It completes the design chapters with estimands, non-compliance, interim analyses, cluster trials and subgroups, and the survival chapters with the restricted mean survival time, the sample size of survival trials, competing risks and the validation of survival models. Each section relies only on the chapters it cites.

8.1 Estimands and intercurrent events

The PICO question of section 1.4 leaves open what happens when patients stop their treatment, take a rescue medication or die before the endpoint is measured. The ICH E9(R1) addendum (2019) asks for an explicit estimand.

The addendum names five strategies for intercurrent events: treatment policy, hypothetical, composite, while on treatment, and principal stratum.

8.2 Non-compliance and the complier average effect

Section 2.7 noted that non-compliance dilutes the intention-to-treat effect. Let Z be the randomized arm and D the treatment received, and suppose that controls cannot obtain the treatment.

Dividing the ITT effect by \pi_c estimates the effect among compliers. This is the instrumental-variable estimator, with randomization as the instrument (Angrist, Imbens and Rubin, 1996).

8.3 Interim analyses and group-sequential designs

Section 2.6 noted that repeated looks at accumulating data are a form of multiplicity: stopping a trial at the first p < 0.05 inflates the type I error.

In practice an independent data monitoring committee performs the interim analyses, so that the investigators stay blind to the accumulating results.

import numpy as np
from scipy import optimize

# K equally spaced looks at the accumulating z-statistic of a trial under H0.
# Z_k = S_k / sqrt(k) where S_k is a sum of k independent N(0, 1) increments.
def crossing_prob(bounds, reps=400_000, seed=3):
    rng = np.random.default_rng(seed)
    K = len(bounds)
    S = np.cumsum(rng.standard_normal((reps, K)), axis=1)
    Z = S / np.sqrt(np.arange(1, K + 1))
    return np.mean((np.abs(Z) > np.asarray(bounds)).any(axis=1))

for K in (1, 2, 3, 5):
    print(f"{K} looks at 1.96: type I error {crossing_prob([1.96] * K):.3f}")

# O'Brien-Fleming boundaries c * sqrt(K / k), with c tuned to a 5 % type I error
K = 5
f = lambda c: crossing_prob(c * np.sqrt(K / np.arange(1, K + 1))) - 0.05
c = optimize.brentq(f, 1.8, 2.4, xtol=1e-3)
print("O'Brien-Fleming, 5 looks:", np.round(c * np.sqrt(K / np.arange(1, K + 1)), 2))
# 1 looks at 1.96: type I error 0.050
# 2 looks at 1.96: type I error 0.083
# 3 looks at 1.96: type I error 0.107
# 5 looks at 1.96: type I error 0.141
# O'Brien-Fleming, 5 looks: [4.56 3.23 2.63 2.28 2.04]

With five looks, the final boundary of 2.04 is barely above 1.96: early stopping costs little power, which is why O’Brien–Fleming boundaries are popular.

8.4 Cluster and crossover trials

Some interventions act on groups: a hygiene programme in a hospital ward, a training of general practitioners. Such trials randomize clusters, and patients of the same cluster resemble each other.

The factor 1 + (m-1)\rho multiplies the sample size of Theorem 2.7. With clusters of m = 20 and \rho = 0.05, it equals 1.95: the trial of Example 2.8 needs 759 patients per arm instead of 389.

In a crossover trial, each patient receives both treatments in random order, separated by a washout period. Each patient is his or her own control, which removes the between-patient variance. The design suits chronic, stable conditions and treatments without lasting (carry-over) effects.

8.5 Subgroup analyses

Trials often report effects in subgroups (men and women, age classes). These analyses are fragile.

The proper analysis tests the interaction, the difference between subgroup effects, for a few pre-specified subgroups. In the ISIS-2 trial (ISIS-2 Collaborative Group, 1988), aspirin appeared useless for patients born under Gemini or Libra, a result reported precisely to show how chance produces such findings.

8.6 Restricted mean survival time

When hazards are not proportional, a single hazard ratio summarizes the difference poorly. A summary that is always interpretable is the area under the survival curve up to a horizon.

The RMST is estimated by the area under the Kaplan–Meier curve. The difference between two arms is the mean event-free time gained over the horizon, in months, a quantity that patients understand (Royston and Parmar, 2013).

8.7 Sample size of a survival trial

For a time-to-event endpoint, the information in a trial is carried by the events, not by the patients.

The expected proportion of events depends on the accrual rate, the follow-up and the control survival; trials therefore fix the number of events and analyse when it is reached.

8.8 Competing risks

A patient waiting for a transplant can die, or receive the transplant; an elderly patient followed for cancer recurrence can die of another cause. When one type of event prevents the others, the events compete.

The quantity 1 - e^{-H_1(t)} describes a hypothetical world in which the competing event has been abolished, not the risk that patients actually face. The real-world risk is estimated by the Aalen–Johansen estimator (Aalen and Johansen, 1978), \hat F_j(t) = \sum_{i:\ t_i \le t} \hat S(t_i^-)\, \frac{d_{j,i}}{n_i}, where \hat S is the all-cause Kaplan–Meier estimator and d_{j,i} counts the events of cause j at t_i.

Two regression models answer two questions. The cause-specific Cox model, a Cox model for h_j with competing events censored, describes the aetiology. The Fine–Gray model (Fine and Gray, 1999) puts proportional hazards on the subdistribution hazard of F_j and is used to predict cumulative incidences (Putter, Fiocco and Geskus, 2007; Austin, Lee and Fine, 2016).

import numpy as np
rng = np.random.default_rng(8)

# Two competing causes with constant cause-specific hazards, plus censoring
h1, h2, n = 0.05, 0.10, 5000
T = rng.exponential(1 / (h1 + h2), n)
cause = np.where(rng.random(n) < h1 / (h1 + h2), 1, 2)
C = rng.uniform(0, 30, n)
X = np.minimum(T, C)
status = np.where(T <= C, cause, 0)                 # 0 = censored

order = np.argsort(X)
X, status = X[order], status[order]
at_risk = n - np.arange(n)                          # no ties with continuous times
S_all = np.cumprod(1 - (status > 0) / at_risk)      # all-cause KM
S_prev = np.r_[1.0, S_all[:-1]]
cif1 = np.cumsum(S_prev * (status == 1) / at_risk)  # Aalen-Johansen, cause 1
one_minus_km = 1 - np.cumprod(1 - (status == 1) / at_risk)

k = np.searchsorted(X, 10.0) - 1
true_cif = h1 / (h1 + h2) * (1 - np.exp(-(h1 + h2) * 10))
print(f"t = 10: 1-KM {one_minus_km[k]:.3f} (target {1-np.exp(-h1*10):.3f}), "
      f"Aalen-Johansen {cif1[k]:.3f} (true CIF {true_cif:.3f})")
# t = 10: 1-KM 0.387 (target 0.393), Aalen-Johansen 0.253 (true CIF 0.259)

With constant hazards 0.05 and 0.10, the true cumulative incidence of cause 1 at 10 years is F_1(10) = \tfrac13\big(1 - e^{-1.5}\big) = 0.259, while one minus Kaplan–Meier targets 1 - e^{-0.5} = 0.393: a 50 % overestimate.

8.9 Validating survival models, external validation and updating

A prognostic model for a time to event predicts \hat S(\tau \mid x) at a horizon \tau, for instance with the Breslow estimator of Definition 5.4. Its validation adapts the measures of ch. 6 to censoring.

• Discrimination: Harrell’s or Uno’s C index (section 6.2), or the time-dependent AUC at \tau, which compares subjects with an event before \tau to those still event-free at \tau (Heagerty, Lumley and Pepe, 2000; Blanche, Dartigues and Jacqmin-Gadda, 2013).
• Calibration: within groups of predicted risk, compare the mean of 1 - \hat S(\tau \mid x) with the Kaplan–Meier estimate of the observed risk at \tau.
• Overall error: the Brier score at \tau with inverse probability of censoring weights (Graf et al., 1999), defined below.

On patients from another time, place or setting, the case mix changes. A narrower range of patients lowers the AUC even for a correct model, and a different baseline risk shifts all predictions. Models are therefore updated in increasing order of complexity: the intercept (or the baseline hazard), then the calibration slope of Definition 6.5, then selected coefficients (Steyerberg, 2019). Royston and Altman (2013) describe these steps for Cox models.

9. Research directions and open problems

This chapter is for self-study. It surveys questions that are open at the time of writing, with entry points into the literature; each section ends with a question that could start a project or an internship.

9.1 Estimands and principal strata

The ICH E9(R1) addendum (section 8.1) made intercurrent events explicit, but several of its strategies target quantities that are not identified without strong assumptions. The principal stratum strategy, for instance, targets the effect among patients who would not experience the intercurrent event under either treatment (Frangakis and Rubin, 2002). Membership of such a stratum is never observed, so estimation relies on untestable assumptions or on sensitivity analyses.

Open question: which estimands for trials with frequent treatment discontinuation and deaths are both clinically meaningful and estimable under assumptions that can be checked?

9.2 Machine learning for time-to-event data

Random survival forests (Ishwaran et al., 2008) relax both the linearity and the proportional hazards assumptions of ch. 5. Neural extensions of the Cox model such as DeepSurv (Katzman et al., 2018) replace the linear predictor by a network but keep proportional hazards. On low-dimensional clinical data, a systematic review by Christodoulou et al. (2019) found no performance benefit of machine learning over logistic regression for binary outcomes, and the gains reported for survival are similarly debated.

The censoring mechanism complicates both training and evaluation. Conformal methods give distribution-free lower prediction bounds on survival times under censoring (Candès, Lei and Ren, 2023), and the IPCW measures of section 8.9 remain the standard for evaluation.

Open question: when do flexible survival models improve calibrated individual predictions, and how should their uncertainty be reported to clinicians?

9.3 Heterogeneous treatment effects

The average effect of ch. 7 hides differences between patients. Causal forests (Wager and Athey, 2018), meta-learners (Künzel et al., 2019) and double machine learning with cross-fitting (Chernozhukov et al., 2018) estimate conditional average treatment effects; Proposition 7.7 explains why products of slowly converging errors remain acceptable.

The PATH statement (Kent et al., 2020) recommends a cautious, risk-based approach in trials: estimate a baseline risk model first, then model the treatment effect across risk strata. Validating a model of individual effects is hard, since individual effects are never observed.

Open question: how can a model of treatment-effect heterogeneity be validated externally, with calibration and discrimination measures analogous to those of ch. 6?

9.4 Real-world evidence and transportability

Target trial emulation (section 7.7) is now applied to large health-care databases. The RCT DUPLICATE initiative compared emulations with the trials they mimicked and reported results that agree for some questions and not for others (Franklin et al., 2021).

When trial results must be carried to a different target population, transportability methods reweight or standardize the trial to the covariate distribution of that population (Dahabreh et al., 2020). Both lines rest on exchangeability assumptions that sensitivity analyses, such as the E-value, can only partly probe.

Open question: which features of a clinical question and of a database predict whether an emulated target trial will reproduce the randomized result?

9.5 Prediction models with artificial intelligence

A review of prediction models for COVID-19 (Wynants et al., 2020) found that almost all were at high risk of bias, mostly because of small samples, poor handling of missing data and no proper validation, the failures that ch. 6 warns against. TRIPOD+AI (Collins et al., 2024) extends the reporting rules to machine-learning models.

Two further problems are specific to deployment. Dataset shift between hospitals or over time degrades calibration first. And a model can be accurate while encoding inequities: Obermeyer et al. (2019) showed that an algorithm that used health-care costs as a proxy for health needs underestimated the needs of Black patients.

Open question: how should a deployed clinical prediction model be monitored and updated so that its calibration stays acceptable in every patient group?

9.6 Adaptive and platform trials

Group-sequential designs (section 8.3) are the simplest adaptive designs. Platform trials evaluate several treatments under a master protocol, adding and dropping arms as evidence accrues (Woodcock and LaVange, 2017). The RECOVERY trial used this design during the COVID-19 pandemic and showed that dexamethasone reduced mortality among hospitalized patients receiving respiratory support (RECOVERY Collaborative Group, 2021). Bayesian adaptive methods (Berry et al., 2010) adjust allocation probabilities or stopping rules as the posterior evolves.

Open question: how should type I error, time trends in the control arm and the sharing of controls between arms be handled when a platform trial runs for years?

Exercises

The exercises are sorted by level: ★ applies the course, ★★ combines results or proves a short statement, ★★★ is open-ended. The Practice column of the session plan says when to attempt each one.

1. ★ Effect measures from a trial. A trial records 12 deaths among 150 treated patients and 24 among 150 controls. Compute the risk difference, the risk ratio and the odds ratio with their 95 % confidence intervals, and the number needed to treat. Which measure would you report to a patient?

   Solution: \hat p_1 = 0.08, \hat p_0 = 0.16. The risk difference is -0.08, 95 % CI [-0.153, -0.007], so the NNT is 12.5: about 13 patients treated to avoid one death.

   By Proposition 0.4, the risk ratio is 0.50, CI [0.26, 0.96], and the odds ratio is 0.457, CI [0.219, 0.951]. For a patient, the absolute measure (8 fewer deaths per 100 treated) is the most informative.

2. ★ What a case-control study can estimate. A case-control study includes 200 cases of a rare cancer and 400 controls; 80 cases and 80 controls were exposed to a solvent. (a) Estimate the odds ratio with its 95 % CI. (b) Why can the risk of cancer among the exposed not be estimated? (c) What does the odds ratio approximate, and why?

   Solution: (a) \mathrm{OR} = (80 \times 320)/(120 \times 80) = 2.67; the variance of its log is 1/80 + 1/120 + 1/80 + 1/320, which gives the CI [1.83, 3.88].

   2. The proportion of cases in the sample, 1/3, was fixed by the investigators. (c) By Theorem 1.4 the sample odds ratio equals the population odds ratio, which approximates the risk ratio because the cancer is rare.

3. ★ Sample size for a continuous endpoint. A trial must detect a reduction of 0.4 percentage points of HbA1c, with a standard deviation of 1.1 in each arm, two-sided \alpha = 0.05 and power 0.80. (a) How many patients per arm? (b) With 15 % dropout? (c) And for a reduction of 0.3 points?

   Solution: (a) Theorem 2.7: n = 2.802^2 \times 2 \times 1.21 / 0.16 = 118.7, so 119 per arm. (b) 118.7/0.85 gives 140 per arm. (c) The factor (0.4/0.3)^2 = 1.78 gives 212 per arm.

4. ★ Simpson’s paradox in an observational study. Among patients with a mild disease, 20/200 treated and 60/400 untreated died; among patients with a severe disease, 120/400 treated and 50/100 untreated died. (a) Compute the crude odds ratio of death, treated versus untreated. (b) Compute the stratum-specific odds ratios and the Mantel–Haenszel odds ratio. (c) Explain the difference.

   Solution: (a) Crude: 140/600 against 110/500, odds ratio 1.08, so the treatment looks slightly harmful.

   2. Mild: \frac{20 \times 340}{180 \times 60} = 0.63; severe: \frac{120 \times 50}{280 \times 50} = 0.43. Mantel–Haenszel: \frac{11.33 + 12}{18 + 28} = 0.51.

   3. Severe patients are treated more often (80 % against 33 %) and die more: severity confounds the crude comparison (confounding by indication, section 3.3).

5. ★ Kaplan–Meier and log-rank by hand. Group A has times 2, 4+, 6, 8, 10+ and group B has times 1, 3, 3, 5+, 7 (months, + for censored). (a) Compute the two Kaplan–Meier curves and medians. (b) Compute O_B, E_B, V and the log-rank statistic.

   Solution: (a) A: \hat S = 0.8 at 2, 0.533 at 6, 0.267 at 8; median 8 months. B: 0.8 at 1, 0.4 at 3, 0 at 7; median 3 months.

   2. The event times 1, 2, 3, 6, 7, 8 give E_B = 0.5 + 0.444 + 1.0 + 0.25 + 0.333 + 0 = 2.53 against O_B = 4, and V = 1.335. The statistic is 1.472^2/1.335 = 1.62, with p = 0.20: no evidence of a difference with ten patients.

6. ★ Reading a Cox model. A Cox model of death gives a treatment coefficient -0.35 (standard error 0.12) and an age coefficient 0.40 per 10 years (standard error 0.08), with age centred at 60. The Breslow baseline cumulative hazard at 5 years is 0.30. (a) Give the treatment hazard ratio and its 95 % CI. (b) What is the hazard ratio of a 70-year-old against a 50-year-old with the same treatment? (c) Predict the 5-year survival of a treated 70-year-old and of an untreated 60-year-old. (d) Which assumption would you check first?

   Solution: (a) e^{-0.35} = 0.70, CI [0.56, 0.89]. (b) e^{2 \times 0.40} = 2.23.

   3. Definition 5.4: \exp(-0.30\, e^{-0.35 + 0.40}) = 0.730 for the treated 70-year-old, and e^{-0.30} = 0.741 for the untreated 60-year-old. (d) Proportional hazards, with Schoenfeld residuals (section 5.3).

7. ★ Discrimination, calibration and net benefit by hand. Eight patients have predicted risks 0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1 and outcomes 1, 1, 0, 1, 0, 0, 1, 0. Compute the AUC, the calibration-in-the-large, the Brier and scaled Brier scores, and the net benefit at t = 0.5.

   Solution: Of the 16 case–control pairs, the case ranks higher in 4 + 4 + 3 + 1 = 12: \mathrm{AUC} = 0.75 (Theorem 6.3). The mean prediction and the observed proportion are both 0.5.

   The Brier score is 1.6/8 = 0.20; predicting 0.5 for everyone scores 0.25, so the scaled Brier score is 0.20. At t = 0.5 four patients are treated, 3 true and 1 false positives: \mathrm{NB} = 3/8 - 1/8 = 0.25, against 0 for treating all or none.

8. ★ Standardization equals weighting on a table. On the kidney-stone data of Example 3.6, (a) compute the propensity of surgery A in each stratum, (b) compute the IPW estimate of the success rate under A for everyone, and (c) explain why it equals the standardized rate of Example 7.3.

   Solution: (a) e = 87/357 = 0.244 for small stones and 263/343 = 0.767 for large ones. (b) (81/0.244 + 192/0.767)/700 = 0.833.

   3. With a saturated propensity model, \sum_{i \in k,\, A_i = 1} Y_i / \hat e_k = n_k\, \hat p_{A,k} in each stratum k, so the IPW estimate is \sum_k (n_k/n)\, \hat p_{A,k}, the standardized rate.

9. ★★ Misclassification and regression to the mean. (a) True one-year risks are 0.10 and 0.20, and the outcome is recorded with \mathrm{Se} = 0.85, \mathrm{Sp} = 0.95 in both arms. Compute the observed risk difference and risk ratio. (b) Prove that non-differential outcome misclassification with \mathrm{Se} + \mathrm{Sp} > 1 moves the risk ratio towards 1. (c) Cholesterol has mean 5.5 mmol/L, standard deviation 1.0 and test–retest correlation 0.8. Patients are enrolled if their value is at least 7.5. What drop is expected at the next visit without treatment?

   Solution: (a) The observed risks are 0.13 and 0.21: \mathrm{RD}^\ast = -0.08 = 0.8 \times (-0.10) and \mathrm{RR}^\ast = 0.62 instead of 0.5.

   2. \mathrm{RR}^\ast = (a p_1 + c)/(a p_0 + c) with a = \mathrm{Se} + \mathrm{Sp} - 1 > 0 and c = 1 - \mathrm{Sp} \ge 0. This is a mediant of p_1/p_0 and c/c = 1, so it lies between them.

   3. For a standard normal, E[Z \mid Z \ge 2] = \varphi(2)/(1 - \Phi(2)) = 2.37. The enrolled mean is 7.87 and the expected retest mean 5.5 + 0.8 \times 2.37 = 7.40 (Proposition 3.4): a drop of 0.47 mmol/L.

10. ★★ The exponential model with censoring. (a) From Proposition 4.6, derive the maximum likelihood estimator of a constant hazard and its observed information. (b) For the ten patients of Example 4.9, give a 95 % CI for \lambda and for the median survival. (c) Prove that \hat\lambda is consistent under independent censoring when the hazard is truly constant.

    Solution: (a) \ell(\lambda) = D \log\lambda - \lambda \sum_i x_i gives \hat\lambda = D/\sum_i x_i and I(\hat\lambda) = D/\hat\lambda^2, so \mathrm{se}(\log\hat\lambda) = 1/\sqrt D.

    2. \hat\lambda = 0.064; the interval \hat\lambda\, e^{\pm 1.96/\sqrt 5} is [0.027, 0.154] per month, and the median \log 2/\lambda lies in [4.5, 26.0] months: five events say little.

    3. By the law of large numbers, \hat\lambda \to E[\delta]/E[X]. With E[X] = \int_0^\infty S_T S_C\, dt and E[\delta] = \int_0^\infty \lambda S_T S_C\, dt, the ratio is \lambda.

11. ★★ Immortal time with administrative censoring (computational). Modify Codes 3.12 and 5.7: deaths at rate \lambda = 0.1 per year, treatment start at rate \mu = 0.3, follow-up stopped at 5 years. Compute (a) the naive “ever versus never” rate ratio, (b) the person-time rate ratio, and (c) a landmark rate ratio at 1 year. Explain each result.

    Hints and expected result: with 2 \times 10^5 simulated subjects, the naive ratio is about 0.30, and the person-time and landmark ratios are both 1.00 up to Monte Carlo error (about ±0.02). The landmark analysis is unbiased here because the treatment has no effect; with a real effect it would misclassify patients who start after the landmark, which dilutes the effect.

12. ★★ Optimism with variable selection (computational). In Code 6.14, add a selection step: keep the predictors whose Wald statistic exceeds 1.96 in the full model, then refit. Compare the optimism-corrected AUC when the bootstrap repeats the selection and when it only refits the selected model, with the AUC on 20 000 new patients.

    Hints and expected result: with the seed of Code 6.14, seven predictors are selected; the apparent AUC is 0.841, the corrected AUC 0.788 when the selection is repeated and 0.817 when it is not, and the external AUC 0.735. Omitting the selection from the bootstrap hides part of the optimism (Algorithm 6.10, step 2).

13. ★★ Horvitz–Thompson, Hájek and stabilized weights. (a) Show that E[A/e(L)] = 1 under positivity. (b) Deduce that the Hájek estimator, the weighted mean of Y among the treated with weights 1/\hat e(L), is consistent for E[Y(1)] under the assumptions of Theorem 7.2. (c) Show that the stabilized weights P(A = 1)/e(L) for the treated have mean 1 among the treated, and explain why they are preferred for marginal structural models.

    Solution: (a) E[A/e(L)] = E\big[E[A \mid L]/e(L)\big] = 1. (b) The Hájek estimator is the ratio of \frac1n\sum A_iY_i/\hat e(L_i), which tends to E[Y(1)] by Theorem 7.4, to \frac1n \sum A_i/\hat e(L_i), which tends to 1 by (a).

    3. E[P(A=1)/e(L) \mid A = 1] = E[A/e(L)] = 1. Stabilized weights keep the size of the pseudo-population equal to the sample size and vary less, which typically narrows the confidence intervals of a marginal structural model that is not saturated (Hernán and Robins, 2020, chapter 12).

14. ★★ Non-compliance and clustering in a trial design. (a) In the trial of Example 2.8, only 60 % of the patients assigned to the treatment take it, and controls have no access. What ITT effect should the design assume, and how many patients per arm are needed (power 0.90)? (b) A cluster trial randomizes 30 general practices per arm, with 25 patients each, and the intracluster correlation is 0.02. What is the design effect, and the effective number of patients per arm?

    Solution: (a) By Proposition 8.3 the ITT difference is 0.6 \times (-0.10) = -0.06, from 0.30 to 0.24. Theorem 2.7 gives 10.51 \times (0.21 + 0.1824)/0.0036 = 1146 per arm, three times more than 389.

    2. 1 + 24 \times 0.02 = 1.48 (Proposition 8.7), so the 750 patients per arm are worth 507 independent patients.

15. ★★ The Aalen–Johansen estimates add up. (a) Prove that the Aalen–Johansen estimates of all causes sum to one minus the all-cause Kaplan–Meier estimate: \sum_j \hat F_j(t) = 1 - \hat S(t). (b) With constant hazards 0.05 and 0.10, compute the sum of the two “one minus Kaplan–Meier” quantities at t = 20, and comment.

    Solution: (a) At each event time t_i, with d_i = \sum_j d_{j,i} events of all causes, \hat S(t_i^-) - \hat S(t_i) = \hat S(t_i^-)\,\frac{d_i}{n_i} = \sum_j \hat S(t_i^-)\,\frac{d_{j,i}}{n_i} . Summing over t_i \le t telescopes to 1 - \hat S(t) on the left and gives \sum_j \hat F_j(t) on the right.

    2. (1 - e^{-1}) + (1 - e^{-2}) = 1.50: the naive “probabilities” of the two causes add up to more than 1, which shows that they are not real-world risks (Proposition 8.14).

16. ★★★ Non-collapsibility of the odds ratio. In a randomized trial, a prognostic factor L with P(L = 1) = 0.5 is independent of treatment. The control risk is 0.2 if L = 0 and 0.6 if L = 1, and the treatment multiplies the odds by 3 in both strata. (a) Compute the marginal risks and the marginal odds ratio. (b) Do the same for the risk difference. (c) Discuss which effect a trial should report, and what happens to a hazard ratio when a prognostic covariate is omitted from a Cox model.

    Hints and expected result: (a) The treated risks are 0.429 and 0.818, so the marginal risks are 0.623 and 0.400 and the marginal odds ratio is 2.48, below the common conditional value 3, without any confounding. (b) The marginal risk difference, 0.223, is the average of the conditional ones (0.229 and 0.218): it is collapsible.

    3. Conditional and marginal odds ratios answer different questions; adjusting for prognostic factors in a logistic model increases power, although the standard error of the log odds ratio grows (Robinson and Jewell, 1991), and it changes the estimand. The same attenuation affects hazard ratios when a prognostic covariate is omitted (Gail, Wieand and Piantadosi, 1984). The choice links to the estimand of section 8.1 and to the open questions of ch. 9.

17. ★★★ A group-sequential survival trial (computational). The trial of Example 8.12 (HR 0.75, power 0.80) adds an interim analysis after 190 of the 380 events, with O’Brien–Fleming boundaries 2.797 and 1.977. (a) Simulate the log-rank statistics at both looks, using Theorem 8.11 and the correlation of Proposition 8.4, and estimate the type I error, the power, the probability of stopping early and the expected number of events. (b) Then simulate full trials with exponential survival and uniform accrual and compare.

    Hints and expected result: in the normal approximation, the type I error is 0.050, the power 0.797 (0.801 for the fixed design), the probability of stopping at the interim 0.21 under the alternative, and the expected number of events 340 instead of 380. Full simulations should agree within Monte Carlo error.

18. ★★★ Doubly robust estimation with machine learning. Replace the parametric models of Code 7.13 by random forests or gradient boosting from scikit-learn, (a) fitted and evaluated on the same data, and (b) with two-fold cross-fitting, as in double machine learning. Over 100 replications with n = 2000, compare the bias, the standard deviation and the coverage of the Wald interval based on the influence function.

    Hints and expected result: the variance is estimated by the empirical variance of the terms in the sum of Proposition 7.7, divided by n. Without cross-fitting, overfitted outcome models leave residuals too small and the coverage falls below 95 %; with cross-fitting it should be close to 95 % (Chernozhukov et al., 2018). Clip the estimated propensity scores, for instance to [0.01, 0.99].

Projects and assessment

The official assessment (master-mva.com, “Mode de validation”) combines a final exam (60%) and a home project with its presentation (40%); the home projects are carried out in groups and presented in the final 3-hour session. The group subjects below are proposals sized for that home project, the solo subjects suit individual practice or a smaller variant, and the written exam is a proposal in the spirit of the final exam. All use public or simulated data and run on a laptop.

Solo projects

1. [Solo · ≈ 20 h] Immortal time bias, reproduced. Paper study of Suissa (2008), “Immortal time bias in pharmaco-epidemiology”. Simulate cohorts in which a drug without effect is started during follow-up, and compare the analyses described in the paper: time-fixed “ever users”, exclusion of the immortal time, landmark, and time-dependent Cox model.
   • Success criterion: the naive rate ratios match Proposition 3.10 within Monte Carlo error, and the 95 % interval of the time-dependent hazard ratio covers 1 in 93–97 % of 500 replications.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters: ch. 3, 5. Evaluation: the rubric below.
2. [Solo · ≈ 20 h] Internal validation strategies compared. Paper study of Steyerberg et al. (2001). On the Wisconsin diagnostic breast cancer data (569 patients, bundled with scikit-learn as load_breast_cancer), draw development samples of 100 patients, fit a penalized logistic model, and estimate its AUC and calibration slope by split-sample, 10-fold cross-validation and bootstrap optimism correction.
   • Success criterion: using the remaining patients as a test set, report the bias and root mean squared error of each method over 200 draws, and check whether the bootstrap has the smallest error, as the paper concludes.
   • Deliverables: code, a 4–6 page report with a decision curve, a 10-minute talk. Chapters: ch. 0, 6. Evaluation: the rubric below.
3. [Solo · ≈ 20 h] Competing risks in primary biliary cholangitis. Use the Mayo Clinic PBC trial data (Fleming and Harrington, 1991, appendix D; distributed as pbc with the R package survival, from which it exports to a CSV file), with death and liver transplantation as competing events.
   • Work: implement the Aalen–Johansen estimator and the cause-specific Cox model from Codes 5.7 and 8.15, with tied times handled, compare with one minus Kaplan–Meier, and interpret the effect of bilirubin on each cause.
   • Success criterion: your curves match an established implementation (survfit of the R package survival, or AalenJohansenFitter of lifelines, which breaks ties by a small jitter) to three decimals, and the overestimation of Proposition 8.14 is quantified.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters: ch. 4, 5, 8. Evaluation: the rubric below.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Emulating a target trial: smoking cessation and weight gain. Use the NHEFS data distributed with the book of Hernán and Robins (2020), which follow smokers of the NHANES I Epidemiologic Follow-up Study from 1971–1975 to 1982.
   • Work: write the target-trial protocol, then estimate the effect of quitting smoking on weight gain by standardization, inverse probability weighting, AIPW and propensity-score matching; check balance; compute an E-value.
   • Success criterion: reproduce the IP-weighted and standardized estimates reported in chapters 12 and 13 of the book to within 0.1 kg, then explain the differences between the four estimators.
   • Split: protocol and data preparation; outcome and propensity models; sensitivity analyses and report. Deliverables: code, a 6–8 page report, a 15-minute group presentation. Chapters: ch. 3, 7.
   • Evaluation: the rubric below, with the presentation given in the final 3-hour session.
2. [Group of 2–3 · ≈ 30 h each] Designing a survival trial when hazards are not proportional. Paper study of Royston and Parmar (2013). Design a two-arm trial with Schoenfeld’s formula, then simulate it under proportional hazards and under a delayed treatment effect, as seen with some immunotherapies.
   • Work: compare the log-rank test, a weighted log-rank test emphasizing late differences, and the difference in restricted mean survival time, with and without an O’Brien–Fleming interim analysis.
   • Success criterion: type I errors within 0.045–0.055 over 10 000 simulated trials; power curves and expected numbers of events for each method and scenario.
   • Split: simulation engine and accrual model; tests and interim analysis; summary of operating characteristics. Deliverables: code, a 6–8 page report, a 15-minute presentation. Chapters: ch. 2, 4, 5, 8.
   • Evaluation: the rubric below, with the presentation given in the final 3-hour session.
3. [Group of 2–3 · ≈ 30 h each] Developing and externally validating a prognostic model in breast cancer. Paper study of Royston and Altman (2013), who developed a Cox model on the Rotterdam breast cancer data and validated it on the German Breast Cancer Study Group trial (Schumacher et al., 1994); both data sets are distributed with the R packages survival and TH.data.
   • Work: develop a Cox model for recurrence-free survival, then validate it externally at 5 years: Harrell’s and Uno’s C, calibration plot, IPCW Brier score, and recalibration of the baseline hazard.
   • Success criterion: the full validation table of the paper’s kind for your model, with bootstrap intervals, and a comparison with a random survival forest.
   • Split: development and internal validation; external validation measures; recalibration, comparison and report. Deliverables: code, a 6–8 page report, a 15-minute presentation. Chapters: ch. 5, 6, 8.
   • Evaluation: the rubric below, with the presentation given in the final 3-hour session.

Written exam (proposal)

Format: 2 hours, one handwritten A4 sheet and a calculator allowed, four exercises.

1. Question (4 points). Define confounding and immortal-time bias. For each, give a clinical example, say in which direction it biases the estimate, and name one design and one analysis that remove it.

   Expected: Definition 3.5 with confounding by indication as example, either direction; Definition 3.9 with the “ever treated” example, bias in favour of treatment (Proposition 3.10). Remedies: randomization or stratification, regression, standardization, IPW for confounding; time-dependent covariates, landmark or target trial emulation with an aligned time zero for immortal time.

2. Question (6 points). Group A has survival times 3, 5, 7+, 9, 11+ and group B has 1, 2, 4, 6+, 8 (months). (a) Compute the Kaplan–Meier curves and the median survival in each group. (b) Compute the log-rank statistic and conclude at the 5 % level. (c) Give one reason why a hazard ratio from a Cox model could be misleading here.

   Expected: (a) A: 0.8, 0.6, 0.3 at 3, 5, 9, median 9; B: 0.8, 0.6, 0.4, 0 at 1, 2, 4, 8, median 4. (b) O_B = 4, E_B = 2.41, V = 1.42, statistic 1.77, p = 0.18: not significant. (c) Ten patients, few events, and nothing to check proportional hazards with.

3. Question (5 points). A trial compares two strategies on mortality at one year, expected to be 30 % with the standard and 20 % with the new strategy. (a) Compute the number of patients per arm for \alpha = 0.05 and power 0.90. (b) The trial is finally randomized by hospital, with 40 patients per hospital and an intracluster correlation of 0.01: correct the sample size. (c) Explain why an interim analysis at half the patients cannot use the threshold 1.96.

   Expected: (a) 389 per arm (Example 2.8). (b) Design effect 1 + 39 \times 0.01 = 1.39, so 541 per arm, about 14 hospitals. (c) Two looks at 1.96 give a type I error of about 0.083 (Proposition 8.4); use, for example, O’Brien–Fleming boundaries.

4. Question (5 points). Under consistency, conditional exchangeability given L and positivity, prove that E[Y(1)] = E[AY/e(L)] with e(L) = P(A = 1 \mid L). Then prove that A \perp L \mid e(L).

   Expected: the three-line proof of Theorem 7.4 with each assumption named where it is used; the balancing property of Theorem 7.6 through the tower property.

Grading rubric (proposal)

References

Official sources (2026–2027):

• Course page “Biostatistics” on master-mva.com: https://www.master-mva.com/cours/biostatistics/ — objectives, organisation of the sessions, assessment and themes. The official pages list no bibliography; the references below are standard works chosen for these notes.
• Université Paris-Saclay, M2 Mathématiques, Vision, Apprentissage, programme table and course detail “Biostatistics” (volume, ECTS, nature of the assessment).

Textbooks:

• S. J. Pocock, Clinical Trials: A Practical Approach, Wiley, 1983.
• L. M. Friedman, C. D. Furberg, D. L. DeMets, D. M. Reboussin, C. B. Granger, Fundamentals of Clinical Trials, 5th ed., Springer, 2015.
• S. Senn, Statistical Issues in Drug Development, 2nd ed., Wiley, 2007.
• C. Jennison, B. W. Turnbull, Group Sequential Methods with Applications to Clinical Trials, Chapman & Hall/CRC, 2000.
• K. J. Rothman, S. Greenland, T. L. Lash, Modern Epidemiology, 3rd ed., Lippincott Williams & Wilkins, 2008.
• J. D. Kalbfleisch, R. L. Prentice, The Statistical Analysis of Failure Time Data, 2nd ed., Wiley, 2002.
• J. P. Klein, M. L. Moeschberger, Survival Analysis: Techniques for Censored and Truncated Data, 2nd ed., Springer, 2003.
• T. R. Fleming, D. P. Harrington, Counting Processes and Survival Analysis, Wiley, 1991.
• P. K. Andersen, Ø. Borgan, R. D. Gill, N. Keiding, Statistical Models Based on Counting Processes, Springer, 1993.
• T. M. Therneau, P. M. Grambsch, Modeling Survival Data: Extending the Cox Model, Springer, 2000.
• E. W. Steyerberg, Clinical Prediction Models: A Practical Approach to Development, Validation, and Updating, 2nd ed., Springer, 2019.
• F. E. Harrell Jr., Regression Modeling Strategies, 2nd ed., Springer, 2015.
• M. A. Hernán, J. M. Robins, Causal Inference: What If, Chapman & Hall/CRC, 2020.
• G. W. Imbens, D. B. Rubin, Causal Inference for Statistics, Social, and Biomedical Sciences, Cambridge University Press, 2015.
• J. Pearl, Causality: Models, Reasoning, and Inference, 2nd ed., Cambridge University Press, 2009.
• A. W. van der Vaart, Asymptotic Statistics, Cambridge University Press, 1998.
• S. M. Berry, B. P. Carlin, J. J. Lee, P. Müller, Bayesian Adaptive Methods for Clinical Trials, Chapman & Hall/CRC, 2010.

Chapters 0–2 (inference, designs, trials):

• B. Woolf, “On estimating the relation between blood group and disease”, Annals of Human Genetics 19(3), 1955.
• D. G. Altman, J. M. Bland, “Absence of evidence is not evidence of absence”, BMJ 311, 1995.
• J. Cornfield, “A method of estimating comparative rates from clinical data. Applications to cancer of the lung, breast, and cervix”, Journal of the National Cancer Institute 11(6), 1951.
• K. F. Schulz, D. G. Altman, D. Moher, for the CONSORT Group, “CONSORT 2010 Statement: updated guidelines for reporting parallel group randomised trials”, BMJ 340, c332, 2010.
• E. von Elm et al., “The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement: guidelines for reporting observational studies”, The Lancet 370(9596), 2007.
• J. Neyman, “On the application of probability theory to agricultural experiments. Essay on principles. Section 9” (1923), translated by D. M. Dabrowska and T. P. Speed, Statistical Science 5(4), 1990.
• D. B. Rubin, “Estimating causal effects of treatments in randomized and nonrandomized studies”, Journal of Educational Psychology 66(5), 1974.
• D. R. Taves, “Minimization: a new method of assigning patients to treatment and control groups”, Clinical Pharmacology & Therapeutics 15(5), 1974.
• S. J. Pocock, R. Simon, “Sequential treatment assignment with balancing for prognostic factors in the controlled clinical trial”, Biometrics 31(1), 1975.
• The Cardiac Arrhythmia Suppression Trial (CAST) Investigators, “Preliminary report: effect of encainide and flecainide on mortality in a randomized trial of arrhythmia suppression after myocardial infarction”, New England Journal of Medicine 321(6), 1989.
• G. Piaggio, D. R. Elbourne, S. J. Pocock, S. J. W. Evans, D. G. Altman, for the CONSORT Group, “Reporting of noninferiority and equivalence randomized trials: extension of the CONSORT 2010 statement”, JAMA 308(24), 2012.
• S. Holm, “A simple sequentially rejective multiple test procedure”, Scandinavian Journal of Statistics 6(2), 1979.
• The Coronary Drug Project Research Group, “Influence of adherence to treatment and response of cholesterol on mortality in the Coronary Drug Project”, New England Journal of Medicine 303(18), 1980.

Chapter 3 (therapeutic evaluation and biases):

• D. Schwartz, J. Lellouch, “Explanatory and pragmatic attitudes in therapeutical trials”, Journal of Chronic Diseases 20(8), 1967.
• C. R. Charig, D. R. Webb, S. R. Payne, J. E. A. Wickham, “Comparison of treatment of renal calculi by open surgery, percutaneous nephrolithotomy, and extracorporeal shockwave lithotripsy”, British Medical Journal 292(6524), 1986.
• E. H. Simpson, “The interpretation of interaction in contingency tables”, Journal of the Royal Statistical Society, Series B 13(2), 1951.
• N. Mantel, W. Haenszel, “Statistical aspects of the analysis of data from retrospective studies of disease”, Journal of the National Cancer Institute 22(4), 1959.
• N. E. Breslow, “Odds ratio estimators when the data are sparse”, Biometrika 68(1), 1981.
• S. Suissa, “Immortal time bias in pharmaco-epidemiology”, American Journal of Epidemiology 167(4), 2008.
• L. E. Lévesque, J. A. Hanley, A. Kezouh, S. Suissa, “Problem of immortal time bias in cohort studies: example using statins for preventing progression of diabetes”, BMJ 340, b5087, 2010.
• N. Mantel, D. P. Byar, “Evaluation of response-time data involving transient states: an illustration using heart-transplant data”, Journal of the American Statistical Association 69(345), 1974.
• J. R. Anderson, K. C. Cain, R. D. Gelber, “Analysis of survival by tumor response”, Journal of Clinical Oncology 1(11), 1983.

Chapters 4–5 (censored data):

• E. L. Kaplan, P. Meier, “Nonparametric estimation from incomplete observations”, Journal of the American Statistical Association 53(282), 1958.
• M. Greenwood, “The natural duration of cancer”, Reports on Public Health and Medical Subjects 33, HMSO, 1926.
• O. Aalen, “Nonparametric inference for a family of counting processes”, Annals of Statistics 6(4), 1978.
• N. Mantel, “Evaluation of survival data and two new rank order statistics arising in its consideration”, Cancer Chemotherapy Reports 50(3), 1966.
• R. Peto, J. Peto, “Asymptotically efficient rank invariant test procedures”, Journal of the Royal Statistical Society, Series A 135(2), 1972.
• D. R. Cox, “Regression models and life-tables”, Journal of the Royal Statistical Society, Series B 34(2), 1972.
• D. R. Cox, “Partial likelihood”, Biometrika 62(2), 1975.
• N. E. Breslow, “Covariance analysis of censored survival data”, Biometrics 30(1), 1974.
• P. K. Andersen, R. D. Gill, “Cox’s regression model for counting processes: a large sample study”, Annals of Statistics 10(4), 1982.
• P. M. Grambsch, T. M. Therneau, “Proportional hazards tests and diagnostics based on weighted residuals”, Biometrika 81(3), 1994.
• M. A. Hernán, “The hazards of hazard ratios”, Epidemiology 21(1), 2010.

Chapter 6 (validation of prediction models):

• G. S. Collins, J. B. Reitsma, D. G. Altman, K. G. M. Moons, “Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD): the TRIPOD statement”, BMJ 350, g7594, 2015.
• G. S. Collins et al., “TRIPOD+AI statement: updated guidance for reporting clinical prediction models that use regression or machine learning methods”, BMJ 385, e078378, 2024.
• E. W. Steyerberg, A. J. Vickers, N. R. Cook, T. Gerds, M. Gonen, N. Obuchowski, M. J. Pencina, M. W. Kattan, “Assessing the performance of prediction models: a framework for traditional and novel measures”, Epidemiology 21(1), 2010.
• J. A. Hanley, B. J. McNeil, “The meaning and use of the area under a receiver operating characteristic (ROC) curve”, Radiology 143(1), 1982.
• F. E. Harrell, R. M. Califf, D. B. Pryor, K. L. Lee, R. A. Rosati, “Evaluating the yield of medical tests”, JAMA 247(18), 1982.
• H. Uno, T. Cai, M. J. Pencina, R. B. D’Agostino, L. J. Wei, “On the C-statistics for evaluating overall adequacy of risk prediction procedures with censored survival data”, Statistics in Medicine 30(10), 2011.
• B. Van Calster, D. J. McLernon, M. van Smeden, L. Wynants, E. W. Steyerberg, “Calibration: the Achilles heel of predictive analytics”, BMC Medicine 17, 230, 2019.
• G. W. Brier, “Verification of forecasts expressed in terms of probability”, Monthly Weather Review 78(1), 1950.
• B. Efron, “Estimating the error rate of a prediction rule: improvement on cross-validation”, Journal of the American Statistical Association 78(382), 1983.
• F. E. Harrell, K. L. Lee, D. B. Mark, “Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors”, Statistics in Medicine 15(4), 1996.
• E. W. Steyerberg, F. E. Harrell, G. J. J. M. Borsboom, M. J. C. Eijkemans, Y. Vergouwe, J. D. F. Habbema, “Internal validation of predictive models: efficiency of some procedures for logistic regression analysis”, Journal of Clinical Epidemiology 54(8), 2001.
• P. Peduzzi, J. Concato, E. Kemper, T. R. Holford, A. R. Feinstein, “A simulation study of the number of events per variable in logistic regression analysis”, Journal of Clinical Epidemiology 49(12), 1996.
• R. D. Riley et al., “Calculating the sample size required for developing a clinical prediction model”, BMJ 368, m441, 2020.
• J. C. van Houwelingen, S. le Cessie, “Predictive value of statistical models”, Statistics in Medicine 9(11), 1990.
• A. J. Vickers, E. B. Elkin, “Decision curve analysis: a novel method for evaluating prediction models”, Medical Decision Making 26(6), 2006.

Chapter 7 (causal inference):

• J. M. Robins, “A new approach to causal inference in mortality studies with a sustained exposure period — application to control of the healthy worker survivor effect”, Mathematical Modelling 7(9–12), 1986.
• D. G. Horvitz, D. J. Thompson, “A generalization of sampling without replacement from a finite universe”, Journal of the American Statistical Association 47(260), 1952.
• P. R. Rosenbaum, D. B. Rubin, “The central role of the propensity score in observational studies for causal effects”, Biometrika 70(1), 1983.
• P. C. Austin, “An introduction to propensity score methods for reducing the effects of confounding in observational studies”, Multivariate Behavioral Research 46(3), 2011.
• J. M. Robins, A. Rotnitzky, L. P. Zhao, “Estimation of regression coefficients when some regressors are not always observed”, Journal of the American Statistical Association 89(427), 1994.
• J. Pearl, “Causal diagrams for empirical research”, Biometrika 82(4), 1995.
• M. A. Hernán, S. Hernández-Díaz, J. M. Robins, “A structural approach to selection bias”, Epidemiology 15(5), 2004.
• M. A. Hernán, J. M. Robins, “Using big data to emulate a target trial when a randomized trial is not available”, American Journal of Epidemiology 183(8), 2016.
• M. A. Hernán, B. C. Sauer, S. Hernández-Díaz, R. Platt, I. Shrier, “Specifying a target trial prevents immortal time bias and other self-inflicted injuries in observational analyses”, Journal of Clinical Epidemiology 79, 2016.
• T. J. VanderWeele, P. Ding, “Sensitivity analysis in observational research: introducing the E-value”, Annals of Internal Medicine 167(4), 2017.

Chapter 8 (complements):

• International Council for Harmonisation, ICH E9(R1) Addendum on Estimands and Sensitivity Analysis in Clinical Trials to the Guideline on Statistical Principles for Clinical Trials, 2019.
• J. D. Angrist, G. W. Imbens, D. B. Rubin, “Identification of causal effects using instrumental variables”, Journal of the American Statistical Association 91(434), 1996.
• S. J. Pocock, “Group sequential methods in the design and analysis of clinical trials”, Biometrika 64(2), 1977.
• P. C. O’Brien, T. R. Fleming, “A multiple testing procedure for clinical trials”, Biometrics 35(3), 1979.
• K. K. G. Lan, D. L. DeMets, “Discrete sequential boundaries for clinical trials”, Biometrika 70(3), 1983.
• ISIS-2 (Second International Study of Infarct Survival) Collaborative Group, “Randomised trial of intravenous streptokinase, oral aspirin, both, or neither among 17,187 cases of suspected acute myocardial infarction: ISIS-2”, The Lancet 332(8607), 1988.
• P. Royston, M. K. B. Parmar, “Restricted mean survival time: an alternative to the hazard ratio for the design and analysis of randomized trials with a time-to-event outcome”, BMC Medical Research Methodology 13, 152, 2013.
• D. A. Schoenfeld, “Sample-size formula for the proportional-hazards regression model”, Biometrics 39(2), 1983.
• O. O. Aalen, S. Johansen, “An empirical transition matrix for non-homogeneous Markov chains based on censored observations”, Scandinavian Journal of Statistics 5(3), 1978.
• J. P. Fine, R. J. Gray, “A proportional hazards model for the subdistribution of a competing risk”, Journal of the American Statistical Association 94(446), 1999.
• H. Putter, M. Fiocco, R. B. Geskus, “Tutorial in biostatistics: competing risks and multi-state models”, Statistics in Medicine 26(11), 2007.
• P. C. Austin, D. S. Lee, J. P. Fine, “Introduction to the analysis of survival data in the presence of competing risks”, Circulation 133(6), 2016.
• P. J. Heagerty, T. Lumley, M. S. Pepe, “Time-dependent ROC curves for censored survival data and a diagnostic marker”, Biometrics 56(2), 2000.
• P. Blanche, J.-F. Dartigues, H. Jacqmin-Gadda, “Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks”, Statistics in Medicine 32(30), 2013.
• E. Graf, C. Schmoor, W. Sauerbrei, M. Schumacher, “Assessment and comparison of prognostic classification schemes for survival data”, Statistics in Medicine 18(17–18), 1999.
• P. Royston, D. G. Altman, “External validation of a Cox prognostic model: principles and methods”, BMC Medical Research Methodology 13, 33, 2013.

Chapter 9 (research directions):

• C. E. Frangakis, D. B. Rubin, “Principal stratification in causal inference”, Biometrics 58(1), 2002.
• H. Ishwaran, U. B. Kogalur, E. H. Blackstone, M. S. Lauer, “Random survival forests”, Annals of Applied Statistics 2(3), 2008.
• J. L. Katzman, U. Shaham, A. Cloninger, J. Bates, T. Jiang, Y. Kluger, “DeepSurv: personalized treatment recommender system using a Cox proportional hazards deep neural network”, BMC Medical Research Methodology 18, 24, 2018.
• E. Christodoulou, J. Ma, G. S. Collins, E. W. Steyerberg, J. Y. Verbakel, B. Van Calster, “A systematic review shows no performance benefit of machine learning over logistic regression for clinical prediction models”, Journal of Clinical Epidemiology 110, 2019.
• E. J. Candès, L. Lei, Z. Ren, “Conformalized survival analysis”, Journal of the Royal Statistical Society, Series B 85(1), 2023.
• S. Wager, S. Athey, “Estimation and inference of heterogeneous treatment effects using random forests”, Journal of the American Statistical Association 113(523), 2018.
• S. R. Künzel, J. S. Sekhon, P. J. Bickel, B. Yu, “Metalearners for estimating heterogeneous treatment effects using machine learning”, Proceedings of the National Academy of Sciences 116(10), 2019.
• V. Chernozhukov, D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, J. Robins, “Double/debiased machine learning for treatment and structural parameters”, The Econometrics Journal 21(1), 2018.
• D. M. Kent et al., “The Predictive Approaches to Treatment effect Heterogeneity (PATH) statement”, Annals of Internal Medicine 172(1), 2020.
• J. M. Franklin et al., “Emulating randomized clinical trials with nonrandomized real-world evidence studies: first results from the RCT DUPLICATE initiative”, Circulation 143(10), 2021.
• I. J. Dahabreh, S. E. Robertson, J. A. Steingrimsson, E. A. Stuart, M. A. Hernán, “Extending inferences from a randomized trial to a new target population”, Statistics in Medicine 39(14), 2020.
• L. Wynants et al., “Prediction models for diagnosis and prognosis of covid-19: systematic review and critical appraisal”, BMJ 369, m1328, 2020.
• Z. Obermeyer, B. Powers, C. Vogeli, S. Mullainathan, “Dissecting racial bias in an algorithm used to manage the health of populations”, Science 366(6464), 2019.
• J. Woodcock, L. M. LaVange, “Master protocols to study multiple therapies, multiple diseases, or both”, New England Journal of Medicine 377(1), 2017.
• The RECOVERY Collaborative Group, “Dexamethasone in hospitalized patients with Covid-19”, New England Journal of Medicine 384(8), 2021.

Projects and exercises:

• M. H. Gail, S. Wieand, S. Piantadosi, “Biased estimates of treatment effect in randomized experiments with nonlinear regressions and omitted covariates”, Biometrika 71(3), 1984.
• L. D. Robinson, N. P. Jewell, “Some surprising results about covariate adjustment in logistic regression models”, International Statistical Review 59(2), 1991.
• W. N. Street, W. H. Wolberg, O. L. Mangasarian, “Nuclear feature extraction for breast tumor diagnosis”, IS&T/SPIE International Symposium on Electronic Imaging: Science and Technology, volume 1905, 1993.
• M. Schumacher et al., “Randomized 2 × 2 trial evaluating hormonal treatment and the duration of chemotherapy in node-positive breast cancer patients”, Journal of Clinical Oncology 12(10), 1994.