G-Estimation of Structural Nested Models

In plain language

G-estimation of structural nested models is a method for measuring how much a treatment that changes over time actually caused an outcome, in situations where a patient's lab value or disease marker both responds to earlier treatment and influences whether the doctor continues that treatment. Standard regression fails here because adjusting for that marker blocks part of the treatment effect while not adjusting leaves confounding unsolved. G-estimation gets around this by modeling the treatment decision process separately, then working backward to find the treatment effect that makes the adjusted outcome statistically independent of the treatment decision. It requires careful pre-specification of how the effect is structured and is harder to implement than conventional regression, but it is one of the few methods that can give a correct answer in this setting.

G-estimation of structural nested models (SNMs)

is one of the three core g-methods James Robins developed to handle the specific failure that defeats ordinary regression: a time-varying confounder that is simultaneously a consequence of prior treatment and a cause of future treatment and of the outcome (treatment-confounder feedback). When LDL falls because a patient took a statin, and that low LDL then drives the decision to keep prescribing, conditioning on LDL in a time-dependent Cox model blocks part of the treatment effect and opens a collider path; not conditioning on it leaves confounding. No single regression can be right. G-methods break the impasse by modeling the treatment process and the counterfactual outcome separately. Where the g-formula plugs in a fully parameterized outcome model and marginal structural models (MSMs) reweight to a pseudo-population, g-estimation targets a structural nested model: a parameterization of the incremental (blip) effect of the treatment received in interval k on the outcome, conditional on the observed history through k.

Core estimand distinction

A structural nested mean model (SNMM) parameterizes γ(k, h; ψ) = E[Y(ā_k, 0) − Y(ā_{k-1}, 0) | H_k = h, A_k] — the contrast, within the stratum defined by history H_k, between receiving the observed treatment in interval k and then nothing afterward, versus receiving nothing from k onward. For survival, the structural nested (accelerated) failure time model and the rank-preserving structural failure time model (RPSFTM) parameterize how treatment in each interval compresses or expands counterfactual survival time. The defining trick of g-estimation: under no-unmeasured-confounding-at-each-interval, the treatment-removed (blipped-down) counterfactual H(ψ) is conditionally independent of treatment A_k given history H_k. You find the ψ that makes the data obey that independence by solving the estimating equation U(ψ) = Σ_k [A_k − E(A_k | H_k)] · H(ψ) · q(H_k) = 0, where E(A_k | H_k) comes from a fitted treatment (propensity) model. This is the crucial contrast with MSMs: SNMs estimate a conditional effect that can directly carry effect modification by time-varying, post-baseline covariates (e.g., how the benefit of continuing a drug depends on current disease activity); MSMs estimate a marginal or baseline-modified effect and integrate that interaction away. The estimand is naturally a per-protocol / sustained-strategy effect — "the effect of taking treatment as assigned versus never" — which maps cleanly onto a target-trial per-protocol question.

Pros, cons, and trade-offs

- vs marginal structural models with IPTW (marginal-structural-models-g-methods): G-estimation inherits a partial double-robustness / null-preservation property — under the sharp null of no effect at any interval, it gives a valid test even if the blip form is wrong — and it does not require positivity as strictly as IPTW, because no extreme weights are formed; an SNM can therefore be far more efficient and stable when some covariate strata deterministically predict treatment. It uniquely estimates effect modification by time-varying covariates. Cost: the blip function and the treatment model must be specified; the estimating equation is solved numerically (grid search or root-finding), interpretation of blip parameters is harder to communicate than a hazard ratio, and off-the-shelf software is thin. Prefer g-estimation/SNMs when effect modification by post-baseline factors is the scientific target, when positivity is borderline and IPTW weights blow up, or when you want the null-robustness of a test. - vs the parametric g-formula: Both handle treatment-confounder feedback, but the g-formula requires correctly modeling the full joint density of all time-varying covariates and the outcome (the g-null paradox lurks when these are misspecified), whereas g-estimation only needs the treatment model plus the blip. Prefer the g-formula when you need absolute risks/survival curves under multiple hypothetical regimes; prefer g-estimation when the contrast of interest is a single sustained-strategy effect and you distrust a fully parametric outcome model. - vs time-dependent / weighted Cox (standard-cox-time-dependent, cox-ph-regression): A naive time-dependent Cox model that conditions on the time-varying confounder is biased under feedback (it blocks the mediated effect and conditions on a collider); g-estimation is consistent. Cost: a much higher technical and communication barrier, and parameters that are structural-model coefficients, not plug-and-play hazard ratios. Prefer standard methods only when the time-varying confounder is not itself affected by prior treatment — then there is no feedback and ordinary adjustment suffices. - vs clone-censor-weight per-protocol (clone-censor-weight-per-protocol): CCW is a design-forward way to get the same sustained-strategy estimand (clone, censor at deviation, weight for informative censoring); g-estimation is a modeling-forward route to it. Prefer CCW for transparency and easy communication of dynamic/grace-period strategies; prefer g-estimation for efficiency and for direct parameterization of effect modification.

