运行超过 10 K 的回归并存储系数

如何解决运行超过 10 K 的回归并存储系数

八年来，我拥有 73000 个社区的人口。我需要为每个提取线性趋势的斜率。所以我需要估计 73000 次 2021-08-08 00:00:00 模型。并将它们的斜率 (b) 存储在列表中。我有专门针对此类问题的 R 包吗？

解决方法

假设您将每个邻域的邻域数据作为列表的单独元素，则可以使用 sapply()。这是假设您想在每个社区运行完全不相关的回归。这是一个两个街区的例子：

nbr_list <- list(n1 = c(1,4,3,5,6),n2 = c(4,6,7,9,9))

get_slope <- function(y){
  slopeOut <- lm(y~I(1:8))$coefficients[2]
  slopeOut }

sapply(nbr_list,FUN=get_slope)
#> n1.I(1:8) n2.I(1:8) 
#> 0.5238095 0.7142857

^{由 reprex package (v0.3.0) 于 2021 年 2 月 23 日创建}

tidyverse 的解决方案是使用包 df 和 library(dplyr) df %>% group_by(Chain,Product) %>% mutate(I = as.integer(row_number() == match(TRUE,Sale > 0))) -> result result 按组运行许多回归并将它们存储到列表列中，然后提取 beta。主要优点是您可以拥有更复杂的模型，并且整个回归输出都存储在同一个数据帧中，您只需要在之后检索它们。最好保留所有输出以备将来使用（例如，如果您只存储斜率并且接下来需要截距，则需要再次重新运行所有回归）。

嵌套解决方案如下：

按社区对您的数据进行分组。
将这些组嵌套到列表列中。
创建一个模型列，将 lm 函数映射到数据列表列。
取消嵌套模型列表列。
将其过滤为想要的预测系数。

broom

或者更简单，没有嵌套：

purrr

我本来打算告诉您使用 lme4::lmList()，但是如果您设置正确，普通的旧式 lm 会快得多。

样本数据

nn <- 73000
ny <- 8
set.seed(101)
dd <- data.frame(
    nbr=rep(1:nn,each=ny),yr=rep(1:ny,nn),y=rnorm(nn*ny)
)

使用 lm()

Y <- matrix(dd$y,nrow=8)  ## set up matrix of response variables
yr <- dd$yr[1:8]          ## predictor variables
system.time(mm2 <- lm(Y~yr))
cc2 <- coef(mm2)
head(cc2["yr",])
## [1]  0.088903516 -0.138058182  0.006353931  0.106540632 -0.021719822
## [6] -0.020777392

这或多或少是瞬时的（在快速机器上

使用 lme4::lmList()

这是我的第一个想法，可能更容易理解（并且会更普遍地工作，例如，如果您对每个组有不同的预测变量），但速度要慢 100-1000 倍。

library(lme4)
system.time(mm <- lmList(y~yr|nbr,data=dd)) ## fit models,34 seconds
system.time(cc <- coef(mm))  ## extract coefs,43 seconds
head(cc$yr)
## [1]  0.088903516 -0.138058182  0.006353931  0.106540632 -0.021719822
## [6] -0.020777392