本章内容
拟合并解释线性模型
检验模型假设
模型选择
8.1 回归的多面性
8.2 OLS回归
基本介绍：
8.2.1 用lm（）拟合回归模型
格式
myfit <- lm(formula,data)
表达式（formula） 形式如下： Y ~ X1 + X2 + ... + Xk
~左边为响应变量，右边为各个预测变量，预测变量之间用+符号分隔。
R表达式中常用的符号
应用于lm（）返回的对象：
8.2.2 简单线性回归
fit <- lm(weight ~ height, data=women) # women是R中的基础数据集
summary(fit)
women$weight
fitted(fit) #给出预测值
residuals(fit)
plot(women$height,women$weight,
xlab="Height (in inches)",
ylab="Weight (in pounds)")
abline(fit) #画出拟合直线
8.2.3 多项式回归：
添加二次项：
fit2 <- lm(weight ~ height + I(height^2), data=women)
I(height^2)表示向预测等式添加一个身高的平方项。I函数将括号的内容看作R的一个常规表 达式。
图像上看也拟合得更好
线性模型与非线性模型
一般来说，n次多项式生成一个n–1个弯曲的曲线。
拟合三次多项式，可用：
fit3 <- lm(weight ~ height + I(height^2) +I(height^3), data=women)
car包中的scatterplot()函数，它可以很容易、方便地绘制 二元关系图。以下代码能生成图8-3所示的图形：
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选项设置点为实心圆（默认为空心圆）。
粗略地看一下图8-3可知，两个变量基本对称， 曲线拟合得比直线更好。
8.2.4多元线性回归
以基础包中的state.x77数据集为例，我们想探究一个州的犯罪率和其他因素的关系，包括 人口、文盲率、平均收入和结霜天数（温度在冰点以下的平均天数）。
因为lm()函数需要一个数据框（state.x77数据集是矩阵），为了以后处理方便，你需要做 如下转化：
states <- as.data.frame(state.x77[,c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
这行代码创建了一个名为states的数据框，包含了我们感兴趣的变量。
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)
当预测变量不止一个时，回归系数的含义为：一个预测变量增加一个单位，其他预测变量保 持不变时，因变量将要增加的数量。
8.2.5 有交互项的多元线性回归
例子
fit <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
summary(fit)
若两个预测 变量的交互项显著，说明响应变量与其中一个预测变量的关系依赖于另外一个预测变量的水平。 因此此例说明，每加仑汽油行驶英里数与汽车马力的关系依车重不同而不同。
通过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不会发生改变。
8.3 回归诊断
在上一节中，你使用lm()函数来拟合OLS回归模型，通过summary()函数获取模型参数和 相关统计量。但是，没有任何输出告诉你模型是否合适，你对模型参数推断的信心依赖于它在多 大程度上满足OLS模型统计假设。虽然在代码清单8-4中summary()函数对模型有了整体的描述， 但是它没有提供关于模型在多大程度上满足统计假设的任何信息。
为什么这很重要？因为数据的无规律性或者错误设定了预测变量与响应变量的关系，都将致 使你的模型产生巨大的偏差。
一方面，你可能得出某个预测变量与响应变量无关的结论，但事实 上它们是相关的另一方面，情况可能恰好相反。当你的模型应用到真实世界中时，预测效果可 能很差，误差显著。
states <- as.data.frame(state.x77[,c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
confint(fit)
结论都只建立在你的数据满足统计假设的前提之上。回归诊断技术向你提供了评价回归模型适用性的必要工具，它能帮助发现并纠正问题。
8.3.1 标准方法
最常见的方法就是对lm()函数返回的对象使用plot()函数，可以生成评价模型拟合情况的四幅图形。
fit <- lm(weight ~ height, data=women)
par(mfrow=c(2,2))
plot(fit)
每张图的对应作用：为理解这些图形，我们来回顾一下OLS回归的统计假设。
1. 正态性：正态Q-Q图
当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正 态分布。“正态Q-Q图”（Normal Q-Q，右上）是在正态分布对应的值下，标准化残差的概 率图。若满足正态假设，那么图上的点应该落在呈45度角的直线上若不是如此，那么 就违反了正态性的假设。
2. 独立性
你无法从这些图中分辨出因变量值是否相互独立，只能从收集的数据中来验证。 上面的例子中，没有任何先验的理由去相信一位女性的体重会影响另外一位女性的体重。 假若你发现数据是从一个家庭抽样得来的，那么可能必须要调整模型独立性的假设。
3. 线性：残差图与拟合图
若因变量与自变量线性相关，那么残差值与预测（拟合）值就没有任何系统关联。换句话说，除了白噪声，模型应该包含数据中所有的系统方差。在“残差图与拟合图” （Residuals vs Fitted，左上）中可以清楚地看到一个曲线关系，这暗示着你可能需要对回 归模型加上一个二次项。模型拟合值与残差之间的散点图，红色线条表示二者关系的平滑曲线。若满足线性假设，二者应该不存在任何趋势性的关系，即红色线条应该与y = 0基本重合。
4. 同方差性：位置尺度图
若满足不变方差假设，那么在“位置尺度图”（Scale-Location Graph，左下） 中，水平线周围的点应该随机分布。该图似乎满足此假设。
离群点：残差与杠杆图（可读性较差且不实用
最后一幅“残差与杠杆图”（Residuals vs Leverage，右下）提供了你可能关注的单个观测点 的信息。从图形可以鉴别出离群点、高杠杆值点和强影响点。
一个观测点是离群点，表明拟合回归模型对其预测效果不佳（产生了巨大的或正或负的 残差）。一个观测点有很高的杠杆值，表明它是一个异常的预测变量值的组合。也就是说，在预 测变量空间中，它是一个离群点。因变量值不参与计算一个观测点的杠杆值。一个观测点是强影响点（influential observation），表明它对模型参数的估计产生的影响过 大，非常不成比例。强影响点可以通过Cook距离即Cook’s D统计量来鉴别。
再看看这个例子二次拟合的诊断图（按前述，二次拟合效果应该更好）
基本符合了线性假设、残差正态性（除了 观测点13）和同方差性（残差方差不变）。观测点15看起来像是强影响点（根据是它有较大的 Cook距离值），删除它将会影响参数的估计。事实上，删除观测点13和15，模型会拟合得会更 好。
剔除点后的模型（删除观测点13和15）
newfit <- lm(weight~ height + I(height^2), data=women[-c(13,15),])
但是对于删除数据，要非常小心，因为本应是你的模型去匹配数据， 而不是反过来
对多元回归也进行一下诊断：
states <- as.data.frame(state.x77[,c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
par(mfrow=c(2,2))
plot(fit)
8.3.2 改进的方法：car包（基于2.x版本）
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")
id.method = "identify"选项能够交互式绘图——
待图形绘制后，用鼠标单击图形内的点，将会标注函数中labels选项的设定值。敲击Esc键，从图 形下拉菜单中选择Stop，或者在图形上右击，都将关闭这种交互模式。此处，我已经鉴定出了Nevada 异常。当simulate=TRUE时，95%的置信区间将会用参数自助法（自助法可参见第12章）生成。除了Nevada，所有的点都离直线很近，并都落在置信区间内，这表明正态性假设符合得很好。但是你也必须关注Nevada，它有一个很大的正残差值（真实值－预测值），表明模型低估了该州 的谋杀率。
改进：
states["Nevada",] # murder=11.5
fitted(fit)["Nevada"] # 3.878958
residuals(fit)["Nevada"] # 7.621042
rstudent(fit)["Nevada"] # 3.542929
可以看到，Nevada的谋杀率是11.5%，而模型预测的谋杀率为3.9%。
可视化误差还有其他方法，比如使用代码清单8-6中的代码。
residplot()函数生成学生化 残差柱状图（即直方图），并添加正态曲线、核密度曲线和轴须图。它不需要加载car包。
绘制学生化残差图的函数
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)
正如你所看到的，除了一个很明显的离群点，误差很好地服从了正态分布
2. 误差的独立性
之前章节提过，判断因变量值（或残差）是否相互独立，最好的方法是依据收集数据方式的 先验知识。
1）时间序列数据
例如，时间序列数据通常呈现自相关性——相隔时间越近的观测相关性大于相隔越远 的观测。
car包提供了一个可做Durbin-Watson检验的函数，能够检测误差的序列相关性。在多元 回归中，使用下面的代码可以做Durbin-Watson检验：
durbinWatsonTest(fit)
p值不显著（p=0.282）说明无自相关性，误差项之间独立。滞后项（lag=1）表明数据集中 每个数据都是与其后一个数据进行比较的。该检验适用于时间独立的数据，对于非聚集型的数据 并不适用。注意，durbinWatsonTest()函数使用自助法（参见第12章）来导出p值。如果添加 了选项simulate=TRUE，则每次运行测试时获得的结果都将略有不同。
3. 线性：成分残差图（偏残差图）
看看因变量与自变量之间是否呈非线性关系，也可以看看是否有不同于已设定线性模型的系统偏差，图形可用car包中的crPlots()函数绘制。
library(car)
crPlots(fit)
结果如图8-11所示。若图形存在非线性，则说明你可能对预测变量的函数形式建模不够充分， 那么就需要添加一些曲线成分，比如多项式项，或对一个或多个变量进行变换（如用log(X)代 替X），或用其他回归变体形式而不是线性回归。
从图8-11中可以看出，成分残差图证实了你的线性假设，线性模型形式对该数据集看似是合 适的。
4. 同方差性
car包提供了两个有用的函数，可以判断误差方差是否恒定。
library(car)
ncvTest(fit)
spreadLevelPlot(fit)
ncvTest()函数生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。若检验显著， 则说明存在异方差性（误差方差不恒定）。
计分检验不显著（p=0.19），说明满足方差不变假设。
spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值 与拟合值的关系。
你还可以通过分布水平图 （图8-12）看到这一点，其中的点在水平的最佳拟合曲线周围呈水平随机分布。若违反了该假设， 你将会看到一个非水平的曲线。
代码结果建议幂次变换（suggested power transformation）的含义 是，经过p次幂（Yp ）变换，非恒定的误差方差将会平稳。此处p=1.2，是指异方差性很不明显，因此建议幂次接近1（不需要进行变换）。
8.3.3 线性模型假设的综合验证：习gvlma包中的gvlma()函数
gvlma()函数由Pena和Slate（2006）编 写，能对线性模型假设进行综合验证，同时还能做偏斜度、峰度和异方差性的评价。换句话说， 它给模型假设提供了一个单独的综合检验（通过/不通过）。
library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)
8.3.4 多重共线性
假设你正在进行一项握力研究，自变量包括DOB（Date Of Birth，出生日期）和年龄。你用 握力对DOB和年龄进行回归，F检验显著，p<0.001。
但是当你观察DOB和年龄的回归系数时，却 发现它们都不显著（也就是说无法证明它们与握力相关）。到底发生了什么呢？
原因是DOB与年龄在四舍五入后相关性极大。回归系数测量的是当其他预测变量不变时，某 个预测变量对响应变量的影响。
那么此处就相当于假定年龄不变，然后测量握力与年龄的关系， 这种问题就称作多重共线性（multicollinearity）。它会导致模型参数的置信区间过大，使单个系数 解释起来很困难。
多重共线性可用统计量VIF（Variance Inflation Factor，方差膨胀因子）进行检测。VIF的平 方根表示变量回归参数的置信区间能膨胀为与模型无关的预测变量的程度（因此而得名）。
car 包中的vif()函数提供VIF值。一般原则下， vif（开根号） >2就表明存在多重共线性问题。
library(car)
vif(fit)
sqrt(vif(fit)) > 2 # problem?
8.4 异常观测值
8.4.1 离群点
离群点是指那些模型预测效果不佳的观测点。它们通常有很大的、或正或负的残差（Yi–Ŷi）。 正的残差说明模型低估了响应值，负的残差则说明高估了响应值。
你已经学习过一种鉴别离群点的方法：图8-9的Q-Q图，落在置信区间带外的点即可被认为是 离群点。另外一个粗糙的判断准则：标准化残差值大于2或者小于–2的点可能是离群点，需要特 别关注
car包也提供了一种离群点的统计检验方法。outlierTest()函数可以求得最大标准化残差 绝对值Bonferroni调整后的p值：
library(car)
outlierTest(fit)
此处，你可以看到Nevada被判定为离群点（p=0.048）。注意，该函数只是根据单个最大（或 正或负）残差值的显著性来判断是否有离群点。若不显著，则说明数据集中没有离群点；若显著， 则你必须删除该离群点，然后再检验是否还有其他离群点存在。
8.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倍的位置。定位函数（locator function）能以交互模式绘图： 单击感兴趣的点，然后进行标注，停止交互时，用户可按Esc键退出，或从图形下拉菜单中选择 Stop，或直接右击图形。此图中，可以看到Alaska和California非常异常，查看它们的预测变量值，与其他48个州进行 比较发现：
Alaska收入比其他州高得多，而人口和温度却很低；California人口比其他州府多得多， 但收入和温度也很高。
高杠杆值点可能是强影响点，也可能不是，这要看它们是否是离群点。
8.4.3 强影响点
即对模型参数估计值影响有些比例失衡的点。例如，若移除模型的一个观测点时 模型会发生巨大的改变，那么你就需要检测一下数据中是否存在强影响点了。
有两种方法可以检测强影响点：Cook距离，或称D统计量，以及变量添加图（added variable plot）。
一般来说，Cook’s D值大于4/(n–k–1)，则表明它是强影响点，其中n为样本量大小，k是预 测变量数目。可通过如下代码绘制Cook’s D图形：
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是强影响点。若删除这些点，将会导致回归模型 截距项和斜率发生显著变化。注意，虽然该图对搜寻强影响点很有用，但我逐渐发现以1为分割 点比4/(n–k–1)更具一般性。若设定D=1为判别标准，则数据集中没有点看起来像是强影响点。Cook’s D图有助于鉴别强影响点，但是并不提供关于这些点如何影响模型的信息。
变量添加图弥补了这个缺陷。对于一个响应变量和k个预测变量，你可以如下图创建k个变量添加图。
所谓变量添加图，即对于每个预测变量Xk，绘制Xk在其他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")
图8-16反映出Nevada和Rhode Island是离群点，New York、California、Hawaii和Washington 有高杠杆值，Nevada、Alaska和Hawaii为强影响点。
8.5 改进措施
有四种方法可以处理违背回归假设的问题：
删除观测点；
变量变换；
添加或删除变量；
使用其他回归方法。
8.5.1 删除观测点（谨慎使用）
删除离群点通常可以提高数据集对于正态假设的拟合度，而强影响点会干扰结果，通常也会 被删除。删除最大的离群点或者强影响点后，模型需要重新拟合。若离群点或强影响点仍然存在， 重复以上过程直至获得比较满意的拟合。
8.5.2 变量变换
当模型不符合正态性、线性或者同方差性假设时，一个或多个变量的变换通常可以改善或调 整模型效果。
当模型违反正态假设时，通常可以对响应变量x尝试某种变换。
# 代码清单8-10 Box-Cox正态变换
library(car)
summary(powerTransform(states$Murder))
当违反了线性假设时，对预测变量y进行变换常常会比较有用。
car包中的boxTidwell()函 数通过获得预测变量幂数的最大似然估计来改善线性关系。
看MLE of lambda：0.87和1.36，还有p-value说明使用：
响应变量变换还能改善异方差性（误差方差非恒定）。在代码清单8-7中，你可以看到car包 中spreadLevelPlot()函数提供的幂次变换应用，
8.5.3 增删变量
有时，添加一个重要变量可以解决我们已经讨论过 的许多问题，删除一个冗余变量也能达到同样的效果。
删除变量在处理多重共线性时是一种非常重要的方法。如果你仅仅是做预测，那么多重共 线性并不构成问题，但是如果还要对每个预测变量进行解释，那么就必须解决这个问题。最常 见的方法就是删除某个存在多重共线性的变量（某个变量 vif（开根号） >2）。
另外一个可用的方法便是 岭回归——多元回归的变体，专门用来处理多重共线性问题。
8.5.4 尝试其他方法
情况方法多重共线性
岭回归删除变量存在离群点/强影响点稳健回归违背了正态性假设非参数模型存在显著非线性非线性回归违背误差独立性 用专门研究误差结构的模型 比如：时间序列模型/多层次回归模型 还可以用广泛应用的广义线性模型
8.6 选择“最佳”的回归模型
是不是所 有的变量都要包括？还是去掉那个对预测贡献不显著的变量？是否需要添加多项式项和/或交互 项来提高拟合度？最终回归模型的选择总是会涉及预测精度（模型尽可能地拟合数据）与模型简 洁度（一个简单且能复制的模型）的调和问题。如果有两个几乎相同预测精度的模型，你肯定喜 欢简单的那个。
8.6.1 模型比较（模型数量少时）：anova()（需要嵌套模型）和aic()（不需要）
anova()函数可以比较两个嵌套模型的拟合优度。
所谓嵌套模型，即它的一 些项完全包含在另一个模型中。在states的多元回归模型中，我们发现Income和Frost的回归系 数不显著，此时你可以检验不含这两个变量的模型与包含这两项的模型预测效果是否一样好
例子：
此时是模型1嵌套在模型2中。anova()函数同时还对是否应该添加Income和Frost到线性模型 中进行了检验。由于检验不显著（p=0.994），我们可以得出结论：不需要将这两个变量添加到线 性模型中，可以将它们从模型中删除。
AIC也可以用来比较模型：考虑了模型的 统计拟合度以及用来拟合的参数数目
（Akaike Information Criterion，赤池信息准则）
它考虑了模型的 统计拟合度以及用来拟合的参数数目。AIC值较小的模型要优先选择，它说明模型用较少的参数 获得了足够的拟合度。该准则可用AIC()函数实现：
此处AIC值表明没有Income和Frost的模型更佳（AIC值更小）。
注意，ANOVA需要嵌套模型，而AIC方法不 需要。
比较两模型相对来说更为直接，但如果有4个、10个或者100个可能的模型该怎么办呢？这便 是下节的主题。
8.6.2 变量选择：逐步回归法（stepwise method）和全子集回归（all-subsets regression）。
1. 逐步回归
逐步回归中，模型会一次添加或者删除一个变量，直到达到某个判停准则为止。
例如，向前 逐步回归（forward stepwise regression）每次添加一个预测变量到模型中，直到添加变量不会使 模型有所改进为止。向后逐步回归（backward stepwise regression）从模型包含所有预测变量开始， 一次删除一个变量直到会降低模型质量为止。而向前向后逐步回归（stepwise stepwise regression， 通常称作逐步回归，以避免听起来太冗长），结合了向前逐步回归和向后逐步回归的方法，变量 每次进入一个，但是每一步中，变量都会被重新评价，对模型没有贡献的变量将会被删除，预测 变量可能会被添加、删除好几次，直到获得最优模型为止。
逐步回归法的实现依据增删变量的准则不同而不同。
MASS包中的stepAIC()函数可以实现 逐步回归模型（向前、向后和向前向后），依据的是精确AIC准则。
例子：向后回归
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, direction="backward")
开始时模型包含4个（全部）预测变量，然后每一步中，AIC列提供了删除一个行中变量后 模型的AIC值，none中的AIC值表示没有变量被删除时模型的AIC。
第一步，Frost被删除，AIC 从97.75降低到95.75第二步，Income被删除，AIC继续下降，成为93.76。然后再删除变量将会 增加AIC，因此终止选择过程。
逐步回归法其实存在争议，虽然它可能会找到一个好的模型，但是不能保证模型就是最佳模 型，因为不是每一个可能的模型都被评价了。为克服这个限制，便有了全子集回归法。
2. 全子集回归
全子集回归是指所有可能的模型都会被检验。分析员可以选择展示所有可能的结 果，也可以展示n个不同子集大小（一个、两个或多个预测变量）的最佳模型。例如，若nbest=2， 先展示两个最佳的单预测变量模型，然后展示两个最佳的双预测变量模型，以此类推，直到包含 所有的预测变量。
全子集回归可用leaps包中的regsubsets()函数实现。你能通过R平方、调整R平方或 Mallows Cp统计量等准则来选择“最佳”模型：
R平方含义是预测变量解释响应变量的程度；调整R平方与之类似，但考虑了模型的参数数 目。R平方总会随着变量数目的增加而增加。当与样本量相比，预测变量数目很大时，容易导致 过拟合。R平方很可能会丢失数据的偶然变异信息，而调整R平方则提供了更为真实的R平方估计。另外，Mallows Cp统计量也用来作为逐步回归的判停规则。广泛研究表明，对于一个好的模型， 它的Cp统计量非常接近于模型的参数数目（包括截距项）。
在代码清单8-14中，我们对states数据进行了全子集回归。结果可用leaps包中的plot() 函数绘制（如图8-17所示），或者用car包中的subsets()函数绘制
第一行中（图底部开始），可以看到含intercept（截距项）和Income 的模型调整R平方为0.33含intercept和Population的模型调整R平方为0.1。跳至第12行，你会看到 含intercept、Population、Illiteracy和Income的模型调整R平方值为0.54，而仅含intercept、Population 和Illiteracy的模型调整R平方为0.55。此处，你会发现含预测变量越少的模型调整R平方越大（对 于非调整的R平方，这是不可能的）。图形表明，双预测变量模型（Population和Illiteracy）是最佳 模型（最高处）
大部分情况中，全子集回归要优于逐步回归，因为考虑了更多模型。但是，当有大量预测变 量时，全子集回归会很慢。一般来说，变量自动选择应该被看做是对模型选择的一种辅助方法， 而不是直接方法。拟合效果佳而没有意义的模型对你毫无帮助，主题背景知识的理解才能最终指 引你获得理想的模型。
8.7 深层次分析：介绍评价模型泛化能力和变量相对重要性的方法
8.7.1 交叉验证（好像就是数挖课要的作业？？？？）
上一节我们学习了为回归方程选择变量的方法。若你最初的目标只是描述性分析，那么只需 要做回归模型的选择和解释。但当目标是预测时，你肯定会问：“这个方程在真实世界中表现如 何呢？”
从定义来看，回归方法本就是用来从一堆数据中获取最优模型参数。
对于OLS回归，通过使 得预测误差（残差）平方和最小和对响应变量的解释度（R平方）最大，可获得模型参数。由于 等式只是最优化已给出的数据，所以在新数据集上表现并不一定好。
在本章开始，我们讨论了一个例子，生理学家想通过个体锻炼的时间和强度、年龄、性别与 BMI来预测消耗的卡路里数。如果用OLS回归方程来拟合该数据，那么仅仅是对一个特殊的观测 集最大化R平方，但是研究员想用该等式预测一般个体消耗的卡路里数，而不是原始数据。你知 道该等式对于新观测样本表现并不一定好，但是预测的损失会是多少呢？你可能并不知道。通过 交叉验证法，我们便可以评价回归方程的泛化能力。
所谓交叉验证，即将一定比例的数据挑选出来作为训练样本，另外的样本作保留样本，先在 训练样本上获取回归方程，然后在保留样本上做预测。由于保留样本不涉及模型参数的选择，该 样本可获得比新数据更为精确的估计
在k重交叉验证中，样本被分为k个子样本，轮流将k–1个子样本组合作为训练集，另外1个子 样本作为保留集。
这样会获得k个预测方程，记录k个保留样本的预测表现结果，然后求其平均值。 （当n是观测总数目，且k为n时，该方法又称作刀切法，jackknifing。）
bootstrap 包中的 crossval() 函数可以实现 k 重交叉验证。
# R平方的k重交叉验证：
# 定义了shrinkage()函数，创建了一个包含预测变量和预测值的矩阵，可获得初始R平方以及交叉验证的R平方。
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")
}
对states数据所有预测变量进行回归，然后再用shrinkage()函数做10重交叉验证：
# 10重交叉验证
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.448）。（注意，观测被随机分配到k个群组中，因此每次运行shrinkage() 函数，得到的结果都会有少许不同。
通过选择有更好泛化能力的模型，还可以用交叉验证来挑选变量。例如，含两个预测变量 （Population和Illiteracy）的模型，比全变量模型R平方减少得更少（0.03 vs. 0.12）
这使得双预测变量模型显得更有吸引力。
其他情况类似，基于大训练样本的回归模型和更接近于感兴趣分布的回归模型，其交叉验证 效果更好。R平方减少得越少，预测则越精确。
8.7.2 相对重要性
不相关：根据预测变量与响应变量的相关系数来进 行排序。
最简单的莫过于比较标准化的回归系数，它表 示当其他预测变量不变时，该预测变量一个标准差的变化可引起的响应变量的预期变化（以标准 差单位度量）。
在进行回归分析前，可用scale()函数将数据标准化为均值为0、标准差为1的数 据集，这样用R回归即可获得标准化的回归系数。
（注意，scale()函数返回的是一个矩阵，而 lm()函数要求一个数据框，你需要用一个中间步骤来转换一下。data.frame()）
例子：
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是最不重要的
还有许多其他方法可定量分析相对重要性。比如，可以将相对重要性看作每个预测变量（本 身或与其他预测变量组合）对R平方的贡献。Ulrike Grömping写的relaimpo包涵盖了一些相对重 要性的评价方法（http://prof.beuth-hochschule.de/groemping/relaimpo/）。
相对权重（relative weight）是一种比较有前景的新方法，它是对所有可能子模型添加一个预 测变量引起的R平方平均增加量的一个近似值（Johnson，2004；Johnson & Lebreton，2004；LeBreton & Tonidandel，2008）。
relweights（）函数
# relweights()函数，计算预测变量的相对权重
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)
}
现将代码清单8-17中的relweights()函数应用到states数据集。
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")
报错了，不知道为啥
根据相对权重法，Illiteracy有最大的相对重要性，余 下依次是Frost、Population和Income。它比标准化回归系数更为直观
8.8 小结
文章链接
发表评论