When to use

Sustained or repeatedly-decided pharmacotherapy where (1) a time-varying covariate (LDL, HbA1c, eGFR, disease-activity score, prior adverse event, adherence) lies on the causal pathway and drives subsequent treatment, (2) the question is a per-protocol / "as-assigned-sustained" effect, and (3) you suspect the effect is modified by an evolving clinical factor. It is the natural estimation engine inside a target-trial emulation when the per-protocol estimand involves time-varying treatment.

When NOT to use — and when it is actively misleading or dangerous

- No treatment-confounder feedback exists. If the time-varying covariate is not affected by prior treatment, g-estimation is needless machinery; a correctly adjusted time-dependent Cox or an ACNU + PS analysis is simpler and equally valid. Reaching for g-estimation here trades clarity for nothing. - The blip and treatment models are both badly specified. Double robustness is partial, not magic: away from the null, a wrong blip and a wrong propensity model yield a confidently wrong, smooth answer with no warning. A "clean" point estimate from g-estimation is not evidence of correctness — this is its most dangerous failure mode, because nothing crashes. - Severe non-positivity in the treatment model. If E(A_k | H_k) is ≈0 or ≈1 in occupied strata, the estimating equation is uninformative there and the solution is driven by a few influential person-intervals; the apparent stability (no extreme weights to flag, unlike IPTW) hides the problem. - Discrete-time approximation of a continuous decision is too coarse. If treatment is re-decided daily but you bin into 90-day intervals, you misclassify the timing of the blip and of confounder feedback; the structural model no longer corresponds to any real intervention. - The outcome model / blip is being chosen by fit to the data. Tuning the blip until the answer looks plausible invalidates the null-robustness and inference. The SNM form must be pre-specified from subject-matter knowledge.

Data-source operational depth

- Claims (FFS or commercial): You must build a long-format, discrete-time person-interval table (one row per person per month/quarter) with a time-varying treatment indicator stitched from `fill_date` + `days_supply` (with explicit stockpiling/grace rules), time-varying confounders proxied from diagnoses/procedures/utilization, and the per-interval propensity model E(A_k | H_k). Failure modes: (1) Medicare Advantage person-time lacks FFS claims, so "no fill" can be missingness rather than non-treatment, corrupting both the treatment indicator and the confounder history that drive the estimating equation — restrict to enrollees with full A/B/D (or commercial pharmacy benefit) and exclude MA-only spans. (2) Lab values are largely absent in pure claims, so the very feedback-confounders that motivate g-estimation (LDL, HbA1c) are unobserved; without lab linkage you cannot satisfy the no-unmeasured-confounding-at-each-interval assumption and the method is being run on a fiction. (3) Differential informative censoring (disenrollment, death) by exposure requires IPCW layered onto g-estimation; in elderly cohorts competing mortality is heavy and must be handled, not ignored. (4) Immortal time sneaks in if interval 0 is mis-set before the first fill. - EHR: Supplies the time-varying labs/vitals that make the feedback structure measurable — its key advantage over claims — but visit-driven, irregular measurement means H_k is observed only when the patient shows up, inducing measurement timing that must be carried-forward or modeled; out-of-system care leaves exposure history incomplete, biasing the treatment model. Prefer EHR linked to claims so the fill history is complete and the confounder history is rich. - Registry: Protocolized repeated measures of treatment, severity, and adjudicated outcomes are an excellent substrate for the H_k history and are often used to validate claims-based g-estimation; the gap is usually complete pharmacy exposure (link to claims) and death (link to a vital-records index). - Linked claims–EHR–vital records: The ideal substrate — claims completeness for exposure, EHR labs for the feedback confounders, vital records for censoring — but linkage selection (only the linkable subset) and order/fill/service date discrepancies must be reconciled before interval boundaries and time zero are fixed.

Worked claims example

