Module 4 · Cox model diagnostics

Lecture 4.1 onward. Reading: Lecture Notes Part 1. Code: R script — Lecture Video 14, R script — Lecture Video 16, R script — Lecture Video 20, R script — Lecture Video 24, R script — Lecture Video 26, R script — Lecture Video 29.

← Back to notes index

Picks up where Notes 3 left off — the Cox model $\lambda(t \mid \mathbf{Z}) = \lambda_0(t)\, e^{\beta^T \mathbf{Z}}$ and the per-event-time partial likelihood contribution $L_j(\beta) = \exp(\beta^T \mathbf{Z}j) / \sum{l \in \mathcal{R}_j} \exp(\beta^T \mathbf{Z}_l)$.

Lecture 4.1

Martingale residuals

In ordinary linear regression a residual is

\[Y_i - \hat{Y}_i\]

— observed minus fitted. Survival data has no single “observed $Y$” to subtract from, because of censoring; the martingale residual is the survival-analysis analogue, built on the same observed-minus-expected idea:

the 0/1 event indicator vs. the expected value of that indicator.

For subject $i$:

\[m_i = \delta_i - \hat{\Lambda}_0(x_i)\, \exp(\hat{\beta}^T \mathbf{Z}_i)\]
Piece Role
$\delta_i$ observed — did subject $i$ actually have the event (1) or were they censored (0)?
$\hat{\Lambda}_0(x_i)\exp(\hat{\beta}^T \mathbf{Z}_i)$ expected — the fitted cumulative hazard for subject $i$ over their observed follow-up $x_i$, i.e. the model’s expected number of events for them

The subtracted term is just the Cox model’s cumulative hazard evaluated at that subject’s own censored time $x_i$ and own covariates: the baseline cumulative hazard $\hat{\Lambda}_0$ scaled by their relative-risk multiplier $\exp(\hat{\beta}^T \mathbf{Z}_i)$. Both pieces are estimated from the fitted model.

Properties

  • Range: $-\infty$ to $1$. The maximum is $1$ — a subject who has the event ($\delta_i = 1$) but whose expected cumulative hazard is essentially $0$, i.e. an event the model did not see coming at all. There is no lower bound: a subject can accumulate arbitrarily large expected hazard (long follow-up, high-risk covariates) and still never have the event, driving $m_i$ far negative.
  • Expected value is $0$. Like ordinary residuals, they are centred at zero under a correctly specified model.

The asymmetry of the range is worth holding onto — martingale residuals are not symmetric like linear-regression residuals, so a plot of them is skewed by construction, not necessarily by model misfit.

Reading the sign

Sign of $m_i$ Meaning
positive the patient died sooner than expected — the event arrived before the model had accumulated much hazard for them
negative the patient survived longer than expected — the model piled up expected hazard that never cashed out in an event

So a large positive residual flags a subject the model treated as low-risk who died anyway; a large negative one flags a subject the model treated as high-risk who kept going.

Deviance residuals

\[d_i = \operatorname{sign}(m_i) \Big\{ -2 \big[\, m_i + \delta_i \log(\delta_i - m_i) \,\big] \Big\}^{1/2}\]

A rescaling of the martingale residual that keeps its sign — $\operatorname{sign}(m_i)$ is pulled out front, so the interpretation above (positive = died sooner than expected, negative = survived longer) carries over unchanged — while the ${-2[\cdot]}^{1/2}$ transform pulls in the long negative tail.

Note what the log term is: since $m_i = \delta_i - \hat{\Lambda}_0(x_i)\exp(\hat{\beta}^T \mathbf{Z}_i)$, the quantity inside it is exactly the expected cumulative hazard,

\[\delta_i - m_i = \hat{\Lambda}_0(x_i)\, \exp(\hat{\beta}^T \mathbf{Z}_i)\]

and for a censored subject ($\delta_i = 0$) the whole log term drops out.

Properties:

  • A symmetric version of the martingale residuals — this is the point of the transform: it fixes the awkward $(-\infty, 1]$ asymmetry.
  • If the fitted model is correct, the deviance residuals are symmetrically distributed with expected value $0$ — which makes them usable the way ordinary regression residuals are: plot them, and look for asymmetry, structure, or outliers as evidence against the model.

The proportional hazards assumption

The Cox model’s central assumption needs checking, and the check comes from pushing the assumption through the $S \leftrightarrow \Lambda$ identities until it lands on something plottable.

Take the simplest case: a single binary covariate $Z \in {0, 1}$, so the hazard ratio for $Z = 1$ vs. $Z = 0$ is $e^{\beta}$. Write $\lambda_1, \lambda_0$ for the two groups’ hazards. The assumption is exactly

\[\frac{\lambda_1(t)}{\lambda_0(t)} = e^{\beta}\]

— constant in $t$. Now carry it along the chain, using $S(t) = e^{-\Lambda(t)}$:

\[\begin{aligned} \frac{\lambda_1(t)}{\lambda_0(t)} = e^{\beta} &\;\Longrightarrow\; \lambda_1(t) = e^{\beta}\, \lambda_0(t) \\[6pt] &\;\Longrightarrow\; \Lambda_1(t) = e^{\beta}\, \Lambda_0(t) \\[6pt] &\;\Longrightarrow\; -\log\{S_1(t)\} = -\log\{S_0(t)\}\, e^{\beta} \\[6pt] &\;\Longrightarrow\; \log\big[-\log\{S_1(t)\}\big] = \log\big[-\log\{S_0(t)\}\big] + \beta \end{aligned}\]

Step by step: multiplying the hazards by a constant multiplies their integrals by the same constant, so the cumulative hazards are in the same proportion $e^{\beta}$. Then $\Lambda(t) = -\log S(t)$ rewrites that in terms of the survivor functions. Finally take logs of both sides — the multiplicative factor $e^{\beta}$ becomes an additive $\beta$.

The function applied to the survivor function on that last line,

\[S \;\longmapsto\; \log\big[-\log S\big]\]

is the complementary log-log transform. The circled ratio at the top is the assumption; everything below it is that assumption re-expressed, and the transform is what turns it into a statement you can check by eye.

The complementary log-log plot

That last line is the diagnostic. Plot $\log[-\log \hat{S}_j(t)]$ against $t$ (or $\log t$) for each group, using the Kaplan–Meier estimates:

  • Under proportional hazards the two curves are parallel, separated by a constant vertical gap of exactly $\beta$ — and that gap is the log hazard ratio.
  • Curves that converge, diverge, or cross are evidence against proportional hazards: the gap is $\beta(t)$, not a single $\beta$.

What to look at: is the gap constant or not? Plotting $\log[-\log \hat{S}(t)]$ against $t$ for both groups, the entire diagnostic reduces to one question — does the vertical distance between the two curves stay the same as $t$ moves?

