嵌套层级数据下，如何修正回归模型的标准误？

你的担心是对的：普通OLS（lm拟合）假设所有观测的残差独立，但同一县内的教区共享未观测的县层面特征，残差会存在聚类相关性；再加上X3、X4是县层面变量，每个县内的取值完全相同，它们的实际变异仅来自40个县之间，而非3000个教区。这两种因素会导致OLS严重低估这两个变量的标准误，以下是三种可行的修正方法：

1. 聚类稳健标准误（最直接的OLS修正） #

不需要改变原模型的系数估计，仅调整标准误以考虑县内的聚类相关性，适合你想保留OLS系数解释、同时修正标准误的场景。

使用sandwich和lmtest包实现：

# 安装依赖包（首次使用时运行）
install.packages(c("sandwich", "lmtest"))
library(sandwich)
library(lmtest)

# 拟合原OLS模型
fit <- lm(Y ~ X1 + X2 + X3 + X4, data = mydata)

# 按县（假设你的县变量名为county）计算聚类稳健标准误
cluster_results <- coeftest(fit, vcov = vcovCL, cluster = ~ county)
print(cluster_results)

这个方法会基于县层面的自由度（40个县）计算X3、X4的标准误，避免因错误使用教区层面自由度导致的低估。

2. 多层线性模型（嵌套数据的专属建模） #

直接建模嵌套结构，将县层面的未观测变异作为随机效应，能同时区分教区和县层面的变异，更贴合数据的实际结构。

使用lme4和lmerTest包实现（lmerTest用于提供p值）：

# 安装依赖包（首次使用时运行）
install.packages(c("lme4", "lmerTest"))
library(lme4)
library(lmerTest)

# 拟合带县层面随机截距的分层模型
hier_fit <- lmer(Y ~ X1 + X2 + X3 + X4 + (1 | county), data = mydata)
summary(hier_fit)

模型中的(1 | county)表示每个县有独立的截距，捕捉县层面的未观测异质性。县层面变量X3、X4的系数标准误会自动基于县之间的变异计算，结果更可靠。如果需要检验县层面的斜率变异，还可以扩展为随机斜率模型（比如(X1 + X2 | county)）。

3. 聚合数据回归（简化的县层面分析） #

将教区层面的数据聚合到县层面，直接用40个县的观测拟合回归，标准误自然反映县层面的变异。但这个方法会丢失教区层面的信息，估计效率较低，适合仅关注县层面效应的场景。

使用dplyr包聚合数据：

# 安装依赖包（首次使用时运行）
install.packages("dplyr")
library(dplyr)

# 按县聚合变量（X3、X4每个县内取值相同，直接取第一个值即可）
county_level_data <- mydata %>%
  group_by(county) %>%
  summarize(
    Y_avg = mean(Y),
    X1_avg = mean(X1),
    X2_avg = mean(X2),
    X3 = first(X3),
    X4 = first(X4)
  )

# 拟合县层面回归
agg_fit <- lm(Y_avg ~ X1_avg + X2_avg + X3 + X4, data = county_level_data)
summary(agg_fit)

内容的提问来源于stack exchange，提问作者Mikhail