8.6 選擇“最佳”的回歸模型

8.6.1 模型比較

用基礎安裝中的anova()函數可以比較兩個嵌套模型的擬合優度。所謂嵌套模型，即它的一

些項完全包含在另一個模型中

用anova()函數比較

> 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 47289.25

2 45289.17 2 0.078505 0.0061 0.9939

AIC（AkaikeInformation Criterion，赤池信息準則）也可以用來比較模型，它考慮了模型的

統計擬合度以及用來擬合的參數數目。AIC值越小的模型要優先選擇，它說明模型用較少的參數

獲得了足夠的擬合度。

> AIC(fit1,fit2)

df AIC

fit1 6 241.6429

fit2 4 237.6565

8.6.2變量選擇

1. 逐步回歸stepwise method

逐步回歸中，模型會一次添加或者刪除一個變量，直到達到某個判停準則為止。向前

逐步回歸（forward stepwise）每次添加一個預測變量到模型中，直到添加變量不會使模型有所改

進為止。向后逐步回歸（backward stepwise）從模型包含所有預測變量開始，一次刪除一個變量

直到會降低模型質量為止。而向前向后逐步回歸（stepwise stepwise，通常稱作逐步回歸

），結合了向前逐步回歸和向后逐步回歸的方法，變量每次進入一個，但是每一步

中，變量都會被重新評價，對模型沒有貢獻的變量將會被刪除，預測變量可能會被添加、刪除好

幾次，直到獲得最優模型為止。。MASS包中的stepAIC()函數可以實現

逐步回歸模型（向前、向后和向前向后），依據的是精確AIC準則。

> 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.658332.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.646588.89 127.311

Call:

lm(formula = Murder ~Population + Illiteracy, data = states)

Coefficients:

(Intercept) Population Illiteracy

1.6515497 0.0002242 4.0807366

2. 全子集回歸

全子集回歸可用leaps包中的regsubsets()函數實現。你能通過R平方、調整R平方或

Mallows Cp統計量等準則來選擇“最佳”模型

> library("leaps", lib.loc="d:/ProgramFiles/R/R-3.1.3/library")

>leaps<-regsubsets(Murder~Population+Illiteracy+Income+Frost,data=states,nbest=4)

> plot(leaps,scal="adjr2")

> library(car)

> subsets(leaps,statistic="cp",main="cpplot for all subsets regression")

> abline(1,1,lty=2,col="red")

8.7 深層次分析

8.7.1 交叉驗證

所謂交叉驗證，即將一定比例的數據挑選出來作為訓練樣本，另外的樣本作保留樣本，先在

訓練樣本上獲取回歸方程，然后在保留樣本上做預測。由于保留樣本不涉及模型參數的選擇，該

樣本可獲得比新數據更為精確的估計。在k 重交叉驗證中，樣本被分為k個子樣本，輪流將k?1個子樣本組合作為訓練集，另外1個子樣本作為保留集。這樣會獲得k 個預測方程，記錄k 個保留樣本的預測表現結果，然后求其平均值。[當n 是觀測總數目，k 為n 時，該方法又稱作刀切法（jackknifing）]bootstrap 包中的crossval() 函數可以實現k 重交叉驗證。

fit<-lm(mpg~hp+wt+hp:wt,data=mtcars)

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")

}