19  G-Methods for Time-Varying Treatments

using DataFrames
using Distributions
using Random
using Statistics
using LinearAlgebra
using GLM
using Printf

Most chapters so far use one treatment decision. Many real problems are longitudinal. Treatment changes over time, covariates respond to past treatment, and those covariates then affect future treatment. Standard regression adjustment can fail in this setting.

Related reading: The g-estimation with time-varying covariates chapter of Topics on Econometrics and Causal Inference derives the identification result and applies the parametric g-formula to a worked binary-treatment example. The LMTP framework for continuous and modified-policy interventions is covered in Longitudinal modified treatment policy (LMTP) and in the Continuous Treatments chapter.

The main tools are:

  1. G-formula: simulate outcomes under a treatment regime.
  2. IPTW: reweight observations by the probability of their treatment history.
  3. Marginal structural models: fit a treatment-history regression in the weighted population.

For longitudinal TMLE, use R’s ltmle package; there is no Julia equivalent yet. The g-formula and IPTW/MSM translate directly to Julia with GLM.jl.

19.1 The time-varying confounding problem

Consider a two-period example:

  • \(L_0\): baseline covariate (e.g. health status at study entry)
  • \(A_0\): treatment at time 0 (e.g. medication)
  • \(L_1\): covariate at time 1, influenced by \(A_0\) (e.g. health status after first dose)
  • \(A_1\): treatment at time 1, depends on \(L_1\)
  • \(Y\): outcome at the end of follow-up

\(L_1\) is the problem variable. It confounds the effect of \(A_1\) on \(Y\), so we would like to adjust for it. But it is also affected by \(A_0\), so adjusting for it blocks part of the effect of \(A_0\). Regression without \(L_1\) is confounded; regression with \(L_1\) over-controls.

19.1.1 Simulating the failure of standard regression

Random.seed!(1)
n = 5000

L0 = randn(n)
A0 = Float64.(rand(n) .< @. 1 / (1 + exp(-(-0.5 + 0.8 * L0))))
L1 = @. 0.5 * L0 + 1.0 * A0 + randn()
A1 = Float64.(rand(n) .< @. 1 / (1 + exp(-(-0.5 + 0.8 * L1))))
Y  = @. 0.5 * A0 + 0.5 * A1 + 0.5 * L1 + 0.5 * L0 + randn()

df = DataFrame(L0 = L0, A0 = A0, L1 = L1, A1 = A1, Y = Y)

# A0 has a direct effect (0.5) AND an indirect effect through L1
# (A0 -> L1 -> Y contributes 1.0 * 0.5 = 0.5), so A0's total effect is 1.0;
# A1's total effect is 0.5; always-vs-never total = 1.5.
@printf("True total effect of always-treated vs never-treated: 1.500\n\n")

naive_no_L1   = lm(@formula(Y ~ A0 + A1 + L0),       df)
naive_with_L1 = lm(@formula(Y ~ A0 + A1 + L0 + L1), df)

@printf("Naive (no L1):   A0=%.3f, A1=%.3f, total=%.3f\n",
        coef(naive_no_L1)[2], coef(naive_no_L1)[3],
        coef(naive_no_L1)[2] + coef(naive_no_L1)[3])
@printf("Naive (with L1): A0=%.3f, A1=%.3f, total=%.3f\n",
        coef(naive_with_L1)[2], coef(naive_with_L1)[3],
        coef(naive_with_L1)[2] + coef(naive_with_L1)[3])
True total effect of always-treated vs never-treated: 1.500

Naive (no L1):   A0=0.944, A1=0.850, total=1.794
Naive (with L1): A0=0.474, A1=0.478, total=0.951

Neither regression recovers the true total effect:

  • Excluding \(L_1\) leaves \(L_1 \to A_1\) confounding that biases \(A_1\).
  • Including \(L_1\) blocks the \(A_0 \to L_1 \to Y\) pathway, so the coefficient on \(A_0\) now measures only the direct effect, not the total effect.

This is exactly the setting g-methods were designed for.

19.2 Parametric g-formula

The g-formula writes the mean outcome under a treatment regime as:

\[ \mathbb{E}[Y(\bar a)] = \sum_{\bar l} \mathbb{E}[Y \mid \bar a, \bar l] \prod_{t=0}^{K} f(l_t \mid \bar a_{t-1}, \bar l_{t-1}), \]

where \(\bar a=(a_0,a_1,\ldots,a_K)\) is a treatment regime and \(\bar l\) is the covariate history. In practice we fit models for the time-varying covariates and the outcome, then simulate the data under the treatment regime.

