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

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

展開(kāi)全文

【科研繪圖點(diǎn)我】【付費(fèi)精品合集】【SEER點(diǎn)我】

今天是各類統(tǒng)計(jì)方法R語(yǔ)言實(shí)現(xiàn)的第八期,，我們主要介紹選擇“最佳”回歸模型與深層次分析,。

選擇“最佳”回歸模型

當(dāng)我們構(gòu)建回歸方程時(shí),，我們一方面需要考慮是否去除不顯著的變量,、是否需要添加交互項(xiàng)和/或多項(xiàng)式項(xiàng),，另一方面還需要權(quán)衡方程的簡(jiǎn)潔性和精度,。

因此,，我們希望最終構(gòu)建的“最佳”回歸方程,，應(yīng)該至少兼有簡(jiǎn)潔和有效兩個(gè)特點(diǎn),。

模型比較方法

方差分析

anova（）函數(shù)可比較兩個(gè)嵌套模型的擬合優(yōu)度。嵌套模型即指一個(gè)模型的一些項(xiàng)完全包含在另一個(gè)模型中,。

states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])

fit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
fit2<-lm(Murder~Population+Illiteracy,data=states)
anova(fit2,fit1)
## Analysis of Variance Table
## 
## Model 1: Murder ~ Population + Illiteracy
## Model 2: Murder ~ Population + Illiteracy + Income + Frost
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     47 289.25                           
## 2     45 289.17  2  0.078505 0.0061 0.9939

p=0.9939,，兩個(gè)模型無(wú)顯著差異，說(shuō)明剔除模型1中的Income,、Frost對(duì)于回歸結(jié)果影響不大,。

AIC法

AIC即Akaike Information Criterion，赤池信息準(zhǔn)則,，考慮了模型的統(tǒng)計(jì)擬合度及用來(lái)擬合的參數(shù)數(shù)目 ,。（不需要嵌套模型）

AIC值越小的模型可優(yōu)先選擇，說(shuō)明模型用較少的參數(shù)獲得了足夠的擬合度

fit1<-lm(Murder~Population+Illiteracy+Income+Frost,
         data=states)
fit2<-lm(Murder~Population+Illiteracy,data=states)
AIC(fit1,fit2)
##      df      AIC
## fit1  6 241.6429
## fit2  4 237.6565

模型1的AIC為241.6429,，模型2的AIC為237.6565,，表面模型2效果更好，與方差分析法結(jié)果一致,。

變量選擇

初步回歸法,、全子集回歸法。

1,、逐步回歸法

模型每次增加或刪除一個(gè)變量,，直到達(dá)到某個(gè)判停準(zhǔn)則（如AIC），可分為向前,、向后和雙向逐步回歸,。

向前：每次添加一個(gè)預(yù)測(cè)變量到模型中，直到添加變量不會(huì)使模型有所改進(jìn)為止,。

向后：從模型包含所有預(yù)測(cè)變量開(kāi)始,，一次刪除一個(gè)變量直到會(huì)降低模型質(zhì)量為止。

雙向：綜合向前向后,。

##向后回歸
library(MASS)
fit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
stepAIC(fit1,direction="backward")
## Start:  AIC=97.75
## Murder ~ Population + Illiteracy + Income + Frost
## 
##              Df Sum of Sq    RSS     AIC
## - Frost       1     0.021 289.19  95.753
## - Income      1     0.057 289.22  95.759
## <none>                    289.17  97.749
## - Population  1    39.238 328.41 102.111
## - Illiteracy  1   144.264 433.43 115.986
## 
## Step:  AIC=95.75
## Murder ~ Population + Illiteracy + Income
## 
##              Df Sum of Sq    RSS     AIC
## - Income      1     0.057 289.25  93.763
## <none>                    289.19  95.753
## - Population  1    43.658 332.85 100.783
## - Illiteracy  1   236.196 525.38 123.605
## 
## Step:  AIC=93.76
## Murder ~ Population + Illiteracy
## 
##              Df Sum of Sq    RSS     AIC
## <none>                    289.25  93.763
## - Population  1    48.517 337.76  99.516
## - Illiteracy  1   299.646 588.89 127.311
## 
## Call:
## lm(formula = Murder ~ Population + Illiteracy, data = states)
## 
## Coefficients:
## (Intercept)   Population   Illiteracy  
##   1.6515497    0.0002242    4.0807366

一開(kāi)始包含四個(gè)變量,，Population + Illiteracy + Income + Frost，最終得到兩個(gè)變量Population + Illiteracy,；AIC從97.75到95.75再到93.76,，逐步降低,。

2、全子集回歸法

