按组进行快速线性回归答案

【问题标题】：Fast linear regression by group按组进行快速线性回归
【发布时间】：2015-06-30 11:46:12
【问题描述】：

我有 500K 个用户，我需要为每个用户计算线性回归 （带截距）。

每个用户大约有 30 条记录。

我尝试使用dplyr 和lm，这太慢了。用户大约 2 秒。

  df%>%                       
      group_by(user_id, add =  FALSE) %>%
      do(lm = lm(Y ~ x, data = .)) %>%
      mutate(lm_b0 = summary(lm)$coeff[1],
             lm_b1 = summary(lm)$coeff[2]) %>%
      select(user_id, lm_b0, lm_b1) %>%
      ungroup()
    )

我尝试使用lm.fit，它已知更快，但它似乎与dplyr 不兼容。

有没有快速按组进行线性回归的方法？

【问题讨论】：

答案可能取决于效率低下的来源。你做了一些分析吗？拆分数据是缓慢的过程，还是拟合模型是缓慢的部分？如果是后者，您可能需要查看例如fastLm 来自 RcppArmadillo 包。
慢的部分是：do(lm = lm(Y ~ x, data = .)) & mutate(lm_b0 = summary(lm)$coeff[1], lm_b1 = summary(lm)$coeff[2])。如何有效地在 50 万用户上使用 fastLm？
@psql 我必须做这个确切的事情，我能想到的最有效的方法就是这样做，它最终花费了大约 2 个小时 dplyr 这是来自 a 的一些基准测试问题stackoverflow.com/questions/29641366/…
在我的情况下大约需要 10 天，所以我真的需要加快这个解决方案！ :) ps：如果结果粗糙，那完全没问题！。
为什么不使用user_id 作为因子和交互项来拟合单个回归？ Y ~ x * user_id + 0 将为每个用户提供不同的斜率和截距（假设 user_id 是一个因素）。使用lm 会快很多，使用fastLm 会更快。

标签： r dplyr lm

【解决方案1】：

您可以只使用基本公式来计算斜率和回归。如果您只关心这两个数字，lm 会做很多不必要的事情。这里我使用data.table 进行聚合，但您也可以在base R 中进行聚合（或dplyr）：

system.time(
  res <- DT[, 
    {
      ux <- mean(x)
      uy <- mean(y)
      slope <- sum((x - ux) * (y - uy)) / sum((x - ux) ^ 2)
      list(slope=slope, intercept=uy - slope * ux)
    }, by=user.id
  ]
)

为 50 万用户生成约 30 个 obs（以秒为单位）：

 user  system elapsed 
 7.35    0.00    7.36

或大约 每位用户 15 微秒。

更新：我最终写了一堆blog posts 也涉及到这一点。

并确认这是否按预期工作：

> summary(DT[user.id==89663, lm(y ~ x)])$coefficients
             Estimate Std. Error   t value  Pr(>|t|)
(Intercept) 0.1965844  0.2927617 0.6714826 0.5065868
x           0.2021210  0.5429594 0.3722580 0.7120808
> res[user.id == 89663]
   user.id    slope intercept
1:   89663 0.202121 0.1965844

数据：

set.seed(1)
users <- 5e5
records <- 30
x <- runif(users * records)
DT <- data.table(
  x=x, y=x + runif(users * records) * 4 - 2, 
  user.id=sample(users, users * records, replace=T)
)

【讨论】：

我开始回答“如果你想要的只是系数......”，但如果这就是你真正想要的，那么这个答案当然是方法。
@BrodieG 为什么在fast_lm 函数中返回一个列表？不能使用矢量？ c(slope, intercept) ?
@Ricol，以便data.table 将它们解释为不同的列。

【解决方案2】：

如果你想要的只是系数，我会使用user_id 作为回归中的一个因素。使用@miles2know 的模拟数据代码（尽管重命名，因为exp() 以外的对象共享该名称对我来说看起来很奇怪）

dat <- data.frame(id = rep(c("a","b","c"), each = 20),
                  x = rnorm(60,5,1.5),
                  y = rnorm(60,2,.2))

mod = lm(y ~ x:id + id + 0, data = dat)

我们不拟合全局截距 (+ 0)，因此每个 id 的截距是 id 系数，并且没有 x 本身，因此 x:id 交互是每个 id 的斜率：

coef(mod)
#      ida      idb      idc    x:ida    x:idb    x:idc 
# 1.779686 1.893582 1.946069 0.039625 0.033318 0.000353

因此，对于id 的a 水平，ida 系数 1.78 是截距，x:ida 系数 0.0396 是斜率。

我将把这些系数收集到数据框的适当列中交给你...

此解决方案应该非常快，因为您不必处理数据帧的子集。使用fastLm 等可能会加快速度。

关于可扩展性的说明：

