Recall that the stabilized IPW estimator for the ATE is \[ \hat{\rho}_{IPW} = \frac{\sum_{i=1}^n \frac{D_i}{\hat{e}(X_i)}Y_i}{\sum_{i=1}^n \frac{D_i}{\hat{e}(X_i)}} - \frac{\sum_{i=1}^n \frac{1-D_i}{1-\hat{e}(X_i)}Y_i}{\sum_{i=1}^n \frac{1-D_i}{1-\hat{e}(X_i)}}, \] and that this can be implementing with weighted least squares regression in which you regress \(Y\) on \(D\), using the inverse of the estimate propensity score \(\hat{e}(X_i)\) as the weights for the treated and the inverse of \(1-\hat{e}(X_i)\) as the weights for the control.
Inference for \(\hat{\rho}_{IPW}\) should account for the sampling variation in the estimated propensity scores, \(\hat{e}(\cdot)\). In the Robust Inference I lecture notes, we showed how to derived the variance when \(\hat{e}(\cdot)\) is estimated using an \(M\)-estimator. Here we will use the bootstrap instead.
We will simulate a population with confounding in \(X\). We do it in a way that makes it so that simple control specifications for \(X\) in a regression would be inadequate, and so the specification challenge is a little harder than usual:
# Make population data
set.seed(11)
N.pop <- 10000
index <- 1:N.pop
X <- .5*exp(rnorm(N.pop))
Y0 <- 50-.05*X^2+(1+.5*(max(X)-X))*rnorm(N.pop)
Y1 <- -1.5 + Y0 + 10*X + .2*X^2+ .25*X*rnorm(N.pop)
e.X <- (1+exp(1-.2*X))^(-1)
D <- rbinom(N.pop, 1, e.X)
Y <- D*Y1 + (1-D)*Y0
rho <- mean(Y1-Y0)
plot(X, Y0, col="blue", pch=19,
cex=.25, ylim=range(c(Y1,Y0)),
main="Potential outcomes")
points(X, Y1, col="red", pch=19, cex=.25)
plot(X, e.X, pch=19, main="Propensity scores")
plot(X, Y, type="n", ,
main="Realized outcomes")
points(X[D==0], Y[D==0], col="blue", pch=19, cex=.25)
points(X[D==1], Y[D==1], col="red", pch=19, cex=.25)
rho
## [1] 7.124273
summary(lm(Y~D))[[4]][,1]
## (Intercept) D
## 49.800748 9.237881
summary(lm(Y~D+X))[[4]][,1]
## (Intercept) D X
## 45.379997 7.790128 5.881640
summary(lm(Y~D*I(X-mean(X))))[[4]][,1]
## (Intercept) D I(X - mean(X)) D:I(X - mean(X))
## 49.7601797 7.3606239 -0.5394921 11.7581545
pop.data <- data.frame(index, Y, D, X)
Next we draw a sample and estimate the effect using simple regression strategies and then IPW:
# Draw a sample:
n.samp <- 1000
samp.i <- sample(index, n.samp)
samp.data <- pop.data[samp.i,]
plot(samp.data$X, samp.data$Y, type="n")
points(samp.data$X[D==0], samp.data$Y[D==0], col="blue")
points(samp.data$X[D==1], samp.data$Y[D==1], col="red")