Chapter 17: Causal Survival Analysis

Published

Last modified: 2026-10-09 09:46:54 (UTC)

📝 Preview Changes: This page has been modified in this pull request (~0% of content changed).
🎨 Highlighting Legend: Modified text (yellow) shows changed words/phrases, added text (green) shows new content, and new sections (blue) highlight entirely new paragraphs.

So far the outcomes have been measured at one fixed time, such as weight gain by 1982 in the smoking cessation example. Many causal questions are instead about how long it takes until an event happens, for example how quitting smoking affects the time until death. Questions of this kind call for survival analysis.

Despite the name, the event need not be death: “survival analysis” (also called “failure time analysis”) covers any time-to-event outcome, whether the event is death, marriage, incarceration, a cancer diagnosis, or a flu infection. What sets these analyses apart is that many individuals will not have had the event when the study ends, so their event times are unknown. This chapter covers the basic methods for the simple case of a treatment that is fixed at baseline.

This chapter is based on Hernán and Robins (2020, chap. 17, pp. 227-240).

R and Stata code for this chapter’s programs (Programs 17.x) is in Tom Palmer’s cibookex-r companion (GPL-3.0), which descends from the authors’ own code.

Key challenge: Survival outcomes involve time, and individuals may be censored before experiencing the event. Estimating hazards is often a convenient way to get at survivals and risks. This chapter shows how IP weighting and the parametric g-formula extend naturally to survival outcomes.

1 17.1 Hazards and Risks (pp. 227-229)


Let \(A\) indicate smoking cessation (1: yes, 0: no) and let \(T\) be the time to death, counted from the start of follow-up. The target is the average causal effect of \(A\) on \(T\), an outcome that can occur at any point after follow-up begins.

Definition 1 (Administrative censoring) Usually only a minority of participants have the event during the study, because follow-up stops for everybody on a fixed date, the administrative end of follow-up. For anyone still event-free on that date we only know that the event time exceeds it: their survival time is administratively censored.

Example 1 (NHEFS: Administrative censoring) In the NHEFS analysis, the population is the 1629 smokers aged 25-74 at baseline who were still alive in 1982. Follow-up runs from January 1, 1983 to December 31, 1992, so everyone has the same administrative censoring time of 120 months. Only 318 of the 1629 died during this period; the other 1311 have administratively censored survival times.

Any survival analysis faces administrative censoring. With staggered entry, where people enter the study on different dates, administrative censoring times differ between individuals even if the study closes on a single date.

Survival analyses may also face other kinds of censoring, such as loss to follow-up or competing events (Fine Point 17.1). Earlier chapters showed how standardization or IP weighting can correct the selection bias such censoring causes, and the same tools carry over to survival outcomes. This chapter deals only with administrative censoring; other censoring usually changes over time, whereas the administrative censoring time is known at baseline, so it is postponed to Part III.

The book simplifies by treating everyone without a confirmed death as alive at the end of follow-up (some deaths may in fact have gone unconfirmed) and by setting aside the issue raised in Fine Point 12.1.

In the NHEFS example, the month of death \(T\) can take values from 1 (January 1983) to 120 (December 1992). \(T\) is known for 102 treated (\(A = 1\)) and 216 untreated (\(A = 0\)) individuals who died during follow-up (\(102 + 216 = 318\)), and is administratively censored (all we know is \(T > 120\)) for the remaining 1311. Therefore we cannot estimate the mean survival \(\operatorname{E}\mathopen{}\left[T\right]\mathclose{}\) as we did with the outcomes of previous chapters; we need measures that accommodate administrative censoring.

1.1 Survival, Risk, and Hazard

Definition 2 (Survival, risk, and hazard of an event time) The survival at month \(k\), \(\Pr[T > k]\), is the share of individuals still alive after \(k\). Computing it for every month up to the administrative end of follow-up, \(k_\text{end} = 120\), and plotting it against time gives the survival curve, which equals \(\Pr[T > 0] = 1\) at \(k = 0\) and can only go down as \(k = 1, 2, \ldots, k_\text{end}\) increases.

The risk, or cumulative incidence, at \(k\) is the complement of the survival:

\[\Pr[T \leq k] = 1 - \Pr[T > k]\]

Its curve begins at \(\Pr[T \leq 0] = 0\) and can only go up.

The hazard at \(k\) is the share of individuals who have the event at \(k\) among those who were still event-free just before \(k\):

\[\Pr[T = k \mid T > k - 1]\]

Remark 1 (Discrete-time hazards). Strictly speaking, the hazard of Definition 2 is a discrete-time hazard, defined for time measured in intervals rather than continuously. Since real data always record time in units such as years, months, or days, we simply call it the hazard.

A natural summary of the treatment effect in a survival analysis contrasts the survival, or the risk, under each treatment level at one or more times \(k\). Survival curves also yield other measures, such as the restricted mean survival time or the years of life lost.

If we suppose (provisionally, until Section 17.4) that quitters (\(A = 1\)) and non-quitters (\(A = 0\)) are marginally exchangeable, we can compare \(\Pr[T > k \mid A = 1]\) with \(\Pr[T > k \mid A = 0]\) for all \(k\) (Figure 17.1 in the book). At 120 months, survival was 76.2% for quitters versus 82.0% for non-quitters. Equivalently, the 120-month risk was \(100\% - 76.2\% = 23.8\%\) among quitters and \(100\% - 82.0\% = 18.0\%\) among non-quitters. A log-rank test comparing the two survival curves gave \(P = 0.005\). These figures are not adjusted for confounding and thus do not have a causal interpretation; see Section 17.4 for the adjusted estimates.

At 120 months the hazard was 0% among quitters (their last death occurred at month 113) and 1/986 = 0.10% among non-quitters; over the 120 months, the hazard curves looked roughly like the letter M.

1.2 Risk vs. Hazard

Remark 2 (Risk versus hazard). Risk and hazard are built differently:

  • Risk: the denominator, everyone at baseline, is the same at every \(k\), and the numerator accumulates all events from baseline through \(k\), so the risk never decreases.
  • Hazard: the denominator, those still alive at \(k\), shrinks over time, and the numerator counts only the events at \(k\), so the hazard can rise or fall.

Many survival analyses summarize the effect by the hazard ratio, treated versus untreated. Fine Point 17.2 explains why that summary is problematic, so the book’s analyses focus on survival and risk instead. Hazards still matter, though: estimating them is often the easiest route to survival and risk.

NoteFine Point 17.1: Competing Events

