R語言中的回歸診斷-car包

如何判斷我們的線性回歸模型是正確的？

1、回歸診斷的基本方法
opar<-par(no.readOnly=TRUE)

fit <- lm(weight ~ height, data = women)
par(mfrow = c(2, 2))
plot(fit)
par(opar)

為理解這些圖形，我們來回顧一下OLS回歸的統(tǒng)計假設。
(1)正態(tài)性（主要使用QQ圖）當預測變量值固定時，因變量成正態(tài)分布，則殘差值也應該是一個均值為0的正態(tài)分布。正態(tài)Q-Q圖（Normal Q-Q，右上）是在正態(tài)分布對應的值下，標準化殘差的概率圖。若滿足正態(tài)假設，那么圖上的點應該落在呈45度角的直線上；若不是如此，那么就違反了正態(tài)性的假設。
(2)獨立性 你無法從這些圖中分辨出因變量值是否相互獨立，只能從收集的數(shù)據(jù)中來驗證。上面的例子中，沒有任何先驗的理由去相信一位女性的體重會影響另外一位女性的體重。假若你發(fā)現(xiàn)數(shù)據(jù)是從一個家庭抽樣得來的，那么可能必須要調整模型獨立性的假設。
(3)線性（使用左上角的圖，該曲線盡量擬合所有點） 若因變量與自變量線性相關，那么殘差值與預測（擬合）值就沒有任何系統(tǒng)關聯(lián)。換句話說，除了白噪聲，模型應該包含數(shù)據(jù)中所有的系統(tǒng)方差。在“殘差圖與擬合圖”Residuals vs Fitted，左上）中可以清楚的看到一個曲線關系，這暗示著你可能需要對回歸模型加上一個二次項。
(4)同方差性（左下角，點隨機分布在曲線的周圍） 若滿足不變方差假設，那么在位置尺度圖（Scale-Location Graph，左下）中，水平線周圍的點應該隨機分布。該圖似乎滿足此假設。最后一幅“殘差與杠圖”（Residuals vs Leverage，右下）提供了你可能關注的單個觀測點的信息。從圖形可以鑒別出離群點、高杠桿值點和強影響點

通過看圖重新修改模型

newfit <- lm(weight ~ height + I(height^2), data = women[-c(13, 15),])
par(mfrow = c(2, 2))
plot(newfit)
par(opar)

2、使用改進的方法進行

主要使用的car包，進行回歸診斷

（1）自變量的正態(tài)分布

qqPlot()函數(shù)提供了更為精確的正態(tài)假設檢驗方法

library(car)

fit <- lm(Murder ~ Population + Illiteracy + Income + 

Frost, data = states)
qqPlot(fit, labels = FALSE, simulate = TRUE, main = "Q-Q Plot")

（2）誤差的獨立性

durbinWatsonTest(fit)

lag Autocorrelation D-W Statistic p-value
   1      -0.2006929      2.317691   0.248
 Alternative hypothesis: rho != 0

(3)線性相關性

crPlots(fit, one.page = TRUE, ask = FALSE)

（4）同方差性

1、car包提供了兩個有用的函數(shù)，可以判斷誤差方差是否恒定。ncvTest()函數(shù)生成一個計分檢驗，零假設為誤差方差不變，備擇假設為誤差方差隨著擬合值水平的變化而變化。

2、spreadLevelPlot()函數(shù)創(chuàng)建一個添加了最佳擬合曲線的散點圖，展示標準化殘差絕對值與擬合值的關系

library(car)
ncvTest(fit)

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 1.746514    Df = 1     p = 0.1863156    

滿足方差不變 p = 0.1863156   
spreadLevelPlot(fit)

3、線性模型假設的綜合驗證

library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)

Call:
lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
    data = states)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7960 -1.6495 -0.0811  1.4815  7.6210 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.235e+00  3.866e+00   0.319   0.7510    
Population  2.237e-04  9.052e-05   2.471   0.0173 *  
Illiteracy  4.143e+00  8.744e-01   4.738 2.19e-05 ***
Income      6.442e-05  6.837e-04   0.094   0.9253    
Frost       5.813e-04  1.005e-02   0.058   0.9541    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared:  0.567,     Adjusted R-squared:  0.5285 
F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08

ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = fit) 
                    Value p-value                Decision
Global Stat        2.7728  0.5965 Assumptions acceptable.
Skewness           1.5374  0.2150 Assumptions acceptable.
Kurtosis           0.6376  0.4246 Assumptions acceptable.
Link Function      0.1154  0.7341 Assumptions acceptable.
Heteroscedasticity 0.4824  0.4873 Assumptions acceptable.

4、多重共線性

如何檢測多重共線性

library(car)

vif(fit)

Population Illiteracy     Income      Frost 
  1.245282   2.165848   1.345822   2.082547 
sqrt(vif(fit)) > 2

Population Illiteracy     Income      Frost 
     FALSE      FALSE      FALSE      FALSE 

如何解決多重共線性？

逐步回歸法（此法最常用的，也最有效）

R語言回歸分析中的異常值點的介紹

（1）離群點

如何識別離群點？

1、Q-Q圖，落在置信區(qū)間帶［－2，2］外的點即可被認為是離群點。

2、一個粗糙的判斷準則：標準化殘差值大于2或者小于2的點可能是離群

3、library(car)
   outlierTest(fit)  顯示離群點

rstudent unadjusted p-value Bonferonni p
Nevada 3.542929         0.00095088     0.047544

（2）高杠桿值點

它們是由許多異常的預測變量值組合起來的，與響應變量值沒有關系

高杠桿值的觀測點可通過帽子統(tǒng)計量（hat statistic）判斷

hat.plot <- function(fit){
    p <- length(coefficients(fit))
    n <- length(fitted(fit))
    plot(hatvalues(fit), main = "Index Plot of Hat Values")
    abline(h = c(2, 3) * p/n, col = "red", lty = 2)
    identify(1:n, hatvalues(fit), names(hatvalues(fit)))
}
hat.plot(fit)

（3）強影響點

強影響點，即對模型參數(shù)估計值影響有些比例失衡的點。例如，若移除模型的一個觀測點時模型會發(fā)生巨大的改變，那么你就需要檢測一下數(shù)據(jù)中是否存在強影響點了

cutoff <- 4/(nrow(states) - length(fit$coefficients) - 2)
plot(fit, which = 4, cook.levels = cutoff)
abline(h = cutoff, lty = 2, col = "red")

4、如何對線性模型進行改進？

1、刪除觀測點；

刪除離群點通?？梢蕴岣邤?shù)據(jù)集對于正態(tài)假設的擬合度，而強影響點會干擾結果，通常也會被刪除。刪除最大的離群點或者強影響點后，模型需要重新擬合
2、變量變換：

Box-Cox正態(tài)變換

library(car)
summary(powerTransform(states$Murder))

library(car)
boxTidwell(Murder ~ Population + Illiteracy, data = states)
3、添加或刪除變量；
4、使用其他回歸方法。