Bang and Robins doubly-robust estimator biased and with large variance?

Question

In their 2005 paper (also see the correction here) Bang and Robins describe a doubly robust estimator of the average treatment effect.

In short, the procedure is:

1. Estimate inverse probability of treatment weights (IPTW). These are $1 / Pr[A=1|L]$ for those with $A = 1$ and $-1 / (1 - Pr[A=1|L])$ for those with $A=0$, for treatment $A$ and vector of confounders $L$.
2. Then, include the IPTW in the conditional mean outcome model as a 'clever covariate': $E[Y|A,L,R]$, where $R$ is the weights defined above.
3. Finally, standardise over the confounder distribution of the sample.

In the words of the Bang and Robins correction:

...see that we must add to the regression the inverse probability of treatment weighted (IPTW) covariate, which is the (estimated) inverse of the PS for treated subjects(Δ= 1)and the inverse of the negative of “1 minus the PS” for untreated subjects(Δ= 0). Other choices can result in inconsistent estimation of the average treatment effect.

I have tried to simulate simple data and test out this estimator under correct specification of the treatment and outcome models. The results I'm getting are confusing.

Here is the simulation code:

library(tidyverse)


dr_sim <- function(SEED){

## Simulate data 

set.seed(SEED)
  
L <- rnorm(1000)
A <- rbinom(1000,1,plogis(-0.5 + 0.25*L))
Y <- 0.5*L + rnorm(1000) # No causal effect of A, confounding by L 

d <- tibble(L,A,Y)

## fit IPTW MSM

ip_mod <- glm(A ~ L, family = binomial, data = d)

ipt_weight <- ifelse(d$A == 1, 
                     1 / predict(ip_mod, type = "response"), 
                     1 / (1 - predict(ip_mod, type = "response")))

msm <- lm(Y ~ A, data = d, weights = ipt_weight)

msm_est <- coef(msm)["A"] # correct

## parametric g-formula

gmod <- lm(Y ~ A + L, data = d)

d1 <- d0 <- d

d0$A <- 0

E_0 <- mean(predict(gmod, newdata = d0))

d1$A <- 1

E_1 <- mean(predict(gmod, newdata = d1))

gform_est <- E_1 - E_0 # correct


## Bang & Robins doubly robust estimator

dr_weight <- ifelse(d$A == 1, 
                    1 / predict(ip_mod, type = "response"),
                    1 / -(1 - predict(ip_mod, type = "response")))

dr_mod <- lm(Y ~ A + L + dr_weight, data = d)

d1 <- d0 <- d

d0$A <- 0

E_0 <- mean(predict(dr_mod, newdata = d0))

d1$A <- 1

E_1 <- mean(predict(dr_mod, newdata = d1))

dr_est <- E_1 - E_0

out <- tibble(msm_est, gform_est, dr_est)

out
}


## Simulate

results <- tibble(seed = 1:1000) |> 
  mutate(out = map(seed, dr_sim)) |> 
  unnest(out)


## results

tibble(
  estimator = c("IPTW","g-formula","doubly-robust"),
  mean_estimate = c(mean(results$msm_est), mean(results$gform_est), mean(results$dr_est)),
  sd_estimate = c(sd(results$msm_est), sd(results$gform_est), sd(results$dr_est)))

And here are the results:

# A tibble: 3 × 3
  estimator     mean_estimate sd_estimate
  <chr>                 <dbl>       <dbl>
1 IPTW               0.000271      0.0628
2 g-formula          0.000285      0.0626
3 doubly-robust      0.0593        1.19  

The IPTW marginal structural model and parametric g-formula estimators are unbiased. The doubly robust robust estimator is not. The magnitude of bias and the extent of the variance seems to depend on the strength of the confounder-treatment relationship. When the confounder strongly affects the treatment, the doubly-robust estimator is unbiased and has low variance (lower than the IPTW estimator), but when it is weak, as above, the variance and bias can be very large.

Can anyone shed any light on this? Am I incorrectly implementing the estimator? Or is there something else I'm missing?

Answer 1

This is a slight error in how you programmed the DR estimator. Although dr_weight is the clever covariate, it is actually a function of the treatment. That means when you do g-computation on the DR outcome model, you need to set the value of the treatment to the specified value not just in the dataset but in the clever covariate, too, rather than leaving the clever covariate as is. Here is what the estimator should look like:

p <- predict(ip_mod, type = "response")
  
dr_mod <- lm(Y ~ A + L + I(A/p - (1-A)/(1-p)), data = d)
  
E_0 <- mean(predict(dr_mod, newdata = d0))
  
E_1 <- mean(predict(dr_mod, newdata = d1))
  
dr_est <- E_1 - E_0

That is, instead of defining the clever covariate outside the model and fixing its values, we let it be a function of the treatment and propensity score, i.e., I(A/p - (1 - A)/(1 - p)). This way, in the g-computation step, the value of A that is set in newdata is inserted into the clever covariate as well.