从许多方面来看，回归分析都是统计学的核心。它其实是一个广义的概念，通指那些用一个或多个预测变量（也称自变量或解释变量）来预测响应变量（也称因变量、效标变量或结果变量）的方法。通常，回归分析可以用来挑选与响应变量相关的解释变量，可以描述两者的关系，也可以生成一个等式，通过解释变量来预测响应变量。本文主要讲述了线性回归的问题，详细描述了线性回归的建模、模型诊断、模型操作及模型可视化等。

《R语言实战--第2版》---------Robert I.Kabacoff

本文内容：

• 回归的多面性
• OLS回归
• 回归诊断
• 异常观测值
• 改进措施

1. 回归的多面性

对于回归模型的拟合，R提供的强大而丰富的功能和选项也同样令人困惑。

在本次分享中，我们的重点是普通最小二乘（OLS）回归法，包括简单线性回归、多项式回归和多元线性回归。

2. OLS回归

OLS回归拟合模型的形式：

我们的目标是通过减少响应变量的真实值与预测值的差值来获得模型参数 （截距项和斜率） 。具体而言，即使得残差平方和最小。

观蒙斋主人：回归分析的基本假设​zhuanlan.zhihu.com

机器学习---最小二乘线性回归模型的5个基本假设（Machine Learning Least Squares Linear Regression Assumptions）​www.cnblogs.com

如果违背了以上假设， 你的统计显著性检验结果和所得的置信区间就很可能不精确了。 注意，OLS回归还假定自变量是固定的且测量无误差，但在实践中通常都放松了这个假设。

2.1 用 lm()拟合回归模型

在R中，拟合线性模型最基本的函数就是lm()，格式为：

myfit <- lm(formula, data)

其中，formula指要拟合的模型形式，data是一个数据框，包含了用于拟合模型的数据。结果对象（本例中是myfit）存储在一个列表中，包含了所拟合模型的大量信息。

表达式（formula）形式如下：

Y ~ X1 + X2 + … + Xk

~左边为响应变量，右边为各个预测变量，预测变量之间用+符号分隔。

除了lm()，表8-3还列出了其他一些对做简单或多元回归分析有用的函数。拟合模型后，将这些函数应用于lm()返回的对象，可以得到更多额外的模型信息。

当回归模型包含一个因变量和一个自变量时， 我们称为简单线性回归。

当只有一个预测变量，但同时包含变量的幂（比如，

）时，我们称为多项式回归。当有不止一个预测变量时，则称为多元线性回归。

2.2 简单线性回归

数据集women提供了15个年龄在30~39岁间女性的身高和体重信息，我们想通过身高来预测体重，获得一个等式可以帮助我们分辨出那些过重或过轻的个体。

fit 

通过输出结果，可以得到预测等式：

因为身高不可能为0， 所以没必要给截距项一个物理解释， 它仅仅是一个常量调整项。 在Pr(>|t|)栏，可以看到回归系数（3.45）显著不为0（p<0.001） ，表明身高每增高1英寸，体重将预期增加3.45磅。R平方项（0.991）表明模型可以解释体重99.1%的方差，它也是实际和预测值之间相关系数的平方（$R^2=r^2_{hat{y}y}$） 。残差标准误（1.53 lbs）则可认为是模型用身高预测体重的平均误差。F统计量检验所有的预测变量预测响应变量是否都在某个几率水平之上。由于简单回归只有一个预测变量，此处F检验等同于身高回归系数的t检验。

2.3多项式回归

fit2 <- lm(weight~height+I(height^2),data=women)
summary(fit2)
plot(women$height,women$weight, xlab="Height (in inches)", ylab="Weight (in lbs)")
lines(women$height,fitted(fit2))

新的预测等式为：

在p<0.001水平下，回归系数都非常显著。模型的方差解释率已经增加到了99.9%。二次项的显著性（t=13.89，p<0.001）表明包含二次项提高了模型的拟合度。从图也可以看出曲线确实拟合得较好。

那线性回归的非线性回归的区别是什么呢？线性模型是指参数线性而不是x线性，同理，非线性模型是指参数非线性而非x非线性，如：

