使用统计包对数百个协变量进行最大似然估计的方法答案

【问题标题】：Approaches for using statistics packages for maximum likelihood estimation for hundreds of covariates使用统计包对数百个协变量进行最大似然估计的方法
【发布时间】：2019-07-19 10:47:50
【问题描述】：

我正在尝试研究最大似然估计的分布，特别是针对大量协变量 p 和高维状态（这意味着 p/n，样本大小为 n，约为 1/5）。我正在生成数据，然后使用statsmodels.api.Logit 将参数拟合到我的模型。

问题是，这似乎只适用于低维状态（如 300 个协变量和 40000 个观察值）。具体来说，我知道已经达到最大迭代次数，对数似然是 inf，即已经发散，并且出现“奇异矩阵”错误。

我不知道如何解决这个问题。最初，当我仍在使用较小的值（比如 80 个协变量，4000 个观察值）时，我偶尔会遇到这个错误，我设置了最多 70 次迭代而不是 35 次。这似乎有帮助。

但是现在它显然无济于事，因为我的对数似然函数正在发散。这不仅仅是在最大迭代次数内不收敛的问题。

很容易回答，这些包根本不适合处理这些数字，但是已经有论文专门研究了这种高维机制，比如here，其中使用了p=800 协变量和n=4000 观察值。

当然，本文使用的是 R 而不是 python。不幸的是我不知道 R。但是我应该认为 python 优化应该具有相当的“质量”？

我的问题：

R 是否比 python statsmodels 更适合处理这种高 p/n 机制中的数据？如果是这样，为什么以及可以使用 R 的技术来修改 python statsmodels 代码？

如何修改我的代码以适用于 p=800 和 n=4000 左右的数字？

【问题讨论】：

这个问题需要code和data（或者意味着生成一些）。也就是说，我记得几天前你的一个类似问题。你的主要目标是什么？这是对通过逻辑回归获得的 MLE 的渐近行为的好奇，还是您对这么多数据有特定的研究问题，您认为逻辑回归可能是解决问题的方法吗？
@N.Wouda 这确实只是对逻辑回归渐近行为的研究。我问过类似的问题，但没有一个是相同的。以前，我什至不确定哪些包最适合用于最大似然估计。你的帖子对此很有帮助。但现在我意识到 statsmodels 包似乎无法在我感兴趣的机制中获得 MLE（大量协变量，少量观察）。然而，我意识到不同的软件包的工作方式非常不同。有些内置了正则化，这在中并不好
我的书房。 statsmodels 似乎没有这个，但它无法获得大维数据的估计。我曾尝试更改逻辑回归模块中的代码，但无济于事。然而，有些论文已经在这个制度下运行了模拟，所以必须用一些包来探索......除非他们只是尝试运行它很多次并且幸运地在
目前，我只是想弄清楚如何让计算工具发挥作用（我的弱项），以便我可以将理论（我的强项）与模拟进行比较。

标签： python optimization statistics statsmodels singular

【解决方案1】：

在您当前使用的代码中（来自其他几个问题），您隐式使用了 Newton-Raphson 方法。这是sm.Logit 模型的默认值。它计算并反转 Hessian 矩阵以加快估计速度，但这对于大型矩阵来说是难以置信的 - 更不用说当矩阵接近奇异时经常导致数值不稳定，正如您已经看到的那样。这在相关的维基百科entry上有简要说明。

您可以通过使用不同的求解器来解决此问题，例如bfgs（或lbfgs），像这样，

model = sm.Logit(y, X)
result = model.fit(method='bfgs')

即使使用n = 10000, p = 2000，这对我来说也运行得很好。

除了估计之外，更成问题的是，您用于生成样本的代码会导致数据受到很大程度的准可分离性影响，在这种情况下，整个 MLE 方法充其量是有问题的。您应该紧急调查这一点，因为它表明您的数据可能不像您希望的那样表现良好。准可分性很好解释here。

【讨论】：

太棒了！ bfgs 真的很好用，我真的很想看到它给我的估计器带来了更大的差异——就像我读过的论文一样。我会试着弄清楚为什么会这样:)