β β Z = 1 Z = 0 t proportional hazards holds log[−log Ŝ(t)]
small gap early gap grows Z = 1 Z = 0 t proportional hazards fails log[−log Ŝ(t)]

Constant gap (top): the vertical distance is the same $\beta$ whether measured early or late. The curves track each other, and that one number summarizes the covariate’s effect for all time — proportional hazards holds.

Non-constant gap (bottom): the curves fan apart — small gap early, large gap late. The separation is $\beta(t)$, a moving target, so no single hazard ratio describes the covariate and the assumption fails.

Note it is only the gap that matters, not the shape of either curve: both curves may bend however they like, as long as they bend together.

Scaled Schoenfeld residuals

A second diagnostic, and the one that reads most directly. Schoenfeld residuals are defined per covariate at each event time; scaling them turns the plot into an estimate of the coefficient as a function of time,

\[\hat{\beta}(t)\]

so the plot’s vertical axis carries the same units as the fitted $\hat{\beta}$ itself. Under proportional hazards $\beta$ does not depend on $t$, so the question is again one of shape: flat, or trending?

0.2 t flat — proportional hazards holds β̂(t)
0.2 t trending — proportional hazards fails β̂(t)
  • Flat (top) — $\hat{\beta}(t)$ hovers around a single level, here $0.2$: the effect is the same early and late, which is precisely the assumption. The level it sits at is the $\hat{\beta}$ the Cox fit reports.
  • Trending (bottom) — $\hat{\beta}(t)$ starts near $0.2$ and climbs: the covariate’s effect strengthens over time. One reported number is then a time-average of something that is genuinely changing, and the proportional hazards assumption is violated.

Same logic as the complementary log-log plot, read off a different axis: there the question was whether the gap between two curves is constant in $t$; here it is whether one curve is constant in $t$.

When proportional hazards fails

The diagnostics only tell you the assumption is broken; three ways to respond:

  • Piecewise constant baseline hazard function — cut the time axis into intervals and let the hazard be constant within each one, fitting a separate effect per interval. The covariate’s effect is then allowed to differ early vs. late instead of being forced into one number.
  • Stratify the baseline hazard — let the offending variable define strata, each with its own unspecified $\lambda_{0j}(t)$, rather than entering the model as a covariate. Nothing is assumed about how the strata relate to each other, so the assumption can no longer be violated by that variable — at the cost of no longer getting a hazard ratio for it.
  • Interaction with time — put the variable’s effect explicitly in terms of $t$, e.g. a $Z \times t$ (or $Z \times \log t$) term, so the coefficient becomes $\beta(t)$ by construction. This models the trend the Schoenfeld plot revealed rather than working around it.
  • Nonparametric kernel smoothing — e.g. a smoothed covariate: rather than forcing a functional form, let the data determine the shape, via splines and in particular cubic splines. The covariate (or its effect) enters through a flexible smooth term instead of a single linear coefficient.

In R

Everything above, run on real data in R script — Lecture Video 14 (survival, plus pharmacoSmoking from asaur and BrainCancer from ISLR2).

Residuals

my.model = coxph(Surv(ttr, relapse) ~ age + grp + longestNoSmoke,
                 data = pharmacoSmoking)
mm = residuals(my.model)                      # martingale (the default)
dd = residuals(my.model, type = "deviance")   # deviance

The residuals are then plotted against each covariate — age, longestNoSmoke, grp — with a lowess curve and 95% band drawn on top by the script’s smoothSEcurve() helper (a loess fit, its fit line, and $\pm\, t_{0.975, df} \times \text{se}$ dashed bands).

This is the second job martingale residuals do, alongside flagging odd subjects: checking the functional form of a covariate. Plotted against a covariate the smooth should be flat around $0$; systematic curvature says the covariate enters the model wrongly. That is exactly what the script demonstrates on longestNoSmoke — a curved smooth, then a refit with the covariate log-transformed:

my.model = coxph(Surv(ttr, relapse) ~ age + grp + log(longestNoSmoke + 1),
                 data = pharmacoSmoking)

(+1 because the variable can be $0$.) One thing to watch in the script as given: the refit block plots against log(longestNoSmoke+1) but then calls smoothSEcurve with log(priorAttempts+1) — the smooth there is drawn over a different variable than the points.

Complementary log-log plot

Built by hand from two Kaplan–Meier fits, one per group, rather than from a model:

surv.p = survfit(Surv(ttr, relapse) ~ 1, subset = (grp == "patchOnly"),
                 data = pharmacoSmoking)
surv.c = survfit(Surv(ttr, relapse) ~ 1, subset = (grp == "combination"),
                 data = pharmacoSmoking)

cll.surv.p = log(-log(surv.p$surv))           # the transform, literally
cll.surv.c = log(-log(surv.c$surv))

plot(log(time.p), cll.surv.p, typ = "l", lwd = 3)
points(log(time.c), cll.surv.c, typ = "l", lwd = 3, col = 2)

Note the horizontal axis is $\log t$, not $t$ — the usual convention, and it does not affect the question being asked, since a constant vertical gap stays constant under any relabelling of the time axis.

Schoenfeld residuals and stratification

my.model = coxph(Surv(time, status) ~ sex + diagnosis + loc + ki + gtv + stereo,
                 data = BrainCancer)
model.test = cox.zph(my.model)

plot(model.test, var = "gtv");  abline(h = 0, col = 2)
plot(model.test, var = "sex");  abline(h = 0, col = 2)
model.test

cox.zph() is the scaled-Schoenfeld diagnostic: plot() draws $\hat{\beta}(t)$ with a confidence band, one panel per var, and printing the object gives a formal test per covariate plus a global test — the null being that $\beta$ does not depend on $t$. The abline(h = 0, col = 2) reference line makes the flat-vs-trending read easy.

And the stratification remedy, letting the baseline hazard differ by sex instead of estimating a hazard ratio for it:

my.model = coxph(Surv(time, status) ~ diagnosis + loc + ki + gtv + stereo +
                   strata(sex), data = BrainCancer)

sex no longer appears among the coefficients in summary() — that is the trade noted above.

Modeling decisions

Reading: Lecture Notes Part 1, “Modeling Decisions” (p. 33). Code: R script — Lecture Video 16.

Choosing among candidate Cox models uses the same toolkit as linear regression — stepwise or regularized variable selection, trees and boosting, and the model fit criteria below. Which criterion applies depends on one question: are the two models nested or not?

Nested models: the likelihood ratio test

Let Model A be the larger model with $a$ degrees of freedom and (partial) log-likelihood $l_a$, and Model B the smaller model nested inside it, with $b$ degrees of freedom and log-likelihood $l_b$. Under $H_0$ that the extra terms in A are all zero,