全子集回歸,，即所有可能的模型都會(huì)被檢驗(yàn),，可選擇展示所有可能的結(jié)果，也可展示n個(gè)不同變量（一個(gè),、兩個(gè)或多個(gè)預(yù)測(cè)變量）的最佳模型

可通過(guò)模型R平方,、調(diào)整R平方或Mallows Cp統(tǒng)計(jì)量等準(zhǔn)則來(lái)選擇“最佳”模型

調(diào)整R平方比R平方增加了模型參數(shù)數(shù)目信息。

##調(diào)整R平方
library(leaps)
leaps<-regsubsets(Murder~Population+Illiteracy+Income+Frost,data=states,nbest=4)
plot(leaps,scale="adjr2")

img

該圖第一行包含兩個(gè)變量Population + Illiteracy,，調(diào)整后的R平方最高,，是最佳模型，和向后回歸結(jié)果一致,。

大部分情況全子集回歸法較逐步回歸法更優(yōu),。

深層次分析

評(píng)價(jià)模型泛化能力和變量相對(duì)重要性的方法，包括交叉驗(yàn)證和相對(duì)重要性,。

1,、交叉驗(yàn)證

交叉驗(yàn)證是模型內(nèi)驗(yàn)證的一種方法，將一定比例的數(shù)據(jù)挑選出來(lái)作為訓(xùn)練樣本,，另外的樣本作為保留樣本,，先在訓(xùn)練樣本上獲取回歸方程，然后在保留樣本上做預(yù)測(cè),。由于保留樣本不涉及模型及參數(shù)的選擇,，該樣本可獲得比新數(shù)據(jù)更為精確的估計(jì)。

最常用的交叉驗(yàn)證方法是k重交叉驗(yàn)證,，在k重交叉驗(yàn)證中,，樣本被平均分為k個(gè)子樣本，輪流將k-1個(gè)子樣本組合作為訓(xùn)練集,，另外1個(gè)子樣本作為保留集,，這樣會(huì)獲得k個(gè)預(yù)測(cè)方程，記錄k個(gè)保留樣本的預(yù)測(cè)表現(xiàn)結(jié)果,，然后求其平均值,。

#10重交叉驗(yàn)證
set.seed(123)
library(bootstrap)
shrinkage<-function(fit,k=10){
  require(bootstrap)
  theta.fit<-function(x,y){lsfit(x,y)}
  theta.predict<-function(fit,x){cbind(1,x)%*%fit$coef}
  x<-fit$model[,2:ncol(fit$model)]
  y<-fit$model[,1]
  results<-crossval(x,y,theta.fit,theta.predict,ngroup=k)
  r2<-cor(y,fit$fitted.values)^2
  r2cv<-cor(y,results$cv.fit)^2
  cat("Original R-square=",r2,"\n")
  cat(k,"Fold Cross-Validated R-square=",r2cv,"\n")
  cat("Change=",r2-r2cv,"\n")
}
fit<-lm(Murder~Population+Income+Illiteracy+Frost,data=states)
shrinkage(fit)
## Original R-square= 0.5669502 
## 10 Fold Cross-Validated R-square= 0.492987 
## Change= 0.0739632
fit2<-lm(Murder~Population+Illiteracy,data=states)
shrinkage(fit2)
## Original R-square= 0.5668327 
## 10 Fold Cross-Validated R-square= 0.5296188 
## Change= 0.03721387

R平方減少越少,，結(jié)果越精確,，可以看到第二個(gè)模型更穩(wěn)健，泛化能力更強(qiáng),。

2,、相對(duì)重要性

評(píng)價(jià)模型中哪個(gè)指標(biāo)更重要，最簡(jiǎn)單的方法是比較標(biāo)準(zhǔn)化的回歸系數(shù),。

#z-score標(biāo)準(zhǔn)化
zstates<-as.data.frame(scale(states))
zfit<-lm(Murder~Population+Income+Illiteracy+Frost,data=zstates)
coef(zfit)
##   (Intercept)    Population        Income    Illiteracy         Frost 
## -2.054026e-16  2.705095e-01  1.072372e-02  6.840496e-01  8.185407e-03

Illiteracy最重要,， Frost最不重要,。

最小二乘回歸小結(jié)

前幾次推文主要介紹了最小二乘回歸，尤其是簡(jiǎn)單線性回歸相關(guān)的分析,，可以看出構(gòu)建一個(gè)回歸模型有許多注意事項(xiàng)和技巧,，當(dāng)然，之前的大多數(shù)方法有一個(gè)前提,，即變量服從正態(tài)分布,，那么假如我們的變量不服從正態(tài)分布怎么辦呢？因此,，之后一段時(shí)間,，我們將介紹重抽樣法、自主法與廣義線性回歸,。

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