Tutorials › Marginal Structural Models › Worked Example With Time-Varying Confounding

The Parametric g-Formula · Tutorial 322 of 501

Worked Example With Time-Varying Confounding

Fit and combine longitudinal models to estimate regime-specific risks when a time-varying confounder is itself affected by prior treatment.

Advanced 12 min read

What You'll Learn

  • Identify how an intermediate covariate can be both affected by earlier treatment and confound later treatment.
  • Fit saturated transition and outcome models from observed cell counts.
  • Derive the parametric g-formula for two fixed treatment decisions.
  • Calculate and compare regime-specific risks by integrating over the right covariate distribution.
  • Distinguish a lack of observed support in a sample from a population positivity violation.

Why the Intermediate Covariate Must Be Simulated Under Each Regime

Suppose treatment is decided at two times. An intermediate health measure \(L_1\) is recorded between decisions. Earlier treatment \(A_0\) can change \(L_1\), while \(L_1\) predicts both later treatment \(A_1\) and the outcome \(Y\). Thus \(L_1\) is a time-varying confounder for the later treatment, but it is also affected by prior treatment. Estimating a risk under a treatment sequence must account for both facts.

The parametric g-formula does this by combining a model for the distribution of \(L_1\) under the first treatment with a model for the outcome conditional on treatment and \(L_1\). The key is to average the outcome predictions over the \(L_1\) distribution generated by the sequence being considered, rather than treating the observed distribution of \(L_1\) as fixed.

Definition: Let \(A_0,A_1\in\{0,1\}\) denote treatment at two decision times, let \(L_1\) be a covariate measured after \(A_0\) and before \(A_1\), and let \(Y\in\{0,1\}\) be the final outcome. For the fixed sequence \((a_0,a_1)\), write \(Y^{a_0,a_1}\) for the outcome that would occur if treatment were set to \(a_0\) and then \(a_1\). The regime-specific risk is \(P(Y^{a_0,a_1}=1)\).

The example below omits baseline covariates to keep the calculations visible. The same logic applies conditional on baseline covariates, followed by averaging over their baseline distribution. We assume consistency, no interference, sequential exchangeability at both treatment decisions, and positivity for the treatment choices and histories needed by the regimes. We also assume the transition and outcome models are correctly specified. These identification and model assumptions are substantive; fitting a model does not make them true.

Fitting the Transition and Outcome Models

Consider a dataset of 1,000 people. For this worked example, use saturated binomial models: one separate transition probability for each value of \(A_0\), and one separate outcome risk for each observed combination of \((A_0,L_1,A_1)\). In these models, the fitted probability in a cell is the number of events in that cell divided by its sample size.

Worked Example: Estimating the Intermediate-Covariate Distribution

Suppose the observed transition counts are as follows. The intermediate covariate is binary, and its fitted transition probabilities are the within-\(A_0\) proportions.

\(A_0\)Number with \(L_1=0\)Number with \(L_1=1\)TotalFitted \(P(L_1=1\mid A_0)\)
0400100500\(100/500=0.20\)
1200300500\(300/500=0.60\)

The fitted transition model therefore gives \(P(L_1=1\mid A_0=0)=0.20\) and \(P(L_1=1\mid A_0=1)=0.60\). The covariate distribution changes substantially with earlier treatment. Under the model, the probability of \(L_1=0\) is the complement in each group: \(0.80\) when \(A_0=0\), and \(0.40\) when \(A_0=1\).

Now suppose the outcome counts, split by \(A_0,L_1,A_1\), are these:

\(A_0\)\(L_1\)\(A_1\)Outcomes \(Y=1\)Cell totalFitted outcome risk
00020200\(20/200=0.10\)
00110200\(10/200=0.05\)
0101550\(15/50=0.30\)
0111050\(10/50=0.20\)
10020100\(20/100=0.20\)
10110100\(10/100=0.10\)
11075150\(75/150=0.50\)
11145150\(45/150=0.30\)

For example, the fitted risk for \((A_0,L_1,A_1)=(1,1,0)\) is \(75/150=0.50\). The saturated outcome model records a fitted risk for every history needed in the calculations below. It does not, by itself, produce a regime-specific risk: the fitted risks still need to be averaged over the appropriate transition distribution.

