Marginal estimates and G-computation

Two promises this page is here to keep

The logistic regression page, when it got to non-collapsibility, said that “if you want an effect measure that can be interpreted at the population level (to compute an NNT, say), report a marginal estimate”. Its table of common mistakes says the same thing again: taking the model’s OR and computing an NNT from it is wrong, and what you need is “standardisation / G-computation”.

Neither sentence says how. This page does.

Same MASS::birthwt, same binary outcome (birth weight under 2500 g), and deliberately the same covariates as that page, so the numbers on the two pages can be compared directly.

Why a conditional OR cannot be turned into an NNT

Exponentiating a logistic coefficient gives an OR, and specifically a conditional OR. Read out in full, it says: “holding age, mother’s weight, race, hypertension history and uterine irritability fixed, smoking multiplies the odds of low birth weight by this much.”

The trouble is the phrase “holding fixed”. The clinical question usually looks like this instead:

If nobody in this population smoked, how much less low birth weight would there be?

Nothing is held fixed in that sentence. The people in the population differ in age, in weight, in race, and the question asks what happens to the average risk across all of them. A conditional OR does not answer it; it answers what happens within one particular combination of covariate values.

Two further things put the NNT out of reach of a conditional OR:

• An NNT is by definition the reciprocal of a risk difference, and an OR is not a risk difference and cannot be converted into one on its own — the same OR corresponds to wildly different risk differences depending on the baseline risk, which is the point of the OR-versus-RR figure on the logistic regression page.
• The OR is non-collapsible, so even if you convert it into some “average” quantity, what comes out is still not the population-level quantity you wanted. The section on non-collapsibility deals with this on its own.

G-computation in three steps

Standardisation applied to a regression model is called G-computation, and it has only three steps:

1. Fit a model, as usual. Here it is an ordinary logistic regression, with nothing special about it.
2. Make two copies of the entire dataset. In the first, set everyone’s exposure to “yes”; in the second, set everyone’s exposure to “no”. Change nothing else. Every mother therefore appears once in each copy, once as a smoker and once as a non-smoker.
3. Predict in each copy and average. The two averages are the two marginal risks; subtract them for the risk difference, divide them for the risk ratio.

Step 2 is the whole method. Leaving every other column alone means the study’s own covariate distribution is being used as the standard population, so what comes out is “this population, fully exposed versus fully unexposed” rather than some hypothetical population.

library(MASS)
data(birthwt, package = "MASS")

bw <- birthwt
bw$race_f  <- factor(bw$race,  levels = 1:3, labels = c("White", "Black", "Other"))
bw$smoke_f <- factor(bw$smoke, levels = 0:1, labels = c("No", "Yes"))

FORM <- low ~ smoke_f + age + lwt + race_f + ht + ui
fit <- glm(FORM, data = bw, family = binomial)

# Steps 2 and 3: two copies, predict in each, average
gcomp <- function(model, data) {
  d1 <- d0 <- data
  d1$smoke_f <- factor("Yes", levels = levels(data$smoke_f))
  d0$smoke_f <- factor("No",  levels = levels(data$smoke_f))
  r1 <- mean(predict(model, d1, type = "response"))
  r0 <- mean(predict(model, d0, type = "response"))
  c(risk1 = r1, risk0 = r0, RD = r1 - r0, RR = r1 / r0, NNT = 1 / (r1 - r0))
}
gcomp(fit, bw)

# The marginal OR comes from the two average risks, not from a coefficient
r <- gcomp(fit, bw)
(r["risk1"] / (1 - r["risk1"])) / (r["risk0"] / (1 - r["risk0"]))
exp(coef(fit))["smoke_fYes"]        # conditional OR -- a different number

# Confidence interval: resample individuals, rerun the whole procedure
set.seed(20260823)
boot <- replicate(2000, {
  ix <- sample(nrow(bw), replace = TRUE)
  d  <- bw[ix, ]
  gcomp(glm(FORM, data = d, family = binomial), d)
})
apply(boot, 1, quantile, c(0.025, 0.975))

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

bw = sm.datasets.get_rdataset("birthwt", "MASS").data
bw["race_f"]  = bw["race"].map({1: "White", 2: "Black", 3: "Other"})
bw["smoke_f"] = bw["smoke"].map({0: "No", 1: "Yes"})

