set.seed(4444)
library(estimatr)

Simulated population

N_clusters <- 100
cluster_size <- 10
N_pop <- N_clusters*cluster_size
cluster_index <- (1:100)%x%rep(1, cluster_size)
unit_index <- rep(1, N_clusters)%x%(1:cluster_size)
between_sd <- 2
cluster_random_effect <- rnorm(N_clusters, sd=between_sd)
cluster_random_effect_l <- cluster_random_effect%x%rep(1, cluster_size)

# This is a cluster-level covariate, 
# constant for units within the same cluster 
x <- (0.25*cluster_random_effect + rnorm(N_clusters, sd=1))%x%rep(1, cluster_size)

# Potential outcomes that exhibit high intra-cluster correlation
within_sd <- 1
y0 <- 0.25*cluster_random_effect_l + 0.25*x + rnorm(N_pop, sd=within_sd)
y1 <- y0 + .1*x + rnorm(N_pop, sd=within_sd)
dt_cluster <- data.frame( id = 1:N_pop,
                          cluster_index,
                          unit_index,
                          x,
                          y0,
                          y1)

Simulation setup:

nsim <- 10000
b_block_2 <- vhat_block_2 <- 
b_block_1 <- vhat_block_1 <-
b_block_nostrat <- vhat_block_nostrat <-
b_cluster <- vhat_cluster <- 
b_unit <- vhat_unit <- 
rep(NA, nsim)

Unit random assignment

I will use lm_robust from the estimatr package (loaded above) because it gives the right standard error.

for(i in 1:nsim){
    dt_cluster$z <- 0
    dt_cluster$z[dt_cluster$id %in% sample(dt_cluster$id, N_pop/2)] <- 1
    dt_cluster$y <- dt_cluster$z*dt_cluster$y1 + (1-dt_cluster$z)*dt_cluster$y0
    fit <- lm_robust(y~z, data=dt_cluster)
    b_unit[i] <- coef(fit)["z"]
    vhat_unit[i] <- diag(vcov(fit))["z"]
}

## Warning: Assigning non-quosure objects to quosure lists is deprecated as of rlang 0.3.0.
## Please coerce to a bare list beforehand with `as.list()`
## This warning is displayed once per session.

Unbiasedness and variance:

mean(y1-y0)

## [1] 0.02531441

mean(b_unit)

## [1] 0.0255967

var(b_unit)

## [1] 0.006822015

mean(vhat_unit)

## [1] 0.007893876

Recall that the standard error estimator here is the estimator of the Neyman variance. This will be conservative, as the results above indicate.

Cluster random assignment

Again, I use lm_robust here because it gives the right standard error.

for(i in 1:nsim){
    dt_cluster$z <- 0
    dt_cluster$z[dt_cluster$cluster_index %in% sample(unique(dt_cluster$cluster_index), N_clusters/2)] <- 1
    dt_cluster$y <- dt_cluster$z*dt_cluster$y1 + (1-dt_cluster$z)*dt_cluster$y0
    fit <- lm_robust(y~z,
                                    data=dt_cluster,
                                    clusters=cluster_index)
    b_cluster[i] <- coef(fit)["z"]
    vhat_cluster[i] <- diag(vcov(fit))["z"]
}

Unbiasedness and variance:

mean(y1-y0)

## [1] 0.02531441

mean(b_cluster)

## [1] 0.02586395

var(b_cluster)

## [1] 0.02474477

mean(vhat_cluster)

## [1] 0.02628188

Design effect relative to unit random assignment:

var(b_cluster)/var(b_unit)

## [1] 3.627194

Block cluster assignment

Note that we need to implement the blocking in the random assignment, and we should also account for the blocking in the estimation. For this particular simulation, the point estimates will be the same even if we do not account for the stratification because we have a balanced design (half treated, half control) with equal sized blocks. But if there is heterogeneity in block sizes or share of units treated per block, then it is important to account for the stratification. Also, your standard error estimates will be too big if you do not account for the stratification. You will see all of this below.

