我想使用 Statsmodels 在 Python 中进行逻辑回归。
X 和 y 各有 750 行，y 是二进制结果，X 中是 10 个特征(包括拦截)。</p>
这是 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: overflow encountered in exp
return 1/(1+np.exp(-X))/home/ipattern/anaconda3/lib/python3.6/site-packages/statsmodels/discrete/discrete_model.py:1264: RuntimeWarning: divide by zero encountered in log
return 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) 某处。
还有伪R-squ。是-inf 和Log-Likelihood 是-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 值很大时，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 的近似值。由于函数以 double 实现的方式，这可能只是一个数值问题。
关于python - Statsmodels 抛出 "overflow in exp"和 "divide by zero in log"警告，伪 R 平方是 -inf，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/46173061/