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 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 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 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個子樣本,輪流將k1個子樣本組合作為訓練集,另外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") }