企业官网建设流程全解析

简介：在回归分析中，离群值往往会让最小二乘估计产生明显偏差，导致模型结论失真。稳健回归通过替换损失函数或迭代加权策略，有效削弱异常点对系数的影响，是数据清洗之外的可靠选择。R语言中MASS和robustbase包提供了rlm与lmrob等函数，分别实现M估计、S估计和MM估计，兼顾稳健性与统计效率。除了系数估计，结合sandwich包计算稳健标准误，或使用covMcd做多元离群点检测，也能在保留原模型的条件下修正推断结论。本文从稳健估计的核心概念讲起，对比不同估计器的崩溃点与适用场景，并用可复现的R代码演示污染数据下的拟合与诊断，帮助数据分析人员在实际项目中合理选择稳健回归方案，避免因少数极端点而误判变量关系。

1. 稳健性估计：为什么最小二乘会被几个离群值牵着鼻子走

拿到一批数据，跑完lm()发现 R² 很高、系数也显著，但一画残差图就冒冷汗：三五个离群值把回归线拽得明显偏离主体。这时候别急着删点，更别轻信"清理异常值"的直觉操作，你应该考虑的是 R 语言里的稳健性估计。稳健性估计的核心是把平方损失换成对异常值不敏感的损失函数，让少数坏点对系数的影响被压到最低，M 估计、S 估计、MM 估计在MASS和robustbase两个包里都有成熟实现。这篇笔记用可复现的 R 代码走完整条路：什么时候该上稳健方法、三种估计器怎么选、稳健标准误和稳健协方差怎么用，以及实操里最容易翻车的五个坑。适合处理回归任务、做实验数据分析、经常被离群值困扰的 R 用户。

2. 稳健性估计的底牌：M 估计、S 估计与 MM 估计怎么选

稳健性估计不是把最小二乘扔掉，而是换一套损失函数和迭代策略。普通最小二乘对离群值脆弱，根源在平方损失：残差每增大一倍，惩罚变成四倍，一个离谱的观测就能把拟合线拉过去。换损失函数之后，同样一份数据，系数估计可能变化很大，这正是稳健性估计存在的意义。但"稳健"也有程度之分，选错了估计器，效果可能比 OLS 还差，所以先把三个关键概念搞清楚。

2.1 影响函数与崩溃点：先看懂稳健性的两个硬指标

判断一个估计器稳不稳健，统计学家看两个数：影响函数和崩溃点。影响函数描述的是"单个观测对估计结果的影响幅度"。OLS 的影响函数是无界的，残差越大，这个点对系数的拉动越没有上限；而稳健估计器的影响函数有界，极端点上再极端，影响也被锁在一个阈值内。崩溃点描述的是"数据里最多能掺入多少比例的坏点，估计器还能给出有效结果"。OLS 的崩溃点是 1/n，意思是一个离群点就可能把结果彻底带偏；中位数的崩溃点是 0.5，哪怕一半数据被污染，估计依然不崩。

这两个指标直接对应选型：如果你的数据里坏点比例可能在 10% 左右，崩溃点至少要有 0.2 才安全；如果你关心的是某个点是否在悄悄拉动系数，影响函数的有界性比崩溃点更值得看。稳健性估计在 R 里的实现，本质都是在"稳健性"和"效率"之间做交易：越稳健，正常数据下的方差就越大。下面这张表是三个主流估计器的对比，也是后面代码里三个函数的选型依据。

2.2 M 估计、S 估计、MM 估计：各自防哪种离群值

M 估计是最直观的稳健回归思路：保留加权最小二乘的框架，但给大残差降权。Huber 的 psi 函数是单调的，残差超过阈值 k 之后，惩罚从平方变成线性，所以它防得住"响应变量方向的离群值"，也就是 y 异常的点。但如果有某个 x 本身离群，形成高杠杆点，Huber 仍然会被该点主导，因为杠杆来自设计矩阵，残差惩罚管不住它。Bisquare 更进一步，残差超过阈值 c 之后权重直接归零，相当于把极端点彻底排除，对高杠杆和 y 方向离群都有一定抵抗力，代价是迭代容易震荡，对初值和尺度估计更敏感。