Question: the per-protocol effect of sustained high-intensity statin therapy versus never on first myocardial infarction over 24 months, in adults initiating a statin in a linked Medicare FFS + lab-feed cohort, where time-varying LDL and a statin-associated myalgia event drive discontinuation (classic treatment-confounder feedback). (1) Build a discrete-time table at monthly resolution: one row per `person_id` per month from the index fill. (2) Time-varying treatment `A_k` = covered by high-intensity statin in month k, derived from `fill_date` + `days_supply` with a 30-day grace period and stockpiling carried forward. (3) History `H_k` = baseline covariates (age, sex, prior CVD, diabetes) plus time-varying LDL (most recent lab carried forward), a myalgia/myopathy diagnosis flag, and prior-month treatment `A_{k-1}` — exactly the post-baseline factors a time-dependent Cox would mishandle. (4) Fit the per-interval treatment model E(A_k | H_k) by pooled logistic regression over person-months. (5) Specify an SNMM blip linear in treatment with an LDL interaction term to let the effect be modified by current LDL, then g-estimate ψ by solving Σ_k [A_k − Ê(A_k|H_k)]·H(ψ) = 0 over a grid (or by root-finding), bootstrapping persons for confidence intervals. (6) Layer IPCW for informative censoring at disenrollment and apply a competing-risks / death handling for the elderly; restrict to full A/B/D person-time and drop MA-only months so "no fill" is a true non-treatment, not missingness. (7) Report the sustained-vs-never effect and its modification by LDL, with sensitivity analyses on grace period, interval width, blip form, and a negative-control outcome to probe residual confounding — never tuning the blip to the result.

Interpreting the output

Using the worked example: the g-estimated blip parameter is psi ≈ −0.04 per treated month, representing the causal effect of one interval of sustained statin therapy on 12-month heart attack probability.

Formal interpretation: The g-estimate of psi ≈ −0.04 is the structural nested model blip-function parameter — the per-interval causal effect of treatment on the subsequent outcome, estimated by finding the value of psi at which the treatment residual (observed treatment minus predicted treatment probability from the propensity nuisance model) is uncorrelated with the blipped-down outcome (observed outcome adjusted to remove the posited treatment effect). Unlike the marginal coefficient from an MSM, psi is a conditional structural parameter describing the treatment effect at the individual level, net of the full time-varying LDL confounder history. It is valid under sequential exchangeability — at each month, the probability of statin receipt given the full observed covariate history is correctly modeled — and consistency: "one month of high-intensity statin" corresponds to a well-defined, replicable intervention. Confidence intervals derive from a bootstrap or the variance of the estimating equation, not from a parametric model formula.

Practical interpretation: Each month of sustained high-intensity statin therapy is estimated to reduce 12-month heart attack probability by approximately 4 percentage points, after correctly accounting for the LDL feedback loop that defeats standard regression. A standard time-fixed regression that controls for monthly LDL would block the very pathway (statin lowers LDL, lower LDL prevents heart attacks) that constitutes the drug's mechanism of benefit; g-estimation avoids this over-adjustment by modeling the treatment decision rather than conditioning on the intermediate confounder in the outcome model.

Worked example

Scenario

A cardiologist starts 200 patients on a statin and follows them for 12 monthly intervals to see who has a heart attack. Each month, the patient's LDL cholesterol is measured. The statin lowers LDL, but low LDL also makes the cardiologist more likely to keep prescribing. We want to know the true effect of sustained statin therapy on heart attack risk. A standard regression that controls for monthly LDL would be wrong — here is why, and how g-estimation fixes it.

Dataset

Three monthly snapshots for one patient illustrating how LDL sits between prior treatment and future treatment.

Steps

• WHY ORDINARY REGRESSION FAILS: In month 2, LDL dropped from 142 to 118 because the statin taken in month 1 worked. That lower LDL then influenced whether the doctor prescribed the statin again in month 3. LDL in month 2 is simultaneously a consequence of month-1 treatment and a cause of month-3 treatment and the eventual outcome.

• If we add monthly LDL as a covariate in a time-dependent regression, we block the path 'statin lowers LDL, lower LDL prevents heart attacks' — we would be subtracting out part of the very benefit we are trying to measure.

• If we instead leave LDL out to avoid blocking that path, the regression is confounded because sicker patients (higher LDL) get more aggressive treatment, which distorts the treatment-outcome comparison.

• Neither choice is correct with a single regression model. This is the core problem g-estimation is designed to solve.

• HOW G-ESTIMATION RECOVERS THE EFFECT: Instead of adjusting for LDL in an outcome model, g-estimation focuses on the treatment decision. At each month, it fits a model asking: given everything known about this patient at this moment (age, sex, prior LDL, prior treatment), what was the probability the doctor would prescribe the statin? The difference between the actual prescription (1 or 0) and that predicted probability is the treatment residual.

• G-estimation then proposes a candidate value for the statin effect — say, psi = -0.04, meaning each month of treatment reduces 12-month heart attack risk by 0.04 on the probability scale.

• It uses that candidate to mathematically remove the treatment effect from the observed outcomes, creating a blipped-down outcome: what each patient's outcome would have looked like if no one had ever received any treatment.

