在Python中,当依赖变量存在较大值时,线性回归会失败。

4
我正在尝试使用Python(使用pandas.stats.api.ols)重写一个预测模型(在Stata中),但是在线性回归方面遇到了问题:由pandas计算的系数和截距与Stata中的不匹配。
调查显示,根本原因可能是因变量的值非常大。我基于以下发现怀疑这一点:
1)我在Python中创建了一个简单的DataFrame,并对其进行了线性回归:
from pandas.stats.api import ols
import pandas as pd
df = pd.DataFrame({"A": [10,20,30,40,50,21], "B": [20, 30, 10, 40, 50,98], "C": [32, 234, 23, 23, 31,21], "D":[12,28,12,98,51,87], "E": [1,8,12,9,12,91]})
ols(y=df['A'], x=df[['B','C', 'D', 'E']])

LR的总结如下:

-------------------------Summary of Regression Analysis-------------------------

Formula: Y ~ <B> + <C> + <D> + <E> + <intercept>

Number of Observations:         6
Number of Degrees of Freedom:   5

R-squared:         0.4627
Adj R-squared:    -1.6865

Rmse:             23.9493

F-stat (4, 1):     0.2153, p-value:     0.9026

Degrees of Freedom: model 4, resid 1

-----------------------Summary of Estimated Coefficients------------------------
      Variable       Coef    Std Err     t-stat    p-value    CI 2.5%   CI 97.5%
--------------------------------------------------------------------------------
             B     0.3212     1.1176       0.29     0.8218    -1.8693     2.5117
             C    -0.0488     0.1361      -0.36     0.7806    -0.3155     0.2178
             D     0.1512     0.4893       0.31     0.8092    -0.8077     1.1101
             E    -0.4508     0.8268      -0.55     0.6822    -2.0713     1.1697
     intercept    20.9222    23.6280       0.89     0.5386   -25.3887    67.2331
---------------------------------End of Summary---------------------------------

我将这个DataFrame保存为Stata .dta文件,并在Stata中运行了LR,如下所示:
 use "/tmp/lr.dta", clear
 reg A B C D E

结果是相同的:
      Source |       SS       df       MS              Number of obs =       6
-------------+------------------------------           F(  4,     1) =    0.22
       Model |  493.929019     4  123.482255           Prob > F      =  0.9026
    Residual |  573.570981     1  573.570981           R-squared     =  0.4627
-------------+------------------------------           Adj R-squared = -1.6865
       Total |      1067.5     5       213.5           Root MSE      =  23.949

------------------------------------------------------------------------------
           A |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           B |   .3211939   1.117591     0.29   0.822    -13.87914    14.52153
           C |  -.0488429   .1360552    -0.36   0.781    -1.777589    1.679903
           D |   .1512067   .4892539     0.31   0.809    -6.065353    6.367766
           E |  -.4508122   .8267897    -0.55   0.682    -10.95617    10.05455
       _cons |    20.9222   23.62799     0.89   0.539    -279.2998    321.1442
------------------------------------------------------------------------------

我在R中尝试了这个方法,结果相同。

2) 但是,如果我在Python中增加了因变量的值:

df = pd.DataFrame({"A": [10.0,20.0,30.0,40.0,50.0,21.0]})
df['B'] = pow(df['A'], 30)
df['C'] = pow(df['A'], 5)
df['D'] = pow(df['A'], 15)
df['E'] = pow(df['A'], 25)

我已经确保这里的所有列都使用float64: df.dtypes A float64 B float64 C float64 D float64 E float64 dtype: object
我得到的结果是:
-------------------------Summary of Regression Analysis-------------------------

Formula: Y ~ <B> + <C> + <D> + <E> + <intercept>

Number of Observations:         6
Number of Degrees of Freedom:   2

R-squared:        -0.7223
Adj R-squared:    -1.1528

Rmse:             21.4390

F-stat (4, 4):    -1.6775, p-value:     1.0000

Degrees of Freedom: model 1, resid 4

-----------------------Summary of Estimated Coefficients------------------------
      Variable       Coef    Std Err     t-stat    p-value    CI 2.5%   CI 97.5%
--------------------------------------------------------------------------------
             B    -0.0000     0.0000      -0.00     0.9973    -0.0000     0.0000
             C     0.0000     0.0000       0.00     1.0000    -0.0000     0.0000
             D     0.0000     0.0000       0.00     1.0000    -0.0000     0.0000
             E     0.0000     0.0000       0.00     0.9975    -0.0000     0.0000
     intercept     0.0000    21.7485       0.00     1.0000   -42.6271    42.6271
---------------------------------End of Summary---------------------------------

但是在Stata中,我得到了非常不同的结果:
      Source |       SS       df       MS              Number of obs =       6
-------------+------------------------------           F(  4,     1) =  237.35
       Model |  1066.37679     4  266.594196           Prob > F      =  0.0486
    Residual |   1.1232144     1   1.1232144           R-squared     =  0.9989
