Doubly robust estimation with AIPW

causal inference

ecology tutorial

Augmented inverse-probability weighting (AIPW) is doubly robust: the causal effect is estimated consistently when either the propensity or the outcome model is right.

Author

Tidy Ecology

Published

2026-07-28

The two previous tutorials each rely on getting one model right. Inverse-probability weighting needs a correct model for the treatment; standardisation needs a correct model for the outcome. In practice you are rarely sure which of the two you can specify well. Augmented inverse-probability weighting, AIPW, combines them and buys a genuine insurance policy: it estimates the effect consistently if either the treatment model or the outcome model is correct, not necessarily both. That property is called double robustness.

The estimator

Take the standardisation estimate of the mean outcome under treatment level a, which is the average of the fitted values m_a, and add a weighted correction built from the residuals. For the treated arm the correction reweights the residual y - m_1 by the inverse propensity, and symmetrically for the control arm. The contrast of the two corrected means is the estimate. In one line per unit,

\[ \psi_i = \left[\frac{D_i\,(Y_i - m_1(X_i))}{\pi(X_i)} + m_1(X_i)\right] - \left[\frac{(1-D_i)\,(Y_i - m_0(X_i))}{1-\pi(X_i)} + m_0(X_i)\right], \]

and the estimated average treatment effect is the mean of psi over the sample. Here pi is the fitted propensity and m_1, m_0 are the outcome-model predictions with treatment set to one and to zero. If the outcome model is perfect, the residuals vanish in expectation and the estimate reduces to standardisation. If the propensity model is perfect, the weighted residual removes whatever bias a wrong outcome model leaves behind. Either way the bias cancels.

expit <- function(z) 1 / (1 + exp(-z))
pa <- c(a0=-0.5, a1=0.6, a2=0.4, a3=0.5, aq=0.35)          # true propensity
pb <- c(b0=-0.6, bD=log(2), b1=0.7, b2=0.5, b3=0.6, bq=0.5) # true outcome

gen <- function(n) {
  x1 <- rnorm(n); x2 <- rnorm(n); x3 <- rbinom(n, 1, 0.5)
  ps <- expit(pa["a0"] + pa["a1"]*x1 + pa["a2"]*x2 + pa["a3"]*x3 + pa["aq"]*x1^2)
  d  <- rbinom(n, 1, ps)
  py <- expit(pb["b0"] + pb["bD"]*d + pb["b1"]*x1 + pb["b2"]*x2 + pb["b3"]*x3 + pb["bq"]*x1^2)
  y  <- rbinom(n, 1, py)
  data.frame(y, d, x1, x2, x3)
}
true_ate <- function(dat) {
  l1 <- pb["b0"]+pb["bD"]*1+pb["b1"]*dat$x1+pb["b2"]*dat$x2+pb["b3"]*dat$x3+pb["bq"]*dat$x1^2
  l0 <- pb["b0"]+pb["bD"]*0+pb["b1"]*dat$x1+pb["b2"]*dat$x2+pb["b3"]*dat$x3+pb["bq"]*dat$x1^2
  as.numeric(mean(expit(l1)) - mean(expit(l0)))
}

Both the propensity and the outcome truly depend on x1 squared. A model with only linear terms is therefore misspecified, which lets us switch each nuisance model between correct and wrong by including or dropping the square.

gcomp <- function(om, dat)
  mean(predict(om, transform(dat, d=1), type="response")) -
  mean(predict(om, transform(dat, d=0), type="response"))

ipw <- function(pm, dat) {                          # normalised (Hajek) ATE
  e <- predict(pm, type="response"); d <- dat$d; y <- dat$y
  sum((d/e)*y)/sum(d/e) - sum(((1-d)/(1-e))*y)/sum((1-d)/(1-e))
}
aipw <- function(pm, om, dat) {
  e  <- predict(pm, type="response"); d <- dat$d; y <- dat$y
  m1 <- predict(om, transform(dat, d=1), type="response")
  m0 <- predict(om, transform(dat, d=0), type="response")
  psi <- (d*(y-m1)/e + m1) - ((1-d)*(y-m0)/(1-e) + m0)
  list(ate = mean(psi), se = sd(psi)/sqrt(length(psi)))
}
fit_all <- function(dat) list(
  om_c = glm(y ~ d + x1 + x2 + x3 + I(x1^2), binomial, dat),   # correct outcome
  om_w = glm(y ~ d + x1 + x2 + x3,           binomial, dat),   # wrong outcome
  pm_c = glm(d ~ x1 + x2 + x3 + I(x1^2),      binomial, dat),  # correct propensity
  pm_w = glm(d ~ x1 + x2 + x3,                binomial, dat))  # wrong propensity

