笔记说明
读《Discovering Statistics Using R》第七章 Regression中的7.7-节做的笔记。使用R进行实操部分来自节。
模型诊断概述
基于样本数据建立模型后有两个问题要考虑：
1.模型对样本数据拟合是否够好？/模型是否被少数样本点所影响？
2.模型是否能外推到其他样本？
这两个问题是有层次的——先1后2.
书中7.7.1对应第1个问题。.7.7.2对应第2个问题。
可以通过考察离群值和影响点来回答第一个问题。
离群值和残差
离群值(outlier)指与数据主体趋势有实质性不同的样本点。在5.8.1节介绍过《Discovering Statistics Using R》笔记6-箱形图和离群值
图7.9展示了线性回归中离群值对建立回归方程的影响：
离群值对回归模型的影响
离群值会影响回归模型的回归系数从而产生系统误差。图中虚线展示的是原始线性模型，实线为加入离群值样本点后拟合得到的线性模型。
残差(residual)指模型预测值与因变量实测值的差值，反映模型的误差。可以通过查看残差值很大的样本点识别离群值。
标准化残差
残差带有与因变量一致的单位，无法跨模型比较残差，也无法判定“残差很大”有统一判定标准。为克服这个困难可使用标准化残差（standardized residuals）:残差值除以残差标准误。
通过标准化，可对不同模型的残差进行比较并设立标准判断残差的可接受范围：按照正态分布样本规律，95%的Z值在-1.95，1.96之间；99%的Z值在-2.58，2.58之间；99.9%的Z值在-3.29，3.29之间。由此产生关于标准化残差的一般规则：
- 标准化残差绝对值大于3.29（可取近似值3）的样本应引起关注。
- 若有超过1%的样本其标准化残差绝对值大于2.58（可取近似值2.5），可认为模型对样本数据拟合较差。
- 若有超过5%的样本其标准化残差绝对值大于1.96（可取近似值2），可认为模型对样本数据拟合较差。
强影响点
删除某个样本后重新建模，看回归系数是否发生明显改变。这类分析可以考察回归模型是否稳定，是否有强影响点给模型带来系统误差。这个分析也能识别离群值。
评估样本对模型影响的统计量有以下几个：
- DFBeta 将某样本剔除前后分别建立模型，两模型的参数(回归系数、截距)之差为DFBeta.
- DFFit 某样本的调整预测值与原始预测值之差为DFFit。调整预测值（adjusted predicted value）是将该样本点剔除后建立新模型，使用新模型计算出的该样本点的预测值。若该样本点非强影响点，DFFit值应较小。
- 学生化残差(studentized residual) 使用调整预测值计算出的残差值除以其标准误极为学生化残差。学生化残差可在不同模型间比较，并服从t分布。
- Cook距离(Cook's distance)Cook距离衡量某样本点对模型的影响，Cook距离>1提示需要关注。
- 帽子值（hat values，也称leverage）leverage的均值定义为（k+1）/n,其中k为自变量个数，n为样本量。
leverage取值范围为0-1。leverage越接近1表示样本点对模型影响大。如果所有样本都不是强影响点，则所有样本的leverage值应该都接近（k+1）/n。
Hoaglin and Welsch (1978)建议leverage值大于2(k+1)/n的样本点视为强影响点；
Stevens (2002)建议leverage值大于3(k+1)/n的样本点视为强影响点. - Covariance ratio(CVR) 协方差比：评估样本影响回归模型参数方差的统计量。该值接近1表示对应样本对模型参数方差的影响很小。
Belsey et al.(1980)建议如果样本的CVR＞1+[3(k+1)/n],剔除该样本会减小模型参数估计的准确性；如果样本的CVR＜1-[3(k+1)/n]，则剔除该样本有助于增加模型参数估计的准确性。（k为自变量个数，n为样本量）
R实操
使用的示例数据为：Album Sales 2.dat
自变量：
- adverts：广告投入费用
- airplay：唱片发布前1周内，唱片中歌曲在广播中播放的次数
- attract：乐队的吸引程度。（打分0-10，10分表示最高）
因变量：
- sales：唱片销量
library(rio)
album2 <- import("data/Album Sales 2.dat")
str(album2)
## 'data.frame': 200 obs. of 4 variables:
## $ adverts: num 10.3 985.7 1445.6 1188.2 574.5 ...
## $ sales : int 330 120 360 270 220 170 70 210 200 300 ...
## $ airplay: int 43 28 35 33 44 19 20 22 21 40 ...
## $ attract: int 10 7 7 7 5 5 1 9 7 7 ...
模型建立：
albumSales.3 <- lm(sales ~ adverts + airplay + attract, data = album2)
上面介绍的关于离群值和强影响点的诊断统计量是对样本的，即每个样本有对应统计量。R中有很多函数计算这类统计量，一般的使用方法为：function(regressionModel)，即只需要将模型带入函数即可。
- 计算离群值相关统计量的函数：
resid（）-计算残差
rstandard()-计算标准化残差
rstudent()-计算学生化残差 - 计算强影响点相关统计量的函数：
cooks.distance()-计算Cook距离
dfbeta()-计算DFBeta
dffits()-计算DFFit
hatvalues()-计算leverage
covratio-计算协方差比
直接运行这些函数，R会在console打印很长的list，并不方便查看。建议将这些统计量存入数据集中。这里只查看前6行数据情况。
# 离群值与强影响点诊断
album2$resid <- resid(albumSales.3)
album2$stz.r<- rstandard(albumSales.3)
album2$stu.r<-rstudent(albumSales.3)
album2$cooks<-cooks.distance(albumSales.3)
album2$dfbeta<-dfbeta(albumSales.3)
album2$dffit<-dffits(albumSales.3)
album2$leverage<-hatvalues(albumSales.3)
album2$covariance.ratios<-covratio(albumSales.3)
head(round(album2, digits = 3))
## adverts sales airplay attract resid stz.r stu.r cooks dfbeta.(Intercept)
## 1 10.256 330 43 10 100.080 2.177 2.199 0.059 -5.422
## 2 985.685 120 28 7 -108.949 -2.323 -2.350 0.011 0.216
## 3 1445.563 360 35 7 68.442 1.469 1.473 0.011 -0.659
## 4 1188.193 270 33 7 7.024 0.150 0.150 0.000 -0.045
## 5 574.513 220 44 5 -5.753 -0.124 -0.123 0.000 -0.149
## 6 568.954 170 19 5 28.905 0.618 0.617 0.001 1.143
## dfbeta.adverts dfbeta.airplay dfbeta.attract dffit leverage covariance.ratios
## 1 -0.002 0.043 0.853 0.489 0.047 0.971
## 2 -0.001 0.003 -0.045 -0.211 0.008 0.920
## 3 0.001 0.013 -0.013 0.214 0.021 0.997
## 4 0.000 0.001 0.000 0.017 0.013 1.033
## 5 0.000 -0.004 0.033 -0.020 0.026 1.048
## 6 0.000 -0.006 -0.125 0.074 0.014 1.027
我们以查看标准化残差为例，根据前面的介绍可知正常情况下约有95%的样本其标准化残差绝对值≤2（更准确值为1.96）。示例数据有200个样本，则预计会有约10个样本的标准化残差绝对值大于2.
使用dplyr包的filter()筛选出标准化残差绝对值大于2的样本：
library(dplyr)
filter(album2, stz.r > 2 | stz.r < -2)
## adverts sales airplay attract resid stz.r stu.r cooks
## 1 10.256 330 43 10 100.07975 2.177404 2.198596 0.058703882
## 2 985.685 120 28 7 -108.94899 -2.323083 -2.349724 0.010889432
## 3 174.093 300 40 7 99.53375 2.130289 2.149882 0.017756472
## 4 102.568 40 25 8 -114.96982 -2.460996 -2.493538 0.024115188
## 5 405.913 190 12 4 97.40266 2.099446 2.118034 0.033159177
## 6 1542.329 190 33 8 -114.12308 -2.455913 -2.488224 0.040415897
## 7 579.321 300 30 7 98.81030 2.104079 2.122816 0.005948358
## 8 56.895 70 37 7 -110.41564 -2.363549 -2.391845 0.022288983
## 9 1000.000 250 5 7 97.28666 2.095399 2.113858 0.031364021
## 10 9.104 120 53 8 -121.32405 -2.628814 -2.669584 0.070765882
## 11 145.585 360 42 8 144.13246 3.093333 3.163622 0.050867000
## 12 785.694 110 20 9 -97.20606 -2.088044 -2.106269 0.025134553
## dfbeta.(Intercept) dfbeta.adverts dfbeta.airplay dfbeta.attract dffit
## 1 -5.4218270671 -0.0016615915 0.0433929166 0.8529699235 0.4892940
## 2 0.2160170176 -0.0008649690 0.0025870806 -0.0450304095 -0.2110983
## 3 -0.2159709599 -0.0010709506 0.0461632913 0.0162397840 0.2689580
## 4 1.1378163541 0.0013393286 0.0132378865 -0.4296547666 -0.3146882
## 5 6.0692407906 -0.0001976727 -0.0376293901 -0.6515969602 0.3674177
## 6 2.9843774733 -0.0022309557 -0.0063243216 -0.2992085381 -0.4073640
## 7 0.0140823505 -0.0001055773 0.0076884972 0.0496393949 0.1556248
## 8 -0.0481327226 0.0014466228 -0.0405291895 -0.0423968247 -0.3021645
## 9 1.0513491554 0.0009966248 -0.0825582156 0.1635148508 0.3573180
## 10 3.0723591414 0.0019761696 -0.1096533563 -0.2810254509 -0.5402885
## 11 -2.8531867723 -0.0017440692 0.0699075972 0.4044744026 0.4613239
## 12 2.8608326140 -0.0003183919 0.0391384577 -0.6261084622 -0.3198451
## leverage covariance.ratios
## 1 0.047190526 0.9712750
## 2 0.008006536 0.9201832
## 3 0.015409738 0.9439200
## 4 0.015677123 0.9145800
## 5 0.029213132 0.9599533
## 6 0.026103520 0.9248580
## 7 0.005345708 0.9365377
## 8 0.015708852 0.9236983
## 9 0.027779409 0.9588774
## 10 0.039348661 0.9203731
## 11 0.020821154 0.8532470
## 12 0.022539842 0.9543502
示例数据中有12个样本（6%）的标准化残差值绝对值大于2。拟合情况较好。
网友评论