car包中的scatterplot()函数，它可以很容易、方便地绘制二元关系图:

library(car)
scatterplot(weight~height,data=women,spread=FALSE, smoother.args=list(lty=2), pch=19, main="Women Age 30-39", xlab="Height (inches)", ylab="Weight (lbs.)")

这个功能加强的图形，既提供了身高与体重的散点图、线性拟合曲线和平滑拟合（loess）曲线，还在相应边界展示了每个变量的箱线图。spread=FALSE选项删除了残差正负均方根在平滑曲线上的展开和非对称信息。smoother.args=list(lty=2)选项设置loess拟合曲线为虚线。pch=19选项设置点为实心圆（默认为空心圆）。粗略地看一下图可知，两个变量基本对称，曲线拟合得比直线更好。

2.4 多元线性回归

当预测变量不止一个时，简单线性回归就变成了多元线性回归，分析也稍微复杂些。以基础包中的state.x77数据集为例， 我们想探究一个州的犯罪率和其他因素的关系， 包括人口、文盲率、平均收入和结霜天数（温度在冰点以下的平均天数） 。多元回归分析中，第一步最好检查一下变量间的相关性。cor()函数提供了二变量之间的相关系数， car包scatterplotMatrix()函数则会生成散点图矩阵 ：

states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
cor(states)
library(car)
scatterplotMatrix(states, spread=FALSE, smoother.args=list(lty=2), main="Scatter Plot Matrix")

scatterplotMatrix()函数默认在非对角线区域绘制变量间的散点图，并添加平滑和线性拟合曲线。对角线区域绘制每个变量的密度图和轴须图。从图中可以看到，谋杀率是双峰的曲线，每个预测变量都一定程度上出现了偏斜。谋杀率随着人口和文盲率的增加而增加，随着收入水平和结霜天数增加而下降。同时，越冷的州府文盲率越低，收入水平越高。

states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
fit <- lm(Murder~Population+Illiteracy+Income+Frost,data=states)
summary(fit)

当预测变量不止一个时，回归系数的含义为：一个预测变量增加一个单位，其他预测变量保持不变时，因变量将要增加的数量。例如本例中，文盲率的回归系数为4.14，表示控制人口、收入和温度不变时，文盲率上升1%，谋杀率将会上升4.14%，它的系数在p<0.001的水平下显著不为0。相反，Frost的系数没有显著不为0（p=0.954），表明当控制其他变量不变时，Frost与Murder不呈线性相关。总体来看，所有的预测变量解释了各州谋杀率57%的方差。

2.5 有交互项的多元线性回归

以mtcars数据框中的汽车数据为例，若你对汽车重量和马力感兴趣，可以把它们作为预测变量，并包含交互项来拟合回归模型。

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

可以看到Pr(>|t|)栏中，马力与车重的交互项是显著的，这意味着什么呢？若两个预测变量的交互项显著，说明响应变量与其中一个预测变量的关系依赖于另外一个预测变量的水平。因此此例说明，每加仑汽油行驶英里数与汽车马力的关系依车重不同而不同。

通过effects包中的effect()函数，你可以用图形展示交互项的结果。格式为：

plot(effect(term, mod,, xlevels), multiline=TRUE)

term即模型要画的项，mod为通过lm()拟合的模型，xlevels是一个列表，指定变量要设定的常量值，multiline=TRUE选项表示添加相应直线。

library(effects)
plot(effect("hp:wt",fit,,list(wt=c(2.2,3.2,4.2))),multiline=TRUE)

从图中可以很清晰地看出，随着车重的增加，马力与每加仑汽油行驶英里数的关系减弱了。当wt=4.2时，直线几乎是水平的，表明随着hp的增加，mpg不会发生改变。

3.回归诊断

# 3.回归诊断

使用lm()函数来拟合OLS回归模型，通过summary()函数获取模型参数和相关统计量。但是，没有任何输出告诉你模型是否合适，你对模型参数推断的信心依赖于它在多大程度上满足OLS模型统计假设。

