筆記說(shuō)明
讀《Discovering Statistics Using R》第七章 Regression中的7.7-節(jié)做的筆記。使用R進(jìn)行實(shí)操部分來(lái)自節(jié)。
模型診斷概述
基于樣本數(shù)據(jù)建立模型后有兩個(gè)問(wèn)題要考慮：
1.模型對(duì)樣本數(shù)據(jù)擬合是否夠好？/模型是否被少數(shù)樣本點(diǎn)所影響？
2.模型是否能外推到其他樣本？
這兩個(gè)問(wèn)題是有層次的——先1后2.
書(shū)中7.7.1對(duì)應(yīng)第1個(gè)問(wèn)題。.7.7.2對(duì)應(yīng)第2個(gè)問(wèn)題。
可以通過(guò)考察離群值和影響點(diǎn)來(lái)回答第一個(gè)問(wèn)題。
離群值和殘差
離群值(outlier)指與數(shù)據(jù)主體趨勢(shì)有實(shí)質(zhì)性不同的樣本點(diǎn)。在5.8.1節(jié)介紹過(guò)《Discovering Statistics Using R》筆記6-箱形圖和離群值
圖7.9展示了線性回歸中離群值對(duì)建立回歸方程的影響：
離群值會(huì)影響回歸模型的回歸系數(shù)從而產(chǎn)生系統(tǒng)誤差。圖中虛線展示的是原始線性模型，實(shí)線為加入離群值樣本點(diǎn)后擬合得到的線性模型。
殘差(residual)指模型預(yù)測(cè)值與因變量實(shí)測(cè)值的差值，反映模型的誤差?？梢酝ㄟ^(guò)查看殘差值很大的樣本點(diǎn)識(shí)別離群值。
標(biāo)準(zhǔn)化殘差
殘差帶有與因變量一致的單位，無(wú)法跨模型比較殘差，也無(wú)法判定“殘差很大”有統(tǒng)一判定標(biāo)準(zhǔn)。為克服這個(gè)困難可使用標(biāo)準(zhǔn)化殘差（standardized residuals）:殘差值除以殘差標(biāo)準(zhǔn)誤。
通過(guò)標(biāo)準(zhǔn)化，可對(duì)不同模型的殘差進(jìn)行比較并設(shè)立標(biāo)準(zhǔn)判斷殘差的可接受范圍：按照正態(tài)分布樣本規(guī)律，95%的Z值在-1.95，1.96之間；99%的Z值在-2.58，2.58之間；99.9%的Z值在-3.29，3.29之間。由此產(chǎn)生關(guān)于標(biāo)準(zhǔn)化殘差的一般規(guī)則：
- 標(biāo)準(zhǔn)化殘差絕對(duì)值大于3.29（可取近似值3）的樣本應(yīng)引起關(guān)注。
- 若有超過(guò)1%的樣本其標(biāo)準(zhǔn)化殘差絕對(duì)值大于2.58（可取近似值2.5），可認(rèn)為模型對(duì)樣本數(shù)據(jù)擬合較差。
- 若有超過(guò)5%的樣本其標(biāo)準(zhǔn)化殘差絕對(duì)值大于1.96（可取近似值2），可認(rèn)為模型對(duì)樣本數(shù)據(jù)擬合較差。
強(qiáng)影響點(diǎn)
刪除某個(gè)樣本后重新建模，看回歸系數(shù)是否發(fā)生明顯改變。這類分析可以考察回歸模型是否穩(wěn)定，是否有強(qiáng)影響點(diǎn)給模型帶來(lái)系統(tǒng)誤差。這個(gè)分析也能識(shí)別離群值。
評(píng)估樣本對(duì)模型影響的統(tǒng)計(jì)量有以下幾個(gè)：
- DFBeta 將某樣本剔除前后分別建立模型，兩模型的參數(shù)(回歸系數(shù)、截距)之差為DFBeta.
- DFFit 某樣本的調(diào)整預(yù)測(cè)值與原始預(yù)測(cè)值之差為DFFit。調(diào)整預(yù)測(cè)值（adjusted predicted value）是將該樣本點(diǎn)剔除后建立新模型，使用新模型計(jì)算出的該樣本點(diǎn)的預(yù)測(cè)值。若該樣本點(diǎn)非強(qiáng)影響點(diǎn)，DFFit值應(yīng)較小。
- 學(xué)生化殘差(studentized residual) 使用調(diào)整預(yù)測(cè)值計(jì)算出的殘差值除以其標(biāo)準(zhǔn)誤極為學(xué)生化殘差。學(xué)生化殘差可在不同模型間比較，并服從t分布。
- Cook距離(Cook's distance)Cook距離衡量某樣本點(diǎn)對(duì)模型的影響，Cook距離>1提示需要關(guān)注。
- 帽子值（hat values，也稱leverage）leverage的均值定義為（k+1）/n,其中k為自變量個(gè)數(shù)，n為樣本量。
leverage取值范圍為0-1。leverage越接近1表示樣本點(diǎn)對(duì)模型影響大。如果所有樣本都不是強(qiáng)影響點(diǎn)，則所有樣本的leverage值應(yīng)該都接近（k+1）/n。
Hoaglin and Welsch (1978)建議leverage值大于2(k+1)/n的樣本點(diǎn)視為強(qiáng)影響點(diǎn)；
Stevens (2002)建議leverage值大于3(k+1)/n的樣本點(diǎn)視為強(qiáng)影響點(diǎn). - Covariance ratio(CVR) 協(xié)方差比：評(píng)估樣本影響回歸模型參數(shù)方差的統(tǒng)計(jì)量。該值接近1表示對(duì)應(yīng)樣本對(duì)模型參數(shù)方差的影響很小。
Belsey et al.(1980)建議如果樣本的CVR＞1+[3(k+1)/n],剔除該樣本會(huì)減小模型參數(shù)估計(jì)的準(zhǔn)確性；如果樣本的CVR＜1-[3(k+1)/n]，則剔除該樣本有助于增加模型參數(shù)估計(jì)的準(zhǔn)確性。（k為自變量個(gè)數(shù)，n為樣本量）
R實(shí)操
使用的示例數(shù)據(jù)為：Album Sales 2.dat
自變量：
- adverts：廣告投入費(fèi)用
- airplay：唱片發(fā)布前1周內(nèi)，唱片中歌曲在廣播中播放的次數(shù)
- attract：樂(lè)隊(duì)的吸引程度。（打分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)
上面介紹的關(guān)于離群值和強(qiáng)影響點(diǎn)的診斷統(tǒng)計(jì)量是對(duì)樣本的，即每個(gè)樣本有對(duì)應(yīng)統(tǒng)計(jì)量。R中有很多函數(shù)計(jì)算這類統(tǒng)計(jì)量，一般的使用方法為：function(regressionModel)，即只需要將模型帶入函數(shù)即可。
- 計(jì)算離群值相關(guān)統(tǒng)計(jì)量的函數(shù)：
resid（）-計(jì)算殘差
rstandard()-計(jì)算標(biāo)準(zhǔn)化殘差
rstudent()-計(jì)算學(xué)生化殘差 - 計(jì)算強(qiáng)影響點(diǎn)相關(guān)統(tǒng)計(jì)量的函數(shù)：
cooks.distance()-計(jì)算Cook距離
dfbeta()-計(jì)算DFBeta
dffits()-計(jì)算DFFit
hatvalues()-計(jì)算leverage
covratio-計(jì)算協(xié)方差比
直接運(yùn)行這些函數(shù)，R會(huì)在console打印很長(zhǎng)的list，并不方便查看。建議將這些統(tǒng)計(jì)量存入數(shù)據(jù)集中。這里只查看前6行數(shù)據(jù)情況。
# 離群值與強(qiáng)影響點(diǎn)診斷
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
我們以查看標(biāo)準(zhǔn)化殘差為例，根據(jù)前面的介紹可知正常情況下約有95%的樣本其標(biāo)準(zhǔn)化殘差絕對(duì)值≤2（更準(zhǔn)確值為1.96）。示例數(shù)據(jù)有200個(gè)樣本，則預(yù)計(jì)會(huì)有約10個(gè)樣本的標(biāo)準(zhǔn)化殘差絕對(duì)值大于2.
使用dplyr包的filter()篩選出標(biāo)準(zhǔn)化殘差絕對(duì)值大于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
示例數(shù)據(jù)中有12個(gè)樣本（6%）的標(biāo)準(zhǔn)化殘差值絕對(duì)值大于2。擬合情況較好。