One dataset

Fit all four nuisance models, then read off standardisation, weighting and AIPW under every combination of correct and wrong.

set.seed(167)
n <- 2000
dat <- gen(n)
f <- fit_all(dat)
e_c <- predict(f$pm_c, type="response")
TA <- true_ate(dat)

g_c <- gcomp(f$om_c, dat); g_w <- gcomp(f$om_w, dat)
i_c <- ipw(f$pm_c, dat);   i_w <- ipw(f$pm_w, dat)
a_cc <- aipw(f$pm_c, f$om_c, dat); a_cw <- aipw(f$pm_c, f$om_w, dat)
a_wc <- aipw(f$pm_w, f$om_c, dat); a_ww <- aipw(f$pm_w, f$om_w, dat)

data.frame(
  estimator = c("g-comp correct","g-comp wrong","IPW correct","IPW wrong",
                "AIPW ok/ok","AIPW ok/wrong","AIPW wrong/ok","AIPW wrong/wrong"),
  estimate  = round(c(g_c,g_w,i_c,i_w,a_cc$ate,a_cw$ate,a_wc$ate,a_ww$ate), 4),
  error     = round(c(g_c,g_w,i_c,i_w,a_cc$ate,a_cw$ate,a_wc$ate,a_ww$ate) - TA, 4))
#>          estimator estimate  error
#> 1   g-comp correct   0.1527 0.0108
#> 2     g-comp wrong   0.1902 0.0483
#> 3      IPW correct   0.1611 0.0192
#> 4        IPW wrong   0.1966 0.0547
#> 5       AIPW ok/ok   0.1564 0.0145
#> 6    AIPW ok/wrong   0.1580 0.0162
#> 7    AIPW wrong/ok   0.1554 0.0135
#> 8 AIPW wrong/wrong   0.1984 0.0566

The true effect on this sample is a risk difference of 0.142. The single-model estimators fail exactly where expected: standardisation with the wrong outcome model gives 0.190 and weighting with the wrong propensity gives 0.197, both well above the truth. AIPW with at least one correct model stays close: 0.156 with both correct, 0.158 with only the propensity correct, 0.155 with only the outcome correct. Only when both models are wrong does AIPW drift, to 0.198. One dataset is suggestive rather than conclusive, so repeat the whole thing many times.

set.seed(16700)
nsim <- 1000
cols <- c("gc_w","ipw_w","aipw_cc","aipw_cw","aipw_wc","aipw_ww")
E <- matrix(NA_real_, nsim, length(cols), dimnames = list(NULL, cols))
tr <- numeric(nsim); cov_cc <- cov_cw <- cov_wc <- 0L
for (s in 1:nsim) {
  ds <- gen(n); ta <- true_ate(ds); tr[s] <- ta
  ff <- fit_all(ds)
  E[s,"gc_w"]  <- gcomp(ff$om_w, ds)
  E[s,"ipw_w"] <- ipw(ff$pm_w, ds)
  acc <- aipw(ff$pm_c, ff$om_c, ds); E[s,"aipw_cc"] <- acc$ate
  acw <- aipw(ff$pm_c, ff$om_w, ds); E[s,"aipw_cw"] <- acw$ate
  awc <- aipw(ff$pm_w, ff$om_c, ds); E[s,"aipw_wc"] <- awc$ate
  E[s,"aipw_ww"] <- aipw(ff$pm_w, ff$om_w, ds)$ate
  if (ta >= acc$ate-1.96*acc$se && ta <= acc$ate+1.96*acc$se) cov_cc <- cov_cc+1L
  if (ta >= acw$ate-1.96*acw$se && ta <= acw$ate+1.96*acw$se) cov_cw <- cov_cw+1L
  if (ta >= awc$ate-1.96*awc$se && ta <= awc$ate+1.96*awc$se) cov_wc <- cov_wc+1L
}
mtruth <- mean(tr)
bias <- colMeans(E) - mtruth
round(rbind(mean = colMeans(E), bias = bias, sd = apply(E, 2, sd)), 4)
#>        gc_w  ipw_w aipw_cc aipw_cw aipw_wc aipw_ww
#> mean 0.1844 0.1875  0.1406  0.1407  0.1406  0.1905
#> bias 0.0443 0.0475  0.0005  0.0007  0.0005  0.0504
#> sd   0.0223 0.0226  0.0221  0.0222  0.0219  0.0227