为什么这很重要？因为数据的无规律性或者错误设定了预测变量与响应变量的关系，都将致使你的模型产生巨大的偏差。一方面，你可能得出某个预测变量与响应变量无关的结论，但事实上它们是相关的；另一方面，情况可能恰好相反。当你的模型应用到真实世界中时，预测效果可能很差，误差显著。

states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
fit <- lm(Murder~Population+Illiteracy+Income+Frost,data=states)
confint(fit)

因为Frost的置信区间包含0，所以貌似可以得出结论：当其他变量不变时，温度的改变与谋杀率无关。(不一定对)不过，你对这些结果的信念，都只建立在你的数据满足统计假设的前提之上。 回归诊断技术向你提供了评价回归模型适用性的必要工具， 它能帮助发现并纠正问题。

3.1 标准方法

R基础安装中提供了大量检验回归分析中统计假设的方法。最常见的方法就是对lm()函数返回的对象使用plot()函数，可以生成评价模型拟合情况的四幅图形。

fit <- lm(weight ~ height, data=women)
par(mfrow=c(2,2))
plot(fit) 

让我们再看看二次拟合的诊断图：

fit2 <- lm(weight ~ height + I(height^2), data=women)
par(mfrow=c(2,2))
plot(fit2)

这第二组图表明多项式回归拟合效果比较理想，基本符合了线性假设、残差正态性（除了观测点13）和同方差性（残差方差不变）。观测点15看起来像是强影响点（根据是它有较大的Cook距离值） ，删除它将会影响参数的估计。事实上，删除观测点13和15，模型会拟合得会更好。

newfit <- lm(weight~ height + I(height^2), data=women[-c(13,15),])
par(mfrow=c(2,2))
plot(newfit)

3.2 改进的方法

car包提供了大量函数，大大增强了拟合和评价回归模型的能力。

(1) 正态性

与基础包中的plot()函数相比，qqPlot()函数提供了更为精确的正态假设检验方法，它画出了在n–p–1个自由度的t分布下的学生化残差（studentized-residual，也称学生化删除残差或折叠化残差）图形，其中n是样本大小，p是回归参数的数目（包括截距项）。

library(car)
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
qqPlot(fit, labels=row.names(states), id.method="identify", simulate=TRUE, main="Q-Q Plot")

qqPlot()函数生成的概率图见图8-9。id.method=“identify”选项能够交互式绘图——待图形绘制后，用鼠标单击图形内的点，将会标注函数中labels选项的设定值。敲击Esc键，从图形下拉菜单中选择Stop， 或者在图形上右击， 都将关闭这种交互模式。 此处， 我已经鉴定出了Nevada异常。当simulate=TRUE时，95%的置信区间将会用参数自助法。

除了Nevada， 所有的点都离直线很近， 并都落在置信区间内， 这表明正态性假设符合得很好。但是你也必须关注Nevada，它有一个很大的正残差值（真实值－预测值） ，表明模型低估了该州的谋杀率。

states["Nevada",]
fitted(fit)["Nevada"]
residuals(fit)["Nevada"]
rstudent(fit)["Nevada"]

可视化误差还有其他方法，residplot()函数生成学生化残差柱状图（即直方图） ，并添加正态曲线、核密度曲线和轴须图。

residplot <- function(fit, nbreaks=10) {  z <- rstudent(fit) hist(z, breaks=nbreaks, freq=FALSE, xlab="Studentized Residual", main="Distribution of Errors") rug(jitter(z), col="brown") curve(dnorm(x, mean=mean(z), sd=sd(z)), add=TRUE, col="blue", lwd=2) lines(density(z)$x, density(z)$y, col="red", lwd=2, lty=2) legend("topright", legend = c( "Normal Curve", "Kernel Density Curve"), lty=1:2, col=c("blue","red"), cex=.7) } residplot(fit) 

正如你所看到的，除了一个很明显的离群点，误差很好地服从了正态分布。虽然Q-Q图已经蕴藏了很多信息，但我总觉得从一个柱状图或者密度图测量分布的斜度比使用概率图更容易。因此为何不一起使用这两幅图呢？

(2) 误差的独立性

