【原】各類統(tǒng)計(jì)方法R語(yǔ)言實(shí)現(xiàn)（七）

醫(yī)科研 2021-01-25

展開全文

今天是各類統(tǒng)計(jì)方法R語(yǔ)言實(shí)現(xiàn)的第七期，我們主要介紹多重共線性、異常觀察值的分析和回歸模型改進(jìn)措施。

多重共線性

多重共線性是指線性回歸模型中的解釋變量之間由于存在強(qiáng)相關(guān)關(guān)系而使模型估計(jì)失真或難以估計(jì)準(zhǔn)確，它會(huì)導(dǎo)致模型參數(shù)的置信區(qū)間過(guò)大，使參數(shù)解釋較困難。

多重共線性可用VIF（Variance Inflation Factor，方差膨脹因子）進(jìn)行檢測(cè)，該指標(biāo)的經(jīng)驗(yàn)判斷方法：VIF在5到10之間：中度共線性。VIF大于10：重度共線性。

多重共線性解決方法

手動(dòng)移除出共線性的自變量
逐步回歸法
增加樣本容量
嶺回歸或lasso回歸
利用因子分析合并變量

#模型擬合
fit<-lm(mpg~hp+wt,data=mtcars)

##展示模型
summary(fit)
##
## Call:
## lm(formula = mpg ~ hp + wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.941 -1.600 -0.182 1.050 5.854
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.22727 1.59879 23.285 < 2e-16 ***
## hp -0.03177 0.00903 -3.519 0.00145 **
## wt -3.87783 0.63273 -6.129 1.12e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.593 on 29 degrees of freedom
## Multiple R-squared: 0.8268, Adjusted R-squared: 0.8148
## F-statistic: 69.21 on 2 and 29 DF, p-value: 9.109e-12
#計(jì)算vif
library(car)
## Loading required package: carData
vif(fit)
## hp wt
## 1.766625 1.766625

不存在多重共線性

異常觀測(cè)值

上次推文已經(jīng)介紹了異常觀測(cè)值主要有三類：離群值點(diǎn)、高杠桿值點(diǎn)、強(qiáng)影響點(diǎn)，具體如下：

a.離群點(diǎn)：擬合回歸模型對(duì)其預(yù)測(cè)效果不佳（即殘差的絕對(duì)值較大）。

b.有高杠桿值的變量表明它是一個(gè)異常的自變量組合。

c.強(qiáng)影響點(diǎn)表明他對(duì)模型參數(shù)的估計(jì)產(chǎn)生的影響過(guò)大。

離群點(diǎn)

之前已經(jīng)介紹在標(biāo)準(zhǔn)化殘差的QQ圖中，偏離其他值的異常點(diǎn)可能是離群點(diǎn)，一般認(rèn)為標(biāo)準(zhǔn)化殘差絕對(duì)值大于2的點(diǎn)為離群點(diǎn)。

接下來(lái)介紹另一種判斷離群值的方法，即使用car包中的outlierTest()函數(shù)。

fit2<-lm(weight ~ height + I(height^2),data = women)
summary(fit2)
##
## Call:
## lm(formula = weight ~ height + I(height^2), data = women)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50941 -0.29611 -0.00941 0.28615 0.59706
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 261.87818 25.19677 10.393 2.36e-07 ***
## height -7.34832 0.77769 -9.449 6.58e-07 ***
## I(height^2) 0.08306 0.00598 13.891 9.32e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3841 on 12 degrees of freedom
## Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994
## F-statistic: 1.139e+04 on 2 and 12 DF, p-value: < 2.2e-16
library(car)
outlierTest(fit2)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 15 2.575781 0.025783 0.38675