The Two-Time g-Formula

For binary \(L_1\) and \(Y\), the identified risk under a fixed sequence is

$$ P(Y^{a_0,a_1}=1) = \sum_{l\in\{0,1\}} P(Y=1\mid A_0=a_0,L_1=l,A_1=a_1) P(L_1=l\mid A_0=a_0). $$

The two factors have different roles. The outcome model predicts risk under the specified treatments and intermediate-covariate value. The transition model supplies the probability of that value under the first treatment in the sequence.

Theorem: Under consistency, no interference, sequential exchangeability at both treatment decisions, and positivity for the relevant treatment and covariate histories, the displayed two-time g-formula identifies \(P(Y^{a_0,a_1}=1)\) from the observed-data distribution.

Proof. By the law of total probability over the potential intermediate covariate \(L_1^{a_0}\),

$$ P(Y^{a_0,a_1}=1) = \sum_{l\in\{0,1\}} P(Y^{a_0,a_1}=1\mid L_1^{a_0}=l) P(L_1^{a_0}=l). $$

By exchangeability at the first decision and positivity, the distribution of \(L_1^{a_0}\) equals the observed distribution of \(L_1\) among those with \(A_0=a_0\): \(P(L_1^{a_0}=l)=P(L_1=l\mid A_0=a_0)\). For any \(l\) with positive probability under this intervention, exchangeability at the second decision identifies the conditional potential-outcome risk as \(P(Y^{a_0,a_1}=1\mid L_1^{a_0}=l)=P(Y=1\mid A_0=a_0,L_1=l,A_1=a_1)\). Consistency links the observed outcomes in those treatment-history cells to the corresponding potential outcomes. Substitution gives the stated formula. \(\square\)

Worked Example: Risks Under Two Fixed Sequences

Use the fitted models to estimate the risk under the sequence \((0,0)\). Under \(A_0=0\), the fitted probabilities for \(L_1=0\) and \(L_1=1\) are \(0.80\) and \(0.20\). The outcome-model risks under \((A_0,A_1)=(0,0)\) are \(0.10\) when \(L_1=0\) and \(0.30\) when \(L_1=1\). Thus

$$ \widehat P(Y^{0,0}=1) =(0.80)(0.10)+(0.20)(0.30) =0.08+0.06 =0.14. $$

The estimated risk under \((0,0)\) is \(0.14\), or 14%. Now consider the sequence \((1,1)\). Under \(A_0=1\), the fitted probabilities for \(L_1=0\) and \(L_1=1\) are \(0.40\) and \(0.60\). The outcome risks under \((A_0,A_1)=(1,1)\) are \(0.10\) and \(0.30\), respectively. Therefore

$$ \widehat P(Y^{1,1}=1) =(0.40)(0.10)+(0.60)(0.30) =0.04+0.18 =0.22. $$

The estimated risk under \((1,1)\) is 22%. These are estimated regime-specific risks, not a risk difference or risk ratio; those contrasts can be formed from the two risks if desired. The difference here is \(0.22-0.14=0.08\).

Worked Example: A Sequence With Different Treatments at the Two Times

The same fitted models can estimate a mixed sequence, such as \((1,0)\). Its intermediate-covariate distribution is still the one under \(A_0=1\): \(P(L_1=0\mid A_0=1)=0.40\) and \(P(L_1=1\mid A_0=1)=0.60\). For \(A_1=0\), the fitted outcome risks under \(A_0=1\) are \(0.20\) at \(L_1=0\) and \(0.50\) at \(L_1=1\). Hence

$$ \widehat P(Y^{1,0}=1) =(0.40)(0.20)+(0.60)(0.50) =0.08+0.30 =0.38. $$

This calculation illustrates why the transition model is indexed by the first treatment in the sequence, while the outcome model uses both treatments. Substituting the \(A_0=0\) transition probabilities here would describe a different distribution of \(L_1\) and would not estimate the risk under \((1,0)\).

Why This Is Not Ordinary Covariate Adjustment