\[2\,(l_a - l_b) \;\sim\; \chi^2_{a-b}\]

so compare $2(l_a - l_b)$ to a $\chi^2$ with $a - b$ degrees of freedom to get a $p$-value. If $p < 0.05$, Model A (the larger one) is the superior model — the extra terms earn their keep.

“Nested” means B is A with some coefficients fixed at $0$: dropping a covariate, or dropping every dummy of a factor at once. The degrees of freedom counted are coefficients, so a three-level factor contributes $2$.

Non-nested models: AIC and BIC

When neither model is a special case of the other the likelihood ratio test does not apply, and an information criterion is used instead. For Model A,

\[\mathrm{AIC} = -2\,l_a + 2a\]

Two competing pieces:

Term Property of the model Direction
$-2 l_a$ goodness of fit — smaller for models that fit the data well rewards fit
$2a$ complexity — the number of parameters, entering as a penalty punishes size

A “good” model fits well (small $-2l_a$) with few parameters (small $2a$), so smaller AIC is better. Unlike the likelihood ratio test there is no $p$-value — it is a ranking, not a hypothesis test.

The Bayesian Information Criterion swaps the penalty:

\[\mathrm{BIC} = -2\,l_a + a\,\log(n)\]

The only difference is that each parameter costs $\log(n)$ rather than $2$. Since $\log(n) > 2$ once $n > 7$, BIC penalizes parameters more heavily than AIC, and so model selection by BIC tends to pick models with fewer parameters than AIC does.

In R

From R script — Lecture Video 16, on pharmacoSmoking again.

Likelihood ratio test, by hand and then via anova():

model.a = coxph(Surv(ttr, relapse) ~ age + grp + employment, data = pharmacoSmoking)
model.b = coxph(Surv(ttr, relapse) ~ age + grp,              data = pharmacoSmoking)

logLik(model.a)   # -375.1393 (df=4)
logLik(model.b)   # -379.2411 (df=2)
lrt = 2*(-375.1393 + 379.2411)          # 2 (l_a - l_b) = 8.2037
pchisq(lrt, df = 2, lower = F)          # 0.0165

anova(model.a, model.b)                 # same thing

model.b is model.a without employment — nested. employment has three levels, hence two dummies (employmentother, employmentpt) and df = 4 - 2 = 2 in the $\chi^2$. The lower = F is what turns the CDF into an upper-tail $p$-value. anova() on two coxph fits prints the same log-likelihoods, $\chi^2$, df and $p$ as an “Analysis of Deviance Table”:

   loglik  Chisq Df Pr(>|Chi|)
1 -375.14
2 -379.24 8.2037  2    0.01654 *

$p = 0.017 < 0.05$, so the larger model with employment is preferred.

AIC and BIC, against a third model that is not nested in model.a:

model.c = coxph(Surv(ttr, relapse) ~ age + grp + gender + longestNoSmoke + race,
                data = pharmacoSmoking)

AIC(model.a); AIC(model.c)   # 758.28   768.49
BIC(model.a); BIC(model.c)   # 768.23   785.91

model.c swaps employment for gender, longestNoSmoke and race — neither model contains the other, so no likelihood ratio test. Both criteria are smaller for model.a, so it wins either way; note the BIC gap is wider than the AIC gap, as expected from the heavier per-parameter penalty on the larger model.c.

Sanity check on the arithmetic: AIC(model.a) $= -2(-375.1393) + 2 \cdot 4 = 758.28$.

Prediction: evaluating predictions

Reading: Lecture Notes Part 1, §1.8.1–1.8.3 (pp. 34–43). Code: R script — Lecture Video 20.

Estimation finds the optimal value of a parameter from data already collected. Prediction uses historical data to compute the random value of future, unseen data. Everything before this point was estimation (plus computing individual predicted survival); this section is about evaluating individual-level predictions.

Loss functions and the 2×2 table

For a general outcome $Y$ with prediction $\hat{Y}$, accuracy is measured by a loss:

\[\frac{1}{n}\sum_{i=1}^{n}(Y_i - \hat{Y}_i)^2 \quad\text{(mean squared error, } L_2\text{)} \qquad\qquad \frac{1}{n}\sum_{i=1}^{n}|Y_i - \hat{Y}_i| \quad\text{(mean absolute error, } L_1\text{)}\]

Smaller is better for both. For a binary outcome with a binary prediction — the diagnostic-test setting, $D^{\pm}$ for true disease status and $T^{\pm}$ for the test result — $L_1$ and $L_2$ coincide and both equal the overall misclassification rate. But the two kinds of error usually carry different costs (a false-positive COVID test means lost wages; a false negative means spreading it), so they are reported separately:

  $D^{+}$ $D^{-}$  
$T^{+}$ $a$ (true positive) $b$ (false positive) $a+b$
$T^{-}$ $c$ (false negative) $d$ (true negative) $c+d$
  $a+c$ $b+d$  
Quantity Definition Estimate
Sensitivity $P(T^{+} \mid D^{+})$ $a/(a+c)$
Specificity $P(T^{-} \mid D^{-})$ $d/(b+d)$
PPV (positive predictive value) $P(D^{+} \mid T^{+})$ $a/(a+b)$
NPV (negative predictive value) $P(D^{-} \mid T^{-})$ $d/(c+d)$
FPR (false positive rate) $P(T^{+} \mid D^{-})$ $b/(b+d)$
FNR (false negative rate) $P(T^{-} \mid D^{+})$ $c/(a+c)$

Sensitivity and specificity are characteristics of the test; sensitivity $+$ FNR $= 1$ and specificity $+$ FPR $= 1$. PPV and NPV additionally depend on the prevalence of disease in the population, so they do not transfer between populations the way sensitivity and specificity do.

Thresholds

Real tests produce a measurement $R$ (antibody count, a cognitive score, a band ratio) and a threshold $\theta$ turns it into $T^{\pm}$. In the notes’ cognitive-functioning example a low score flags impairment, so $R \le \theta \Rightarrow T^{+}$ and $R > \theta \Rightarrow T^{-}$, giving

\[\text{Sens}(\theta) = P(R \le \theta \mid D^{+}), \qquad \text{Spec}(\theta) = P(R > \theta \mid D^{-})\]

Every $\theta$ yields its own 2×2 table and its own (sensitivity, specificity) pair. At the extremes, $\theta = 0$ calls nobody positive (sensitivity $0$, specificity $1$) and $\theta = 100$ calls everybody positive (sensitivity $1$, specificity $0$). The class exercise (18 subjects, 10 $D^{+}$ and 8 $D^{-}$) walks through $\theta = 10$ and $\theta = 50$ by hand.

The ROC curve and AUC