A competing event (Section 8.5) is one, usually death, that makes the event of interest impossible afterwards: someone who dies of cancer can no longer have a stroke. The central choice is whether to treat competing events as non-administrative censoring.

  • Treated as censoring: the analysis tries to mimic a world in which death from other causes has been eliminated, or made unrelated to stroke risk factors. The resulting estimand is hard to interpret and may not be meaningful (Chapter 8), and the censoring can create selection bias even under the null, which then has to be removed (for instance by IP weighting) using measured risk factors for stroke.

  • Not treated as censoring: the event time is set to infinity for those who die, so they have zero probability of stroke between death and the end of follow-up. A non-null effect on stroke is then hard to interpret, because it could arise solely because treatment changes the risk of death, which precludes stroke.

Other options have drawbacks too. A composite outcome (stroke or death) removes the competing event but answers a different question, and a non-null estimate could reflect an effect on either component. Restricting attention to the principal stratum of people who would survive under either treatment targets a local effect of the kind seen in Chapter 16, which is hard both to interpret and to estimate.

No statistical technique resolves the problem, because competing events make the causal estimand itself unclear. Young et al. (2019) review the available approaches and their difficulties, and Stensrud et al. (2020, 2021) propose separable effects (see Technical Point 23.3 of the book).

2 17.2 From Hazards to Risks (pp. 229-232)


Definition 3 (Person and person-time data layouts) A survival dataset can be laid out in two ways.

  • One row per person: the layout used so far in the book, called the “wide” format when treatments and confounders vary over time. In the previous section the data had 1629 rows.
  • One row per person-time: person 1 at \(k = 0\), then person 1 at \(k = 1\), and so on until that person’s follow-up stops, then person 2, and so on. This “long” layout is used for most of this chapter and for all of Part III’s time-varying treatments. In the smoking cessation data it has 176,764 rows, one per person-month.

2.1 The Event Indicator \(D_k\)

Definition 4 (Event indicator) The long format records survival through a time-varying event indicator \(D_k\): for each person and month \(k\), \(D_k = 1\) if \(T \leq k\) and \(D_k = 0\) if \(T > k\).

The book’s Figure 17.2 is a causal diagram with treatment \(A\), event indicators \(D_1\) and \(D_2\), and \(U\), usually unmeasured, for whatever makes a person more susceptible to the event; susceptibility can itself change over time (\(U_0, U_1, \ldots\)).

Each person-time row at \(k\) stores the next month’s indicator \(D_{k+1}\): the \(k = 0\) row holds \(D_1\) (1 for a death in month 1), the \(k = 1\) row holds \(D_2\), and so on. A person’s rows stop at the first one with \(D_{k+1} = 1\), or else at month 119. Inclusion requires being alive at month 0, so \(D_0 = 0\) for everyone.

Using the time-varying outcome variable \(D_k\), we can define:

Definition 5 (Survival, risk, and hazard via the event indicator) These are the quantities of Definition 2, rewritten in terms of \(D_k\):

  • Survival at \(k\): \(\Pr[D_k = 0] = \Pr[T > k]\)
  • Risk at \(k\): \(\Pr[D_k = 1] = \Pr[T \leq k]\)
  • Hazard at \(k\): \(\Pr[D_k = 1 \mid D_{k-1} = 0] = \Pr[T = k \mid T > k - 1]\), since \(D_{k-1} = 0\) says \(T > k - 1\), and, with \(T\) measured in whole months (Remark 1), adding \(D_k = 1\) says exactly that \(T = k\)

At \(k = 1\) hazard and risk coincide, since everyone is alive at \(k = 0\) by construction.

2.2 From Hazards to Survival

Surviving through \(k\) means surviving each interval up to \(k\) in turn, so the survival at \(k\) is the product of the interval-specific conditional survivals.

Proposition 1 (Survival as a product of conditional survivals) Suppose \(D_0 = 0\) for everyone and that \(D_{m-1} = 1\) implies \(D_m = 1\) for every \(m \geq 1\). Then for every \(k \geq 1\) such that the conditioning events have positive probability,

\[\Pr[D_k = 0] = \prod_{m=1}^k \Pr[D_m = 0 \mid D_{m-1} = 0]\]

Proof. An event that has happened stays happened, so \(D_m = 0\) implies \(D_{m-1} = 0\). Hence the event \(\{D_k = 0\}\) equals \(\{D_1 = 0, D_2 = 0, \ldots, D_k = 0\}\). By the chain rule of probability, its probability is \(\Pr[D_1 = 0] \times \prod_{m=2}^k \Pr[D_m = 0 \mid D_{m-1} = 0, \ldots, D_1 = 0]\). Again because \(D_{m-1} = 0\) implies \(D_1 = \cdots = D_{m-2} = 0\), each conditioning event reduces to \(D_{m-1} = 0\). Finally, \(D_0 = 0\) for everyone, so \(\Pr[D_1 = 0] = \Pr[D_1 = 0 \mid D_0 = 0]\), which is the \(m = 1\) factor.

Each factor is one minus the hazard at \(m\), so the hazards through \(k\) determine the survival at \(k\), and hence the risk at \(k\) as well.

Example 2 (Survival through two months) For \(k = 2\), \(\Pr[D_2 = 0] = \Pr[D_1 = 0] \times \Pr[D_2 = 0 \mid D_1 = 0]\): survive the first month, then survive the second given survival of the first.

2.3 Nonparametric and Parametric Estimation

Definition 6 (Kaplan-Meier estimator) Nonparametrically, the hazard \(\Pr[D_k = 1 \mid D_{k-1} = 0]\) is estimated by the number of events in interval \(k\) divided by the number of people alive at the end of interval \(k - 1\). Plugging these estimates into the product gives the Kaplan-Meier (product-limit) estimator of the survival \(\Pr[D_k = 0]\), which the book used for Figure 17.1.

It estimates the survival curve well when the total number of events is reasonably large. Each interval, however, usually contains few or no events, so the individual hazard estimates are erratic, and estimating the hazard at a specific \(k\) may require parametric smoothing (Chapter 11 and Fine Point 17.3).

Example 3 (NHEFS: A logistic hazards model) Parametrically, a convenient choice is a logistic regression for \(\Pr[D_{k+1} = 1 \mid D_k = 0]\), fit to the person-time rows of those still alive at \(k\). For the treated and the untreated in our example:

\[\operatorname{logit}\,\Pr[D_{k+1} = 1 \mid D_k = 0, A] = \theta_{0,k} + \theta_1 A + \theta_2 (A \times k) + \theta_3 (A \times k^2)\]

