估算R中的OLS模型，包含数百万个观测值和数千个变量

Question

我试图估计一个大的OLS回归与大约100万观测和~50,000变量使用biglm。

我计划使用每个大约100个观察的块来运行每个估计。我用一个小样本测试了这个策略，它运行良好。

但是，对于真实数据，我在尝试定义biglm函数的公式时遇到“错误：protect（）：保护堆栈溢出”。

我已经尝试过了：

• 以-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个变量的平均效果。因此，简化模型选择较少的变量并不是我现在正在寻找的解决方案。

Answer 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
}