我正在尝试在 R 中运行固定效应回归。当我在没有应用固定效应因子的情况下运行线性模型时,模型工作得很好。但是当我应用因子 - 这是用户 ID 的数字代码时,我收到以下错误:
Error in rep.int(c(1, numeric(n)), n - 1L) : cannot allocate vector of length 1055470143
我不确定错误是什么意思,但我担心这可能是在 R 中正确编码变量的问题。
我正在尝试在 R 中运行固定效应回归。当我在没有应用固定效应因子的情况下运行线性模型时,模型工作得很好。但是当我应用因子 - 这是用户 ID 的数字代码时,我收到以下错误:
Error in rep.int(c(1, numeric(n)), n - 1L) : cannot allocate vector of length 1055470143
我不确定错误是什么意思,但我担心这可能是在 R 中正确编码变量的问题。
我认为这是更多的统计问题和更少的编程问题,原因有两个:
首先,我不确定您使用的是横截面数据还是面板数据。如果您使用横截面数据,控制 30000 个人是没有意义的(当然,它们会增加变化)。
其次,如果您使用的是面板数据,那么plm
R 中的 package 等很好的包可以进行这种计算。
一个例子:
set.seed(42)
DF <- data.frame(x=rnorm(1e5),id=factor(sample(seq_len(1e3),1e5,TRUE)))
DF$y <- 100*DF$x + 5 + rnorm(1e5,sd=0.01) + as.numeric(DF$id)^2
fit <- lm(y~x+id,data=DF)
R 会话需要将近 2.5 GB 的 RAM(如果您添加操作系统所需的 RAM,这比许多 PC 都可用)并且需要一些时间才能完成。结果很没用。
如果您没有遇到 RAM 限制,您可能会受到向量长度的限制(例如,如果您有更多的因子级别),特别是如果您使用旧版本的 R。
怎么了?
第一步lm
是使用函数创建设计矩阵model.matrix
。这是一个较小的示例,说明因子会发生什么:
model.matrix(b~a,data=data.frame(a=factor(1:5),b=2))
# (Intercept) a2 a3 a4 a5
# 1 1 0 0 0 0
# 2 1 1 0 0 0
# 3 1 0 1 0 0
# 4 1 0 0 1 0
# 5 1 0 0 0 1
# attr(,"assign")
# [1] 0 1 1 1 1
# attr(,"contrasts")
# attr(,"contrasts")$a
# [1] "contr.treatment"
看看 n 个因子水平如何导致 n-1 个虚拟变量?如果您有许多因子水平和许多观察值,则该矩阵会变得很大。
你该怎么办?
我很确定,您应该使用混合效果模型。有两个重要的包在 R 中实现线性混合效果模型,包 nlme 和更新的包 lme4。
library(lme4)
fit.mixed <- lmer(y~x+(1|id),data=DF)
summary(fit.mixed)
Linear mixed model fit by REML
Formula: y ~ x + (1 | id)
Data: DF
AIC BIC logLik deviance REMLdev
1025277 1025315 -512634 1025282 1025269
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 8.9057e+08 29842.472
Residual 1.3875e+03 37.249
Number of obs: 100000, groups: id, 1000
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.338e+05 9.437e+02 353.8
x 1.000e+02 1.180e-01 847.3
Correlation of Fixed Effects:
(Intr)
x 0.000
这需要很少的 RAM,计算速度很快,并且是一个更正确的模型。
看看随机截距是如何解释大部分方差的?
所以,你需要研究混合效应模型。有一些不错的出版物,例如Baayen、Davidson、Bates (2008),解释了如何使用 lme4。