# Implementing analytic block-cluster-robust standard errors

How do we calculate block-cluster-robust SEs for the average-treatment effect? (Note, I do not want block bootstrap. I want the analytic estimate, calculated with block-population-weighted block-level SE estimates.)

This is for a research design with blocks, clusters within blocks, and we want to use Eicker-Huber-White robust SEs.

To calculate the block SEs we need to calculate the SE within each block, and weight by the share the observations in each block.

Below you'll see a function that calculates cluster-robust SEs.

The first problem is how to integrate the blocking adjustment into the function, but I cannot figure it out. At present the function outputs a covariance matrix, and only calculate SEs later in the coeftest() function, which prevents us from calculating SEs by block.

A second related question, I find no resources discussing estimation of SEs that are blocked and clustered and robust. Why? Is there any reason why I am not finding resources? Is there any reason to avoid estimating block-cluster-robust SEs?

  remove(list = ls())

require(sandwich, quietly = TRUE)
require(lmtest, quietly = TRUE)
require(tidyverse)

set.seed(42)

N <- 560
k <- 56

data <- data.frame(id = 1:N)

# Simulate data with outcome, treatment, block, and cluster
data <-
data %>%
mutate(y1 = rnorm(n = N),
z = rep(x = c(1,0), each = 10, times = k/2),
block = rep(x = c(1,0), each = N/2),
cluster = rep(seq(1:k), each = 10))

#write your own function to return variance covariance matrix under clustered SEs
get_CL_vcov<-function(model, cluster){
#calculate degree of freedom adjustment
M <- length(unique(cluster))
N <- length(cluster)
K <- model$rank dfc <- (M/(M-1))*((N-1)/(N-K)) #calculate the uj's uj <- apply(estfun(model),2, function(x) tapply(x, cluster, sum)) #use sandwich to get the var-covar matrix vcovCL <- dfc*sandwich(model, meat=crossprod(uj)/N) return(vcovCL) } # Define a model m1<-lm(y1 ~ z, data=data) #call our new function and save the var-cov matrix output in an object m1.vcovCL <- get_CL_vcov(m1, data$cluster)

#the regular OLS standard errors
coeftest(m1)

#the clustered standard errors by indicating the correct var-covar matrix
coeftest(m1, m1.vcovCL)