R Example 9: Inverse propensity score weighting

We still use the lalonde data from the MatchIt package and use the propensity score model that we found out in R example 7.

library("MatchIt")
data("lalonde")

model <- glm(treat ~ . , data = lalonde[, -9], family = "binomial")
eps <- predict(model, type = "response")

lalonde

1. Esimate the inverse probability weights

• Visualize the weights

The figures are not very informative as there are units with extremely large weights. We can perfrom log transformation to make the histogram more informative.

## Calculate the raw weights (before normalization)
## Weights to estimate ATE
n.treated <- sum(lalonde$treat == 1)
n.control <- sum(lalonde$treat == 0)
weights <- ifelse(lalonde$treat == 1, 1/eps, 1/(1 - eps))

## Check the weights histogram
library(ggplot2)
temp.data <- data.frame(weights = weights, treated = as.factor(lalonde$treat))
ggplot(temp.data, aes(x = weights, fill = treated, color = treated)) + 
  geom_histogram(alpha = 0.5, position = "identity") + 
  xlab("Weights") 

## We can perform log transformation to make the histogram more informative
ggplot(temp.data, aes(x = weights, fill = treated, color = treated)) + 
  geom_histogram(alpha = 0.5, position = "identity") +  scale_x_log10() + 
  xlab("Weights") 

• Visualize the estimated propensity scores

temp.data.eps <- data.frame(eps = eps, treated = as.factor(lalonde$treat))
ggplot(temp.data.eps, aes(x = eps, fill = treated, color = treated)) + 
  geom_histogram(alpha = 0.5, position = "identity") + xlim(c(0, 1)) +
  ggtitle("Histogram of eps before trimming")

We likely will do some trimming to make sure overlapping assumption is satisfied and probably will also improve covariates balancing.

• Check covariate balancing.

We can check covariate balancing from the love plot after inverse probability weighting. Here, we need to write our own function to draw the love plot.

## Need to change race (categorical) into indicators (numerical)
lalonde$black <- lalonde$race == "black"
lalonde$hispan <- lalonde$race == "hispan"
lalonde$white <- lalonde$race == "white"

## Draw love plot
love.plot = function(cov, treat,  ## cov is the matrix of covariates and treat is a vector of treatment assignment
                     weights = rep(1, length(treat)),
                     plot = F) 
{
    
    ## mean with normalized weights \sum w_i x_i / (\sum w_i)
  treat.means <- colSums(cov[treat == 1,] * weights[treat == 1])/sum(weights[treat == 1])
  treat.var <- colSums(t(t(cov[treat == 1,]) - treat.means)^2 *
                          weights[treat == 1])/sum(weights[treat == 1])
  
  control.means <- colSums(cov[treat == 0,] * weights[treat == 0])/sum(weights[treat == 0])
  control.var <- colSums(t(t(cov[treat == 0,]) - control.means)^2 *
                          weights[treat == 0])/sum(weights[treat == 0])
  
  ## the standardized mean differences for every covariate
  smd <- (treat.means - control.means)/sqrt((treat.var + control.var)/2)
  names(smd) <- colnames(cov)
  
  if (plot == T) {
    plot.data <- data.frame(smd = smd, covariates = names(smd))
    range <- max(abs(smd))
    ggplot(plot.data) + geom_point(aes(x = as.numeric(smd), y = covariates)) +
      geom_vline(xintercept = 0) + xlim(-range, range) +
      labs(x = 'Standardized Difference in Means')
  }
  return(smd)
}

raw.smd <- love.plot(lalonde[, c(2:3, 5:8, 10:12)], lalonde$treat)
weighted.smd <- love.plot(lalonde[, c(2:3, 5:8, 10:12)], lalonde$treat, weights = weights)


plot.data <- data.frame(smd = c(raw.smd, weighted.smd), 
                        covariates = c(names(raw.smd), names(weighted.smd)),
                        category = c(rep("Original", length(raw.smd)), rep("IPW", length(weighted.smd))))
range <- max(abs(plot.data$smd))

ggplot(plot.data) + geom_point(aes(x = as.numeric(smd), y = covariates, color = category)) +
      geom_vline(xintercept = c(-0.1, -0.05, 0, 0.05, 0.1),
                 linetype = c("solid", "dashed", "solid", "dashed", "solid")) + 
      xlim(-range, range) +
      labs(x = 'Standardized Difference in Means')

Covariate balancing is greatly improved after IPW. There is still room for further improvement.

2. Trimming

2.1 Check weights and covariates balancing after trimming

• Perform trimming and check the histogram of weights

rm.idx <- which(eps < 0.1 | eps > 0.9)

temp.data.trimmed <- temp.data[-rm.idx, ]
lalonde.trimmed <- lalonde[-rm.idx, ]

## Check the histogram of weights again

ggplot(temp.data.trimmed, aes(x = weights, fill = treated, color = treated)) + 
  geom_histogram(alpha = 0.5, position = "identity") +  scale_x_log10() + 
  xlab("Weights") + ggtitle("Histogram of weights after trimming")

• Check covariate balancing after trimming.

raw.smd <- love.plot(lalonde.trimmed[, c(2:3, 5:8, 10:12)], lalonde.trimmed$treat)
weighted.smd <- love.plot(lalonde.trimmed[, c(2:3, 5:8, 10:12)], lalonde.trimmed$treat, 
                          weights = weights[-rm.idx])