S 估计的思路不同，它不去拟合系数，而是直接最小化残差的稳健尺度（比如残差 MAD），天然具备 0.5 的崩溃点，但效率偏低，正常数据下估计方差比 OLS 大不少。MM 估计是两者结合：先用 S 估计得到一个稳健的初始尺度和系数初值，再在这个初值基础上做一步 M 估计，把效率拉回到 95%。这个"先保稳健、再保效率"的两段式设计，让它成为目前最常用的高崩溃点稳健回归方案，也就是robustbase::lmrob的默认行为。选型时我一般按这个逻辑走：怀疑只是 y 方向有少量离群值，用rlm(psi=psi.huber)足够；数据污染比例高、或者 x 里也有异常点，直接上lmrob；只有需要多元协方差或离群点检测时才单独考虑covMcd。

2.3 用一段 R 代码判断数据是否需要稳健性估计

在动手换模型之前，先做两个便宜又快的检查。第一个是看 OLS 残差的 SD 和 MAD 比值：MAD 本身是稳健的尺度估计，如果两者比值明显大于 1.5，说明残差分布有重尾，离群值在拖大标准差。第二个更直接：用残差分位数把两侧各 5% 的极端观测剔除后再拟合一次，看系数变化了多少。

library(MASS) library(robustbase) # 假设 dat 是你的数据框，至少包含 x 和 y 两列 # 这里用后面第 3 节会构造的模拟数据做示范 fit0 <- lm(y ~ x, data = dat) r <- resid(fit0) # 检查 1：OLS 标准差与稳健标准差 MAD 的比值 ratio <- sd(r) / mad(r) cat("SD / MAD =", round(ratio, 2), "\n") # 检查 2：剔除残差两端各 5% 后，系数变化幅度 q <- quantile(r, probs = c(0.05, 0.95)) keep <- which(r >= q[1] & r <= q[2]) fit_trim <- lm(y ~ x, data = dat, subset = keep) cat("OLS 斜率:", coef(fit0)["x"], "\n") cat("剔除极端残差后斜率:", coef(fit_trim)["x"], "\n")

这段代码的逻辑是：mad(r)用中位数绝对偏差估计残差的稳健尺度，不受少数极端残差影响，所以sd(r)/mad(r)越大，说明极端残差对普通标准差的贡献越大。第二次拟合的subset = keep只保留残差在 5% 到 95% 分位数之间的观测，如果斜率变化超过 15% 到 20%，说明那 10% 的数据正在实质性地拉动系数，这时候就值得换稳健估计。注意这只是一种快速体检，不是正式检验，但它能帮你在没跑稳健回归之前就判断方向。

3. 用 rlm 与 lmrob 跑通稳健回归：最小命令与参数调优

理论看完了，现在动手。为了让你能原样复现，先构造一份带污染的仿真数据，然后分别用MASS::rlm和robustbase::lmrob拟合，对比系数、标准误和残差表现。这里的所有代码都假设你已经装好了依赖包：install.packages(c("MASS", "robustbase", "sandwich", "lmtest", "boot"))一次装完。换用你自己的数据时，把dat换成read.csv()读进来的数据框即可，列名对应改一下。

3.1 构造一份带离群值的实验数据：set.seed 之后的完整代码

为了模拟真实场景，生成 200 个样本，x 服从均值 10、标准差 3 的正态分布，y 的生成机制是y = 2 + 1.5x + ε，其中 ε 服从均值为 0、标准差为 2 的正态分布。然后随机挑 20 个样本，在 y 上额外叠加一个标准差 25 的噪声，模拟"部分样本被污染"的情况。

