R Example 7: Subclassification

We still use the lalonde data, perform trimming as we did before. We can also remove units whose eps are outside of the range the control (for a treated unit) or outside of the range of the treated (for a control units)

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

model <- glm(treat ~ re74 + race + married + I(re74^2) + re74:race, data = lalonde, family = "binomial")
eps <- predict(model, type = "response")

rm.idx <- which(eps < 0.1 | eps > 0.9)
idx.treated <- which((lalonde$treat == 1) & (eps > max(eps[lalonde$treat == 0])))
idx.control <- which((lalonde$treat == 0) & (eps < min(eps[lalonde$treat == 1])))
rm.idx <- union(rm.idx, c(idx.treated, idx.control))

lalonde <- lalonde[-rm.idx, ]
eps <- eps[-rm.idx]
lalonde

• Now, we perform sequencing splitting to find the subclasses

## First, create a data frame to store the current grouping information
temp <- data.frame(e = eps, treat = lalonde$treat,
                           b = 1)

t.max <-  1.96
# check whether t. stat is above t.max for first iteration
t.stat <- t.test(x = temp$e[temp$treat == 1], y =  temp$e[temp$treat == 0], var.equal=T)$statistic
condition = t.stat > t.max
# minimum n.t or n.c in each block
size <- 3
size.new <- 20

# continue until all t statistics are below t.max
set.seed(0)
while(condition)
{
  # calculate size in each group
  # we want to not split a block if it is too small
  b <- max(temp$b)
  ignore <- sapply(1:b, function(j) {
     n.t <- sum(temp$treat == 1 & temp$b == j)
     n.c <- sum(temp$treat == 0 & temp$b == j)
     return(n.t < size | n.c < size | (n.t + n.c < size.new * 2))
  })
  
  # split unbalanced blocks into more blocks
  split <- which((abs(t.stat) > t.max) & (!ignore))
  
  if(length(split) == 0)
    break

  
  ## we need to keep a current copy of the block information and which block to ignore as block assignments are going to change later
  b.current <- temp$b
  
  for(j in split)
  {
    
    cutoff <- median(temp$e[b.current == j])
    ## We split units into two new blocks
    ## extract the index of units belonging to each new stratum
    idx.s <- which(b.current == j & temp$e < cutoff)
    idx.l <- which(b.current == j & temp$e > cutoff)
    ## randomly put half of the ties into one category
    idx.e <- sample(which(temp$e == cutoff & b.current == j))
    n.tie <- length(idx.e)
    if (n.tie >= 1) {
      if (n.tie > 1)
        idx.s <- c(idx.s, idx.e[1:round(n.tie/2)])
      idx.l <- c(idx.l, idx.e[(round(n.tie/2)+ 1):n.tie])
    }
      
    ## we split only when new stratum has at least size number of control/treated units
    if (sum(temp$treat[idx.s]==1) > size && sum(temp$treat[idx.s]==0) > size && 
        sum(temp$treat[idx.l]==1) > size && sum(temp$treat[idx.l]==0) > size) {
      # anything above the current will have to be moved up 1
      temp$b[b.current > j] <- temp$b[b.current > j] + 1
      ## the upper new stratum will also have the block idx added by 1
      temp$b[idx.l] <- temp$b[idx.l] + 1
    }
    ## We don't do anything if we do not want to split
  }
  
   # calculate t statistic for each block
  b <- max(temp$b)
  t.stat <- sapply(1:b, function(j) {
    t.test(x = temp$e[temp$treat == 1 & temp$b == j], 
                        y = temp$e[temp$treat == 0 & temp$b == j], var.equal=T)$statistic
  })
  
  ## Update condition
  # check whether ANY blocks are above t.max
  # AND are not too small
  condition <- any(abs(t.stat) > t.max)
}

lalonde$blocks <- temp$b
table(lalonde[, c("treat", "blocks")])

     blocks
treat   1   2   3   4   5
    0  72  36  18   7  67
    1   7  12   8  19 130

## check the range of estimated propensity scores
sapply(1:max(lalonde$blocks), function(j) range(temp$e[temp$b == j]))

          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.1011357 0.1451386 0.2376577 0.4672575 0.5555612
[2,] 0.1451386 0.2376577 0.4567325 0.5555612 0.7043876

• Next, we treat the data as from a stratified randomized experiment and use Neyman’s estimator

## CI from Neyman
est_per_block <- sapply(1:max(lalonde$blocks), function(idx) {
  Yt <- lalonde$re78[lalonde$blocks == idx & lalonde$treat == 1]
  Yc <- lalonde$re78[lalonde$blocks == idx & lalonde$treat == 0]
  
  tau_hat <- mean(Yt) - mean(Yc)

  st <- sd(Yt)
  sc <- sd(Yc)
  varhat_tauhat <- st^2/length(Yt) + sc^2/length(Yc)
  return(c(tau_hat, varhat_tauhat))
})
colnames(est_per_block)<- 1:max(lalonde$blocks)
rownames(est_per_block) <- c("mean", "var")
est_per_block

               1          2            3             4           5
mean    2252.462   -1206.69     3752.993     -980.4199    1318.621
var  7019811.047 2558729.60 11917739.160 12153089.3809 1021087.344

wts <- table(lalonde$blocks)/length(lalonde$blocks)

tau_hat <- sum(est_per_block[1, ] * wts)
sd_tau_hat <- sqrt(sum(est_per_block[2, ] * wts^2))
  

result <- c(tau_hat, sd_tau_hat, c(tau_hat- 1.96 * sd_tau_hat, tau_hat + 1.96 * sd_tau_hat))
names(result) <- c("est", "sd", "CI_lower", "CI_upper")
result

      est        sd  CI_lower  CI_upper 
1201.8050  864.2807 -492.1852 2895.7952 

This result is very close to what we get from IPW.