参考：
1、https://blog.csdn.net/daunxx/article/details/51858208
2、https://stats.idre.ucla.edu/r/dae/robust-regression/
1、数据准备与数据理解
数据源：https://stats.idre.ucla.edu/stat/data/crime.dta
> library(pacman)
> p_load(foreign,dplyr)
> crime <- read.dta("./crime.dta")
> str(crime)
## 'data.frame': 51 obs. of 9 variables:
## $ sid : num 1 2 3 4 5 6 7 8 9 10 ...
## $ state : chr "ak" "al" "ar" "az" ...
## $ crime : int 761 780 593 715 1078 567 456 686 1206 723 ...
## $ murder : num 9 11.6 10.2 8.6 13.1 ...
## $ pctmetro: num 41.8 67.4 44.7 84.7 96.7 ...
## $ pctwhite: num 75.2 73.5 82.9 88.6 79.3 ...
## $ pcths : num 86.6 66.9 66.3 78.7 76.2 ...
## $ poverty : num 9.1 17.4 20 15.4 18.2 ...
## $ single : num 14.3 11.5 10.7 12.1 12.5 ...
## - attr(*, "datalabel")= chr "crime data from agresti & finlay - 1997"
## - attr(*, "time.stamp")= chr "24 Feb 2001 14:23"
## - attr(*, "formats")= chr [1:9] "%9.0g" "%9s" "%8.0g" "%9.0g" ...
## - attr(*, "types")= int [1:9] 102 130 105 102 102 102 102 102 102
## - attr(*, "val.labels")= chr [1:9] "" "" "" "" ...
## - attr(*, "var.labels")= chr [1:9] "" "" "violent crime rate" "murder rate" ...
## - attr(*, "version")= int 6
数据集共51行，9个变量。查看每个变量的注释：
> attr(crime,"var.labels")
## [1] "" "" "violent crime rate"
## [4] "murder rate" "pct metropolitan" "pct white"
## [7] "pct hs graduates" "pct poverty" "pct single parent"
sid：州ID
state：州名
crime：暴力犯罪率
murder：谋杀率
pctmetro：生活在大都市地区的人口比例
pctwhite：白人人口比例
pcths：高中及以上学历人口比例
poverty：生活在贫困线以下的人口比例
single：单亲家庭比例
2、最小二乘回归
我们只研究暴力犯罪率与贫困线以下人口比例和单亲家庭之间的关系。
首先使用OLS回归并进行一些诊断，作为稳健回归的参考依据。
> fit.ols <- lm(crime ~ poverty + single,data = crime)
> summary(fit.ols)
##
## Call:
## lm(formula = crime ~ poverty + single, data = crime)
##
## Residuals:
## Min 1Q Median 3Q Max
## -811.14 -114.27 -22.44 121.86 689.82
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1368.189 187.205 -7.308 2.48e-09 ***
## poverty 6.787 8.989 0.755 0.454
## single 166.373 19.423 8.566 3.12e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 243.6 on 48 degrees of freedom
## Multiple R-squared: 0.7072, Adjusted R-squared: 0.695
## F-statistic: 57.96 on 2 and 48 DF, p-value: 1.578e-13
> par(mfrow=c(2,2))
> plot(fit.ols)
残差（Residual）: 基于回归方程的预测值与观测值的差。
离群点（Outlier）: 线性回归（linear regression）中的离群点是指对应残差较大的观测值。也就是说，当某个观测值与基于回归方程的预测值相差较大时，该观测值即可视为离群点。 离群点的出现一般是因为样本自身较为特殊或者数据录入错误导致的，当然也可能是其他问题。
杠杆率（Leverage）: 当某个观测值所对应的预测值为极端值时，该观测值称为高杠杆率点。杠杆率衡量的是独立变量对自身均值的偏离程度。高杠杆率的观测值对于回归方程的参数有重大影响。
影响力点：（Influence）: 若某观测值的剔除与否，对回归方程的系数估计有显著影响，则该观测值是具有影响力的，称为影响力点。影响力是高杠杆率和离群情况引起的。
Cook距离（Cook's distance）: 综合了杠杆率信息和残差信息的统计量。
从上面的图中我们可以看到9、25、51在四个图中都出现了，说明它们对我们的模型可能是有问题的。我们已经确定这些数据点不是数据输入错误，也不是来自于与其他大多数数据不同的人群，所以也没有足够理由剔除它们。
先观察下这些观测来自于哪些州：
> crime[c(9,25,51),c(1,2)]
## sid state
## 9 9 fl
## 25 25 ms
## 51 51 dc
fl、ms、dc州要么有高杠杆率，要么有高残差。我们再看看Cook's distance有相对较大值的观测值：
> # 计算cook距离
> d <- cooks.distance(fit.ols)
> # 残差
> r <- resid(fit.ols)
> # 合并到原始数据
> df <- cbind(crime,d,r)
> # cook距离的cut-off点一般是4/n
> df[d>4/nrow(crime),]
## sid state crime murder pctmetro pctwhite pcths poverty single d
## 1 1 ak 761 9.0 41.8 75.2 86.6 9.1 14.3 0.1254750
## 9 9 fl 1206 8.9 93.0 83.5 74.4 17.8 10.6 0.1425891
## 25 25 ms 434 13.5 30.7 63.3 64.3 24.7 14.7 0.6138721
## 51 51 dc 2922 78.5 100.0 31.8 73.1 26.4 22.1 2.6362519
## r
## 1 -311.7055
## 9 689.8233
## 25 -811.1373
## 51 434.1663
dc州的cook距离是较大的。现在我们看看残差，因为残差的符号是不重要的，我们新生成一列为残差的绝对值，并按大小排序：
> df <- df %>% mutate(rabs=abs(r)) %>% arrange(-rabs)
> head(df,10)
## sid state crime murder pctmetro pctwhite pcths poverty single d
## 1 25 ms 434 13.5 30.7 63.3 64.3 24.7 14.7 0.61387212
## 2 9 fl 1206 8.9 93.0 83.5 74.4 17.8 10.6 0.14258909
## 3 51 dc 2922 78.5 100.0 31.8 73.1 26.4 22.1 2.63625193
## 4 46 vt 114 3.6 27.0 98.4 80.8 10.0 11.0 0.04271548
## 5 26 mt 178 3.0 24.0 92.6 81.0 14.9 10.8 0.01675501
## 6 21 me 126 1.6 35.7 98.5 78.8 10.7 10.6 0.02233128
## 7 31 nj 627 5.3 100.0 80.8 76.7 10.9 9.6 0.02229184
## 8 14 il 960 11.4 84.0 81.0 76.2 13.6 11.5 0.01265689
## 9 1 ak 761 9.0 41.8 75.2 86.6 9.1 14.3 0.12547500
## 10 20 md 998 12.7 92.8 68.9 78.4 9.7 12.0 0.03569623
## r rabs
## 1 -811.1373 811.1373
## 2 689.8233 689.8233
## 3 434.1663 434.1663
## 4 -415.7843 415.7843
## 5 -351.7679 351.7679
## 6 -341.9864 341.9864
## 7 324.0288 324.0288
## 8 322.5949 322.5949
## 9 -311.7055 311.7055
## 10 303.8792 303.8792
总之，我们没有充分的理由剔除这些离群点和高杠杆率点，但是它们又对模型产生了较大的影响，如下图：
数据中存在一个异常点，如果不剔除该点，适用OLS方法来做回归的话，那么就会得到图中红色的那条线；如果将这个异常点剔除掉的话，那么就可以得到图中蓝色的那条线。显然，蓝色的线比红色的线对数据有更强的解释性，这就是OLS在做回归分析时候的弊端。
3、稳健回归
稳健回归（Robust regression），就是当最小二乘法遇到上述的异常点的时候，用于代替最小二乘法的一个算法。当然，稳健回归还可以用于异常点检测，或者是找出那些对模型影响最大的样本点。
稳健回归在剔除离群点或者高杠杆率点和保留离群点或高杠杆率点并像最小二乘法那样平等使用各点之间找到了一个折中方案。其在估计回归参数时，根据观测值的稳健情况对观测值进行赋权。简而言之，稳健回归是加权最小二乘回归。
MASS 包中的 rlm命令提供了不同形式的稳健回归拟合方式。接下来，以基于Huber方法和bisquare方法下的M估计为例来进行演示。
Huber方法下，残差较小的观测值被赋予的权重为1，残差较大的观测值的权重随着残差的增大而递减。而bisquare方法下，所有的非0残差所对应观测值的权重都是递减的。
首先演示一下Huber方法。演示过程中，重点关注IRLS过程得出的权重结果。
> p_load(MASS)
> fit.huber <- rlm(crime ~ poverty + single,data = crime)
> summary(fit.huber)
##
## Call: rlm(formula = crime ~ poverty + single, data = crime)
## Residuals:
## Min 1Q Median 3Q Max
## -846.09 -125.80 -16.49 119.15 679.94
##
## Coefficients:
## Value Std. Error t value
## (Intercept) -1423.0373 167.5899 -8.4912
## poverty 8.8677 8.0467 1.1020
## single 168.9858 17.3878 9.7186
##
## Residual standard error: 181.8 on 48 degrees of freedom
> hweight <- data.frame(state=crime$state,
+ resid=fit.huber$residuals,
+ weight=fit.huber$w)
> hweight <- arrange(hweight,weight)
> head(hweight,15)
## state resid weight
## 1 ms -846.08536 0.2889618
## 2 fl 679.94327 0.3595480
## 3 vt -410.48310 0.5955740
## 4 dc 376.34468 0.6494131
## 5 mt -356.13760 0.6864625
## 6 me -337.09622 0.7252263
## 7 nj 331.11603 0.7383578
## 8 il 319.10036 0.7661169
## 9 ak -313.15532 0.7807432
## 10 md 307.19142 0.7958154
## 11 ma 291.20817 0.8395172
## 12 la -266.95752 0.9159411
## 13 al 105.40319 1.0000000
## 14 ar 30.53589 1.0000000
## 15 az -43.25299 1.0000000
从结果中可以粗略的看到，当绝对残差下降时，权重增加。在OLS回归中，所有观测的权重都是1。因此，在稳健回归中权重接近于1的情况越多，OLS和稳健回归的结果就越接近。
> p_load(ggplot2)
> ggplot(hweight,aes(abs(resid),weight)) +
+ geom_point(size=1) +
+ theme_bw()
接下来使用bisquare方法：
> fit.bisquare <- rlm(crime ~ poverty + single,
+ data = crime,psi=psi.bisquare)
> summary(fit.bisquare)
##
## Call: rlm(formula = crime ~ poverty + single, data = crime, psi = psi.bisquare)
## Residuals:
## Min 1Q Median 3Q Max
## -905.59 -140.97 -14.98 114.65 668.38
##
## Coefficients:
## Value Std. Error t value
## (Intercept) -1535.3338 164.5062 -9.3330
## poverty 11.6903 7.8987 1.4800
## single 175.9303 17.0678 10.3077
##
## Residual standard error: 202.3 on 48 degrees of freedom
> biweight <- data.frame(state=crime$state,
+ resid=fit.bisquare$residuals,
+ weight=fit.bisquare$w)
> biweight <- arrange(biweight,weight)
> head(biweight,15)
## state resid weight
## 1 ms -905.5931 0.007652565
## 2 fl 668.3844 0.252870542
## 3 vt -402.8031 0.671495418
## 4 mt -360.8997 0.731136908
## 5 nj 345.9780 0.751347695
## 6 la -332.6527 0.768938330
## 7 me -328.6143 0.774103322
## 8 ak -325.8519 0.777662383
## 9 il 313.1466 0.793658594
## 10 md 308.7737 0.799065530
## 11 ma 297.6068 0.812596833
## 12 dc 260.6489 0.854441716
## 13 wy -234.1952 0.881660897
## 14 ca 201.4407 0.911713981
## 15 ga -186.5799 0.924033113
> ggplot(bihweight,aes(abs(resid),weight)) +
+ geom_point(size=1) +
+ theme_bw()
> coef <- data.frame(OLS=fit.ols$coefficients,
> HUBER=fit.huber$coefficients,
> BISQUARE=fit.bisquare$coefficients)
> coef
## OLS HUBER BISQUARE
## (Intercept) -1368.188661 -1423.037337 -1535.33376
## poverty 6.787359 8.867678 11.69033
## single 166.372670 168.985787 175.93032
可以看到，ms州的权重显著低于Huber方法中计算的权重，并且三种方法估计出的回归参数也相差甚大。通常，当稳健回归跟OLS回归的分析结果相差较大时，数据分析者采用稳健回归较为明智。稳健回归和OLS回归的分析结果的较大差异通常暗示着离群点对模型参数产生了较大影响。
所有的回归方法都有自己的长处和缺陷，稳健回归也不例外。稳健回归中，Huber方法的软肋在于无法很好的处理极端离群点，而bisquare方法的软肋在于回归结果不易收敛，以至于经常有多个最优解。
除此之外，两种方法得出的参数结果也极为不同，尤其是single变量的系数和截距项(intercept)。不过，一般而言无需关注截距项，除非事先已经对预测变量进行了中心化，此时截距项才显的有些用处。再有，变量poverty的系数在两种方法下都不显著，而变量single则刚好相反，都较为显著。