It is rare to get high sensitivity and high specificity; the receiver operating characteristic (ROC) curve shows the trade-off by plotting 1 − specificity ($x$) against sensitivity ($y$) for every candidate $\theta$. The diagonal through the origin with slope 1 is a useless test — equivalent to a coin flip. The point closest to the upper-left corner is sometimes taken as the “optimal” $\theta$, though the real optimum depends on the costs of the two errors.

The area under the curve (AUC) collapses this to a single number:

  • an overall summary of accuracy across all $\theta$ — the degree of separability between $D^{+}$ and $D^{-}$;
  • $1$ is perfect; the diagonal has AUC $0.5$ (useless); higher is better;
  • it is the same thing as the C-statistic (the Concordance line in summary(coxph)).

From a diagnostic test to a model

With a logistic regression (or any model) giving $p_i = \hat{P}(D^{+} \mid \mathbf{Z}_i)$, nothing changes except that the thresholds are probability cutoffs $\theta \in [0,1]$: $p_i \ge \theta \Rightarrow T^{+}$. Note the inequality has flipped relative to the cognitive example — here higher $p$ means more likely $D^{+}$ — which is exactly the direction issue in the R code below. Formally,

\[\text{Spec}(\theta) = P(p < \theta \mid D^{-}), \qquad \text{Sens}(\theta) = P(p \ge \theta \mid D^{+}), \qquad \text{ROC}(u) = \text{Sens}\{\text{Spec}^{-1}(1-u)\}, \qquad \text{AUC} = \int \text{ROC}(u)\,du\]

Estimation. Let $X_i$ ($i = 1,\dots,m$) be the values of $p$ among the $D^{+}$ and $Y_j$ ($j = 1,\dots,n$) among the $D^{-}$, $N = m + n$. Then

\[\widehat{\text{sens}}(\theta) = \frac{1}{m}\sum_{i=1}^{m} I(X_i \ge \theta), \qquad \widehat{\text{spec}}(\theta) = \frac{1}{n}\sum_{j=1}^{n} I(Y_j < \theta)\]

and the AUC has the interpretation

\[\text{AUC} = P(X > Y)\]

— the probability that a randomly chosen $D^{+}$ subject is ranked as more likely diseased than a randomly chosen $D^{-}$ subject. Its nonparametric estimate is the Mann–Whitney $U$-statistic over all $m \times n$ pairs,

\[\widehat{\text{AUC}} = \frac{1}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n} \psi(X_i, Y_j), \qquad \psi(X,Y) = \begin{cases} 1 & X > Y \\ \tfrac{1}{2} & X = Y \\ 0 & X < Y \end{cases}\]

(the $\tfrac12$ handles ties). Parametric alternatives exist — the binormal approach assumes both $X$ and $Y$ are normal — but the course uses the empirical estimate.

Brier score. With $A_i \in {0,1}$ the true status and $p_i$ the predicted probability,

\[\text{BS} = E\{(p_i - A_i)^2\}, \qquad \widehat{\text{BS}} = \frac{1}{N}\sum_{i=1}^{N}(p_i - A_i)^2\]

— the $L_2$ loss again, now between a probability and a 0/1 outcome. Smaller is better. It becomes the tool of choice when a held-out set is too small for an AUC (see LOOCV below).

Net reclassification improvement (NRI) is an alternative for incremental value when classifying into more than two categories (Pencina et al. 2008); it can be used for two categories but in practice the AUC dominates.

Incremental value

The incremental value of a predictor $Z_b$ is the AUC gain from including it:

\[IV(Z_b) = \text{AUC}\{p(Z_b)\} - \text{AUC}\{p(-Z_b)\}\]

where $p(Z_b)$ is the prediction from a model using all predictors and $p(-Z_b)$ from the same model without $Z_b$. Why care: $Z_b$ is typically costly or invasive to obtain — a genetic score, stool DNA test, breast biopsy, wearables data — versus cheap predictors like age, gender, race/ethnicity, employment, routine labs. If it adds nothing on top of what you already have, it is not worth measuring.

Two things to hold onto:

  • A predictor can be significantly associated with the outcome in a logistic regression and still have little incremental value (Pepe et al. 2004). A small $p$-value for the coefficient is not the same claim as an improvement in the ROC curve.
  • If you are arguing a new predictor is useful, you must show incremental value.

In R

From R script — Lecture Video 20. Two packages: tdROC (for calc.AUC() when you already have sensitivities and specificities) and pROC (for everything from raw data).

The class exercise ROC curve, from the four hand-computed (sens, spec) pairs:

sens = c(1, 0.9, 0.7, 0)
spec = c(0, 0.375, 0.75, 1)
plot(1-spec, sens, typ = "o", pch = 20, lwd = 2, main = "ROC curve",
     xlab = "1-Specificity", ylab = "Sensitivity")
polygon(c(1-spec, 1), c(sens, 0), col = rgb(0.51, 0.44, 1, alpha = 0.2), border = 0)
abline(a = 0, b = 1, lwd = 2, lty = 2)             # the coin-flip line

library(tdROC)
calc.AUC(sens = sens, spec = spec)                 # 0.74375

AUC from raw data — the same 18 subjects, scores in pred:

library(pROC)
outcome = c(rep(1,10), rep(0,8))
pred = c(2,5,7,9,10,7,4,15,35,51,  9,7,12,26,49,87,99,60)

roc.output = roc(response = outcome, predictor = pred)
# Setting levels: control = 0, case = 1
# Setting direction: controls > cases

Read those two messages. direction: controls > cases means roc() has decided that higher values of the predictor indicate controls — correct here, since a low cognitive score means impairment, but pROC picks the direction from a median comparison and does not always default to >. Always set direction explicitly: direction = "<" when higher values mean cases (the usual situation for a fitted probability).

cbind(roc.output$thresholds, roc.output$sensitivities, roc.output$specificities)
roc.output$auc                                     # 0.79375
calc.AUC(sens = roc.output$sensitivities, spec = roc.output$specificities)   # matches
plot(roc.output)

roc(cases = pred[1:10], controls = pred[11:18])$auc          # same, given as two groups
roc(response = outcome, predictor = pred,
    smooth = TRUE, smooth.method = "binormal")               # parametric smoothed curve

The thresholds table is the machine version of the class exercise — every midpoint between adjacent scores is a $\theta$, and the rows for $\theta = 11$ and $\theta = 50$ reproduce the hand calculations (sens $0.7$/spec $0.75$ and sens $0.9$/spec $0.375$). The AUC of $0.794$ is higher than the class exercise’s $0.744$ only because the curve now has all 16 thresholds rather than four points joined by straight lines.

From a logistic model. The general recipe — fit, predict(type = "response"), then roc() on the fitted probabilities with direction = "<":