A conventional regression that includes \(L_1\) can estimate conditional outcome associations, but a conditional risk is not automatically a marginal risk under a fixed sequence. The g-formula averages those conditional predictions using the intermediate-covariate distribution that would arise under the regime. This matters precisely because earlier treatment can change \(L_1\).

For example, if one incorrectly used the \(A_0=0\) transition probabilities to calculate the risk under \((1,1)\), the result would be \((0.80)(0.10)+(0.20)(0.30)=0.14\), rather than the correct model-based calculation of \(0.22\). The outcome risks are the same as in the \((1,1)\) calculation, but the covariate weights describe the wrong first treatment. Conversely, simply conditioning on \(L_1\) and reporting separate risks by \(L_1\) does not answer the marginal regime-risk question.

The g-formula is one approach among several. Inverse-probability weighting represents a regime through weights based on treatment probabilities, while the parametric g-formula models covariate transitions and outcomes and then averages predictions. Ordinary regression adjustment is not a substitute for this marginalization when the covariate is treatment-affected. The g-formula’s advantage is that it makes the implied covariate distribution explicit; its limitation is that errors in the transition or outcome models can distort the estimate.

Diagnostics and Failure Modes

Check the fitted models as conditional models, not only through the final risk. Compare predicted and observed transition proportions across treatment histories, and compare predicted and observed outcomes across the histories used in the outcome model. For a saturated binary model, the fitted cell risks reproduce the cell proportions, but sparse cells can still make those proportions unstable. With less saturated models, inspect whether predicted probabilities are plausible and whether the model captures important differences across histories.

Also inspect support. A cell with very few observations can make its estimated risk unstable. A cell with no observations signals a lack of empirical support and prevents direct estimation of that cell’s risk from this sample. It does not by itself establish a population positivity violation: a population positivity violation occurs when the relevant population treatment probability is zero. Sparse observed data can be consistent with a small but positive population probability. Either way, extrapolating a model into unsupported histories requires caution.

A further limitation is that this simple example assumes no loss to follow-up and uses binary \(L_1\). In an applied analysis, missingness, additional covariates, or more complex longitudinal histories require corresponding models and assumptions. The introductory material on parametric g-computation covers the general simulation and integration framework; here the focus is the particular role of a treatment-affected confounder in a two-decision calculation.

Derived result: If each relevant history has positive probability, the sample is independent and identically distributed, and the saturated cell risks are estimated by observed proportions, then the plug-in g-formula estimate converges in probability to the identified risk.

Proof. For each finite cell, its sample proportion converges in probability to its population probability by the law of large numbers. If a cell probability is positive, its sample count divided by \(n\) converges to a positive limit, so the corresponding conditional sample proportion (events divided by cell count) converges in probability to the conditional cell probability. The estimated transition probabilities and outcome risks therefore converge jointly in probability to their population values. The g-formula estimate is a finite sum of products of these probabilities. Addition and multiplication are continuous, so the continuous mapping theorem gives convergence in probability of that sum to the population g-formula risk. By the identification result, this limit is \(P(Y^{a_0,a_1}=1)\). \(\square\)

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. In the two-decision setup, why must the distribution used to average over \(L_1\) depend on \(a_0\)?
  2. Suppose \(P(L_1=1\mid A_0=1)=0.70\), and the outcome risks under a fixed sequence are \(0.20\) at \(L_1=0\) and \(0.50\) at \(L_1=1\). Calculate the g-formula risk.
  3. Does a treatment-history cell with no observations in a finite sample prove a population positivity violation?
  4. What are the two conditional models combined in the two-time parametric g-formula?

Solutions. (1) Earlier treatment can affect \(L_1\), so the covariate distribution under the regime must be \(P(L_1\mid A_0=a_0)\), not a distribution taken from a different first treatment. (2) The complementary probability for \(L_1=0\) is \(0.30\), so the risk is \((0.30)(0.20)+(0.70)(0.50)=0.06+0.35=0.41\). (3) No. It indicates a lack of empirical support and prevents direct estimation from that sample; population positivity fails only if the relevant population treatment probability is zero. (4) A transition model for \(L_1\) given \(A_0\), and an outcome model for \(Y\) given \(A_0,L_1,A_1\).

Clinical Biostats