我只是在@nrussell 的模拟全尺寸数据上尝试过这个，但遇到了内存分配问题。根据您拥有的内存量，它可能无法一次性使用，但您可能可以分批使用用户 ID。他的回答和我的回答结合起来可能是最快的——或者 nrussell 的可能更快——将用户 ID 因子扩展到数千个虚拟变量可能计算效率不高，因为我一直在等待超过一个现在只需几分钟即可运行 5000 个用户 ID。

【讨论】：

sparse.model.matrix()（来自Matrix 包）和lm.fit 可能值得考虑。我很好奇lme4::lmer 会如何解决这个问题...
或者也许将用户分成一组并使用parallel或其他多核工具？

【解决方案3】：

更新： 正如 Dirk 所指出的，通过直接指定 x 和 Y 而不是使用基于公式的接口 fastLm，可以大大改进我的原始方法，这会产生（相当大的）处理开销。为了比较，使用原始的全尺寸数据集，

R> system.time({
  dt[,c("lm_b0", "lm_b1") := as.list(
    unname(fastLm(x, Y)$coefficients))
    ,by = "user_id"]
})
#  user  system elapsed 
#55.364   0.014  55.401 
##
R> system.time({
  dt[,c("lm_b0","lm_b1") := as.list(
    unname(fastLm(Y ~ x, data=.SD)$coefficients))
    ,by = "user_id"]
})
#   user  system elapsed 
#356.604   0.047 356.820

这个简单的改变产生了大约 6.5 倍的加速。

[原来的做法]

可能还有一些改进的余地，但以下在运行 64 位 R 的 Linux VM（2.6 GHz 处理器）上花费了大约 25 分钟：

library(data.table)
library(RcppArmadillo)
##
dt[
  ,c("lm_b0","lm_b1") := as.list(
    unname(fastLm(Y ~ x, data=.SD)$coefficients)),
  by=user_id]
##
R> dt[c(1:2, 31:32, 61:62),]
   user_id   x         Y     lm_b0    lm_b1
1:       1 1.0 1674.8316 -202.0066 744.6252
2:       1 1.5  369.8608 -202.0066 744.6252
3:       2 1.0  463.7460 -144.2961 374.1995
4:       2 1.5  412.7422 -144.2961 374.1995
5:       3 1.0  513.0996  217.6442 261.0022
6:       3 1.5 1140.2766  217.6442 261.0022

数据：

dt <- data.table(
  user_id = rep(1:500000,each=30))
##
dt[, x := seq(1, by=.5, length.out=30), by = user_id]
dt[, Y := 1000*runif(1)*x, by = user_id]
dt[, Y := Y + rnorm(
  30, 
  mean = sample(c(-.05,0,0.5)*mean(Y),1), 
  sd = mean(Y)*.25), 
  by = user_id]

【讨论】：

真的需要 25 分钟，而不是几秒钟吗？
@BrodieG 是的，我也对此感到惊讶。我假设这是因为我在资源有限的虚拟机上运行它，而不是在适当的（4+ 核心、8 GB RAM 等）台式机上运行。顺便说一句，我喜欢你跳出框框思考的方法，+1
这是一种巧妙的方法。
新手错误 :) 解析公式比 fastLm 节省的时间多。使用显式 y 和 X 作为向量和矩阵重试。
添加公式界面是喜忧参半。我们希望它是为了界面的一致性，而不是为了扰乱性能……顺便说一句，对答案的更新很好。

【解决方案4】：

您可以像这样使用 data.table 来尝试一下。我刚刚创建了一些玩具数据，但我想 data.table 会有所改进。这是相当快的。但这是一个相当大的数据集，因此也许可以在较小的样本上对这种方法进行基准测试，看看速度是否要好得多。祝你好运。


    library(data.table)

    exp <- data.table(id = rep(c("a","b","c"), each = 20), x = rnorm(60,5,1.5), y = rnorm(60,2,.2))
    # edit: it might also help to set a key on id with such a large data-set
    # with the toy example it would make no diff of course
    exp <- setkey(exp,id)
    # the nuts and bolts of the data.table part of the answer
    result <- exp[, as.list(coef(lm(y ~ x))), by=id]
    result
       id (Intercept)            x
    1:  a    2.013548 -0.008175644
    2:  b    2.084167 -0.010023549
    3:  c    1.907410  0.015823088

【讨论】：

【解决方案5】：

使用 Rfast 的示例。

假设单个响应和 500K 预测变量。

y <- rnorm(30)
x <- matrnorm(500*1000,30)
system.time( Rfast::univglms(y, x,"normal") )  ## 0.70 seconds

假设 500K 响应变量和单个预测变量。

system.time( Rfast::mvbetas(x,y) )  ## 0.60 seconds

注意：上述时间将在不久的将来减少。

【讨论】：