Averaged over 1000 datasets, the mean true effect is 0.140. The two misspecified single-model estimators carry a bias of +0.044 (wrong outcome model) and +0.047 (wrong propensity). Every AIPW variant with one correct model is effectively unbiased: +0.0005, +0.0007 and +0.0005. The one AIPW that fails is the one with both models wrong, at +0.050.

lab <- c(gc_w="g-computation (wrong outcome)", ipw_w="IPW (wrong propensity)",
         aipw_ww="AIPW (both wrong)", aipw_cc="AIPW (both correct)",
         aipw_cw="AIPW (propensity ok, outcome wrong)", aipw_wc="AIPW (outcome ok, propensity wrong)")
ord <- c("gc_w","ipw_w","aipw_ww","aipw_cc","aipw_cw","aipw_wc")
b <- data.frame(key=ord, bias=bias[ord], se=apply(E,2,sd)[ord]/sqrt(nsim))
b$label <- factor(lab[ord], levels=lab[ord])
b$kind  <- ifelse(abs(b$bias) < 0.01, "consistent", "biased")
ggplot(b, aes(bias, label, colour=kind)) +
  geom_vline(xintercept=0, colour=te$label, linewidth=0.5) +
  geom_pointrange(aes(xmin=bias-1.96*se, xmax=bias+1.96*se), linewidth=0.9, size=0.6) +
  scale_colour_manual(values=c(consistent=te$mown, biased=te$abandoned), guide="none") +
  scale_x_continuous(limits=c(-0.01, 0.06)) +
  labs(x="mean bias in the estimated ATE (over 1000 datasets)", y=NULL,
       title="AIPW is unbiased when either nuisance model is correct",
       subtitle="Only when both the propensity and the outcome model are wrong does the bias return") +
  theme_te()

A horizontal interval plot. Three AIPW variants with at least one correct model sit on the zero line. Standardisation with a wrong outcome model, weighting with a wrong propensity, and AIPW with both wrong sit near a bias of 0.05. — Figure 1: Mean bias in the estimated average treatment effect over 1000 datasets, with a tight Monte Carlo interval. AIPW is centred on zero whenever at least one nuisance model is correct; the bias returns only when both are wrong, matching the two misspecified single-model estimators.

Why one correct model suffices

The augmentation term is a product of two errors: the treatment-model error and the outcome-model error. If the propensity is right, the inverse-probability weighting makes the outcome-model error average to zero, so a wrong outcome model does no harm. If the outcome model is right, its residuals already have mean zero within covariate strata, so a wrong propensity multiplies something that is already centred. The bias survives only when both pieces are wrong at once, because then neither factor is there to cancel it. That is the whole content of double robustness.

Uncertainty and coverage

Because AIPW is an average of the per-unit quantity psi, its standard error is just the standard deviation of psi divided by the square root of the sample size, the influence-function estimator. It performs well here.

c(both_correct = cov_cc, ps_ok_out_wrong = cov_cw, out_ok_ps_wrong = cov_wc) / nsim
#>    both_correct ps_ok_out_wrong out_ok_ps_wrong 
#>           0.938           0.955           0.941

Interval coverage is close to the nominal 95%: 93.8% with both models correct, 95.5% with only the propensity correct, and 94.1% with only the outcome correct. The plug-in standard error stays honest even when one nuisance model is misspecified.

