R、SciKitLearn 和 Statsmodels 提供不同的多重线性回归拟合

Question

我正在尝试建立一个多元线性回归模型来预测一个名为 CTRRMAX 的量。这个 MLR 模型纯粹是一个预测模型 - 我不需要了解预测变量对结果变量的重要性或关系。我只是希望模型能够解释尽可能多的方差。

我学习了如何在 R 中进行回归分析，因此我在那里进行了大部分转换和特征选择。当我尝试将其嵌入到我通常使用的 Python 代码中时，scikitlearn 为我提供了一个完全不同的模型，statsmodels 也是如此。我不明白为什么我会从所有这些不同的包中得到不同的 R 平方、系数和截距。它们不应该都一样吗？

我已经（尽我所能）提供了不同包之间的代码。

R代码：

final_MLRmodel = lm(ctrrmax~cape+area+newrhsurf+newlcl+mucape+ecape+shear1+newpsurf+newmaxcape+newrhcolint+newshear3)
summary(final_MLRmodel)

输出给我 R2=0.4593 和 adjR2=0.3999

ScikitLearn 代码：

from sklearn.linear_model import LinearRegression

df = pd.DataFrame([CAPE, AREA, newRHsurf, newLCL, MUCAPE, ECAPE, shear1, newPsurf, newmaxCAPE, newrh_colint, newshear3, CTRRMAX])
df = df.transpose()

X = df[df.columns[0:-1]]
y = df[df.columns[-1]]

lm = LinearRegression()
model = lm.fit(X, y)

y_pred = model.predict(X)
R2 = model.score(X,y)
print('R2: ', model.score(X,y))

n = len(CTRRMAX)
p=np.shape(X)[1]
adjr2 = 1-(1-R2)*(n-1)/(n-p-1)
print('adjR2: ', adjr2)

输出给出 R2=0.1855 和 adjR2=0.0959。

统计模型：

import statsmodels.api as sm
X1 = sm.add_constant(X)
mod = sm.OLS(y,X1)
res = mod.fit()
res.summary()

输出给出 R2=-7.283 和 adjR2=-7.283。

我已经确认了一百万次，即使在转换之后，所有模型都在读取完全相同的数据。我想这可能是多重共线性的问题，但如果我只想要一个预测模型，那么我认为我不必关心多重共线性？任何帮助将不胜感激。

这个模型显然也需要一个数据集来运行。这是我的第一个堆栈溢出问题，我不知道如何附加用于将数据从 Python 传输到 R 的 .csv 文件。我在这里仅发布了每个变量的前 20 个数据点，在希望不需要整个数据集。

CTRRMAX = np.array(0.65109549, 0.57082317, 0.66105562, 0.65956423, 0.62482739,
       0.66715294, 0.58653562, 0.67174625, 0.64250845, 0.48642454,
       0.77204295, 0.66184076, 0.68038761, 0.73653454, 0.60328276,
       0.62336814, 0.64258676, 0.65066922, 0.82871873, 0.58368109)
CAPE = np.array(2160., 3030., 2090., 1420., 2540.,  470., 1730., 2310., 2010.,
       5020., 1490., 2830.,  510., 1310., 3800., 1700., 1660., 1990.,
        130., 3080.)
AREA = np.array(266., 175., 122., 161., 165., 225.,  86., 129., 295., 316.,  44.,
       116., 244., 253., 118., 138., 128.,  69., 118., 117.)
newRHsurf = np.array(0.37423277, -0.38957954, -0.00965352, -0.74678448,  0.01749659,
        0.78598837,  0.40301939, -0.99305805,  0.57022447,  0.99996833,
       -0.40080868,  0.59670993,  0.8575545 , -0.26668987,  0.81043055,
       -0.15807503,  0.7954259 , -0.86284556, -0.23866762,  0.1054419)
newLCL = np.array(-0.99553022,  0.09343053,  0.80974088, -0.99510076, -0.95194702,
        0.16104252,  0.8352307 ,  0.99146248, -0.83244151, -0.25468422,
       -0.78162311,  0.13424665, -0.55608987, -0.31989455, -0.4858867 ,
        0.78814695, -0.3772229 , -0.14639698, -0.49059305,  0.9434852)