判断因变量值（或残差）是否相互独立，最好的方法是依据收集数据方式的先验知识。例如，时间序列数据通常呈现自相关性——相隔时间越近的观测相关性大于相隔越远的观测。car包提供了一个可做Durbin-Watson检验的函数，能够检测误差的序列相关性。

durbinWatsonTest(fit)  

p值不显著（p=0.282）说明无自相关性，误差项之间独立。滞后项（lag=1）表明数据集中每个数据都是与其后一个数据进行比较的。该检验适用于时间独立的数据，对于非聚集型的数据并不适用。

(3) 线性

通过成分残差图（component plus residual plot）也称偏残差图（partial-residual-plot），你可以看看因变量与自变量之间是否呈非线性关系，也可以看看是否有不同于已设定线性模型的系统偏差，图形可用car包中的crPlots()函数绘制。

创建变量X的成分残差图，每幅图都是

的直线。

library(car)
crPlots(fit) 

若图形存在非线性， 则说明你可能对预测变量的函数形式建模不够充分，那么就需要添加一些曲线成分，比如多项式项，或对一个或多个变量进行变换（如用log(X)代替X），或用其他回归变体形式而不是线性回归。

(4) 同方差性

car包提供了两个有用的函数，可以判断误差方差是否恒定。ncvTest()函数生成一个计分检验， 零假设为误差方差不变， 备择假设为误差方差随着拟合值水平的变化而变化。 若检验显著，则说明存在异方差性（误差方差不恒定）。spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系。

library(car)
ncvTest(fit)
spreadLevelPlot(fit)

可以看到，计分检验不显著（p=0.19） ，说明满足方差不变假设。你还可以通过分布水平图（图8-12）看到这一点，其中的点在水平的最佳拟合曲线周围呈水平随机分布。若违反了该假设，你将会看到一个非水平的曲线。代码结果建议幂次变换（suggested-power-transformation）的含义是，经过p次幂（

）变换，非恒定的误差方差将会平稳。例如，若图形显示出了非水平趋势，建议幂次转换为0.5，在回归等式中用

代替Y，可能会使模型满足同方差性。若建议幂次为0，则使用对数变换。对于当前例子，异方差性很不明显，因此建议幂次接近1（不需要进行变换）。

3.3 线性模型假设的综合验证

最后，让我们一起学习gvlma包中的gvlma()函数。gvlma()函数由Pena和Slate（2006）编写，能对线性模型假设进行综合验证，同时还能做偏斜度、峰度和异方差性的评价。换句话说，它给模型假设提供了一个单独的综合检验（通过/不通过） 。

library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)

从输出项（Global Stat中的文字栏）我们可以看到数据满足OLS回归模型所有的统计假设（p=0.597） 。若Decision下的文字表明违反了假设条件（比如p<0.05） ，你可以使用前几节讨论的方法来判断哪些假设没有被满足。

3.4 多重共线性

假设你正在进行一项握力研究，自变量包括DOB（Date Of Birth，出生日期）和年龄。你用握力对DOB和年龄进行回归，F检验显著，p<0.001。但是当你观察DOB和年龄的回归系数时，却发现它们都不显著（也就是说无法证明它们与握力相关） 。到底发生了什么呢？ 原因是DOB与年龄在四舍五入后相关性极大。回归系数测量的是当其他预测变量不变时，某个预测变量对响应变量的影响。那么此处就相当于假定年龄不变，然后测量握力与年龄的关系，这种问题就称作多重共线性（multicollinearity）。它会导致模型参数的置信区间过大，使单个系数解释起来很困难。多重共线性可用统计量VIF（V ariance Inflation Factor，方差膨胀因子）进行检测。VIF的平方根表示变量回归参数的置信区间能膨胀为与模型无关的预测变量的程度（因此而得名） 。car包中的vif()函数提供VIF值。一般原则下，

就表明存在多重共线性问题。

library(car)
vif(fit)
sqrt(vif(fit)) >2

结果表明预测变量不存在多重共线性问题。

4.异常观测值

一个全面的回归分析要覆盖对异常值的分析，包括离群点、高杠杆值点和强影响点。这些数据点需要更深入的研究，因为它们在一定程度上与其他观测点不同，可能对结果产生较大的负面影响。