# Step 1: Fit nuisance models
mod_L1 = lm(@formula(L1 ~ L0 + A0), df)
mod_Y  = lm(@formula(Y  ~ L0 + A0 + L1 + A1), df)

# Step 2: Simulate counterfactual outcomes under a regime (a0, a1).
# The L1 distribution under A0 = a0 differs from the observed L1.
function g_compute(a0::Real, a1::Real; n_mc::Int = 100)
    sigma_L1 = sqrt(sum(abs2.(residuals(mod_L1))) / dof_residual(mod_L1))
    Y_sim = zeros(n)
    for _ in 1:n_mc
        df_a0 = DataFrame(L0 = L0, A0 = fill(a0, n))
        L1_sim = predict(mod_L1, df_a0) .+ sigma_L1 .* randn(n)
        df_sim = DataFrame(L0 = L0, A0 = fill(a0, n),
                           L1 = L1_sim, A1 = fill(a1, n))
        Y_sim .+= predict(mod_Y, df_sim)
    end
    Y_sim ./= n_mc
    return mean(Y_sim)
end

Random.seed!(99)
E_Y11 = g_compute(1.0, 1.0)
E_Y00 = g_compute(0.0, 0.0)
@printf("G-formula estimate (always vs never): %.3f  (true = 1.500)\n",
        E_Y11 - E_Y00)
G-formula estimate (always vs never): 1.483  (true = 1.500)

The important step is simulating \(L_1\) under the counterfactual treatment. We do not condition on the observed \(L_1\), because observed \(L_1\) is partly a consequence of observed treatment.

The g-formula can handle other regimes too. For example, treat only at \(t=0\):

Random.seed!(99)
E_Y10 = g_compute(1.0, 0.0)
@printf("Effect of treating only at t=0: %.3f  (true = 1.0)\n",
        E_Y10 - E_Y00)
Effect of treating only at t=0: 1.005  (true = 1.0)

19.3 IPTW for time-varying treatments

IPTW reweights observations by the inverse probability of their observed treatment history. The stabilized weight is:

\[ SW_i = \prod_{t=0}^{K} \frac{f(A_t \mid \bar A_{t-1})}{f(A_t \mid \bar A_{t-1}, \bar L_{t})}. \]

The numerator uses treatment history only. The denominator also conditions on the covariate history up to and including \(L_t\) — the covariates measured just before \(A_t\). (Skipping \(L_t\) would leave the \(L_t \to A_t\) confounding intact in the pseudo-population.)

# Numerator: P(A0) and P(A1 | A0)
n0 = glm(@formula(A0 ~ 1),  df, Binomial(), LogitLink())
n1 = glm(@formula(A1 ~ A0), df, Binomial(), LogitLink())

# Denominator: P(A0 | L0) and P(A1 | A0, L0, L1)
d0 = glm(@formula(A0 ~ L0),           df, Binomial(), LogitLink())
d1 = glm(@formula(A1 ~ A0 + L0 + L1), df, Binomial(), LogitLink())

p_n0 = predict(n0); p_n1 = predict(n1)
p_d0 = predict(d0); p_d1 = predict(d1)

# Probability of the OBSERVED treatment at each time
w0_num = ifelse.(df.A0 .== 1, p_n0, 1 .- p_n0)
w0_den = ifelse.(df.A0 .== 1, p_d0, 1 .- p_d0)
w1_num = ifelse.(df.A1 .== 1, p_n1, 1 .- p_n1)
w1_den = ifelse.(df.A1 .== 1, p_d1, 1 .- p_d1)

sw = (w0_num .* w1_num) ./ (w0_den .* w1_den)
@printf("Weight summary: min=%.3f, mean=%.3f, max=%.3f\n",
        minimum(sw), mean(sw), maximum(sw))

# Trim extreme weights at 99th percentile
sw_trim = min.(sw, quantile(sw, 0.99))
Weight summary: min=0.284, mean=1.000, max=13.489
5000-element Vector{Float64}:
 1.2320285379453222
 0.7918791756556957
 1.4947259119545413
 1.2957761198231386
 0.47742706181142075
 2.761282241893212
 0.9322305437789774
 0.758079673579434
 0.4446543045790307
 1.821169577919855
 1.0703793185873152
 1.0814119421933566
 0.8862804877411058
 ⋮
 1.1814410813547576
 0.9407057320445324
 0.6072055546765932
 3.0826283622213997
 0.6439023204591842
 1.2723812770139957
 1.700832121460976
 1.390763055139526
 1.1802015422293588
 1.7280793712231706
 0.612366299960335
 2.803538167023195

19.4 Marginal structural models