可以看出15號(hào)點(diǎn)p=0.38675，不顯著，表明沒有離群點(diǎn)。注意：outlierTest（）函數(shù)是根據(jù)單個(gè)最大殘差值(絕對(duì)值）的顯著性來(lái)判斷是否有離群點(diǎn)，若不顯著，則說(shuō)明數(shù)據(jù)集中沒有離群點(diǎn)，若顯著，則必須刪除該離群點(diǎn)，然后再檢驗(yàn)是否還有其他離群點(diǎn)存在。

高杠桿值點(diǎn)

有高杠桿值的變量表明它是一個(gè)異常的自變量組合，即由許多異常的自變量組合起來(lái)的異常點(diǎn)，與因變量值沒有關(guān)系。

高杠桿值的觀測(cè)點(diǎn)可通過(guò)帽子統(tǒng)計(jì)量（hat statistic）判斷。對(duì)于一個(gè)給定的數(shù)據(jù)集，帽子均值為p/n，其中p是模型估計(jì)的參數(shù)數(shù)目（包含截距項(xiàng)），n是樣本量。

一般來(lái)說(shuō)，若觀測(cè)點(diǎn)的帽子值大于帽子均值的2或3倍，則可認(rèn)定為高杠桿值點(diǎn)。

hat.plot<-function(fit){
p<-length(coefficients(fit))
n<-length(fitted(fit))
plot(hatvalues(fit),main="Index Plot of Hat Values",ylim=c(0,3*p/n)+0.2)

abline(h=c(2,3)*p/n,col="red",lty=2)
identify(1:n,hatvalues(fit),names(hatvalues(fit)))
}
hat.plot(fit2)

## integer(0)

簡(jiǎn)化代碼

##簡(jiǎn)化代碼
hat<-hatvalues(fit2)
hat_mean<-mean(hat)
plot(hat,ylim=c(0,3*hat_mean))
abline(h=c(2,3)*mean(hatvalues(fit2)) , col="red",lty=2)

水平的兩根紅線表示帽子均值的2和3倍,可以看出1號(hào)點(diǎn)和15號(hào)點(diǎn)超過(guò)了2倍但沒到3倍。

強(qiáng)影響點(diǎn)

表明某點(diǎn)對(duì)模型參數(shù)的估計(jì)產(chǎn)生的影響過(guò)大，即移除該點(diǎn)，模型會(huì)發(fā)生巨大的變化。

檢測(cè)方法：

Cook距離，或稱為D統(tǒng)計(jì)量：Cook’s D值大于4/(n-k-1)，則表明它是強(qiáng)影響點(diǎn)，其中n為樣本量大小，k是預(yù)測(cè)變量數(shù)目（有助于鑒別強(qiáng)影響點(diǎn)，但并不提供關(guān)于這些點(diǎn)如何影響模型的信息）

變量添加圖（added variable plot）：對(duì)于每個(gè)預(yù)測(cè)變量Xk，繪制Xk在其他k-1個(gè)預(yù)測(cè)變量上回歸的殘差值相對(duì)于響應(yīng)變量在其他k-1個(gè)預(yù)測(cè)變量上回歸的殘差值的關(guān)系圖。

#Cook距離
cutoff<-4/(nrow(women)-length(fit2$coefficients)-2)
plot(fit2,which=4,cook.levels=cutoff)
abline(h=cutoff,lty=2,col="red")

紅線表示4/(n-k-1)，可以發(fā)現(xiàn)15號(hào)cook距離最大，與上次推文結(jié)果一致。

#變量添加圖
library(car)
avPlots(fit2,ask=FALSE,onepage=TRUE,id.method="identify")

對(duì)于此圖，可以想象去掉某一個(gè)點(diǎn)之后，直線擬合是否會(huì)有大范圍變動(dòng)，此處15號(hào)點(diǎn)的影響在所有點(diǎn)中算是比較大的了。

結(jié)果整合

car包中的influencePlot（）函數(shù)

hat<-hatvalues(fit2)
hat_mean<-mean(hat)

library(car)
influencePlot(fit2,id.method="identify",main="Influence Plot",
sub="Circle size if proportional to Cook's distance",
xlim=c(0,3*hat_mean))

## StudRes Hat CookD
## 1 -0.3527249 0.4647059 0.03883656
## 2 -1.5156988 0.2680672 0.25310078
## 13 -1.5312900 0.1656755 0.13956756
## 15 2.5757809 0.4647059 1.30646190

本質(zhì)上是將三個(gè)值繪制在一張圖里。

縱坐標(biāo)超過(guò)2或小于-2的點(diǎn)可被認(rèn)為是離群點(diǎn)，水平軸超過(guò)2倍或3倍帽子值均值的點(diǎn)有高杠桿值。圓圈大小與影響成比例，圓圈很大的點(diǎn)可能是對(duì)模型估計(jì)造成的不成比例影響的強(qiáng)影響點(diǎn)。

回歸模型改進(jìn)措施

主要有四種方法：

（1）刪除異常值.

（2）變量變換。

（3）添加或刪除變量。

（4）使用其他回歸方法。

刪除異常值

通常刪除離群點(diǎn)和強(qiáng)影響點(diǎn)，直到擬合較滿意。

當(dāng)然刪除要謹(jǐn)慎，可以探究產(chǎn)生異常值的原因。

變量變換

可以嘗試各類變換方法，使變量滿足正態(tài)性、線性、同方差性假設(shè)，可以嘗試之前各類統(tǒng)計(jì)方法R語(yǔ)言實(shí)現(xiàn)（四）介紹的方法，但是變量變換之后需要有具體意義。

添加或刪除變量

可以刪除多重共線性的變量（根據(jù)VIF），也可以嶺回歸或lasso回歸。

使用其他回歸方法

存在離群點(diǎn)或強(qiáng)影響點(diǎn)，可使用穩(wěn)健回歸代替最小二乘回歸。

違背了正態(tài)性假設(shè)，可以用非參數(shù)回歸模型。

存在顯著非線性，可以使用非線性模型。

違背了誤差獨(dú)立性假設(shè)，可以使用專門研究誤差結(jié)構(gòu)的模型，如時(shí)間序列模型或多層次回歸模型。

最后，還能依據(jù)數(shù)據(jù)的分布形式選擇不同的廣義線性模型。

好了，今天的R語(yǔ)言實(shí)現(xiàn)統(tǒng)計(jì)方法系列推文暫時(shí)告一段落，我們下次再見吧！小伙伴們?nèi)绻惺裁唇y(tǒng)計(jì)上的問(wèn)題，或者如果想要學(xué)習(xí)什么方面的生物信息內(nèi)容，可以在微信群或者知識(shí)星球提問(wèn)，沒準(zhǔn)哪天的推文就是專門解答你的問(wèn)題哦！