如何让 group_by 和 lm 变快？
How to make group_by and lm fast?
这是一个示例。
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的速度而是很多subject和subsetting.
理论上
首先，您可以。
其次，显式数据拆分并不是分组回归的唯一方法（或推荐方法）。参见 and 。因此，将模型构建为 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)，这样第一列是截距，第二列是斜率。
我的想法：我们可能需要更好的包来进行大回归
在这里使用 glm4 并没有比 更有吸引力，因为它基本上也只是告诉你回归系数。它也比 data.table 方法慢一点，因为它计算拟合值和残差。
简单回归分析不需要适当的模型拟合例程。我有一些关于如何在这种回归上做花哨的事情的答案，比如 ，其中还给出了如何计算所有统计数据的详细信息。但是当我写这个答案时，我正在考虑一些关于大回归问题的一般性问题。我们可能需要更好的包，否则个案编码就没完没了。
回复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.
更新：现在可以使用快速解决方案
查看 。
由于 OP 似乎只寻找 beta，这里有一种使用 data.table 包来仅计算 beta 的方法。公式见参考资料。
dt[, sumx := sum(day), by=.(subject)][,
denom := sum(day^2) - sumx^2 / .N, by=.(subject)]
dt[, lapply(.SD, function(y) (sum(day*y) - (sumx[1L] * sum(y))/.N) / denom[1L]),
by=.(subject),
.SDcols = paste0("y", 1:4)]
数据：
library(data.table)
set.seed(0L)
nSubj <- 200e3
nr <- 1e6
dt <- data.table(
subject = rep(1:nSubj, each=5),
day = 3:7,
y1 = runif(nr),
y2 = rpois(nr, 3),
y3 = rnorm(nr),
y4 = rnorm(nr, 1, 5))
dt2 <- copy(dt)
参考：
Inference in Linear Regression
在对线性回归的智力理解上，李哲元的解法最好(+1)
添加一些时间：
system.time({
dt[, lapply(.SD, function(y) cov(x,y) / var(x) ),
by=.(subject),
.SDcols=paste0("y", 1:4)]
})
user system elapsed
73.96 0.00 74.15
system.time({
dt2[, sumx := sum(day), by=.(subject)][,
denom := sum(day^2) - sumx^2 / .N, by=.(subject)]
dt2[, lapply(.SD, function(y) (sum(day*y) - (sumx[1L] * sum(y))/.N) / denom[1L]),
by=.(subject),
.SDcols = paste0("y", 1:4)]
})
user system elapsed
2.14 0.00 2.14
这是一个示例。
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的速度而是很多subject和subsetting.
理论上
首先，您可以
其次，显式数据拆分并不是分组回归的唯一方法（或推荐方法）。参见 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)，这样第一列是截距，第二列是斜率。
我的想法：我们可能需要更好的包来进行大回归
在这里使用 glm4 并没有比 data.table 方法慢一点，因为它计算拟合值和残差。
简单回归分析不需要适当的模型拟合例程。我有一些关于如何在这种回归上做花哨的事情的答案，比如
回复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.
更新：现在可以使用快速解决方案
查看
由于 OP 似乎只寻找 beta，这里有一种使用 data.table 包来仅计算 beta 的方法。公式见参考资料。
dt[, sumx := sum(day), by=.(subject)][,
denom := sum(day^2) - sumx^2 / .N, by=.(subject)]
dt[, lapply(.SD, function(y) (sum(day*y) - (sumx[1L] * sum(y))/.N) / denom[1L]),
by=.(subject),
.SDcols = paste0("y", 1:4)]
数据：
library(data.table)
set.seed(0L)
nSubj <- 200e3
nr <- 1e6
dt <- data.table(
subject = rep(1:nSubj, each=5),
day = 3:7,
y1 = runif(nr),
y2 = rpois(nr, 3),
y3 = rnorm(nr),
y4 = rnorm(nr, 1, 5))
dt2 <- copy(dt)
参考： Inference in Linear Regression
在对线性回归的智力理解上，李哲元的解法最好(+1)
添加一些时间：
system.time({
dt[, lapply(.SD, function(y) cov(x,y) / var(x) ),
by=.(subject),
.SDcols=paste0("y", 1:4)]
})
user system elapsed
73.96 0.00 74.15
system.time({
dt2[, sumx := sum(day), by=.(subject)][,
denom := sum(day^2) - sumx^2 / .N, by=.(subject)]
dt2[, lapply(.SD, function(y) (sum(day*y) - (sumx[1L] * sum(y))/.N) / denom[1L]),
by=.(subject),
.SDcols = paste0("y", 1:4)]
})
user system elapsed
2.14 0.00 2.14