model = glm(outcome ~ z1 + z2, family = "binomial", data = mydata)
preds = predict(model, type = "response")
roc_object = roc(mydata$outcome, preds, direction = "<")
roc_object$auc

Real data: aSAH (ships with pROC) — 113 patients with aneurysmal subarachnoid haemorrhage (Turck et al. 2010); predictors age, gender and two biomarkers s100b and ndka; outcome Good vs. Poor.

aSAH$outcome.binary = 1*(aSAH$outcome == "Good")       # be explicit about what 1 means
model = glm(outcome.binary ~ age + gender + s100b + ndka, family = "binomial", data = aSAH)
preds = predict(model, type = "response")
roc_object = roc(aSAH$outcome.binary, preds, direction = "<")
roc_object$auc                                          # 0.8096
plot(roc_object)

Incremental value — refit without each biomarker in turn and difference the AUCs:

model.without.ndka = glm(outcome.binary ~ age + gender + s100b, family = "binomial", data = aSAH)
preds.without.ndka = predict(model.without.ndka, type = "response")
roc_object.without.ndka = roc(aSAH$outcome.binary, preds.without.ndka, direction = "<")
roc_object$auc - roc_object.without.ndka$auc            # 0.8096 - 0.7824 = 0.027

plot(roc_object, lwd = 2)
plot(roc_object.without.ndka, add = TRUE, col = 2, lwd = 2)
legend("bottomright", lty = c(1,1), col = c(1,2), c("ROC with ndka", "ROC without ndka"))
Model AUC $IV$
age + gender + s100b + ndka 0.810 —
drop ndka 0.782 0.027
drop s100b 0.709 0.100

So s100b carries real incremental value and ndka very little — consistent with the summary(model) table, where s100b is strongly significant ($p = 0.0001$) and ndka is not ($p = 0.08$). But per the Pepe point above, that agreement is not guaranteed: the coefficient test and the AUC difference answer different questions.

Prediction: cross-validation

Reading: Lecture Notes Part 1, §1.8.4 (pp. 43–49), following James, Witten, Hastie & Tibshirani, An Introduction to Statistical Learning. Code: R script — Lecture Video 24.

Everything in the previous section built a model and evaluated it on the same data. That is inappropriate — a kind of cheating: of course the tool performs well on the data used to make it. Ideally you build on training data and evaluate on totally separate test (validation) data.

Terminology, with “statistical learning method” standing for whatever model/method/tool produces the predictions:

  • Test error rate — the average error from predicting a new observation, one not used in training. This is the rate you actually care about; a method is warranted if its test error is low.
  • Training error rate — the error on the observations used to fit. Often quite different from the test error, and in particular it can dramatically underestimate it.

Usually there is only one dataset, so the test error is estimated by holding out a subset of the observations from the fit and predicting them.

Single holdout

Randomly split the data into a training set and a test set; fit on the first, score (e.g. AUC) on the second. Simple, but two drawbacks:

  1. Highly variable — do it twice and there is no guarantee the two estimates are close.
  2. Trained on less data — methods perform worse with fewer observations, so the estimate may over-estimate the test error.

Cross-validation addresses both.

Leave-one-out cross-validation (LOOCV)

Hold out a single observation $(\mathbf{Z}_1, Y_1)$, fit on the remaining $N-1$, predict $\hat{Y}_1$, and record the error $TE_1$ (AUC, $L_1$, $L_2$, …). Since the held-out point was not used in fitting, $TE_1$ is approximately unbiased for the test error — but a poor estimate, being based on one observation. So repeat for every observation and average:

\[\frac{1}{N}\sum_{i=1}^{N} TE_i\]

Advantages over a single split: far less bias (each fit uses $N-1$ observations, almost the whole dataset, instead of roughly half) and no randomness — LOOCV always gives the same answer, since there is nothing to randomize. Disadvantage: $N$ model fits, which can be expensive when $N$ is large or each fit is slow.

$k$-fold cross-validation

Randomly divide the observations into $k$ folds of roughly equal size. Hold out fold 1, fit on the other $k-1$, compute $TE_1$ on the held-out fold; repeat with each fold as the test set; average:

\[\frac{1}{k}\sum_{j=1}^{k} TE_j\]

LOOCV is the special case $k = N$. In practice $k = 5$ or $k = 10$, the advantage being mainly computational — $5$ or $10$ fits instead of $N$. Because the fold assignment is random, different runs give different answers.

Generalized cross-validation

Like the single split, but repeated: pick a proportion $h$ of the data as training (e.g. $h = 0.5$, or $0.75$), fit, score the rest as $TE_1$; repeat $M$ times ($M = 100$ or $500$) with fresh random subsets — which, unlike folds, are not mutually exclusive — and average:

\[\frac{1}{M}\sum_{m=1}^{M} TE_m\]

This lets you tune how much data goes to training vs. test and generally has less bias than $k$-fold, at more computational cost than $k$-fold (though usually less than LOOCV).

Two further points

  • With real data the true test error is unknown, so the accuracy of a CV estimate cannot be checked; with simulated data it can. Sometimes the actual level of the estimated test error matters (how will this tool do on new data?); sometimes only the location of the minimum across methods or tuning levels matters, and the absolute value does not.
  • CV is not specific to binary outcomes — continuous outcomes and multi-category (misclassification-rate) outcomes work the same way.

A note on AUC variance

For a standard error or confidence interval on an AUC, the asymptotic theory is DeLong, DeLong & Clarke-Pearson (1988); pROC::ci.auc() implements it (method = "delong") and a bootstrap alternative (method = "bootstrap"). Caveat: neither accounts for $p$ itself being an estimate with its own variability. To capture that you must bootstrap the entire procedure — (a) build the model, (b) evaluate it — which ci.auc does not do and has to be coded by hand.

In R

From R script — Lecture Video 24, continuing with aSAH. Reference point — the training performance from the last section:

model = glm(outcome.binary ~ age + gender + s100b + ndka, family = "binomial", data = aSAH)
preds = predict(model, type = "response")
roc(aSAH$outcome.binary, preds, direction = "<")$auc       # 0.8096 — training AUC

Single split. sample() picks half the row indices; setdiff() gives the rest. The one new ingredient is predict(..., newdata = test.data) — score rows the model never saw:

total.n = length(aSAH$age)
training = sample(total.n, round(0.5*total.n), replace = F)
test = setdiff(1:total.n, training)
training.data = aSAH[training,];  test.data = aSAH[test,]

model = glm(outcome.binary ~ age + gender + s100b + ndka, family = "binomial", data = training.data)
preds = predict(model, newdata = test.data, type = "response")
roc(test.data$outcome.binary, preds, direction = "<")$auc  # an estimate of test performance

Run it twice and the AUC moves — that is drawback 1 above.

