# Comparison of IPTW and regression adjustment in causal inference

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
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
$$$$


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.

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()`).
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.