我正在将 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
------------------------------------------------------------------------------
最佳答案
预测变量是 time
的 1 次... 5 次方,它在 1627 年和 2007 年之间变化(大概是一个日历年,这无关紧要)。即使使用现代软件,改变时间原点以减少数值应变也是明智的,例如使用 (time
- 1800) 的幂。
无论如何,重做回归表明 Stata 将第一个预测变量丢弃为共线。 Python 和 R 会发生什么?这些是对数字棘手挑战的不同 react 。
(拟合五次多项式很少有科学值(value),但这在这里可能无关紧要。基于 2 到 5 次方的拟合曲线对于这些数据效果不是很好,这些数据看起来很经济。更有意义的是前 5 个残差都是正数,但并非所有残差都是正数!)
关于python - Python(和 R)和 Stata 中的线性回归之间的区别,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/32280444/