G-Computation and the Parametric G-Formula

In plain language

G-computation estimates how much a treatment changes an outcome by fitting a statistical model that predicts the outcome from patient characteristics and treatment status, then using that model to predict each patient's result twice — once as if treated and once as if untreated — and averaging across everyone in the study. The difference between the two averages is a treatment effect expressed as a risk difference or risk ratio on a plain probability scale, a number a clinician or payer can act on directly. This predict-everyone-twice-and-average approach sidesteps a problem specific to logistic regression where odds ratios cannot be simply averaged across patient groups or turned directly into risk reductions, even when there is no confounding at all. Confidence intervals come from bootstrapping — repeating the whole procedure hundreds of times on resampled data to quantify the spread of plausible answers.

The standardization idea — fit, predict, average

G-computation (also called model-based standardization or outcome regression standardization) estimates the causal effect of a treatment by exploiting the counterfactual structure embedded in any outcome regression model. The procedure has three steps: (1) fit an outcome model — logistic for binary endpoints, Poisson or negative-binomial for counts, Cox for time-to-event, linear for continuous scores — on the observed cohort, including treatment and all measured confounders as predictors; (2) duplicate the dataset twice, setting treatment = 1 in the first copy and treatment = 0 in the second, and generate model predictions for every patient under each assignment regardless of what they actually received; (3) average the predicted outcomes under treatment across all patients, average the predicted outcomes under no treatment, and contrast the two averages as a risk difference, risk ratio, or rate ratio. The result is the population-averaged causal effect estimate — the difference in average outcomes you would expect if, contrary to fact, the entire cohort had been treated versus if none had been.

The intuition maps directly onto direct standardization or post-stratification: segment the cohort into covariate cells, estimate the expected outcome in each cell under each treatment from the fitted model, and re-weight by the cell sizes in the target population. G-computation performs this standardization continuously over the full covariate distribution of the outcome model rather than over a small number of discrete strata. The "G" in g-computation refers to Robins's general product formula for the joint density of a time-varying treatment history — the same mathematical foundation underlying the entire g-methods family: IPTW-estimated marginal structural models (MSMs), g-estimation of structural nested models, and the parametric g-formula.

Why g-computation returns marginal effects on any scale — and how it resolves noncollapsibility

Logistic regression reports a conditional odds ratio (OR): the effect within strata of the adjustment covariates, on the odds scale. Noncollapsibility means that even in the complete absence of confounding, the conditional OR is not a weighted average of stratum-specific ORs; it changes as more covariates enter the model not because the truth changed but because the reference point of "holding covariates fixed" shifted. The fully-adjusted OR from a logistic model cannot be directly converted to a population-level risk difference or risk ratio without additional steps, and comparing ORs across studies or populations requires careful accounting for the baseline risk. The same noncollapsibility problem affects Cox hazard ratios: a fully-adjusted HR changes as the covariate set expands, even without additional confounding control, and is therefore not directly comparable to the marginal effect a payer needs.

G-computation solves both problems by standardizing over the observed covariate distribution: it produces a marginal risk difference or risk ratio from the same logistic or Cox model that would otherwise report only a conditional OR or HR. The marginal risk difference is scale-commensurable with baseline risk, directly interpretable for number-needed-to-treat (NNT) calculations, and appropriate input for cost-effectiveness models that require absolute event probabilities. The analyst need never report the OR as the primary causal effect estimate — g-computation extracts the population-relevant quantity from the same fitted model.

Inference via the bootstrap

G-computation is a plug-in estimator: point estimates are computed by substituting fitted model parameters into the standardization formula. Analytic standard errors via the delta method exist for simple parametric models but are rarely used in RWE applications. The nonparametric bootstrap is the standard approach because it (a) propagates model-fitting uncertainty into the final interval without assumptions about the sampling distribution of the marginal contrast, (b) works identically for any outcome scale — RD, RR, NNT — without re-deriving variance formulas, and (c) accommodates multi-model or survival-model g-computation settings where analytic variances are intractable. The procedure resamples the cohort with replacement, re-fits the outcome model, and re-applies the standardization in each resample; the 2.5th and 97.5th percentiles of the bootstrap distribution form the 95% confidence interval. For large datasets (n > 100,000), 500 to 1,000 resamples are typically sufficient; for small samples, bias-corrected and accelerated (BCa) intervals should be preferred over the naive percentile method to improve coverage.

The time-varying generalization: the parametric g-formula

The point-treatment version above assumes treatment is decided once at baseline. When treatment is time-varying — a patient starts, stops, escalates, or switches over follow-up — and post-baseline health measures are simultaneously confounders of future treatment and mediators of past treatment (treatment-confounder feedback), neither standard regression nor point-treatment g-computation can give an unbiased causal estimate. Conditioning on a post-baseline variable that is on the causal path from past treatment to the outcome blocks part of the true effect; omitting it leaves future treatment confounded. No regression specification escapes both horns.

