如何从Sklearn的HuberRegressor鲁棒回归中获取r值与p值

首先明确：Sklearn内置的HuberRegressor接口没有直接输出r值、p值的参数，需要手动计算，且计算得到的指标统计含义和普通最小二乘（OLS）的对应指标存在差异，仅可作为显著性参考，不能完全套用OLS的显著性判断逻辑。

1. 计算r值 #

你之前通过stats.linregress得到的r_value为皮尔逊相关系数，要和该指标对齐，可直接计算Huber回归预测值与真实值的皮尔逊相关系数，也可根据需求计算衡量模型解释力的决定系数R²：

import numpy as np
from scipy import stats
from sklearn.metrics import r2_score

# 提取真实值，用训练好的Huber模型生成预测值
y_true = mn_all_df['Y'].to_numpy()
y_pred = huber.predict(mn_all_df['X'].to_numpy().reshape(-1, 1))

# 计算与linregress输出对齐的皮尔逊r值
r_value = np.corrcoef(y_true, y_pred)[0, 1]

# 可选：计算决定系数R²
r2 = r2_score(y_true, y_pred)

2. 计算p值 #

Huber回归没有标准的p值闭式解，工业界和学术界常用大样本下的Wald检验近似计算，实现逻辑如下：

# 获取模型拟合的核心参数
coef = huber.coef_[0]
epsilon = huber.epsilon
scale = huber.scale_
X_arr = mn_all_df['X'].to_numpy()
n_sample = len(X_arr)
n_param = 2 # 本场景为单变量回归，参数为截距+自变量系数两个

# 构造带截距项的自变量矩阵
X_mat = np.vstack([np.ones(n_sample), X_arr]).T
# 计算残差与Huber权重
residuals = y_true - y_pred
huber_weights = np.where(np.abs(residuals) < epsilon * scale,
                         1,
                         epsilon * scale / np.abs(residuals))
W = np.diag(huber_weights)

# 用三明治估计法计算系数的鲁棒协方差矩阵
bread = np.linalg.inv(X_mat.T @ W @ X_mat)
meat = X_mat.T @ np.diag(residuals ** 2) @ X_mat
cov_mat = bread @ meat @ bread

# 提取自变量系数的标准误，构造t统计量
coef_se = np.sqrt(np.diag(cov_mat))[1]
t_stat = coef / coef_se
# 计算双尾检验p值，自由度为样本量减参数个数
p_value = 2 * (1 - stats.t.cdf(np.abs(t_stat), df = n_sample - n_param))

注意事项 #

• 上述p值为大样本下的近似结果，样本量过小时偏差会明显升高
• 若对p值精度要求极高，可采用Bootstrap抽样法多次拟合Huber回归，基于系数的经验分布计算p值，仅计算成本会更高
• 你之前拟合时传入的Y为二维数组，HuberRegressor的fit方法默认接收一维Y，可调整为huber.fit(X, y.ravel())避免接口警告

内容的提问来源于stack exchange，提问作者Dong-gyun Kim