4.1 离群点

离群点是指那些模型预测效果不佳的观测点。它们通常有很大的、或正或负的残差$Y_i-hat{Y}_i$,正的残差说明模型低估了响应值，负的残差则说明高估了响应值。

Q-Q图， 落在置信区间带外的点即可被认为是离群点。另外一个粗糙的判断准则：标准化残差值大于2或者小于–2的点可能是离群点，需要特别关注。

car包也提供了一种离群点的统计检验方法。 outlierTest()函数可以求得最大标准化残差绝对值Bonferroni调整后的p值：

library(car)
outlierTest(fit)

此处，你可以看到Nevada被判定为离群点（p=0.048）。注意，该函数只是根据单个最大（或正或负）残差值的显著性来判断是否有离群点。若不显著，则说明数据集中没有离群点； 若显著，则你必须删除该离群点，然后再检验是否还有其他离群点存在。

4.2 高杠杆点

高杠杆值观测点，即与其他预测变量有关的离群点。换句话说，它们是由许多异常的预测变量值组合起来的，与响应变量值没有关系。

高杠杆值的观测点可通过帽子统计量（hat statistic）判断。对于一个给定的数据集，帽子均值为p/n，其中p是模型估计的参数数目（包含截距项） ，n是样本量。一般来说，若观测点的帽子值大于帽子均值的2或3倍，就可以认定为高杠杆值点。

hat.plot <- function(fit) {  p <- length(coefficients(fit)) n <- length(fitted(fit)) plot(hatvalues(fit), main="Index Plot of Hat Values") abline(h=c(2,3)*p/n, col="red", lty=2) identify(1:n, hatvalues(fit), names(hatvalues(fit))) }
hat.plot(fit) 

水平线标注的即帽子均值2倍和3倍的位置。此图中，可以看到Alaska和California非常异常，查看它们的预测变量值，与其他48个州进行比较发现： Alaska收入比其他州高得多， 而人口和温度却很低； California人口比其他州府多得多，但收入和温度也很高。高杠杆值点可能是强影响点，也可能不是，这要看它们是否是离群点。

4.3 强影响点

强影响点，即对模型参数估计值影响有些比例失衡的点。例如，若移除模型的一个观测点时模型会发生巨大的改变，那么你就需要检测一下数据中是否存在强影响点了。

有两种方法可以检测强影响点：Cook距离，或称D统计量，以及变量添加图（added variable plot）。一般来说，Cook’sD值大于4/(n–k–1)，则表明它是强影响点，其中n为样本量大小，k是预测变量数目。

cutoff <- 4/(nrow(states)-length(fit$coefficients)-2)
plot(fit, which=4, cook.levels=cutoff)
abline(h=cutoff, lty=2, col="red")

通过图形可以判断Alaska、Hawaii和Nevada是强影响点。若删除这些点，将会导致回归模型截距项和斜率发生显著变化。注意，该图对搜寻强影响点很有用。Cook’sD图有助于鉴别强影响点，但是并不提供关于这些点如何影响模型的信息。变量添加图弥补了这个缺陷。对于一个响应变量和k个预测变量，你可以如下图创建k个变量添加图。所谓变量添加图，即对于每个预测变量，绘制

在其他k–1个预测变量上回归的残差值相对于响应变量在其他k–1个预测变量上回归的残差值的关系图。car包中的avPlots()函数可提供变量添加图：

library(car)
avPlots(fit, ask=FALSE, id.method="identify")

图中的直线表示相应预测变量的实际回归系数。你可以想象删除某些强影响点后直线的改变，以此来估计它的影响效果。例如，来看左下角的图（“Murder | others” vs. “Income | others”） ，若删除点Alaska，直线将往负向移动。事实上，删除Alaska，Income的回归系数将会从0.000 06变为–0.000 85。利用car包中的influencePlot()函数，你还可以将离群点、杠杆值和强影响点的信息整合到一幅图形中：

library(car)
influencePlot(fit, id.method="identify", main="Influence Plot", sub="Circle size is proportional to Cook's distance")  

