# Modeling choices in matched, cluster-randomized experiments

Suppose an experiment with a matched, cluster-randomized design. Units are nested in groups, groups are matched through pairwise matching via Mahalanobis distance, and a treatment is randomly assigned to one group within the matched pairs.

1. It is often proposed to include only variables for the matched pairs in the model, but is it best practice to instead include only the variables used in the matching algorithm? Bruhn and McKenzie (2009) have advocated the latter, but perhaps one should include both (presumably this is not perfectly co-linear as the matching function is a non-linear function of the input variables).

2. I have seen many researchers incorporate these designs using statistical software intended for survey data analysis. For example, researchers will consider matched pairs (or blocks in block-randomized cases) as strata and declare the pairs as strata in the software, while declaring clustered SEs by declaring clusters. However, (a) what these packages are actually estimating when strata are declared and is it equivalent to including the indicator variables? (b) Are more efficient estimators than these weighted-least-squares estimators?

Here's a working example comparing a survey package in R with strata declared (i.e. matched pairs) and OLS with indicators. As we see, they produce different results.

# Load packages
library(dplyr)
library(survey)

# Create 3-strata dataset
set.seed(123)
strata <- factor(c(1,2,3,1,2,3))
y <- rnorm(6,0,1)
x <- seq(1:6)
dat <- data.frame(strata, y, x)

# Declare survey design with strata, no clusters
des <- survey::svydesign(ids = ~0, strata = strata, data = dat)

# Estimate survey linear model with strata
m_svy <- survey::svyglm(y ~ x, design = des) %>% summary

# Estimate model with dummies for each stratum
m_dummy <- lm(y ~ x + strata, data = dat) %>% summary

# Pull out estimates
svy <- m_svy$coefficients[ , 1:2] %>% round(2) %>% as.data.frame dum <- m_dummy$coefficients[1:2,1:2] %>% round(2) %>% as.data.frame

# Produce estimate table
mtab <- dplyr::bind_rows(svy, dum)
rownames(mtab) <- c("Intercept_svy","X_svy","Intercept_dum","X_dum")

# Check outcomes
mtab