• If psi is the correct treatment effect, that blipped-down outcome should be statistically uncorrelated with the treatment residuals computed in step 5 — because once you remove the treatment effect, the leftover variation in outcomes should no longer track who got treated.

• The algorithm tests candidate values of psi — conceptually sweeping across a range — and finds the one value where the correlation between treatment residuals and blipped-down outcomes equals zero. That is the g-estimate.

• The result is a consistent estimate of the per-protocol effect of sustained statin therapy versus never treating, even though LDL was both affected by earlier statin use and drove later statin decisions.

Result

The g-estimated blip psi represents the causal effect of one interval of statin treatment on heart attack risk, adjusted for the full time-varying confounder history in a way that neither blocks the mediated pathway nor leaves confounding unaddressed. In the statin cohort the method would report something like psi = -0.04 per treated month (with a bootstrap confidence interval), interpreted as: each month of sustained high-intensity statin therapy reduces 12-month heart attack probability by approximately 4 percentage points, after correctly accounting for the LDL feedback loop that defeats standard regression.

Runnable example

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy.optimize import brentq

TX_FORMULA = "A_k ~ age + sex + ldl_k + myalgia_k + A_prev + k"

def g_estimate_snmm(df: pd.DataFrame, outcome_col: str = "Y") -> float:
    """Return psi: the per-interval additive blip on a standardized SNMM.

    df is long (one row per person-interval). A_k is the interval treatment; outcome_col is the
    person-level final outcome repeated across that person's rows.
    """
    d = df.copy()

    # (1) Per-interval treatment (propensity) model E(A_k | H_k), pooled over person-intervals.
    tx = smf.logit(TX_FORMULA, data=d).fit(disp=False)
    d["a_resid"] = d["A_k"] - tx.predict(d)          # A_k - Ehat(A_k | H_k)

    # Treatment dose carried by each person = sum of A_k; blip removes psi per treated interval.
    person = d.groupby("person_id").agg(dose=("A_k", "sum"),
                                        Y=(outcome_col, "first")).reset_index()
    d = d.merge(person[["person_id", "dose"]], on="person_id")

    def estimating_eq(psi: float) -> float:
        # H(psi) = blipped-down outcome; under correct psi it is independent of A_k given H_k,
        # i.e. the residual-weighted sum is zero in expectation.
        h = d[outcome_col] - psi * d["dose"]
        return float(np.sum(d["a_resid"] * h))

    # Solve U(psi) = 0 by bracketing root-find (widen bracket if needed for your outcome scale).
    return brentq(estimating_eq, -50.0, 50.0)

def bootstrap_ci(df: pd.DataFrame, n_boot: int = 500, seed: int = 1) -> tuple[float, float]:
    """Person-level bootstrap of psi (resample clusters, not rows)."""
    rng = np.random.default_rng(seed)
    ids = df["person_id"].unique()
    est = []
    for _ in range(n_boot):
        draw = rng.choice(ids, size=len(ids), replace=True)
        boot = pd.concat([df[df["person_id"] == i] for i in draw], ignore_index=True)
        est.append(g_estimate_snmm(boot))
    return float(np.percentile(est, 2.5)), float(np.percentile(est, 97.5))

library(gesttools)
library(data.table)

# gesttools expects an ordered long data.frame keyed by id and time, with outcome and exposure named.
dat <- FormatData(
  data       = as.data.frame(df),
  idvar      = "person_id",
  timevar    = "k",
  An         = "A_k",                 # time-varying exposure
  Cn         = NA,                    # add a censoring indicator name for IPCW
  outcomevar = "Y",
  varying    = c("ldl_k", "myalgia_k", "A_prev")
)

fit <- gestSingle(
  data          = dat,
  idvar         = "person_id",
  timevar       = "k",
  Yn            = "Y",
  An            = "A_k",
  Cn            = NA,
  outcomemodels = list("Y ~ A_k + age + sex + ldl_k"),  # blip; add A_k:ldl_k for effect modification
  propensitymodel = "A_k ~ age + sex + ldl_k + myalgia_k + A_prev + k",
  type          = 1                    # 1 = SNMM blip constant over time
)
print(fit$psi)        # the g-estimated blip parameter(s)

# --- Transparent fallback: solve sum_k (A_k - Ehat) * (Y - psi*dose) = 0 by hand ---
setDT(df)
df[, a_resid := A_k - predict(glm(A_k ~ age + sex + ldl_k + myalgia_k + A_prev + k,
                                  family = binomial, data = df), type = "response")]
df[, dose := sum(A_k), by = person_id]
U <- function(psi) df[, sum(a_resid * (Y - psi * dose))]
psi_hat <- uniroot(U, interval = c(-50, 50))$root