纵坐标超过+2或小于–2的州可被认为是离群点，水平轴超过0.2或0.3的州有高杠杆值（通常为预测值的组合）。圆圈大小与影响成比例，圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影响点。图反映出Nevada和Rhode Island是离群点，New York、California、Hawaii和Washington有高杠杆值，Nevada、Alaska和Hawaii为强影响点。

5.改进措施

我们已经花费了不少篇幅来学习回归诊断，你可能会问： “如果发现了问题，那么能做些什么呢？”有四种方法可以处理违背回归假设的问题：

• 删除观测点；
• 变量变换；
• 添加或删除变量；
• 使用其他回归方法。

5.1 删除观测点

删除离群点通常可以提高数据集对于正态假设的拟合度，而强影响点会干扰结果，通常也会被删除。删除最大的离群点或者强影响点后，模型需要重新拟合。若离群点或强影响点仍然存在，重复以上过程直至获得比较满意的拟合。

不过在其他情况下，所收集数据中的异常点可能是最有趣的东西。发掘为何该观测点不同于其他点，有助于你更深刻地理解研究的主题，或者发现其他你可能没有想过的问题。

5.2 变量变换

当模型不符合正态性、线性或者同方差性假设时，一个或多个变量的变换通常可以改善或调整模型效果。

car包中的powerTransform()函数通过λ的最大似然估计来正态化变量

。

library(car)
summary(powerTransform(states$Murder))

结果表明，你可以用Murder0.6来正态化变量Murder。由于0.6很接近0.5，你可以尝试用平方根变换来提高模型正态性的符合程度。但在本例中，

的假设也无法拒绝（p=0.145） ，因此没有强有力的证据表明本例需要变量变换，这与图8-9的Q-Q图结果一致。

当违反了线性假设时，对预测变量进行变换常常会比较有用。car包中的boxTidwell()函数通过获得预测变量幂数的最大似然估计来改善线性关系。

library(car)
boxTidwell(Murder~Population+Illiteracy,data=states)

结果显示，使用变换Population0.87和Illiteracy1.36能够大大改善线性关系。但是对Population（p=0.75）和Illiteracy（p=0.54）的计分检验又表明变量并不需要变换。与成分残差图是一致的。

5.3 增删变量

删除变量在处理多重共线性时是一种非常重要的方法。如果你仅仅是做预测，那么多重共线性并不构成问题，但是如果还要对每个预测变量进行解释，那么就必须解决这个问题。最常见的方法就是删除某个存在多重共线性的变量（某个变量

） 。另外一个可用的方法便是岭回归——多元回归的变体，专门用来处理多重共线性问题。

5.4 尝试其他方法

正如刚才提到的，处理多重共线性的一种方法是拟合一种不同类型的模型（本例中是岭回归）。其实，如果存在离群点和/或强影响点，可以使用稳健回归模型替代OLS回归。如果违背了正态性假设，可以使用非参数回归模型。如果存在显著的非线性，能尝试非线性回归模型。如果违背了误差独立性假设，还能用那些专门研究误差结构的模型，比如时间序列模型或者多层次回归模型。最后，你还能转向广泛应用的广义线性模型，它能适用于许多OLS回归假设不成立的情况。

5.5 变量选择

从大量候选变量中选择最终的预测变量有以下两种流行的方法：逐步回归法（stepwise method）和全子集回归（all-subsets regression） 。

(1) 逐步回归

逐步回归中，模型会一次添加或者删除一个变量，直到达到某个判停准则为止。例如，向前逐步回归（forward-stepwise-regression）每次添加一个预测变量到模型中，直到添加变量不会使模型有所改进为止。 向后逐步回归（backward-stepwise-regression）从模型包含所有预测变量开始，一次删除一个变量直到会降低模型质量为止。而向前向后逐步回归（stepwise stepwise regression，通常称作逐步回归，以避免听起来太冗长），结合了向前逐步回归和向后逐步回归的方法，变量每次进入一个，但是每一步中，变量都会被重新评价，对模型没有贡献的变量将会被删除，预测变量可能会被添加、删除好几次，直到获得最优模型为止。

library(MASS)
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
stepAIC(fit)