FORM = "low ~ C(smoke_f) + age + lwt + C(race_f) + ht + ui"
fit = smf.glm(FORM, data=bw, family=sm.families.Binomial()).fit()

def gcomp(model, d):
    r1 = model.predict(d.assign(smoke_f="Yes")).mean()
    r0 = model.predict(d.assign(smoke_f="No")).mean()
    return dict(risk1=r1, risk0=r0, RD=r1 - r0, RR=r1 / r0, NNT=1 / (r1 - r0))

print(gcomp(fit, bw))

rng = np.random.default_rng(20260823)
boot = []
for _ in range(2000):
    d = bw.iloc[rng.integers(0, len(bw), len(bw))]
    f = smf.glm(FORM, data=d, family=sm.families.Binomial()).fit()
    g = gcomp(f, d)
    boot.append([g["RD"], g["RR"]])
print(np.percentile(np.array(boot), [2.5, 97.5], axis=0))

What comes out

The NNT reads as “for every this many mothers in this population who stop smoking, one fewer low-birth-weight baby”. There is no “with the other variables held fixed” in that sentence, because G-computation never held anything fixed.

Non-collapsibility: two ORs that disagree, and neither is wrong

One model, one set of people, one contrast — and two odds ratios:

• Conditional OR 2.794
• Marginal OR 2.419

The gap is 0.375, a ratio of 1.155. The conditional OR sits further from 1 — and it always sits on that side.

Why does this happen? Because the odds ratio is a non-collapsible measure: the average of the subgroup odds ratios is not the odds ratio of the average. Non-linearity alone is not the reason — a log link is non-linear too, and its rate ratio collapses perfectly well.

Working out the odds for a group of people at different risks and then averaging those odds is simply not the same operation as averaging their risks first and converting afterwards. The predicted risks this model assigns are spread widely — from 6.9% to 90.7% under the smoking scenario, and from 2.6% to 77.7% under the non-smoking one — and the wider that spread, the further apart the two ORs sit.

Why the confidence interval has to come from a bootstrap

The marginal risk difference is a non-linear function of every coefficient and every covariate, not one coefficient. The model output contains no standard error for it, and the standard-error column in summary() is not it either.

The standard approach is a non-parametric bootstrap: resample individuals, refit the model, and run the whole G-computation again, 2000 times, then take percentiles. “Run the whole thing again” is the load-bearing part — resampling without refitting amounts to pretending the coefficients are known.

So what happens if you cut corners? Here are three ways of getting an interval for the risk difference from the same data:

The NNT interval is lopsided — never write it as a plus-or-minus

An NNT is the reciprocal of a risk difference, and taking a reciprocal turns “close to zero” into “close to infinity”. The lower limit of the risk difference here is 4.2 percentage points, not far from zero, which pushes the upper limit of the NNT out to 21.7:

• Point estimate 5.2
• 95% CI 2.9 to 21.7

The lower limit sits only 2.3 below the point estimate; the upper limit sits 16.5 above it.

The predicted-probability figure is a standardised curve, not a partial plot

The right panel is easy to mistake for an ordinary partial dependence plot — “hold the other variables at their means and let weight vary”. It is not. Each point on those curves is computed by setting everyone’s weight to that one value, leaving every other column as observed, predicting, and averaging. It is G-computation applied a second time, now to mother’s weight.

The difference is not pedantry:

• Holding covariates at their means invents an “average patient” who may not exist: average age, average weight, and a race that is some decimal fraction. The more covariates there are, and the more of them are categorical, the more absurd that fictional person becomes.
• Standardisation invents nobody. It keeps the whole population’s heterogeneity and replaces only the one variable being varied.

The vertical distance between the two curves is the marginal risk difference at that level of weight. That distance shrinks in absolute terms as weight increases — from about 22.0 percentage points at the left-hand end to about 9.6 at the right.

Which one to report

Common mistakes

Further reading

• Logistic regression and the odds ratio — where the conditional OR comes from, and the first appearance of non-collapsibility
• Adjusted risk ratios — the other way around the OR: estimate an RR directly
• IPTW and propensity score weighting — another route to a marginal estimate, by weighting rather than by duplicating the data

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B2-08-marginal.R

Report a problem on this page