R和Python中线性回归的差异[关闭]

问题描述 投票：2回答：2

我试图将线性回归R结果与python的结果相匹配

匹配每个自变量的系数

以下是代码：

数据已上传。

https://www.dropbox.com/s/oowe4irm9332s78/X.csv?dl=0

https://www.dropbox.com/s/79scp54unzlbwyk/Y.csv?dl=0

R代码：

#define pathname = " "

X <- read.csv(file.path(pathname,"X.csv"),stringsAsFactors = F)
Y <- read.csv(file.path(pathname,"Y.csv"),stringsAsFactors = F)

d1 = Y$d1

reg_data <- cbind(d1, X)
head(reg_data)

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

names(which(is.na(reg_model$coefficients)))

reg_model$coefficients

R结果

> summary(reg_model)

Call:
lm(formula = d1 ~ ., data = reg_data)

Residuals:
ALL 60 residuals are 0: no residual degrees of freedom!

Coefficients: (14 not defined because of singularities)
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -67.37752         NA      NA       NA
v1            1.30214         NA      NA       NA
v2           -2.93118         NA      NA       NA
v3            7.58902         NA      NA       NA
v4           11.88570         NA      NA       NA
v5            1.60622         NA      NA       NA
v6            3.71528         NA      NA       NA
v7           -9.34627         NA      NA       NA
v8           -3.84694         NA      NA       NA
v9           -2.51332         NA      NA       NA
v10           4.22403         NA      NA       NA
v11          -9.70126         NA      NA       NA
v12                NA         NA      NA       NA
v13           4.67276         NA      NA       NA
v14          -6.57924         NA      NA       NA
v15          -3.68065         NA      NA       NA
v16           5.25168         NA      NA       NA
v17          14.60444         NA      NA       NA
v18          16.00679         NA      NA       NA
v19          24.79622         NA      NA       NA
v20          13.85774         NA      NA       NA
v21           2.16022         NA      NA       NA
v22         -36.65361         NA      NA       NA
v23           2.26554         NA      NA       NA
v24                NA         NA      NA       NA
v25                NA         NA      NA       NA
v26           7.00981         NA      NA       NA
v27           0.88904         NA      NA       NA
v28           0.34400         NA      NA       NA
v29          -5.27597         NA      NA       NA
v30           5.21034         NA      NA       NA
v31           6.79640         NA      NA       NA
v32           2.96346         NA      NA       NA
v33          -1.52702         NA      NA       NA
v34          -2.74632         NA      NA       NA
v35          -2.36952         NA      NA       NA
v36          -7.76547         NA      NA       NA
v37           2.19630         NA      NA       NA
v38           1.63336         NA      NA       NA
v39           0.69485         NA      NA       NA
v40           0.37379         NA      NA       NA
v41          -0.09107         NA      NA       NA
v42           2.06569         NA      NA       NA
v43           1.57505         NA      NA       NA
v44           2.70535         NA      NA       NA
v45           1.17634         NA      NA       NA
v46         -10.51141         NA      NA       NA
v47          -1.15060         NA      NA       NA
v48           2.87353         NA      NA       NA
v49           3.37740         NA      NA       NA
v50          -5.89816         NA      NA       NA
v51           0.85851         NA      NA       NA
v52           3.73929         NA      NA       NA
v53           4.93265         NA      NA       NA
v54           3.45650         NA      NA       NA
v55           0.12382         NA      NA       NA
v56          -0.21171         NA      NA       NA
v57           4.37199         NA      NA       NA
v58           3.21456         NA      NA       NA
v59           0.09012         NA      NA       NA
v60          -0.85414         NA      NA       NA
v61          -3.29856         NA      NA       NA
v62           4.38842         NA      NA       NA
v63                NA         NA      NA       NA
v64                NA         NA      NA       NA
v65                NA         NA      NA       NA
v66                NA         NA      NA       NA
v67                NA         NA      NA       NA
v68                NA         NA      NA       NA
v69                NA         NA      NA       NA
v70                NA         NA      NA       NA
v71                NA         NA      NA       NA
v72                NA         NA      NA       NA
v73                NA         NA      NA       NA

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:    NaN 
F-statistic:   NaN on 59 and 0 DF,  p-value: NA

Python代码：

Y = pd.read_csv(pathname+"Y.csv")
X = pd.read_csv(pathname+"X.csv")

lr = LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
lr.fit(X, Y['d1'])

(list(zip(lr.coef_, X)))

lr.intercept_

Python结果：

intercept = 29.396033164254106

