企业官网建设流程全解析

简介：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̄)' (X'X)^{-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(X'X)^{-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_i

Cook 距离是回归诊断中使用频率最高的影响度量指标，其公式为：

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% 分位数作为阈值，这在p=2时通常比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（也常称为 Tukey's 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 标准化残差最大的点，但权重不会降到 0；Bisquare 方法则可能将极端残差点权重直接置零。两个模型的系数估计差异反映了稳健回归的“折中”程度。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()的参数选择需要结合数据特征和业务需求。以下是我常用的参数设置参考表：

这个表的核心逻辑是：异常点越多、越极端，越倾向于使用分解值更高的估计方法和更激进的权重函数。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 相关验证降权逻辑是否与诊断结论一致。只有这三步全部完成，我才敢把模型结果写进分析报告。希望帮到你。

本文还有配套的精品资源，点击获取