dd <- rbind(data.frame(ate=E[,"aipw_cc"], est="AIPW (both correct)"),
            data.frame(ate=E[,"gc_w"],    est="g-computation (wrong outcome)"),
            data.frame(ate=E[,"ipw_w"],   est="IPW (wrong propensity)"))
dd$est <- factor(dd$est, levels=c("AIPW (both correct)","g-computation (wrong outcome)","IPW (wrong propensity)"))
ggplot(dd, aes(ate, fill=est, colour=est)) +
  geom_density(alpha=0.32, linewidth=0.7) +
  geom_vline(xintercept=mtruth, linetype="dashed", colour=te$ink) +
  annotate("text", x=mtruth, y=Inf, label="true ATE", vjust=1.6, hjust=1.08, size=3, colour=te$ink) +
  scale_fill_manual(values=c("AIPW (both correct)"=te$mown,
                             "g-computation (wrong outcome)"=te$abandoned,
                             "IPW (wrong propensity)"=te$grazed), name=NULL) +
  scale_colour_manual(values=c("AIPW (both correct)"=te$forest,
                               "g-computation (wrong outcome)"=te$abandoned,
                               "IPW (wrong propensity)"=te$low), name=NULL) +
  labs(x="estimated ATE (risk difference)", y="density",
       title="One correct model is enough for AIPW to centre on the truth",
       subtitle="Both misspecified single-model estimators are shifted away from the true ATE") +
  theme_te() + theme(legend.position="top")

Three overlaid density curves. The AIPW curve is centred on a dashed line at the true effect near 0.14. The wrong-outcome standardisation curve and the wrong-propensity weighting curve are separate peaks to the right, near 0.19. — Figure 2: Sampling distribution of the estimated average treatment effect over 1000 datasets. AIPW with both models correct is centred on the true effect, while both misspecified single-model estimators are shifted to the right.

What it does not buy you

Double robustness is insurance against one wrong model, not against both: the both-wrong AIPW above is as biased as the misspecified single-model fits, so the property is not a substitute for thinking about either model. It also does not repair poor overlap. The correct propensities in this sample run from 0.137 to 0.997, so a few units carry large weights, and near-violations of positivity inflate the variance whichever estimator you use. When both nuisance models are correct AIPW is also efficient, reaching the smallest possible large-sample variance, which is a second reason to prefer it. Fitting the nuisance models with flexible machine learning brings its own subtlety and needs sample splitting (cross-fitting), the route taken by targeted learning and double machine learning, beyond the scope here.

References

Robins J, Rotnitzky A, Zhao L 1994. Journal of the American Statistical Association 89(427):846-866 (10.2307/2290910).

Bang H, Robins J 2005. Biometrics 61(4):962-973 (10.1111/j.1541-0420.2005.00377.x).

Kang J, Schafer J 2007. Statistical Science 22(4):523-539 (10.1214/07-STS227).

Funk M, Westreich D, Wiesen C, Sturmer T, Davidian M, Cole S 2011. American Journal of Epidemiology 173(7):761-767 (10.1093/aje/kwq439).

Larsen A, Meng K, Kendall B 2019. Methods in Ecology and Evolution 10(7):924-934 (10.1111/2041-210X.13190).

Arif S, MacNeil M 2022. Ecosphere 13(5):e4009 (10.1002/ecs2.4009).

Hernan M, Robins J 2020. Causal Inference: What If. Chapman and Hall/CRC (ISBN 978-1-4200-7616-5).