plot.data <- data.frame(smd = c(raw.smd, weighted.smd), 
                        covariates = c(names(raw.smd), names(weighted.smd)),
                        category = c(rep("Original", length(raw.smd)), rep("IPW", length(weighted.smd))))
range <- max(abs(plot.data$smd))
ggplot(plot.data) + geom_point(aes(x = as.numeric(smd), y = covariates, color = category)) +
      geom_vline(xintercept = c(-0.1, -0.05, 0, 0.05, 0.1),
                 linetype = c("solid", "dashed", "solid", "dashed", "solid")) + 
      xlim(-range, range) +
      labs(x = 'Standardized Difference in Means')

Though covariate balanicng of the raw data gets better, covariates balancing with IPW does not improve.

2.2 Refit the propensity score model after trimming

We do this because the covariate balancing the previous step is not satisfying.

## Refit the propensity score model
model.trimmed <- glm(treat ~ . , data = lalonde.trimmed[, -(9:12)], family = "binomial")
eps.trimmed <- predict(model.trimmed, type = "response")
weights.trimmed <- ifelse(lalonde.trimmed$treat == 1, 1/eps.trimmed, 1/(1 - eps.trimmed))
temp.data.trimmed$weights <- weights.trimmed

## Check the histogram of weights again

ggplot(temp.data.trimmed, aes(x = weights, fill = treated, color = treated)) + 
  geom_histogram(alpha = 0.5, position = "identity") +  scale_x_log10() + 
  xlab("Weights") + ggtitle("Histogram of weights after trimming")

• Check covariate balancing again.

raw.smd <- love.plot(lalonde.trimmed[, c(2:3, 5:8, 10:12)], lalonde.trimmed$treat)
weighted.smd <- love.plot(lalonde.trimmed[, c(2:3, 5:8, 10:12)], lalonde.trimmed$treat, 
                          weights = weights.trimmed)
plot.data <- data.frame(smd = c(raw.smd, weighted.smd), 
                        covariates = c(names(raw.smd), names(weighted.smd)),
                        category = c(rep("Original", length(raw.smd)), rep("IPW", length(weighted.smd))))
range <- max(abs(plot.data$smd))
ggplot(plot.data) + geom_point(aes(x = as.numeric(smd), y = covariates, color = category)) +
      geom_vline(xintercept = c(-0.1, -0.05, 0, 0.05, 0.1),
                 linetype = c("solid", "dashed", "solid", "dashed", "solid")) + 
      xlim(-range, range) +
      labs(x = 'Standardized Difference in Means')

3. Estimate the ATE using IPW

• Estimate ATE by weighted least square regression

lm.result <- lm(re78 ~ treat, weights = weights.trimmed, data = lalonde.trimmed)
summary(lm.result)

Call:
lm(formula = re78 ~ treat, data = lalonde.trimmed, weights = weights.trimmed)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
-16478  -7345  -2741   3852  61182 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5402.7      537.4  10.053   <2e-16 ***
treat         1155.1      756.9   1.526    0.128    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9888 on 339 degrees of freedom
Multiple R-squared:  0.006823,  Adjusted R-squared:  0.003894 
F-statistic: 2.329 on 1 and 339 DF,  p-value: 0.1279

• Statistical inference treating the weights as fixed.

We use the Sandwich estimator allowing heteroscedastic noise levels between the treated and control grousp

library(sandwich)
tau_hat <- lm.result$coefficients[2]
SE <- sqrt(diag(vcovHC(lm.result, type = "HC2")))[2]

## get the 95% CI
result <- c(tau_hat, SE, c(tau_hat- 1.96 * SE, tau_hat + 1.96 * SE))
names(result) <- c("est", "sd", "CI_lower", "CI_upper")
result

      est        sd  CI_lower  CI_upper 
1155.0913  827.0484 -465.9235 2776.1062 

• Statistical inference using bootstrap

IPW_estimator <- function(W, Y, X) {
  ## Estimate propensity score
  model <- glm(W ~ X , family = "binomial")
  eps <- predict(model, type = "response")
  
  ## Calculate the weights
  weights <- ifelse(W == 1, 1/eps, 1/(1 - eps))
  
  ## Calculate weighted mean difference between treated and control group
  est <- lm(Y ~ W, weights = weights)$coef[2]
  return(est)
}

X <- model.matrix(model.trimmed)
Y <- lalonde.trimmed$re78
W <- lalonde.trimmed$treat

IPW_bootstrap <- function(W, Y, X, n.boot = 200){
  est <- IPW_estimator(W, Y, X)
  IPWboot <- sapply(1:n.boot, function(i) {
    id.boot <- sample(1:length(W), replace = T)
    IPW_estimator(W[id.boot], Y[id.boot], X[id.boot, ])
  })
  return(c(est, sd(IPWboot)))
}

SE_boostrap <- IPW_bootstrap(W, Y, X, 5000)[2]
result_bootstrap <- c(tau_hat, SE_boostrap, 
                      c(tau_hat- 1.96 * SE_boostrap, tau_hat + 1.96 * SE_boostrap))
names(result_bootstrap) <- c("est", "sd (bootstrap)", "CI_lower (bootstrap)", "CI_upper (bootstrap)")
result_bootstrap

                 est       sd (bootstrap) CI_lower (bootstrap) CI_upper (bootstrap) 
           1155.0913             869.8960            -549.9049            2860.0875