最小二乘法（OLS）是很常用的線性回歸。
本文介紹的IRLS是其變化版。
對數據中異常值的處理會有很大提升。
簡單搜了一下，網上對該方法還沒有中文的說明，也可能是我沒有找到。
幾個基本概念：
Residual：殘差，預測值(基於回歸方程)與實際觀測值之間的差值。
Outlier：在線性回歸中，離群值是具有較大殘差的觀測值。
Leverage：在預測變量上具有極值的觀測值是具有高槓桿的點。槓桿是衡量一個自變量偏離其均值的程度。高槓桿點對回歸係數的估計有很大的影響。
Influence：如果移除觀測結果會使回歸係數的估計發生很大的變化，那麼該觀測結果就是有影響的。影響力可以被認為是槓桿和離群值的產物。
Cook’s distance：測量槓桿信息和殘差的方法。
關於IRLS：
rlm屬於穩健回歸（Robust regression）的一個方法。
穩健回歸可以用在任何使用最小二乘回歸的情況下。在擬合最小二乘回歸時，我們可能會發現一些異常值或高槓桿數據點。已經確定這些數據點不是數據輸入錯誤，也不是來自另一個群落。所以我們沒有令人信服的理由將它們排除在分析之外。
穩健回歸可能是一種好的策略，它是在將這些點完全從分析中排除；和包括所有數據點；以及在OLS回歸中平等對待所有數據點之間的妥協。他可以個給每個樣本一個權重，離群值權重低一些，正常值權重高一些，進行校正。
rlm可在MASS包實現。
1#搞一個數據
2cdata <- read.dta("https://stats.idre.ucla.edu/stat/data/crime.dta")
3
4#先用OLS試試
5ols <- lm(crime ~ poverty + single, data = cdata)
6opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
7plot(ols, las = 1)
1#從殘差結果可知，9, 25, 51 是異常值。
2#然後計算Cook’s distance.一般將高於4/n的值為異常高的值。
3d1 <- cooks.distance(ols)
4r <- stdres(ols)
5a <- cbind(cdata, d1, r)
6a[d1 > 4/51, ]
7
8 sid state crime murder pctmetro pctwhite pcths poverty single d1 r
91 1 ak 761 9.0 41.8 75.2 86.6 9.1 14.3 0.1254750 -1.397418
109 9 fl 1206 8.9 93.0 83.5 74.4 17.8 10.6 0.1425891 2.902663
1125 25 ms 434 13.5 30.7 63.3 64.3 24.7 14.7 0.6138721 -3.562990
1251 51 dc 2922 78.5 100.0 31.8 73.1 26.4 22.1 2.6362519 2.616447
13
14#可以看到4個樣本的cook值很高，高於我們設的閾值。
15
16#接下來用rlm試試~
17#默認的權重算法為Huber方法~
18rr.huber <- rlm(crime ~ poverty + single, data = cdata)
19#將權重和殘差排個序輸出
20hweights <- data.frame(state = cdata$state, resid = rr.huber$resid, weight = rr.huber$w)
21hweights2 <- hweights[order(rr.huber$w), ]
22hweights2[1:15, ]
23
24 state resid weight
2525 ms -846.08536 0.2889618
269 fl 679.94327 0.3595480
2746 vt -410.48310 0.5955740
2851 dc 376.34468 0.6494131
2926 mt -356.13760 0.6864625
3021 me -337.09622 0.7252263
3131 nj 331.11603 0.7383578
3214 il 319.10036 0.7661169
331 ak -313.15532 0.7807432
3420 md 307.19142 0.7958154
3519 ma 291.20817 0.8395172
3618 la -266.95752 0.9159411
372 al 105.40319 1.0000000
38
39#明顯的看到，殘差越高的樣本權重越低。
40
41#然後再換一種權重算法bisquare
42rr.bisquare <- rlm(crime ~ poverty + single, data=cdata, psi = psi.bisquare)
43biweights <- data.frame(state = cdata$state, resid = rr.bisquare$resid, weight = rr.bisquare$w)
44biweights2 <- biweights[order(rr.bisquare$w), ]
45biweights2[1:15, ]
46
47 state resid weight
4825 ms -905.5931 0.007652565
499 fl 668.3844 0.252870542
5046 vt -402.8031 0.671495418
5126 mt -360.8997 0.731136908
5231 nj 345.9780 0.751347695
5318 la -332.6527 0.768938330
5421 me -328.6143 0.774103322
551 ak -325.8519 0.777662383
5614 il 313.1466 0.793658594
5720 md 308.7737 0.799065530
5819 ma 297.6068 0.812596833
5951 dc 260.6489 0.854441716
6050 wy -234.1952 0.881660897
615 ca 201.4407 0.911713981
6210 ga -186.5799 0.924033113
63
64#這個權重和上面的方法差別就很大了
綜上，rlm是比OLS更好的方法。
但是巨大的差異表明模型參數受到異常值的高度影響。
不同的權重算法各有優點和缺點。Huber可能會難以處理嚴重的異常值，而bisquare可能會難以收斂或產生多個解決方案。
Reference:
https://stats.idre.ucla.edu/r/dae/robust-regression/