--- title: "Doubly robust estimation with AIPW" description: "Augmented inverse-probability weighting (AIPW) is doubly robust: the causal effect is estimated consistently when either the propensity or the outcome model is right." date: "2026-07-28 10:00" categories: [causal inference, R, ecology tutorial] image: thumbnail.png image-alt: "Mean bias of six estimators over 1000 datasets, showing AIPW centred on zero whenever at least one nuisance model is correct and biased only when both are wrong." --- ```{r} #| label: setup #| include: false knitr::opts_chunk$set(collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE, dev = "png", dpi = 110, fig.width = 8, fig.height = 4.4) library(ggplot2) te <- list(ink="#16241d", body="#2c3a31", forest="#275139", label="#46604a", paper="#f5f4ee", line="#dad9ca", faint="#5d6b61", grazed="#cda23f", mown="#2f8f63", abandoned="#b5534e", high="#1d5b4e", low="#c9b458") theme_te <- function(base_size = 12) { theme_minimal(base_size = base_size) + theme(plot.background=element_rect(fill=te$paper,colour=NA), panel.background=element_rect(fill=te$paper,colour=NA), panel.grid.major=element_line(colour=te$line,linewidth=0.3), panel.grid.minor=element_blank(), axis.text=element_text(colour=te$faint), axis.title=element_text(colour=te$body), plot.title=element_text(colour=te$ink,face="bold"), plot.subtitle=element_text(colour=te$body), legend.text=element_text(colour=te$body), legend.title=element_text(colour=te$body), plot.margin=margin(10,12,10,10)) } expit <- function(z) 1 / (1 + exp(-z)) ``` The two previous tutorials each rely on getting one model right. [Inverse-probability weighting](../propensity-scores-and-ipw/) needs a correct model for the treatment; [standardisation](../g-computation-and-standardisation/) needs a correct model for the outcome. In practice you are rarely sure which of the two you can specify well. Augmented inverse-probability weighting, AIPW, combines them and buys a genuine insurance policy: it estimates the effect consistently if either the treatment model or the outcome model is correct, not necessarily both. That property is called double robustness. ## The estimator Take the standardisation estimate of the mean outcome under treatment level `a`, which is the average of the fitted values `m_a`, and add a weighted correction built from the residuals. For the treated arm the correction reweights the residual `y - m_1` by the inverse propensity, and symmetrically for the control arm. The contrast of the two corrected means is the estimate. In one line per unit, $$ \psi_i = \left[\frac{D_i\,(Y_i - m_1(X_i))}{\pi(X_i)} + m_1(X_i)\right] - \left[\frac{(1-D_i)\,(Y_i - m_0(X_i))}{1-\pi(X_i)} + m_0(X_i)\right], $$ and the estimated average treatment effect is the mean of `psi` over the sample. Here `pi` is the fitted propensity and `m_1`, `m_0` are the outcome-model predictions with treatment set to one and to zero. If the outcome model is perfect, the residuals vanish in expectation and the estimate reduces to standardisation. If the propensity model is perfect, the weighted residual removes whatever bias a wrong outcome model leaves behind. Either way the bias cancels. ```{r} #| label: dgp expit <- function(z) 1 / (1 + exp(-z)) pa <- c(a0=-0.5, a1=0.6, a2=0.4, a3=0.5, aq=0.35) # true propensity pb <- c(b0=-0.6, bD=log(2), b1=0.7, b2=0.5, b3=0.6, bq=0.5) # true outcome gen <- function(n) { x1 <- rnorm(n); x2 <- rnorm(n); x3 <- rbinom(n, 1, 0.5) ps <- expit(pa["a0"] + pa["a1"]*x1 + pa["a2"]*x2 + pa["a3"]*x3 + pa["aq"]*x1^2) d <- rbinom(n, 1, ps) py <- expit(pb["b0"] + pb["bD"]*d + pb["b1"]*x1 + pb["b2"]*x2 + pb["b3"]*x3 + pb["bq"]*x1^2) y <- rbinom(n, 1, py) data.frame(y, d, x1, x2, x3) } true_ate <- function(dat) { l1 <- pb["b0"]+pb["bD"]*1+pb["b1"]*dat$x1+pb["b2"]*dat$x2+pb["b3"]*dat$x3+pb["bq"]*dat$x1^2 l0 <- pb["b0"]+pb["bD"]*0+pb["b1"]*dat$x1+pb["b2"]*dat$x2+pb["b3"]*dat$x3+pb["bq"]*dat$x1^2 as.numeric(mean(expit(l1)) - mean(expit(l0))) } ``` Both the propensity and the outcome truly depend on `x1` squared. A model with only linear terms is therefore misspecified, which lets us switch each nuisance model between correct and wrong by including or dropping the square. ```{r} #| label: estimators gcomp <- function(om, dat) mean(predict(om, transform(dat, d=1), type="response")) - mean(predict(om, transform(dat, d=0), type="response")) ipw <- function(pm, dat) { # normalised (Hajek) ATE e <- predict(pm, type="response"); d <- dat$d; y <- dat$y sum((d/e)*y)/sum(d/e) - sum(((1-d)/(1-e))*y)/sum((1-d)/(1-e)) } aipw <- function(pm, om, dat) { e <- predict(pm, type="response"); d <- dat$d; y <- dat$y m1 <- predict(om, transform(dat, d=1), type="response") m0 <- predict(om, transform(dat, d=0), type="response") psi <- (d*(y-m1)/e + m1) - ((1-d)*(y-m0)/(1-e) + m0) list(ate = mean(psi), se = sd(psi)/sqrt(length(psi))) } fit_all <- function(dat) list( om_c = glm(y ~ d + x1 + x2 + x3 + I(x1^2), binomial, dat), # correct outcome om_w = glm(y ~ d + x1 + x2 + x3, binomial, dat), # wrong outcome pm_c = glm(d ~ x1 + x2 + x3 + I(x1^2), binomial, dat), # correct propensity pm_w = glm(d ~ x1 + x2 + x3, binomial, dat)) # wrong propensity ``` ## One dataset Fit all four nuisance models, then read off standardisation, weighting and AIPW under every combination of correct and wrong. ```{r} #| label: estimates set.seed(167) n <- 2000 dat <- gen(n) f <- fit_all(dat) e_c <- predict(f$pm_c, type="response") TA <- true_ate(dat) g_c <- gcomp(f$om_c, dat); g_w <- gcomp(f$om_w, dat) i_c <- ipw(f$pm_c, dat); i_w <- ipw(f$pm_w, dat) a_cc <- aipw(f$pm_c, f$om_c, dat); a_cw <- aipw(f$pm_c, f$om_w, dat) a_wc <- aipw(f$pm_w, f$om_c, dat); a_ww <- aipw(f$pm_w, f$om_w, dat) data.frame( estimator = c("g-comp correct","g-comp wrong","IPW correct","IPW wrong", "AIPW ok/ok","AIPW ok/wrong","AIPW wrong/ok","AIPW wrong/wrong"), estimate = round(c(g_c,g_w,i_c,i_w,a_cc$ate,a_cw$ate,a_wc$ate,a_ww$ate), 4), error = round(c(g_c,g_w,i_c,i_w,a_cc$ate,a_cw$ate,a_wc$ate,a_ww$ate) - TA, 4)) ``` The true effect on this sample is a risk difference of `r sprintf("%.3f", TA)`. The single-model estimators fail exactly where expected: standardisation with the wrong outcome model gives `r sprintf("%.3f", g_w)` and weighting with the wrong propensity gives `r sprintf("%.3f", i_w)`, both well above the truth. AIPW with at least one correct model stays close: `r sprintf("%.3f", a_cc$ate)` with both correct, `r sprintf("%.3f", a_cw$ate)` with only the propensity correct, `r sprintf("%.3f", a_wc$ate)` with only the outcome correct. Only when both models are wrong does AIPW drift, to `r sprintf("%.3f", a_ww$ate)`. One dataset is suggestive rather than conclusive, so repeat the whole thing many times. ```{r} #| label: mc set.seed(16700) nsim <- 1000 cols <- c("gc_w","ipw_w","aipw_cc","aipw_cw","aipw_wc","aipw_ww") E <- matrix(NA_real_, nsim, length(cols), dimnames = list(NULL, cols)) tr <- numeric(nsim); cov_cc <- cov_cw <- cov_wc <- 0L for (s in 1:nsim) { ds <- gen(n); ta <- true_ate(ds); tr[s] <- ta ff <- fit_all(ds) E[s,"gc_w"] <- gcomp(ff$om_w, ds) E[s,"ipw_w"] <- ipw(ff$pm_w, ds) acc <- aipw(ff$pm_c, ff$om_c, ds); E[s,"aipw_cc"] <- acc$ate acw <- aipw(ff$pm_c, ff$om_w, ds); E[s,"aipw_cw"] <- acw$ate awc <- aipw(ff$pm_w, ff$om_c, ds); E[s,"aipw_wc"] <- awc$ate E[s,"aipw_ww"] <- aipw(ff$pm_w, ff$om_w, ds)$ate if (ta >= acc$ate-1.96*acc$se && ta <= acc$ate+1.96*acc$se) cov_cc <- cov_cc+1L if (ta >= acw$ate-1.96*acw$se && ta <= acw$ate+1.96*acw$se) cov_cw <- cov_cw+1L if (ta >= awc$ate-1.96*awc$se && ta <= awc$ate+1.96*awc$se) cov_wc <- cov_wc+1L } mtruth <- mean(tr) bias <- colMeans(E) - mtruth round(rbind(mean = colMeans(E), bias = bias, sd = apply(E, 2, sd)), 4) ``` Averaged over `r nsim` datasets, the mean true effect is `r sprintf("%.3f", mtruth)`. The two misspecified single-model estimators carry a bias of `r sprintf("%+.3f", bias["gc_w"])` (wrong outcome model) and `r sprintf("%+.3f", bias["ipw_w"])` (wrong propensity). Every AIPW variant with one correct model is effectively unbiased: `r sprintf("%+.4f", bias["aipw_cc"])`, `r sprintf("%+.4f", bias["aipw_cw"])` and `r sprintf("%+.4f", bias["aipw_wc"])`. The one AIPW that fails is the one with both models wrong, at `r sprintf("%+.3f", bias["aipw_ww"])`. ```{r} #| label: fig-bias #| fig-cap: "Mean bias in the estimated average treatment effect over 1000 datasets, with a tight Monte Carlo interval. AIPW is centred on zero whenever at least one nuisance model is correct; the bias returns only when both are wrong, matching the two misspecified single-model estimators." #| fig-alt: "A horizontal interval plot. Three AIPW variants with at least one correct model sit on the zero line. Standardisation with a wrong outcome model, weighting with a wrong propensity, and AIPW with both wrong sit near a bias of 0.05." lab <- c(gc_w="g-computation (wrong outcome)", ipw_w="IPW (wrong propensity)", aipw_ww="AIPW (both wrong)", aipw_cc="AIPW (both correct)", aipw_cw="AIPW (propensity ok, outcome wrong)", aipw_wc="AIPW (outcome ok, propensity wrong)") ord <- c("gc_w","ipw_w","aipw_ww","aipw_cc","aipw_cw","aipw_wc") b <- data.frame(key=ord, bias=bias[ord], se=apply(E,2,sd)[ord]/sqrt(nsim)) b$label <- factor(lab[ord], levels=lab[ord]) b$kind <- ifelse(abs(b$bias) < 0.01, "consistent", "biased") ggplot(b, aes(bias, label, colour=kind)) + geom_vline(xintercept=0, colour=te$label, linewidth=0.5) + geom_pointrange(aes(xmin=bias-1.96*se, xmax=bias+1.96*se), linewidth=0.9, size=0.6) + scale_colour_manual(values=c(consistent=te$mown, biased=te$abandoned), guide="none") + scale_x_continuous(limits=c(-0.01, 0.06)) + labs(x="mean bias in the estimated ATE (over 1000 datasets)", y=NULL, title="AIPW is unbiased when either nuisance model is correct", subtitle="Only when both the propensity and the outcome model are wrong does the bias return") + theme_te() ``` ## Why one correct model suffices The augmentation term is a product of two errors: the treatment-model error and the outcome-model error. If the propensity is right, the inverse-probability weighting makes the outcome-model error average to zero, so a wrong outcome model does no harm. If the outcome model is right, its residuals already have mean zero within covariate strata, so a wrong propensity multiplies something that is already centred. The bias survives only when both pieces are wrong at once, because then neither factor is there to cancel it. That is the whole content of double robustness. ## Uncertainty and coverage Because AIPW is an average of the per-unit quantity `psi`, its standard error is just the standard deviation of `psi` divided by the square root of the sample size, the influence-function estimator. It performs well here. ```{r} #| label: coverage c(both_correct = cov_cc, ps_ok_out_wrong = cov_cw, out_ok_ps_wrong = cov_wc) / nsim ``` Interval coverage is close to the nominal 95%: `r sprintf("%.1f%%", 100*cov_cc/nsim)` with both models correct, `r sprintf("%.1f%%", 100*cov_cw/nsim)` with only the propensity correct, and `r sprintf("%.1f%%", 100*cov_wc/nsim)` with only the outcome correct. The plug-in standard error stays honest even when one nuisance model is misspecified. ```{r} #| label: fig-dist #| fig-cap: "Sampling distribution of the estimated average treatment effect over 1000 datasets. AIPW with both models correct is centred on the true effect, while both misspecified single-model estimators are shifted to the right." #| fig-alt: "Three overlaid density curves. The AIPW curve is centred on a dashed line at the true effect near 0.14. The wrong-outcome standardisation curve and the wrong-propensity weighting curve are separate peaks to the right, near 0.19." dd <- rbind(data.frame(ate=E[,"aipw_cc"], est="AIPW (both correct)"), data.frame(ate=E[,"gc_w"], est="g-computation (wrong outcome)"), data.frame(ate=E[,"ipw_w"], est="IPW (wrong propensity)")) dd$est <- factor(dd$est, levels=c("AIPW (both correct)","g-computation (wrong outcome)","IPW (wrong propensity)")) ggplot(dd, aes(ate, fill=est, colour=est)) + geom_density(alpha=0.32, linewidth=0.7) + geom_vline(xintercept=mtruth, linetype="dashed", colour=te$ink) + annotate("text", x=mtruth, y=Inf, label="true ATE", vjust=1.6, hjust=1.08, size=3, colour=te$ink) + scale_fill_manual(values=c("AIPW (both correct)"=te$mown, "g-computation (wrong outcome)"=te$abandoned, "IPW (wrong propensity)"=te$grazed), name=NULL) + scale_colour_manual(values=c("AIPW (both correct)"=te$forest, "g-computation (wrong outcome)"=te$abandoned, "IPW (wrong propensity)"=te$low), name=NULL) + labs(x="estimated ATE (risk difference)", y="density", title="One correct model is enough for AIPW to centre on the truth", subtitle="Both misspecified single-model estimators are shifted away from the true ATE") + theme_te() + theme(legend.position="top") ``` ## What it does not buy you Double robustness is insurance against one wrong model, not against both: the both-wrong AIPW above is as biased as the misspecified single-model fits, so the property is not a substitute for thinking about either model. It also does not repair poor overlap. The correct propensities in this sample run from `r sprintf("%.3f", min(e_c))` to `r sprintf("%.3f", max(e_c))`, so a few units carry large weights, and near-violations of positivity inflate the variance whichever estimator you use. When both nuisance models are correct AIPW is also efficient, reaching the smallest possible large-sample variance, which is a second reason to prefer it. Fitting the nuisance models with flexible machine learning brings its own subtlety and needs sample splitting (cross-fitting), the route taken by targeted learning and double machine learning, beyond the scope here. ## References Robins J, Rotnitzky A, Zhao L 1994. Journal of the American Statistical Association 89(427):846-866 (10.2307/2290910). Bang H, Robins J 2005. Biometrics 61(4):962-973 (10.1111/j.1541-0420.2005.00377.x). Kang J, Schafer J 2007. Statistical Science 22(4):523-539 (10.1214/07-STS227). Funk M, Westreich D, Wiesen C, Sturmer T, Davidian M, Cole S 2011. American Journal of Epidemiology 173(7):761-767 (10.1093/aje/kwq439). Larsen A, Meng K, Kendall B 2019. Methods in Ecology and Evolution 10(7):924-934 (10.1111/2041-210X.13190). Arif S, MacNeil M 2022. Ecosphere 13(5):e4009 (10.1002/ecs2.4009). Hernan M, Robins J 2020. Causal Inference: What If. Chapman and Hall/CRC (ISBN 978-1-4200-7616-5). ## Related tutorials - [G-computation and standardisation](../g-computation-and-standardisation/) - [Propensity scores and inverse-probability weighting](../propensity-scores-and-ipw/) - [Propensity score matching](../propensity-score-matching/) - [Sensitivity to unmeasured confounding](../sensitivity-to-unmeasured-confounding/)