最小二乘法（OLS）是很常用的线性回归。
本文介绍的IRLS是其变化版。
对数据中异常值的处理会有很大提升。
简单搜了一下，网上对该方法还没有中文的说明，也可能是我没有找到。
几个基本概念：
Residual：残差，预测值(基于回归方程)与实际观测值之间的差值。 Outlier：在线性回归中，离群值是具有较大残差的观测值。 Leverage：在预测变量上具有极值的观测值是具有高杠杆的点。杠杆是衡量一个自变量偏离其均值的程度。高杠杆点对回归系数的估计有很大的影响。 Influence：如果移除观测结果会使回归系数的估计发生很大的变化，那么该观测结果就是有影响的。影响力可以被认为是杠杆和离群值的产物。 Cook’s distance：测量杠杆信息和残差的方法。
关于IRLS： rlm属于稳健回归（Robust regression）的一个方法。 稳健回归可以用在任何使用最小二乘回归的情况下。在拟合最小二乘回归时，我们可能会发现一些异常值或高杠杆数据点。已经确定这些数据点不是数据输入错误，也不是来自另一个群落。所以我们没有令人信服的理由将它们排除在分析之外。 稳健回归可能是一种好的策略，它是在将这些点完全从分析中排除；和包括所有数据点；以及在OLS回归中平等对待所有数据点之间的妥协。他可以个给每个样本一个权重，离群值权重低一些，正常值权重高一些，进行校正。 rlm可在MASS包实现。
#搞一个数据
cdata <- read.dta("https://stats.idre.ucla.edu/stat/data/crime.dta")
#先用OLS试试
ols <- lm(crime ~ poverty single, data = cdata)
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(ols, las = 1)
#从残差结果可知，9, 25, 51 是异常值。
#然后计算Cook’s distance.一般将高于4/n的值为异常高的值。
d1 <- cooks.distance(ols)
r <- stdres(ols)
a <- cbind(cdata, d1, r)
a[d1 > 4/51, ]
sid state crime murder pctmetro pctwhite pcths poverty single d1 r
1 1 ak 761 9.0 41.8 75.2 86.6 9.1 14.3 0.1254750 -1.397418
9 9 fl 1206 8.9 93.0 83.5 74.4 17.8 10.6 0.1425891 2.902663
25 25 ms 434 13.5 30.7 63.3 64.3 24.7 14.7 0.6138721 -3.562990
51 51 dc 2922 78.5 100.0 31.8 73.1 26.4 22.1 2.6362519 2.616447
#可以看到4个样本的cook值很高，高于我们设的阈值。
#接下来用rlm试试~
#默认的权重算法为Huber方法~
rr.huber <- rlm(crime ~ poverty single, data = cdata)
#将权重和残差排个序输出
hweights <- data.frame(state = cdata$state, resid = rr.huber$resid, weight = rr.huber$w)
hweights2 <- hweights[order(rr.huber$w), ]
hweights2[1:15, ]
state resid weight
25 ms -846.08536 0.2889618
9 fl 679.94327 0.3595480
46 vt -410.48310 0.5955740
51 dc 376.34468 0.6494131
26 mt -356.13760 0.6864625
21 me -337.09622 0.7252263
31 nj 331.11603 0.7383578
14 il 319.10036 0.7661169
1 ak -313.15532 0.7807432
20 md 307.19142 0.7958154
19 ma 291.20817 0.8395172
18 la -266.95752 0.9159411
2 al 105.40319 1.0000000
#明显的看到，残差越高的样本权重越低。
#然后再换一种权重算法bisquare
rr.bisquare <- rlm(crime ~ poverty single, data=cdata, psi = psi.bisquare)
biweights <- data.frame(state = cdata$state, resid = rr.bisquare$resid, weight = rr.bisquare$w)
biweights2 <- biweights[order(rr.bisquare$w), ]
biweights2[1:15, ]
state resid weight
25 ms -905.5931 0.007652565
9 fl 668.3844 0.252870542
46 vt -402.8031 0.671495418
26 mt -360.8997 0.731136908
31 nj 345.9780 0.751347695
18 la -332.6527 0.768938330
21 me -328.6143 0.774103322
1 ak -325.8519 0.777662383
14 il 313.1466 0.793658594
20 md 308.7737 0.799065530
19 ma 297.6068 0.812596833
51 dc 260.6489 0.854441716
50 wy -234.1952 0.881660897
5 ca 201.4407 0.911713981
10 ga -186.5799 0.924033113
#这个权重和上面的方法差别就很大了
综上，rlm是比OLS更好的方法。
但是巨大的差异表明模型参数受到异常值的高度影响。 不同的权重算法各有优点和缺点。Huber可能会难以处理严重的异常值，而bisquare可能会难以收敛或产生多个解决方案。
Reference: https://stats.idre.ucla.edu/r/dae/robust-regression/
