R语言稳健回归实战：从lm到rlm的异常值诊断与处理

简介R语言稳健性估计实例分析资源面向数据分析、统计建模及回归诊断学习者。压缩包共1个pptx文件大小仅716KB以幻灯形式系统展示线性回归诊断与稳健回归的完整思路。内容从lm()基础拟合与plot()四联诊断图出发逐步讲解残差、异常点、高杠杆点与强影响点的判别方法涵盖学生化残差、帽子矩阵及Cook距离等关键指标的计算与应用同时理清三类特殊点的联系与区别在此基础上引入Huber与Bisquare两种M估计稳健回归方法通过实例演示如何在异常值存在时进行加权迭代获得更可靠的参数估计。整体框架紧凑适合教学演示、课后复习或项目参考能帮助读者快速构建回归稳健性分析的知识体系。目前已有1129人学习对于需要处理含异常值数据的分析人员具有较高参考价值。1. R 语言稳健性估计从 lm() 到 rlm() 的完整实例分析做回归分析时我经常碰到一种场景数据里混进了几个“不老实”的点普通最小二乘回归OLS的结果被它们牵着鼻子走模型系数变得面目全非。R 语言里处理这类问题有一套成熟的工具链从lm()拟合、plot(lm.fit1)出四张诊断图到cooks.distance()计算 Cook 距离再到稳健回归中的 Huber 和 Bisquare M 估计每一步都有对应的函数和判断标准。这篇文章围绕一套完整的 R 实例分析展开包含可直接运行的 R 代码和一份 crime 数据集的分析流程适合正在做回归诊断、异常值处理或需要提高模型稳健性的数据分析师和统计专业学生。你将看到普通残差、学生化残差、杠杆率、Cook 距离这几个概念如何串成一条识别异常点的完整链路以及rlm()在 Huber 和 Bisquare 两种权重函数下的实际表现——这些内容在多数教材里只讲公式很少告诉你参数怎么选、输出怎么读、哪些“经验分界点”其实有争议。文章会以一份真实可跑通的 R 代码为主线把每个函数的作用、每段输出的含义、每个阈值的由来都拆开讲清楚。2. 从普通残差到学生化残差异常点的识别逻辑与帽子矩阵2.1 普通残差为什么不能直接用方差不等齐问题任何一本回归分析教材都会告诉你残差是观测值Y与预测值Ŷ的差表达式为e Y - Ŷ。但实际用 R 做诊断时直接比较普通残差的大小是有问题的。问题出在方差上普通残差的方差不是常数它依赖于帽子矩阵的对角线元素h_ii具体形式是Var(e_i) σ²(1 - h_ii)。这意味着什么不同观测点的残差天然具有不同的方差如果直接比较e_i的绝对值大小那些h_ii较大的点即远离自变量均值的点残差方差更小同样的偏差会被放大从而被误判为异常点。我一般会在 R 里这样获取普通残差# 读取数据并拟合普通线性回归模型 c1 - read.csv(E:/RData/20170917.csv) attach(c1) lm.fit1 - lm(Weight ~ Height, data c1) # 提取普通残差和拟合值 resid_ols - resid(lm.fit1) fitted_ols - fitted(lm.fit1) # 查看前六个残差 head(resid_ols)这段代码中resid()函数提取 OLS 回归的普通残差fitted()提取模型对每个样本的预测值。attach(c1)把数据框的列变量直接暴露到工作环境中方便后续直接引用Weight和Height但要注意使用后建议用detach(c1)释放避免变量名冲突。plot(lm.fit1)是诊断的第一道工序它一次生成四幅图残差对拟合值图、残差的正态 Q-Q 图、标准化残差绝对值平方根对拟合值图、Cook 距离图。这里面第三幅图横轴是拟合值纵轴是sqrt(|standardized residuals|)主要用来检查方差齐性。如果你看到散点呈现漏斗形分布说明方差不稳定这时普通残差的可比性进一步下降。2.2 帽子矩阵与杠杆率h_ii 如何刻画点的“偏远程度”杠杆率衡量的是自变量X对自身均值的偏异程度。公式为h_ii (1/n) (X_i - X̄) (XX)^{-1} (X_i - X̄)从公式可以直接读出两层含义第一项1/n是基础杠杆所有点共享第二项是第i个点到样本中心X̄的 Mahalanobis 距离。在样本空间中h_ii较大的点位于自变量空间的边缘它们可能把回归线拉向自己对回归系数的 LS 估计影响可能很大。在 R 中提取杠杆率有很多路径常见做法是# 通过 lm.influence 获取帽子矩阵对角线元素 H - hatvalues(lm.fit1) # 查看杠杆率最高的几个样本 head(sort(H, decreasing TRUE), 5) # 结合模型矩阵手动计算杠杆率 X - model.matrix(lm.fit1) H_manual - diag(X %*% solve(t(X) %*% X) %*% t(X))代码中hatvalues()返回帽子矩阵的对角线元素是官方推荐做法。model.matrix()提取设计矩阵X包括截距列和自变量列然后用矩阵运算手动复现X(XX)^{-1}X的对角线。手动计算的目的是验证对帽子矩阵的理解实际项目中直接用hatvalues()即可。注意1/n这一项说明即使所有自变量都等于均值杠杆率也至少是1/n所以看杠杆率时不要只看绝对值还要结合2p/n或3p/n这类经验阈值判断。2.3 学生化残差的计算公式拆解与 R 实现由于普通残差存在方差不齐的问题需要标准化后比较。学生化残差的形式是r_i e_i / (s * sqrt(1 - h_ii))其中s是剩余标准差h_ii是帽子矩阵对角线元素。从公式可以看出学生化残差同时考虑了残差本身的偏差程度和杠杆率的影响。h_ii越大分母越小同一个残差对应的学生化残差越大。在 R 中可以直接用rstandard()或rstudent()得到内部学生化残差和外部学生化残差# 内部学生化残差使用当前模型的误差方差估计 r_int - rstandard(lm.fit1) # 外部学生化残差删除第i个点后重新估计误差方差 r_ext - rstudent(lm.fit1) # 判断哪些点超过阈值 outlier_flag - abs(r_ext) 3 sum(outlier_flag)rstandard()计算时使用包含所有样本的误差方差估计rstudent()则对每个点执行“删除一个样本后再估计方差”的策略对异常点更敏感。经验上外部学生化残差绝对值大于 3 的点值得高度关注。代码中最后一行的sum()统计异常点数量方便批量筛查。2.4 避坑学生化残差与普通残差的三个典型误用现象直接比较普通残差的大小把e_i最大的几个点当作异常点结果剔除后模型反而变得更差某些正常点被误删。原因普通残差方差不齐h_ii较大的点天然残差方差更小同样的偏离程度表现为更大的e_i导致高杠杆点被优先标记为异常点而真正的离群点可能因为杠杆率低被漏掉。解决用rstandard()或rstudent()替代普通残差。学生化残差分母中加入了sqrt(1 - h_ii)修正了方差不等的影响。我在实际项目中基本只用rstudent()它对单个异常点更敏感。现象abs(r_ext) 2标记出大量点把阈值放宽到 2 后异常点比例超过 10%模型被削掉太多样本。原因样本量较大时学生化残差的分布接近t分布在n 50时约 5% 的点可能超过 2但这不代表它们是异常点。阈值设置过松会把正常波动当成异常。解决以abs(r_ext) 3作为首要关注线同时结合 Cook 距离判断强影响性。不要只依据单一指标删点应该综合残差、杠杆率、Cook 距离三维度。现象删除了所有学生化残差超阈值的点后重新拟合发现删点前模型系数还在合理范围删点后某个自变量变得不显著或系数符号反转。原因一个点既是异常点又是强影响点时它对系数的拉动作用可能掩盖了其他点的模式。盲目删除所有异常点可能破坏本来稳定的数据结构。解决先看 Cook 距离优先关注“影响大”的点而不是“偏差大”的点。异常点不一定有强影响高杠杆点也不一定是强影响点需要区分对待。3. Cook 距离与强影响点综合杠杆率和残差的判断标准3.1 Cook 距离公式拆解为什么它同时包含 h_ii 和 r_iCook 距离是回归诊断中使用频率最高的影响度量指标其公式为D_i (r_i² / p) * (h_ii / (1 - h_ii))其中r_i是第i个点的学生化残差p是模型中参数个数含截距h_ii是杠杆率。这个结构很有深意第一项r_i² / p度量残差偏离程度第二项h_ii / (1 - h_ii)是杠杆率的单调变换。两个因子相乘意味着一个点只有同时具备“残差大”和“杠杆高”两个特征时Cook 距离才会显著。单纯残差大但杠杆低或者杠杆高但残差小D_i都不会太大。这与强影响点的定义高度吻合强影响点是指剔除后对回归系数估计有显著效应的观测值。在 R 中的计算方式非常直接# 使用基本包的 cooks.distance 函数 d1 - cooks.distance(ols) # 查看 Cook 距离最大的样本 which.max(d1) # 结合学生化残差和杠杆率构成诊断矩阵 r - stdres(ols) h - hatvalues(ols) # 输出高杠杆、高残差、高 Cook 距离的样本 diag_matrix - data.frame( id 1:nrow(cdata), cook_d round(d1, 4), std_resid round(r, 3), leverage round(h, 4) ) head(diag_matrix[order(-diag_matrix$cook_d), ], 10)代码中cooks.distance()返回每个样本的 Cook 距离stdres()提取标准化残差hatvalues()提取杠杆率。构建的数据框把三个核心诊断量并列展示按 Cook 距离降序排列后可以直观看到哪些点对模型影响最大。这里的ols是之前lm(crime ~ poverty single, data cdata)的拟合结果在 UCLA 的crime.dta数据集上运行分析crime与poverty、single两个自变量的关系。3.2 经验分界点 4/n 的由来与争议Cook 距离的判断阈值在学术界一直存在争议。最常用的经验分界点是4/n其中n是样本量。在 R 中筛选强影响点的标准写法是# 按 4/n 阈值筛选强影响点 n - nrow(cdata) influential - cdata[d1 4 / n, ] influential # 同时也可以参考 F 分布的分位数 qf_threshold - qf(0.5, df1 2, df2 n - 2) influential_f - cdata[d1 qf_threshold, ]4/n是经验法则来源于 Cook 距离与 F 分布近似关系中取F(0.5, p, n-p)的近似结果。另一种做法是用qf(0.5, p, n-p)直接计算 F 分布 50% 分位数作为阈值这在p2时通常比4/n略宽松。实际问题中我一般两种都跑一遍把落在两个阈值之间但又不算极端的样本标记为“重点关注”。3.3 实际分析crime 数据集中第 9、25、51 号样本的处理在 crime 数据集上运行plot(ols, las 1)会生成四张诊断图。从残差图和 Cook 距离图可以清晰看到第 9、25、51 号观测值位于边缘位置。进一步用数值确认# 查看这3个样本的具体诊断值 target_ids - c(9, 25, 51) diag_matrix[target_ids, ] # 输出这些样本的原始数据 cdata[target_ids, ]输出的诊断矩阵显示这三个点的 Cook 距离都超过了4/51的阈值标准化残差绝对值也偏高。此时面临一个典型决策场景如果直接采用 OLS你可能会倾向于删除这三行数据再重新拟合但如果删除后模型系数变化巨大说明这些点是强影响点但未必是“错误数据”。稳健回归提供了第三条路不剔除样本而是降低它们的权重。3.4 避坑Cook 距离阈值的两个常见翻车现场现象使用4/n阈值筛出 5 个强影响点全部删除后重新拟合发现某个自变量系数符号反向拟合优度下降。原因强影响点不一定都是“坏点”。如果这个点代表了真实存在的特殊子群体比如高收入低犯罪率的城市删除它会让模型丧失对这类群体的解释能力。4/n是经验阈值样本量小或自变量维度高时容易误判。解决先记录强影响点对应的实际业务含义再决定是否删除。通常我会保留这些样本改用稳健回归或加权回归让数据自己决定权重。现象plot(lm.fit1)四张图中 Cook 距离图看起来没有超过红虚线但手工计算cooks.distance()却发现值超过4/n两套结果不一致。原因plot()函数绘制的 Cook 距离图纵轴范围可能被自动缩放红虚线是 R 根据 Cook 距离分布计算的可视化阈值而不是严格的4/n边界。两种呈现逻辑不同导致肉眼判断与数值判断冲突。解决以cooks.distance()的数值结果为准plot()图只作为初步筛查。数值筛选后用identify()或which()定位具体样本ID再回到业务层面判断。4. rlm() 实现稳健回归Huber 与 Bisquare 两种 M 估计的完整实战4.1 为什么选择 rlm最小二乘在异常点面前的两个困境最小二乘估计的目标是使残差平方和最小这意味着一个大残差点会以平方级别拉动回归线。面对异常点和高杠杆点时OLS 有两个困境第一如果异常点来自数据录入错误理论上应该剔除但数据分析者很难有充分证据证明“这个点一定是错的”第二如果异常点来自另一个总体或特殊子群体直接删除会造成样本选择偏差。稳健回归的思路是在“完全剔除”与“一视同仁”之间折中对残差较大的观测值赋予较低权重对正常样本保留高权重。rlm()是 MASS 包中的核心函数实现了 M 估计的迭代重复加权最小二乘算法。其基本流程是先用 OLS 得到初始残差根据残差大小计算观测权重再用加权最小二乘更新系数然后重新计算残差和权重迭代直到收敛。权重函数的选择决定了稳健性的具体形式。4.2 Huber 方法的权重函数与参数选择Huber 方法的权重函数是分段函数w(e) 1当|e| cw(e) c / |e|当|e| c其中c是截断常数R 中默认取1.345。这意味着残差在阈值内的观测获得权重 1残差超过阈值的观测权重随残差增大而递减。Huber 估计对中等程度的异常值表现稳健同时保留了较高的统计效率。在 R 中的用法# 加载 MASS 包 library(MASS) # Huber 方法的 M 估计 rr.huber - rlm(crime ~ poverty single, data cdata) # 查看模型摘要 summary(rr.huber) # 查看每个观测的最终权重 weights_huber - rr.huber$w head(sort(weights_huber, decreasing FALSE), 10)summary(rr.huber)输出与lm()类似包含系数估计和t值但注意这里不展示 F 统计量和 R²因为迭代加权过程让这些统计量的解释变得复杂。rr.huber$w保存了每个观测的最终权重权重最小的点就是被降权最厉害的点。4.3 Bisquare 方法的权重函数与参数选择Bisquare也常称为 Tukeys biweight方法的权重函数是w(e) (1 - (e/c)²)²当|e| cw(e) 0当|e| c与 Huber 方法不同Bisquare 给所有非零残差的观测都赋予递减权重残差超过c的观测权重直接归零。R 中默认c 4.685。这意味着 Bisquare 比 Huber 更“激进”它可以完全剔除极端异常点的影响而 Huber 对极端残差仍然保留c/|e|的微小权重。# Bisquare 方法的 M 估计 rr.bisq - rlm(crime ~ poverty single, data cdata, method MM) # 或者显式指定 psi 函数为 bisquare rr.bisq2 - rlm(crime ~ poverty single, data cdata, psi psi.bisquare) # 查看权重分布 summary(rr.bisq$w)代码中method MM表示使用 MM 估计它结合了高分解值和高效率特性是处理强影响点时的推荐选择。psi psi.bisquare显式指定所用的psi函数MASS 包中内置了psi.huber和psi.bisquare。MM 估计在初始化阶段使用高分解值的估计方法然后进入 Bisquare 迭代比默认的 M 估计更稳健。4.4 权重结果对比同一批样本在两种方法下的待遇差异将两种方法的权重提取出来对比是理解稳健回归最直观的方式# 合并两种权重进行对比 weight_compare - data.frame( id 1:nrow(cdata), huber_w round(rr.huber$w, 4), bisq_w round(rr.bisq$w, 4), std_resid_ols round(stdres(ols), 3) ) # 查看权重最低的10个样本 head(weight_compare[order(weight_compare$huber_w), ], 10) # 计算两种权重与 OLS 标准化残差的相关性 cor(weight_compare$huber_w, abs(weight_compare$std_resid_ols)) cor(weight_compare$bisq_w, abs(weight_compare$std_resid_ols))通常你会发现Huber 方法中权重最小的点对应原始 OLS 标准化残差最大的点但权重不会降到 0Bisquare 方法则可能将极端残差点权重直接置零。两个模型的系数估计差异反映了稳健回归的“折中”程度。Huber 适合你怀疑异常点有少量信息但不愿完全放弃的场景Bisquare 适合你认为部分点真的来自其他总体的场景。4.5 避坑rlm() 使用中的四个高频报错与处理现象rlm()运行后提示convergence相关警告或者迭代次数未达到默认上限就停止结果似乎仍未稳定。原因M 估计的迭代是从 OLS 初始值开始的如果初始模型中有极端强影响点权重函数可能在某些点产生周期性振荡迭代难以收敛。默认最大迭代次数可能不足。解决增加迭代次数或调整初始值。可以传入maxit 100参数也可以先利用lm()拟合后剔除极端 Cook 距离点再用剩余样本的系数作为初值。现象rlm(crime ~ poverty single, data cdata)报错提示variable lengths differ或者NA/NaN/Inf in foreign function call。原因数据中存在缺失值。rlm()默认使用na.omit处理缺失值但部分情况下数据框中的NA会在权重计算中引发错误。解决拟合前手动执行cdata - na.omit(cdata)同时检查是否存在Inf值。如果某个自变量的分布严重偏态考虑先做对数变换再进入模型。现象拟合成功但summary(rr.huber)输出的系数与lm()差别不大怀疑稳健回归没有起作用。原因数据集中本身没有严重的异常点或高杠杆点稳健回归和 OLS 自然结果接近。这不是 bug而是正常现象。稳健回归的价值在数据“脏”的时候才体现。解决在拟合前先画出散点图或执行诊断矩阵确认数据中确实存在候选异常点。如果诊断结果表明数据干净直接报 OLS 结果即可。现象Bisquare 方法拟合后大量观测权重为 0模型的有效样本量大幅下降标准误增大。原因psi.bisquare的默认截断常数c 4.685对应的残差阈值是在正态误差假设下确定的如果数据中存在多个相互靠近的异常点遮蔽效应可能导致过多样本被降权。解决改用method MM提高分解值或者适当调大c值比如psi psi.bisquare, c 5.5。但注意调大c会降低稳健性需要权衡。5. 完整 R 代码实战从 OLS 诊断到稳健回归的参数对比5.1 数据读取与模型拟合的完整流程结合前文提到的crime.dta数据集完整流程从读取外文格式数据开始。R 中读取 Stata 格式数据需要使用foreign包# 加载所需包 require(foreign) require(MASS) # 读取 Stata 格式数据 cdata - read.dta(https://stats.idre.ucla.edu/stat/data/crime.dta) # 查看数据结构 str(cdata) names(cdata) # 拟合普通最小二乘回归 ols - lm(crime ~ poverty single, data cdata) # 输出模型摘要 summary(ols)read.dta()是读取 Stata 数据文件的标准函数其网络路径直接加载数据。str(cdata)查看各变量的类型和取值分布确保crime、poverty、single都是数值型。summary(ols)输出的系数表中需要重点关注poverty和single的估计值及显著性。5.2 四图诊断与数值诊断的配合诊断不能只依赖plot()生成的图形还需要数值输出来确定具体样本编号。完整流程如下# 四图诊断 opar - par(mfrow c(2, 2), oma c(0, 0, 1.1, 0)) plot(ols, las 1) # 计算 Cook 距离和标准化残差 d1 - cooks.distance(ols) r - stdres(ols) h - hatvalues(ols) # 构建诊断矩阵 a - cbind(cdata, d1, r, h) # 按 4/n 阈值筛选 n - nrow(cdata) a[d1 4 / n, ]代码中par(mfrow c(2, 2))将图形区域分割成 2x2 的网格四张诊断图依次排列。cbind()将原始数据与三个诊断量合并成新数据框方便筛选和查看。a[d1 4 / n, ]筛选出 Cook 距离超阈值的全部样本输出包括原始变量和诊断量可以直接对照样本 ID 查看业务含义。5.3 稳健回归与 OLS 系数对比表# OLS 系数 coef_ols - coef(ols) # Huber 稳健回归系数 coef_huber - coef(rr.huber) # Bisquare 稳健回归系数 coef_bisq - coef(rr.bisq) # 合并结果生成对比表 compare_table - data.frame( OLS round(coef_ols, 4), Huber round(coef_huber, 4), Bisquare round(coef_bisq, 4) ) print(compare_table)对比表的价值在于直观展示三种方法对同一批数据的系数估计差异。如果 Huber 和 Bisquare 的系数与 OLS 明显不同说明异常点对 OLS 的拉动效应已经被稳健回归修正如果三者结果接近说明数据本身质量较好。同时可以对比标准误# 对比标准误 se_ols - summary(ols)$coefficients[, 2] se_huber - summary(rr.huber)$coefficients[, 2] se_bisq - summary(rr.bisq)$coefficients[, 2] cbind(OLS_se se_ols, Huber_se se_huber, Bisquare_se se_bisq)5.4 参数选择建议不同场景下的 c 值与 method 设置rlm()的参数选择需要结合数据特征和业务需求。以下是我常用的参数设置参考表数据特征methodpsi 函数c 值理由基本干净偶发小异常Mpsi.huber1.345保留效率只修正重尾存在若干个孤立异常值Mpsi.bisquare4.685对极端残差直接归零异常点较多或聚集成簇MMpsi.bisquare4.685高分解值抗遮蔽效应高杠杆点与异常并存MMpsi.huber3.0杠杆点需要更渐进地降权大样本追求效率Mpsi.huber1.5放宽阈值减少有效样本损失这个表的核心逻辑是异常点越多、越极端越倾向于使用分解值更高的估计方法和更激进的权重函数。method MM比默认的M估计多一个高分解值初始化步骤能有效抵抗多个异常点相互遮蔽的情况。5.5 避坑稳健回归结果解读中的三个常见错误现象用summary(rr.huber)中的 R² 与 OLS 的 R² 比较认为稳健回归拟合效果“更好”或“更差”。原因rlm()的输出并不包含与传统 OLS 直接可比的 R²。迭代加权过程中使用的权重改变了目标函数R² 不再具有“解释方差比例”的标准含义。解决比较模型时使用系数大小、标准误、残差的稳健性和预测效果不要用 R² 作为主要判据。现象把 Huber 和 Bisquare 的权重当作样本质量的绝对评分权重低的样本被认为“一定有问题”。原因权重反映的是“在当前模型设定下这个样本对回归拟合的影响相对较小”不直接等同于“这个样本是错误的”。一个在业务上重要但偏离主趋势的样本权重可能被压低但它仍然包含真实信息。解决将低权重样本单独输出到业务层面验证是否符合预期。若符合业务逻辑应保留在数据集中甚至可以考虑单独建模。现象直接引用rr.huber$w中的权重进行二次加权分析没有意识到权重是在拟合后固定的。原因rlm()的权重是迭代收敛后的产物它们依赖于最终系数估计。换个模型设定权重会完全改变不能当作外生变量使用。解决除非在做敏感性分析否则不要在后续分析中直接使用rlm()的权重作为通用样本权重。如果需要稳定的加权方案应基于领域知识预先定义权重。6. 杠杆率、Cook 距离与权重的联动验证一个手工计算技巧验证稳健回归是否“做对了事”有一个很实用的技巧把手动计算的杠杆率、Cook 距离与rlm()输出的权重放到同一个数据框里用相关性检验判断降权是否准确瞄准了最需要降权的样本。具体做法是计算每个样本的 Cook 距离或者杠杆率与其在稳健回归中权重的 Spearman 相关如果降权逻辑正确高 Cook 距离的样本应该获得低权重。这样做的价值在于它用数据验证了“权重函数是否真的在折中处理极端点”而不是只看系数差异。具体验证代码如下# 计算三个诊断量 h - hatvalues(ols) d1 - cooks.distance(ols) r - abs(stdres(ols)) # 提取两种稳健回归的权重 w_huber - rr.huber$w w_bisq - rr.bisq$w # 构建验证数据框 verify_df - data.frame( leverage h, cook_d d1, abs_stdres r, w_huber w_huber, w_bisq w_bisq ) # 计算 Spearman 相关系数 cor_leverage_huber - cor(verify_df$cook_d, verify_df$w_huber, method spearman) cor_leverage_bisq - cor(verify_df$cook_d, verify_df$w_bisq, method spearman) # 输出相关系数 cat(Cook距离与Huber权重的Spearman相关:, cor_leverage_huber, \n) cat(Cook距离与Bisque权重的Spearman相关:, cor_leverage_bisq, \n) # 找出权重最低但 Cook 距离不高的样本检查是否有异常降权 low_w_but_low_cook - verify_df[ verify_df$w_huber quantile(verify_df$w_huber, 0.1) verify_df$cook_d quantile(verify_df$cook_d, 0.5), ] print(low_w_but_low_cook)这种方法在 Huber 下通常表现出高度负相关因为 Huber 的权重直接由残差大小决定而 Cook 距离的主要驱动因子恰恰是学生化残差但在 Bisquare 下由于权重函数在阈值处截断相关可能变弱。这解释了为什么 Bisquare 对极端点的处理更彻底对中间型异常点的降权却可能更温和。验证完成后把注意力放回业务层面。我通常会在输出结果时保留三样东西OLS 残差的散点图、稳健回归权重的分布直方图、以及按权重排序的前十个样本的业务标签。这三样配合能有效回答“为什么某个样本被降权”以及“这个降权是否合理”。这是我自己比较习惯的一种做法。从那以后我每次做稳健回归都会强制走一遍这个流程先用plot()和cooks.distance()做诊断确认异常点和高杠杆点的位置再用rlm()配合 Huber 或 Bisquare 权重跑一遍最后用 Spearman 相关验证降权逻辑是否与诊断结论一致。只有这三步全部完成我才敢把模型结果写进分析报告。希望帮到你。本文还有配套的精品资源点击获取