Comparison of IPTW and regression adjustment in causal inference

Question

Please see the reproducible R code in the end. The simulated data is from section 4.1 in the paper "ipw: An R Package for Inverse Probability Weighting", where we have measurements in 1000 individuals on a continuous confounder L, a dichotomous treatment A and a continuous outcome Y, which is the simplest case. The simulated data is generated as:

$L \sim N(10, 5); \ logitPr(A = 1) = −10 + L; \ Y = 10A + 0.5L + N(−10, 5)$.

I am comparing four methods to estimate the average treatment effect (in terms of mean squared error and coverage probability):

1. Regression of Y on A and L;
2. Parametric G-formula estimate with 95% bootstrap CI;
3. IPTW estimate using unstablized weight from "ipw" package and then estimate the ATE using svyglm function from "survey" package (this is the same method used in Section 4.1 in the paper "ipw: An R Package for Inverse Probability Weighting");
4. IPTW estimate by hand with 95% bootstrap CI.

The coverage probabilities are 0.95, 0.96, 0.58, 0.89 and MSE over 100 simulated datasets are 0.1813837, 0.1839691, 3.8638934, 3.6837547, respectively, for these four methods. So there are substantial bias using IPTW .

Q1: Although Methods 3 and 4 produce the same weight (see head(simdat) below), they do not produce the same ATE, is it because Method 3 with svyglm function is a parametric approach while method 4 is non-parametric ?

Q2: Is it possible for IPTW approach (e.g., methods 3 and 4) to match the bias with regression methods (methods 1 and 2) ?

Q3: The coverage probability for method 3 is very low (0.58) and has substantial bias compared to regression methods, does this mean we should not use this approach in practice ? (the authors mentioned in the paper that " we can fit the marginal structural model (MSM), estimating the marginal causal effect of A on Y , using a robust standard error estimate from the survey package") A broader question would be, under what circumstances IPTW approach (methods 3 and 4) is advantageous to regression methods ?

library(survey)
library(ipw)
library(boot)
set.seed(100)
Nsim <- 100
count <-MSE <- rep(0,4)

g.comp=function(data,indices)       # Use parametric G-formula to estimate the ATE
{
  dat=data[indices,]
  glm1  <- glm(y ~ l, data=dat[dat$a==1,])
  glm2  <- glm(y ~ l, data=dat[dat$a==0,])
  Y.1 = predict(glm1, newdata=data.frame(l=simdat$l), type="response")
  Y.0 = predict(glm2, newdata=data.frame(l=simdat$l), type="response")
  mean((Y.1) - mean(Y.0))
}

iptw.w = function(data,indices)     # Define the IPTW function to estimate the ATE
{
  dat=data[indices,]
  mean(dat$IPTW*as.numeric(dat$a==1)*dat$y) - mean(dat$IPTW*as.numeric(dat$a==0)*dat$y)
}

n <- 1000
for (i in 1:Nsim){
simdat <- data.frame(l = rnorm(n, 10, 5))
a.lin <- simdat$l - 10
pa <- exp(a.lin)/(1 + exp(a.lin))
simdat$a <- rbinom(n, 1, prob = pa)
simdat$y <- 10*simdat$a + 0.5*simdat$l + rnorm(n, -10, 5)

m_ps <- glm(a~ l,
            family = binomial(), data = simdat)

pr_score <- predict(m_ps, type = "response")

#inverse weights from the propensity score
simdat$IPTW <-simdat$a/pr_score + (1-simdat$a)/(1-pr_score )

### regression estimate of the coefficient
res1 <- lm(y ~ a+l,   data = simdat)
MSE[1] <- (coef(res1)[2] - 10)^2 + MSE[1]
int1 <-  confint.lm(res1)[2,]

### ATE estimate using g-formula, 95% bootstrap CI 
MSE[2] <- (g.comp(simdat,indices=1:nrow(simdat))  - 10)^2 + MSE[2]
boot.out=boot(simdat,g.comp,200)  
res2 <- boot.ci(boot.out,type="norm",conf=0.95) 
int2 <- c(res2$normal[2],res2$normal[3])

### IPW estimate, an example from section 4.1 in the paper "ipw: An R Package for Inverse Probability Weighting"
temp <- ipwpoint(exposure = a, family = "binomial", link = "logit",
                # numerator = ~ 1,
                 denominator = ~ l, data = simdat)
simdat$usw <- temp$ipw.weights
msm <- (svyglm(y ~ a, design = svydesign(~ 1, weights = ~ usw,
                                          data = simdat)))
int3 <- confint(msm)[2,]
MSE[3] <- (coef(msm)[2] - 10)^2 + MSE[3]

### IPTW estimate by hand
est <- mean(simdat$IPTW*as.numeric(simdat$a==1)*simdat$y) - mean(simdat$IPTW*as.numeric(simdat$a==0)*simdat$y)
MSE[4] <- (est  - 10)^2 + MSE[4]
boot.out4=boot(simdat,iptw.w,200)  
res4 <- boot.ci(boot.out4,type="norm",conf=0.95) 
int4 <- c(res4$normal[2],res4$normal[3])

if (int1[1]< 10 & int1[2] > 10) count[1] <- count[1] + 1
if (int2[1]< 10 & int2[2] > 10) count[2] <- count[2] + 1
if (int3[1]< 10 & int3[2] > 10) count[3] <- count[3] + 1
if (int4[1]< 10 & int4[2] > 10) count[4] <- count[4] + 1
}
> count/Nsim
[1] 0.95 0.96 0.58 0.89
> MSE/Nsim
[1] 0.1813837 0.1839691 3.8638934 3.6837547
#### IPTW by hand matches the unstablized weight from ipw package
> head(simdat)
          l a         y     IPTW      usw
