在 R 中估计具有数百万个观测值和数千个变量的 OLS 模型
Posted
技术标签:
【中文标题】在 R 中估计具有数百万个观测值和数千个变量的 OLS 模型【英文标题】:Estimating an OLS model in R with million observations and thousands of variables 【发布时间】:2019-04-04 20:05:04 【问题描述】:我正在尝试使用 biglm 估计具有约 100 万个观察值和约 50,000 个变量的大型 OLS 回归。
我计划使用每个大约 100 个观察值的块来运行每个估计。我用一个小样本测试了这个策略,效果很好。
但是,在尝试为 biglm 函数定义公式时,我收到了“错误:protect():保护堆栈溢出”的真实数据。
我已经试过了:
从 R 开始 --max-ppsize=50000
设置选项(表达式 = 50000)
但错误仍然存在
我在 Windows 上工作并使用 Rstudio
# create the sample data frame (In my true case, I simply select 100 lines from the original data that contains ~1,000,000 lines)
DF <- data.frame(matrix(nrow=100,ncol=50000))
DF[,] <- rnorm(100*50000)
colnames(DF) <- c("y", paste0("x", seq(1:49999)))
# get names of covariates
my_xvars <- colnames(DF)[2:( ncol(DF) )]
# define the formula to be used in biglm
# HERE IS WHERE I GET THE ERROR :
my_f <- as.formula(paste("y~", paste(my_xvars, collapse = " + ")))
编辑 1: 我练习的最终目标是估计所有 50,000 个变量的平均效果。因此,简化模型选择更少的变量并不是我现在正在寻找的解决方案。
【问题讨论】:
您希望从 50,000 个变量模型中得到什么?在建模之前可能值得研究降维技术。 您可能想要使用 OLS 以外的其他方法,例如可以比 OLS 更好地应对许多(相关)预测变量的惩罚回归 sthda.com/english/articles/37-model-selection-essentials-in-r/… 我练习的最终目标确实是估计因变量中 50,000 个协变量中每一个的平均影响。因此,用更少的变量来简化模型并不是一个理想的解决方案。我已经使用大块的行来分割估计。如果有任何方法也可以使用列子集来分离估计,然后能够恢复对每个变量的影响的无偏估计,那就太棒了,但是我还找不到类似的东西。 这似乎无效,因为每个变量的系数还取决于模型中的其他变量。因此,与所有变量的一个模型的平均系数相比,2 个不同模型的平均系数可能会有很大差异。 @igoR87 你是绝对正确的,特征空间分区似乎比样本空间分区复杂得多。我从一个杜克大学的学生那里找到了下面的论文,该论文为这类问题提出了一套解决方案。但是,在尝试实现(和理解)类似的东西之前,我希望也许已经有一个更容易获得的解决方案dukespace.lib.duke.edu/dspace/bitstream/handle/10161/12912/… 【参考方案1】:第一个瓶颈(我不能保证不会有其他瓶颈)在于公式的构建。 R 无法从文本中构造出这么长的公式(细节太难看,现在无法探索)。下面我展示了biglm
代码的破解版本,它可以直接获取模型矩阵X
和响应变量y
,而不是使用公式来构建它们。 然而:下一个瓶颈是内部函数biglm:::bigqr.init()
(在biglm
内部调用)尝试分配大小为choose(nc,2)=nc*(nc-1)/2
的数字向量(其中nc
是列数. 当我尝试使用 50000 列时,我得到了
错误:无法分配大小为 9.3 Gb 的向量
(nc
为 25000 时需要 2.3Gb)。下面的代码在nc <- 10000
时在我的笔记本电脑上运行。
我对这种方法有一些注意事项:
由于上述问题,除非您有至少 10G 的内存,否则您将无法处理具有 50000 列的问题。biglm:::update.biglm
必须以并行方式修改(这应该不会太难)
我不知道 p>>n 问题(适用于拟合初始块的级别)是否会咬你。在下面运行我的示例时(10 行,10000 列),除了 10 个参数之外,所有参数都是NA
。我不知道这些NA
值是否会污染结果,从而导致连续更新失败。如果是这样,我不知道是否有办法解决这个问题,或者它是否是基本的(所以你至少需要nr>nc
来进行初始拟合)。 (直接做一些小实验看看有没有问题,但我已经在这上面花了太长时间了……)
不要忘记,使用这种方法,您必须在模型矩阵中显式添加一个截距列(例如,X <- cbind(1,X)
,如果需要的话。
示例(先将底部的代码保存为my_biglm.R
):
nr <- 10
nc <- 10000
DF <- data.frame(matrix(rnorm(nr*nc),nrow=nr))
respvars <- paste0("x", seq(nc-1))
names(DF) <- c("y", respvars)
# illustrate formula problem: fails somewhere in 15000 < nc < 20000
try(reformulate(respvars,response="y"))
source("my_biglm.R")
rr <- my_biglm(y=DF[,1],X=as.matrix(DF[,-1]))
my_biglm <- function (formula, data, weights = NULL, sandwich = FALSE,
y=NULL, X=NULL, off=0)
if (!is.null(weights))
if (!inherits(weights, "formula"))
stop("`weights' must be a formula")
w <- model.frame(weights, data)[[1]]
else w <- NULL
if (is.null(X))
tt <- terms(formula)
mf <- model.frame(tt, data)
if (is.null(off <- model.offset(mf)))
off <- 0
mm <- model.matrix(tt, mf)
y <- model.response(mf) - off
else
## model matrix specified directly
if (is.null(y)) stop("both y and X must be specified")
mm <- X
tt <- NULL
qr <- biglm:::bigqr.init(NCOL(mm))
qr <- biglm:::update.bigqr(qr, mm, y, w)
rval <- list(call = sys.call(), qr = qr, assign = attr(mm,
"assign"), terms = tt, n = NROW(mm), names = colnames(mm),
weights = weights)
if (sandwich)
p <- ncol(mm)
n <- nrow(mm)
xyqr <- bigqr.init(p * (p + 1))
xx <- matrix(nrow = n, ncol = p * (p + 1))
xx[, 1:p] <- mm * y
for (i in 1:p) xx[, p * i + (1:p)] <- mm * mm[, i]
xyqr <- update(xyqr, xx, rep(0, n), w * w)
rval$sandwich <- list(xy = xyqr)
rval$df.resid <- rval$n - length(qr$D)
class(rval) <- "biglm"
rval
【讨论】:
非常感谢!解决了公式限制。但是,正如您所料,我遇到了 50k vars 的 RAM 限制。关于您的第 1 点和第 3 点:-我没有 10G 可用-关于块内的 p>n,我用我的数据测试了问题的简化版本(1k 行,200 个变量,50 行/块),并返回了分块估计与一次估计的结果相同。但是,我的真实数据是稀疏的(> 80% 零)。然后我用密集矩阵对其进行了测试,结果非常不同,分块估计返回了很多 NA 我将寻找解决我的问题的替代解决方案以上是关于在 R 中估计具有数百万个观测值和数千个变量的 OLS 模型的主要内容,如果未能解决你的问题,请参考以下文章
如何使用超过 5000 万个观测值的样本计算具有固定效应的 logit 模型的边际效应
在一个线程上删除具有数百万个字符串的大型哈希图会影响另一个线程的性能
升级到 xCode 4 后,我收到数百万个“此类不符合键值编码”错误