-------------+------------------------------           Adj R-squared =  0.9947
       Total |      1067.5     5       213.5           Root MSE      =  1.0598

------------------------------------------------------------------------------
           A |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           B |  -1.45e-45   2.32e-46    -6.24   0.101    -4.40e-45    1.50e-45
           C |   2.94e-06   3.67e-07     8.01   0.079    -1.72e-06    7.61e-06
           D |  -3.86e-21   6.11e-22    -6.31   0.100    -1.16e-20    3.90e-21
           E |   4.92e-37   7.88e-38     6.24   0.101    -5.09e-37    1.49e-36
       _cons |   9.881129    1.07512     9.19   0.069    -3.779564    23.54182
------------------------------------------------------------------------------

在R中的结果与Stata相一致: lm(formula = A ~ B + C + D + E, data = stata)

Residuals:
         1          2          3          4          5          6 
-1.757e-01  8.211e-01  1.287e-03 -1.269e-06  1.289e-09 -6.467e-01 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  9.881e+00  1.075e+00   9.191    0.069 .
B           -1.449e-45  2.322e-46  -6.238    0.101  
C            2.945e-06  3.674e-07   8.015    0.079 .
D           -3.855e-21  6.106e-22  -6.313    0.100  
E            4.919e-37  7.879e-38   6.243    0.101  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 1.06 on 1 degrees of freedom
Multiple R-squared:  0.9989,  Adjusted R-squared:  0.9947 
F-statistic: 237.3 on 4 and 1 DF,  p-value: 0.04864

因此,据我所知,在这里 pandas 存在一些问题。请问有人能够提供帮助和建议吗?

仔细看一下这个代码:df['B'] = pow(df['A'], 30),它将 B 替换为 A**30。类似的代码行还有 CDE,分别用 A 的不同幂次替换了每个变量的原始数据。这是一个打字错误吗?你是否想要新的 B 等于旧的 B **30?你在 Stata 和 R 中使用了什么?它们匹配吗? - Paul
@Paul 不,这不是打字错误。正如我所提到的,我在Python中将数据保存为 .dta 文件,然后在R和Stata中加载它们。所以数据是匹配的。 - Jackie Wong
我改用了sklearn.linear_model.LinearRegression,因为它支持归一化处理,现在的结果与R和Stata相匹配。 - Jackie Wong
1个回答

2

我认为这是由于Python中的相对精度问题(不仅限于Python,大多数其他编程语言也是如此,例如C++)。np.finfo(float).eps给出了2.2204460492503131e-16,因此当你尝试进行任何基本操作,比如+ - * /时,小于eps*max_value_of_your_data的所有内容都将被视为实际上是0。例如,1e117 + 1e100 == 1e117返回True,因为1e100/1e117 = 1e-17 < eps。现在看看你的数据。

# your data
# =======================
print(df)

    A           B          C           D           E
0  10  1.0000e+30     100000  1.0000e+15  1.0000e+25
1  20  1.0737e+39    3200000  3.2768e+19  3.3554e+32
2  30  2.0589e+44   24300000  1.4349e+22  8.4729e+36
3  40  1.1529e+48  102400000  1.0737e+24  1.1259e+40
4  50  9.3132e+50  312500000  3.0518e+25  2.9802e+42
5  21  4.6407e+39    4084101  6.8122e+19  1.1363e+33

当考虑相对精度时,

# ===================================================
import numpy as np

np.finfo(float).eps # 2.2204460492503131e-16

df[df < df.max().max()*np.finfo(float).eps] = 0
df

   A           B  C  D           E
0  0  0.0000e+00  0  0  0.0000e+00
1  0  1.0737e+39  0  0  0.0000e+00
2  0  2.0589e+44  0  0  8.4729e+36
3  0  1.1529e+48  0  0  1.1259e+40
4  0  9.3132e+50  0  0  2.9802e+42
5  0  4.6407e+39  0  0  0.0000e+00

所以y(A)没有任何变化,这就是为什么statsmodels返回所有0系数的原因。提醒一下,在运行回归之前,始终最好将数据标准化。


1
回复:“在运行回归之前,规范化数据通常是一个好的做法。” http://stats.stackexchange.com/questions/29781/when-should-you-center-your-data-when-should-you-standardize - Brandon Bertelsen
@jianxun-li:非常感谢!这正是我怀疑的。R或Stata是否对浮点数有特殊处理,还是在运行LR之前隐式地进行归一化? - Jackie Wong
@brandon-bertelsen:感谢提供链接。我会尝试理解它。很抱歉我不是一个统计学家,对此了解甚少。 - Jackie Wong
@JackieWong 我不确定那些领域特定的统计软件如何处理这种相对精度问题。但正如你所提到的,首先归一化RHS变量,然后在回归估计后通过标准误差进行缩放似乎是一种有效的方法。 - Jianxun Li

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接