To account for the blocking in the analysis, you use block fixed-effects. You could just insert fixed effects for each block. This is okay if the treatment assignment probabilities are uniform across all blocks. That condition holds in our simulation, so it is fine. But it may not always be the case — for example if you have different numbers of units or clusters per block. In those cases, there are two potential estimators that are design-consistent. The first includes the block fixed effects and also weights by the inverse of the probabilty of assignment. See either my notes or Gerber and Green (2012, Ch. Ch. 4). Another approach (that is algebraically equivalent to the weighted fixed effects estimator) is the “centered interaction”" specification, which includes the block fixed effects as well as the interactions between the centered block fixed effects and the treatment variable. This can be implemented using the lm_lin function in the estimatr package.

The code below covers all of these estimators, even though in our simulation they don’t make a difference for the point estimates. They do make a difference for the standard error/variance estimates though, as you will see.

N_blocks <- 10
clusters_per_block <- N_clusters/N_blocks
dt_cluster <- dt_cluster[order(dt_cluster$x),]
dt_cluster$block <- (1:N_blocks%x%rep(1, clusters_per_block))%x%rep(1, cluster_size)


for(i in 1:nsim){
    dt_cluster$z <- 0
    for(j in unique(dt_cluster$block)){
        clus_treat <- sample(unique(dt_cluster[dt_cluster$block==j,"cluster_index"]), 
                            clusters_per_block/2)
        dt_cluster[dt_cluster$block==j&dt_cluster$cluster_index%in% clus_treat,"z"] <- 1
    }
    dt_cluster$y <- dt_cluster$z*dt_cluster$y1 + (1-dt_cluster$z)*dt_cluster$y0
    # This specification is ignoring the stratification like most of you did in your simulations
    fit <- lm_robust(y~z,
                     data=dt_cluster,
                     clusters = cluster_index)
    b_block_nostrat[i] <- coef(fit)["z"]
  vhat_block_nostrat[i] <- diag(vcov(fit))["z"]
  # Here is how you would construct the inverse probability weights for a 
    # block-cluster randomized experiment:
    z_aggregated <- aggregate(dt_cluster[,c("x","z","block")],
                            by=list(cluster=dt_cluster$cluster_index),
                            mean)
  z_aggregated <- z_aggregated[order(z_aggregated$x),]
  p_treatment <- tapply(z_aggregated$z,
                        z_aggregated$block,
                        mean)
  dt_cluster$ipw <- NA
  for(k in unique(dt_cluster$block)){
    dt_cluster[dt_cluster$block==k,"ipw"]  <- dt_cluster[dt_cluster$block==k,
                                                         "z"]*(1/p_treatment[as.numeric(names(p_treatment))==k]) +
                                              (1-dt_cluster[dt_cluster$block==k,
                                                            "z"])*(1/(1-p_treatment[as.numeric(names(p_treatment))==k]))
  }

    # Weighted Block FE (note that the weights in our simulation are uniform, so not actually needed,
  # but I am including them for completeness)
  fit <- lm_robust(y~z+as.factor(block),
                         data=dt_cluster,
                                         clusters=cluster_index,
                                         weights=ipw)
    b_block_1[i] <- coef(fit)["z"]
    vhat_block_1[i] <- diag(vcov(fit))["z"]
    # Block FE with centered interaction specification
    fit <- lm_lin(y~z,
                  covariates=~as.factor(block),
                  clusters=cluster_index,
                  data=dt_cluster)
    b_block_2[i] <- coef(fit)["z"]
    vhat_block_2[i] <- diag(vcov(fit))["z"]
}

Unbiasedness and variance

mean(y1-y0)

## [1] 0.02531441

mean(b_block_nostrat)

## [1] 0.02310357

mean(b_block_1)

## [1] 0.02310357

mean(b_block_2)

## [1] 0.02310357

var(b_block_nostrat)

## [1] 0.01298488

var(b_block_1)

## [1] 0.01298488

var(b_block_2)

## [1] 0.01298488

mean(vhat_block_nostrat)

## [1] 0.02641759

mean(vhat_block_1)

## [1] 0.01390732

mean(vhat_block_2)

## [1] 0.0138043

You can see that the estimator without the fixed effect does not return the right variance estimate (and so the estimated standard error would be way too large too).

Design effect relative to cluster random assignment and unit random assignment

var(b_block_2)/var(b_cluster)

## [1] 0.5247525

var(b_block_2)/var(b_unit)

## [1] 1.903379

Design Effects Exercise

Cyrus Samii

9/23/2020

Simulated population

Unit random assignment

Cluster random assignment

Block cluster assignment