简介一份聚焦R语言稳健性估计的实例分析演示文稿面向统计学、数据分析及建模学习者。内容围绕线性回归模型稳健性问题系统讲解异常点、高杠杆点与强影响点的定义和区别并结合R代码展示lm()拟合、plot()四幅诊断图、resid()残差提取、cooks.distance()距离计算等操作同时给出学生化残差绝对值大于3、Cook距离大于4/(n-p-1)等识别准则。PPT中还包括杠杆率、学生化残差、Cook距离等关键公式的推导以及Huber和Bisquare两种M估计方法的原理介绍说明稳健回归本质上是加权最小二乘通过迭代赋权降低异常值影响帮助读者在数据存在异常值时选用稳健回归替代最小二乘法提高预测准确性与模型泛化能力。资源共1个pptx文件压缩包大小716KB按PPT页顺序组织适合课堂展示或自学查阅。当前已有1129人学习下载内容精炼且可直接用于教学演示。1. 当异常值把回归线拉弯稳健性估计到底解决什么问题做回归分析时最怕的不是模型不显著而是数据里藏着几个异常值把一条本该平直的拟合线硬生生拉弯。稳健性估计要解决的就是这个问题——在不提前删除任何观测的前提下通过加权手段把异常点和高杠杆点的影响压下去让回归系数更贴近多数数据的真实趋势。这套 R 语言实例包含完整的回归诊断流程残差图、杠杆率、Cook 距离和稳健回归代码Huber 与 Bisquare 两种权重函数适合经常和数据打架的统计建模从业者也适合刚学会lm()却不知道怎么处理异常值的 R 初学者。下面我会从回归诊断讲起逐步落到rlm()的实现细节和踩坑记录。2. 回归诊断三件套残差、高杠杆与强影响点的判定界线2.1 普通残差为什么不能直接当异常点判据普通残差的定义很直白e_i Y_i - Ŷ_i也就是观测值与模型预测值的差。在 R 里用resid(lm.fit1)就能取出来。但这里有个隐蔽的坑直接拿e_i的大小比异常点是不可靠的。原因是残差的方差不是常数它等于Var(e_i) σ²(1 - h_ii)其中h_ii是帽子矩阵对角线上的第 i 个元素。也就是说位于自变量空间边缘的点残差方差反而更小。两个绝对值相同的残差一个出现在数据中心一个出现在数据边缘它们对回归诊断的意义完全不同。正文中那句“直接比较 e_i 是不适宜的”说的就是这件事。学生化残差就是把这个差异修正掉公式为r_i e_i / (s * sqrt(1 - h_ii))其中s是剩余标准差h_ii是杠杆率。这个统计量把残差、杠杆率、数据尺度三者揉在一起是判断异常点更合理的基准。R 里可以用rstandard()或stdres()拿到它后一个来自 MASS 包。2.2 帽子矩阵与杠杆率离中心越远杠杆越大杠杆率的概念来自帽子矩阵H X(XX)^{-1}X。为什么叫帽子矩阵因为Ŷ HY这个矩阵给 Y 戴了一顶“帽子”。对角线的第 i 个元素h_ii就是第 i 个观测的杠杆率正文给出了它的展开形式h_ii 1/n (x_i - x̄)(XX)^{-1}(x_i - x̄)这个式子的含义比帽子矩阵本身更直观第一项1/n是截距项带来的基础杠杆第二项刻画了x_i偏离样本中心x̄的距离。离中心越远h_ii越大这个点就越可能把回归线拉向自己。全体观测的杠杆率之和等于模型参数个数 p自变量数加截距项所以平均杠杆是p/n。经验上h_ii超过2p/n或3p/n就需要留意但这不是硬性阈值更多是用来看排序。需要注意高杠杆点并不等于坏点。如果它的 Y 值恰好落在由其他多数点决定的回归趋势上它对系数估计的影响很小。高杠杆只是说明它在 X 空间里离群是否真的有影响要看残差是否也大。2.3 Cook 距离把杠杆和残差焊在一起Cook 距离的综合性强得多它的公式是D_i [h_ii / (p(1 - h_ii))] * r_i²看起来复杂其实结构很清楚前半部分是杠杆项后半部分是学生化残差平方。这两项一乘高杠杆且高残差的点Cook 距离会特别突出而单纯杠杆高但残差小或者残差大但杠杆低Cook 距离都不容易超标。关于判定阈值正文里提到了两个标准。一个是用D_i 4/n作为经验分界其中 n 是观测个数这是 UCLA 教程和很多回归诊断教材里推荐的做法灵敏度较高另一个是教科书里常见的D_i 1这个阈值非常保守通常只有极端强影响点才能超过。实际项目中我一般先按4/n筛一遍再人工检查 Top 10 的观测是否在业务上说得通而不是机械地看有没有超线。2.4 异常点、高杠杆点、强影响点之间的纠缠关系这三类点在概念上相关但不等价。异常点是从 Y 空间判别的看的是残差高杠杆点是从 X 空间判别的看的是杠杆率强影响点是从结果判别的看的是剔除后系数变化幅度。正文里用 B 点和 C 点说明得很清楚B 点残差大但 X 值在中心区域杠杆低它是个异常点但对回归系数影响有限C 点 X 值离中心远残差又大杠杆和残差双重突出才是真正能改变回归方程的强影响点。实际诊断中的正确顺序是先看残差找出异常点候选再看杠杆率找出高杠杆候选最后用 Cook 距离把两者交汇的点捞出来。捞出来之后不要急着下结论还要结合业务背景确认这些观测是数据录入错误还是真实存在的边缘情况。这决定了下一步是修正还是做稳健回归。3. 用 lm() 做回归诊断四张图与 Cook 距离的实操流程3.1 载入数据与 OLS 基线正文里给的示例数据是 UCLA 教学常用的crime.dta包含犯罪率、贫困率、单亲家庭比例等变量共 51 个观测美国 50 个州加华盛顿特区。读入用foreign包的read.dta()如果需要联网可以先下载到本地再读取。library(foreign) library(MASS) # 如果网络不通可先下载 crime.dta 到本地再运行这一行 cdata - read.dta(https://stats.idre.ucla.edu/stat/data/crime.dta) str(cdata) summary(cdata)read.dta()是读 Stata 格式数据的老牌函数使用前必须加载foreign包。MASS包是给后面的rlm()和stdres()做准备。str()和summary()用来检查数据结构和缺失情况跑任何回归前都建议先过一遍这两步。然后拟合普通最小二乘基线模型ols - lm(crime ~ poverty single, data cdata) summary(ols)这里的crime是因变量poverty和single是两个自变量。lm()做 OLS 估计summary()输出系数、标准误、t 值和 F 检验。保留这份摘要里的残差标准误和 R²后面跟rlm()对比时用得上。3.2 plot(lm) 的四张图怎么读opar - par(mfrow c(2, 2), oma c(0, 0, 1.1, 0)) plot(ols, las 1) par(opar)plot()函数会对lm对象一次性生成四幅诊断图mfrow c(2, 2)把它们排成 2×2 网格oma留出顶部空间。las 1让坐标轴文字横排避免图例标签挤在一起。四幅图分别对应残差对拟合值、残差的正态 Q-Q 图、标准化残差绝对值平方根对拟合值、Cook 距离图。第一幅图重点看残差是否围绕 0 随机分布出现喇叭形说明有异方差出现曲线说明可能有遗漏的非线性项。第二幅图看残差是否近似正态样本点偏离直线越远正态性假设越可疑。第三幅图进一步检查方差是否稳定。第四幅图直接按观测序号画出 Cook 距离超过阈值线的点就是强影响点候选。在这份crime数据上第 9、25、51 号观测在图中明显突出。这三条正是正文里点名需要重点检查的候选项。3.3 Cook 距离与标准化残差联合筛选诊断图用来做视觉检查定量筛选还需要直接算统计量d1 - cooks.distance(ols) r - stdres(ols) a - cbind(cdata, d1, r) a[d1 4/51, ]这段代码做了三件事。第一行cooks.distance()计算每个观测的 Cook 距离第二行stdres()计算学生化残差注意这个函数来自 MASS 包加载顺序不能错第三行用cbind()把原始数据、Cook 距离、学生化残差拼在一起最后按d1 4/51筛选。阈值4/51就是正文里说的经验分界n 为 51 个观测4/n ≈ 0.0784。这个阈值筛出的观测基本覆盖了图上突出的强影响点。3.4 筛出强影响点后先别动手删很多初学者走到这一步看到 Cook 距离超线的观测第一反应是“删掉重跑”。这个操作在业务上非常危险——你删掉的可能是真实存在但边缘的客户群体或地区样本。正文的思路是先用 OLS 的残差、拟合值、Cook 距离和杠杆率做完完整诊断再转向稳健回归而不是直接删点。正确的做法是先把强影响点列出来逐一核对原始数据是否录错。如果确认是真实观测就应该用稳健回归作为最终模型。所谓稳健回归本质是在“完全保留异常值”和“完全删除异常值”之间做一个加权折中下一章详细展开。4. 用 rlm() 做稳健回归Huber 与 Bisquare 的实现与参数4.1 M 估计为什么能抗异常值OLS 的目标是最小化残差平方和即Σe_i²。这个目标函数的问题在于残差大一点平方后就成了主导项。一个离群点的残差可能是正常点残差的 10 倍平方后贡献就是 100 倍回归系数被它牵着走。M 估计把目标函数换成Σρ(e_i / s)其中ρ是损失函数它对大残差的惩罚不再按平方增长。Huber 损失函数在残差较小时保留二次形式残差超过阈值后转为线性增长这样大残差虽然仍有影响但不会主导整个目标函数。Bisquare 更极端超过阈值的残差其观测的权重直接归零相当于完全排除。因为权重的计算依赖残差残差又依赖权重所以 M 估计不能一步到位必须用迭代重复加权最小二乘IRLS求解。R 里的rlm()内部就是这么做的先跑一遍 OLS 得到初始残差计算权重再用加权最小二乘更新系数算新残差再算新权重循环直到收敛。4.2 Huber 和 Bisquare 权重函数的核心差异两种权重函数的核心参数差异可以放在一起对比。Huber 的调节参数是k 1.345当标准化残差|u| ≤ k时权重为 1当|u| k时权重按k/|u|递减。Huber 权重永远不会归零无论残差多大观测点都保留着部分影响力。Bisquare 的调节参数是c 4.685标准化的残差绝对值在c以内时权重按[1 - (u/c)²]²递减一旦超过c权重直接清零。这意味着 Bisquare 会把极端异常点彻底排除在模型之外。项目HuberBisquare常用参数k 1.345c 4.685残差较小时权重 1权重递减残差较大时权重 k/u极端点的地位保留部分影响完全排除实际分析中Huber 更稳健受异常值影响小但对极端异常值仍有敏感度Bisquare 对极端值更无情但对初始值和尺度估计更敏感。如果数据里异常值特别极端建议先跑 Huber 看权重分布再跑 Bisquare 对比系数。4.3 rlm() 代码与参数解读加载完MASS包后rlm()可以直接用公式接口拟合# Huber 方法rlm 默认 rr.huber - rlm(crime ~ poverty single, data cdata) summary(rr.huber) # Bisquare 方法显式指定权重函数和最大迭代次数 rr.bisquare - rlm(crime ~ poverty single, data cdata, psi psi.bisquare, maxit 100) summary(rr.bisquare)第一条rlm()使用默认的 Huber 权重不需要额外参数。第二条将psi参数设为psi.bisquare切换为 Bisquare 方法。maxit 100用来控制迭代上限实际中如果数据特别脏默认的 20 次迭代可能不够。看summary()输出时注意两个细节。第一是系数对比 OLS 的结果稳健回归的系数通常有变化尤其是poverty和single的估计值可能向 0 缩一些也可能反过来变大具体方向取决于强影响点从哪个方向拉动回归线。第二是标准误稳健回归的标准误通常比 OLS 大这是合理的因为它没有假装那些异常值不存在而是保留了数据本身的不确定性。4.4 从权重向量看数据rlm()对象里保存了最后一次迭代的权重可以直接用$w取出来w.huber - rr.huber$w w.bisquare - rr.bisquare$w # 看一下各自权重最小的 3 个观测 head(sort(w.huber), 3) head(sort(w.bisquare), 3)Huber 的权重大多为 1只有少数观测的权重低于 1且不会出现 0。Bisquare 的权重是连续递减的极端点直接是 0。把这两个向量和原始数据的poverty、single列对照能看到被压权的观测在自变量空间中确实处于边缘位置。这里有一个正文提到的关键点在 Bisquare 方法下所有非零残差对应观测的权重都是递减的所以 Bisquare 对异常值的惩罚更均匀Huber 对“中间偏大”的残差则采取了一种更宽容的姿态。选哪种没有绝对答案我的习惯是两种都跑如果系数差异不大就选解释上更简单的那一个如果差异大说明异常值确实在影响结论需要进一步核查数据源。5. 稳健回归避坑指南五个常见问题与排查方法5.1 rlm() 报“不收敛”或迭代超限现象rlm()运行后提示类似 “failed to converge in 20 iterations”。原因这是最常见的问题。rlm()默认最大迭代次数是 20 次当数据中同时存在多个强影响点时权重和残差互相迭代来回振荡20 次内达不到收敛容差。并不是模型坏了而是迭代不够。解决在rlm()里显式调大迭代上限比如maxit 100。如果调到 100 仍不收敛需要检查初始值是否被异常值污染得太过分可以先手动剔除 Cook 距离最大的 1 个点再跑一次看权重是否稳定。5.2 Huber 权重几乎全是 1稳健回归结果和 OLS 一模一样现象rr.huber$w里绝大多数权重都等于 1系数和lm()几乎没有差别。原因Huber 的阈值k 1.345对应标准正态分布的约 81% 分位点也就是说只有当标准化残差超过 1.345 时权重才会开始减少。如果异常值虽然大但标准化后没有越过这个门槛Huber 就识别不出来。这种情况并不罕见因为学生化残差的分母里包含了可能的异常值贡献方差一大阈值就被稀释了。解决换 Bisquare 权重函数跑一遍或者把 Huber 参数k调小。rlm()支持psi psi.huber, k 1这样的写法。但调参要谨慎阈值越小正常点被误伤的概率越高。我的做法是先用默认 Huber 跑通流程再用 Bisquare 看权重分布最后回头判断哪些点真正需要压权。5.3 Cook 距离阈值之争4/n、4/(n-p-1)、还是 1现象同一份数据用D_i 4/n筛出 9 个点用D_i 1只筛出 1 个点。文献里还能找到4/(n-p-1)这个版本三套标准结果差距很大。原因Cook 距离没有精确的抽样分布所有阈值都是近似经验法则。D_i 1非常保守适合小样本4/n灵敏度高样本越大越严格4/(n-p-1)是 F 分布近似下的简化版性质和4/n类似但更贴近统计理论。解决不要纠结选哪个阈值把 Cook 距离按从大到小排个序看前 5 到前 10 名是谁。阈值只是筛子业务判断才是最终依据。如果某个观测排第二但没超过4/n它同样值得调查。正文中也明确提到“此处的判断标准有争议”所以把 Cook 距离当成排序工具而不是判定工具是更稳的用法。5.4 高杠杆点被当成强影响点误伤现象hatvalues(ols)里某个观测的杠杆率明显超过2p/n看起来像个“问题点”但 Cook 距离并不高删除前后回归系数变化很小。原因这个观测可能既远离数据中心又恰好落在其他多数点决定的趋势线上。也就是说它在 X 空间里离群但在 Y 方向上没有偏离——回归线本来就在它脚底下。正文里 B 点就是这么个情况从残差角度看是异常点从影响角度看不是。解决先画散点图或半正态概率图把 X 空间边缘点标注出来再结合残差看。只有当杠杆率高且学生化残差也大时才是真正值得警惕的强影响点。单纯杠杆高但不影响系数不需要做任何处理。5.5 手动删除异常点后再回归后果比想象严重现象为了“让模型干净”先删掉几个 Cook 距离超线的观测再跑lm()结果 R² 上升、p 值变漂亮看起来数据被“修正”了。原因手动删点相当于把模型的不确定性“修”掉了。你是在看完全部残差和诊断图之后才决定删除哪些点这本身已经是利用结果做决策再用这份数据做统计推断标准误会被低估置信区间和 p 值都会失真。解决能不动手就不动手。诊断发现异常值后先核实数据录入是否有误如果数据正确直接用rlm()做稳健回归。稳健回归把“保留”和“删除”做成了权重折中异常值不会被完全排除但也不至于拉着系数跑偏。汇报结果时把权重最低的几个观测列进附录说明它们对模型的潜在影响这才是更诚实的做法。6. 进阶验证把权重和残差摊开对比 OLS 与稳健回归到底差在哪6.1 权重可视化看出哪些点被压权了rlm()跑完后我习惯把权重画出来看一眼这比任何统计量都直观par(mfrow c(1, 2)) plot(cdata$poverty, rr.bisquare$w, xlab poverty, ylab Bisquare weight, main Bisquare weights by poverty) plot(cdata$single, rr.bisquare$w, xlab single, ylab Bisquare weight, main Bisquare weights by single)左图看权重随poverty的变化右图看随single的变化。如果某个区域出现明显的低权重点说明这些观测在对应自变量维度上处于边缘位置且残差大正是它们在拉动 OLS 的系数。6.2 系数对比表稳健回归改变了什么用一行代码把三套系数放在一起coef_compare - rbind( OLS coef(ols), Huber coef(rr.huber), Bisquare coef(rr.bisquare) ) round(coef_compare, 4)如果三行系数差别不大说明这份数据的异常值虽然存在但影响有限如果某个变量的系数在 OLS 和稳健回归之间明显移动说明强影响点确实在单方向拉动估计值。这个对比结果值得写进分析报告比单独列某个统计量有说服力得多。6.3 残差分布对比稳健回归是否改善了拟合boxplot(resid(ols), resid(rr.huber), resid(rr.bisquare), names c(OLS, Huber, Bisquare), ylab Residuals)稳健回归的残差分布通常更紧凑——因为极端残差对应的权重被压低了但它们并没有从数据集里消失。箱线图会清楚显示这一点。如果 Huber 和 Bisquare 的残差分布差别很大说明存在极端异常值Bisquare 把它们的权重清零后残差整体大幅收窄如果差别不大说明数据中的异常值不算极端。最后把我自己的习惯分享给你从那以后我每次跑完summary(lm())之后的下一步操作永远是cooks.distance()加rlm()哪怕最终结论没有变化也要走一遍才放心。数据里有没有藏着影响结论的点不到权重摊开那一刻谁也不知道。如果你手头的回归模型刚好有几个不省心的观测不妨把这份资源里的代码复制出来跑一遍让稳健回归替你回答“该不该删、怎么处理权重”。希望帮到你。本文还有配套的精品资源点击获取