筆記說明
讀《Discovering Statistics Using R》第七章 Regression中的7.7-節(jié)做的筆記竞惋。使用R進(jìn)行實操部分來自節(jié)姜胖。
模型診斷概述
基于樣本數(shù)據(jù)建立模型后有兩個問題要考慮:
1.模型對樣本數(shù)據(jù)擬合是否夠好?/模型是否被少數(shù)樣本點所影響?
2.模型是否能外推到其他樣本?
這兩個問題是有層次的——先1后2.
書中7.7.1對應(yīng)第1個問題。.7.7.2對應(yīng)第2個問題氢烘。
可以通過考察離群值和影響點來回答第一個問題。
離群值和殘差
離群值(outlier)指與數(shù)據(jù)主體趨勢有實質(zhì)性不同的樣本點家厌。在5.8.1節(jié)介紹過《Discovering Statistics Using R》筆記6-箱形圖和離群值
圖7.9展示了線性回歸中離群值對建立回歸方程的影響:
離群值會影響回歸模型的回歸系數(shù)從而產(chǎn)生系統(tǒng)誤差播玖。圖中虛線展示的是原始線性模型,實線為加入離群值樣本點后擬合得到的線性模型饭于。
殘差(residual)指模型預(yù)測值與因變量實測值的差值蜀踏,反映模型的誤差维蒙。可以通過查看殘差值很大的樣本點識別離群值果覆。
標(biāo)準(zhǔn)化殘差
殘差帶有與因變量一致的單位颅痊,無法跨模型比較殘差,也無法判定“殘差很大”有統(tǒng)一判定標(biāo)準(zhǔn)局待。為克服這個困難可使用標(biāo)準(zhǔn)化殘差(standardized residuals):殘差值除以殘差標(biāo)準(zhǔn)誤斑响。
通過標(biāo)準(zhǔn)化,可對不同模型的殘差進(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)化殘差絕對值大于3.29(可取近似值3)的樣本應(yīng)引起關(guān)注饲漾。
- 若有超過1%的樣本其標(biāo)準(zhǔn)化殘差絕對值大于2.58(可取近似值2.5),可認(rèn)為模型對樣本數(shù)據(jù)擬合較差候醒。
- 若有超過5%的樣本其標(biāo)準(zhǔn)化殘差絕對值大于1.96(可取近似值2)能颁,可認(rèn)為模型對樣本數(shù)據(jù)擬合較差杂瘸。
強(qiáng)影響點
刪除某個樣本后重新建模倒淫,看回歸系數(shù)是否發(fā)生明顯改變。這類分析可以考察回歸模型是否穩(wěn)定败玉,是否有強(qiáng)影響點給模型帶來系統(tǒng)誤差敌土。這個分析也能識別離群值。
評估樣本對模型影響的統(tǒng)計量有以下幾個:
- DFBeta 將某樣本剔除前后分別建立模型运翼,兩模型的參數(shù)(回歸系數(shù)返干、截距)之差為DFBeta.
- DFFit 某樣本的調(diào)整預(yù)測值與原始預(yù)測值之差為DFFit。調(diào)整預(yù)測值(adjusted predicted value)是將該樣本點剔除后建立新模型血淌,使用新模型計算出的該樣本點的預(yù)測值矩欠。若該樣本點非強(qiáng)影響點,DFFit值應(yīng)較小悠夯。
- 學(xué)生化殘差(studentized residual) 使用調(diào)整預(yù)測值計算出的殘差值除以其標(biāo)準(zhǔn)誤極為學(xué)生化殘差癌淮。學(xué)生化殘差可在不同模型間比較,并服從t分布沦补。
- Cook距離(Cook's distance)Cook距離衡量某樣本點對模型的影響乳蓄,Cook距離>1提示需要關(guān)注。
- 帽子值(hat values夕膀,也稱leverage)leverage的均值定義為(k+1)/n,其中k為自變量個數(shù)虚倒,n為樣本量美侦。
leverage取值范圍為0-1。leverage越接近1表示樣本點對模型影響大魂奥。如果所有樣本都不是強(qiáng)影響點菠剩,則所有樣本的leverage值應(yīng)該都接近(k+1)/n。
Hoaglin and Welsch (1978)建議leverage值大于2(k+1)/n的樣本點視為強(qiáng)影響點耻煤;
Stevens (2002)建議leverage值大于3(k+1)/n的樣本點視為強(qiáng)影響點. - Covariance ratio(CVR) 協(xié)方差比:評估樣本影響回歸模型參數(shù)方差的統(tǒng)計量赠叼。該值接近1表示對應(yīng)樣本對模型參數(shù)方差的影響很小。
Belsey et al.(1980)建議如果樣本的CVR>1+[3(k+1)/n],剔除該樣本會減小模型參數(shù)估計的準(zhǔn)確性违霞;如果樣本的CVR<1-[3(k+1)/n]嘴办,則剔除該樣本有助于增加模型參數(shù)估計的準(zhǔn)確性。(k為自變量個數(shù)买鸽,n為樣本量)
R實操
使用的示例數(shù)據(jù)為:Album Sales 2.dat
自變量:
- adverts:廣告投入費用
- airplay:唱片發(fā)布前1周內(nèi)涧郊,唱片中歌曲在廣播中播放的次數(shù)
- 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)
上面介紹的關(guān)于離群值和強(qiáng)影響點的診斷統(tǒng)計量是對樣本的妆艘,即每個樣本有對應(yīng)統(tǒng)計量。R中有很多函數(shù)計算這類統(tǒng)計量看幼,一般的使用方法為:function(regressionModel)
批旺,即只需要將模型帶入函數(shù)即可。
- 計算離群值相關(guān)統(tǒng)計量的函數(shù):
resid()
-計算殘差
rstandard()
-計算標(biāo)準(zhǔn)化殘差
rstudent()
-計算學(xué)生化殘差 - 計算強(qiáng)影響點相關(guān)統(tǒng)計量的函數(shù):
cooks.distance()
-計算Cook距離
dfbeta()
-計算DFBeta
dffits()
-計算DFFit
hatvalues()
-計算leverage
covratio
-計算協(xié)方差比
直接運行這些函數(shù)诵姜,R會在console打印很長的list汽煮,并不方便查看。建議將這些統(tǒng)計量存入數(shù)據(jù)集中棚唆。這里只查看前6行數(shù)據(jù)情況暇赤。
# 離群值與強(qiáng)影響點診斷
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)化殘差絕對值≤2(更準(zhǔn)確值為1.96)宵凌。示例數(shù)據(jù)有200個樣本鞋囊,則預(yù)計會有約10個樣本的標(biāo)準(zhǔn)化殘差絕對值大于2.
使用dplyr
包的filter()
篩選出標(biāo)準(zhǔn)化殘差絕對值大于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個樣本(6%)的標(biāo)準(zhǔn)化殘差值絕對值大于2。擬合情況較好瞎惫。