Python和R之间的线性回归系数之间的差异

Posted 2023-02-23

【中文标题】Python和R之间的线性回归系数之间的差异

【英文标题】：Difference between Linear Regression Coefficients between Python and R

【发布时间】：2017-09-17 09:31:03

【问题描述】：

我正在尝试在 Python 中运行我已经在 R 中完成的线性回归，以便找到系数为 0 的变量。我遇到的问题是 R 中的线性回归为低方差的列返回 NA，而 scikit 学习回归返回系数。在 R 代码中，我通过将带有 NA 的变量保存为线性回归的输出来查找并保存这些变量，但我似乎无法找到一种在 python 中模仿这种行为的方法。我正在使用的代码可以在下面找到。

R 代码：

a <- c(23, 45, 546, 42, 68, 15, 47)
b <- c(1, 2, 4, 6, 34, 2, 8)
c <- c(22, 33, 44, 55, 66, 77, 88)
d <- c(1, 1, 1, 1, 1, 1, 1)
e <- c(1, 1, 1, 1, 1, 1, 1.1)
f <- c(1, 1, 1, 1, 1, 1, 1.01)
g <- c(1, 1, 1, 1, 1, 1, 1.001)

df <- data.frame(a, b, c, d, e, f, g)
var_list = c('b', 'c', 'd', 'e', 'f', 'g')

target <- temp_dsin.df$a
reg_data <- cbind(target, df[, var_list])


if (nrow(reg_data) < length(var_list))
  message(paste0('    WARNING: Data set is rank deficient. Result may be doubtful'))

reg_model <- lm(target ~ ., data = reg_data)

print(reg_model$coefficients)

#store the independent variables with 0 coefficients
zero_coef_IndepVars.v <- names(which(is.na(reg_model$coefficients)))

print(zero_coef_IndepVars.v)

Python 代码：

import pandas as pd
from sklearn import linear_model

a = [23, 45, 546, 42, 68, 15, 47]
b = [1, 2, 4, 6, 34, 2, 8]
c = [22, 33, 44, 55, 66, 77, 88]
d = [1, 1, 1, 1, 1, 1, 1]
e = [1, 1, 1, 1, 1, 1, 1.1]
q = [1, 1, 1, 1, 1, 1, 1.01]
f = [1, 1, 1, 1, 1, 1, 1.001]


df = pd.DataFrame('a': a,
                             'b': b,
                             'c': c,
                             'd': d,
                             'e': e,
                             'f': q,
                             'g': f)


var_list = ['b', 'c', 'd', 'e', 'f', 'g']

# build linear regression model and test for linear combination
target = df['a']
reg_data = pd.DataFrame()
reg_data['a'] = target
train_cols = df.loc[:,df.columns.str.lower().isin(var_list)]


if reg_data.shape[0] < len(var_list):
    print('    WARNING: Data set is rank deficient. Result may be doubtful')

# Create linear regression object
reg_model = linear_model.LinearRegression()

# Train the model using the training sets
reg_model.fit(train_cols , reg_data['a'])

print(reg_model.coef_)

R 的输出：

(Intercept)           b           c           d           e           f           g 
 537.555988   -0.669253   -1.054719          NA -356.715149          NA          NA 

> print(zero_coef_IndepVars.v)
[1] "d" "f" "g"

Python 的输出：

           b             c   d               e              f            g
[-0.66925301   -1.05471932   0.   -353.1483504   -35.31483504   -3.5314835]

如您所见，“b”、“c”和“e”列的值很接近，但“d”、“f”和“g”列的值却大不相同。对于这个示例回归，我想返回 ['d', 'f', 'g'] 因为它们的输出是来自 R 的 NA。问题是 sklearn 线性回归为 col 'd' 返回 0，而它返回-35.31 用于 col 'f' 和 -3.531 用于 col 'g'。

有谁知道 R 如何决定是返回 NA 还是返回值/如何将这种行为实现到 Python 版本中？了解差异的来源可能会帮助我在 python 中实现 R 行为。我需要 python 脚本的结果与 R 输出完全匹配。

【问题讨论】：

请注意他们是NA，而不是NaN。

@SinanÜnür 所以你认为在 R 线性回归中存在共线性检查吗？我认为它是这样的，这就是为什么我以我的方式选择我的数据，但我需要在 python 中复制这种行为。

