多库线性回归不一致问题解析

当使用 scikit-learn、statsmodels 或 R 进行多元线性回归时，若输入特征量纲差异极大（如某特征达 10¹⁸ 级），会导致矩阵病态、浮点精度损失，从而产生截然不同的 R²、系数与截距——这并非算法差异，而是数值计算稳定性问题。

在构建纯预测型多元线性回归模型时，我们期望不同工具（R、scikit-learn、statsmodels）在相同数据上给出高度一致的结果——毕竟它们都求解同一数学问题：最小化残差平方和 $\min_{\beta} |y - X\beta|^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）→ 检查特征量纲 → 必要时标准化 → 再建模”的流程，即可消除此类陷阱，确保模型结果可靠、可复现、可迁移。