R + FE 回归具有许多固定效应，我提取了 fe 系数和标准误差

Question

我有一个包含大约 200 万个人的数据集。数据集的简化版本包含三种类型的变量：（1）因变量（y），（2）协变量（x1，x2），（3）区域固定效应（reg，大约 200 万）。

（对于区域级别的第二次回归）我需要提取区域固定效应及其标准误。主要问题是我的电脑卡住了，显然是因为内存不足（我有 16GB，4 核 macbook pro），而目前我无法使用更强大的机器。

使用的策略：

1. 简单的lm(y ~ x1 + x2 + reg,data=DT) 卡住了

   • 我已经尝试减少区域假人的数量（例如，在不同地区投放人员），当我有大约 500 个假人时，这很有效。尽管如此，这意味着要删除我一半的数据集。

2. 使用 lm_robust：

• lm_robust(y ~ x1 + x2, fixed_effects= ~ reg, se_type='stata',data = DT) 有效，比较快，给出固定效应（fixed_effects）的系数，但不给出标准误差
• lm_robust(y ~ x1 + x2 + reg, se_type='stata',data = DT) 给了我固定效应系数和标准误差，但是当我使用超过 500 个区域时会卡住

1. 使用 felm 或 feols

   • felm(y ~ x1 + x2 | reg, data=DT)
   • feols(y ~ x1 + x2 | reg, data=DT)

   两种方法都有效，速度非常快，但输出中没有 FE 系数和标准误差

我考虑过的一个策略是使用 felm/feols 并使用输出“residuals”（整个系统的残差，有假人）和“r.residuals”（没有假人的预测产生的残差）退出区域残差（即对于区域 i，hat(reg)_i = (y-hat(b1)*x1 - hat(b2)*x2 - hat(reg_i)) - (y-hat(b1)*x1 - hat( b2)*x2). 问题是我不会得到标准错误. 一种解决方案是使用引导来计算标准错误. 这应该可以工作(我还没有尝试过), 但是好像太复杂了。

Answer 1

正如我所建议的，我最终通过自举计算了标准误差。 我的数据包含四种类型的变量：（1）因变量（y），（2）协变量（x1，x2），（3）区域固定效应（reg，约200万）和组变量（clus）。

我遇到的另一个复杂情况是我需要计算聚类标准误差（clus 是组变量）。这是我使用的代码。

# Construct variable with all the cluster groups. This will be used later to create the bootstrap samples
DT[, `:=` (reg = as.factor(reg),clus = as.factor(clus)) ]
clusters <- as.data.table(table(DT$clus))

# To speed up the process I parallelize my code
cores <- detectCores() -1
registerDoParallel(cores)
b <- 7000   # Number of bootstrap repetitions

set.seed(123456)
foreach (i=1:b, .combine = "rbind") %dopar% {
  # Create a sample with replacement, where the sample is constructed taken into account the clustered nature of the data
  units <- sample(clusters[[1]], size = nrow(clusters), replace=T)
  DT_b <- DT_res[as.character(units), on = "occind",allow.cartesian=TRUE]
  
  # Use felm to perform the regression. This way of doing it absorbs the fe, and thus can deal with a lot of fixed effects
  A <- felm(y ~ x1 + x2 | reg | 0 | clus , data=DT_b)

  # Extract the fixed effects using getfe
  A_fe <- as.data.table(getfe(A)[c("idx","effect")])
  setnames(A_fe,c(1,2),c("rn","coef"))
  A <- as.data.table(A[c("coefficients")],keep.rownames=TRUE)
  setnames(A,c(1,2),c("rn","coef"))

  A <- rbind(A,A_fe,fill=TRUE)
} -> coef.parallel

coef.parallel 数据集如下所示

rn       coef
x1       -0.123
x2       -0.321
reg1     0.0012
reg2     0.003
...
reg2000  -1.2
...
...
x1       0.124
x2       0.521
reg1     0.032
reg2     0.03
...
reg2000  -0.9

下一步就是简单地计算每个变量的平均值 (rn) 和标准误差

# Compute the bootstrap estimate and clustered standard errors
estimation.bootstrap.parallelized <- coef.parallel[,.(coef.b=mean(coef),se.b=sd(coef,na.rm=TRUE)),by="rn"]

estimation.bootstrap.parallelized 数据集如下所示

rn       coef.b   se.b
x1       -0.208   0.0032
x2       -0.321   0.123
reg1     0.0012   0.00041
reg2     0.003    0.00021
...
reg2000  -0.945     0.32