with a time-varying intercept \(\theta_{0,k} = \theta_0 + \theta_4 k + \theta_5 k^2\). The treatment-by-time product terms let the hazard ratio change over follow-up (Technical Point 17.1 explains why a logistic model approximates a hazards model).

Multiplying the fitted values of one minus the hazard over time, separately for each treatment group, gives estimates of the survival \(\Pr[D_{k+1} = 0 \mid A = a]\). These curves (the book’s Figure 17.4) are smoothed versions of the Kaplan-Meier curves of Figure 17.1.

WarningThe Hazards Model Must Be Correct

This procedure is valid only if the hazards model is correctly specified, which looks plausible here because the parametric and nonparametric survival curves nearly coincide. Bootstrapping individuals gives 95% confidence intervals for the survival estimates.

People contribute many rows each, yet an ordinary logistic regression program reports correct standard errors as long as the hazards model is correct. Other link functions, such as the probit, would also work.

NoteFine Point 17.2: The Hazards of Hazard Ratios

Two features of the hazard ratio complicate its use as a causal effect measure.

It changes over time. Hazards vary with \(k\), and so does their ratio. Published analyses nonetheless often report one hazard ratio, typically from a Cox model with no treatment-by-time interaction, which forces the ratio to be constant. That single number is a weighted average of the time-specific ratios, which makes it hard to interpret. When the event is rare and the only censoring is administrative at a common time \(k_\text{end}\), the weight for time \(k\) is proportional to the number of untreated events at \(k\) (formally, the conditional density at \(k\) of \(T\) given \(A = 0\) and \(T < k_\text{end}\)). Being an average, it can equal 1 even when the survival curves differ. Survival and risk contrasts, by contrast, refer to a stated period, such as a 5-year survival difference or a 120-month risk ratio.

Even time-specific ratios are hard to interpret causally. Imagine a treatment that kills every high-risk person by time \(k\) and does nothing to anyone else. At \(k + 1\) the hazard ratio compares treated and untreated survivors of \(k\): only low-risk people remain among the treated, while the untreated survivors include both high-risk and low-risk people. The hazard ratio at \(k + 1\) is therefore below 1, although treatment helps nobody.

The explanation is selection bias from conditioning on survival, a variable affected by treatment. The hazard at time 2 is \(\Pr[D_2 = 1 \mid D_1 = 0, A]\), and conditioning on the collider \(D_1\) generally opens the path \(A \to D_1 \leftarrow U \to D_2\), creating an association between \(A\) and \(D_2\) among those with \(D_1 = 0\). This bias does not arise if treatment has no effect on the event, so that the survival curves of treated and untreated coincide (no arrows from \(A\) into the event indicators; the book’s Figure 17.3 includes such arrows). Hernán (2010) gives an example.

NoteFine Point 17.3: Models for Survival Analysis

All survival methods must handle failure times that are censored by the administrative end of follow-up.

  • Nonparametric methods, such as Kaplan-Meier curves, make no assumption about the distribution of the censored failure times.
  • Parametric models posit a distribution (exponential, Weibull, and so on) for the failure times or the hazards; the logistic hazards model of the main text is one example.
  • Semiparametric models, such as the accelerated failure time (AFT) model and the Cox proportional hazards model, leave the shape of the baseline hazard (the hazard when all covariates equal zero) unspecified, but restrict in advance how the hazard at other covariate values relates to it.

For applied survival analysis see Hosmer, Lemeshow, and May (2008); for more formal treatments see Kalbfleisch and Prentice (2002) and Fleming and Harrington (2005).

NoteTechnical Point 17.1: Approximating the Hazard Ratio via a Logistic Model

The discrete-time hazard ratio \(\frac{\Pr[D_{k+1}=1 \mid D_k=0,A=1]}{\Pr[D_{k+1}=1 \mid D_k=0,A=0]}\) equals \(\operatorname{exp}\mathopen{}\left\{\alpha_1\right\}\mathclose{}\) at all times \(k+1\) in the hazards model \(\Pr[D_{k+1} = 1 \mid D_k = 0, A] = \Pr[D_{k+1} = 1 \mid D_k = 0, A = 0] \times \operatorname{exp}\mathopen{}\left\{\alpha_1 A\right\}\mathclose{}\).

Taking logs on both sides of the hazards model gives

\[\log \Pr[D_{k+1} = 1 \mid D_k = 0, A] = \alpha_{0,k} + \alpha_1 A, \quad \text{where } \alpha_{0,k} \stackrel{\text{def}}{=}\log \Pr[D_{k+1} = 1 \mid D_k = 0, A = 0].\]

Suppose the hazard at \(k + 1\) is small, i.e., \(\Pr[D_{k+1} = 1 \mid D_k = 0, A] \approx 0\). Then one minus the hazard, \(\Pr[D_{k+1} = 0 \mid D_k = 0, A]\), is close to 1, so the hazard is approximately equal to the odds:

\[\Pr[D_{k+1} = 1 \mid D_k = 0, A] \approx \frac{\Pr[D_{k+1} = 1 \mid D_k = 0, A]}{\Pr[D_{k+1} = 0 \mid D_k = 0, A]}.\]

Taking logs of both sides and substituting the log-linear hazards model:

\[\operatorname{logit}\,\Pr[D_{k+1} = 1 \mid D_k = 0, A] = \log \frac{\Pr[D_{k+1} = 1 \mid D_k = 0, A]}{\Pr[D_{k+1} = 0 \mid D_k = 0, A]} \approx \log \Pr[D_{k+1} = 1 \mid D_k = 0, A] = \alpha_{0,k} + \alpha_1 A\]

So when the hazard at \(k + 1\) is near zero, \(\theta_1\) in the logistic model \(\operatorname{logit}\,\Pr[D_{k+1} = 1 \mid D_k = 0, A] = \theta_{0,k} + \theta_1 A\) approximates the log hazard ratio \(\alpha_1\) (Thompson 1977). A common rule of thumb deems the approximation adequate if \(\Pr[D_{k+1} = 1 \mid D_k = 0, A] < 0.1\) at every \(k\).

The rare-event condition can nearly always be met by choosing a short enough time unit: years might do for lung cancer, whereas days might be needed for the common cold. Shorter units mean more person-time rows for fitting the logistic model.

3 17.3 Why Censoring Matters (pp. 232-234)


In the NHEFS example the only censoring is administrative, at the same time \(k_\text{end} = 120\) for everyone. Then the hazard-based procedure is more than needed: \(\Pr[D_{k+1} = 0 \mid A = a]\) can be estimated directly as the proportion of people with treatment \(a\) still alive at \(k + 1\), or with a separate logistic model for \(\Pr[D_{k+1} = 0 \mid A]\) at each \(k = 0, 1, \ldots, k_\text{end}\).

