今天是各類統(tǒng)計(jì)方法R語言實(shí)現(xiàn)的第八期，我們主要介紹選擇“最佳”回歸模型與深層次分析。
選擇“最佳”回歸模型
當(dāng)我們構(gòu)建回歸方程時(shí)，我們一方面需要考慮是否去除不顯著的變量、是否需要添加交互項(xiàng)和/或多項(xiàng)式項(xiàng)，另一方面還需要權(quán)衡方程的簡潔性和精度。
因此，我們希望最終構(gòu)建的“最佳”回歸方程，應(yīng)該至少兼有簡潔和有效兩個(gè)特點(diǎn)。
模型比較方法
方差分析
anova（）函數(shù)可比較兩個(gè)嵌套模型的擬合優(yōu)度。嵌套模型即指一個(gè)模型的一個(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è)模型無顯著差異，說明剔除模型1中的Income、Frost對(duì)于回歸結(jié)果影響不大。
AIC法
AIC即Akaike Information Criterion，赤池信息準(zhǔn)則，考慮了模型的統(tǒng)計(jì)擬合度及用來擬合的參數(shù)數(shù)目 。（不需要嵌套模型）
AIC值越小的模型可優(yōu)行選擇，說明模型用較少的參數(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ù)測變量到模型中，直到添加變量不會(huì)使模型有所改進(jìn)為止。
向后：從模型包含所有預(yù)測變量開始，一次刪除一個(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
一開始包含四個(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ù)測變量）的最佳模型
可通過模型R平方、調(diào)整R平方或Mallows Cp統(tǒng)計(jì)量等準(zhǔn)則來選擇“最佳”模型
調(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")
該圖第一行包含兩個(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ù)挑選出來作為訓(xùn)練樣本，另外的樣本作為保留樣本，先在訓(xùn)練樣本上獲取回歸方程，然后在保留樣本上做預(yù)測。 由于保留樣本不涉及模型及參數(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ù)測方程，記錄k個(gè)保留樣本的預(yù)測表現(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)更重要，最簡單的方法是比較標(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é)
前幾次推文主要介紹了最小二乘回歸，尤其是簡單線性回歸相關(guān)的分析，可以看出構(gòu)建一個(gè)回歸模型有許多注意事項(xiàng)和技巧，當(dāng)然，之前的大多數(shù)方法有一個(gè)前提，即變量服從正態(tài)分布，那么假如我們的變量不服從正態(tài)分布怎么辦呢？因此，之后一段時(shí)間，我們將介紹重抽樣法、自主法與廣義線性回歸。
好了，今天的R語言實(shí)現(xiàn)統(tǒng)計(jì)方法系列推文暫時(shí)告一段落，我們下次再見吧！ 小伙伴們?nèi)绻惺裁唇y(tǒng)計(jì)上的問題，或者如果想要學(xué)習(xí)什么方面的生物信息內(nèi)容，可以在微信群或者知識(shí)星球提問，沒準(zhǔn)哪天的推文就是專門解答你的問題哦！
