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

Posted 2023-02-16

【中文标题】R中的大型固定效应二项式回归

【英文标题】：Large fixed effects binomial regression in R

【发布时间】：2015-04-22 10:20:12

【问题描述】：

我需要对一个相对较大的数据框进行逻辑回归，该数据框包含 480.000 个条目和 3 个固定效应变量。固定效应 var A 有 3233 级，var 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 package 的 felm 回归函数，在运行回归之前减去固定效应。这完美地工作，并允许我在短短几分钟内将上述回归作为正常的线性模型运行。但是，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 函数使用“Gauss-Seidel”算法执行与lfe 包类似的固定效应减法，用于逻辑回归。不幸的是，该软件包不再被开发。对 R 来说相对较新，并且对统计没有深入的经验，我无法理解输出，也不知道如何以一种可以给我正常的“效果大小”、“模型拟合”、“ glm 回归摘要提供的显着性区间”指标。

我向包的作者发送了一条消息。他们友好地回应如下：

该包不提供与 glm 对象相同格式的输出。然而，你 可以很容易地计算出大部分拟合统计量（ 估计，拟合优度）给定当前输出（在 CRAN 中 版本，我相信当前输出是估计的向量 系数，以及相关的标准误差向量；同样的 协方差分量，但如果你不需要担心它们 是没有随机效应的拟合模型）。只需要注意 用于计算标准误差的协方差矩阵是 与相关的精度矩阵的对角块的逆 Gauss-Seidel 算法，因此他们倾向于低估 联合可能性的标准误。我不维护 包裹不再，我没有时间进入具体 细节;包背后的开创性理论可以在 paper referenced in the manual，其他都需要解决 由你用笔和纸 :)。

如果有人可以解释如何“轻松计算大多数拟合统计数据”，让没有任何统计学知识的人可以理解（可能是不可能的），或者提供 R 代码来举例说明如何做到这一点我将非常感激！

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

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

我很乐意就上述任何问题或完全不同的方法在 R 中运行具有内存限制的多个大型固定效应的逻辑回归提供建议。

【问题讨论】：

是来自Matrix 包的?sparse.model.matrix，还是来自MatrixModels 包的model.Matrix，有用吗？

【参考方案1】：

我同意任何人（我猜是@Ben Bolker？）建议您使用MatrixModels 中的glm4 函数。首先，如果您使用 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 类似乎没有摘要方法。虽然，插槽似乎有明显的命名，并且不应该花很长时间来获得例如估计器的标准误差、计算方差估计等。

【讨论】：

这是我经历过的最快的一个

【参考方案2】：

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

【参考方案3】：

看看

glmmbootglmmML

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)

计算时间差“巨大”！

【讨论】：

这看起来很棒。唯一剩下的问题：您知道如何将其用于多个集群吗？ IE。如果我有 3 个不同的固定效果，那将是 4 个“组”集群。如何在公式中定义这些单独的集群？编辑：仅仅通过使用 paste(A,B) 或 factor(A):factor(B) 从不同的固定效果创建“新”因子不起作用/不会产生与正常 glm 相同的结果。

@Phil 感谢您提出原始问题 - 似乎与我现在正在经历的过程相同。您最终发现如何定义多个集群了吗？