set.seed(2024) n <- 200 x <- rnorm(n, mean = 10, sd = 3) e <- rnorm(n, mean = 0, sd = 2) y <- 2 + 1.5 * x + e # 注入 10% 的离群值，方向随机 n_out <- 20 idx <- sample(1:n, n_out) y[idx] <- y[idx] + rnorm(n_out, mean = 0, sd = 25) dat <- data.frame(x = x, y = y) # 快速查看污染效果 summary(lm(y ~ x, data = dat))

这里set.seed(2024)固定随机数种子，保证你得到的dat和我完全一致，这是所有可复现分析的第一步。sample(1:n, n_out)随机挑出 20 个索引，rnorm(n_out, mean = 0, sd = 25)给这些样本的 y 加上大幅噪声，标准差 25 是正常误差标准差 2 的十几倍，足以把 OLS 拉偏。污染比例 10% 是一个很典型的场景：低于 5% 时 OLS 通常还能扛，超过 15% 后 Huber 也开始吃力，10% 正好能区分出不同估计器的差异。

3.2 MASS::rlm：Huber 与 Bisquare 的最小拟合命令

rlm是最容易上手的稳健回归函数，默认用 Huber 的 psi 函数。下面这段代码同时拟合 OLS、Huber、Bisquare 三个模型，并对比斜率。

fit_ols <- lm(y ~ x, data = dat) fit_huber <- rlm(y ~ x, data = dat, psi = psi.huber, k = 1.345) fit_bisq <- rlm(y ~ x, data = dat, psi = psi.bisquare, c = 4.685) # 汇总三个模型的 x 系数 coefs <- data.frame( method = c("OLS", "rlm_Huber", "rlm_Bisquare"), intercept = c(coef(fit_ols)[1], coef(fit_huber)[1], coef(fit_bisq)[1]), slope = c(coef(fit_ols)[2], coef(fit_huber)[2], coef(fit_bisq)[2]) ) print(coefs) # 对比残差尺度 cat("OLS 残差 SD:", sd(resid(fit_ols)), "\n") cat("Huber 残差 MAD:", mad(resid(fit_huber)), "\n")

真值斜率是 1.5，截距是 2。污染数据下，OLS 的斜率会被那 20 个离群点拽着偏离 1.5 不少，Huber 和 Bisquare 会明显更靠近真值。k = 1.345是 Huber 的经典阈值，对应 95% 渐近效率：再调小会更稳健但方差更大，调大则越来越接近 OLS。c = 4.685是 Bisquare 的经典阈值，含义是残差超过 4.685 倍尺度后权重归零。rlm默认用 MAD 作为尺度估计，迭代次数上限在maxit参数里，默认 20 次一般够用，但如果数据污染复杂，可以调到maxit = 100。需要提醒的是，rlm返回的是加权最小二乘的最终结果，summary(fit_huber)里的标准误是基于迭代权重算的，和真正意义上的稳健标准误不完全是一回事，但作为近似参考足够。

3.3 robustbase::lmrob：MM 估计的默认设置与 tuning.psi

如果污染比例超过 10%，或者你怀疑 x 里有高杠杆点，rlm的 Huber 版本会显得"太软"，这时候直接上lmrob。

fit_mm <- lmrob(y ~ x, data = dat) cat("lmrob 斜率:", coef(fit_mm)["x"], "\n") # 显式调整效率参数 fit_mm_95 <- lmrob(y ~ x, data = dat, tuning.psi = 1.5484) cat("tuning.psi=1.5484 斜率:", coef(fit_mm_95)["x"], "\n") fit_mm_85 <- lmrob(y ~ x, data = dat, tuning.psi = 2.5) cat("tuning.psi=2.5 斜率:", coef(fit_mm_85)["x"], "\n")

