柚子快報(bào)激活碼778899分享：R語言統(tǒng)計(jì)分析——回歸診斷2

http://yzkb.51969.com/

參考資料：R語言實(shí)戰(zhàn)【第2版】

????????R語言的car包提供的大量函數(shù)，大大增強(qiáng)了擬合和評(píng)價(jià)回歸模型的能力。如下：

函數(shù)目的qqPlot()分位數(shù)比較圖durbinWatsonTest()對(duì)誤差自相關(guān)性做Durbin-Watson檢驗(yàn)crPlots()成分與殘差圖ncvTest()對(duì)非恒定的誤差方差做得分檢驗(yàn)spreadLevelPlot()分散水平檢驗(yàn)outlierTest() Bonferroni離群點(diǎn)檢驗(yàn) avPlots()添加的變量圖形influencePlot()回歸影響圖scatterplot()增強(qiáng)的散點(diǎn)圖scatterplotMatrix()增強(qiáng)的散點(diǎn)圖矩陣vif()方差膨脹因子

? ? ? ? 另外，gvlma包提供了對(duì)所有線性模型假設(shè)進(jìn)行檢驗(yàn)的方法。

1、正態(tài)性

? ? ? ? 與基礎(chǔ)包中的plot()函數(shù)相比，qqPlot()函數(shù)提供了更為精確的正態(tài)假設(shè)檢驗(yàn)方法，它畫出了在n-p-1個(gè)自由度的圖分布下的學(xué)生化殘差（studentized residual，也稱學(xué)生化刪除殘差或折疊化殘差）圖形，其中n是樣本大小，p是回歸參數(shù)的數(shù)目（包括截距項(xiàng)）。如下：

# 加載car包

library(car)

# 獲取數(shù)據(jù)

states<-as.data.frame(state.x77[,c("Murder","Population",

"Illiteracy","Income","Frost")])

# 擬合多元線性模型

fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)

# 繪制Q-Q圖

qqPlot(fit,

simulate=TRUE,

main="Q-Q Plot")

? ? ? ? ?當(dāng)simulate=TRUE時(shí)，將會(huì)使用參數(shù)自助法生成95%的置信區(qū)間。

? ? ? ? 從上圖中可以看出，除了Nevada，所有的點(diǎn)都離直線很近，并都落在置信區(qū)間內(nèi)，這表明正態(tài)性假設(shè)符合得很好。作為一個(gè)異常值，我們也應(yīng)該關(guān)注Nevada，它有一個(gè)很大的正殘差值（真實(shí)值-預(yù)測(cè)值），表明模型低估了該地的謀殺率。

# 查看Nevada數(shù)據(jù)

states["Nevada",]

# 查看Nevada的擬合數(shù)據(jù)（預(yù)測(cè)值）

fitted(fit)["Nevada"]

# 查看Nevada的殘差

residuals(fit)["Nevada"]

# 查看Nevada的學(xué)生化殘差

rstudent(fit)["Nevada"]

注：學(xué)生化殘差是殘差除以它的標(biāo)準(zhǔn)差后得到的數(shù)值，用以直觀地判斷誤差項(xiàng)服從正態(tài)分布這一假定是否成立，若假定成立，學(xué)生化殘差的分布也應(yīng)服從正態(tài)分布?！緦W(xué)生化殘差_百度百科】

? ? ? ? 可視化誤差還有其他方法，比如使用residplot()函數(shù)生成學(xué)生化殘差柱狀圖（即直方圖），并添加正態(tài)曲線、核密度曲線和軸須圖。

resiplot<-function(fit,nbreaks=10){

z<-rstudent(fit)

hist(z,breaks=nbreaks,freq=FALSE,

xlab="Studentized Residual",

main="Distribution of Error")

rug(jitter(z),col="brown")

curve(dnorm(x,mean=mean(z),sd=sd(z)),

add=TRUE,col="blue",lwd=2)

lines(density(z)$x,density(z)$y,

col="red",lwd=2,lty=2)

legend("topright",

legend=c("Normal Curve","Kernel Density Curve"),

lty=1:2,col=c("blue","red"),cex=0.7)

}

resiplot(fit)

? ? ? ? 如上圖所示，除了一個(gè)明顯的離群點(diǎn)，誤差基本上都很好的服從了正態(tài)分布。

