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