开始时模型包含4个（全部）预测变量，然后每一步中，AIC列提供了删除一个行中变量后模型的AIC值，<none>中的AIC值表示没有变量被删除时模型的AIC。第一步，Frost被删除，AIC从97.75降低到95.75；第二步，Income被删除，AIC继续下降，成为93.76。然后再删除变量将会增加AIC，因此终止选择过程。

逐步回归法其实存在争议，虽然它可能会找到一个好的模型，但是不能保证模型就是最佳模型，因为不是每一个可能的模型都被评价了。为克服这个限制，便有了全子集回归法。

(2) 全子集回归法

全子集回归是指所有可能的模型都会被检验。全子集回归可用leaps包中的regsubsets()函数实现。你能通过R平方、调整R平方或Mallows Cp统计量等准则来选择“最佳”模型。

library(leaps)
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) leaps <-regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data=states, nvmax =4)
summary.leaps <- summary(leaps)
names(summary.leaps)
summary.leaps$bic
par(mfrow=c(2,2))
plot(summary.leaps$rss ,xlab="Number of Variables ",ylab="RSS",type="l")
plot(summary.leaps$adjr2 ,xlab="Number of Variables ",ylab="Adjusted RSq",type="l")
points (2,summary.leaps$adjr2[2], col="red",cex=2,pch=20)
plot(summary.leaps$cp ,xlab="Number of Variables ",ylab="Cp",type='l')
points (2,summary.leaps$cp [2],col="red",cex=2,pch=20)
plot(summary.leaps$bic ,xlab="Number of Variables ",ylab="BIC",type='l')
points(2,summary.leaps$bic [2],col="red",cex=2,pch =20)
coef(leaps ,2)

5.6 相对重要性

我们一直都有一个疑问： “哪些变量对预测有用呢？”但你内心真正感兴趣的其实是：“哪些变量对预测最为重要？”潜台词就是想根据相对重要性对预测变量进行排序。这个问题很有实际用处。例如，假设你能对团队组织成功所需的领导特质依据相对重要性进行排序，那么就可以帮助管理者关注他们最需要改进的行为。若预测变量不相关，过程就相对简单得多，你可以根据预测变量与响应变量的相关系数来进行排序。但大部分情况中，预测变量之间有一定相关性，这就使得评价变得复杂很多。评价预测变量相对重要性的方法一直在涌现。

(1) 标准化回归系数

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

states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
zstates <- as.data.frame(scale(states))
zfit <- lm(Murder~Population + Income + Illiteracy + Frost, data=zstates)
coef(zfit)

此处可以看到， 当其他因素不变时， 文盲率一个标准差的变化将增加0.68个标准差的谋杀率。根据标准化的回归系数，我们可认为Illiteracy是最重要的预测变量，而Frost是最不重要的。

(2) 相对权重法

相对权重（relative weight）是一种比较有前景的新方法，它是对所有可能子模型添加一个预测变量引起的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$vectors ev <- svd$values delta <- diag(sqrt(ev)) lambda <- evec %*% delta %*% t(evec) lambdasq <- lambda ^ 2 beta <- solve(lambda) %*% rxy rsquare <- colSums(beta ^ 2) rawwgt <- lambdasq %*% beta ^ 2 import <- (rawwgt / rsquare) * 100 import <- as.data.frame(import) row.names(import) <- names(fit$model[2:nvar]) names(import) <- "Weights" import <- import[order(import),1, drop=FALSE] 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") 

可以看到各个预测变量对模型方差的解释程度（R平方=0.567），Illiteracy解释了59%的R平方，Frost解释了20.79%，等等。根据相对权重法，Illiteracy有最大的相对重要性，余下依次是Frost、Population和Income。相对重要性的测量（特别是相对权重方法）有广泛的应用，它比标准化回归系数更为直观。

6.结语

在统计中， 回归分析是许多方法的一个总称。相信你已经看到，它是一个交互性很强的方法，包括拟合模型、检验统计假设、修正数据和模型，以及为达到最终结果的再拟合等过程。从很多角度来看，获得模型的最终结果不仅是一种科学，也是一种艺术和技巧。