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):
- Regression of Y on A and L;
- Parametric G-formula estimate with 95% bootstrap CI;
- 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");
- 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
```