R中的大固定效应二项回归



我需要在一个相对较大的数据框架上运行逻辑回归,该数据框架包含480.000个条目和3个固定的影响变量。固定效果变量A有3233个级别,变量B有2326个级别,var C有811个级别。所以我总共有6370个固定效果。数据是横截面的。如果我不能使用正常的glm函数运行此回归,因为回归矩阵对我的内存来说太大了(我得到消息"Error: cannot allocate vector of size 22.9 Gb")。我正在寻找在我的Macbook Air(OS X 10.9.5 8GB RAM)上运行此回归的替代方法。我还可以访问一台16GB RAM的服务器。

我尝试过用几种不同的方法来解决这个问题,但到目前为止,没有一种能带来令人满意的结果:

lfe/felm:使用lfe包的felm回归函数,该函数在运行回归之前减去固定效应。这非常有效,使我能够在几分钟内将上述回归作为一个正常的线性模型运行。然而,lfe不支持逻辑回归和glms。因此,felm很好地了解了模型是否适合不同的模型,但对最终的逻辑回归模型不起作用。

biglm/biglm:我考虑过使用bigglm将我的函数分解为更易于管理的块。然而,一些来源(例如链接1、链接2、链接3)提到,为了实现这一点,因子水平需要在块之间保持一致,即每个块必须包含每个因子变量的每个因子中的至少一个。因子A和因子B包含只出现一次的级别,所以我不能将集合拆分为级别一致的不同块。如果我删除10个固定影响因素A和8个因素B(一个小变化),我将只剩下4个以上级别的因素,将我的数据划分为4个块将使其更易于管理。然而,我仍然需要弄清楚如何对我的df进行排序,以确保我的480.000个条目被排序为4个块,其中3个因子中的每个因子的每个因子级别至少出现一次。

GlmmGS/glmgs:具有相同名称的包中的glmmgs函数使用"高斯-塞德尔"算法执行固定效果减法,类似于逻辑回归的lfe包。不幸的是,该方案已不再开发。由于对R相对较新,对统计学没有深入的经验,我无法理解输出,也不知道如何将其转换为glm回归摘要提供的正常"效应大小"、"模型拟合"、"显著性区间"指标。

我给包的作者发了一条信息。他们友好地回应如下:

该包不提供与glm对象相同格式的输出。然而,你可以很容易地计算大多数拟合统计数据(估计,拟合优度)版本,我认为当前输出是系数和相关的标准误差矢量;相同协方差分量,但如果您是没有随机效应的拟合模型)。只需注意用于计算标准误差的协方差矩阵是与相关联的精度矩阵的对角块的逆高斯-塞德尔算法,因此他们倾向于低估联合似然的标准误差。我没有维护包裹再长了,我没有时间进入具体细节;该方案背后的开创性理论可以在手册中提到的论文,其他的都需要解决由您用笔和纸:)。

如果有人能解释如何"轻松计算大多数拟合统计数据",让没有受过统计学教育的人能够理解(可能是不可能的),或者提供R代码来展示如何做到这一点,我将不胜感激!

革命分析:我在一台模拟Mac上Windows7的虚拟机上安装了革命分析企业。该程序有一个名为RxLogit的函数,该函数针对大型逻辑回归进行了优化。使用RxLogit函数,我得到了the error (Failed to allocate 326554568 bytes. Error in rxCall("RxLogit", params) : bad allocation),因此该函数似乎也遇到了内存问题。然而,该软件使我能够在分布式计算集群上运行回归。因此,我可以通过在具有大量内存的集群上购买计算时间来"解决问题"。然而,我想知道革命分析程序是否提供了任何我不知道的公式或方法,可以让我进行某种类似lfe的固定效果减法运算或类似bigglm的分块运算,将因素考虑在内。

MatrixModels/glm4:有人建议我使用具有sparse = TRUE属性的MatrixModels包的glm4函数来加快计算速度。如果我在所有固定效应的情况下运行glm4回归,我会得到"Error in Cholesky(crossprod(from), LDL = FALSE) : internal_chm_factor: Cholesky factorization failed"错误。如果我只在固定效应变量B或a和C的情况下执行它,计算会起作用并返回"glpModel"对象。与glmmGS一样,我在将输出转换为对我来说有意义的形式时遇到了一些问题,因为标准的summary()方法似乎不起作用。

我很乐意就上述任何问题或完全不同的方法获得建议,以在记忆受限的R中运行具有多个大固定效应的逻辑回归。

查看

glmmboot{glmmML}

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

还有一份Brostrom和Holmberg的漂亮文件(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)

计算时间差是"巨大的"!

对于子孙后代,我还想推荐speedglm包,我发现它在尝试对大型数据集执行逻辑回归时很有用。它似乎占用了大约一半的内存,并且完成得比glm快得多。

我同意任何人(我想是@Ben Bolker?)建议您使用MatrixModels中的glm4函数。首先,如果使用sparse参数,它可以解决内存问题。一个包含480.000个条目和6370个固定效果的密集设计矩阵将需要6371*480.000*8=24.64640.000个字节。然而,你的设计矩阵将是非常稀疏的(很多零),所以如果你使用稀疏的设计矩阵,你可以使用更小的(内存中的)设计矩阵。其次,您可以利用稀疏性来更快地进行估计。

至于选项,快速搜索显示speedglm也有sparse参数,尽管我还没有尝试过。无论你使用什么方法,我的关键是它都应该使用你的设计矩阵是稀疏的,以减少计算时间和内存需求。

您得到的错误(Error in Cholesky(crossprod(from), LDL = FALSE) : internal_chm_factor: Cholesky factorization failed" error)可能是因为您的设计矩阵是奇异的。在这种情况下,您的问题没有唯一的解决方案,一些选项是合并一些组级别,使用惩罚或随机效应模型。

您说得对,glpModel类似乎没有一个summary方法。尽管如此,这些槽似乎有明显的名称,你应该不会花很长时间来获得,例如,你的估计器上的标准误差,计算方差估计等。

您可能需要羊驼软件包中的函数feglm。它允许同时运行logit和probit。

语法与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)

最新更新