
当使用 scikit-learn、statsmodels 或 r 进行多元线性回归时,若输入特征量纲差异极大(如某特征达 10¹⁸ 级),会导致矩阵病态、浮点精度损失,从而产生截然不同的 r²、系数与截距——这并非算法差异,而是数值计算稳定性问题。
在构建纯预测型多元线性回归模型时,我们期望不同工具(R、scikit-learn、statsmodels)在相同数据上给出高度一致的结果——毕竟它们都求解同一数学问题:最小化残差平方和 $min_{beta} |y – Xbeta|^2$。然而,实际中却常出现 R² 差异巨大(如 R 报 0.46,sklearn 报 0.19,statsmodels 甚至报负值)、系数符号/量级完全相悖的现象。根本原因往往不是模型逻辑错误,而是数值不稳定性(numerical instability),其典型诱因是特征尺度严重失衡。
以提问者数据为例,newrh_colint 的取值在 $10^{18}$ 数量级,而其余特征(如 CAPE、shear1)均在 $10^0 sim 10^3$ 范围内。这种跨越 18 个数量级的差异,会使设计矩阵 $X$ 的条件数(condition number)急剧升高,逼近或超过双精度浮点数(float64)的有效位数极限(约 16 位十进制精度)。此时,正规方程 $(X^top X)^{-1}X^top y$ 中的矩阵求逆或 QR 分解过程极易引入显著舍入误差,最终导致:
- sklearn 的 LinearRegression(底层调用 LAPACK dgelsd)因 SVD 截断容差失效而拟合失准;
- statsmodels 的 OLS 因 X.T @ X 严重失真,计算出荒谬的负 R²(本质是残差平方和大于总平方和);
- R 的 lm() 内部使用更稳健的 QR 分解与条件数检测机制,相对更能容忍一定尺度差异,因而结果更可信。
✅ 快速验证与修复方法:
对异常大尺度特征进行简单缩放,例如将 newrh_colint 除以 $10^{18}$:
import numpy as np import pandas as pd from sklearn.linear_model import LinearRegression import statsmodels.api as sm # 构建原始 DataFrame(同提问者) df = pd.DataFrame([ CAPE, AREA, newRHsurf, newLCL, MUCAPE, ECAPE, shear1, newpsurf, newmaxCAPE, newrh_colint, newshear3, CTRRMAX ]).T df.columns = ['CAPE', 'AREA', 'newRHsurf', 'newLCL', 'MUCAPE', 'ECAPE', 'shear1', 'newPsurf', 'newmaxCAPE', 'newrh_colint', 'newshear3', 'CTRRMAX'] # 关键修复:缩放极端特征(10^18 → ~1) df['newrh_colint'] *= 1e-18 X = df.iloc[:, :-1] y = df['CTRRMAX'] # sklearn 拟合(修复后) lm = LinearRegression() model_sk = lm.fit(X, y) print(f"sklearn R²: {model_sk.score(X, y):.4f}") # ≈ 0.458(与 R 接近) # statsmodels 拟合(修复后) X_sm = sm.add_constant(X) res_sm = sm.OLS(y, X_sm).fit() print(f"statsmodels R²: {res_sm.rsquared:.4f}") # ≈ 0.459
? 推荐生产级解决方案:标准化(Standardization)
为避免手动缩放、提升泛化性与可复现性,应统一使用 StandardScaler 对所有特征进行零均值、单位方差变换:
from sklearn.preprocessing import StandardScaler from sklearn.pipeline import Pipeline # 构建带标准化的流水线 pipeline = Pipeline([ ('scaler', StandardScaler()), ('lr', LinearRegression()) ]) model_pipe = pipeline.fit(X, y) print(f"Pipeline R²: {model_pipe.score(X, y):.4f}") # 同样收敛至 ≈0.459
⚠️ 重要注意事项:
- 标准化会改变系数解释:StandardScaler 后的系数 $hat{beta}j$ 表示「当第 $j$ 个特征增加 1 个标准差时,目标变量的平均变化量」,而非原始单位变化;若需原始尺度系数,可通过 `scaler.scale和scaler.mean_` 反向转换;
- 永远不要对目标变量 y 标准化(除非明确需要),否则 score() 返回的是标准化后的 R²,无法直接对比;
- 在交叉验证中,务必在每折内独立拟合 StandardScaler(Pipeline 自动保证),防止数据泄露;
- 若后续需特征重要性排序,建议使用标准化后的系数绝对值,或更稳健的 PermutationImportance。
总结而言,多库结果不一致极少源于“算法不同”,绝大多数情况是数据预处理缺失导致的数值灾难。坚持“先探索性数据分析(EDA)→ 检查特征量纲 → 必要时标准化 → 再建模”的流程,即可消除此类陷阱,确保模型结果可靠、可复现、可迁移。