Python（和R）和Stata中的线性回归之间的区别

Question

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

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 会发生什么？这些是对数字棘手挑战的不同反应。

（拟合五次多项式很少有科学价值，但这在这里可能无关紧要。基于 2 到 5 次方的拟合曲线不适用于这些看起来经济的数据。更有意义的是前 5 个残差都是正数，但并非全部都是这样！）

Answer 2

这是一个通配符问题。在您的 Stata 代码中，time* 将匹配 time2, time3... 但不匹配 time。如果将Python 代码更改为lr = linear_regression(df, 'growth', 'time2 time3 time4 time5')，它将得到完全相同的结果。

编辑

出现Stata 删除了第一个自变量。拟合可以如下可视化：

lr1 = linear_regression(df, 'growth', 'time time2 time3 time4 time5')
lr2 = linear_regression(df, 'growth', 'time2 time3 time4 time5')
pred_x1 = ((np.linspace(1620, 2000)[..., np.newaxis]**np.array([1,2,3,4,5]))*lr1.coef_).sum(1)+lr1.intercept_
pred_x2 = ((np.linspace(1620, 2000)[..., np.newaxis]**np.array([2,3,4,5]))*lr2.coef_).sum(1)+lr2.intercept_
plt.plot(np.linspace(1620, 2000), pred_x1, label='Python/R fit')
plt.plot(np.linspace(1620, 2000), pred_x2, label='Stata fit')
plt.plot(df.time, df.growth, '+', label='Data')
plt.legend(loc=0)

还有残差平方和：

In [149]:

pred1 = (df.time.values[..., np.newaxis]**np.array([1,2,3,4,5])*lr1.coef_).sum(1)+lr1.intercept_
pred2 = (df.time.values[..., np.newaxis]**np.array([2,3,4,5])*lr2.coef_).sum(1)+lr2.intercept_
print 'Python fit RSS',((pred1 - df.growth.values)**2).sum()
print 'Stata fit RSS',((pred2 - df.growth.values)**2).sum()
Python fit RSS 7.2062436549
Stata fit RSS 28.0291464826