lmrob默认使用 Tukey Bisquare 作为 psi 函数，tuning.psi = 1.5484对应 95% 效率，这是默认值，通常不用改。tuning.psi调大（比如 2.5）相当于放松对离群值的惩罚，效率提高但稳健性下降；调小则反过来。lmrob的summary()会给出各系数的稳健标准误、t 值和 p 值，比rlm的输出更完整。另外lmrob自带plot()方法，可以画残差 vs 拟合值、标准化残差 vs 稳健距离等诊断图，跑完建议看一眼plot(fit_mm)，特别是第 4 张图的"稳健距离 vs 残差"，能直接标出哪些点是高杠杆或离群。实际项目中我通常把lmrob的结果当成基准答案，再用rlm或 OLS 做敏感性对比，三个结果放在一起，比单独看任何一个都更能说明数据质量。

4. 稳健标准误与稳健协方差：不换模型也能修正推断

稳健回归换的是系数估计，但有时你并不想换模型，比如你用的是一个比较复杂的 GLM 或固定效应模型，重写稳健版本成本太高；或者你只关心系数显著性和置信区间，不关心系数本身被拉动。这时候有一个便宜得多的方案：保留原模型，只替换标准误。这就是稳健标准误（也叫异方差稳健标准误）的用武之地。R 语言里最常用的是sandwich包。另一条线是多元场景：当你有一堆变量、想看哪些观测是多元离群点时，稳健协方差矩阵估计covMcd比逐个变量做单变量检验可靠得多。

4.1 sandwich 包：OLS 系数搭配 vcovHC 稳健标准误

回归模型的经典假设是误差同方差。一旦数据里有离群值或异方差，OLS 系数虽然可能还大致可靠，但标准误会失真，导致 t 检验和 p 值一起失真。vcovHC的作用是不动系数，只换协方差矩阵的估计方式。

library(sandwich) library(lmtest) fit_ols <- lm(y ~ x, data = dat) # 默认 HC3 稳健标准误 coeftest(fit_ols, vcov = vcovHC(fit_ols, type = "HC3")) # 对比普通标准误 coeftest(fit_ols)

coeftest来自lmtest包，它的作用是把回归模型的系数、稳健标准误、t 值、p 值重新组装成一张熟悉的表。vcovHC的type参数控制小样本校正方式：HC0是最基础的 sandwich 估计，不校正；HC1做自由度校正，等价于很多软件里的"稳健标准误"；HC2考虑杠杆值；HC3对高杠杆点惩罚最重，是sandwich包的默认值，样本量不大时最推荐。如果你的数据是时间序列，误差还可能自相关，vcovHAC或者NeweyWest是替代方案，原理类似，只是把相关结构也纳入估计。需要明确一点：vcovHC只改变推断，不改变系数，离群值照旧拉动 OLS 系数。所以这个方案适合"系数不是重点、显著性才是重点"的场景。

4.2 covMcd 稳健协方差：多元离群点检测与马氏距离

单变量离群可以用箱线图找，但多元离群是"组合异常"：每个变量单独看都正常，组合起来却偏离主体。这时要用马氏距离，而马氏距离依赖协方差矩阵，普通协方差矩阵又被离群点污染，形成死循环。covMcd的思路是挑出数据中"最紧凑"的一半（或指定比例）子集，用这部分干净数据估计中心和协方差，再做全样本的马氏距离检验。

library(robustbase) # 用多个变量做稳健协方差估计，这里取 x, y 两个演示 mcd <- covMcd(dat, alpha = 0.75) # 稳健中心与稳健协方差 mcd$center mcd$cov # 计算全样本的稳健马氏距离 md_rob <- mahalanobis(dat, center = mcd$center, cov = mcd$cov) # 卡方分布 97.5% 分位数作为离群阈值 cutoff <- qchisq(0.975, df = ncol(dat)) flag <- md_rob > cutoff cat("被标记为多元离群点的数量:", sum(flag), "\n")

alpha = 0.75表示用 75% 的最紧凑子集估计协方差，alpha 越小越稳健、但估计方差越大，默认 0.5 时崩溃点最高。mahalanobis计算出每个点到稳健中心的距离，在多元正态假设下近似服从卡方分布，自由度等于变量数，所以用qchisq(0.975, df = ncol(dat))做阈值。这里的df要写ncol(dat)，不是样本量，这是最常见的错误。被标记出来的观测就是同时影响均值和协方差的可疑点，拿到这些点之后回到rlm或lmrob里做二次确认。