Things change with staggered entry.

Definition 7 (Staggered entry) Under staggered entry, people enter the study on different dates while the study closes on one date. Each person’s administrative censoring time is the closing date minus their entry date, so it now varies across individuals.

3.1 The Censoring Indicator

Definition 8 (Time-varying censoring indicator) The time-varying censoring indicator \(C_k\) equals 0 if the person’s administrative end of follow-up lies at or beyond \(k\) and 1 otherwise. In the person-time format, the row for an individual at time \(k\) includes \(C_{k+1}\).

The NHEFS dataset omits this variable because \(C_{k+1} = 0\) for everyone at all \(k\) before 120 months; with individual-specific administrative censoring, \(C_{k+1}\) switches from 0 to 1 at different times for different people.

Definition 9 (Survival without censoring) The target is the survival curve we would see if no one were censored before \(k_\text{end}\), that is, \(\Pr[D_k = 0 \mid A = a]\) with \(D_k\) observed even after censoring. In the notation of Chapter 12 this is \(\Pr[D^{\bar{c}=0}_k = 0 \mid A = a]\), where \(\bar{c} = (c_1, c_2, \ldots, c_{k_\text{end}})\) and the superscript \(\bar{c} = 0\) denotes no censoring; the superscript is dropped when the meaning is clear.

To keep things simple, suppose entry dates are effectively random, as they would be if no variable showed a secular trend. Censoring times, and so \(C\), are then independent of both treatment and time of death.

3.2 Why Naive Estimation Fails

WarningUncensored Survivors Do Not Estimate Survival

The proportion of people with treatment \(a\) who are both alive and uncensored through \(k\) does not estimate \(\Pr[D_k = 0 \mid A = a]\). It estimates the joint probability \(\Pr[C_{k+1} = 0, D_{k+1} = 0 \mid A = a]\) instead.

Example 4 (Naive estimation under censoring) To see why, consider a study with \(k_\text{end} = 2\) in which:

  • \(\Pr[C_1 = 0] = 1\): nobody is censored by \(k = 1\)
  • \(\Pr[D_1 = 0 \mid C_0 = 0] = 0.9\): 90% of individuals survive through \(k = 1\)
  • \(\Pr[C_2 = 0 \mid D_1 = 0, C_1 = 0] = 0.5\): a random half of survivors is censored by \(k = 2\)
  • \(\Pr[D_2 = 0 \mid C_2 = 0, D_1 = 0, C_1 = 0] = 0.9\): 90% of the remaining individuals survive through \(k = 2\)

The fraction of uncensored survivors at \(k = 2\) is \(1 \times 0.9 \times 0.5 \times 0.9 = 0.405\). Without censoring, that is, with \(\Pr[C_2 = 0 \mid D_1 = 0, C_1 = 0] = 1\), the survival would be \(1 \times 0.9 \times 1 \times 0.9 = 0.81\).

3.3 Correct Estimation Under As-If Randomly Assigned Censoring

Proposition 2 (Survival under as-if random censoring) Fix a treatment level \(a\) with \(\Pr[A = a] > 0\) and a time \(k \in \{1, \ldots, k_\text{end}\}\). Suppose the uncensored event indicators start at zero and never revert, that is, \(D^{\bar{c}=0}_0 = 0\) for everyone and \(D^{\bar{c}=0}_{m-1} = 1\) implies \(D^{\bar{c}=0}_m = 1\) for every \(m \geq 1\). Suppose also that, for every \(m = 1, \ldots, k\):

  • Censoring as if randomly assigned: given \(A\) and survival without censoring through \(m - 1\), being censored by \(m\) is independent of the event at \(m\) without censoring, \(C_m \perp\!\!\!\perp D^{\bar{c}=0}_m \mid A = a, D^{\bar{c}=0}_{m-1} = 0\);
  • Positivity: \(\Pr[C_m = 0 \mid A = a, D^{\bar{c}=0}_{m-1} = 0] > 0\) and \(\Pr[D^{\bar{c}=0}_{m-1} = 0 \mid A = a] > 0\);
  • Consistency: for anyone with \(C_m = 0\), \(D_m = D^{\bar{c}=0}_m\) and \(D_{m-1} = D^{\bar{c}=0}_{m-1}\). This is plausible because \(C_m = 0\) implies \(C_1 = \cdots = C_m = 0\): by Definition 8, a follow-up that ends at or beyond \(m\) also ends at or beyond every earlier month. So these people actually experienced the no-censoring regime through \(m\), and their observed \(D_{m-1}\) and \(D_m\) are the values that regime would produce.

Then the survival without censoring (Definition 9) at \(k\) is

\[\Pr[D^{\bar{c}=0}_k = 0 \mid A = a] = \prod_{m=1}^k \Pr[D_m = 0 \mid D_{m-1} = 0, C_m = 0, A = a]\]

Proof. Applying Proposition 1 to the uncensored event indicators \(D^{\bar{c}=0}_m\), within the stratum \(A = a\), gives

\[\Pr[D^{\bar{c}=0}_k = 0 \mid A = a] = \prod_{m=1}^k \Pr[D^{\bar{c}=0}_m = 0 \mid D^{\bar{c}=0}_{m-1} = 0, A = a]\]

Fix one factor \(m\). By the independence assumption (with positivity ensuring the added conditioning event has positive probability), adding \(C_m = 0\) to the conditioning event leaves the factor unchanged: it equals \(\Pr[D^{\bar{c}=0}_m = 0 \mid D^{\bar{c}=0}_{m-1} = 0, C_m = 0, A = a]\). Among the people with \(C_m = 0\), consistency replaces \(D^{\bar{c}=0}_m\) by \(D_m\) and \(D^{\bar{c}=0}_{m-1}\) by \(D_{m-1}\), so the factor equals \(\Pr[D_m = 0 \mid D_{m-1} = 0, C_m = 0, A = a]\), which involves only observed data. Multiplying the \(k\) factors gives the result.

The range \(k = 1, \ldots, k_\text{end}\) in Proposition 2 covers the same times as \(D_{k+1}\) for \(k = 0, \ldots, k_\text{end} - 1\) elsewhere in this chapter.

Definition 10 (Cause-specific hazard) The cause-specific hazard at \(k + 1\) is \(\Pr[D_{k+1} = 1 \mid D_k = 0, C_{k+1} = 0, A = a]\).

Estimation proceeds as before, except that the quantity estimated nonparametrically or with a logistic model is the cause-specific hazard (Definition 10).

WarningRandom Censoring Is Often Hard to Defend

