Statsmodels 抛出“exp 溢出”和“log 中除以零”警告，伪 R 平方为 -inf答案

【问题标题】：Statsmodels throws "overflow in exp" and "divide by zero in log" warnings and pseudo-R squared is -infStatsmodels 抛出“exp 溢出”和“log 中除以零”警告，伪 R 平方为 -inf
【发布时间】：2018-02-20 17:16:59
【问题描述】：

我想在 Python 中使用 Statsmodels 进行逻辑回归。

X 和 y 各有 750 行，y 是二元结果，X 中是 10 个特征（包括 intecept）。

这里是 X 的前 12 行（最后一列是截距）：

      lngdp_      lnpop    sxp      sxp2    gy1    frac  etdo4590  geogia  \
0   7.367709  16.293980  0.190  0.036100 -1.682   132.0         1   0.916   
1   7.509883  16.436258  0.193  0.037249  2.843   132.0         1   0.916   
2   7.759187  16.589224  0.269  0.072361  4.986   132.0         1   0.916   
3   7.922261  16.742384  0.368  0.135424  3.261   132.0         1   0.916   
4   8.002359  16.901037  0.170  0.028900  1.602   132.0         1   0.916   
5   7.929126  17.034786  0.179  0.032041 -1.465   132.0         1   0.916   
6   6.594413  15.627563  0.360  0.129600 -9.321  4134.0         0   0.648   
7   6.448889  16.037861  0.476  0.226576 -2.356  3822.0         0   0.648   
8   8.520786  16.919334  0.048  0.002304  2.349   434.0         1   0.858   
9   8.637107  16.991980  0.050  0.002500  2.326   434.0         1   0.858   
10  8.708144  17.075489  0.042  0.001764  1.421   465.0         1   0.858   
11  8.780480  17.151779  0.080  0.006400  1.447   496.0         1   0.858   

    peace  intercept  
0    24.0        1.0  
1    84.0        1.0  
2   144.0        1.0  
3   204.0        1.0  
4   264.0        1.0  
5   324.0        1.0  
6     1.0        1.0  
7    16.0        1.0  
8   112.0        1.0  
9   172.0        1.0  
10  232.0        1.0  
11  292.0        1.0

这是我的代码：

import statsmodels.api as sm

logit = sm.Logit(y, X, missing='drop')
result = logit.fit()
print(result.summary())

这是输出：

     Optimization terminated successfully.

     Current function value: inf

     Iterations 9
/home/ipattern/anaconda3/lib/python3.6/site-packages/statsmodels/discrete/discrete_model.py:1214： RuntimeWarning：exp 中遇到溢出
返回 1/(1+np.exp(-X))

/home/ipattern/anaconda3/lib/python3.6/site-packages/statsmodels/discrete/discrete_model.py:1264： RuntimeWarning：在日志中遇到除以零
返回 np.sum(np.log(self.cdf(q*np.dot(X,params))))

                           Logit Regression Results                           
==============================================================================
Dep. Variable:                  warsa   No. Observations:                  750
Model:                          Logit   Df Residuals:                      740
Method:                           MLE   Df Model:                            9
Date:                Tue, 12 Sep 2017   Pseudo R-squ.:                    -inf
Time:                        11:16:58   Log-Likelihood:                   -inf
converged:                       True   LL-Null:                   -4.6237e+05
                                        LLR p-value:                     1.000
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
lngdp_        -0.9504      0.245     -3.872      0.000      -1.431      -0.469
lnpop          0.5105      0.128      3.975      0.000       0.259       0.762
sxp           16.7734      5.206      3.222      0.001       6.569      26.978
sxp2         -23.8004     10.040     -2.371      0.018     -43.478      -4.123
gy1           -0.0980      0.041     -2.362      0.018      -0.179      -0.017
frac          -0.0002    9.2e-05     -2.695      0.007      -0.000   -6.76e-05
etdo4590       0.4801      0.328      1.463      0.144      -0.163       1.124
geogia        -0.9919      0.909     -1.091      0.275      -2.774       0.790
peace         -0.0038      0.001     -3.808      0.000      -0.006      -0.002
intercept     -3.4375      2.486     -1.383      0.167      -8.310       1.435
==============================================================================

底部的系数、std err、p 值等都是正确的（我知道这一点是因为我有“解决方案”）。

但是你可以看到Current function value is inf 我认为这是错误的。

我收到两个警告。显然 statsmodels 确实是 np.exp(BIGNUMBER)，例如np.exp(999) 和 np.log(0) 某处。

还有Pseudo R-squ. is -inf 和Log-Likelihood is -inf，我认为不应该是-inf。

那我做错了什么？

编辑：

X.describe():

           lngdp_       lnpop         sxp        sxp2         gy1  \
count  750.000000  750.000000  750.000000  750.000000  750.000000   
mean     7.766948   15.702191    0.155329    0.043837    1.529772   
std      1.045121    1.645154    0.140486    0.082838    3.546621   
min      5.402678   11.900227    0.002000    0.000004  -13.088000   
25%      6.882694   14.723123    0.056000    0.003136   -0.411250   
50%      7.696212   15.680984    0.111000    0.012321    1.801000   
75%      8.669355   16.652981    0.203000    0.041209    3.625750   
max      9.851826   20.908354    0.935000    0.874225   14.409000   

              frac    etdo4590      geogia       peace  intercept  