MUCAPE = np.array(3231.9094098 , 2722.89503908, 2434.98996265, 2483.8618223 ,
       1534.22930423, 2274.16312798, 1758.04659665, 2013.51926427,
       1852.74436288, 3849.84084006, 1538.14989687, 2040.98226786,
       1049.41283138, 1623.10054915, 4055.06571398, 2386.93537232,
       2421.36562901, 3093.63218601, 1886.79025799, 3033.74984899)
ECAPE = np.array(544.31677561, 628.24486587, 198.8664456 ,  27.78930588,
       726.19478451,  28.43822235, 180.08772011, 267.41910501,
       205.58738719, 282.76765822, 140.42604251, 453.39778625,
       110.01754724, 173.23091256, 546.4515334 ,  28.88800818,
        32.49516751, 194.49062916,  47.18208684, 150.4724515)
shear1 = np.array(0.38087247, 1.57931439, 1.39215283, 1.10672164, 3.02872954,
       1.80076502, 0.9350142 , 1.9858369 , 0.56996634, 3.39201731,
       1.86925209, 1.74190136, 0.40244472, 2.16668611, 2.02204818,
       1.50670497, 1.07298262, 0.74519335, 1.04231893, 1.52101836)
newPsurf = np.array(0.12141161,  0.99513855, -0.47831087,  0.37072577, -0.92267433,
        0.30656576, -0.22574095, -0.99994693, -0.45593397,  0.95358944,
       -0.9355281 , -0.67766388, -0.84796369,  0.69423688,  0.99605423,
        0.86181924, -0.30228107, -0.3340524 ,  0.02340334,  0.99402057)
newmaxCAPE = np.array(0.59753941,  0.63973417, -0.88611046,  0.97564725, -0.56573558,
        0.99420651,  0.69055665,  0.18575401,  0.94000581,  0.32183912,
        0.98409834, -0.99429748,  0.96010656,  0.7027577 ,  0.97306961,
        0.9949444 ,  0.99237209, -0.97644958,  0.04583703,  0.99022841)
newrh_colint = np.array(4.35888606e+18, 3.47399496e+18, 3.33191607e+18, 3.72920308e+18,
       3.61121221e+18, 3.17535430e+18, 3.59729890e+18, 2.86802689e+18,
       3.44806618e+18, 5.08126931e+18, 2.34607596e+18, 3.28603005e+18,
       2.03731238e+18, 2.78882246e+18, 4.56923554e+18, 3.21237248e+18,
       3.41041888e+18, 3.53788584e+18, 2.80536323e+18, 4.32299498e+18)
newshear3 = np.array(0.44118218, -0.99932764,  0.42862348,  0.65427145, -0.94152187,
        0.86240296,  0.38829232, -0.33614615,  0.83542138, -0.46978768,
        0.6141936 , -0.61328793,  0.77510501,  0.8593079 , -0.03691697,
        0.99619598,  0.55651457,  0.97479817,  0.99382273,  0.23179711)

Answer 1

我怀疑您遇到了变量缩放问题。

newrh_colint

变量比其他变量大 10^18 倍。

如果我加载你的变量来预测，像这样：

df = pd.DataFrame([CAPE, AREA, newRHsurf, newLCL, MUCAPE, ECAPE, shear1, newPsurf, newmaxCAPE, newrh_colint, newshear3, CTRRMAX])
df = df.transpose()
df[9] *= 1e-18

然后 newrh_colint 的大小与其他变量大致相同，并且我得到的 R 平方更类似于 R 生成的值。

为什么这很重要？理论上，浮点数在 10^18 时的相对精度与在 10^0 时的相对精度一样高，并且将 X 变量之一乘以常数不会改变结果。但是，我怀疑 newrh_colint 列比其他列大得多，以至于结果几乎是奇异矩阵。

如果您想要一个自动缩放变量以使它们具有相同大小的解决方案，您可以使用StandardScaler。