先容几个线性回归（linear regression）中的术语：
残差（Residual）: 基于回归方程的预测值与视察值的差。
离群点（Outlier）: 线性回归（linear regression）中的离群点是指对应残差较大的视察值。也就是说，当某个视察值与基于回归方程的预测值相差较大时，该视察值即可视为离群点。 离群点的呈现一般是因为样本自身较为非凡可能数据录入错误导致的，虽然也大概是其他问题。
杠杆率（Leverage）: 当某个视察值所对应的预测值为极度值时，该视察值称为高杠杆率点。杠杆率权衡的是独立变量对自身均值的偏异水平。高杠杆率的视察值对付回归方程的参数有重大影响。
影响力点：（Influence）: 若某视察值的剔除与否，对回归方程的系数预计有显著相应，则该视察值是具有影响力的，称为影响力点。影响力是高杠杆率和离群环境引起的。
Cook间隔（Cook’s distance）: 综合了杠杆率信息和残差信息的统计量。
利用最小二乘回归时，有时候会碰着离群点和高杠杆率点。此时，若认定离群点可能高杠杆率点的呈现并非因为数据录入错误可能该该视察值来自别的一个总体的话，利用最小二乘回归会变得很棘手，因为数据阐明者因为没有充实的来由剔除离群点和高杠杆率。此时稳健回归是个极佳的替代方案。稳健回归在剔除离群点可能高杠杆率点和保存离群点或高杠杆率点并像最小二乘法那样平等利用各点之间找到了一个折中。其在预计回归参数时，按照视察值的稳健环境对视察值举办赋权。简而言之，稳健回归是加权最小二乘回归，或称文艺最小二乘回归。
MASS 包中的 rlm呼吁提供了差异形式的稳健回归拟合方法。接下来，以基于Huber要领和bisquare要领下的M预计为例来举办演示。这是两种更为根基的M预计要领。在M预计中，要做的工作是在满意约束 ∑ i=1 n wi (yi-x’b) xi’ = 0 时，求出使得 ∑wi2ei2 最小的参数。由于权重的预计依赖于残差，而残差的预计又反过来依赖于权重，因此需用迭代反复加权最小二乘（ Iteratively Reweighted Least Squares ，IRLS）来预计参数。举例，第j次迭代获得的系数矩阵为: Bj = [ X’ Wj-1 X ] -1 X’ Wj-1 Y ，这里下脚标暗示求解进程中的迭代次数，而不是凡是的行标可能列标，一连这一进程，直到功效收敛为止。Huber要领下，残差较小的视察值被赋予的权重为1，残差较大的视察值的权重跟着残差的增大而递减，详细函数为：w(e)={1 for |e|<=kk|e| for |e|>k . 而bisquare要领下，所有的非0残差所对应视察值的权重都是递减的。
数据描写
下面用到的数据是Alan Agresti和Barbara Finlay所著的《Statistical Methods for Social Sciences》(Third Edition，Prentice Hall, 1997)中的crime数据集。该数据集的变量别离是：
state id (sid),
state name (state),
violent crimes per 100,000 people (crime),
murders per 1,000,000 (murder),
the percent of the population living in metropolitan areas (pctmetro),
the percent of the population that is white (pctwhite),
percent of population with a high school education or above (pcths),
percent of population living under poverty line (poverty),
percent of population that are single parents (single).
该数据集共有51个视察值。接下来用数据会合的poverty 和 single 变量来预测 crime.
library(foreign) cdata <- read.dta(“http://www.ats.ucla.edu/stat/data/crime.dta”) summary(cdata)
sid state crime murder
Min. : 1.0 Length:51 Min. : 82 Min. : 1.60
1st Qu.:13.5 Class :character 1st Qu.: 326 1st Qu.: 3.90
Median :26.0 Mode :character Median : 515 Median : 6.80
Mean :26.0 Mean : 613 Mean : 8.73
3rd Qu.:38.5 3rd Qu.: 773 3rd Qu.:10.35
Max. :51.0 Max. :2922 Max. :78.50
pctmetro pctwhite pcths poverty single
Min. : 24.0 Min. :31.8 Min. :64.3 Min. : 8.0 Min. : 8.4
1st Qu.: 49.5 1st Qu.:79.4 1st Qu.:73.5 1st Qu.:10.7 1st Qu.:10.1
Median : 69.8 Median :87.6 Median :76.7 Median :13.1 Median :10.9
Mean : 67.4 Mean :84.1 Mean :76.2 Mean :14.3 Mean :11.3
3rd Qu.: 84.0 3rd Qu.:92.6 3rd Qu.:80.1 3rd Qu.:17.4 3rd Qu.:12.1
Max. :100.0 Max. :98.5 Max. :86.6 Max. :26.4 Max. :22.1
稳健回归
先对数据举办OLS回归，重点调查回归功效中的残差、拟合值、Cook间隔和杠杆率。
ols <- lm(crime ~ poverty + single, data = cdata) summary(ols)
Call:
lm(formula = crime ~ poverty + single, data = cdata)
Residuals:
Min 1Q Median 3Q Max
-811.1 -114.3 -22.4 121.9 689.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1368.19 187.21 -7.31 2.5e-09 ***
poverty 6.79 8.99 0.76 0.45
single 166.37 19.42 8.57 3.1e-11 ***
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
Residual standard error: 244 on 48 degrees of freedom
Multiple R-squared: 0.707, Adjusted R-squared: 0.695
F-statistic: 58 on 2 and 48 DF, p-value: 1.58e-13
opar <- par(mfrow = c(2,2),oma = c(0, 0, 1.1, 0)) plot(ols, las = 1)
从图上看出，第 9, 第25和第5个视察值大概是离群点，看看这些视察值所属的是哪些州。
cdata[c(9, 25, 51), 1:2]
sid state
9 9 fl
25 25 ms
51 51 dc
可以揣摩，DC, Florida 和Mississippi这三个处所所对应的视察值大概具有较大的残差可能杠杆率。下面调查一下Cook间隔较大的视察值有哪些。在判定Cook间隔巨细的时候，凡是回收过的履历分界点是Cook间隔序列的4/n处，个中n是视察值的个数。
library(MASS) 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.1255 -1.397
9 9 fl 1206 8.9 93.0 83.5 74.4 17.8 10.6 0.1426 2.903
25 25 ms 434 13.5 30.7 63.3 64.3 24.7 14.7 0.6139 -3.563
51 51 dc 2922 78.5 100.0 31.8 73.1 26.4 22.1 2.6363 2.616
原来该当先删除DC所对应的视察值，因为DC对应的并不是州。然而，由于DC所对应的Cook间隔较大保存DC有助于我们举办调查。下面生成一个absr1变量, 其对应的为残差序列的值，取出残差值较大的视察值：
rabs <- abs(r) a <- cbind(cdata, d1, r, rabs) asorted <- a[order(-rabs), ] asorted[1:10, ]
sid state crime murder pctmetro pctwhite pcths poverty single d1 r rabs
25 25 ms 434 13.5 30.7 63.3 64.3 24.7 14.7 0.61387 -3.563 3.563
9 9 fl 1206 8.9 93.0 83.5 74.4 17.8 10.6 0.14259 2.903 2.903
51 51 dc 2922 78.5 100.0 31.8 73.1 26.4 22.1 2.63625 2.616 2.616
46 46 vt 114 3.6 27.0 98.4 80.8 10.0 11.0 0.04272 -1.742 1.742
26 26 mt 178 3.0 24.0 92.6 81.0 14.9 10.8 0.01676 -1.461 1.461
21 21 me 126 1.6 35.7 98.5 78.8 10.7 10.6 0.02233 -1.427 1.427
1 1 ak 761 9.0 41.8 75.2 86.6 9.1 14.3 0.12548 -1.397 1.397
31 31 nj 627 5.3 100.0 80.8 76.7 10.9 9.6 0.02229 1.354 1.354
14 14 il 960 11.4 84.0 81.0 76.2 13.6 11.5 0.01266 1.338 1.338
20 20 md 998 12.7 92.8 68.9 78.4 9.7 12.0 0.03570 1.287 1.287
此刻转向稳健回归。再提示一下，稳健回归是通过迭代反复加权最小二乘（iterated re-weighted least squares ，IRLS)来完成的。其对应的R函数是MASS包中的rlm。IRLS对应的有多个权重函数（ weighting functions），首先演示一下Huber要领。 演示进程中，重点存眷IRLS进程得出的权重功效。
rr.huber <- rlm(crime ~ poverty + single, data = cdata) summary(rr.huber)
Call: rlm(formula = crime ~ poverty + single, data = cdata)
Residuals:
Min 1Q Median 3Q Max
-846.1 -125.8 -16.5 119.2 679.9
Coefficients:
Value Std. Error t value
(Intercept) -1423.037 167.590 -8.491
poverty 8.868 8.047 1.102
single 168.986 17.388 9.719
Residual standard error: 182 on 48 degrees of freedom
hweights <- data.frame(state = cdata$state, resid = rr.huber$resid, weight = rr.huber$w) hweights2 <- hweights[order(rr.huber$w), ] hweights2[1:15, ]
25 ms -846.09 0.2890
9 fl 679.94 0.3595
46 vt -410.48 0.5956
51 dc 376.34 0.6494
26 mt -356.14 0.6865
21 me -337.10 0.7252
31 nj 331.12 0.7384
14 il 319.10 0.7661
1 ak -313.16 0.7807
20 md 307.19 0.7958
19 ma 291.21 0.8395
18 la -266.96 0.9159
2 al 105.40 1.0000
3 ar 30.54 1.0000
4 az -43.25 1.0000
容易看出来，视察值的残差值越大，其被赋予的权重越小。功效表白Mississippi所对应的视察值被赋予的权重是最小的，其次是Florida所对应的视察值，而所有未被展示的视察值的权重皆为1。由于OLS回归中所有视察值的权重都为1，因此，稳健回归中权重为1的视察值越多，则稳健回归于OLS回归的阐明功效越临近。
接下来，用bisquare要领来举办稳健回归进程。
rr.bisquare <- rlm(crime ~ poverty + single, data=cdata, psi = psi.bisquare) summary(rr.bisquare)
Call: rlm(formula = crime ~ poverty + single, data = cdata, psi = psi.bisquare)
Residuals:
Min 1Q Median 3Q Max
-906 -141 -15 115 668
Coefficients:
Value Std. Error t value
(Intercept) -1535.334 164.506 -9.333
poverty 11.690 7.899 1.480
single 175.930 17.068 10.308
Residual standard error: 202 on 48 degrees of freedom 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.6 0.007653
9 fl 668.4 0.252871
46 vt -402.8 0.671495
26 mt -360.9 0.731137
31 nj 346.0 0.751348
18 la -332.7 0.768938
21 me -328.6 0.774103
1 ak -325.9 0.777662
14 il 313.1 0.793659
20 md 308.8 0.799066
19 ma 297.6 0.812597
51 dc 260.6 0.854442
50 wy -234.2 0.881661
5 ca 201.4 0.911714
10 ga -186.6 0.924033
与Huber要领对比，bisquare要领下的 Mississippi视察值被赋予了极小的权重，而且两种要领预计出的回归参数也相差甚大。凡是，当稳健回归跟OLS回归的阐明功效相差较大时，数据阐明者回收稳健回归较为明智。稳健回归和OLS回归的阐明功效的较大差别凡是体现着离群点对模子参数发生了较大影响。所有的要领都有优点和软肋，稳健回归也不破例。稳健回归中，Huber要领的软肋在于无法很好的而处理惩罚极度离群点，而bisquare要领的软肋在于回归功效不易收敛，以至于常常有多个最优解。
除此之外，两种要领得出的参数功效极为差异，尤其是single变量的系数和截距项(intercept)。不外，一般而言无需存眷截距项，除非事先已经对预测变量举办了中心化，此时截距项才显的有些用处。再有， 变量 poverty的系数在两种要领下都不显著，而变量 single则恰好相反，都较为显著。