The parametric g-formula (Robins 1986) resolves this by modeling the joint distribution of the time-varying confounders, treatment, and outcome at each interval, and then Monte Carlo simulating the counterfactual history forward under each hypothetical treatment strategy. At each time step, simulated covariate values are drawn from the fitted covariate models conditional on the simulated — not the observed — history, generating potential-outcome trajectories under the target strategy without the confounding present in the natural course. The result is the full distribution of the counterfactual outcome under the hypothetical strategy, including absolute cumulative incidence curves and risk differences at any follow-up horizon. The parametric g-formula is uniquely suited to dynamic treatment regimes (e.g., "treat whenever LDL exceeds 130 mg/dL"), a capability that IPTW-estimated MSMs handle less naturally.

Relationship to IPW-MSMs and doubly-robust methods

G-computation and IPTW-estimated MSMs are dual approaches to the same causal estimand: g-computation models the outcome process (the conditional distribution of Y given treatment and confounders), while IPTW models the treatment process (the probability of observed treatment given confounders). Each is consistent if its respective model is correctly specified; neither has double robustness. The comparison with TMLE is sharper: TMLE is formally g-computation plus a targeting (fluctuation) step driven by the propensity score so the final estimate solves the efficient influence-curve equation and achieves double robustness. In practice, g-computation with a well-specified outcome model and bootstrap inference is a transparent, defensible primary analysis; TMLE or AIPW is the natural doubly-robust sensitivity check, especially when high-dimensional confounding motivates Super Learner nuisance fits.

G-null paradox

Under a sharp null (treatment truly has no effect) with treatment-confounder feedback, a parametric g-formula can be mathematically guaranteed to be misspecified — the g-null paradox. When the null is plausible and feedback is present, g-estimation of structural nested models is robust to this paradox; IPTW-MSMs and TMLE are unaffected because they model the treatment process rather than the covariate evolution conditionally on treatment.

Identifying assumptions

G-computation requires three assumptions that data cannot confirm: (1) exchangeability — no unmeasured confounding; all variables that jointly predict treatment assignment and the outcome are measured and included in the outcome model; (2) positivity — every patient in the target population has a positive probability of receiving each treatment level conditional on measured covariates; and (3) consistency — the treatment in the data corresponds to a well-defined, replicable intervention with no interference between patients. Modeling does not buy identification: a g-computation estimate from a confounded observational dataset is a biased quantity, not a causal effect. The identifying assumptions are identical to those for IPTW, TMLE, and any other standard observational causal method — the approach to estimation changes, but the assumptions needed to claim causation do not.

Pros, cons, and trade-offs

Pros: G-computation is the most flexible route to marginal effects on any scale from any outcome model — risk difference, risk ratio, NNT, expected count difference — without additional software or re-derived variance formulas. It naturally handles treatment-covariate interactions because predictions are generated under each patient's own covariate profile, not at artificially fixed covariate means. The parametric g-formula handles dynamic strategies and treatment-confounder feedback, producing absolute cumulative incidence curves that IPTW-MSMs cannot easily generate. Bootstrap CIs are straightforward and generalizable to any contrast. The output — a marginal risk difference — is the number a payer, HTA body, or clinical guideline committee needs for formulary and value decisions.

Cons: G-computation is consistent only if the outcome model is correctly specified — there is no robustness fallback if the model is wrong, unlike doubly-robust estimators such as TMLE. The parametric g-formula requires fitting separate models for each time-varying confounder, multiplying model-specification risk with each additional variable. Bootstrap inference is computationally expensive for large cohorts or complex g-formula runs. G-computation is less familiar than propensity-score matching in some regulatory and payer review contexts.

Trade-offs vs IPTW: G-computation avoids positivity diagnostics and extreme-weight truncation because it extrapolates from the outcome model rather than reweighting. However, this also means it extrapolates into covariate regions where no patients were treated — a risk if the model is misspecified in those regions. Use IPTW as a cross-check when the treatment process is well-modeled but the outcome model is uncertain; use g-computation when the outcome model is credible and extreme weights would be a concern.

When to use

Use g-computation when the primary deliverable is a marginal risk difference, risk ratio, or rate ratio from an observational cohort: comparative effectiveness research where absolute benefit estimates are required for formulary or guideline decisions; target-trial emulation where counterfactual risks under sustained strategies must be reported on the absolute scale; any analysis where logistic or Poisson regression is the underlying model but the conditional OR or IRR is not the quantity the audience needs. Use the parametric g-formula for time-varying treatment with confounder feedback, dynamic treatment rules, or when absolute cumulative incidence curves under competing strategies are required.

When NOT to use

