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