基于Gauss-Seidel算法的R高维回归固定效应估计报错排查

问题背景 #

拥有包含20年信息的匹配雇主-雇员数据集（约200万个体、40万家企业），目标是通过工资方程，用Gauss-Seidel算法实现精确最小二乘估计，得到个体、企业、职业三类固定效应向量。因Stata收敛过慢转用R编写函数，修复缺失值导致的模型矩阵与标识符不匹配问题后，出现新错误：

Error in h(simpleError(msg, call)) :  error in evaluating the argument 'x' in selecting a method for function 'mean': non-conformable arguments

核心错误分析与修复步骤 #

1. 修正beta更新的逻辑错误 #

原代码在个体/企业/职业循环中错误地用lm(y_i - X_i %*% beta ~ 0)$coefficients更新beta，这行代码是对无自变量的常数回归，仅返回因变量均值，导致beta维度从协变量列数变为1，后续X_i %*% beta必然出现矩阵维度不兼容。

Gauss-Seidel的正确逻辑是：更新某一类固定效应时，固定当前beta和另外两类固定效应，计算当前组的固定效应；仅在所有固定效应更新完成后，统一用残差回归协变量更新beta。

2. 修复固定效应的索引匹配问题 #

原代码用df$ntrab == i匹配个体，但i是固定效应数组的序号（1、2...），而非ntrab的实际标识符，若标识符不是连续整数，会导致无法匹配到对应观测，甚至生成空向量引发计算错误。

正确做法是将标识符转为因子，利用因子水平索引匹配：

# 转换为因子，确保索引与固定效应数组位置对应
df$ntrab = factor(df$ntrab)
df$emp_id = factor(df$emp_id)
df$occupation4 = factor(df$occupation4)

# 按因子水平数初始化固定效应
worker_fe = rep(0, nlevels(df$ntrab))
firm_fe = rep(0, nlevels(df$emp_id))
occupation_fe = rep(0, nlevels(df$occupation4))

3. 确保矩阵维度一致性 #

• 生成协变量矩阵时显式去除截距，避免多余列干扰：

  X = model.matrix(~ p_age2 + antig + antig2 + firm_age + lfprod - 1, df)
  

• 删除循环内的beta更新操作，仅保留所有固定效应更新完成后的beta统一更新步骤。

4. 彻底处理缺失值 #

提前删除所有含缺失值的观测，确保后续生成的模型矩阵、标识符向量行数完全一致：

df = na.omit(df)

修正后的代码框架示例 #

estimate_fixed_effects = function(df, max_iter = 1000, tol = 1e-6) {
  # 预处理：删除缺失值，转换标识符为因子
  df = na.omit(df)
  df$ntrab = factor(df$ntrab)
  df$emp_id = factor(df$emp_id)
  df$occupation4 = factor(df$occupation4)
  
  y = df$real_lrganho
  # 协变量矩阵：无截距
  X = model.matrix(~ p_age2 + antig + antig2 + firm_age + lfprod - 1, df)
  beta = rep(0, ncol(X))
  
  # 初始化固定效应（对应因子水平顺序）
  worker_fe = rep(0, nlevels(df$ntrab))
  firm_fe = rep(0, nlevels(df$emp_id))
  occupation_fe = rep(0, nlevels(df$occupation4))
  
  for (iter in 1:max_iter) {
    beta_old = beta
    
    # 更新个体固定效应：控制beta、企业和职业固定效应
    for (i in 1:nlevels(df$ntrab)) {
      idx = which(df$ntrab == levels(df$ntrab)[i])
      y_i = y[idx]
      X_i = X[idx, ]
      other_effects = firm_fe[df$emp_id[idx]] + occupation_fe[df$occupation4[idx]]
      resid_i = y_i - X_i %*% beta - other_effects
      worker_fe[i] = mean(resid_i)
    }
    
    # 更新企业固定效应：控制beta、个体和职业固定效应
    for (j in 1:nlevels(df$emp_id)) {
      idx = which(df$emp_id == levels(df$emp_id)[j])
      y_j = y[idx]
      X_j = X[idx, ]
      other_effects = worker_fe[df$ntrab[idx]] + occupation_fe[df$occupation4[idx]]
      resid_j = y_j - X_j %*% beta - other_effects
      firm_fe[j] = mean(resid_j)
    }
    
    # 更新职业固定效应：控制beta、个体和企业固定效应
    for (k in 1:nlevels(df$occupation4)) {
      idx = which(df$occupation4 == levels(df$occupation4)[k])
      y_k = y[idx]
      X_k = X[idx, ]
      other_effects = worker_fe[df$ntrab[idx]] + firm_fe[df$emp_id[idx]]
      resid_k = y_k - X_k %*% beta - other_effects
      occupation_fe[k] = mean(resid_k)
    }
    
    # 更新beta：控制所有固定效应后回归协变量
    total_fe = worker_fe[df$ntrab] + firm_fe[df$emp_id] + occupation_fe[df$occupation4]
    beta = lm(y - total_fe ~ X - 1)$coefficients
    
    # 收敛判断
    if (max(abs(beta - beta_old)) < tol) {
      cat("收敛于迭代次数：", iter, "\n")
      break
    }
  }
  
  return(list(
    beta = beta,
    worker_fe = data.frame(ntrab = levels(df$ntrab), fe = worker_fe),
    firm_fe = data.frame(emp_id = levels(df$emp_id), fe = firm_fe),
    occupation_fe = data.frame(occupation4 = levels(df$occupation4), fe = occupation_fe)
  ))
}

额外优化建议 #

• 对于百万级数据集，三重循环效率极低，可改用data.table的分组向量化操作替代循环，大幅提升速度：

  library(data.table)
  setDT(df)
  # 示例：向量化更新个体固定效应
  df[, total_other := firm_fe[emp_id] + occupation_fe[occupation4]]
  df[, resid := y - X %*% beta - total_other]
  worker_fe = df[, mean(resid), by = ntrab]$V1
  

• 可直接使用专门的固定效应估计包（如lfe），其内部实现了更高效的迭代算法，稳定性和速度远优于手动编写的Gauss-Seidel。

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