1  4.955049 0 -5.393570 1.011515 1.011515
2 11.521077 0 -1.469263 5.066890 5.066890
3 17.228639 1  1.024713 1.001499 1.001499
4 15.842916 1 11.506791 1.005172 1.005172
5 10.047000 1  2.190818 1.917832 1.917832
6 19.214193 1 13.470528 1.000254 1.000254
```

Answer 1

This is an interesting question. I will say this simulation is biased to favor regression methods because the predictors explain a huge amount of variance in the treatment, meaning there is incomplete overlap and extremely severe imbalance, which are the situations in which IPW fares the worst. Take a look at the distribution of the covariate in the last dataset in the simulation:

There is clearly so little overlap and the true propensity scores for a large number of units (i.e., anyone outside the area of overlap) are near 0 or 1. This means that any method to estimate the ATE will necessarily require extrapolation. IPW cannot do extrapolation; indeed, this is its strength. Regression gets the answer right but only because the extrapolation is correct because the outcome model is linear; if the outcome model outside the area of overlap were not linear, the extrapolation would be incorrect and regression would be severely biased. Only using IPW allows you to see that the effect cannot be estimated without extrapolation.

In order for IPW to yield the correct answer, it needs to achieve balance. And here, it does not. Again, using the final dataset, the unadjusted standardized mean difference (SMD) is 2.18, and the weighted SMD is .66. SMDs need to be less than .1 to even be plausibly considered balanced. So, IPW failed to balance the covariates (because it is impossible for it to do so given the initial imbalance), and you should not attempt to estimate the effect from an imbalanced dataset.

Note that the phenomenon observed here is basically identical to those observed in this question and this question. If you want to designed an unbiased simulation to accurately compare the performance of two methods, you can't do so in a way that is unrealistic and automatically guarantees one of the methods to fail.

Let's go through your questions:

Q1) You use two different estimators. syvglm() uses what's known as a Hajek estimator, and your manual approach uses a Horvitz-Thompson estimator. You should always use Hajek estimators (even though in this case the MSE was large for the Hajek estimator). The way covariate balance is assessed and the way to incorporate an outcome model into an IPW weighted estimator is more in line with Hajek estimators. In more realistic scenarios, Hajek estimators perform better.

The HT estimator is, as you wrote, $$ \tau_{HT} = \frac{1}{n}\sum_{i=1}^n{A_i w_i Y_i} - \frac{1}{n}\sum_{i=1}^n{(1-A_i) w_i Y_i} $$ The Hajek estimator is $$ \tau_{Hajek} = \frac{\sum_{i=1}^n{A_i w_i Y_i}}{\sum_{i=1}^n{A_i w_i}} - \frac{\sum_{i=1}^n{(1-A_i) w_i Y_i}}{\sum_{i=1}^n{(1-A_i) w_i}} $$ The weighted least square estimator (what svyglm() does with a continuous outcome) is equal to the Hajek estimator (so you will get the same answer using lm()).

Q2 and Q3) These two questions are essentially the same, which is under what circumstances is it better to use regression vs. IPW? Rather than answer that explicitly, I'll describe the features of this simulation design that make it unfairly biased toward regression.

1. The outcome model is linear and correctly specified. When the outcome model is incorrectly specified but the treatment model is correctly specified, IPW will generally do better than regression in terms of bias.
2. There is severe imbalance. IPW does not perform well in cases of severe imbalance because the weights will be extreme. This doesn't affect regression much. Other methods of weighting are less affected by this; for example, overlap weighting, where the weights are $1-p$ for the treated units and $p$ for the control units (where $p$ is the propensity score), does not have this problem and can yield the right answer without extrapolating.
3. There is no effect modification. Regression assumes the relationship between the covariates and each potential outcome is the same. Doing g-computation by fitting separate models relaxes this requirement. IPW always averages over any possible effect modification, so it doesn't matter whether you think there is or isn't any. This flexibility comes at the cost of precision
4. There is very little overlap. I mentioned this before; the lack of overlap yields the significant imbalance observed and requires extrapolation to correctly estimate the effect. Regression does extrapolation and in this case gets it right because the outcome model is linear even beyond the area of overlap. IPW can only interpolate. Overlap weights restrict the estimation to the area of overlap, so they do not suffer from this problem.

In realistic settings, the true outcome model is likely not linear, so you can't rely on a linear model to get it right. You can use a machine learning method to do g-computation, but then you can't do valid inference. You can augment the g-computation estimator with a doubly-robust method that incorporates IPW using augmented IPW or targeted minimum loss-based estimation (TMLE), which de-biases the treatment effect estimate and provides a pathway for valid inference. Or, you can use a flexible machine learning method or optimization-based method to estimate the weights without requiring the outcome model to be correct. You can augment the weighted estimator using outcome regression using augmented IPW or TMLE as previously mentioned, which improves the precision and makes inference more valid. See my (cheeky) post here about when to use matching (e.g., weighting) vs. outcome regression models to estimate effects.