OLS与自助聚类标准误线性模型的ANOVA分析实现求助

问题背景 #

需要对比标准OLS回归模型和通过fwildclusterboot包boottest函数得到的聚类自助稳健标准误模型的方差分析结果，但直接对自助得到的系数矩阵调用anova()函数无法运行，用户尝试的错误代码如下：

# "lm" is my linear model
lm_coefnames <- c("treatment", "var1", "var2", "Intercept")
set.seed(2020)

boot_lm <- matrix(NA, length(lm_coefnames), 4)
for (i in 1:length(lm_coefnames)){
  
  boot_lm[i, ] <- as.numeric(summary(boottest(lm, clustid = "cluster_variable", 
                           param = lm1_coefnames[i], B = 9999))[1, c(2,4:6)])
}

anova(summary(boot_lm))

核心问题：anova()函数仅支持接收统计模型对象（如lm的输出结果），无法直接处理自助抽样得到的系数矩阵，因此直接传入矩阵会触发报错。

可行解决思路 #

思路1：用boottest直接做联合F检验（对应ANOVA整体显著性） #

ANOVA的核心是检验模型整体或部分系数的联合显著性，boottest支持直接对多个参数执行联合检验，无需逐个提取系数后再处理：

set.seed(2020)
# 检验除截距外的所有自变量联合显著性（对应ANOVA的模型整体检验）
boot_anova <- boottest(
  lm, 
  clustid = "cluster_variable", 
  param = c("treatment", "var1", "var2"),  # 指定需联合检验的参数
  B = 9999,
  test = "F"  # 选择F检验，匹配ANOVA的联合显著性逻辑
)

# 查看聚类自助稳健的ANOVA结果
summary(boot_anova)

该结果直接对应标准OLS模型anova(lm)中的整体F检验，但使用的是聚类自助稳健的标准误与p值。

思路2：构造自助样本模型集合，批量生成ANOVA结果 #

如果需要对比每个系数的ANOVA分解结果（类似anova(lm)的逐行输出），可以通过循环生成每个自助样本的回归模型，再批量提取ANOVA统计量：

set.seed(2020)
# 获取聚类ID的唯一值
clusters <- unique(lm$model$cluster_variable)
B <- 9999

# 存储每次自助抽样的ANOVA结果
anova_results <- list()

for (b in 1:B) {
  # 对聚类进行放回式自助抽样
  boot_clusters <- sample(clusters, replace = TRUE)
  # 构造自助数据集：按抽样的聚类ID合并原始数据
  boot_data <- do.call(rbind, lapply(boot_clusters, function(c) lm$model[lm$model$cluster_variable == c, ]))
  # 拟合OLS模型
  boot_model <- lm(formula(lm), data = boot_data)
  # 存储当前自助样本的ANOVA结果
  anova_results[[b]] <- anova(boot_model)
}

# 整理目标变量的ANOVA统计量分布，以treatment为例
treatment_f_values <- sapply(anova_results, function(x) x["treatment", "F value"])
treatment_p_values <- sapply(anova_results, function(x) x["treatment", "Pr(>F)"])

# 计算统计量的均值、分位数，与原OLS的ANOVA结果对比
mean(treatment_f_values)
quantile(treatment_f_values, c(0.025, 0.975))

这种方法模拟了聚类自助下的ANOVA统计量分布，可直接与标准OLS的ANOVA结果做分布层面的对比。

思路3：手动计算稳健ANOVA统计量 #

ANOVA的F统计量可通过系数估计值和稳健方差-协方差矩阵手动计算：

1. 用sandwich包获取聚类稳健的方差-协方差矩阵：

library(sandwich)
# 生成聚类稳健的VCV矩阵
vcv_cluster <- vcovCL(lm, cluster = lm$model$cluster_variable)

2. 构造假设矩阵，手动计算F统计量（以检验treatment、var1、var2联合为0为例）：

# 构造假设矩阵：排除截距行，检验前3个系数联合为0
hyp_mat <- diag(length(coef(lm)))[-4, ]
beta_hat <- coef(lm)

# 计算F统计量
f_stat <- t(hyp_mat %*% beta_hat) %*% solve(hyp_mat %*% vcv_cluster %*% t(hyp_mat)) %*% (hyp_mat %*% beta_hat) / nrow(hyp_mat)

这种方法灵活性强，适合自定义的ANOVA检验场景，p值可通过自助法或渐近分布进一步计算。

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