r - 如何使group_by和lm快速？

Question

这是一个示例。

df <- tibble(
      subject = rep(letters[1:7], c(5, 6, 7, 5, 2, 5, 2)),
      day = c(3:7, 2:7, 1:7, 3:7, 6:7, 3:7, 6:7),
      x1 = runif(32), x2 = rpois(32, 3), x3 = rnorm(32), x4 = rnorm(32, 1, 5))

df %>%
  group_by(subject) %>%
  summarise(
    coef_x1 = lm(x1 ~ day)$coefficients[2],
    coef_x2 = lm(x2 ~ day)$coefficients[2],
    coef_x3 = lm(x3 ~ day)$coefficients[2],
    coef_x4 = lm(x4 ~ day)$coefficients[2])

此数据很小，因此性能不是问题。

但是我的数据是如此之大，大约有1,000,000行和200,000个主题，并且此代码非常慢。

我认为原因不是lm的速度，而是很多主题和设置了的原因。

Answer 1

理论上

首先，您可以fit a linear model with multiple LHS。

其次，显式数据拆分不是分组回归的唯一方法(或推荐的方法)。参见R regression analysis: analyzing data for a certain ethnicity和R: build separate models for each category。因此，将模型构建为cbind(x1, x2, x3, x4) ~ day * subject，其中subject是一个因子变量。

最后，由于您具有许多因子级别并使用大型数据集，因此lm是不可行的。考虑将speedglm::speedlm与sparse = TRUE一起使用，或MatrixModels::glm4与sparse = TRUE一起使用。

事实上
speedlm和glm4均未处于积极开发中。 (在我看来)它们的功能是原始的。
speedlm和glm4都不支持多个LHS作为lm。因此，您需要将4个单独的模型x1 ~ day * subject改为x4 ~ day * subject。

这两个软件包在sparse = TRUE后面具有不同的逻辑。

speedlm首先使用标准model.matrix.default构造一个密集的设计矩阵，然后使用is.sparse来检查其是否稀疏。如果为TRUE，则后续计算可以使用稀疏方法。

glm4使用model.Matrix构造设计矩阵，并且可以直接构建稀疏矩阵。

因此，在这个稀疏性问题中speedlm和lm一样糟糕并不奇怪，并且glm4是我们真正想要使用的。
glm4没有用于分析拟合模型的完整有用的通用函数集。您可以通过coef，fitted和residuals提取系数，拟合值和残差，但必须自己计算所有统计信息(标准误差，t统计量，F统计量等)。对于那些非常了解回归理论的人来说，这并不是什么大不了的事情，但是仍然很不方便。
glm4仍然希望您使用最佳的模型公式，以便可以构造最稀疏的矩阵。传统的~ day * subject确实不是一个好方法。以后我可能应该对此问题进行问答。基本上，如果您的公式具有截距并且将因素进行了对比，那么您将失去稀疏性。这是我们应该使用的一种:~ 0 + subject + day:subject。

用glm4进行测试

## use chinsoon12's data in his answer
set.seed(0L)
nSubj <- 200e3
nr <- 1e6
DF <- data.frame(subject = gl(nSubj, 5),
                 day = 3:7,
                 y1 = runif(nr), 
                 y2 = rpois(nr, 3), 
                 y3 = rnorm(nr), 
                 y4 = rnorm(nr, 1, 5))

library(MatrixModels)
fit <- glm4(y1 ~ 0 + subject + day:subject, data = DF, sparse = TRUE)

在我的1.1GHz Sandy Bridge笔记本电脑上大约需要6到7秒。让我们提取其系数:

b <- coef(fit)

head(b)
#  subject1   subject2   subject3   subject4   subject5   subject6 
# 0.4378952  0.3582956 -0.2597528  0.8141229  1.3337102 -0.2168463 

tail(b)
#subject199995:day subject199996:day subject199997:day subject199998:day 
#      -0.09916175       -0.15653402       -0.05435883       -0.02553316 
#subject199999:day subject200000:day 
#       0.02322640       -0.09451542 

您可以执行B <- matrix(b, ncol = 2)，以使第一列为intercept，第二列为slope。

我的想法:对于大的回归，我们可能需要更好的软件包

在这里使用glm4并没有比chinsoon12's data.table solution更具吸引力，因为它基本上也只是告诉您回归系数。它也比data.table方法慢一点，因为它会计算拟合值和残差。

简单回归分析不需要适当的模型拟合例程。我对如何在这种回归上做一些花哨的事情有一些答案，例如Fast pairwise simple linear regression between all variables in a data frame，其中也给出了如何计算所有统计信息的详细信息。但是，当我写这个答案时，我在思考关于大型回归问题的一般性问题。我们可能需要更好的程序包，否则将无法进行个案编码。

回复OP

speedglm::speedlm(x1 ~ 0 + subject + day:subject, data = df, sparse = TRUE) gives Error: cannot allocate vector of size 74.5 Gb

是的，因为它的sparse逻辑不好。

MatrixModels::glm4(x1 ~ day * subject, data = df, sparse = TRUE) gives Error in Cholesky(crossprod(from), LDL = FALSE) : internal_chm_factor: Cholesky factorization failed

这是因为某些subject只有一个基准。您至少需要两个数据才能拟合一条线。这是一个示例(在密集设置中):

dat <- data.frame(t = c(1:5, 1:9, 1),
                  f = rep(gl(3,1,labels = letters[1:3]), c(5, 9, 1)),
                  y = rnorm(15))

f的级别“c”仅具有一个基准/行。

X <- model.matrix(~ 0 + f + t:f, dat)
XtX <- crossprod(X)
chol(XtX)
#Error in chol.default(XtX) : 
#  the leading minor of order 6 is not positive definite

Cholesky分解无法解决等级不足的模型。如果使用lm的QR因式分解，我们将看到NA系数。

lm(y ~ 0 + f + t:f, dat)
#Coefficients:
#      fa        fb        fc      fa:t      fb:t      fc:t  
# 0.49893   0.52066  -1.90779  -0.09415  -0.03512        NA  

我们只能估算出水平“c”的截距，而不是斜率。

请注意，如果使用data.table解决方案，则在计算此级别的斜率时最终会得到0 / 0，最终结果是NaN。

更新:现已提供快速解决方案

checkout Fast group-by simple linear regression and general paired simple linear regression。