4.3 用稳健回归还是只换稳健标准误：判断流程

我给一个可执行的判断顺序。第一，系数本身是否敏感：用第 2.3 节的方法对比 OLS 和剔除极端值后的系数，如果系数变化大，说明离群点在实质性拉动结果，必须上稳健回归。第二，你是不是只关心显著性：如果研究问题就是检验某个变量系数是否为 0，而系数方向稳定，那就用vcovHC修正推断，成本最低。第三，模型复杂度：如果是混合效应模型、结构方程这类没有现成稳健版本的模型，sandwich的思路可以扩展到clubSandwich等扩展包；如果模型可以换成线性形式，lmrob仍然是最省事的方案。第四，样本量：小样本下稳健回归的标准误本身方差就大，与其只换标准误，不如同时报告 OLS、rlm、lmrob 三列系数，让读者自己看稳定性。这套流程比"一律上稳健回归"或者"一律只换标准误"都更靠谱。

5. 避坑：稳健性估计在 R 里的五个常见问题

稳健性估计的坑不在理论多深，而在包的行为、默认参数和边界条件上。下面五条都是实际运行中反复遇到的，按"现象 → 原因 → 解决"写清楚。

5.1 lmrob 反复报 scale 不收敛

现象：lmrob()跑完弹出警告，类似 "Scale estimation did not converge"，甚至直接报Error in lmrob.S(...)。原因：lmrob内部的 S 估计在高污染比例、高杠杆点集群或者 x 变量量纲差异极大的情况下，迭代找不到稳定的残差尺度解。解决：先把数值变量做标准化，尤其是量纲差几个数量级的特征；再显式提高迭代上限，lmrob(..., control = lmrob.control(maxit = 500, maxit.scale = 500))；如果还不行，把setting参数从默认改成"KS2014"或"S"试一组初始值策略。换数据比换参数更优先：先画出杠杆图，确认是不是有 x 方向极端点集群。

5.2 默认 rlm 参数对高杠杆点无效

现象：rlm()用了默认的 Huber psi，结果系数和 OLS 差不多，离群值依然把模型拉着跑。原因：Huber 的 psi 是单调函数，只限制大残差的惩罚增长，不限制高杠杆点的影响力，而且默认k = 1.345在污染较重时阈值偏大。解决：改用psi = psi.bisquare, c = 4.685，让超过阈值的点权重归零；同时检查 x 的杠杆，hatvalues(fit_ols)找出高杠杆点，放进lmrob里对比。如果 bisquare 迭代震荡，把maxit调到 100，并换scale.est = "Huber"试试。

5.3 把稳健标准误当成稳健回归

现象：用了vcovHC之后，系数完全没有变化，有人就以为数据没问题，直接下了结论。原因：没分清"稳健推断"和"稳健估计"。vcovHC只换协方差矩阵，系数还是 OLS 估计出来的，离群值照样在拉动系数。解决：先跑第 2.3 节的系数敏感性检查；如果系数差异大，先上rlm或lmrob，再用vcovHC作为附加结果。正确的组合方式是：稳健回归 + 稳健标准误，而不是二选一。也可以在论文或报告里同时报告两套系数，让审阅者看到差异本身。

5.4 小样本直接套 MM 估计

现象：n = 30 甚至 n = 20，跑lmrob得到的结果和rlm差很多，标准误也忽大忽小。原因：崩溃点 0.5 是渐近性质，小样本下稳健估计的方差会明显膨胀，S 估计初始化阶段对随机种子很敏感，同样的数据换一个 seed 结果就不一样。解决：n < 50 时优先用rlm(psi = psi.huber)这种参数较少的估计器，少一个 S 估计阶段就少一层随机性；如果必须用lmrob，把setting = "KS2014"固定下来，并且跑一个多重种子稳定性检查，sapply(1:20, function(s){ set.seed(s); coef(lmrob(y ~ x, dat))["x"] })，看结果波动范围。波动大就说明你的样本撑不起稳健估计。