count   750.000000  750.000000  750.000000  750.000000      750.0  
mean   1812.777333    0.437333    0.600263  348.209333        1.0  
std    1982.106029    0.496388    0.209362  160.941996        0.0  
min      12.000000    0.000000    0.000000    1.000000        1.0  
25%     176.000000    0.000000    0.489250  232.000000        1.0  
50%     864.000000    0.000000    0.608000  352.000000        1.0  
75%    3375.000000    1.000000    0.763000  472.000000        1.0  
max    6975.000000    1.000000    0.971000  592.000000        1.0

logit.loglikeobs(result.params)：

array([ -4.61803704e+01,  -2.26983454e+02,  -2.66741244e+02,
        -2.60206733e+02,  -4.75585266e+02,  -1.76454554e+00,
        -4.86048292e-01,  -8.02300533e-01,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,  -6.02780923e+02,
        -4.12209348e+02,  -6.42901288e+02,  -6.94331125e+02,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf,
                   -inf,             -inf,             -inf, ...

(logit.exog * np.array(result.params)).min(0):

array([ -9.36347474,   6.07506083,   0.03354677, -20.80694575,
        -1.41162588,  -1.72895247,   0.        ,  -0.9631801 ,
        -2.23188846,  -3.4374963 ])

数据集：

X：https://pastebin.com/VRNSepBg

y：https://pastebin.com/j2Udyc7m

【问题讨论】：

@Phillip 如果数据集是公开的，那么您可以将案例发布到 statsmodels 问题跟踪器 github.com/statsmodels/statsmodels/issues 。我猜 Logit 需要使预测的概率远离 0 和 1，而这目前还没有完成。但是，我怀疑在这样的极端情况下，对数似然 llf 将由选择作为削波阈值的任何内容支配。
好的，我会记住的，谢谢。

标签： python statistics logistic-regression statsmodels

【解决方案1】：

我很惊讶它在这种情况下仍然收敛。

当 x 值很大时，使用 Logit 或 Poisson 中使用的 exp 函数可能会出现收敛问题。这通常可以通过重新调整回归量来避免。

但是，在这种情况下，我的猜测是 x 中的异常值。第 6 列的值类似于 4134.0，而其他列的值要小得多。

您可以检查每个观测值 logit.loglikeobs(result.params) 的对数似然，以查看哪些观测值可能会导致问题，其中 logit 是引用模型的名称

每个预测变量的贡献也可能有所帮助，例如

np.argmax(np.abs(logit.exog * result.params), 0)

或

(logit.exog * result.params).min(0)

如果只是一个或几个观察结果，那么删除它们可能会有所帮助。重新调整 exog 很可能对此无济于事，因为在收敛时，它只会通过重新调整估计系数来补偿。

还要检查是否存在编码错误或较大的值作为缺失值的占位符。

编辑

鉴于 loglikeobs 中 -inf 的数量似乎很大，我认为可能存在比异常值更根本的问题，因为 Logit 模型不是该数据集正确指定的最大似然模型。

一般来说有两种可能（因为我还没有看到数据集）：

完美分离：Logit 假设预测的概率远离零和一。在某些情况下，解释变量或它们的组合可以完美地预测因变量。在这种情况下，参数要么没有被识别，要么变为正或负无穷大。实际参数估计取决于优化的收敛标准。 Statsmodels Logit 检测到一些情况，然后引发 PerfectSeparation 异常，但它没有检测到所有部分分离的情况。

Logit 或 GLM-Binomial 属于单参数线性指数族。在这种情况下，参数估计仅取决于指定的均值函数和隐含方差。它不需要正确指定似然函数。因此，即使给定数据集的似然函数不正确，也可以获得良好（一致）的估计。在这种情况下，解是一个准最大似然估计，但对数似然值是无效的。

这可能会产生这样的效果，即收敛性和数值稳定性方面的结果取决于如何处理边缘或极端情况的计算细节。在某些情况下，Statsmodels 正在裁剪这些值以使它们远离界限，但还不是无处不在。

困难在于弄清楚如何处理数值问题，并避免在底层模型不适合或与数据不兼容时在不警告用户的情况下返回“某些”数字。

也许llf = -inf 在这种情况下是“正确”的答案，任何有限数都只是-inf 的近似值。也许这只是一个数值问题，因为函数是以双精度实现的方式。

【讨论】：

感谢您的帮助！我编辑了我的帖子并添加了您所说的内容。在带有 glm() 函数的 R 中，一切看起来都很正常。没有警告等。例如，伪 R 平方不是 -inf。
您也可以尝试 statmodels GLM Binomial，它与 Logit 的二元 endog 模型相同，但使用不同的数值方法。在非边缘和非极端情况下，这两个模型产生相同的结果。 GLM 也存在偏差（但可能会遇到相同的溢出问题）。
我在问题的底部添加了数据集（如果您想尝试一下）。我尝试了 statsmodels GLM Binomial，它工作正常！正如你所说，很难找到 statsmodels 的 Logit 函数有什么问题，所以我认为使用 GLM Binomial 函数。另外，llf = -inf 是正确答案是什么意思？