[(-2.4463986167806304, 'v1'),
 (-1.6293010275307021, 'v2'),
 (0.89089949009506508, 'v3'),
 (-3.1021251646895251, 'v4'),
 (-1.7707078771936109, 'v5'),
 (-2.0474705122225636, 'v6'),
 (-1.5537181337496202, 'v7'),
 (-1.6391241229716156, 'v8'),
 (-1.2981646048517046, 'v9'),
 (0.89221826294889328, 'v10'),
 (-0.56694104645951571, 'v11'),
 (2.042810365310288e-14, 'v12'),
 (-2.0312478672439052, 'v13'),
 (-1.5617121392788413, 'v14'),
 (0.4583365939498274, 'v15'),
 (0.8840538748922967, 'v16'),
 (-5.5952681002058871, 'v17'),
 (2.4937042448512892, 'v18'),
 (0.45806845189176543, 'v19'),
 (-1.1648810657830406, 'v20'),
 (-1.7800004329275585, 'v21'),
 (-5.0132817522704816, 'v22'),
 (3.6862778096189266, 'v23'),
 (2.7533531010703882e-14, 'v24'),
 (1.2150003225741557e-14, 'v25'),
 (0.94669823515018103, 'v26'),
 (-0.3082823207975679, 'v27'),
 (0.53619247380957358, 'v28'),
 (-1.1339902793546781, 'v29'),
 (1.9657159583080186, 'v30'),
 (-0.63200501460653324, 'v31'),
 (1.4741013580918978, 'v32'),
 (-2.4448418291953313, 'v33'),
 (-2.0787115960875036, 'v34'),
 (0.22492914212063603, 'v35'),
 (-0.75136276693004922, 'v36'),
 (1.2838658951186761, 'v37'),
 (0.5816277993227944, 'v38'),
 (-0.11270569554555088, 'v39'),
 (-0.13430982360936233, 'v40'),
 (-3.3189296496897662, 'v41'),
 (-0.452575588270415, 'v42'),
 (6.1329755709937519, 'v43'),
 (0.61559185634634817, 'v44'),
 (-1.206315459828555, 'v45'),
 (-3.7452010299772009, 'v46'),
 (-1.1472174665136678, 'v47'),
 (2.8960489381172652, 'v48'),
 (0.0090220136972478659, 'v49'),
 (-5.264918363314754, 'v50'),
 (1.2194758337662015, 'v51'),
 (2.78655271320092, 'v52'),
 (3.106513852668896, 'v53'),
 (3.5181252502607929, 'v54'),
 (-0.34426523770507278, 'v55'),
 (-0.48792823932479878, 'v56'),
 (0.12284460490031779, 'v57'),
 (1.6860388628044991, 'v58'),
 (1.2823067194737174, 'v59'),
 (2.8352263554153665, 'v60'),
 (-1.304336378501032, 'v61'),
 (0.55226132316435139, 'v62'),
 (1.5416988124754771, 'v63'),
 (-0.2605804175310813, 'v64'),
 (1.2489066081702334, 'v65'),
 (-0.44469553013696161, 'v66'),
 (-1.4102990055550157, 'v67'),
 (3.8150423259022639, 'v68'),
 (0.12039684410168072, 'v69'),
 (-1.340699466779357, 'v70'),
 (1.7066389124439212, 'v71'),
 (0.50470752944860442, 'v72'),
 (1.0024872633969766, 'v73')]

但它不匹配。请帮助。

注意：它匹配以下示例

http://davidcoallier.com/blog/linear-regression-from-r-to-python/

python r scikit-learn linear-regression lm
2个回答
10
投票

tl;博士如果你想在Python中复制R的策略你可能不得不自己实现它，因为R做了一些其他地方没有广泛使用的聪明的东西。

作为参考（因为到目前为止仅在评论中提到），这是一个不适合/等级不足的拟合，当有比预测更多的预测变量时，它总会发生（p> n：在这种情况下p = 73，n = 61），并且通常当存在许多分类响应和/或实验设计以某种方式受到限制时。处理这些情况以获得完全意味着任何事情的答案通常需要仔细考虑和/或高级技术（例如，惩罚性回归：参见Wikipedia article on linear regression中的套索和岭回归的参考）。

处理这种情况最天真的方法是把它全部扔进标准的线性代数中，希望没有什么破坏得太厉害，这显然是python的statsmodels包所做的：来自pitfalls document：

  • 排名不足的矩阵不会引起错误。
  • 几乎完美的多重共线性或病态设计矩阵的情况可能会产生数值不稳定的结果。如果这不是期望的行为，则用户需要手动检查矩阵的等级或条件数。