LOOCV. An AUC on one observation makes no sense, so the per-observation error is the Brier score term $(p_i - A_i)^2$:

test.error.vector = vector(length = total.n)
for(i in 1:total.n){
    training.data = aSAH[-i,];  test.data = aSAH[i,]
    model = glm(outcome.binary ~ age + gender + s100b + ndka, family = "binomial", data = training.data)
    preds = predict(model, newdata = test.data, type = "response")
    test.error.vector[i] = (preds - test.data$outcome.binary)^2
}
mean(test.error.vector)                                    # 0.1842 — everyone gets this

No sample() anywhere, so the answer is the same on every machine.

5-fold CV, back to AUC. Fold membership is assigned by sample(1:5, total.n, replace = T) — one of several ways to do it, giving folds of approximately equal size:

k = 5
test.error.vector = vector(length = k)
assign = sample(1:5, total.n, replace = T)
for(i in 1:k){
    test = which(assign == i);  training = setdiff(1:total.n, test)
    training.data = aSAH[training,];  test.data = aSAH[test,]
    model = glm(outcome.binary ~ age + gender + s100b + ndka, family = "binomial", data = training.data)
    preds = predict(model, newdata = test.data, type = "response")
    test.error.vector[i] = roc(test.data$outcome.binary, preds, direction = "<")$auc
}
mean(test.error.vector)                                    # varies by run

Everyone gets a different answer here. The script suggests swapping the roles — train on one fold, test on four — to watch model building struggle with 1/5 of the data.

Generalized CV with $h = 0.75$, $M = 100$: the single-split code inside a loop, with round(h*total.n) training rows each time:

M = 100;  h = 0.75
test.error.vector = vector(length = M)
for(i in 1:M){
    training = sample(total.n, round(h*total.n), replace = F)
    test = setdiff(1:total.n, training)
    ...                                                    # fit / predict / roc as above
    test.error.vector[i] = roc_object$auc
}
mean(test.error.vector)

Expect glm.fit: fitted probabilities numerically 0 or 1 occurred warnings along the way — with 113 rows and a strong biomarker, some training subsets separate almost perfectly.

AUC confidence interval on a single split:

ci.auc(roc_object, method = "delong")
ci.auc(roc_object, method = "bootstrap")   # 2000 stratified bootstrap replicates

The two agree closely, and both are wide — roughly $\pm 0.1$ around the point estimate with ~56 test rows — and both ignore the model-building variability, as noted above.

Class activity: medicaldata::opt. A randomized trial of periodontal treatment in pregnancy (823 women, 4 centres); outcome low birthweight by the WHO cutoff:

library(medicaldata)
opt.data = opt[!is.na(opt$Birthweight), ]
opt.data$lowbw = 1*(opt.data$Birthweight < 2500)
  1. Logistic model: lowbw ~ Group + Age + Black + Nat.Am + Asian + Education + Hypertension + Diabetes.
  2. Training AUC — predict() then roc(..., direction = "<"). Note that with missing predictor values predict() only returns rows that are complete, so the response vector and the predictions will not line up until missingness is dealt with.
  3. 5-fold CV AUC — the aSAH loop with the formula and data swapped.
  4. Incremental value under CV — the script picks BL.Anti.inf (anti-inflammatory use in the 6 months before baseline; no missing data). The point of set.seed(50) before each of the two 5-fold runs is that both models are scored on the same folds, so the difference in CV AUCs is a fair comparison rather than fold noise.
Quantity Value
training AUC, base model 0.608
5-fold CV AUC, base model (seed 50) 0.562
5-fold CV AUC, + BL.Anti.inf (seed 50) 0.578
incremental value 0.016

The gap between the training AUC and the CV AUC is the underestimate of test error that motivates this whole section; the incremental value of the biomarker is small.

Prediction of time-to-event outcomes

Reading: Lecture Notes Part 1, §1.8.5 (pp. 49–54). Code: R script — Lecture Video 26.

Back to survival. The prediction is now $p = P(T \le t)$ for a chosen horizon $t$, with $\hat{p}_i$ coming from a Cox model, and the question is how to evaluate it when the outcome is a censored event time rather than a 0/1. The answer is time-dependent sensitivity, specificity and AUC — the same quantities as before with $D^{\pm}$ replaced by “had the event by $t$” vs. “survived past $t$.”

Three definitions

Who counts as a case and who as a control depends on how you treat the time axis.

Definition Cases ($D^{+}$) Controls ($D^{-}$) A subject’s role
Cumulative / dynamic $T \le t$ $T > t$ control while $t < T_i$, then case for all $t \ge T_i$
Incident / dynamic $T = t$ $T > t$ control early, case only at $T_i$ — mirrors the partial likelihood’s risk sets
Incident / static $T = t$ $T > t^{}$ for a fixed long-term $t^{}$ fixed; cases stratified by event time, controls event-free through $(0, t^{*})$
\[\text{Spec}(\theta, t) = P(p < \theta \mid T > t), \qquad \text{Sens}(\theta, t) = P(p \ge \theta \mid T \le t) \qquad \text{(cumulative/dynamic)}\]

The course uses cumulative/dynamic: at each $t$ the whole population is split by status at $t$, which is what matters clinically when a specific horizon (5-year survival, 30-day mortality) is the question. Then

\[\text{ROC}(u, t) = \text{Sens}\{\text{Spec}^{-1}(1-u),\, t\}, \qquad \text{AUC}(t) = \int \text{ROC}(u, t)\, du\]

and the AUC keeps its pairwise interpretation,

\[\text{AUC}(t) = P(p_i > p_j \mid T_i \le t,\ T_j > t), \quad i \ne j\]

— given one subject who died before $t$ and one who survived past it, the probability $p$ ranks the one who died as higher risk.

Estimation without censoring

If nobody were censored before $t$, the empirical estimates are exactly the earlier ones with $I(T_i \le t)$ selecting cases and $I(T_j > t)$ selecting controls:

\[\widehat{\text{sens}}(\theta, t) = \frac{\sum_i I(p_i \ge \theta) I(T_i \le t)}{\sum_i I(T_i \le t)}, \qquad \widehat{\text{spec}}(\theta, t) = \frac{\sum_i I(p_i < \theta) I(T_i > t)}{\sum_i I(T_i > t)}\] \[\widehat{\text{AUC}}(t) = \frac{\sum_i \sum_j I(p_i > p_j)\, I(T_i \le t)\, I(T_j > t)} {\sum_i \sum_j I(T_i \le t)\, I(T_j > t)}\]

Estimation with censoring: IPCW

In practice we observe $X_i = \min(T_i, C_i)$, not $T_i$, and three situations arise:

  • $X_i > t$ — then $T_i > t$ for sure: a control.
  • $X_i \le t$ and $\delta_i = 1$ — then $T_i \le t$ for sure: a case.
  • $X_i \le t$ and $\delta_i = 0$ — censored before $t$: status unknown.