2、誤差的獨(dú)立性

? ? ? ? 判斷因變量值（或殘差）是否相互獨(dú)立，最好的方法是依據(jù)收集數(shù)據(jù)方式的先驗(yàn)知識(shí)。例如，時(shí)間序列數(shù)據(jù)通常呈現(xiàn)自相關(guān)——相隔時(shí)間越近的觀測(cè)相關(guān)性大于相隔越遠(yuǎn)的觀測(cè)。car包提供了一個(gè)可做Durbin-Watson檢驗(yàn)的函數(shù)，能夠檢測(cè)誤差的序列相關(guān)性。繼續(xù)使用上面的案例進(jìn)行Durbin-Watson檢驗(yàn)：

durbinWatsonTest(fit)

? ? ? ? 結(jié)果中，p值不顯著（p=0.242）說明無自相關(guān)性，誤差之間相互獨(dú)立。滯后項(xiàng)（lag=1）表明數(shù)據(jù)集中每個(gè)數(shù)據(jù)都是與其中后一個(gè)數(shù)據(jù)進(jìn)行比較的。該檢驗(yàn)適用于時(shí)間獨(dú)立的數(shù)據(jù)，對(duì)于非聚集型的數(shù)據(jù)并不使用。注意，durbinWatsonTest()函數(shù)使用的自助法來導(dǎo)出p值。如果添加了選項(xiàng)simulate=TRUE，則每次運(yùn)行測(cè)試時(shí)獲得的結(jié)果都將略有不同。

3、線性

? ? ? ? 通過成分殘差圖（component plus residual plot）也稱偏殘差圖（partial residual plot），我們可以查看因變量與自變量之間是否呈非線性關(guān)系，也可以查看是否有不同于已設(shè)定線性模型的系統(tǒng)偏差，圖形可用car包中的crPlots()函數(shù)繪制。繼續(xù)使用上面的案例數(shù)據(jù)，如下：

# 導(dǎo)入car包

library(car)

# 繪制成分殘差圖

crPlots(fit)

? ? ? ? 如果成分殘差圖中的圖形存在非線性，則說明我們對(duì)自變量的函數(shù)形式建模不夠充分。那么就需要添加一些曲線成分，如多項(xiàng)式，或?qū)σ粋€(gè)或多個(gè)變量進(jìn)行變換（如log(x)代替x），或用其他回歸變體形式而不是線性回歸。

? ? ? ? 從上圖中可以看出，4個(gè)自變量的成分殘差圖可以認(rèn)為是線性的，尤其是前兩個(gè)自變量。

4、同方差性

? ? ? ? car包提供了兩個(gè)有用的函數(shù)，可以判斷誤差方差是否恒定。

????????ncvTest()函數(shù)生成一個(gè)計(jì)分檢驗(yàn)，零假設(shè)為誤差方差不變，備擇假設(shè)為誤差方差隨著擬合值水平的變化而變化。若檢驗(yàn)顯著，則說明存在異方差性（誤差方差不恒定）。

? ? ? ? spreadLevelPlot()函數(shù)創(chuàng)建一個(gè)添加了最佳擬合曲線的散點(diǎn)圖，展示標(biāo)準(zhǔn)化殘差決對(duì)值與擬合值得關(guān)系。

# 導(dǎo)入car包

library(car)

ncvTest(fit)

spreadLevelPlot(fit)

? ? ? ? 根據(jù)ncvTest()的計(jì)分檢驗(yàn)結(jié)果顯示不顯著（p=0.186），說明滿足方差不變假設(shè)。?

? ? ? ? spreadLevelPlot()函數(shù)作圖顯示，數(shù)據(jù)點(diǎn)在水平的最佳擬合曲線（虛線）周圍呈水平隨機(jī)分布。如果違反了該假設(shè)，我們將看到一個(gè)非水平的曲線。而運(yùn)行代碼結(jié)果顯示“Suggested power transformation: 1.209626”，建議進(jìn)行1.2次冪的冪次轉(zhuǎn)換。因接近1，所以不需要進(jìn)行變換。

柚子快報(bào)激活碼778899分享：R語言統(tǒng)計(jì)分析——回歸診斷2

http://yzkb.51969.com/

精彩文章