Do not use point-treatment g-computation when treatment is time-varying with treatment-confounder feedback — that setting requires the parametric g-formula or an IPTW-MSM. Do not interpret a g-computation confidence interval from a misspecified or confounded outcome model as a valid causal uncertainty range — model diagnostics and doubly-robust sensitivity analyses are mandatory before claiming causal inference. Do not omit the outcome model specification, covariate set, and target population from reporting — these choices define the estimand and must be pre-specified. Do not use g-computation to rescue analyses with substantial unmeasured confounding; modeling does not buy causal identification.

Interpreting the output

In the worked example, the g-computation risk difference is 0.18 - 0.28 = -0.10.

(1) Formal interpretation. The g-computation estimand is the population-averaged causal risk difference — the average difference in 1-year stroke probability if, contrary to fact, the entire cohort had received the drug versus if none had received it, averaging over the observed distribution of disease severity in the cohort. The estimate of -0.10 is a marginal (population-averaged) effect, not the conditional log-odds ratio that the underlying logistic model produces within each severity stratum. It is consistent under three untestable assumptions: exchangeability (no unmeasured confounding given severity), positivity (all patients had a non-zero probability of receiving each assignment), and consistency (the drug is well-defined). A bootstrap 95% confidence interval means that across repeated hypothetical samples from the same data-generating process, the constructed intervals would contain the true marginal risk difference in 95% of replications.

(2) Practical interpretation. If this blood-pressure drug were given to every patient in the cohort instead of none, the outcome model predicts that average 1-year stroke risk would fall from 28% to 18% — an absolute reduction of 10 percentage points. Translating to number-needed-to-treat: treating 10 patients prevents, on average, one stroke. This number is directly usable in a budget-impact model or cost-effectiveness analysis in a way that a logistic OR cannot be without additional computation. The causal claim is conditional on the assumption that disease severity captures all confounding between drug assignment and stroke risk in this dataset.

Worked example

Scenario

A pharmacoepidemiologist wants to estimate the causal effect of a new blood-pressure drug on 1-year stroke risk in a cohort of 1,000 patients. Patients have one of two disease-severity profiles: mild (stratum A, 60% of the cohort) and severe (stratum B, 40%). A logistic regression outcome model is fit on the observed data with treatment and severity as predictors. The analyst applies g-computation: predict each patient's 1-year stroke probability under treatment (treated=1) and under no treatment (treated=0), then average both sets of predictions over the full cohort and contrast the averages.

Dataset

Stratum-level summary of fitted 1-year stroke risks from the logistic outcome model, by disease severity. These are the model-predicted risks that g-computation standardizes over the cohort; the 60%/40% split reflects the observed stratum proportions.

Steps

• Step 1 — Fit the outcome model: run logistic regression of stroke (1=yes, 0=no) on treatment (1=yes, 0=no) and severity. The fitted model produces a predicted stroke probability for every patient under any treatment value.

• Step 2 — Set everyone to treated: copy the dataset and set treated=1 for every patient. Score with the fitted model. Stratum A patients (mild) get predicted risk 0.10; stratum B patients (severe) get 0.30.

• Step 3 — Set everyone to untreated: copy the dataset and set treated=0 for every patient. Score with the fitted model. Stratum A gets 0.20; stratum B gets 0.40.

• Step 4 — Standardize the treated risks over the cohort proportions: 0.60.10 + 0.40.30 = 0.06 + 0.12 = 0.18

• Step 5 — Standardize the untreated risks over the cohort proportions: 0.60.20 + 0.40.40 = 0.12 + 0.16 = 0.28

• Step 6 — Contrast the standardized risks: 0.18 - 0.28 = -0.10

• Step 7 — Bootstrap for the 95% CI: resample the full cohort with replacement 1,000 times, re-fit the outcome model and re-apply Steps 2-6 in each resample; the 2.5th and 97.5th percentiles of the resulting distribution of risk differences give the confidence interval.

Result

standardized_risk_treated = 0.60.10 + 0.40.30 = 0.18. standardized_risk_untreated = 0.60.20 + 0.40.40 = 0.28. RD = 0.18 - 0.28 = -0.10. The drug is estimated to lower 1-year stroke risk by 10 percentage points on average across the cohort — the population-averaged marginal risk difference, not the conditional odds ratio that the logistic model reports by default.

Runnable example

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

