R语言统计分析——回归模型深层次分析

本文主要是介绍R语言统计分析——回归模型深层次分析，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

参考资料：R语言实战【第2版】

本文主要讨论回归模型的泛化能力和变量相对重要性的方法。

1、交叉验证

从定义上看，回归方法就是从一堆数据中获取最优模型参数。对于OLS（普通最小二乘）回归，通过使得预测误差（残差）平方和最小和对响应变量的解释度（R平方）最大，可获得模型参数。由于等式只是最优化已给出的数据，所以在新数据集上表现并不一定好。

通过交叉验证法，我们可以评价回归方程的泛化能力。所谓交叉验证，就是将一定比例的数据挑选出来作为驯良样本，另外的样本作为保留样本，先在训练样本上获取回归方程，然后在保留样本上做预测。由于保留样本不涉及模型参数的选择，该样本可获得比新数据更为精确的估计。

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

bootstrap包中的crossval()函数可以实现k重交叉验证。如下：

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)^2r2cv<-cor(y,results$cv.fit)^2cat("Original R-square=",r2,"\n")cat(k,"Fold Cross-Validated R-square=",r2cv,"\n")cat("Change=",r2-r2cv,"\n")
}
states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
fit<-lm(Murder~Population+Income+Illiteracy+Frost,data=states)
shrinkage(fit)

有结果可知，基于初始用样本的R平法（0.567）过于乐观，对新数据更好的方差解释率估计是交叉验证后的R平方（0.476）。（注意，由于观测被随机分配到k个群组中，因此每次运行shrinkage()函数，得到的结果都会有些许不同）

2、相对重要性

我们根据线性回归模型做出来预测方程后，还有一个问题需要关注：哪些变量对于我们的预测来说更为重要？

若预测变量不相关，过程就像对简单得多，我们可以根据预测变量与响应变量的相关系数来进行排序。但大部分情况中，预测变量之间有一定相关性，这就使得评价变得复杂很多。

评价预测变量的相对重要性最简单的方法就是比较标准化的回归系数，它表示当其他预测变量不变时，该预测变量一个标准差的变化可引起的响应变量的预期变化（以标准差单位度量）。在进行回归分析前，可用scale()函数将数据标准化为均值为0、标准差为1的数据集，这样用R回归即可获得标准化的回归系数。（注意，scale()函数返回的是一个矩阵，而lm()函数要求的是一个数据框，我们需要一个中间步骤来转换一下。）如下：

states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
# 对数据进行标准化
zstates<-as.data.frame(scale(states))
# 查看标准化后的数据
head(zstates)
# 对标准化后的数据进行拟合
zfit<-lm(Murder~Population+Income+Illiteracy+Frost,data=zstates)
# 查看回归系数
coef(zfit)

有结果可知：当其他以因素不变时，Illiteracy一个标准差的变化将增加0.68个标准差的谋杀率。根据标准化的回归系数，我们可以认为Illiteracy是最重要的预测变量。

还有其他方法可定量分析预测变量的相对重要性。比如，可以将相对重要性看作每个预测变量（本身或与其他预测变量组合）对R平方的贡献。相对权重是对所有可能子模型添加一个预测变量引起的R平方平均增加量的一个近似值。如下：

# 编写相对权重函数
relweights<-function(fit,...){R <- cor(fit$model)nvar <- ncol(R)rxx <- R[2:nvar, 2:nvar]rxy <- R[2:nvar, 1]svd <- eigen(rxx)evec <- svd$vectorsev <- svd$valuesdelta <- diag(sqrt(ev))lambda <- evec %*% delta %*% t(evec)lambdasq <- lambda^2beta <- solve(lambda) %*% rxyrsquare <- colSums(beta ^ 2)rawwgt <- lambdasq %*% beta ^ 2import <- (rawwgt / rsquare) * 100import <- as.data.frame(import)row.names(import) <- names(fit$model[2:nvar])names(import) <- "Weights"dotchart(import$Weights, labels=row.names(import),xlab="% of R-Square", pch=19,main="Relative Importance of Predictor Variables",sub=paste("Total R-Square=", round(rsquare, digits=3)),...)return(import)
}
# 应用相对权重函数
states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
relweights(fit,col="blue")