python - Python(和 R)和 Stata 中的线性回归之间的区别

Question

我正在将 Stata 模型移植到 Python，使用相同的输入数据(可用 @ https://drive.google.com/file/d/0B8PLy9yAUHvlcTI1SG5sdzdnaWc/view?usp=sharing)看到 Python 和 Stata 的线性回归结果不同

Stata代码如下:

reg growth time*
predict ghat
predict resid, residuals

结果是(前 5 行):

. list growth ghat resid

     +----------------------------------+
     |    growth       ghat       resid |
     |----------------------------------|
  1. | 2.3527029   2.252279    .1004239 |
  2. |  2.377728   2.214551     .163177 |
  3. | 2.3547957   2.177441     .177355 |
  4. | 3.0027488   2.140942    .8618064 |
  5. | 3.0249328    2.10505    .9198825 |

在 Python 中，代码是:

import pandas as pd
from sklearn.linear_model import LinearRegression


def linear_regression(df, dep_col, indep_cols):
  lf = LinearRegression(normalize=True)
  lf.fit(df[indep_cols.split(' ')], df[dep_col])

  return lf

df = pd.read_stata('/tmp/python.dta')
lr = linear_regression(df, 'growth', 'time time2 time3 time4 time5')

df['ghat'] = lr.predict(df['time time2 time3 time4 time5'.split(' ')])
df['resid'] = df.growth - df.ghat

df.head(5)['growth ghat resid'.split(' ')]

结果是:

     growth      ghat     resid
0  2.352703  3.026936 -0.674233
1  2.377728  2.928860 -0.551132
2  2.354796  2.833610 -0.478815
3  3.002749  2.741135  0.261614
4  3.024933  2.651381  0.373551

我也在 R 中尝试过，并得到了与在 Python 中相同的结果。我无法找出根本原因:是因为 Stata 中使用的算法有点不同吗？我从源代码可以看出sklearn使用的是普通最小二乘法，但不知道Stata中的那个。

有人可以在这里提供建议吗？

------------ 编辑 1 ------------

我尝试将 Stata 中的数据类型指定为 double，但 Stata 仍然产生与使用 float 相同的结果。生成的Stata代码如下:

gen double growth = .
foreach lag in `lags' {
    replace growth = ma_${metric}_per_`group' / l`lag'.ma_${metric}_per_`group' - 1 if nlag == `lag' & in_sample
}

gen double time = day - td(01jan2010) + 1
forvalues i = 2/5 {
    gen double time`i' = time^`i'
}

------------ 编辑 2 ------------

已确认 Stata 确实由于共线性而删除了 time 变量。由于我们的 Stata 代码启用 quiet 模型来抑制不需要的消息，因此之前没有看到该消息。根据我的调查，这不能在 Stata 中禁用。因此看来我还需要检测共线性并删除 Python 中的共线列。

. reg growth time*,
note: time omitted because of collinearity

      Source |       SS       df       MS              Number of obs =     381
-------------+------------------------------           F(  4,   376) =  126.10
       Model |  37.6005042     4  9.40012605           Prob > F      =  0.0000
    Residual |  28.0291465   376  .074545602           R-squared     =  0.5729
-------------+------------------------------           Adj R-squared =  0.5684
       Total |  65.6296507   380  .172709607           Root MSE      =  .27303

------------------------------------------------------------------------------
      growth |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        time |          0  (omitted)
       time2 |  -.0098885   .0009231   -10.71   0.000    -.0117037   -.0080734
       time3 |   .0000108   1.02e-06    10.59   0.000     8.77e-06    .0000128
       time4 |  -4.40e-09   4.20e-10   -10.47   0.000    -5.22e-09   -3.57e-09
       time5 |   6.37e-13   6.15e-14    10.35   0.000     5.16e-13    7.58e-13
       _cons |   3322.727   302.7027    10.98   0.000     2727.525     3917.93
------------------------------------------------------------------------------

Answer 1

预测变量是 time 的 1 次... 5 次方，它在 1627 年和 2007 年之间变化(大概是一个日历年，这无关紧要)。即使使用现代软件，改变时间原点以减少数值应变也是明智的，例如使用 (time - 1800) 的幂。

无论如何，重做回归表明 Stata 将第一个预测变量丢弃为共线。 Python 和 R 会发生什么？这些是对数字棘手挑战的不同 react 。

(拟合五次多项式很少有科学值(value)，但这在这里可能无关紧要。基于 2 到 5 次方的拟合曲线对于这些数据效果不是很好，这些数据看起来很经济。更有意义的是前 5 个残差都是正数，但并非所有残差都是正数!)