Dropping the third group loses information at best and induces bias at worst. The inverse probability of censoring weighting (IPCW) estimator instead gives the subjects whose status is known extra weight to stand in for those whose status is not. Define $G(t) = P(C > t)$, the survival function of the censoring time, estimated by Kaplan–Meier with the event indicator flipped ($1 - \delta_i$). Each subject’s weight is

\[\widehat{W}_i = \frac{I(X_i > t)}{\hat{G}(t)} + \frac{I(X_i \le t)\, \delta_i}{\hat{G}(X_i)}\]
Subject $\widehat{W}_i$ Reasoning
$X_i > t$ $1 / \hat{G}(t)$ made it past $t$ uncensored — “extra credit” inversely proportional to how likely that was
$X_i \le t$, $\delta_i = 1$ $1 / \hat{G}(X_i)$ had the event at $X_i$, so was uncensored up to $X_i$; weight by the probability of that
$X_i \le t$, $\delta_i = 0$ $0$ status unknown — excluded from the sums, but still used in estimating $\hat{G}(\cdot)$

The estimates are the empirical ones with $T \to X$ and a weight in every sum:

\[\widehat{\text{sens}}(\theta, t) = \frac{\sum_i \widehat{W}_i\, I(p_i \ge \theta)\, I(X_i \le t)}{\sum_i \widehat{W}_i\, I(X_i \le t)}, \qquad \widehat{\text{spec}}(\theta, t) = \frac{\sum_i \widehat{W}_i\, I(p_i < \theta)\, I(X_i > t)}{\sum_i \widehat{W}_i\, I(X_i > t)}\] \[\widehat{\text{AUC}}(t) = \frac{\sum_i \sum_j \widehat{W}_i \widehat{W}_j\, I(p_i > p_j)\, I(X_i \le t)\, I(X_j > t)} {\sum_i \sum_j \widehat{W}_i \widehat{W}_j\, I(X_i \le t)\, I(X_j > t)}\]

Everything said about incremental value and cross-validation carries over unchanged; variances are usually obtained by bootstrap or wild bootstrap rather than analytically. References: Heagerty & Zheng (2005), Heagerty, Lumley & Pepe (2000), Uno et al. (2007).

In R

From R script — Lecture Video 26. BrainCancer again, with a 50/50 split (set.seed(100)), a Cox model fit on the training half, and predictions for $t = 30$ months on the test half via the Video 12 recipe — basehaz(centered = FALSE), predict(type = "lp", reference = "zero"), $\hat{S}(30 \mid Z) = \exp{-\hat{\Lambda}_0(30)\, e^{\hat\beta^T Z}}$ — so that pred.event = 1 - pred.surv is the $\hat{p}_i = \hat{P}(T \le 30 \mid Z_i)$ to evaluate.

IPCW by hand, exactly the formulas above with $\theta = 0.5$:

test.data$censor = 1 - test.data$status                 # flip the indicator ...
km = survfit(Surv(test.data$time, test.data$censor) ~ 1) # ... KM for G(t) = P(C > t)
G.t  = approx(km$time, km$surv, xout = t,              method = "constant")$y
G.Xi = approx(km$time, km$surv, xout = test.data$time, method = "constant")$y

test.data$Wi = 0
test.data$Wi[test.data$time >  t]                        = 1/G.t
test.data$Wi[test.data$time <= t & test.data$status == 1] = 1/G.Xi[test.data$time <= t & test.data$status == 1]

sens = sum(Wi*(pred.event >= theta)*(time <= t)) / sum(Wi*(time <= t))   # 0.3507
spec = sum(Wi*(pred.event <  theta)*(time >  t)) / sum(Wi*(time >  t))   # 0.7059

With timeROC, which does the same thing (weighting = "marginal" is the KM estimate of $G$) and returns PPV/NPV as well:

library(timeROC)
SeSpPPVNPV(theta, test.data$time, test.data$status, pred.event,
           cause = 1, weighting = "marginal", times = 30)

     Cases Survivors Censored Se (%) Sp (%) PPV (%) NPV (%)
t=30    13        17       13  35.07  70.59   41.43    64.7

Same 35.07 / 70.59 as by hand; the Censored column is the 13 test subjects with $\widehat{W}_i = 0$. times = c(10, 30, 50) gives one row per horizon.

roc.object = timeROC(test.data$time, test.data$status, pred.event,
                     cause = 1, weighting = "marginal", times = c(10, 30, 50), ROC = TRUE)
roc.object            # AUC(10) = 85.3%,  AUC(30) = 68.8%,  AUC(50) = 62.0%
plot(roc.object, time = 30)
plotAUCcurve(roc.object)

timeROC prints AUC (%) per time; roc.object$AUC is a vector whose first element is for $t = 0$ (always NA), hence $AUC[2] in the loops below. Discrimination falls as the horizon lengthens — the model separates early deaths well, later ones less so.

Class activity: Paquid (ships with timeROC): 2561 elderly subjects, time to dementia, with DSST (Digit Symbol Substitution Test) and MMSE (Mini-Mental State) scores at baseline. status is 0 / 1 / 2 = censored / dementia / death without dementia, so the first move is status.use = 1*(status == 1) — treating death as censoring, the “simplistic” competing-risks approach of §1.10.

  1. Training AUC at $t = 7$ years — fit Surv(time, status.use) ~ DSST + MMSE on everyone, predict $\hat{p}_i(7)$ for everyone, timeROC(..., times = 7).
  2. Test AUC by generalized CV — $h = 0.9$, $M = 100$, set.seed(50): the Cox fit, basehaz, prediction, and timeROC all move inside the loop and run on the held-out 10%.
  3. Incremental value of DSST — same loop, now fitting both DSST + MMSE and MMSE alone on each training draw and scoring both on the same test draw.
Quantity AUC(7)
training, DSST + MMSE 0.781
generalized CV, DSST + MMSE 0.779
generalized CV, MMSE only 0.712
incremental value of DSST 0.068

With $n = 2561$ and two predictors the optimism is negligible (0.781 vs. 0.779); the DSST adds real discrimination on top of the MMSE. Both models were evaluated on identical splits because the seed is set once and the two fits sit inside one loop.

Time-varying covariates

Reading: Lecture Notes Part 1, §1.9 (pp. 55–60). Code: R script — Lecture Video 29.

So far every covariate in $\lambda(t \mid Z) = \lambda_0(t)\exp{\beta^T Z}$ was measured at baseline. Many are not: pollution exposure, blood glucose, medication adherence, alcohol use, or — the running example — whether a patient has received a heart transplant, which is 0 for everyone at enrolment and switches to 1 for some during follow-up. These are time-varying (time-dependent) covariates.

