Большие фиксированные эффекты биномиальной регрессии в R

Question

Мне нужно запустить логистическую регрессию на относительно большом фрейме данных с 480000 записей с 3 фиксированными переменными эффекта. Фиксированный эффект вар A имеет 3233 уровня, вар B имеет 2326 уровней, вар C имеет 811 уровней. В общем, у меня есть 6370 фиксированных эффектов. Данные являются поперечными. Если я не могу запустить эту регрессию с помощью нормального glm функция, потому что матрица регрессии кажется слишком большой для моей памяти (я получаю сообщение " Error: cannot allocate vector of size 22.9 Gb "). Я ищу альтернативные способы запуска этой регрессии на моем Macbook Air (OS X 10.9.5 8 ГБ ОЗУ). У меня также есть доступ к серверу с 16 ГБ ОЗУ.

Я пытался решить проблему несколькими различными способами, но пока ни один из них не привел к удовлетворительным результатам:

lfe / felm: использование функции регрессии фелма lfe пакет, который вычитает фиксированные эффекты перед запуском регрессии. Это прекрасно работает и позволило мне запустить описанную выше регрессию как обычную линейную модель всего за несколько минут. Тем не мение, lfe не поддерживает логистические регрессии и glms. Таким образом, felm отлично подходил для того, чтобы получить представление о подгонке модели для разных моделей, но не работает для окончательных моделей логистической регрессии.

biglm / bigglm: я думал об использовании bigglm разбить мою функцию на более управляемые куски. Тем не менее, несколько источников (например, link1, link2, link3) упоминают, что для того, чтобы это работало, уровни факторов должны быть согласованы между порциями, то есть каждый блок должен содержать по меньшей мере один из каждого фактора каждой переменной фактора. Факторы A и B содержат уровни, которые появляются только один раз, поэтому я не могу разбить наборы на разные куски с одинаковыми уровнями. Если я удаляю 10 факторов с фиксированным эффектом A и 8 факторов B (незначительное изменение), у меня останутся только факторы с 4+ уровнями, и разделение моих данных на 4 части сделает их намного более управляемыми. Однако тогда мне все еще нужно выяснить, как отсортировать мою df таким образом, чтобы гарантировать, что мои 480.000 записей будут отсортированы по 4 блокам, в которых каждый факторный уровень каждого из 3 факторов появляется хотя бы один раз.

GlmmGS / glmgs: glmmgs Функция в пакете с тем же именем выполняет вычитание с фиксированными эффектами, как lfe пакет для логистических регрессий с использованием алгоритма Гаусса-Зейделя. К сожалению, пакет больше не разрабатывается. Будучи относительно новым для R и не имея глубокого опыта со статистикой, я не могу понять результат и не имею представления о том, как преобразовать его таким образом, чтобы получить нормальный "размер эффекта", "подбор модели", " "Интервал значимости" - индикаторы, которые дают резюме регрессии ГЛМ

Я отправил сообщение авторам пакета. Они любезно ответили следующим образом:

Пакет не обеспечивает вывод в том же формате объекта glm. Тем не менее, вы можете легко вычислить большую часть статистики подгонки (стандартная ошибка оценок, достоверность подгонки), учитывая текущий вывод (в версии CRAN, я считаю, что текущий вывод представляет собой вектор оценки коэффициентов и связанный вектор стандартных ошибок, то же самое для ковариационных компонентов, но вам не нужно беспокоиться о них, если вы подходите модели без случайных эффектов). Только имейте в виду, что ковариационные матрицы, используемые для вычисления стандартных ошибок, являются обратными диагональным блокам прецизионной матрицы, связанной с алгоритмом Гаусса-Зейделя, и поэтому они имеют тенденцию недооценивать стандартные ошибки совместной вероятности. Я больше не поддерживаю пакет, и у меня нет времени вдаваться в конкретные детали; основополагающая теория, лежащая в основе пакета, может быть найдена в статье, на которую есть ссылки в руководстве, все остальное вы должны решить с помощью ручки и бумаги:).

Если кто-то может объяснить, как "легко рассчитать большую часть подходящей статистики" таким образом, чтобы кто-то без какого-либо образования в области статистики мог это понять (возможно, было бы невозможно) или предоставить код R, который показывает на примере того, как это можно сделать, я был бы Весьма признателен!

