我正在尝试使用 biglm 估计具有约 100 万个观测值和约 50,000 个变量的大型 OLS 回归。
我计划使用每个大约 100 个观察的块来运行每个估计。我用一个小样本测试了这个策略,效果很好。
但是,使用真实数据时,我在尝试定义 biglm 函数的公式时收到“错误:保护():保护堆栈溢出”。
我已经试过了:
以 --max-ppsize=50000 开始 R
设置选项(表达式 = 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 个变量的平均影响。因此,简化模型选择更少的变量并不是我现在正在寻找的解决方案。
第一个瓶颈(我不能保证不会有其他瓶颈)是在公式的构建中。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
。
我对这种方法有一些警告:
biglm:::update.biglm
会以平行的方式(这不应该是太硬)进行修改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
}
本文收集自互联网,转载请注明来源。
如有侵权,请联系[email protected] 删除。
我来说两句