The model

Define, for the transplant example,

\[Z_1(t) = \begin{cases} 0 & \text{no transplant by time } t \\ 1 & \text{transplant by time } t \end{cases}\]

and write

\[\lambda(t \mid Z) = \lambda_0(t)\exp\{\beta_1 Z_1(t) + \beta_2 Z_2 + \beta_3 Z_3 + \cdots\} \;=\; \lambda_0(t)\exp\{\beta^T Z(t)\}\]

with some components of $Z(t)$ varying and others fixed. Two things to notice:

  • It is still a proportional hazards model — $\beta$ does not depend on time; the covariate values do. Interpretation of $e^{\beta_1}$ follows the usual rubric.
  • Estimation uses the same partial likelihood, with the covariates evaluated at the current event time inside each risk set:
\[L^{partial}(\beta) = \prod_{j=1}^{K} \frac{\exp\{\beta^T z_j(t_{(j)})\}}{\sum_{l \in R_j} \exp\{\beta^T z_l(t_{(j)})\}}\]

As a subject moves through time, each time they are compared to the person who had the event, their covariate value must be the current one, $z_l(t_{(j)})$.

Setting up the data: splitting person-time

Operationally this means one row per interval of constant covariate value. Patient 1 in the notes’ table had a transplant at day 11 and died at day 57, so $Z_1(t) = 0$ for $t < 11$ and $1$ for $t \ge 11$, and the record becomes

ID Age Prior Transplant Start End Status
1 42.5 0 0 0 11 0
1 42.5 0 1 11 57 1

Time-invariant covariates repeat; the event indicator is 1 only on the final row. Once the data are in (start, stop, event) form, standard software fits the model with no further modification.

Do the multiple rows per subject mean correlated data? No. The split is a programming trick: the likelihood at any event time uses exactly one row per subject — the one whose interval contains that time. Two genuine exceptions need a cluster variance: subjects with multiple events, and subjects appearing in overlapping intervals (almost always a data error).

The Stanford heart transplant fallacy

Clark et al. (1971) plotted Kaplan–Meier curves for transplanted vs. non-transplanted patients and concluded transplants dramatically prolonged life. The problem: receipt of a transplant is time-dependent — patients had to survive long enough to receive one. The comparison only shows that patients who live longer live longer, a tautology. Two distinct errors:

  1. Transplant was modelled as a baseline binary covariate instead of a time-varying one — fixable by the person-time split above, after which you can talk about the hazard.
  2. Talking about survival for the transplant group is prediction of the future using data from the future. PREDICTION IS BROKEN with time-varying covariates: you cannot predict survival from covariate values that have not happened yet.

Gail (1972) made the reassessment; the lesson is the rule above.

Landmark prediction

If the goal is prediction or a statement about survival, pick a landmark time $t_0$ and split the sample:

  • Group A — event before $t_0$ ($T \le t_0$). Discarded.
  • Group B — survived past $t_0$ ($T > t_0$). The model is built here only.

The baseline moves to $t_0$, so everything known at $t_0$ — including a time-varying covariate’s value up to then — becomes an ordinary baseline covariate for Group B. Heart transplant data, $t_0 = 30$ days: 79 of 103 patients survived 30 days; of those, 33 had a transplant within 30 days and 46 had not (30 of whom got one later, but still count as “no transplant by day 30”). The resulting transplant30 is fixed for every 30-day survivor.

Prediction and all accuracy measures are then conditional on $T > t_0$:

\[\text{Spec}(\theta, t, t_0) = P(p < \theta \mid T > t), \qquad \text{Sens}(\theta, t, t_0) = P(p \ge \theta \mid T \le t,\ T > t_0), \qquad p = P(T \le t \mid T > t_0)\]

Two disadvantages: how to choose $t_0$, and nothing can be said about Group A. (Dynamic prediction — van Houwelingen & Putter 2011 — extends the idea to a sequence of landmarks.)

In R

From R script — Lecture Video 29, on survival::jasa — the Stanford data: 103 patients, futime/fustat, baseline age and surgery (prior bypass), and transplant with wait.time (NA for the 34 never transplanted).

The problematic analysis, reproducing Clark et al.:

km = survfit(Surv(futime, fustat) ~ transplant, data = jasa)
plot(km, col = c(1,2), lwd = 2)
survdiff(Surv(futime, fustat) ~ transplant, data = jasa)   # Chisq = 33.2, p = 8e-09

A hugely significant benefit of transplant — the tautology.

Landmark at 30 days. Keep the 79 patients with futime >= 30, define the fixed covariate, and compare:

groupb = jasa[jasa$futime >= 30,]
groupb$transplant30 = 1*(groupb$transplant == 1 & groupb$wait.time < 30)
survdiff(Surv(futime, fustat) ~ transplant30, data = groupb)   # Chisq = 0, p = 0.8
coxph(Surv(futime, fustat) ~ transplant30 + age + surgery, data = groupb)
Covariate HR $p$
transplant30 0.96 0.88
age 1.04 0.03
surgery 0.44 0.05

Conditional on surviving 30 days, transplant makes no difference; age and prior surgery do. The & with wait.time < 30 works because R’s NA & FALSE is FALSE, so the never-transplanted get 0 (not NA).

Time-varying covariate with tmerge. Back to all 103, this time letting transplant switch on at wait.time:

jasa2 = jasa
jasa2$id = 1:nrow(jasa)
jasa2$futime[jasa2$futime == 0] = 0.01            # tmerge refuses a zero-length interval
jasa2 = tmerge(jasa2, jasa2, id = id,
               death  = event(futime, fustat),    # the outcome
               transpl = tdc(wait.time))          # tdc = time-dependent covariate

tmerge builds the (tstart, tstop) rows. Patient 4 (transplant at 35, died at 38) and patient 7 (transplant at 50, died at 674) each become two rows:

   fustat futime wait.time transplant id tstart tstop death transpl
4       1     38        35          1  4      0    35     0       0
5       1     38        35          1  4     35    38     1       1
8       1    674        50          1  7      0    50     0       0
9       1    674        50          1  7     50   674     1       1

— exactly the person-time split from the notes. The model then takes the counting-process Surv(tstart, tstop, event) form:

coxph(Surv(tstart, tstop, death) ~ transpl + age + surgery, data = jasa2)
Covariate HR $p$
transpl 0.92 0.80
age 1.03 0.02
surgery 0.46 0.03

n = 169 rows for 103 patients, 75 events. Properly modelled, the hazard of death after transplant is 0.92 times the hazard before, nowhere near significant — the headline benefit of the naive analysis was entirely immortal-time bias. The two approaches agree, and both only license statements about the hazard; the landmark version is the one that could also be used for prediction, conditional on 30-day survival.