Revolution Analytics: Я установил Revolution Analytics Enterprise на виртуальной машине, которая имитирует Windows 7 на моем Mac. Программа имеет функцию под названием RxLogit это оптимизировано для больших логистических регрессий. С использованием RxLogit функция, которую я получаю the error (Failed to allocate 326554568 bytes. Error in rxCall("RxLogit", params) : bad allocation), так что эта функция, кажется, слишком сталкивается с проблемами памяти. Однако программное обеспечение позволяет мне выполнять регрессию на кластере распределенных вычислений. Так что я мог бы просто "убить проблему", купив вычислительное время в кластере с большим объемом памяти. Тем не менее, мне интересно, если программа революционной аналитики предоставляет какие-либо формулы или методы, которые я не знаю, которые позволили бы мне делать какие-то lfe -подобная операция вычитания или bigglm -подобная операция чанкинга, учитывающая факторы.

MatrixModels / glm4: один человек предложил мне использовать glm4 функция MatrixModels пакет с sparse = TRUE атрибут для ускорения расчета. Если я бегу glm4 регрессия со всеми фиксированными эффектами я получаю "Error in Cholesky(crossprod(from), LDL = FALSE) : internal_chm_factor: Cholesky factorization failed "ошибка. Если я запускаю его только с фиксированными переменными эффекта B ИЛИ A и C, расчет работает и возвращает "glpModel" объект. Как с glmmGS У меня есть некоторые проблемы в превращении этого вывода в форму, которая имеет смысл для меня, так как стандарт summary() метод, похоже, не работает на нем.

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

Answer 1

Проверять, выписываться

glmmboot{glmmML}

http://cran.r-project.org/web/packages/glmmML/glmmML.pdf

Есть также хороший документ Брострома и Холмберга ( http://cran.r-project.org/web/packages/eha/vignettes/glmmML.pdf)

Вот пример из их документа:

dat <- data.frame(y = rbinom(5000, size = 1, prob = 0.5),
               x = rnorm(5000), group = rep(1:1000, each = 5))
fit1 <- glm(y ~ factor(group) + x, data = dat, family = binomial)

require(glmmML)
fit2 <- glmmboot(y ~ x, cluster = group,data = dat)

Разница во времени вычислений "огромна"!

Answer 2

Для потомков, я также хотел бы порекомендовать пакет speedglm, который я нашел полезным при попытке выполнить логистическую регрессию на больших наборах данных. Кажется, он использует примерно вдвое меньше памяти и заканчивается быстрее, чем glm,

Answer 3

Я согласен с тем, кто (@Ben Bolker, я полагаю?) Предложил вам использовать glm4 функция от MatrixModels, Во-первых, это решает проблему с памятью, если вы используете sparse аргумент. Матрица плотного проектирования с 480 000 записей и 6370 фиксированных эффектов потребует 6371 * 480 000 * 8 = 24,464,640,000 байтов. Тем не менее, ваша матрица дизайна будет очень разреженной (много нулей), поэтому вы можете использовать гораздо меньшую (в памяти) матрицу дизайна, если будете использовать разреженную. Во-вторых, вы можете использовать разреженность, чтобы ускорить оценку.

Что касается параметров, быстрый поиск показывает, что speedglm также имеет sparse аргумент хотя я не пробовал. Ключевым моментом для любого метода, которым вы заканчиваете, является то, что он должен использовать то, что ваша матрица проектирования разрежена как для сокращения времени вычислений, так и для уменьшения требований к памяти.

Вы получаете ошибку (Error in Cholesky(crossprod(from), LDL = FALSE) : internal_chm_factor: Cholesky factorization failed" error) вероятно, потому что ваша матрица дизайна является единственной. В этом случае ваша проблема не имеет уникального решения, и некоторые варианты состоят в том, чтобы объединить некоторые из групповых уровней, использовать модель штрафов или случайных эффектов.

Вы правы в том, что, похоже, не существует метода glpModel учебный класс. Тем не менее, кажется, что слоты имеют очевидные имена, и вам не понадобится много времени, чтобы получить, например, стандартные ошибки в вашем оценщике, вычислить оценку отклонения и т. Д.

Answer 4

Функция feglmиз пакета альпака может быть то, что вам нужно. Он позволяет запускать как логит, так и пробит.

Синтаксис такой же, как в lfe::felm()

См. этот пример, взятый из документации:

      # Generate an artificial data set for logit models
library(alpaca)
data <- simGLM(1000L, 20L, 1805L, model = "logit")
# Fit 'feglm()'
mod <- feglm(y ~ x1 + x2 + x3 | i + t, data)
summary(mod)