簡介R語言穩(wěn)健性估計(jì)實(shí)例分析資源面向數(shù)據(jù)分析、統(tǒng)計(jì)建模及回歸診斷學(xué)習(xí)者。壓縮包共1個(gè)pptx文件大小僅716KB以幻燈形式系統(tǒng)展示線性回歸診斷與穩(wěn)健回歸的完整思路。內(nèi)容從lm()基礎(chǔ)擬合與plot()四聯(lián)診斷圖出發(fā)逐步講解殘差、異常點(diǎn)、高杠桿點(diǎn)與強(qiáng)影響點(diǎn)的判別方法涵蓋學(xué)生化殘差、帽子矩陣及Cook距離等關(guān)鍵指標(biāo)的計(jì)算與應(yīng)用同時(shí)理清三類特殊點(diǎn)的聯(lián)系與區(qū)別在此基礎(chǔ)上引入Huber與Bisquare兩種M估計(jì)穩(wěn)健回歸方法通過實(shí)例演示如何在異常值存在時(shí)進(jìn)行加權(quán)迭代獲得更可靠的參數(shù)估計(jì)。整體框架緊湊適合教學(xué)演示、課后復(fù)習(xí)或項(xiàng)目參考能幫助讀者快速構(gòu)建回歸穩(wěn)健性分析的知識體系。目前已有1129人學(xué)習(xí)對于需要處理含異常值數(shù)據(jù)的分析人員具有較高參考價(jià)值。1. R 語言穩(wěn)健性估計(jì)從 lm() 到 rlm() 的完整實(shí)例分析做回歸分析時(shí)我經(jīng)常碰到一種場景數(shù)據(jù)里混進(jìn)了幾個(gè)“不老實(shí)”的點(diǎn)普通最小二乘回歸OLS的結(jié)果被它們牽著鼻子走模型系數(shù)變得面目全非。R 語言里處理這類問題有一套成熟的工具鏈從lm()擬合、plot(lm.fit1)出四張?jiān)\斷圖到cooks.distance()計(jì)算 Cook 距離再到穩(wěn)健回歸中的 Huber 和 Bisquare M 估計(jì)每一步都有對應(yīng)的函數(shù)和判斷標(biāo)準(zhǔn)。這篇文章圍繞一套完整的 R 實(shí)例分析展開包含可直接運(yùn)行的 R 代碼和一份 crime 數(shù)據(jù)集的分析流程適合正在做回歸診斷、異常值處理或需要提高模型穩(wěn)健性的數(shù)據(jù)分析師和統(tǒng)計(jì)專業(yè)學(xué)生。你將看到普通殘差、學(xué)生化殘差、杠桿率、Cook 距離這幾個(gè)概念如何串成一條識別異常點(diǎn)的完整鏈路以及rlm()在 Huber 和 Bisquare 兩種權(quán)重函數(shù)下的實(shí)際表現(xiàn)——這些內(nèi)容在多數(shù)教材里只講公式很少告訴你參數(shù)怎么選、輸出怎么讀、哪些“經(jīng)驗(yàn)分界點(diǎn)”其實(shí)有爭議。文章會以一份真實(shí)可跑通的 R 代碼為主線把每個(gè)函數(shù)的作用、每段輸出的含義、每個(gè)閾值的由來都拆開講清楚。2. 從普通殘差到學(xué)生化殘差異常點(diǎn)的識別邏輯與帽子矩陣2.1 普通殘差為什么不能直接用方差不等齊問題任何一本回歸分析教材都會告訴你殘差是觀測值Y與預(yù)測值?的差表達(dá)式為e Y - ?。但實(shí)際用 R 做診斷時(shí)直接比較普通殘差的大小是有問題的。問題出在方差上普通殘差的方差不是常數(shù)它依賴于帽子矩陣的對角線元素h_ii具體形式是Var(e_i) σ2(1 - h_ii)。這意味著什么不同觀測點(diǎn)的殘差天然具有不同的方差如果直接比較e_i的絕對值大小那些h_ii較大的點(diǎn)即遠(yuǎn)離自變量均值的點(diǎn)殘差方差更小同樣的偏差會被放大從而被誤判為異常點(diǎn)。我一般會在 R 里這樣獲取普通殘差# 讀取數(shù)據(jù)并擬合普通線性回歸模型 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) # 查看前六個(gè)殘差 head(resid_ols)這段代碼中resid()函數(shù)提取 OLS 回歸的普通殘差fitted()提取模型對每個(gè)樣本的預(yù)測值。attach(c1)把數(shù)據(jù)框的列變量直接暴露到工作環(huán)境中方便后續(xù)直接引用Weight和Height但要注意使用后建議用detach(c1)釋放避免變量名沖突。plot(lm.fit1)是診斷的第一道工序它一次生成四幅圖殘差對擬合值圖、殘差的正態(tài) Q-Q 圖、標(biāo)準(zhǔn)化殘差絕對值平方根對擬合值圖、Cook 距離圖。這里面第三幅圖橫軸是擬合值縱軸是sqrt(|standardized residuals|)主要用來檢查方差齊性。如果你看到散點(diǎn)呈現(xiàn)漏斗形分布說明方差不穩(wěn)定這時(shí)普通殘差的可比性進(jìn)一步下降。2.2 帽子矩陣與杠桿率h_ii 如何刻畫點(diǎn)的“偏遠(yuǎn)程度”杠桿率衡量的是自變量X對自身均值的偏異程度。公式為h_ii (1/n) (X_i - X?) (XX)^{-1} (X_i - X?)從公式可以直接讀出兩層含義第一項(xiàng)1/n是基礎(chǔ)杠桿所有點(diǎn)共享第二項(xiàng)是第i個(gè)點(diǎn)到樣本中心X?的 Mahalanobis 距離。在樣本空間中h_ii較大的點(diǎn)位于自變量空間的邊緣它們可能把回歸線拉向自己對回歸系數(shù)的 LS 估計(jì)影響可能很大。在 R 中提取杠桿率有很多路徑常見做法是# 通過 lm.influence 獲取帽子矩陣對角線元素 H - hatvalues(lm.fit1) # 查看杠桿率最高的幾個(gè)樣本 head(sort(H, decreasing TRUE), 5) # 結(jié)合模型矩陣手動計(jì)算杠桿率 X - model.matrix(lm.fit1) H_manual - diag(X %*% solve(t(X) %*% X) %*% t(X))代碼中hatvalues()返回帽子矩陣的對角線元素是官方推薦做法。model.matrix()提取設(shè)計(jì)矩陣X包括截距列和自變量列然后用矩陣運(yùn)算手動復(fù)現(xiàn)X(XX)^{-1}X的對角線。手動計(jì)算的目的是驗(yàn)證對帽子矩陣的理解實(shí)際項(xiàng)目中直接用hatvalues()即可。注意1/n這一項(xiàng)說明即使所有自變量都等于均值杠桿率也至少是1/n所以看杠桿率時(shí)不要只看絕對值還要結(jié)合2p/n或3p/n這類經(jīng)驗(yàn)閾值判斷。2.3 學(xué)生化殘差的計(jì)算公式拆解與 R 實(shí)現(xiàn)由于普通殘差存在方差不齊的問題需要標(biāo)準(zhǔn)化后比較。學(xué)生化殘差的形式是r_i e_i / (s * sqrt(1 - h_ii))其中s是剩余標(biāo)準(zhǔn)差h_ii是帽子矩陣對角線元素。從公式可以看出學(xué)生化殘差同時(shí)考慮了殘差本身的偏差程度和杠桿率的影響。h_ii越大分母越小同一個(gè)殘差對應(yīng)的學(xué)生化殘差越大。在 R 中可以直接用rstandard()或rstudent()得到內(nèi)部學(xué)生化殘差和外部學(xué)生化殘差# 內(nèi)部學(xué)生化殘差使用當(dāng)前模型的誤差方差估計(jì) r_int - rstandard(lm.fit1) # 外部學(xué)生化殘差刪除第i個(gè)點(diǎn)后重新估計(jì)誤差方差 r_ext - rstudent(lm.fit1) # 判斷哪些點(diǎn)超過閾值 outlier_flag - abs(r_ext) 3 sum(outlier_flag)rstandard()計(jì)算時(shí)使用包含所有樣本的誤差方差估計(jì)rstudent()則對每個(gè)點(diǎn)執(zhí)行“刪除一個(gè)樣本后再估計(jì)方差”的策略對異常點(diǎn)更敏感。經(jīng)驗(yàn)上外部學(xué)生化殘差絕對值大于 3 的點(diǎn)值得高度關(guān)注。代碼中最后一行的sum()統(tǒng)計(jì)異常點(diǎn)數(shù)量方便批量篩查。2.4 避坑學(xué)生化殘差與普通殘差的三個(gè)典型誤用現(xiàn)象直接比較普通殘差的大小把e_i最大的幾個(gè)點(diǎn)當(dāng)作異常點(diǎn)結(jié)果剔除后模型反而變得更差某些正常點(diǎn)被誤刪。原因普通殘差方差不齊h_ii較大的點(diǎn)天然殘差方差更小同樣的偏離程度表現(xiàn)為更大的e_i導(dǎo)致高杠桿點(diǎn)被優(yōu)先標(biāo)記為異常點(diǎn)而真正的離群點(diǎn)可能因?yàn)楦軛U率低被漏掉。解決用rstandard()或rstudent()替代普通殘差。學(xué)生化殘差分母中加入了sqrt(1 - h_ii)修正了方差不等的影響。我在實(shí)際項(xiàng)目中基本只用rstudent()它對單個(gè)異常點(diǎn)更敏感。現(xiàn)象abs(r_ext) 2標(biāo)記出大量點(diǎn)把閾值放寬到 2 后異常點(diǎn)比例超過 10%模型被削掉太多樣本。原因樣本量較大時(shí)學(xué)生化殘差的分布接近t分布在n 50時(shí)約 5% 的點(diǎn)可能超過 2但這不代表它們是異常點(diǎn)。閾值設(shè)置過松會把正常波動當(dāng)成異常。解決以abs(r_ext) 3作為首要關(guān)注線同時(shí)結(jié)合 Cook 距離判斷強(qiáng)影響性。不要只依據(jù)單一指標(biāo)刪點(diǎn)應(yīng)該綜合殘差、杠桿率、Cook 距離三維度。現(xiàn)象刪除了所有學(xué)生化殘差超閾值的點(diǎn)后重新擬合發(fā)現(xiàn)刪點(diǎn)前模型系數(shù)還在合理范圍刪點(diǎn)后某個(gè)自變量變得不顯著或系數(shù)符號反轉(zhuǎn)。原因一個(gè)點(diǎn)既是異常點(diǎn)又是強(qiáng)影響點(diǎn)時(shí)它對系數(shù)的拉動作用可能掩蓋了其他點(diǎn)的模式。盲目刪除所有異常點(diǎn)可能破壞本來穩(wěn)定的數(shù)據(jù)結(jié)構(gòu)。解決先看 Cook 距離優(yōu)先關(guān)注“影響大”的點(diǎn)而不是“偏差大”的點(diǎn)。異常點(diǎn)不一定有強(qiáng)影響高杠桿點(diǎn)也不一定是強(qiáng)影響點(diǎn)需要區(qū)分對待。3. Cook 距離與強(qiáng)影響點(diǎn)綜合杠桿率和殘差的判斷標(biāo)準(zhǔn)3.1 Cook 距離公式拆解為什么它同時(shí)包含 h_ii 和 r_iCook 距離是回歸診斷中使用頻率最高的影響度量指標(biāo)其公式為D_i (r_i2 / p) * (h_ii / (1 - h_ii))其中r_i是第i個(gè)點(diǎn)的學(xué)生化殘差p是模型中參數(shù)個(gè)數(shù)含截距h_ii是杠桿率。這個(gè)結(jié)構(gòu)很有深意第一項(xiàng)r_i2 / p度量殘差偏離程度第二項(xiàng)h_ii / (1 - h_ii)是杠桿率的單調(diào)變換。兩個(gè)因子相乘意味著一個(gè)點(diǎn)只有同時(shí)具備“殘差大”和“杠桿高”兩個(gè)特征時(shí)Cook 距離才會顯著。單純殘差大但杠桿低或者杠桿高但殘差小D_i都不會太大。這與強(qiáng)影響點(diǎn)的定義高度吻合強(qiáng)影響點(diǎn)是指剔除后對回歸系數(shù)估計(jì)有顯著效應(yīng)的觀測值。在 R 中的計(jì)算方式非常直接# 使用基本包的 cooks.distance 函數(shù) d1 - cooks.distance(ols) # 查看 Cook 距離最大的樣本 which.max(d1) # 結(jié)合學(xué)生化殘差和杠桿率構(gòu)成診斷矩陣 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()返回每個(gè)樣本的 Cook 距離stdres()提取標(biāo)準(zhǔn)化殘差hatvalues()提取杠桿率。構(gòu)建的數(shù)據(jù)框把三個(gè)核心診斷量并列展示按 Cook 距離降序排列后可以直觀看到哪些點(diǎn)對模型影響最大。這里的ols是之前l(fā)m(crime ~ poverty single, data cdata)的擬合結(jié)果在 UCLA 的crime.dta數(shù)據(jù)集上運(yùn)行分析crime與poverty、single兩個(gè)自變量的關(guān)系。3.2 經(jīng)驗(yàn)分界點(diǎn) 4/n 的由來與爭議Cook 距離的判斷閾值在學(xué)術(shù)界一直存在爭議。最常用的經(jīng)驗(yàn)分界點(diǎn)是4/n其中n是樣本量。在 R 中篩選強(qiáng)影響點(diǎn)的標(biāo)準(zhǔn)寫法是# 按 4/n 閾值篩選強(qiáng)影響點(diǎn) n - nrow(cdata) influential - cdata[d1 4 / n, ] influential # 同時(shí)也可以參考 F 分布的分位數(shù) qf_threshold - qf(0.5, df1 2, df2 n - 2) influential_f - cdata[d1 qf_threshold, ]4/n是經(jīng)驗(yàn)法則來源于 Cook 距離與 F 分布近似關(guān)系中取F(0.5, p, n-p)的近似結(jié)果。另一種做法是用qf(0.5, p, n-p)直接計(jì)算 F 分布 50% 分位數(shù)作為閾值這在p2時(shí)通常比4/n略寬松。實(shí)際問題中我一般兩種都跑一遍把落在兩個(gè)閾值之間但又不算極端的樣本標(biāo)記為“重點(diǎn)關(guān)注”。3.3 實(shí)際分析crime 數(shù)據(jù)集中第 9、25、51 號樣本的處理在 crime 數(shù)據(jù)集上運(yùn)行plot(ols, las 1)會生成四張?jiān)\斷圖。從殘差圖和 Cook 距離圖可以清晰看到第 9、25、51 號觀測值位于邊緣位置。進(jìn)一步用數(shù)值確認(rèn)# 查看這3個(gè)樣本的具體診斷值 target_ids - c(9, 25, 51) diag_matrix[target_ids, ] # 輸出這些樣本的原始數(shù)據(jù) cdata[target_ids, ]輸出的診斷矩陣顯示這三個(gè)點(diǎn)的 Cook 距離都超過了4/51的閾值標(biāo)準(zhǔn)化殘差絕對值也偏高。此時(shí)面臨一個(gè)典型決策場景如果直接采用 OLS你可能會傾向于刪除這三行數(shù)據(jù)再重新擬合但如果刪除后模型系數(shù)變化巨大說明這些點(diǎn)是強(qiáng)影響點(diǎn)但未必是“錯(cuò)誤數(shù)據(jù)”。穩(wěn)健回歸提供了第三條路不剔除樣本而是降低它們的權(quán)重。3.4 避坑Cook 距離閾值的兩個(gè)常見翻車現(xiàn)場現(xiàn)象使用4/n閾值篩出 5 個(gè)強(qiáng)影響點(diǎn)全部刪除后重新擬合發(fā)現(xiàn)某個(gè)自變量系數(shù)符號反向擬合優(yōu)度下降。原因強(qiáng)影響點(diǎn)不一定都是“壞點(diǎn)”。如果這個(gè)點(diǎn)代表了真實(shí)存在的特殊子群體比如高收入低犯罪率的城市刪除它會讓模型喪失對這類群體的解釋能力。4/n是經(jīng)驗(yàn)閾值樣本量小或自變量維度高時(shí)容易誤判。解決先記錄強(qiáng)影響點(diǎn)對應(yīng)的實(shí)際業(yè)務(wù)含義再決定是否刪除。通常我會保留這些樣本改用穩(wěn)健回歸或加權(quán)回歸讓數(shù)據(jù)自己決定權(quán)重?，F(xiàn)象plot(lm.fit1)四張圖中 Cook 距離圖看起來沒有超過紅虛線但手工計(jì)算cooks.distance()卻發(fā)現(xiàn)值超過4/n兩套結(jié)果不一致。原因plot()函數(shù)繪制的 Cook 距離圖縱軸范圍可能被自動縮放紅虛線是 R 根據(jù) Cook 距離分布計(jì)算的可視化閾值而不是嚴(yán)格的4/n邊界。兩種呈現(xiàn)邏輯不同導(dǎo)致肉眼判斷與數(shù)值判斷沖突。解決以cooks.distance()的數(shù)值結(jié)果為準(zhǔn)plot()圖只作為初步篩查。數(shù)值篩選后用identify()或which()定位具體樣本ID再回到業(yè)務(wù)層面判斷。4. rlm() 實(shí)現(xiàn)穩(wěn)健回歸Huber 與 Bisquare 兩種 M 估計(jì)的完整實(shí)戰(zhàn)4.1 為什么選擇 rlm最小二乘在異常點(diǎn)面前的兩個(gè)困境最小二乘估計(jì)的目標(biāo)是使殘差平方和最小這意味著一個(gè)大殘差點(diǎn)會以平方級別拉動回歸線。面對異常點(diǎn)和高杠桿點(diǎn)時(shí)OLS 有兩個(gè)困境第一如果異常點(diǎn)來自數(shù)據(jù)錄入錯(cuò)誤理論上應(yīng)該剔除但數(shù)據(jù)分析者很難有充分證據(jù)證明“這個(gè)點(diǎn)一定是錯(cuò)的”第二如果異常點(diǎn)來自另一個(gè)總體或特殊子群體直接刪除會造成樣本選擇偏差。穩(wěn)健回歸的思路是在“完全剔除”與“一視同仁”之間折中對殘差較大的觀測值賦予較低權(quán)重對正常樣本保留高權(quán)重。rlm()是 MASS 包中的核心函數(shù)實(shí)現(xiàn)了 M 估計(jì)的迭代重復(fù)加權(quán)最小二乘算法。其基本流程是先用 OLS 得到初始?xì)埐罡鶕?jù)殘差大小計(jì)算觀測權(quán)重再用加權(quán)最小二乘更新系數(shù)然后重新計(jì)算殘差和權(quán)重迭代直到收斂。權(quán)重函數(shù)的選擇決定了穩(wěn)健性的具體形式。4.2 Huber 方法的權(quán)重函數(shù)與參數(shù)選擇Huber 方法的權(quán)重函數(shù)是分段函數(shù)w(e) 1當(dāng)|e| cw(e) c / |e|當(dāng)|e| c其中c是截?cái)喑?shù)R 中默認(rèn)取1.345。這意味著殘差在閾值內(nèi)的觀測獲得權(quán)重 1殘差超過閾值的觀測權(quán)重隨殘差增大而遞減。Huber 估計(jì)對中等程度的異常值表現(xiàn)穩(wěn)健同時(shí)保留了較高的統(tǒng)計(jì)效率。在 R 中的用法# 加載 MASS 包 library(MASS) # Huber 方法的 M 估計(jì) rr.huber - rlm(crime ~ poverty single, data cdata) # 查看模型摘要 summary(rr.huber) # 查看每個(gè)觀測的最終權(quán)重 weights_huber - rr.huber$w head(sort(weights_huber, decreasing FALSE), 10)summary(rr.huber)輸出與lm()類似包含系數(shù)估計(jì)和t值但注意這里不展示 F 統(tǒng)計(jì)量和 R2因?yàn)榈訖?quán)過程讓這些統(tǒng)計(jì)量的解釋變得復(fù)雜。rr.huber$w保存了每個(gè)觀測的最終權(quán)重權(quán)重最小的點(diǎn)就是被降權(quán)最厲害的點(diǎn)。4.3 Bisquare 方法的權(quán)重函數(shù)與參數(shù)選擇Bisquare也常稱為 Tukeys biweight方法的權(quán)重函數(shù)是w(e) (1 - (e/c)2)2當(dāng)|e| cw(e) 0當(dāng)|e| c與 Huber 方法不同Bisquare 給所有非零殘差的觀測都賦予遞減權(quán)重殘差超過c的觀測權(quán)重直接歸零。R 中默認(rèn)c 4.685。這意味著 Bisquare 比 Huber 更“激進(jìn)”它可以完全剔除極端異常點(diǎn)的影響而 Huber 對極端殘差仍然保留c/|e|的微小權(quán)重。# Bisquare 方法的 M 估計(jì) rr.bisq - rlm(crime ~ poverty single, data cdata, method MM) # 或者顯式指定 psi 函數(shù)為 bisquare rr.bisq2 - rlm(crime ~ poverty single, data cdata, psi psi.bisquare) # 查看權(quán)重分布 summary(rr.bisq$w)代碼中method MM表示使用 MM 估計(jì)它結(jié)合了高分解值和高效率特性是處理強(qiáng)影響點(diǎn)時(shí)的推薦選擇。psi psi.bisquare顯式指定所用的psi函數(shù)MASS 包中內(nèi)置了psi.huber和psi.bisquare。MM 估計(jì)在初始化階段使用高分解值的估計(jì)方法然后進(jìn)入 Bisquare 迭代比默認(rèn)的 M 估計(jì)更穩(wěn)健。4.4 權(quán)重結(jié)果對比同一批樣本在兩種方法下的待遇差異將兩種方法的權(quán)重提取出來對比是理解穩(wěn)健回歸最直觀的方式# 合并兩種權(quán)重進(jìn)行對比 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) ) # 查看權(quán)重最低的10個(gè)樣本 head(weight_compare[order(weight_compare$huber_w), ], 10) # 計(jì)算兩種權(quán)重與 OLS 標(biāo)準(zhǔn)化殘差的相關(guān)性 cor(weight_compare$huber_w, abs(weight_compare$std_resid_ols)) cor(weight_compare$bisq_w, abs(weight_compare$std_resid_ols))通常你會發(fā)現(xiàn)Huber 方法中權(quán)重最小的點(diǎn)對應(yīng)原始 OLS 標(biāo)準(zhǔn)化殘差最大的點(diǎn)但權(quán)重不會降到 0Bisquare 方法則可能將極端殘差點(diǎn)權(quán)重直接置零。兩個(gè)模型的系數(shù)估計(jì)差異反映了穩(wěn)健回歸的“折中”程度。Huber 適合你懷疑異常點(diǎn)有少量信息但不愿完全放棄的場景Bisquare 適合你認(rèn)為部分點(diǎn)真的來自其他總體的場景。4.5 避坑rlm() 使用中的四個(gè)高頻報(bào)錯(cuò)與處理現(xiàn)象rlm()運(yùn)行后提示convergence相關(guān)警告或者迭代次數(shù)未達(dá)到默認(rèn)上限就停止結(jié)果似乎仍未穩(wěn)定。原因M 估計(jì)的迭代是從 OLS 初始值開始的如果初始模型中有極端強(qiáng)影響點(diǎn)權(quán)重函數(shù)可能在某些點(diǎn)產(chǎn)生周期性振蕩迭代難以收斂。默認(rèn)最大迭代次數(shù)可能不足。解決增加迭代次數(shù)或調(diào)整初始值?？梢詡魅雖axit 100參數(shù)也可以先利用lm()擬合后剔除極端 Cook 距離點(diǎn)再用剩余樣本的系數(shù)作為初值?，F(xiàn)象rlm(crime ~ poverty single, data cdata)報(bào)錯(cuò)提示variable lengths differ或者NA/NaN/Inf in foreign function call。原因數(shù)據(jù)中存在缺失值。rlm()默認(rèn)使用na.omit處理缺失值但部分情況下數(shù)據(jù)框中的NA會在權(quán)重計(jì)算中引發(fā)錯(cuò)誤。解決擬合前手動執(zhí)行cdata - na.omit(cdata)同時(shí)檢查是否存在Inf值。如果某個(gè)自變量的分布嚴(yán)重偏態(tài)考慮先做對數(shù)變換再進(jìn)入模型?，F(xiàn)象擬合成功但summary(rr.huber)輸出的系數(shù)與lm()差別不大懷疑穩(wěn)健回歸沒有起作用。原因數(shù)據(jù)集中本身沒有嚴(yán)重的異常點(diǎn)或高杠桿點(diǎn)穩(wěn)健回歸和 OLS 自然結(jié)果接近。這不是 bug而是正?，F(xiàn)象。穩(wěn)健回歸的價(jià)值在數(shù)據(jù)“臟”的時(shí)候才體現(xiàn)。解決在擬合前先畫出散點(diǎn)圖或執(zhí)行診斷矩陣確認(rèn)數(shù)據(jù)中確實(shí)存在候選異常點(diǎn)。如果診斷結(jié)果表明數(shù)據(jù)干凈直接報(bào) OLS 結(jié)果即可?，F(xiàn)象Bisquare 方法擬合后大量觀測權(quán)重為 0模型的有效樣本量大幅下降標(biāo)準(zhǔn)誤增大。原因psi.bisquare的默認(rèn)截?cái)喑?shù)c 4.685對應(yīng)的殘差閾值是在正態(tài)誤差假設(shè)下確定的如果數(shù)據(jù)中存在多個(gè)相互靠近的異常點(diǎn)遮蔽效應(yīng)可能導(dǎo)致過多樣本被降權(quán)。解決改用method MM提高分解值或者適當(dāng)調(diào)大c值比如psi psi.bisquare, c 5.5。但注意調(diào)大c會降低穩(wěn)健性需要權(quán)衡。5. 完整 R 代碼實(shí)戰(zhàn)從 OLS 診斷到穩(wěn)健回歸的參數(shù)對比5.1 數(shù)據(jù)讀取與模型擬合的完整流程結(jié)合前文提到的crime.dta數(shù)據(jù)集完整流程從讀取外文格式數(shù)據(jù)開始。R 中讀取 Stata 格式數(shù)據(jù)需要使用foreign包# 加載所需包 require(foreign) require(MASS) # 讀取 Stata 格式數(shù)據(jù) cdata - read.dta(https://stats.idre.ucla.edu/stat/data/crime.dta) # 查看數(shù)據(jù)結(jié)構(gòu) str(cdata) names(cdata) # 擬合普通最小二乘回歸 ols - lm(crime ~ poverty single, data cdata) # 輸出模型摘要 summary(ols)read.dta()是讀取 Stata 數(shù)據(jù)文件的標(biāo)準(zhǔn)函數(shù)其網(wǎng)絡(luò)路徑直接加載數(shù)據(jù)。str(cdata)查看各變量的類型和取值分布確保crime、poverty、single都是數(shù)值型。summary(ols)輸出的系數(shù)表中需要重點(diǎn)關(guān)注poverty和single的估計(jì)值及顯著性。5.2 四圖診斷與數(shù)值診斷的配合診斷不能只依賴plot()生成的圖形還需要數(shù)值輸出來確定具體樣本編號。完整流程如下# 四圖診斷 opar - par(mfrow c(2, 2), oma c(0, 0, 1.1, 0)) plot(ols, las 1) # 計(jì)算 Cook 距離和標(biāo)準(zhǔn)化殘差 d1 - cooks.distance(ols) r - stdres(ols) h - hatvalues(ols) # 構(gòu)建診斷矩陣 a - cbind(cdata, d1, r, h) # 按 4/n 閾值篩選 n - nrow(cdata) a[d1 4 / n, ]代碼中par(mfrow c(2, 2))將圖形區(qū)域分割成 2x2 的網(wǎng)格四張?jiān)\斷圖依次排列。cbind()將原始數(shù)據(jù)與三個(gè)診斷量合并成新數(shù)據(jù)框方便篩選和查看。a[d1 4 / n, ]篩選出 Cook 距離超閾值的全部樣本輸出包括原始變量和診斷量可以直接對照樣本 ID 查看業(yè)務(wù)含義。5.3 穩(wěn)健回歸與 OLS 系數(shù)對比表# OLS 系數(shù) coef_ols - coef(ols) # Huber 穩(wěn)健回歸系數(shù) coef_huber - coef(rr.huber) # Bisquare 穩(wěn)健回歸系數(shù) coef_bisq - coef(rr.bisq) # 合并結(jié)果生成對比表 compare_table - data.frame( OLS round(coef_ols, 4), Huber round(coef_huber, 4), Bisquare round(coef_bisq, 4) ) print(compare_table)對比表的價(jià)值在于直觀展示三種方法對同一批數(shù)據(jù)的系數(shù)估計(jì)差異。如果 Huber 和 Bisquare 的系數(shù)與 OLS 明顯不同說明異常點(diǎn)對 OLS 的拉動效應(yīng)已經(jīng)被穩(wěn)健回歸修正如果三者結(jié)果接近說明數(shù)據(jù)本身質(zhì)量較好。同時(shí)可以對比標(biāo)準(zhǔn)誤# 對比標(biāo)準(zhǔn)誤 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 參數(shù)選擇建議不同場景下的 c 值與 method 設(shè)置rlm()的參數(shù)選擇需要結(jié)合數(shù)據(jù)特征和業(yè)務(wù)需求。以下是我常用的參數(shù)設(shè)置參考表數(shù)據(jù)特征methodpsi 函數(shù)c 值理由基本干凈偶發(fā)小異常Mpsi.huber1.345保留效率只修正重尾存在若干個(gè)孤立異常值Mpsi.bisquare4.685對極端殘差直接歸零異常點(diǎn)較多或聚集成簇MMpsi.bisquare4.685高分解值抗遮蔽效應(yīng)高杠桿點(diǎn)與異常并存MMpsi.huber3.0杠桿點(diǎn)需要更漸進(jìn)地降權(quán)大樣本追求效率Mpsi.huber1.5放寬閾值減少有效樣本損失這個(gè)表的核心邏輯是異常點(diǎn)越多、越極端越傾向于使用分解值更高的估計(jì)方法和更激進(jìn)的權(quán)重函數(shù)。method MM比默認(rèn)的M估計(jì)多一個(gè)高分解值初始化步驟能有效抵抗多個(gè)異常點(diǎn)相互遮蔽的情況。5.5 避坑穩(wěn)健回歸結(jié)果解讀中的三個(gè)常見錯(cuò)誤現(xiàn)象用summary(rr.huber)中的 R2 與 OLS 的 R2 比較認(rèn)為穩(wěn)健回歸擬合效果“更好”或“更差”。原因rlm()的輸出并不包含與傳統(tǒng) OLS 直接可比的 R2。迭代加權(quán)過程中使用的權(quán)重改變了目標(biāo)函數(shù)R2 不再具有“解釋方差比例”的標(biāo)準(zhǔn)含義。解決比較模型時(shí)使用系數(shù)大小、標(biāo)準(zhǔn)誤、殘差的穩(wěn)健性和預(yù)測效果不要用 R2 作為主要判據(jù)?，F(xiàn)象把 Huber 和 Bisquare 的權(quán)重當(dāng)作樣本質(zhì)量的絕對評分權(quán)重低的樣本被認(rèn)為“一定有問題”。原因權(quán)重反映的是“在當(dāng)前模型設(shè)定下這個(gè)樣本對回歸擬合的影響相對較小”不直接等同于“這個(gè)樣本是錯(cuò)誤的”。一個(gè)在業(yè)務(wù)上重要但偏離主趨勢的樣本權(quán)重可能被壓低但它仍然包含真實(shí)信息。解決將低權(quán)重樣本單獨(dú)輸出到業(yè)務(wù)層面驗(yàn)證是否符合預(yù)期。若符合業(yè)務(wù)邏輯應(yīng)保留在數(shù)據(jù)集中甚至可以考慮單獨(dú)建?！，F(xiàn)象直接引用rr.huber$w中的權(quán)重進(jìn)行二次加權(quán)分析沒有意識到權(quán)重是在擬合后固定的。原因rlm()的權(quán)重是迭代收斂后的產(chǎn)物它們依賴于最終系數(shù)估計(jì)。換個(gè)模型設(shè)定權(quán)重會完全改變不能當(dāng)作外生變量使用。解決除非在做敏感性分析否則不要在后續(xù)分析中直接使用rlm()的權(quán)重作為通用樣本權(quán)重。如果需要穩(wěn)定的加權(quán)方案應(yīng)基于領(lǐng)域知識預(yù)先定義權(quán)重。6. 杠桿率、Cook 距離與權(quán)重的聯(lián)動驗(yàn)證一個(gè)手工計(jì)算技巧驗(yàn)證穩(wěn)健回歸是否“做對了事”有一個(gè)很實(shí)用的技巧把手動計(jì)算的杠桿率、Cook 距離與rlm()輸出的權(quán)重放到同一個(gè)數(shù)據(jù)框里用相關(guān)性檢驗(yàn)判斷降權(quán)是否準(zhǔn)確瞄準(zhǔn)了最需要降權(quán)的樣本。具體做法是計(jì)算每個(gè)樣本的 Cook 距離或者杠桿率與其在穩(wěn)健回歸中權(quán)重的 Spearman 相關(guān)如果降權(quán)邏輯正確高 Cook 距離的樣本應(yīng)該獲得低權(quán)重。這樣做的價(jià)值在于它用數(shù)據(jù)驗(yàn)證了“權(quán)重函數(shù)是否真的在折中處理極端點(diǎn)”而不是只看系數(shù)差異。具體驗(yàn)證代碼如下# 計(jì)算三個(gè)診斷量 h - hatvalues(ols) d1 - cooks.distance(ols) r - abs(stdres(ols)) # 提取兩種穩(wěn)健回歸的權(quán)重 w_huber - rr.huber$w w_bisq - rr.bisq$w # 構(gòu)建驗(yàn)證數(shù)據(jù)框 verify_df - data.frame( leverage h, cook_d d1, abs_stdres r, w_huber w_huber, w_bisq w_bisq ) # 計(jì)算 Spearman 相關(guān)系數(shù) 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) # 輸出相關(guān)系數(shù) cat(Cook距離與Huber權(quán)重的Spearman相關(guān):, cor_leverage_huber, \n) cat(Cook距離與Bisque權(quán)重的Spearman相關(guān):, cor_leverage_bisq, \n) # 找出權(quán)重最低但 Cook 距離不高的樣本檢查是否有異常降權(quán) 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 下通常表現(xiàn)出高度負(fù)相關(guān)因?yàn)?Huber 的權(quán)重直接由殘差大小決定而 Cook 距離的主要驅(qū)動因子恰恰是學(xué)生化殘差但在 Bisquare 下由于權(quán)重函數(shù)在閾值處截?cái)嘞嚓P(guān)可能變?nèi)?。這解釋了為什么 Bisquare 對極端點(diǎn)的處理更徹底對中間型異常點(diǎn)的降權(quán)卻可能更溫和。驗(yàn)證完成后把注意力放回業(yè)務(wù)層面。我通常會在輸出結(jié)果時(shí)保留三樣?xùn)|西OLS 殘差的散點(diǎn)圖、穩(wěn)健回歸權(quán)重的分布直方圖、以及按權(quán)重排序的前十個(gè)樣本的業(yè)務(wù)標(biāo)簽。這三樣配合能有效回答“為什么某個(gè)樣本被降權(quán)”以及“這個(gè)降權(quán)是否合理”。這是我自己比較習(xí)慣的一種做法。從那以后我每次做穩(wěn)健回歸都會強(qiáng)制走一遍這個(gè)流程先用plot()和cooks.distance()做診斷確認(rèn)異常點(diǎn)和高杠桿點(diǎn)的位置再用rlm()配合 Huber 或 Bisquare 權(quán)重跑一遍最后用 Spearman 相關(guān)驗(yàn)證降權(quán)邏輯是否與診斷結(jié)論一致。只有這三步全部完成我才敢把模型結(jié)果寫進(jìn)分析報(bào)告。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取