Chapter 17: Causal Survival Analysis
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.
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
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
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)
2.1 The Event Indicator \(D_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:
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.
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.
2.3 Nonparametric and Parametric Estimation
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).
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.
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.
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).
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.
3.1 The Censoring Indicator
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.
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
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.
3.3 Correct Estimation Under As-If Randomly Assigned Censoring
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.
Estimation proceeds as before, except that the quantity estimated nonparametrically or with a logistic model is the cause-specific hazard (Definition 10).
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)
4.1 Counterfactual Survivals
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.
In the NHEFS, the estimated weights averaged 1, as they should, and ranged from 0.33 to 4.21.
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\).
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
Chapter 12 called a model that conditions on all of \(L\) a faux marginal structural model.
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)
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
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
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.
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.
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:
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.
| 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.
6.5 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).
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
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.
7 Summary
Key concepts from this chapter:
- Survival analysis: Methods for time-to-event outcomes that accommodate administrative censoring
- 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]\)
- Person-time data format: Enables parametric modeling of time-varying hazards via logistic regression
- Why censoring matters: Staggered entry creates individual-specific censoring times; naive estimation is biased
- Counterfactual survivals: The target is \(\Pr[D^a_{k+1} = 0]\) under treatment \(a\) and no censoring
- IP weighting of marginal structural models: Estimates causal survival curves by weighting by \(SW^A\)
- Parametric g-formula: Standardizes conditional survivals to the covariate distribution
- 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
- Reporting only a single (time-averaged) hazard ratio from a Cox model and interpreting it causally
- Treating competing events as independent censoring (see Fine Point 17.1)
- Ignoring time-varying confounding (addressed in Part III)
- 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.