Random censoring is often hard to defend. Under staggered entry, censoring time depends on the calendar date of entry (late entrants are followed for less time), and calendar time may also predict the outcome, so the analysis must adjust for calendar time at baseline. Other baseline prognostic factors may also differ between treatment groups and need adjustment.

The rest of the chapter adds adjustment for baseline confounders with g-methods, and Part III extends it to treatments and confounders that vary over time.

4 17.4 IP Weighting of Marginal Structural Models (pp. 234-235)


Example 5 (NHEFS: Confounding by age) Without exchangeability between treated and untreated, comparing their survival curves does not estimate a causal effect. In the smoking cessation data, 120-month survival was lower among quitters than among non-quitters (76.2% versus 82.0%), which need not mean that quitting raises mortality: quitters tend to be older, and older people die more often. Because quitters include more older people, confounding by age makes quitting look harmful.

4.1 Counterfactual Survivals

Definition 11 (Counterfactual survival) Let \(D^{a,\bar{c}=0}_k\) be the counterfactual death indicator at \(k\) under treatment \(a\) and no censoring. Since every analysis in this chapter targets survival without censoring, we abbreviate it to \(D^a_k\). From here on we also drop \(C_k = 0\) from the conditioning event of the hazard, as if everyone shared one administrative censoring time, as in the NHEFS. We want to compare the counterfactual survivals \(\Pr[D^{a=1}_{k+1} = 0]\) and \(\Pr[D^{a=0}_{k+1} = 0]\) under treatment for all (\(a = 1\)) and under no treatment for all (\(a = 0\)):

\[\Pr[D^{a=1}_{k+1} = 0] \text{ vs. } \Pr[D^{a=0}_{k+1} = 0] \quad \text{for } k = 0, 1, \ldots, k_\text{end} - 1\]

With confounding, the observed survivals \(\Pr[D_{k+1} = 0 \mid A = 1]\) and \(\Pr[D_{k+1} = 0 \mid A = 0]\) do not in general estimate this contrast.

4.2 Estimation via IP Weighting

Assume that treated and untreated are exchangeable conditional on \(L\), which, as in Chapters 12 to 15, contains sex, age, race, education, smoking intensity and duration, weight, recreational exercise, and physical activity in daily life. Assume positivity and consistency as well. IP weighted survival curves are then estimated in two steps.

Algorithm 1 (IP weighted survival curves) Step 1: Estimate each person’s stabilized IP weight \(SW^A\) as in Chapter 12. Fit a logistic model for \(\Pr[A = 1 \mid L]\). The denominator is the fitted \(\mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1 \mid L\right]\mathclose{}\) for the treated and \(1 - \mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1 \mid L\right]\mathclose{}\) for the untreated; the numerator is \(\mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1\right]\mathclose{}\) for the treated and \(1 - \mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1\right]\mathclose{}\) for the untreated, with \(\mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1\right]\mathclose{}\) estimated nonparametrically or by a logistic model without covariates. In the weighted pseudo-population \(L\) and \(A\) are independent, so \(L\) no longer confounds.

Step 2: In the person-time data, fit the hazards model

\[\operatorname{logit}\,\Pr[D^a_{k+1} = 1 \mid D^a_k = 0] = \beta_{0,k} + \beta_1 a + \beta_2 (a \times k) + \beta_3 (a \times k^2)\]

weighting each person by the estimated \(SW^A\). The weighted fit estimates the parameters of a marginal structural logistic model: the hazards over time that would have occurred had the whole study population been treated (\(a = 1\)), and had it been untreated (\(a = 0\)). (The book writes this model for the complementary probability \(\Pr[D^a_{k+1} = 0 \mid D^a_k = 0]\); because \(\operatorname{logit}(1 - p) = -\operatorname{logit}(p)\), the two forms are equivalent with the signs of the \(\beta\)’s reversed.)

The fitted counterfactual hazards \(\Pr[D^a_{k+1} = 1 \mid D^a_k = 0]\) are turned into survival estimates \(\Pr[D^a_{k+1} = 0] = \prod_{m=0}^k (1 - \Pr[D^a_{m+1} = 1 \mid D^a_m = 0])\) for \(a = 1\) and \(a = 0\) (the book’s Figure 17.6).

In the NHEFS, the estimated weights averaged 1, as they should, and ranged from 0.33 to 4.21.

Example 6 (NHEFS: IP weighted survival curves) Estimated 120-month survival: 80.7% if everyone had quit and 80.5% if no one had quit, a difference of 0.2% (95% confidence interval \(-4.1\%\) to \(3.7\%\), from 500 bootstrap samples). The curve under quitting lay below the one under not quitting for most of follow-up, but the difference was never larger in magnitude than 1.4 percentage points (the most extreme value was \(-1.4\%\); 95% confidence interval \(-3.4\%\) to \(0.7\%\)). After IP weighting for \(L\), then, the data show little evidence that quitting affects mortality at any point during follow-up.

WarningTwo Models Must Be Correct

These IP weighted estimates are valid only if the model for treatment and the marginal hazards model are both correctly specified.

5 17.5 The Parametric G-Formula (pp. 236-237)


IP weighting is one way to estimate the marginal survival curves under exchangeability, positivity, and consistency. Another is standardization with parametric models, the parametric g-formula.

5.1 Identification via Standardization

The counterfactual survival \(\Pr[D^a_{k+1} = 0]\) is an average of the survivals at \(k + 1\) conditional on \(L\) and on \(A = a\), weighted by the distribution of \(L\).

Theorem 1 (Identification of counterfactual survival by standardization) Let \(L\) be discrete. Under conditional exchangeability given \(L\), positivity, and consistency, the counterfactual survival equals the standardized survival

\[\Pr[D^a_{k+1} = 0] = \sum_l \Pr[D_{k+1} = 0 \mid L = l, A = a]\, \Pr[L = l]\]

Proof. The argument is the one that identifies the standardization formula in chapter 2, applied to survival. Step by step:

\[ \begin{aligned} \Pr[D^a_{k+1} = 0] &= \sum_l \Pr[D^a_{k+1} = 0 \mid L = l]\, \Pr[L = l] && \text{(law of total probability)}\\ &= \sum_l \Pr[D^a_{k+1} = 0 \mid L = l, A = a]\, \Pr[L = l] && \text{(conditional exchangeability; positivity)}\\ &= \sum_l \Pr[D_{k+1} = 0 \mid L = l, A = a]\, \Pr[L = l] && \text{(consistency)} \end{aligned} \]

5.2 Estimation

Algorithm 2 (Parametric g-formula for survival curves) Estimation again takes two steps.

