如何在R的lm固定效应回归对象中彻底移除冗余因子变量？

核心思路是先将因变量和感兴趣的核心自变量对所有冗余控制变量（固定效应、地区特定趋势）做残差化，再用残差化后的变量拟合回归，这样模型仅包含你关心的10个系数，大幅提升coeftest()结合vcovHAC的运行速度，同时保留原模型的核心系数估计和残差特性。

具体步骤与代码示例 #

假设你的原模型公式是：

formula_full <- Y ~ X1 + X2 + ... + X10 + factor(FE) + factor(distr):trend
model_full <- lm(formula_full, data = your_data)

1. 分离核心变量与控制变量
   提取你关心的核心自变量（X1到X10），以及所有控制变量（固定效应、地区趋势）：

   # 核心自变量矩阵（去除截距项）
   X_main <- model.matrix(~ X1 + X2 + ... + X10 - 1, data = your_data)
   # 控制变量矩阵（包含固定效应和地区特定趋势，去除截距项）
   X_controls <- model.matrix(~ factor(FE) + factor(distr):trend - 1, data = your_data)
   # 因变量
   Y <- your_data$Y
   

2. 对变量做残差化（去除控制变量的影响）
   计算投影矩阵，将Y和核心变量投影到控制变量的正交补空间：

   # 用qr.solve提升大矩阵下的计算效率
   M <- diag(nrow(your_data)) - X_controls %*% qr.solve(t(X_controls) %*% X_controls) %*% t(X_controls)
   # 残差化后的因变量和核心自变量
   Y_resid <- as.vector(M %*% Y)
   X_main_resid <- as.data.frame(M %*% X_main)
   

3. 拟合精简模型
   用残差化后的变量拟合回归，此时模型仅包含10个核心系数：

   model_slim <- lm(Y_resid ~ ., data = X_main_resid)
   

4. 计算HAC标准误
   现在用coeftest()结合vcovHAC处理精简模型，速度会大幅提升：

   library(lmtest)
   library(sandwich)
   coeftest(model_slim, vcov = vcovHAC)
   

关键说明 #

• 这个方法得到的核心系数估计值与原lm()模型完全一致，因为固定效应模型中，控制变量的投影操作不会改变核心变量的系数结果。
• 精简模型的残差与原模型的残差完全相同，你可以直接将其与原数据合并做诊断：

  your_data$residuals <- residuals(model_slim)
  

替代方案：修改原lm模型对象（需谨慎） #

如果你坚持要修改原模型对象，需要同步修改多个关键元素（仅推荐熟悉lm对象结构的用户尝试）：

# 定义你关心的系数名称
keep_coefs <- c("X1", "X2", ..., "X10")
# 获取系数索引
keep_idx <- which(names(model_full$coefficients) %in% keep_coefs)

# 修改模型对象的核心元素
model_slim <- model_full
model_slim$coefficients <- model_slim$coefficients[keep_idx]
model_slim$rank <- length(keep_idx)
model_slim$qr$qr <- model_slim$qr$qr[1:model_slim$rank, 1:model_slim$rank]
model_slim$terms <- delete.response(model_slim$terms)
model_slim$terms <- terms(paste("~", paste(keep_coefs, collapse = "+")))
attr(model_slim$terms, "intercept") <- as.integer(any(names(model_slim$coefficients) == "(Intercept)"))

# 此时coeftest可以正常处理
coeftest(model_slim, vcov = vcovHAC)

注意：这种方法容易遗漏模型对象的其他关联元素，可能导致后续分析（如预测）出错，因此优先推荐残差化方法。

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