下一个最好的事情（当存在较小程度的共线性时是合理的）是在进行线性代数时明显地转动，即重新排列计算问题，以便可以省略共线部分。这就是R的作用;为了做到这一点，R代码的作者必须修改标准的LINPACK例程。

underlying code有这样的评论/解释：

c dqrdc2使用户主变换来计算qr c由p矩阵x对n进行分解。有限的专栏 c基于缩减列的2范数的旋转策略 c将具有接近零范数的列移动到右边缘 c x矩阵。这个策略意味着顺序一个 c自由度效应可以以自然的方式计算。

我非常担心以这种方式修改linpack代码。 c如果你是计算线性代数大师，你真的 c了解如何解决这个问题请随意 c建议改进此代码。

这段代码（和评论）自1998年以来一直在R的代码库中;我很想知道最初是谁写的（基于comments further down in the code，它似乎是Ross Ihaka？），但是在代码的历史回溯到1998年的代码重组之后我遇到了麻烦。（更多的挖掘表明这个代码有自从其记录历史开始以来，它一直存在于R的代码库中，即该文件是在SVN修订版2 1997-09-18中添加的，直到很久之后才被修改。）

MartinMächler最近（2016年10月25日，here）添加了有关?qr的更多信息，因此这些信息实际上可以在下一版R的文档中找到...

如果您知道如何将已编译的FORTRAN代码与Python代码链接（我没有），那么编译src/appl/dqrdc2.f并将lm.fit的内容翻译成Python将非常容易：这是lm.fit的核心，减去错误检查和其他处理...

z <- .Call(C_Cdqrls, x, y, tol, FALSE)
coef <- z$coefficients
pivot <- z$pivot
r1 <- seq_len(z$rank)
dn <- colnames(x)
nmeffects <- c(dn[pivot[r1]], rep.int("", n - z$rank))
r2 <- if (z$rank < p) 
    (z$rank + 1L):p
else integer()
if (is.matrix(y)) { ## ... 
}  else {
    coef[r2] <- NA
    if (z$pivoted) 
        coef[pivot] <- coef
    names(coef) <- dn
    names(z$effects) <- nmeffects
}
z$coefficients <- coef
2
投票

作为Ben Bolker答案的补充。

主要问题是统计软件包应该对病态问题做些什么。据我所知，在处理奇异和几乎单一的设计矩阵时，包之间存在很大差异。唯一确定性的方法是用户明确选择变量选择或惩罚算法。

如果可以清楚地识别排名，那么结果仍然是确定性的，但是根据所选择的奇点策略而变化。 Stata和R drop变量，Stata在列出变量时按顺序丢弃它们，即丢弃最后一个共线变量。我不知道哪个变量R掉了。 statsmodels使用对称处理变量并通过使用基于奇异值分解的广义逆pinv来丢弃奇异值。这相当于PCA降低等级回归或岭回归的微小惩罚。

结果不是确定性的，即取决于线性代数包，并且如果数值噪声影响等级或共线变量的识别，则在不同的计算机上甚至可能不相同。具体地，显示QR和pinv / SVD的等级需要阈值，低于该阈值，等级被识别为减少。在statsmodels中，相对条件数的numpy默认阈值约为1e-15。因此，如果奇异值小于此值，我们将得到正则化解。但在某些情况下，阈值可能太低，而“非奇异”解决方案则受到无法复制的数值噪声的支配。我想对于显示QR或其他纯粹数值解决共线性问题的任何等级都是如此。 （Robustness issue of statsmodel Linear regression (ols) - Python https://github.com/statsmodels/statsmodels/issues/2628 和相关的 statmodels in python package, How exactly duplicated features are handled?）

关于排名揭示，转动QR：

当大多数statsmodels.OLS被编写时，它在scipy中不可用，但它已经在scipy中作为pivoting关键字提供了一段时间，但尚未作为statsmodels的选项添加.OLS https://github.com/scipy/scipy/pull/44

我对它作为一种默认解决方案持怀疑态度，因为我猜测，如果没有验证它，数值噪声会影响哪些变量被转动。然后，它将不再是确定性的，将丢弃哪些变量。在我看来，变量选择应该是用户有意识的选择，而不是纯粹的数字选择。

（免责声明：我是一名statsmodels维护者）

编辑：

该问题在示例中使用了scikit-learn。

据我所知，scikit-learn使用不同的LAPACK函数来解决线性回归问题，但它也基于像statsmodels这样的奇异值分解。但是，scikit-learn目前使用相应scipy函数的默认阈值，该阈值小于statsmodels使用的numpy函数。 例如how does sklearn do Linear regression when p >n?

因此，我希望scikit-learn和statsmodels在单数情况下具有相同的结果，但结果在一些近似单一情况下会有所不同。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.