Оценка модели OLS в R с миллионами наблюдений и тысячами переменных - PullRequest
Новая
       33

Оценка модели OLS в R с миллионами наблюдений и тысячами переменных

1 голос
/ 04 апреля 2019

Я пытаюсь оценить большую регрессию МНК с помощью ~ 1 миллиона наблюдений и ~ 50 000 переменных с использованием biglm.

Я планирую провести каждую оценку, используя куски приблизительно по 100 наблюдений каждый. Я проверил эту стратегию с небольшим образцом, и она работала нормально.

Однако с реальными данными я получаю сообщение «Ошибка: защита (): переполнение стека защиты» при попытке определить формулу для функции biglm.

Я уже пробовал:

  • начиная с 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 переменных. Поэтому упрощение модели выбора меньшего количества переменных не является решением, которое я сейчас ищу.

1 Ответ

0 голосов
/ 05 апреля 2019

Первое узкое место (я не могу гарантировать, что других не будет) - это построение формулы.R не может построить длинную формулу из текста (детали слишком уродливы, чтобы их изучать прямо сейчас).Ниже я показываю взломанную версию кода biglm, которая может принимать матрицу модели X и переменную ответа y напрямую, а не использовать формулу для их построения. Однако : следующее узкое место заключается в том, что внутренняя функция biglm:::bigqr.init(), которая вызывается внутри biglm, пытается выделить числовой вектор размером choose(nc,2)=nc*(nc-1)/2 (где nc - количество столбцов.Когда я пытаюсь использовать 50000 столбцов, я получаю

Ошибка: невозможно выделить вектор размером 9,3 ГБ

(требуется 2,3 ГБ, если nc равно 25000).ниже работает на моем ноутбуке, когда nc <- 10000.

У меня есть несколько предостережений об этом подходе:

  • вы не сможете обработать пробелм с 50000 столбцами, если у вас нетпо крайней мере 10 ГБ памяти из-за описанной выше проблемы.
  • 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
}
...