---
title: "Chapter 21: G-Methods for Time-Varying Treatments"
format:
html: default
revealjs:
output-file: 21-g-methods-time-varying-slides.html
pdf:
output-file: 21-g-methods-time-varying-handout.pdf
docx:
output-file: 21-g-methods-time-varying.docx
preview-changed: true
---
{{< include ../latex-macros/macros.qmd >}}
Chapter 20 showed that, with treatment-confounder feedback,
traditional adjustment methods give a non-null estimate
even when a time-varying treatment has no effect at all.
This chapter shows that the **g-methods** --
the g-formula, IP weighting, g-estimation, and doubly robust versions that combine them --
recover the correct (null) answer in the same dataset.
Throughout, we compare **static** strategies under the identifiability conditions of Chapter 19:
sequential exchangeability, positivity, and consistency.
::: {.notes}
This chapter is based on @hernan2020causal [Chapter 21, pp. 277-303].
**Roadmap**: each method was introduced earlier for a time-fixed treatment:
IP weighting of marginal structural models (Chapter 12),
the g-formula (Chapter 13),
and g-estimation of structural nested models (Chapter 14)
(the book's introduction to Chapter 21 says Chapter 15 [@hernan2020causal, p. 277], but g-estimation is the subject of Chapter 14).
All three g-methods rest on the same identifying assumptions
but model different parts of the data distribution:
the g-formula models the outcome and covariates,
IP weighting models treatment,
and g-estimation models treatment plus a structural model for the effect of each "blip" of treatment.
Section 21.5 adds censoring, and Section 21.6 relates the g-formula to other identifying formulas.
:::
## 21.1 The G-Formula for Time-Varying Treatments (pp. 277-281)
---
The running example is the sequentially randomized experiment of Chapter 20
(Table 20.1, reproduced in the book as Table 21.1),
with treatments $A_0, A_1$, a confounder $L_1$ measured between them,
and an outcome $Y$.
There is no $L_0$ in this study.
::: {#tbl-seq-rand}
| $N$ | $A_0$ | $L_1$ | $A_1$ | Mean $Y$ |
|------:|:---:|:---:|:---:|---:|
| 2400 | 0 | 0 | 0 | 84 |
| 1600 | 0 | 0 | 1 | 84 |
| 2400 | 0 | 1 | 0 | 52 |
| 9600 | 0 | 1 | 1 | 52 |
| 4800 | 1 | 0 | 0 | 76 |
| 3200 | 1 | 0 | 1 | 76 |
| 1600 | 1 | 1 | 0 | 44 |
| 6400 | 1 | 1 | 1 | 44 |
The sequentially randomized experiment of Table 21.1 [@hernan2020causal, p. 277];
32,000 individuals in total.
:::
---
For the effect of the time-fixed treatment $A_1$ alone,
the g-formula is ordinary standardization:
$\E{Y^{a_1}} = \sum_{l_1} \E{Y \mid A_1 = a_1, L_1 = l_1} \Pr[L_1 = l_1]$.
For the joint treatment $(A_0, A_1)$, the weights become the distribution of $L_1$
**given past treatment**.
::: {#def-g-formula-tv}
## The G-Formula for Time-Varying Treatments
For two time points, the g-formula for the static strategy $(a_0, a_1)$ is
$$\sum_{l_1} \E{Y \mid A_0 = a_0, A_1 = a_1, L_1 = l_1} \, f(l_1 \mid a_0).$$
For $K+1$ time points $k = 0, \ldots, K$ and strategy $\bar{a}$, it is
$$\sum_{\bar{l}} \E{Y \mid \bar{A} = \bar{a}, \bar{L} = \bar{l}}
\prod_{k=0}^{K} f\paren{l_k \mid \bar{a}_{k-1}, \bar{l}_{k-1}},$$
where the sum is over all covariate histories $\bar{l}$.
Under sequential exchangeability for $Y^{\bar{a}}$ given $(\bar{L}_k, \bar{A}_{k-1})$ at each $k$,
positivity, and consistency, the g-formula equals $\E{Y^{\bar{a}}}$.
:::
::: {.notes}
Robins (1986, 1987) introduced the g-formula for time-varying treatments,
as cited in @hernan2020causal [p. 278].
Every factor is conditional on prior treatment and covariate history;
this conditioning is unnecessary in the time-fixed case,
where treatment and confounders are measured at a single time.
::: {.callout-warning title="Positivity Is Needed to Compute the G-Formula"}
**Positivity**: the g-formula is computable only if,
for every $l_1$ with $f(l_1 \mid a_0) \neq 0$,
some individuals have $(A_0 = a_0, A_1 = a_1, L_1 = l_1)$
(the definition of positivity in Technical Point 19.2).
:::
:::
::: {.callout-note title="Fine Point 21.1: Treatment and Covariate History"}
The relevant history at time $k$ is the set of treatments and confounders
needed for conditional exchangeability of $A_k$.
Usually this is the chronological past,
but confounders can in principle lie in the temporal future of treatment (Fine Point 7.4),
and adjusting for some past variables can create M-bias (Figure 7.4)
[@hernan2020causal, p. 280].
:::
---
::: {#exm-g-formula-tv}
## G-Formula in the Sequentially Randomized Experiment
**Never treat**, $(a_0, a_1) = (0, 0)$.
Among the 16,000 individuals with $A_0 = 0$,
$2400 + 1600 = 4000$ have $L_1 = 0$ and $2400 + 9600 = 12000$ have $L_1 = 1$, so
$f(L_1 = 0 \mid A_0 = 0) = 0.25$ and $f(L_1 = 1 \mid A_0 = 0) = 0.75$.
Then
$$
\begin{aligned}
\E{Y^{a_0=0, a_1=0}}
&= \E{Y \mid 0, 0, L_1 = 0} \times 0.25 + \E{Y \mid 0, 0, L_1 = 1} \times 0.75 \\
&= 84 \times 0.25 + 52 \times 0.75 \\
&= 21 + 39 = 60.
\end{aligned}
$$
**Always treat**, $(a_0, a_1) = (1, 1)$.
Among the 16,000 individuals with $A_0 = 1$,
$4800 + 3200 = 8000$ have $L_1 = 0$, so $f(L_1 = 0 \mid A_0 = 1) = 0.5$, and
$$
\begin{aligned}
\E{Y^{a_0=1, a_1=1}}
&= 76 \times 0.5 + 44 \times 0.5 \\
&= 38 + 22 = 60.
\end{aligned}
$$
The g-formula estimate of $\E{Y^{1,1}} - \E{Y^{0,0}}$ is $60 - 60 = 0$, the correct null value
[@hernan2020causal, p. 278].
:::
### The G-Formula as a Simulation
::: {#rem-g-formula-simulation}
## The g-formula simulates the counterfactual joint distribution
Under sequential exchangeability for $Y$ and $\bar{L}$ jointly,
the g-formula simulates the joint distribution of the counterfactuals $(Y^{\bar{a}}, \bar{L}^{\bar{a}})$
that would have been observed had everybody followed strategy $\bar{a}$.
On the tree graph (Figures 21.1-21.2), this means building a new tree in which
everyone receives $a_0$ and $a_1$ with probability 1,
while $\Pr[L_1 = l_1 \mid A_0 = a_0]$ and $\E{Y \mid A_0 = a_0, A_1 = a_1, L_1 = l_1}$
keep their observed values.
:::
::: {.notes}
::: {.callout-warning title="Two Cautions About the G-Formula"}
Two cautions [@hernan2020causal, p. 279]:
1. **The value of the g-formula depends on what is in $L$.**
If we wrongly omit $L_1$ (believing it is not a confounder),
the g-formula reduces to $\E{Y \mid A_0 = a_0, A_1 = a_1}$,
which has no causal interpretation when $L_1$ affects $A_1$ (Figure 20.8).
2. **The g-formula can be causal even when its components are not.**
Under Figure 20.9, where only static sequential exchangeability holds,
the g-formula including $L_1$ still identifies $\E{Y^{\bar{a}}}$,
yet neither $\Pr[L_1 = l_1 \mid A_0 = a_0]$ nor $\E{Y \mid A_0 = a_0, A_1 = a_1, L_1 = l_1}$
equals its counterfactual analogue.
In a sequentially randomized trial (Figures 20.1 and 20.2)
the components are causal as well:
$\Pr[L_1 = l_1 \mid A_0 = a_0] = \Pr[L_1^{a_0} = l_1]$ and
$\E{Y \mid A_0 = a_0, A_1 = a_1, L_1 = l_1} = \E{Y^{a_0, a_1} \mid L_1^{a_0} = l_1}$,
so the g-formula is $\sum_{l_1} \E{Y^{a_0, a_1} \mid L_1^{a_0} = l_1} \Pr[L_1^{a_0} = l_1] = \E{Y^{a_0, a_1}}$.
:::
:::
### Estimation: the Plug-In (Parametric) G-Formula
::: {#def-plug-in-g-formula}
## Plug-in and parametric g-formula
With many confounders or time points, the components must be estimated, for example
- a linear regression model for $\E{Y \mid \bar{A} = \bar{a}, \bar{L} = \bar{l}}$, and
- logistic regression models for the discrete confounders $L_k$, $k \neq 0$
(the distribution of $L_0$ can be estimated without a model, as in Section 13.3).
Plugging the estimates into the formula gives the **plug-in g-formula**;
when the estimates come from parametric models, the **parametric g-formula**.
:::
::: {.notes}
::: {#def-g-formula-random-strategy}
## The g-formula for a random strategy
**Random strategies**: under sequential exchangeability the g-formula also handles a random strategy $f^{int}$,
such as "at each $k$, independently treat with probability 0.3"
(so $f^{int}(1 \mid \bar{a}_{k-1}, \bar{l}_k) = 0.3$).
The general expression is
$$\sum_{\bar{a}, \bar{l}} \E{Y \mid \bar{A} = \bar{a}, \bar{L} = \bar{l}}
\prod_{k=0}^{K} f\paren{l_k \mid \bar{a}_{k-1}, \bar{l}_{k-1}}
\prod_{k=0}^{K} f^{int}\paren{a_k \mid \bar{a}_{k-1}, \bar{l}_k}.$$
Replacing $f^{int}$ by the observed $f(a_k \mid \bar{a}_{k-1}, \bar{l}_k)$ gives the observed mean of $Y$.
For a deterministic strategy, $f^{int}$ is 1 for the mandated treatment values and 0 otherwise,
so the $f^{int}$ factors and the sum over $\bar{a}$ can be dropped.
:::
:::
::: {.callout-note title="Technical Point 21.1: The g-formula Density"}
For a static strategy $\bar{a}$,
the g-formula density of $(Y, \bar{L})$ at $(y, \bar{l})$ is
$f(y \mid \bar{a}_K, \bar{l}_K) \prod_{k=0}^{K} f(l_k \mid \bar{a}_{k-1}, \bar{l}_{k-1})$,
and integrating out $\bar{l}$ gives the g-formula density of $Y$.
For a deterministic dynamic strategy $g$, replace each $a_k$
by the value $a_k^g$ that $g$ assigns given the simulated history.
:::
::: {.notes}
::: {.callout-tip title="Software for the Parametric G-Formula"}
**Software**: the `gfoRmula` R package (Lin et al. 2019) is on CRAN,
and a GFORMULA SAS macro is available on GitHub [@hernan2020causal, p. 281].
:::
:::
## 21.2 IP Weighting for Time-Varying Treatments (pp. 282-285)
---
For the time-fixed treatment $A_1$ alone (Chapter 12),
the weights are $W^{A_1} = 1/f(A_1 \mid L_1)$ or $SW^{A_1} = f(A_1)/f(A_1 \mid L_1)$.
For a time-varying treatment, the denominator becomes
each individual's probability of the **treatment history** they received,
given their treatment and covariate history.
::: {#def-ipw-tv}
## IP Weights for Time-Varying Treatments
With $A_{-1} \equiv 0$ by definition, the **nonstabilized** IP weights are
$$W^{\bar{A}} = \prod_{k=0}^{K} \frac{1}{f\paren{A_k \mid \bar{A}_{k-1}, \bar{L}_k}},$$
and the **stabilized** IP weights are
$$SW^{\bar{A}} = \prod_{k=0}^{K} \frac{f\paren{A_k \mid \bar{A}_{k-1}}}{f\paren{A_k \mid \bar{A}_{k-1}, \bar{L}_k}}.$$
For $K = 1$: $W^{\bar{A}} = \dfrac{1}{f(A_0 \mid L_0)} \times \dfrac{1}{f(A_1 \mid A_0, L_0, L_1)}$.
:::
::: {#prp-ipw-tv-identification}
## IP weighting identifies the counterfactual mean
Under the identifiability conditions (sequential exchangeability, positivity, and consistency),
$\E{Y^{a_0, a_1}}$ equals the mean
$\Expfc{ps}{Y \mid A_0 = a_0, A_1 = a_1}$
in the pseudo-population created by either set of weights in @def-ipw-tv.
:::
::: {.notes}
As in Technical Point 12.2, the pseudo-population mean equals
$\E{W^{\bar{A}} Y I(A_0 = a_0, A_1 = a_1)}$
(nonstabilized) or the ratio
$\E{SW^{\bar{A}} Y I(\cdot)} / \E{SW^{\bar{A}} I(\cdot)}$ (stabilized),
whether or not sequential exchangeability holds [@hernan2020causal, p. 282].
In the pseudo-population, the probability of treatment at each $k$ is
the constant $1/2$ (nonstabilized weights)
or depends at most on past treatment (stabilized weights),
so sequential *unconditional* exchangeability holds there,
and $\E{Y^{\bar{a}}} - \E{Y^{\bar{a}'}} = \Expfc{ps}{Y \mid \bar{A} = \bar{a}} - \Expfc{ps}{Y \mid \bar{A} = \bar{a}'}$.
In a true sequentially randomized trial the treatment probabilities are known by design,
so nonstabilized-weight estimates are unbiased;
in observational studies they must be estimated,
for example by a pooled logistic model for $\Pr[A_k = 1 \mid \bar{A}_{k-1}, \bar{L}_k]$
that includes functions of time $k$.
::: {.callout-warning title="The Denominator Model Must Be Correct"}
A misspecified denominator model biases the estimate;
a misspecified numerator model for $f(A_k \mid \bar{A}_{k-1})$ does not.
:::
:::
---
::: {#exm-ipw-tv}
## IP Weighting in the Sequentially Randomized Experiment
There is no $L_0$, so the denominator is $f(A_0) f(A_1 \mid A_0, L_1)$, with $f(A_0 = 0) = 16000/32000 = 0.5$.
For the never-treated rows of @tbl-seq-rand:
- $(A_0, L_1, A_1) = (0, 0, 0)$: $f(A_1 = 0 \mid A_0 = 0, L_1 = 0) = 2400/4000 = 0.6$,
so $W^{\bar{A}} = 1/(0.5 \times 0.6) = 10/3$
and the row contributes $2400 \times 10/3 = 8000$ pseudo-individuals.
- $(0, 1, 0)$: $f(A_1 = 0 \mid A_0 = 0, L_1 = 1) = 2400/12000 = 0.2$,
so $W^{\bar{A}} = 1/(0.5 \times 0.2) = 10$
and the row contributes $2400 \times 10 = 24000$.
The 32,000 pseudo-individuals with $(A_0, A_1) = (0, 0)$ give
$$
\begin{aligned}
\Expfc{ps}{Y \mid A_0 = 0, A_1 = 0}
&= 84 \times \frac{8000}{32000} + 52 \times \frac{24000}{32000} \\
&= 21 + 39 = 60.
\end{aligned}
$$
The same calculation for $(1, 1)$ also gives 60, so the IP weighted effect estimate is 0.
The whole pseudo-population has $4 \times 32000 = 128000$ individuals,
the 32,000 study individuals times the 4 static strategies
[@hernan2020causal, pp. 282-283].
:::
::: {.notes}
::: {#rem-g-formula-ipw-equal}
## Equal nonparametric estimates are not evidence of causality
The nonparametric g-formula and IP weighted estimates are exactly equal.
That equality has nothing to do with causality:
it holds even when the identifiability conditions fail
and neither estimate has a causal interpretation.
Stabilized weights give the same estimate of 0 here (check for yourself)
[@hernan2020causal, p. 283].
:::
::: {.callout-tip title="Estimate the Counterfactual Mean Both Ways"}
**Comparing the two parametric estimators**: if the parametric g-formula and the parametric IP weighted estimates
differ by more than sampling variability (quantified by bootstrapping the difference),
then at least one set of models is misspecified,
whether or not the identifiability assumptions hold.
Similar estimates do not prove the models are correct,
because both may be biased in the same direction.
The book recommends estimating $\E{Y^{\bar{a}}}$ both ways and,
if they differ substantially by a prespecified criterion, revisiting the models
[@hernan2020causal, p. 284].
:::
:::
### Marginal Structural Models
::: {#def-msm-tv}
## Marginal structural mean model for a time-varying treatment
There are far more strategies $\bar{a}$ than can be estimated one at a time,
so we pool information with a **marginal structural mean model**, for example
$$\E{Y^{\bar{a}}} = \beta_0 + \beta_1 \, \mathrm{cum}(\bar{a}), \qquad \mathrm{cum}(\bar{a}) = \sum_{k=0}^{K} a_k.$$
This model is unsaturated: one unknown per strategy on the left, two parameters on the right.
The average causal effect is $\E{Y^{\bar{a}}} - \E{Y^{\bar{a} = \bar{0}}} = \beta_1 \times \mathrm{cum}(\bar{a})$.
:::
---
::: {#prp-msm-tv-wls}
## IP weighted least squares is consistent for the marginal structural model
Suppose the marginal structural model of @def-msm-tv is correctly specified,
and the identifiability conditions (sequential exchangeability, positivity, and consistency) hold.
Let $\hb_1$ be the coefficient of $\mathrm{cum}(\bar{A})$ in the fit of the model
$\E{Y \mid \bar{A}} = \theta_0 + \theta_1 \, \mathrm{cum}(\bar{A})$
by weighted least squares with weights $SW^{\bar{A}}$ or $W^{\bar{A}}$ of @def-ipw-tv
(known, or estimated from a correctly specified treatment model).
Then $\hb_1$ is consistent for the causal $\beta_1$,
which in general differs from the associational $\theta_1$ of the unweighted data.
:::
::: {.notes}
::: {.callout-tip title="Variance Estimation and Choice of Weights"}
**Variance**: use the nonparametric bootstrap, the analytic variance,
or a conservative 95% interval from the robust variance estimator.
For a non-saturated model, intervals are typically narrower with $SW^{\bar{A}}$ than with $W^{\bar{A}}$,
so $SW^{\bar{A}}$ is preferred [@hernan2020causal, p. 285].
:::
::: {.callout-tip title="Checking the Marginal Structural Model"}
**Checking the MSM**: fitting
$$\E{Y \mid \bar{A}} = \theta_0 + \theta_1 \mathrm{cum}(\bar{A}) + \theta_2 \mathrm{cum}_{-5}(\bar{A}) + \theta_3 \mathrm{cum}(\bar{A})^2$$
with the same weights, where $\mathrm{cum}_{-5}$ is cumulative treatment in the final 5 months,
gives a 2-degree-of-freedom Wald test of $\theta_2 = \theta_3 = 0$,
a test that the original MSM is correctly specified.
In practice one may use other summaries of treatment history and flexible functions such as cubic splines.
:::
:::
::: {#def-msm-effect-modification-tv}
## Marginal structural model with a baseline effect modifier
**Effect modification**: for a dichotomous *baseline* variable $V$ in $L_0$,
a marginal structural model with effect modification by $V$ is
$\E{Y^{\bar{a}} \mid V} = \beta_0 + \beta_1 \mathrm{cum}(\bar{a}) + \beta_2 V + \beta_3 \mathrm{cum}(\bar{a}) V$.
Its parameters are estimated by fitting the associational model
$\E{Y \mid \bar{A}, V} = \theta_0 + \theta_1 \mathrm{cum}(\bar{A}) + \theta_2 V + \theta_3 \mathrm{cum}(\bar{A}) V$
by weighted least squares with $W^{\bar{A}}$ or, better,
$SW^{\bar{A}}(V) = \prod_{k=0}^{K} f(A_k \mid \bar{A}_{k-1}, V) / f(A_k \mid \bar{A}_{k-1}, \bar{L}_k)$.
:::
---
::: {#prp-msm-effect-modification-tv-wls}
## IP weighted least squares is consistent for the effect-modification model
Suppose the marginal structural model of @def-msm-effect-modification-tv is correctly specified,
the identifiability conditions (sequential exchangeability, positivity, and consistency) hold,
and $V$ is a baseline variable in $L_0$.
Let $\hb_1$ and $\hb_3$ be the coefficients of $\mathrm{cum}(\bar{A})$ and $\mathrm{cum}(\bar{A}) V$
in the fit of the associational model of @def-msm-effect-modification-tv
by weighted least squares with weights $SW^{\bar{A}}(V)$ of @def-msm-effect-modification-tv
or $W^{\bar{A}}$ of @def-ipw-tv
(known, or estimated from a correctly specified model for treatment in the denominator).
Then $\hb_1$ and $\hb_3$ are consistent for the causal $\beta_1$ and $\beta_3$.
:::
::: {.notes}
::: {.callout-warning title="Only Baseline Variables as Effect Modifiers"}
With treatment-confounder feedback, $V$ may include only baseline variables;
if $V$ included any $L_k$ with $k > 0$,
the weighted-fit coefficients of $\mathrm{cum}(\bar{A})$ and $\mathrm{cum}(\bar{A}) V$ could be non-null even though $\beta_1 = \beta_3 = 0$ under the null
[@hernan2020causal, p. 285].
:::
:::
::: {.callout-note title="Technical Point 21.2: IP Weighting for Dynamic Treatment Strategies"}
For a deterministic dynamic strategy $g$,
the g-formula equals the mean of $Y$ among pseudo-population members who follow $g$
in the pseudo-population created by **nonstabilized** weights.
Stabilized weights cannot be used, because their numerator depends on $A$.
:::
::: {.callout-note title="Technical Point 21.3: The g-null Paradox"}
Under the sharp null with treatment-confounder feedback,
the g-formula does not depend on $\bar{a}$,
but non-saturated parametric models for $\E{Y \mid \bar{A}, \bar{L}}$ and $f(l_k \mid \cdot)$
with variation-independent parameters cannot all be correctly specified
when $L_k$ has a discrete component (Robins and Wasserman 1997).
The estimated g-formula may therefore falsely reject the null even in a sequentially randomized experiment.
In practice the bias seems small relative to random variability.
MSMs and structural nested mean models do not suffer from this paradox:
both are correctly specified under the null whatever functional form we choose
[@hernan2020causal, p. 286].
:::
## 21.3 A Doubly Robust Estimator for Time-Varying Treatments (pp. 286-288)
---
::: {#def-doubly-robust-tv}
## Doubly robust estimator
IP weighting needs a correct model for treatment given confounders;
the g-formula needs a correct model for the outcome given treatment and confounders.
A **doubly robust** estimator is consistent if *either* model is correct,
without our knowing which,
giving "two chances to get it right" [@hernan2020causal, p. 286].
:::
### Review: a Doubly Robust Plug-In Estimator for a Time-Fixed Treatment
::: {#alg-dr-time-fixed}
## Doubly robust plug-in estimator for a time-fixed treatment
For a binary $A$, binary $Y$, and many confounders $L$, estimate $\E{Y^a}$ in three steps:
1. Fit a treatment model and compute $\hat{f}(a \mid L) = \hPr{A = a \mid L}$.
2. Among individuals with $A = a$, fit by maximum likelihood an outcome model
that includes $\hat{W}^a = 1/\hat{f}(a \mid L)$ as a covariate, such as
$b(a, L; \theta) = \expit(\theta_{a,0} + \theta_{a,1} L + \theta_{a,2} \hat{W}^a)$.
3. Average the predictions $b(a, L; \hth)$ over **all** individuals, treated and untreated.
The difference of the step-3 averages for $a = 1$ and $a = 0$
is a doubly robust estimator of $\E{Y^{a=1}} - \E{Y^{a=0}}$.
:::
### Extension to Time-Varying Treatments
::: {#alg-dr-time-varying}
## Doubly robust plug-in estimator for a time-varying treatment
To estimate $\E{Y^{\bar{a} = \bar{1}}}$ ("always treated"):
1. **Treatment model**: fit a model $\pi_k(\bar{L}_k; \alpha)$ for
$\Pr[A_k = 1 \mid \bar{A}_{k-1} = \bar{1}_{k-1}, \bar{L}_k]$, pooled over persons and times,
using only individuals treated through $k - 1$.
For those treated through $m$, compute the time-varying weights
$\hat{W}^{\bar{1}_m} = \prod_{k=0}^{m} 1/\hat{\pi}_k$.
2. **Sequential outcome models**: for $m = K, K-1, \ldots, 0$,
among individuals treated through $m$, fit a model $b_m(\bar{L}_m; \beta_m)$
that includes $\hat{W}^{\bar{1}_m}$ as a covariate.
The time-$K$ model has dependent variable $Y$;
the time-$m$ model, $m < K$, has as dependent variable the predictions
$\hat{B}_{m+1} = b_{m+1}(\bar{L}_{m+1}; \hb_{m+1})$ from the time-$(m+1)$ model.
3. **Average**: estimate $\E{Y^{\bar{a} = \bar{1}}}$ by the sample average of $\hat{B}_0$ over all individuals.
Repeat with $\bar{a} = \bar{0}$ for "never treated";
the difference estimates the average causal effect.
:::
::: {.notes}
::: {.callout-tip title="Choosing the Sequential Outcome Models"}
For a binary $Y$, the time-$m$ model can be logistic,
$\expit(\gamma_m X_m + \varsigma_m \hat{W}^{\bar{1}_m})$ with $X_m$ a function of $\bar{L}_m$:
although $\hat{B}_{m+1}$ is not 0 or 1, it lies in $[0, 1]$.
For a continuous $Y$, use a linear model [@hernan2020causal, p. 288].
:::
This estimator is due to Bang and Robins (2005), building on Robins (2000).
It is a targeted minimum loss-based estimator (TMLE) in the terminology of van der Laan and coauthors
[@hernan2020causal, p. 287].
:::
---
::: {#prp-dr-k-plus-2-robust}
## The time-varying doubly robust estimator is $K + 2$ robust
Suppose each sequential outcome model $b_m$ of @alg-dr-time-varying
is fit by solving the score equations of a generalized linear model with a canonical link
(logistic for a binary $Y$, linear for a continuous $Y$).
For a binary $Y$ and $m < K$, the response $\hat{B}_{m+1}$ lies in $[0, 1]$ rather than $\{0, 1\}$,
so the logistic fit is a quasi-likelihood fit.
Suppose also that each $b_m$ includes the weight $\hat{W}^{\bar{1}_m}$ as a covariate.
Suppose further that positivity holds for the strategy $\bar{a} = \bar{1}$.
Then the estimator of @alg-dr-time-varying is consistent for the g-formula functional of @def-g-formula-tv
at $\bar{a} = \bar{1}$ if
(i) the outcome models $b_m$ are correct for all $m$, or
(ii) the treatment models $\pi_k$ are correct for all $k$.
In fact it is **$K + 2$ robust**: it is also consistent if,
for some $m \in \{0, \ldots, K - 1\}$,
the treatment model is correct for times $0$ to $m$
and the outcome model is correct for times $m + 1$ to $K$
[@hernan2020causal, p. 288].
If, in addition, sequential exchangeability and consistency hold,
that functional equals $\E{Y^{\bar{a} = \bar{1}}}$,
so the estimator is consistent for the counterfactual mean.
:::
::: {.callout-note title="Technical Point 21.4: A K+2 Robust Augmented IP Weighted Estimator"}
For $K = 1$, let $\hat{b}_1(L_0, L_1) = \hE{Y \mid L_0, A_0 = 1, L_1, A_1 = 1}$
and $\hat{b}_0(L_0) = \hE{\hat{b}_1(L_0, L_1) \mid A_0 = 1, L_0}$
(the iterated conditional expectation, ICE, plug-in estimator of the g-formula averages $\hat{b}_0$),
and let $\hat{\pi}_0, \hat{\pi}_1$ estimate $\Pr[A_0 = 1 \mid L_0]$ and $\Pr[A_1 = 1 \mid L_0, L_1, A_0 = 1]$.
The augmented IP weighted estimator of Robins et al. (1994) is the sample average of
$$\hat{U}_{TR} = \hat{b}_0(L_0)
+ \frac{A_0 A_1}{\hat{\pi}_0 \hat{\pi}_1}\paren{Y - \hat{b}_1(L_0, L_1)}
+ \frac{A_0}{\hat{\pi}_0}\paren{\hat{b}_1(L_0, L_1) - \hat{b}_0(L_0)}.$$
It is consistent if (a) $\hat{\pi}_0, \hat{\pi}_1$ are consistent
(the augmentation terms then average to 0, leaving the IP weighted estimator),
(b) both outcome regressions are consistent (the last two terms then average to 0, leaving the ICE estimator), or
(c) $\hat{b}_1$ and $\hat{\pi}_0$ are consistent (Molina et al. 2017):
hence "triply", that is $K + 2$, robust.
A modification due to Tchetgen Tchetgen (2009) is $2^{K+1}$ (quadruply) robust.
Technical Point 21.5 shows that including the terms $A_0 A_1/(\hat{\pi}_0 \hat{\pi}_1)$ and $A_0/\hat{\pi}_0$
as covariates in the sequential logistic outcome models makes both correction terms average to exactly 0,
so the plug-in estimate stays in $[0, 1]$;
the estimator in the main text is an instance of this construction.
:::
::: {.notes}
::: {.callout-warning title="Fact-Check: A Term in the Book's First Display of the Estimator"}
*Fact-check note*: the book's first display of $\hat{U}_{TR}$ [@hernan2020causal, p. 288] writes the $\hat{b}_1$ term as
$-\{A_0 A_1/(\hat{\pi}_0 \hat{\pi}_1) - 1\}\hat{b}_1$,
which disagrees with the rearranged form the same Technical Point gives (the one shown here).
Expanding the rearranged form term by term gives
$$\hat{b}_0 + \frac{A_0 A_1}{\hat{\pi}_0 \hat{\pi}_1}\paren{Y - \hat{b}_1} + \frac{A_0}{\hat{\pi}_0}\paren{\hat{b}_1 - \hat{b}_0}
= \frac{A_0 A_1 Y}{\hat{\pi}_0 \hat{\pi}_1}
- \cb{\frac{A_0 A_1}{\hat{\pi}_0 \hat{\pi}_1} - \frac{A_0}{\hat{\pi}_0}}\hat{b}_1
- \cb{\frac{A_0}{\hat{\pi}_0} - 1}\hat{b}_0,$$
so the "$-1$" inside the braces on $\hat{b}_1$ should be $-A_0/\hat{\pi}_0$.
The two versions differ by $(1 - A_0/\hat{\pi}_0)\,\hat{b}_1$,
which does not average to 0 in general,
so only the corrected form supports the three robustness arguments.
The book's third rearrangement (used for the Molina et al. argument) agrees with the corrected form.
:::
:::
::: {.callout-note title="Fine Point 21.2: Representations of the g-formula"}
The same g-formula can be written as
standardization (the *joint density modeling* estimator of Section 21.1),
as iterated conditional expectations (the ICE estimator, the basis of this section),
or as an IP weighted mean (Section 21.2).
Without parametric restrictions the three coincide, but in practice each suggests a different estimator
[@hernan2020causal, p. 291].
:::
::: {.callout-note title="Technical Point 21.6: A Multiply Robust Estimator"}
This technical point gives a multiply robust plug-in estimator of $\E{Y^g}$
for static, dynamic, or random strategies $g$, based on Rotnitzky et al. (2017);
without the weight covariate, the same algorithm is the ICE estimator of the g-formula.
The book expects such estimators, fit with machine learning and sample splitting,
to become more common as software matures, especially for failure-time outcomes.
:::
## 21.4 G-Estimation for Time-Varying Treatments (pp. 289-296)
---
::: {#def-snmm-tv}
## Structural Nested Mean Model
A **structural nested mean model** (SNMM) for a time-varying treatment
specifies, for each $k = 0, \ldots, K$,
$$\E{Y^{\bar{a}_{k-1}, a_k, \underline{0}_{k+1}} - Y^{\bar{a}_{k-1}, \underline{0}_k}
\mid \bar{L}_k^{\bar{a}_{k-1}} = \bar{l}_k, \bar{A}_{k-1} = \bar{a}_{k-1}, A_k = a_k}
= a_k \, \gamma_k\paren{\bar{a}_{k-1}, \bar{l}_k, \beta},$$
where $(\bar{a}_{k-1}, a_k, \underline{0}_{k+1})$ follows $\bar{a}_{k-1}$ through $k - 1$, gives $a_k$ at $k$,
and gives no treatment from $k + 1$ to $K$;
$\gamma_k(\bar{a}_{k-1}, \bar{l}_k, \psi^\dagger)$ is a known function of past treatment, covariate history, and a parameter value $\psi^\dagger$,
with $\gamma_k(\bar{a}_{k-1}, \bar{l}_k, 0) = 0$;
and $\beta$ is the unknown true value of that parameter.
:::
::: {#rem-snmm-last-blip}
## What an SNMM models
An SNMM models the effect on the mean of $Y$ of a **last blip** of treatment $a_k$ at time $k$,
as a function of past treatment and covariate history.
Because $\gamma_k(\cdot, 0) = 0$, the null hypothesis of no effect corresponds to $\beta = 0$.
:::
---
::: {#exm-snmm-saturated-two-times}
## A saturated SNMM for two time points
An SNMM (@def-snmm-tv) has **one equation per time point**.
For @tbl-seq-rand ($K = 1$), a saturated additive SNMM is
$$
\begin{aligned}
k = 0: &\quad \E{Y^{a_0, a_1 = 0} - Y^{a_0 = 0, a_1 = 0} \mid A_0 = a_0} = \beta_0 a_0, \\
k = 1: &\quad \E{Y^{a_0, a_1} - Y^{a_0, a_1 = 0} \mid L_1^{a_0} = l_1, A_0 = a_0, A_1^{a_0} = a_1}
= a_1 \paren{\beta_{11} + \beta_{12} l_1 + \beta_{13} a_0 + \beta_{14} a_0 l_1}.
\end{aligned}
$$
The first equation is the effect of treating at time 0 and never again;
the second, of treating at time 1 and never again.
Under sequential exchangeability, $A_0 = a_0$ and $A_1^{a_0} = a_1$ can be dropped from the conditioning events.
:::
::: {.notes}
Both equations are saturated:
$\beta_0$ covers the single possible history at time 0,
and the four $\beta_{1j}$ cover the four histories $(A_0, L_1)$ at time 1.
The effect of $a_1$ is
- $\beta_{11}$ when $A_0 = 0, L_1^{a_0} = 0$;
- $\beta_{11} + \beta_{12}$ when $A_0 = 0, L_1^{a_0} = 1$;
- $\beta_{11} + \beta_{13}$ when $A_0 = 1, L_1^{a_0} = 0$;
- $\beta_{11} + \beta_{12} + \beta_{13} + \beta_{14}$ when $A_0 = 1, L_1^{a_0} = 1$.
The two equations nested within each other are why the model is called *nested*
[@hernan2020causal, p. 291].
:::
---
::: {#def-rank-preserving-snmm-tv}
## Additive rank-preserving model for two time points
As in Chapter 14, g-estimation starts from the additive **rank-preserving** model for each individual $i$:
$$
\begin{aligned}
Y_i^{a_0, 0} &= Y_i^{0, 0} + \psi_0 a_0, \\
Y_i^{a_0, a_1} &= Y_i^{a_0, 0} + \psi_{11} a_1 + \psi_{12} a_1 L_{1,i}^{a_0} + \psi_{13} a_1 a_0 + \psi_{14} a_1 a_0 L_{1,i}^{a_0}.
\end{aligned}
$$
The second equation is only *locally* rank-preserving:
ranks are preserved among individuals with the same $a_0$ and $l_1$.
:::
::: {.callout-warning title="Rank Preservation Is Implausible but Not Required"}
Rank preservation is biologically implausible,
but it is not needed.
G-estimation (the Chapter 14 method, extended here to time-varying treatment) is consistent for the SNMM parameters $\beta$ of @def-snmm-tv
even if the rank-preserving model is wrong,
provided the SNMM is correctly specified,
sequential exchangeability for $Y$ and consistency hold,
and the model for treatment given past treatment and covariates is correctly specified.
:::
::: {.notes}
The proof is in Robins (1994).
To fit an unsaturated SNMM by g-estimation, positivity is not required
[@hernan2020causal, p. 292].
:::
---
::: {#def-candidate-counterfactuals-tv}
## Candidate counterfactuals for two time points
By consistency ($Y^{A_0, A_1} = Y$ and $L_1^{A_0} = L_1$),
the model can be written in terms of observed data.
Replacing $\psi$ by a candidate value $\psi^\dagger$ defines the **candidate counterfactuals**
$$
\begin{aligned}
H_1(\psi^\dagger) &= Y - \paren{\psi_{11}^\dagger A_1 + \psi_{12}^\dagger A_1 L_1 + \psi_{13}^\dagger A_1 A_0 + \psi_{14}^\dagger A_1 A_0 L_1}, \\
H_0(\psi^\dagger) &= H_1(\psi^\dagger) - \psi_0^\dagger A_0.
\end{aligned}
$$
:::
::: {#rem-candidate-counterfactuals-true-psi}
## Candidate counterfactuals at the true parameter
When $\psi^\dagger = \psi$, $H_k(\psi^\dagger)$ equals the counterfactual $Y^{\bar{A}_{k-1}, \bar{0}_k}$:
observed treatment through $k - 1$ and no treatment afterwards.
Sequential exchangeability then says $H_k(\psi)$ is independent of $A_k$ given the past.
:::
---
::: {#exm-g-estimation-saturated}
## G-Estimation in the Sequentially Randomized Experiment (Fine Point 21.3)
**Time 1.**
Within each stratum of $(A_0, L_1)$, the mean of $H_1(\psi)$ must not depend on $A_1$.
The means of $H_1(\psi)$ by row of @tbl-seq-rand are $84$ and $84 - \psi_{11}$ in stratum $(0, 0)$,
$52$ and $52 - \psi_{11} - \psi_{12}$ in $(0, 1)$,
$76$ and $76 - \psi_{11} - \psi_{13}$ in $(1, 0)$, and
$44$ and $44 - \psi_{11} - \psi_{12} - \psi_{13} - \psi_{14}$ in $(1, 1)$
(for $A_1 = 0$ and $A_1 = 1$ respectively).
Equating within strata:
$$
\begin{aligned}
84 = 84 - \psi_{11} &\implies \psi_{11} = 0, \\
52 = 52 - \psi_{11} - \psi_{12} = 52 - \psi_{12} &\implies \psi_{12} = 0, \\
76 = 76 - \psi_{11} - \psi_{13} = 76 - \psi_{13} &\implies \psi_{13} = 0, \\
44 = 44 - \psi_{11} - \psi_{12} - \psi_{13} - \psi_{14} = 44 - \psi_{14} &\implies \psi_{14} = 0.
\end{aligned}
$$
**Time 0.**
With all $\psi_{1j} = 0$, $H_1(\psi) = Y$ and $H_0(\psi) = Y - \psi_0 A_0$.
Exchangeability $Y^{0,0} \ind A_0$ requires equal means of $H_0(\psi)$ in the two arms of 16,000 each:
$$
\begin{aligned}
A_0 = 0: &\quad 84 \times 0.25 + 52 \times 0.75 = 60, \\
A_0 = 1: &\quad (76 - \psi_0) \times 0.5 + (44 - \psi_0) \times 0.5 = 60 - \psi_0,
\end{aligned}
$$
so $60 = 60 - \psi_0$ and $\psi_0 = 0$.
All parameters are 0, so $\E{Y^g} = 60$ under every static or dynamic strategy $g$,
in agreement with the g-formula and IP weighting
[@hernan2020causal, pp. 292-293].
:::
### General Structural Nested Mean Models
::: {#def-candidate-counterfactuals-snmm}
## Candidate counterfactuals for a general SNMM
The candidate counterfactuals generalize to
$$H_k(\psi^\dagger) = Y - \sum_{j=k}^{K} A_j \, \gamma_j\paren{\bar{A}_{j-1}, \bar{L}_j, \psi^\dagger}.$$
:::
::: {.notes}
::: {#prp-snmm-h-mean}
## Mean of the candidate counterfactuals at the true parameter
Suppose the SNMM of @def-snmm-tv is correctly specified and consistency holds.
Even when local rank preservation is false (as it essentially always is when treatment has an effect),
$\E{H_k(\beta) \mid \bar{A}_k, \bar{L}_k} = \E{Y^{\bar{A}_{k-1}, \underline{0}_k} \mid \bar{A}_k, \bar{L}_k}$
and $\E{H_0(\beta)} = \E{Y^{\bar{0}}}$.
Hence, if $\hb$ is consistent for $\beta$,
$\E{Y^{\bar{0}}}$ is consistently estimated by the sample average of $H_0(\hb)$
[@hernan2020causal, p. 294].
:::
::: {#exm-blip-functions}
## Unsaturated blip functions
Unsaturated examples of the blip function:
$\beta_1$ (same effect for all histories and times),
$\beta_1 + \beta_2 k$ (effect varies linearly with time),
and $\beta_1 + \beta_2 k + \beta_3 a_{k-1} + \beta_4 l_k + \beta_5 l_k a_{k-1}$
(modification by the most recent treatment and covariate).
:::
:::
::: {.callout-note title="Technical Point 21.7: Marginal Structural Models and Structural Nested Models"}
An SNMM is a semiparametric MSM
if and only if $\gamma_k$ does not depend on $\bar{l}_k$;
the implied MSM is $\E{Y^{\bar{a}}} = \alpha_0 + \sum_{k=0}^{K} a_k \gamma_k(\bar{a}_{k-1}, \beta)$ with $\alpha_0 = \E{Y^{\bar{0}}}$.
Such an SNMM additionally assumes no effect modification by past covariates.
If that assumption holds, g-estimation is more efficient than IP weighting;
if it fails, g-estimation is biased while IP weighting of the MSM is not:
a variance-bias trade-off.
:::
### G-Estimation Algorithm (One Parameter)
::: {#alg-g-estimation-one-parameter}
## G-estimation of a one-parameter SNMM
Take $\gamma_k(\bar{a}_{k-1}, \bar{l}_k, \psi^\dagger) = \psi^\dagger$ (true value $\beta_1$), so $H_k(\psi^\dagger) = Y - \sum_{j=k}^{K} A_j \psi^\dagger$.
1. For each $\psi^\dagger$ on a grid from $\psi_{low}$ to $\psi_{up}$ (e.g., steps of 0.1),
compute $H_k(\psi^\dagger)$ for every person and time.
2. Fit the pooled logistic model
$\logit \Pr[A_k = 1 \mid H_k(\psi^\dagger), \bar{L}_k, \bar{A}_{k-1}] = \alpha_0 + \alpha_1 H_k(\psi^\dagger) + \alpha_2 W_k$,
where $W_k = w_k(\bar{L}_k, \bar{A}_{k-1})$ and each person contributes $K + 1$ rows.
3. The g-estimate $\hb$ is the $\psi^\dagger$ with $\hat{\alpha}_1$ closest to 0;
the 95% confidence interval is the set of $\psi^\dagger$ with $P > 0.05$ for $\alpha_1 = 0$.
:::
---
::: {#rem-g-estimating-equation}
## G-estimation as solving an estimating equation
Equivalently, $\hb$ is the $\psi^\dagger$ at which the score test of $\alpha_1 = 0$ has $P = 1$,
that is, the solution of
$$\sum_{i=1}^{N} \sum_{k=0}^{K} \cb{A_{i,k} - \expit\paren{\hat{\alpha}_0 + \hat{\alpha}_2 W_{i,k}}} H_{i,k}(\psi^\dagger) = 0,$$
with $\hat{\alpha}_0, \hat{\alpha}_2$ from the model fit with $\alpha_1 = 0$.
:::
::: {#prp-g-estimator-consistent}
## Consistency of the g-estimator
Let $\hb$ be the solution of the estimating equation of @rem-g-estimating-equation
(equivalently, the limit of the grid-search estimate of @alg-g-estimation-one-parameter as the grid step goes to 0).
Then $\hb$ is consistent for the true $\beta_1$ if
(i) the one-parameter SNMM of @alg-g-estimation-one-parameter,
$\gamma_k(\bar{a}_{k-1}, \bar{l}_k, \beta) = \beta_1$, is correct,
(ii) sequential exchangeability for $Y$ holds,
(iii) consistency holds, $Y^{\bar{A}} = Y$,
(iv) the treatment model $\logit \Pr[A_k = 1 \mid \bar{L}_k, \bar{A}_{k-1}] = \alpha_0 + \alpha_2 W_k$ is correct, and
(v) $H_k(\psi^\dagger)$ enters the treatment model linearly (Technical Point 14.2).
:::
::: {.notes}
::: {#rem-g-estimation-vector-beta}
## G-estimation of a vector parameter
**Vector $\beta$**: for the five-parameter blip
$\beta_0 + \beta_1 k + \beta_2 a_{k-1} + \beta_3 l_k + \beta_4 l_k a_{k-1}$,
the treatment model needs five terms involving $H_k(\psi^\dagger)$, e.g.
$\alpha_0 + H_k(\psi^\dagger)(\alpha_1 + \alpha_2 k + \alpha_3 A_{k-1} + \alpha_4 L_k + \alpha_5 L_k A_{k-1}) + \alpha_6 W_k$,
and $\hb$ solves the corresponding five-dimensional estimating equation.
The choice of these terms affects the width of the confidence interval, not consistency.
A grid search is impractical beyond about two dimensions
($20^5$ points for 20 values per component),
but for an SNMM linear in $\beta$ the estimating equation has a closed-form solution
(Technical Point 21.8, which also gives a $2^{K+1}$ multiply robust version)
[@hernan2020causal, pp. 295-296].
:::
:::
### From $\hb$ to Counterfactual Means
::: {#rem-counterfactual-means-after-g-estimation}
## Counterfactual means after g-estimation
Given $\hb$, how to estimate a counterfactual mean depends on the strategy and the blip function:
- $\E{Y^{\bar{0}}}$: the sample average of $H_0(\hb)$.
- If there is no effect modification by past covariates,
$\gamma_k(\bar{a}_{k-1}, \bar{l}_k, \beta) = \gamma_k(\bar{a}_{k-1}, \beta)$, and for a static $\bar{a}$
$$\hE{Y^{\bar{a}}} = \hE{Y^{\bar{0}}} + \sum_{k=0}^{K} a_k \gamma_k(\bar{a}_{k-1}, \hb).$$
- If $\gamma_k$ depends on $\bar{l}_k$, or $g$ is dynamic,
simulate $L_k$ from a model for $f(l_k \mid \bar{a}_{k-1}, \bar{l}_{k-1})$ (Technical Point 21.9).
:::
::: {.callout-note title="Technical Point 21.9: Estimating a Counterfactual Mean After g-estimation of a Structural Nested Mean Model"}
Estimate $\E{Y^{\bar{0}}}$ by the average of $H_0(\hb)$;
fit a model for $f(l_k \mid \bar{a}_{k-1}, \bar{l}_{k-1})$;
simulate covariate histories under $g$;
and add to $\hE{Y^{\bar{0}}}$ the average over simulations of
$\sum_{j=0}^{K} a_{v,j} \gamma_j(\bar{a}_{v,j-1}, \bar{l}_{v,j}, \hb)$.
The result is consistent for $\E{Y^g}$ if the covariate model and the SNMM are correct
and either the treatment model or the outcome model $\E{Y^{\bar{A}_{k-1}, \underline{0}_k} \mid \bar{L}_k, \bar{A}_{k-1}}$ is correct.
Under the null, $\gamma_j(\cdot, \hb)$ converges to 0 whenever $\hb$ is consistent for $\beta = 0$,
even when the covariate model is wrong.
So, under the identifiability conditions, this procedure **preserves the null**, unlike the parametric g-formula,
provided $\Pr[A_k = 1 \mid \bar{L}_k, \bar{A}_{k-1}]$ is known (as in a sequentially randomized experiment)
or either the treatment model or the outcome model is correct at each $k$
[@hernan2020causal, p. 297].
:::
::: {.callout-note title="Technical Point 21.13: Formal Definition of a General Structural Nested Mean Model"}
This technical point generalizes SNMMs to blips followed by an arbitrary strategy $g$ rather than "never treat"
(Robins 2004), and notes identification under a time-varying parallel-trends assumption
instead of sequential exchangeability [@hernan2020causal, p. 303].
:::
## 21.5 Censoring Is a Time-Varying Treatment (pp. 297-299)
---
Part II treated censoring $C$ as time-fixed.
::: {#def-time-varying-censoring-tv}
## Time-varying censoring
More realistically, censoring is time-varying: $C_1, C_2, \ldots, C_{K+1}$,
where $C_m = 0$ if the individual is still uncensored at $m$ and 1 otherwise.
- Censoring is **monotone**: $C_m = 0$ implies $C_1 = \cdots = C_{m-1} = 0$.
- $C_0 = 0$ for everybody, since otherwise they would not be in the study.
- After censoring, treatments, confounders, and outcomes are unobserved,
so analyses are restricted to uncensored person-times.
:::
---
::: {#def-g-formula-censoring}
## The g-formula with time-varying censoring
With time-varying censoring, the g-formula of @def-g-formula-tv for strategy $\bar{a}$ becomes
$$\sum_{\bar{l}} \E{Y \mid \bar{C} = \bar{0}, \bar{A} = \bar{a}, \bar{L} = \bar{l}}
\prod_{k=0}^{K} f\paren{l_k \mid \bar{c}_k = \bar{0}, \bar{a}_{k-1}, \bar{l}_{k-1}}.$$
:::
::: {#prp-g-formula-censoring-identifies}
## The censored g-formula identifies the uncensored counterfactual mean
If sequential exchangeability, positivity, and consistency hold
with $A_m$ replaced by the joint treatment $(A_m, C_{m+1})$ at every $m$,
then the g-formula of @def-g-formula-censoring equals $\E{Y^{\bar{a}, \bar{c} = \bar{0}}}$:
the mean outcome had everyone followed $\bar{a}$ and nobody been lost to follow-up.
:::
::: {.notes}
The superscript $\bar{c} = \bar{0}$ makes explicit the contrast
many investigators have in mind when they speak of "the effect of treatment"
[@hernan2020causal, p. 298].
::: {.callout-warning title="Conditioning on Being Uncensored Can Create Selection Bias"}
Conditioning on being uncensored induces selection bias under the null
when $C$ is a collider on a path between $A$ and $Y$, or a descendant of one (Chapter 8).
:::
:::
### IP Weighting with Censoring
::: {#alg-ipw-censoring-tv}
## IP weighting with time-varying censoring
Fit, for example, $\E{Y \mid \bar{A}, \bar{C} = \bar{0}} = \theta_0 + \theta_1 \mathrm{cum}(\bar{A})$
in the pseudo-population created by $W^{\bar{A}} \times W^{\bar{C}}$ or $SW^{\bar{A}} \times SW^{\bar{C}}$, where
$$W^{\bar{C}} = \prod_{k=1}^{K+1} \frac{1}{\Pr\sb{C_k = 0 \mid C_{k-1} = 0, \bar{A}_{k-1}, \bar{L}_{k-1}}},
\qquad
SW^{\bar{C}} = \prod_{k=1}^{K+1} \frac{\Pr\sb{C_k = 0 \mid C_{k-1} = 0, \bar{A}_{k-1}}}{\Pr\sb{C_k = 0 \mid C_{k-1} = 0, \bar{A}_{k-1}, \bar{L}_{k-1}}}.$$
Numerator and denominator are estimated with separate logistic models.
:::
---
::: {#rem-censoring-pseudo-population}
## The pseudo-population created by censoring weights
- **Nonstabilized** weights replace censored individuals by copies of uncensored ones
with the same treatment and covariate history:
the pseudo-population has the size of the study population **before** censoring,
and censoring is abolished.
- **Stabilized** weights keep the size **after** censoring
(the same proportion censored at each $k$),
but make censoring occur at random with respect to $\bar{L}_k$:
selection, but no selection bias.
- With nonstabilized weights the pseudo-population has no arrows from $L_k$ or $A_k$ into later $C_m$ ($m > k$).
With stabilized weights it has no arrows from $L_k$ into later $C_m$,
but censoring may still depend on past treatment,
because the numerator conditions on $\bar{A}_{k-1}$
(as for $SW^C$ in Chapter 12).
Either way, censoring no longer depends on the measured covariate history,
so IP weighting can estimate the joint effect of $(\bar{A}, \bar{C})$
even when $\bar{L}$ is affected by prior treatment.
:::
::: {#rem-g-estimation-censoring}
## G-estimation with censoring
**G-estimation with censoring**: first create, with nonstabilized weights $W^{\bar{C}}$,
a pseudo-population in which nobody is censored, then g-estimate in it.
:::
::: {.notes}
::: {.callout-tip title="Check the Mean of the Stabilized Weights"}
After fitting the censoring models, compute the sample mean of the estimated $SW^{\bar{C}}$.
Correctly specified models give a mean close to 1
[@hernan2020causal, p. 298],
so for $SW^{\bar{C}}$ a mean well away from 1 suggests that the numerator or the denominator model is misspecified.
The same check applies to $SW^{\bar{A}}$ with a narrower reading:
given a correct denominator $f_k$, taking iterated expectations over $a_k$ at each time $k$
gives $\E{\prod_k n_k / f_k} = 1$ for any proper conditional numerator $n_k$,
so for $SW^{\bar{A}}$ a mean away from 1 points to the denominator (treatment) model,
or to extreme weights from near-violations of positivity,
not to the numerator.
:::
::: {.callout-warning title="Fact-Check: Stabilized Weights Keep Arrows from Past Treatment"}
*Fact-check note*: the book states that "regardless of the type of IP weights used",
the pseudo-population has no arrows from $L_k$ and $A_k$ into future $C_m$ [@hernan2020causal, p. 299].
For stabilized weights this holds for $L_k$ only:
their numerator $\Pr[C_k = 0 \mid C_{k-1} = 0, \bar{A}_{k-1}]$ keeps the dependence of censoring on past treatment,
which is harmless because the analysis conditions on $\bar{A}$.
:::
:::
::: {.callout-note title="Technical Point 21.10: Survival Analysis with Time-Varying Treatments"}
To estimate the risk $\Pr[D_{k+1}^{\bar{a}, \bar{c} = \bar{0}} = 1]$,
either use the g-formula with models for the discrete-time hazards
and for the confounder densities, conditional on being event-free and uncensored
(e.g., pooled logistic models, as in Chapter 17),
or fit a pooled logistic hazards model in which each person-time receives the time-varying weight
$W_k^{\bar{A}} \times W_k^{\bar{C}}$ (or its stabilized version);
the latter estimates a marginal structural pooled logistic model (Robins 1998)
[@hernan2020causal, p. 299].
:::
## 21.6 The Big G-Formula (pp. 300-303)
---
This part of the book relies on sequential exchangeability given the measured $L$
and identification by the g-formula,
because few realistic longitudinal analyses use other identifying conditions
(such as the front door criterion).
But all identifying formulas are mathematically linked to the g-formula
based on **all** variables, measured and unmeasured.
::: {#def-big-g-formula}
## The Big G-Formula
Given a causal DAG with treatment $A$, outcome $Y$, measured covariates $L$, and unmeasured variables $U$,
let $X = (L, U)$.
Every parent of a treatment node is in $A$ or $X$, so $X$ ensures sequential exchangeability.
The **big g-formula** is the g-formula with $L$ replaced by $X$.
If positivity held, it would identify $\E{Y^g}$ under any strategy $g$
[@hernan2020causal, p. 300].
:::
::: {#def-factuals}
## Factuals
*Factuals* are variables of the actual world, in contrast to counterfactuals;
some factuals, like $U$, are unavailable for analysis because they were not measured
[@hernan2020causal, p. 300].
:::
::: {#rem-big-g-formula-not-computable}
## The big g-formula cannot be computed directly
The big g-formula is a function of the distribution of the factuals $(A, L, Y, U)$ (@def-factuals) alone,
but it includes unmeasured $U$, so it cannot be computed directly.
:::
### Reducing the Big G-Formula to Observed Data
::: {#rem-reducing-big-g-formula}
## When the big g-formula reduces to an observed-data formula
The key question: can the big g-formula be rewritten as a functional
of the distribution of the observed $(A, L, Y)$ only?
If so, the new formula reproduces the big g-formula and can be used in data analysis.
- Under the **front door** criterion (Figure 7.14), the big g-formula for $\E{Y^a}$
reduces to the front door formula (Technical Point 21.11).
- In general, whether such a reduction exists, and what the formula is,
using only the d-separations implied by the DAG,
was settled by Tian and Pearl (2002), Shpitser and Pearl (2006), and Huang and Valtorta (2006).
:::
::: {.callout-warning title="A Reduction Is Causal Only If the DAG Is"}
These are purely mathematical questions about distributions obeying d-separation;
they acquire a causal meaning only if the DAG is truly causal,
which we can never know for certain in an observational study.
:::
---
::: {#exm-front-door-big-g}
## The Front Door Formula from the Big G-Formula (Technical Point 21.11)
Under Figure 7.14 ($A \to M \to Y$, with $U$ a common cause of $A$ and $Y$),
the big g-formula for $\Pr[Y^a = y]$ is
$$\sum_m \sum_u \Pr[Y = y \mid M = m, A = a, U = u] \Pr[M = m \mid A = a, U = u] \Pr[U = u].$$
Step 1: by $A \ind Y \mid M, U$, drop $A = a$ from the first factor;
by $U \ind M \mid A$, replace $\Pr[M = m \mid A = a, U = u]$ by $\Pr[M = m \mid A = a]$;
and write $\Pr[U = u] = \sum_{a'} \Pr[U = u \mid A = a'] \Pr[A = a']$:
$$\sum_m \Pr[M = m \mid A = a] \sum_u \Pr[Y = y \mid M = m, U = u] \sum_{a'} \Pr[U = u \mid A = a'] \Pr[A = a'].$$
Step 2: by $U \ind M \mid A$, $\Pr[U = u \mid A = a'] = \Pr[U = u \mid M = m, A = a']$,
and by $A \ind Y \mid M, U$, $\Pr[Y = y \mid M = m, U = u] = \Pr[Y = y \mid M = m, A = a', U = u]$:
$$\sum_m \Pr[M = m \mid A = a] \sum_{a'} \cb{\sum_u \Pr[Y = y \mid M = m, A = a', U = u] \Pr[U = u \mid M = m, A = a']} \Pr[A = a'].$$
Step 3: the braces sum over $u$ to $\Pr[Y = y \mid M = m, A = a']$, giving the front door formula
$$\Pr[Y^a = y] = \sum_m \Pr[M = m \mid A = a] \sum_{a'} \Pr[Y = y \mid M = m, A = a'] \Pr[A = a'],$$
which involves observed variables only.
This proof does not require that counterfactuals $Y^m$ exist
[@hernan2020causal, p. 301].
:::
::: {.notes}
::: {#rem-front-door-other-proofs}
## Other proofs of the front door formula
Technical Point 21.11 also gives a coupling argument,
and Technical Point 21.12 a proof based on a SWIG property (Shpitser et al. 2022):
if a fixed treatment node $a_m$ is d-separated from $B^a$ given $C^a$ on the SWIG,
then $\Pr[B^a = b \mid C^a = c]$ does not depend on $a_m$
[@hernan2020causal, pp. 301-302].
:::
:::
## Summary
---
- In @tbl-seq-rand, the treatment has no effect;
the g-formula, IP weighting, and g-estimation all give $\E{Y^{1,1}} - \E{Y^{0,0}} = 60 - 60 = 0$,
where traditional methods failed.
- **G-formula**: standardize to the distribution of each $L_k$ given past treatment and covariates;
estimate with parametric models (the parametric g-formula), which is subject to the g-null paradox.
- **IP weighting**: $W^{\bar{A}}$ or $SW^{\bar{A}}$, a product over time,
fit marginal structural models by weighted regression;
$SW^{\bar{A}}$ preferred for unsaturated models.
- **Doubly robust**: sequential outcome regressions with the time-varying IP weight as a covariate,
consistent if either set of models is correct, and in fact $K + 2$ robust.
- **G-estimation**: one SNMM equation per time;
find $\psi^\dagger$ making $H_k(\psi^\dagger)$ independent of $A_k$ given the past.
- **Censoring** is handled like treatment: the joint strategy $(\bar{a}, \bar{c} = \bar{0})$ and weights $W^{\bar{C}}$ or $SW^{\bar{C}}$.
- The **big g-formula** uses all variables, measured and unmeasured;
other identifying formulas, such as the front door formula, are reductions of it to observed data.
## References
---
::: {#refs}
:::