Step 1: Estimate the conditional survivals \(\Pr[D_{k+1} = 0 \mid L = l, A = a]\) from the administratively censored data with a parametric hazards model like that of Section 17.2, now including the covariates \(L\), for example:

\[\operatorname{logit}\,\Pr[D_{k+1} = 1 \mid D_k = 0, L, A] = \theta_{0,k} + \theta_1 A + \theta_2 A \times k + \theta_3 A \times k^2 + \boldsymbol{\theta}_L^\top L\]

If this model is correct, it estimates the hazards \(\Pr[D_{k+1} = 1 \mid D_k = 0, L, A]\) for every combination of treatment and covariates. The product of one minus these hazards,

\[\prod_{m=0}^k \Pr[D_{m+1} = 0 \mid D_m = 0, L = l, A = a],\]

is the conditional survival \(\Pr[D_{k+1} = 0 \mid L = l, A = a]\). Under conditional exchangeability given \(L\), it is also the survival had everyone with \(L = l\) received treatment \(a\):

\[\Pr[D_{k+1} = 0 \mid L = l, A = a] = \Pr[D^a_{k+1} = 0 \mid L = l]\]

The model can thus produce, say, survival curves with and without treatment for 61-year-old white men with a college education who exercise little. We want marginal curves, though, not conditional ones.

Step 2: Average the conditional survivals over the distribution of \(L\), that is, standardize them, using the procedure of Section 13.3 (expand the dataset, fit the outcome model, predict, and average). The procedure also handles continuous components of \(L\), for which the sum over \(l\) becomes an integral. The book’s Figure 17.7 shows the resulting curves.

Chapter 12 called a model that conditions on all of \(L\) a faux marginal structural model.

Example 7 (NHEFS: Parametric g-formula survival curves) The estimated curve under quitting lay below the one under not quitting throughout follow-up, with a difference never larger in magnitude than 2.0 percentage points (the most extreme value was \(-2.0\%\); 95% confidence interval \(-5.6\%\) to \(1.8\%\)). Estimated 120-month survival was 80.4% under quitting and 80.6% under not quitting; the book reports the difference as 0.2% (95% confidence interval from \(-4.6\%\) to \(4.1\%\)), i.e., a magnitude of 0.2 percentage points (\(80.4\% - 80.6\% = -0.2\%\)). As with IP weighting, standardization for \(L\) yields little evidence of an effect of quitting on mortality during follow-up.

TipWhy the Two Sets of Curves Differ

The IP weighted and g-formula curves are close but not identical, because they depend on different modeling assumptions: IP weighting needs correct models for treatment and for the unconditional hazards, while the parametric g-formula needs a correct model for the conditional hazards.

6 17.6 G-Estimation of Structural Nested Models (pp. 237-240)


Remark 3 (Structural nested models do not give survivals). The methods so far contrast survivals or risks under different values of \(A\), computing survival from hazards fitted by logistic regression. That route works for IP weighting of marginal structural models and for the parametric g-formula, but not for g-estimation of structural nested models.

Structural nested models (Chapter 14) describe conditional causal contrasts, such as a difference or ratio of covariate-specific counterfactual means, rather than the separate quantities being contrasted, so they cannot give us survivals or hazards.

They cannot even approximate a hazard ratio, since structural nested logistic models are hard to extend to treatments that vary over time (Technical Point 14.1).

6.1 Structural Nested Models for Risks, Survivals, and Survival Times

Remark 4 (Structural nested models for survival outcomes). What a structural nested model can do is model ratios:

  • a structural nested log-linear model for the ratio of risks (cumulative incidences) under different treatments, the structural nested cumulative failure time (CFT) model (Technical Point 17.2), suited to rare failures because a log-linear model does not keep risks below 1;
  • for common failures, the analogous model for the ratio of survivals, the structural nested cumulative survival time (CST) model, suited, for the same reason, to settings where survival is rare;
  • more generally, a model for the ratio of survival times under different treatments, the accelerated failure time (AFT) model.
NoteTechnical Point 17.2: Structural Nested CFT and CST Models

With a treatment fixed at baseline, a (non-nested) structural CFT model relates the counterfactual risk under treatment \(a\) to that under treatment 0, given \(A\) and \(L\):

\[\frac{\Pr[D^a_k = 1 \mid L, A]}{\Pr[D^{a=0}_k = 1 \mid L, A]} = \operatorname{exp}\mathopen{}\left\{\gamma_k(L, A; \psi)\right\}\mathclose{}\]

with \(\gamma_k(L, A; \psi)\) a function of treatment and covariates indexed by a parameter \(\psi\), possibly a vector. \(\operatorname{exp}\mathopen{}\left\{\gamma_k(L, A; \psi)\right\}\mathclose{}\) has to equal 1 when \(A = 0\), since both sides then refer to the same treatment, and when treatment does not affect the outcome at \(k\). For instance, with \(\gamma_k(L, A; \psi) = \psi A\), \(\psi = 0\) means no effect, \(\psi < 0\) benefit, and \(\psi > 0\) harm.

A structural CST model is the same with survivals in place of risks:

\[\frac{\Pr[D^a_k = 0 \mid L, A]}{\Pr[D^{a=0}_k = 0 \mid L, A]} = \operatorname{exp}\mathopen{}\left\{\gamma_k(L, A; \psi)\right\}\mathclose{}\]

The two differ only in whether risk or survival is modeled multiplicatively, yet \(\gamma_k\) has a different meaning in each, because a multiplicative model for risk is not multiplicative for survival unless the effect is null. In continuous time, structural CST models and structural additive hazards models coincide (Tchetgen Tchetgen et al., 2015).

Structural CFT models need a specific rare-failure assumption at every value of \(L\). When it holds, they have an edge over AFT models: their unbiased estimating equations are differentiable in the parameters and therefore easy to solve (Page 2005; Picciotto et al. 2012). With time-varying treatments, they fall within the class of multivariate structural nested mean models (Robins 1994).

6.2 Structural Nested AFT Models

Definition 12 (Structural nested AFT model) Let \(T^a_i\) be individual \(i\)’s counterfactual survival time under treatment \(a\). The ratio \(T^{a=1}_i / T^{a=0}_i\) measures the effect of \(A\) on \(i\)’s survival time: above 1, treatment lengthens survival (benefit); below 1, it shortens survival (harm); equal to 1, no effect. Assume for now that the effect is identical for everyone. This gives the structural nested AFT model

\[T^a_i / T^{a=0}_i = \operatorname{exp}\mathopen{}\left\{-\psi_1 a\right\}\mathclose{}\]