【参考方案1】：

这是实现方式的不同。 R 中的lm 使用基于 QR 分解的底层 C 代码。模型矩阵被分解为正交矩阵 Q 和三角矩阵 R。这导致了其他人所说的“共线性检查”。 R 不检查这一点，QR 分解的性质确保最小共线变量在拟合算法中获得“优先级”。

有关线性回归上下文中 QR 分解的更多信息： https://www.stat.wisc.edu/~larget/math496/qr.html

来自 sklearn 的代码基本上是 numpy.linalg.lstsq 的包装，它最小化了欧几里得二次范数。如果您的模型是Y = AX，它会最小化||Y - AX||^2。这是一种不同的（计算上不太稳定的）算法，它没有 QR 分解的副作用。

个人注意：如果您想在经过验证和测试的计算框架中稳健地拟合模型并坚持使用 Python，请寻找基于 QR 或 SVD 的线性回归实现。软件包 scikit-learn 或 statsmodels（截至 2017 年 4 月 22 日仍处于测试阶段）应该可以帮助您。

【讨论】：

谢谢大家，这非常有帮助！我会看一下python中是否存在任何使用QR分解实现的东西，或者如果那不起作用，则坚持使用R。

个人注意：如果你想从回归中随机删除一些变量，那么使用 R。如果你想要一个 SVD/pinv 正则化解决方案，那么使用 Python scikit-learn 或 statsmodels。如果您都不想要，请清理您的数据并自己选择变量。

@user333700 QR 分解没有随机性，如问题所示。但无论如何，感谢您提供指向基于 SVD 的方法的指针。记录一下：我们部门一半的研究是用 Python 完成的。 Python 是一门很棒的语言。但是对于一个简单的线性回归，我们现在实际上已经提到了 5 个包（numpy、scipy、sklearn、scikit-learn 和 statsmodels）。我会在我的课堂统计计算中坚持使用 R：少做解释，标准工具是一个稳定的工具。我猜每个人都有自己的想法。

另请参阅***.com/questions/40935624/… 中的答案以获取更多详细信息。 （即使现在可以使用 scipy.linalg 来实现，statsmodels 也不太可能通过旋转 QR 来获得变量选择。）

我用这个例子玩了一下：如果我们通过添加 1e-6 来更改d，即d <- c(1, 1, 1, 1, 1, 1, 1 + 1e-6)，那么d 被保留，e 被删除。如果我添加 1e-7 而不是 1e-6，那么 d 仍会像原来的情况一样被删除。对于数据集略有不同的重复用例，这对我来说看起来很“随机”（脆弱）。

【参考方案2】：

我猜数据不够。这是 statsmodel 的结果：

import statsmodels.formula.api as smf
lm = smf.ols(formula='a ~ b + c + d + e + f + g', data=df).fit()
lm.summary()

给：

OLS Regression Results
Dep. Variable:	a	R-squared:	0.038
Model:	OLS	Adj. R-squared:	-0.923
Method:	Least Squares	F-statistic:	0.03993
Date:	Fri, 21 Apr 2017	Prob (F-statistic):	0.987
Time:	22:29:16	Log-Likelihood:	-46.059
No. Observations:	7	AIC:	100.1
Df Residuals:	3	BIC:	99.90
Df Model:	3		
Covariance Type:	nonrobust		
coef	std err	t	P>|t|	[95.0% Conf. Int.]
Intercept	151.5350	1065.536	0.142	0.896	-3239.476 3542.545
b	-0.6693	10.324	-0.065	0.952	-33.526 32.188
c	-1.0547	6.412	-0.164	0.880	-21.462 19.352
d	151.5350	1065.536	0.142	0.896	-3239.476 3542.545
e	-368.1353	3862.592	-0.095	0.930	-1.27e+04 1.19e+04
f	99.5679	574.110	0.173	0.873	-1727.506 1926.642
g	146.3383	1016.341	0.144	0.895	-3088.111 3380.788
Omnibus:	nan	Durbin-Watson:	2.447
Prob(Omnibus):	nan	Jarque-Bera (JB):	4.545
Skew:	1.797	Prob(JB):	0.103
Kurtosis:	4.632	Cond. No.	1.34e+18

OLS 提供了几个线索，表明这个线性问题是病态的。