优化在 R lm() 模型中运行固定效果的时间答案

【问题标题】：Optimize the time to run fixed effects in an R lm() model优化在 R lm() 模型中运行固定效果的时间
【发布时间】：2020-08-02 06:28:42
【问题描述】：

我正在尝试运行一个回归模型，其中包含美国城市的固定效应。我有超过 10,000,000 百万行和 600 个城市。下面的代码有效，但它真的很慢。当为具有很多级别的变量包含一个因子时，有什么方法可以更快地运行模型。

x <- data.frame(
    a = sample( 1:1000, 1000000 , replace=T),
    cityfips = sample( 1:250, 1000000 , replace=T),
    d = sample( 1:4, 1000000 , replace=T)
)

system.time(a1 <- lm( a~cityfips+d  , x ) )
system.time(a2 <- lm( a~as.factor(cityfips) + d  , x ) )




> system.time(a1 <- lm( a~cityfips+d  , x ) )
   user  system elapsed 
   0.22    0.00    0.22 
> system.time(a2 <- lm( a~as.factor(cityfips) + d  , x ) )
   user  system elapsed 
  95.65    0.97   96.62 
> system.time(a3 <- slm( a~as.factor(cityfips) + d  , x ) )
   user  system elapsed 
   4.58    2.06    6.65

【问题讨论】：

嗯，我实际上看到了您可能遇到的问题.. 是 model.matrix 正在杀死您的 lm，而重复的问题实际上并没有解决这个问题
您可以尝试使用稀疏矩阵，web.mit.edu/r/current/lib/R/library/SparseM/html/slm.html，但似乎它非常具有实验性
@StupidWolf - 做到了，我更新了给你归属的问题，希望我能投票给你。谢谢
不用担心...好吧，我检查了一下，好像系数是一样的。是的，您需要小心，配件的完成方式略有不同。你可以阅读更多econ.uiuc.edu/~roger/research/sparse/SparseM.pdf
一般来说，您不应该编辑您的问题以包含答案。如果现有答案未涵盖答案的某些方面，您可以发布自己的答案。

标签： r lm

【解决方案1】：

当你有这么多因素时，在lm() 中构造model.matrix 将占用大部分时间，一种方法是使用glmnet 中的sparseMatrix，并且有两个包，sparseM，@987654324 @ 允许 lm 进入 sparseMatrix：

set.seed(111)
x <- data.frame(
    a = sample( 1:1000, 1000000 , replace=T),
    cityfips = sample( 1:250, 1000000 , replace=T),
    d = sample( 1:4, 1000000 , replace=T)
)

library(SparseM)
library(MatrixModels)
library(Matrix)

system.time(f_lm <- lm( a~as.factor(cityfips) + d  , x )  )
   user  system elapsed 
 75.720   2.494  79.365  
system.time(f_sparseM <- slm(a~as.factor(cityfips) + d  , x ))
   user  system elapsed 
  5.373   3.952  10.646
system.time(f_modelMatrix <- glm4(a~as.factor(cityfips) + d  ,data=x,sparse=TRUE))
   user  system elapsed 
  1.878   0.335   2.219

我能找到的最接近的是MatrixModels 中的 glm4，但您可以在下面看到系数与使用 lm 的拟合相同：

all.equal(as.numeric(f_sparseM$coefficients),as.numeric(f_lm$coefficients))
[1] TRUE
all.equal(as.numeric(f_lm$coefficients),as.numeric(coefficients(f_modelMatrix)))
[1] TRUE

除了MatrixModels 中的glm4 之外，另一个选择是使用 lm.fit（正如@BenBolker 所指出的那样：

lm.fit(x=Matrix::sparse.model.matrix(~as.factor(cityfips) + d,data=x),y=x$a)

这会给你一个像lm.fit()这样的列表，你不能应用summary()等函数。

这两个软件包的作者都警告说它是实验性的，因此与 stats::lm 相比可能仍然存在一些差异，请注意检查。

【讨论】：

你也可以使用Matrix::sparse.model.matrix并将结果传递给lm.fit()
感谢@BenBolker 指出。我正在看它，但不确定 OP 需要 lm 对象做什么。例如，如果您需要执行 summary(..) 它不适用于 lm.fit() 或 MatrixModels:::lm.fit.sparse()
但是是的，你是对的，也许我把它包括在答案中

【解决方案2】：

查看lfe package。我还没有深入研究算法的细节，但至少根据我的经验，它在很短的时间内产生了与lm() 完全相同的结果。

作为奖励，它可以很容易地对标准错误进行聚类，因此您无需在之后进行任何聚类和/或三明治估计器业务，尽管这样做的语法有点不寻常

【讨论】：