我正在将 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
------------------------------------------------------------------------------
最佳答案
预测变量是 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/