def g_computation(df: pd.DataFrame, formula: str, n_boot: int = 1000, seed: int = 42) -> dict:
    """
    Point-treatment g-computation with bootstrap 95% CI.

    Parameters
    ----------
    df      : DataFrame with outcome, treatment indicator, and confounders
    formula : statsmodels formula, e.g. 'event ~ treated + severity'
    n_boot  : number of bootstrap resamples (500+ for large datasets)
    seed    : random seed for reproducibility

    Returns
    -------
    dict with risk_treated, risk_untreated, RD, RR, CI_95
    """
    rng = np.random.default_rng(seed)

    def _one_run(data: pd.DataFrame):
        m = smf.logit(formula, data=data).fit(disp=False)
        d1 = data.copy(); d1["treated"] = 1   # everyone treated
        d0 = data.copy(); d0["treated"] = 0   # everyone untreated
        r1 = m.predict(d1).mean()              # marginal risk under treatment
        r0 = m.predict(d0).mean()              # marginal risk under no treatment
        return r1, r0, r1 - r0

    r1, r0, rd = _one_run(df)

    boot_rds = []
    n = len(df)
    for _ in range(n_boot):
        idx = rng.integers(0, n, size=n)
        _, _, rd_b = _one_run(df.iloc[idx].reset_index(drop=True))
        boot_rds.append(rd_b)

    ci = np.percentile(boot_rds, [2.5, 97.5])
    rr = r1 / r0 if r0 > 0 else float("nan")
    return {"risk_treated": r1, "risk_untreated": r0, "RD": rd, "RR": rr, "CI_95": ci.tolist()}

# ── Synthetic cohort: severity confounds treatment and stroke outcome ──
rng0 = np.random.default_rng(7)
n = 400
severity = rng0.binomial(1, 0.40, n)                               # 40% severe
p_treat  = 1 / (1 + np.exp(-(-0.5 + 1.2 * severity)))             # severe patients more likely treated
treated  = rng0.binomial(1, p_treat)
p_event  = 1 / (1 + np.exp(-(-2.5 + 1.5 * treated + 2.0 * severity)))
event    = rng0.binomial(1, p_event)
df = pd.DataFrame({"event": event, "treated": treated, "severity": severity})

result = g_computation(df, "event ~ treated + severity", n_boot=1000)
print(f"Risk (treated):    {result['risk_treated']:.3f}")
print(f"Risk (untreated):  {result['risk_untreated']:.3f}")
print(f"Risk difference:   {result['RD']:.3f}")
print(f"Risk ratio:        {result['RR']:.3f}")
print(f"95% bootstrap CI:  ({result['CI_95'][0]:.3f}, {result['CI_95'][1]:.3f})")

# For the time-varying parametric g-formula with treatment-confounder feedback, use
# R's gfoRmula package (McGrath et al. 2020). Python ports exist but gfoRmula is the
# production-grade implementation; see marginal-structural-models-g-methods entry.

# ── Synthetic cohort: severity (0=mild, 1=severe) confounds treatment and stroke outcome ──
set.seed(7)
n        <- 400
severity <- rbinom(n, 1, 0.40)
treated  <- rbinom(n, 1, plogis(-0.5 + 1.2 * severity))
event    <- rbinom(n, 1, plogis(-2.5 + 1.5 * treated + 2.0 * severity))
df       <- data.frame(event = event, treated = treated, severity = severity)

# ── Manual g-computation: fit, score two counterfactual datasets, average ──
gcomp_once <- function(data) {
  m  <- glm(event ~ treated + severity, data = data, family = binomial)
  d1 <- data; d1$treated <- 1          # everyone treated
  d0 <- data; d0$treated <- 0          # everyone untreated
  r1 <- mean(predict(m, newdata = d1, type = "response"))
  r0 <- mean(predict(m, newdata = d0, type = "response"))
  c(risk_treated = r1, risk_untreated = r0, RD = r1 - r0, RR = r1 / r0)
}
point_est <- gcomp_once(df)

# ── Bootstrap 95% percentile CI ──
B        <- 1000
boot_rds <- replicate(B, gcomp_once(df[sample(nrow(df), replace = TRUE), ])["RD"])
ci       <- quantile(boot_rds, c(0.025, 0.975))

cat(sprintf("Risk (treated):   %.3f\n",  point_est["risk_treated"]))
cat(sprintf("Risk (untreated): %.3f\n",  point_est["risk_untreated"]))
cat(sprintf("Risk difference:  %.3f  95%% CI (%.3f, %.3f)\n",
            point_est["RD"], ci[1], ci[2]))
cat(sprintf("Risk ratio:       %.3f\n",  point_est["RR"]))

# ── One-liner with stdReg (marginal risks + bootstrap CI via stdGlm) ──
# library(stdReg)
# fit <- glm(event ~ treated + severity, data = df, family = binomial)
# std <- stdGlm(fit = fit, data = df, X = "treated")
# summary(std, CI.type = "plain")   # reports marginal risks and RD with bootstrap CI

# ── Parametric g-formula for time-varying treatment (treatment-confounder feedback) ──
# library(gfoRmula)   # McGrath et al. 2020 — Patterns 1(3):100008
# See marginal-structural-models-g-methods entry for a full gfoRmula example with
# time-varying confounder models and dynamic treatment strategies.