我有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.
是否有快速的方法按组进行线性回归？
Bro*_*ieG 20
您可以使用基本公式计算斜率和回归. lm如果您关心的是这两个数字,那么会做很多不必要的事情.这里我data.table用于聚合,但你也可以在基数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
]
)
产量为500K用户〜每个30个盲(秒):
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)
)
Gre*_*gor 11
如果您只想要系数,我只会将其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)
我们不适用全局intercept(+ 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
因此,对于a水平id,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倍的加速.
[原创方法]
可能还有一些改进空间,但是在运行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]
你可以尝试使用像这样的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