各类统计方法R语言实现（八）

【科研绘图点我】【付费精品合集】【SEER点我】

今天是各类统计方法R语言实现的第八期，我们主要介绍选择“最佳”回归模型与深层次分析。

选择“最佳”回归模型

当我们构建回归方程时，我们一方面需要考虑是否去除不显著的变量、是否需要添加交互项和/或多项式项，另一方面还需要权衡方程的简洁性和精度。

因此，我们希望最终构建的“最佳”回归方程，应该至少兼有简洁和有效两个特点。

模型比较方法

方差分析

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     47 289.25                           
## 2     45 289.17  2  0.078505 0.0061 0.9939


p=0.9939，两个模型无显著差异，说明剔除模型1中的Income、Frost对于回归结果影响不大。

AIC法

AIC即Akaike Information Criterion，赤池信息准则，考虑了模型的统计拟合度及用来拟合的参数数目 。（不需要嵌套模型）

AIC值越小的模型可优先选择，说明模型用较少的参数获得了足够的拟合度

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效果更好，与方差分析法结果一致。

变量选择

初步回归法、全子集回归法。

1、逐步回归法

模型每次增加或删除一个变量，直到达到某个判停准则（如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.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


一开始包含四个变量，Population + Illiteracy + Income + Frost，最终得到两个变量Population + Illiteracy；AIC从97.75到95.75再到93.76，逐步降低。

2、全子集回归法

全子集回归，即所有可能的模型都会被检验，可选择展示所有可能的结果，也可展示n个不同变量（一个、两个或多个预测变量）的最佳模型

可通过模型R平方、调整R平方或Mallows Cp统计量等准则来选择“最佳”模型

调整R平方比R平方增加了模型参数数目信息。

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


img

该图第一行包含两个变量Population + Illiteracy，调整后的R平方最高，是最佳模型，和向后回归结果一致。

大部分情况全子集回归法较逐步回归法更优。

深层次分析

评价模型泛化能力和变量相对重要性的方法，包括交叉验证和相对重要性。

1、交叉验证

交叉验证是模型内验证的一种方法，将一定比例的数据挑选出来作为训练样本，另外的样本作为保留样本，先在训练样本上获取回归方程，然后在保留样本上做预测。由于保留样本不涉及模型及参数的选择，该样本可获得比新数据更为精确的估计。

最常用的交叉验证方法是k重交叉验证，在k重交叉验证中，样本被平均分为k个子样本，轮流将k-1个子样本组合作为训练集，另外1个子样本作为保留集，这样会获得k个预测方程，记录k个保留样本的预测表现结果，然后求其平均值。

#10重交叉验证
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平方减少越少，结果越精确，可以看到第二个模型更稳健，泛化能力更强。

2、相对重要性

评价模型中哪个指标更重要，最简单的方法是比较标准化的回归系数。

#z-score标准化
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最不重要。

最小二乘回归小结

前几次推文主要介绍了最小二乘回归，尤其是简单线性回归相关的分析，可以看出构建一个回归模型有许多注意事项和技巧，当然，之前的大多数方法有一个前提，即变量服从正态分布，那么假如我们的变量不服从正态分布怎么办呢？因此，之后一段时间，我们将介绍重抽样法、自主法与广义线性回归。

好了，今天的R语言实现统计方法系列推文暂时告一段落，我们下次再见吧！小伙伴们如果有什么统计上的问题，或者如果想要学习什么方面的生物信息内容，可以在微信群或者知识星球提问，没准哪天的推文就是专门解答你的问题哦！