in which \(\operatorname{exp}\mathopen{}\left\{-\psi_1\right\}\mathclose{}\) is the factor by which treatment stretches or shrinks each person’s survival time:

  • \(\psi_1 < 0\): treatment increases survival time
  • \(\psi_1 > 0\): treatment decreases survival time
  • \(\psi_1 = 0\): treatment does not affect survival time

The minus sign keeps the usual convention that positive parameters mean harm and negative ones benefit. The “nested” part of the name only matters for time-varying treatments (Chapter 21).

If the effect varies with covariates \(L\), a more general structural AFT model is \(T^a_i / T^{a=0}_i = \operatorname{exp}\mathopen{}\left\{-\psi_1 a - \psi_2 a L_i\right\}\mathclose{}\), with \(\psi_1\) and the vector \(\psi_2\) common to all individuals.

Proposition 3 (The AFT model in terms of observed survival times) Under consistency, the structural AFT model \(T^a_i / T^{a=0}_i = \operatorname{exp}\mathopen{}\left\{-\psi_1 a - \psi_2 a L_i\right\}\mathclose{}\) (a generalization of Definition 12) implies a model in terms of the observed survival time:

\[T^{a=0}_i = T_i \operatorname{exp}\mathopen{}\left\{\psi_1 A_i + \psi_2 A_i L_i\right\}\mathclose{}\]

Proof. Multiplying both sides by \(T^{a=0}_i \operatorname{exp}\mathopen{}\left\{\psi_1 a + \psi_2 a L_i\right\}\mathclose{}\) gives the equivalent form

\[T^{a=0}_i = T^a_i \operatorname{exp}\mathopen{}\left\{\psi_1 a + \psi_2 a L_i\right\}\mathclose{} \quad \text{for all individuals } i\]

By consistency, \(T^{A_i}_i = T_i\), so setting \(a = A_i\) gives the result.

The parameters of this AFT model can be estimated by g-estimation, modified to allow for administrative censoring.

WarningDeterministic, Rank-Preserving Models Are Implausible

As a description of reality the structural AFT model is implausible on two counts. It is deterministic: it says that each person’s untreated survival time \(T^{a=0}\) follows exactly from \(T\), \(A\), and \(L\). It is rank-preserving: if \(i\) would die before \(j\) when both are untreated (\(T^{a=0}_i < T^{a=0}_j\)), then \(i\) would also die before \(j\) when both are treated (\(T^{a=1}_i < T^{a=1}_j\)). Rank preservation is not believable, and methods that require it are generally best avoided (Chapter 14). The book uses a rank-preserving model anyway because it makes g-estimation easier to explain, and the procedure is identical with or without rank preservation. Robins (1997b) describes structural nested AFT models that are neither deterministic nor rank-preserving.

6.3 G-Estimation Without Censoring

Take the simpler rank-preserving model \(T^{a=0}_i = T_i \operatorname{exp}\mathopen{}\left\{\psi A_i\right\}\mathclose{}\). If no one were administratively censored, so that \(T\) were known for all, g-estimation would follow Section 14.5:

Algorithm 3 (G-estimation of a structural AFT model without censoring)  

  1. For many candidate values \(\psi^\dagger\), compute \(H_i(\psi^\dagger) = T_i \operatorname{exp}\mathopen{}\left\{\psi^\dagger A_i\right\}\mathclose{}\).
  2. The g-estimate of \(\psi\) is the \(\psi^\dagger\) at which \(H_i(\psi^\dagger)\) is unrelated to \(A\) in a logistic model for \(\Pr[A = 1]\) that includes \(H_i(\psi^\dagger)\) and the confounders \(L\).
TipUse a Directed Search

Statistical software offers directed search methods, such as the Nelder-Mead simplex, that need far less computation than a grid search.

6.4 Selection Bias from Administrative Censoring

With administrative censoring at time \(K\), \(H_i(\psi^\dagger)\) is unavailable for anyone whose \(T_i\) is unknown, and limiting g-estimation to those with observed survival times (\(T_i \leq K\)) introduces selection bias. Consider a 60-month randomized experiment with just three types of people, whose survival times with and without treatment appear in Table 1 (the book’s Table 17.1). Type 3 has the best prognosis and type 1 the worst. Randomization makes the expected share of each type the same in both arms.

Table 1: Survival times (months) under treatment and no treatment for each individual type in a 60-month experiment.
Type \(T^{a=0}\) \(T^{a=1}\)
1 36 24
2 72 48
3 108 72

In this example treatment shortens survival time by the same factor for every type: \(24/36 = 48/72 = 72/108 = 2/3\). In terms of the AFT model \(T^{a=1}/T^{a=0} = \operatorname{exp}\mathopen{}\left\{-\psi_1\right\}\mathclose{}\), this means \(\operatorname{exp}\mathopen{}\left\{-\psi_1\right\}\mathclose{} = 2/3\), so \(\psi_1 = \log(3/2) \approx 0.41 > 0\), a harmful effect.

Example 8 (Selection bias from administrative censoring) Follow-up ends at \(K = 60\) months, so by Table 1:

  • type 1 deaths are observed in either arm (36 and 24 are both below 60);
  • type 3 deaths are censored in either arm (108 and 72 are both above 60);
  • type 2 deaths are observed only in the treated arm (48 is below 60, 72 is not).

Among people with observed death times, then, the arms differ in composition: the treated include types 1 and 2, the untreated only type 1. Exchangeability fails, the untreated look worse off on average, and treatment appears more beneficial than it is. This selection bias (Chapter 8) occurs whenever treatment changes survival time at all.

6.5 Artificial Censoring

Definition 13 (Artificial censoring) To prevent the bias, the analysis must keep only people whose death would have been observed by \(K\) under either treatment, those with \(T^{a=0}_i \leq K\) and \(T^{a=1}_i \leq K\); in the example, every type 2 person is dropped. So besides the administratively censored (\(T_i > K\)), some people with observed deaths (\(T_i \leq K\)) are removed because under the other treatment their death would have come after \(K\). Removing such uncensored people is known as artificial censoring.

Let \(\Delta(\psi)\) indicate inclusion (1) or exclusion (0); g-estimation then uses \(\Delta(\psi^\dagger)\) in place of \(H(\psi^\dagger)\) (Technical Point 17.3).

NoteTechnical Point 17.3: Artificial Censoring