5.5 不固定随机种子，结果无法复现

现象：同一份数据，两次运行lmrob()得到略微不同的系数，covMcd()的标记点也不一样。原因：两个函数内部都有随机初始值或重抽样步骤，lmrob的 S 估计初始化、covMcd的 C-step 迭代都依赖随机子集。解决：在脚本开头固定set.seed(42)，并且把 R 版本、packageVersion("robustbase")记录在分析文档里。这不是玄学，是稳健估计在计算层面自带的行为，任何人都应该这样处理。报告结果的时候附上 seed 和包版本，别人才能在你的数字上做进一步对比。

6. 验证稳健性估计：污染度扫描与自举置信区间

稳健性估计的结论不能只靠一次拟合就信。我习惯做两件事：一是污染度扫描，把数据被人为污染的比例从 0% 扫到 30%，看 OLS、rlm、lmrob 的斜率估计各自怎么偏离真值；二是自举法，对稳健估计的系数分布做非参数 bootstrap，得到不依赖分布假设的置信区间。

library(boot) # 污染度扫描：每个污染水平下重复 200 次，比较斜率偏差 sweep_contam <- function(c_rate, n = 200, reps = 200) { b_ols <- b_rlm <- b_mm <- numeric(reps) for (i in 1:reps) { x <- rnorm(n, 10, 3); e <- rnorm(n, 0, 2) y <- 2 + 1.5 * x + e if (c_rate > 0) { n_out <- round(n * c_rate) idx <- sample(n, n_out) y[idx] <- y[idx] + rnorm(n_out, 0, 25) } d <- data.frame(x, y) b_ols[i] <- coef(lm(y ~ x, d))["x"] b_rlm[i] <- coef(rlm(y ~ x, d, psi = psi.bisquare, c = 4.685))["x"] b_mm[i] <- coef(lmrob(y ~ x, d))["x"] } data.frame(c_rate = c_rate, ols_bias = mean(b_ols) - 1.5, rlm_bias = mean(b_rlm) - 1.5, mm_bias = mean(b_mm) - 1.5) } results <- do.call(rbind, lapply(seq(0, 0.3, 0.05), sweep_contam)) print(results, digits = 3)

这段代码的核心是三层循环：外层扫污染比例，中层设置随机种子对应的仿真，内层拟合三个模型。results里每一行是某个污染水平下 200 次重复的平均偏差。你会看到 OLS 的偏差随污染水平快速变大，Huber 在 10% 内还行、超过 15% 开始吃力，lmrob一直压得很低。这就是稳健性估计"值不值得上"的最直观证据：用你自己数据的污染比例预估，看三种方法在该点的偏差差距。

自举法给系数一个更可信的区间。boot包的做法是把重抽样写成一个函数，返回值是rlm的斜率：

boot_rlm <- function(data, indices) { d <- data[indices, ] fit <- rlm(y ~ x, data = d, psi = psi.bisquare, c = 4.685) unname(coef(fit)["x"]) } set.seed(42) boot_res <- boot(dat, boot_rlm, R = 999) boot.ci(boot_res, type = "perc", index = 2)

boot(dat, boot_rlm, R = 999)对dat做 999 次有放回重抽样，每次重抽样都重新拟合一次rlm。boot.ci的type = "perc"给出百分位置信区间，index = 2对应斜率的位置。如果bca区间可用也可以试，但需要足够多的重抽样次数，不然会警告。拿到这个区间之后，和 OLS 的普通置信区间、lmrob的稳健标准误区间放在一起，你就能看出：稳健估计不只是把系数拉回真值，它的不确定区间往往也更诚实。这也是我每次跑稳健回归的收尾动作——先证明它比 OLS 稳，再给出它的区间。希望帮到你。

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