With stabilized weights, the marginal structural model is a weighted regression of \(Y\) on treatment history:

df.sw = sw_trim
msm_fit = lm(@formula(Y ~ A0 + A1), df, wts = df.sw)

println(coeftable(msm_fit))
@printf("\nMSM total effect (always vs never): %.3f  (true = 1.500)\n",
        coef(msm_fit)[2] + coef(msm_fit)[3])

# For comparison: the untrimmed weights
df.sw_u = sw
msm_untrimmed = lm(@formula(Y ~ A0 + A1), df, wts = df.sw_u)
@printf("Untrimmed MSM total effect:         %.3f\n",
        coef(msm_untrimmed)[2] + coef(msm_untrimmed)[3])
StatsBase.CoefTable(Any[[-0.04989628698526234, 1.0455936766250338, 0.5477314823226674], [0.028823347044086425, 0.03993523083912306, 0.03916614787775476], [-1.7311066237014057, 26.18223695355993, 13.984818829573053], [0.0834956511758907, 1.5200922797061257e-141, 1.3152335380742825e-43], [-0.10640295676856143, 0.9673027377450147, 0.47094829050850884], [0.006610382798036753, 1.1238846155050528, 0.624514674136826]], ["Coef.", "Std. Error", "t", "Pr(>|t|)", "Lower 95%", "Upper 95%"], ["(Intercept)", "A0", "A1"], 4, 3)

MSM total effect (always vs never): 1.593  (true = 1.500)
Untrimmed MSM total effect:         1.486

The weighted population is constructed so treatment is no longer associated with past covariates. That is why the MSM targets the total effect. Note the trimmed estimate sits a little above the truth while the untrimmed one is closer: the weight models here are correctly specified, so the untrimmed weights are consistent, and capping the top 1% of weights reintroduces a bit of confounding bias in exchange for lower variance. Trimming is a bias-variance trade, not a free lunch.

19.5 Hernán-style g-estimation: when to use which

Estimator Strengths Weaknesses
G-formula Handles arbitrary regimes; transparent Sensitive to model misspecification on \(L_t\) and \(Y\)
IPTW + MSM Simple, intuitive; valid with correct propensity model Variance large with extreme weights

In practice, reporting both is useful. Agreement is reassuring; disagreement suggests model misspecification in at least one of them.

For longitudinal TMLE, use R’s ltmle package. A Julia port would be useful but does not exist yet.

19.6 A robustness check: g-formula with covariates only

Even without IPTW + MSM, the g-formula is a useful check on standard regression. If standard regression and the g-formula disagree, that gap is a warning about time-varying confounding.

Random.seed!(123)
g_total = E_Y11 - E_Y00
msm_total = coef(msm_fit)[2] + coef(msm_fit)[3]

@printf("%-30s %.3f\n", "True total effect:", 1.5)
@printf("%-30s %.3f  (wrong)\n", "Standard regression (no L1):",
        coef(naive_no_L1)[2] + coef(naive_no_L1)[3])
@printf("%-30s %.3f  (wrong)\n", "Standard regression (with L1):",
        coef(naive_with_L1)[2] + coef(naive_with_L1)[3])
@printf("%-30s %.3f  (correct)\n", "G-formula:", g_total)
@printf("%-30s %.3f  (correct)\n", "MSM (IPTW):", msm_total)
True total effect:             1.500
Standard regression (no L1):   1.794  (wrong)
Standard regression (with L1): 0.951  (wrong)
G-formula:                     1.483  (correct)
MSM (IPTW):                    1.593  (correct)

The standard regressions are wrong in opposite directions; the g-methods both recover the truth. When standard and g-formula estimates diverge, the divergence is the bias.

19.7 When to reach for these methods

Use g-methods when:

  1. There is a time-varying confounder \(L_t\) that is affected by past treatment \(A_{t-1}\) and itself affects subsequent treatment \(A_t\).
  2. The research question concerns a sustained or dynamic treatment regime, not a single decision at a single time.

This includes many clinical follow-up studies and longitudinal studies of education or labor-market policies.

For a one-shot treatment, the methods from the Estimation and Nonparametric Causal Methods chapters are sufficient. The added complexity of g-methods is only justified when the time-varying confounder problem is present.

19.8 Summary

  • Standard regression fails when a time-varying confounder is itself affected by past treatment.
  • G-formula simulates outcomes under treatment regimes.
  • IPTW + MSM reweights the population and then fits a weighted treatment history model.
  • Longitudinal TMLE is best handled in R’s ltmle package for now.
  • Disagreement across estimators is useful information, not just a nuisance.