Define \(K(\psi)\) as the smallest untreated survival time compatible with a person who actually died at time \(K\). For a dichotomous treatment, \(K(\psi) = \inf\{K \operatorname{exp}\mathopen{}\left\{\psi A\right\}\mathclose{}\}\), the infimum taken over the treatment values \(A \in \{0, 1\}\). Evaluating the two candidates \(K\operatorname{exp}\mathopen{}\left\{\psi \times 0\right\}\mathclose{} = K\) and \(K \operatorname{exp}\mathopen{}\left\{\psi \times 1\right\}\mathclose{} = K\operatorname{exp}\mathopen{}\left\{\psi\right\}\mathclose{}\):

  • if treatment contracts survival time (\(\psi > 0\)), \(K\operatorname{exp}\mathopen{}\left\{\psi\right\}\mathclose{} > K\), so \(K(\psi) = K\);
  • if treatment expands survival time (\(\psi < 0\)), \(K\operatorname{exp}\mathopen{}\left\{\psi\right\}\mathclose{} < K\), so \(K(\psi) = K \operatorname{exp}\mathopen{}\left\{\psi\right\}\mathclose{}\);
  • if treatment has no effect (\(\psi = 0\)), both candidates equal \(K\), so \(K(\psi) = K\).

Every administratively censored person (\(T > K\)) gets \(\Delta(\psi) = 0\), since their survival time exceeds \(K\) under the treatment they actually received, so \(H(\psi) \geq K(\psi)\). Some people with observed deaths (\(T \leq K\)) also get \(\Delta(\psi) = 0\); they are the artificially censored, removed to prevent selection bias.

Because \(\Delta(\psi)\) depends only on \(H(\psi)\) and \(K\), and every function of \(H(\psi)\) at the true \(\psi\) is independent of \(A\) given \(L\) under conditional exchangeability, g-estimation can use \(\Delta(\psi)\) in place of \(H(\psi)\), that is, in Technical Point 14.2’s estimating equation. Hernán et al. (2005, Appendix) give practical details.

6.6 NHEFS Results and Practical Obstacles

Example 9 (NHEFS: G-estimation of a structural AFT model) Applied to the NHEFS with the rank-preserving model \(T^{a=0}_i = T_i \operatorname{exp}\mathopen{}\left\{\psi A_i\right\}\mathclose{}\), g-estimation gave \(\hat{\psi} = -0.047\) (95% confidence interval: \(-0.223\) to \(0.333\)). The quantity \(\operatorname{exp}\mathopen{}\left\{-\hat{\psi}\right\}\mathclose{} = \operatorname{exp}\mathopen{}\left\{0.047\right\}\mathclose{} \approx 1.05\) is the ratio of the median survival time had everyone received \(a = 1\) to the median survival time had everyone received \(a = 0\), so quitting appears to have little effect on time to death.

The point estimate of \(\psi\) is where the estimating function of Technical Point 17.3 reaches its minimum, and the 95% confidence limits are where it equals 3.84, the 0.95 quantile of a \(\chi^2\) distribution with one degree of freedom.

Survival analysis with instrumental variables, described by Tchetgen Tchetgen et al. (2015) and Robins (1997b), runs into problems like the ones described here for structural nested models.

Remark 5 (Why structural nested AFT models are rarely used). Structural nested models, AFT models included, are rarely used in practice (Chapter 14). Reasons include a shortage of user-friendly software and, for survival outcomes, the reliance on search algorithms that may not find a unique solution, which gets worse as parameters are added. Published applications therefore mostly use simple AFT models with no effect modification by covariates. That is a pity, since subject-matter knowledge, such as a biological mechanism, maps more naturally onto the parameters of a structural AFT model than onto those of a structural hazards model, particularly when the AFT model is neither deterministic nor rank-preserving.

7 Summary


Key concepts from this chapter:

  1. Survival analysis: Methods for time-to-event outcomes that accommodate administrative censoring
  2. Survival, risk, and hazard: Related but distinct measures; \(\Pr[D_k = 0] = \prod_{m=1}^k \Pr[D_m = 0 \mid D_{m-1} = 0]\)
  3. Person-time data format: Enables parametric modeling of time-varying hazards via logistic regression
  4. Why censoring matters: Staggered entry creates individual-specific censoring times; naive estimation is biased
  5. Counterfactual survivals: The target is \(\Pr[D^a_{k+1} = 0]\) under treatment \(a\) and no censoring
  6. IP weighting of marginal structural models: Estimates causal survival curves by weighting by \(SW^A\)
  7. Parametric g-formula: Standardizes conditional survivals to the covariate distribution
  8. G-estimation of structural nested AFT models: Models the ratio of survival times; requires artificial censoring to handle administrative censoring

Causal estimands for survival outcomes:

  • Causal survival difference: \(\Pr[D^{a=1}_{k} = 0] - \Pr[D^{a=0}_{k} = 0]\)
  • Causal risk difference: \(\Pr[D^{a=1}_{k} = 1] - \Pr[D^{a=0}_{k} = 1]\)

NHEFS smoking cessation example (120-month results):

  • Unadjusted: 76.2% survival among quitters vs. 82.0% among non-quitters (confounded)
  • IP weighting: 80.7% under cessation vs. 80.5% under continued smoking; difference 0.2%
  • Parametric g-formula: 80.4% under cessation vs. 80.6% under continued smoking; difference \(-0.2\%\) (reported in the book as 0.2%)
  • G-estimation of a structural AFT model: \(\hat{\psi} = -0.047\) (95% CI \(-0.223\) to \(0.333\)); survival time ratio \(\operatorname{exp}\mathopen{}\left\{-\hat{\psi}\right\}\mathclose{} \approx 1.05\)

None of the three adjusted analyses found much evidence that quitting smoking affects mortality.

Why not just use hazard ratios? (Fine Point 17.2):

  • Hazard ratios vary over time; a “single” hazard ratio from Cox regression is a weighted average
  • Even time-specific hazard ratios have built-in selection bias (conditioning on survival is conditioning on a collider)
  • Risk and survival differences/ratios are better-defined causal estimands
WarningCommon Mistakes in Survival Analysis
  1. Reporting only a single (time-averaged) hazard ratio from a Cox model and interpreting it causally
  2. Treating competing events as independent censoring (see Fine Point 17.1)
  3. Ignoring time-varying confounding (addressed in Part III)
  4. Restricting the analysis to uncensored individuals when estimating AFT models, which introduces selection bias

Looking ahead: Part III extends these ideas to time-varying treatments, where survival analysis becomes even more complex. The g-formula and IP weighting remain powerful tools, but they must account for the time-varying nature of both treatment and confounders.

8 References


Hernán, Miguel A, and James M Robins. 2020. Causal Inference: What If. Chapman & Hall/CRC. https://miguelhernan.org/whatifbook.
Back to top