如何快速制作group_by和lm？

这是一个样本.

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速度,而是很多主题和子集.

理论上

首先,您可以使用多个LHS拟合线性模型.

其次,显式数据拆分不是逐组回归的唯一方式(或推荐方式).参见R回归分析:分析某个种族的数据和R:为每个类别建立单独的模型.因此,建立你的模型作为cbind(x1, x2, x3, x4) ~ day * subject其中subject的一个因素变量.

最后,由于您有许多因子级别并且使用大型数据集,因此lm是不可行的.考虑使用speedglm::speedlmwith sparse = TRUE或MatrixModels::glm4with 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),所以第一列是截距,第二列是斜率.

我的想法:我们可能需要更好的大回归包

glm4在这里使用并不比chinsoon12的data.table解决方案具有吸引力的优势,因为它基本上只是告诉你回归系数.它也比data.table方法慢一点,因为它计算拟合值和残差.

简单回归分析不需要适当的模型拟合程序.关于如何在这种回归中做一些奇特的东西,我有一些答案,比如数据框中所有变量之间的快速成对简单线性回归,其中还给出了有关如何计算所有统计数据的详细信息.但是当我写这个答案时,我正在考虑关于大回归问题的一般事情.我们可能需要更好的软件包,否则就无法完成个案编码.

回复OP

speedglm::speedlm(x1 ~ 0 + subject + day:subject, data = df, sparse = TRUE)给出错误:无法分配大小为74.5 Gb的向量

是的,因为它的sparse逻辑很糟糕.

MatrixModels::glm4(x1 ~ day * subject, data = df, sparse = TRUE)在Cholesky中给出错误(crossprod(from),LDL = FALSE):internal_chm_factor:Cholesky分解失败

这是因为你只有一个数据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))

级别"c" f只有一个基准/行.

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因子分解无法解决排名不足的模型.如果我们使用lmQR分解,我们将看到一个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.

更新:快速解决方案现已推出

查看快速分组 - 简单线性回归和一般配对简单线性回归.