```
set.seed(4444)
library(estimatr)
```

```
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)
```

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.

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`

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`