按组快速线性回归
Fast linear regression by group
我有 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 不兼容。
有没有快速的分组线性回归的方法？
您可以像这样使用 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
如果你想要的只是系数，我会使用 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 即可几分钟。
更新：
正如德克所指出的，我原来的方法可以通过直接指定 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 倍的加速。
[原做法]
可能还有一些改进的余地，但以下内容在 Linux VM（2.6 GHz 处理器），运行 64 位 R 上花费了大约 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]
您可以只使用计算斜率和回归的基本公式。如果您只关心这两个数字，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 微秒。
并确认这是否按预期工作：
> 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)
)
一个使用 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
注：以上时间近期会减少。
我有 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 不兼容。
有没有快速的分组线性回归的方法？
您可以像这样使用 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
如果你想要的只是系数，我会使用 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 即可几分钟。
更新：
正如德克所指出的，我原来的方法可以通过直接指定 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 倍的加速。
[原做法]
可能还有一些改进的余地，但以下内容在 Linux VM（2.6 GHz 处理器），运行 64 位 R 上花费了大约 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]
您可以只使用计算斜率和回归的基本公式。如果您只关心这两个数字，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 微秒。
并确认这是否按预期工作：
> 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)
)
一个使用 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
注：以上时间近期会减少。