可以发现残差图和Adjusted R2的提升都甚微，因此没有必要作非线性转换。
4、模型分析
（1）模型比较
前面只是简单得用Adjusted R2来比较模型，《R语言实战》里介绍了可以用方差分析来比较嵌套模型（即它的一些项完全包含在另一个模型中）有没有显著性差异。方差分析的思想是：如果线性模型y~x1+x2+x3与y~x1+x2没有显著性差异，若同时x3变量对模型也不显著，那就没必要加上变量x3。下面进行试验：
aovfit1 <- lm(Murder~Population+Illiteracy+Income+Frost,data=states)
aovfit2 <- lm(Murder~Population+Illiteracy,data=states)
anova(aovfit1,aovfit2)
Analysis of Variance Table
Model 1: Murder ~ Population + Illiteracy + Income + Frost
Model 2: Murder ~ Population + Illiteracy
Res.Df RSS Df Sum of Sq F Pr(>F)
1 45 289.17
2 47 289.25 -2 -0.078505 0.0061 0.9939
summary(aovfit1)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.235e+00 3.866e+00 0.319 0.7510
Population 2.237e-04 9.052e-05 2.471 0.0173 *
Illiteracy 4.143e+00 8.744e-01 4.738 2.19e-05 ***
Income 6.442e-05 6.837e-04 0.094 0.9253
Frost 5.813e-04 1.005e-02 0.058 0.9541
Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared: 0.567, Adjusted R-squared: 0.5285
F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08
summary(aovfit2)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.652e+00 8.101e-01 2.039 0.04713 *
Population 2.242e-04 7.984e-05 2.808 0.00724 **
Illiteracy 4.081e+00 5.848e-01 6.978 8.83e-09 ***
Residual standard error: 2.481 on 47 degrees of freedom
Multiple R-squared: 0.5668, Adjusted R-squared: 0.5484
F-statistic: 30.75 on 2 and 47 DF, p-value: 2.893e-09
Income和Frost两个变量不显著，两个模型之间没有显著性差异，就可以不加这两个变量。删去这两个不显著的变量后，R2略微减少，Adjusted R2增大，这也符合二者的定义。
《R语言实战》里还介绍到了用AIC(Akaike Information Criterion，赤池信息准则)值来比较模型，AIC值越小的模型优先选择，原理不明。
aovfit1 <- lm(Murder~Population+Illiteracy+Income+Frost,data=states)
aovfit2 <- lm(Murder~Population+Illiteracy,data=states)
AIC(aovfit1,aovfit2)
df AIC
aovfit1 6 241.6429
aovfit2 4 237.6565
第二个模型AIC值更小，因此选第二个模型（真是简单粗暴）。注：ANOVA需限定嵌套模型，AIC则不需要。可见AIC是更简单也更实用的模型比较方法。
（2）变量选择
这里的变量选择与最开始的变量选择同也不同，虽然是一回事，但一开始是一个粗略的变量的选择，主要是为了构建模型；这里则要进行细致的变量选择来调整模型。
前面提到的向前或向后选择或者是同时向前向后选择变量都是逐步回归法。MASS包中的stepAIC() 函数可以实现逐步回归模型（向前、向后和向前向后），依据的是精确AIC准则。以下实例是向后回归法：
library(MASS)
aovfit1 <- lm(Murder~Population+Illiteracy+Income+Frost,data=states)
stepAIC(aovfit1,direction = "backward") # “forward”为向前选择，"backward"为向后选择，"both"为混合选择
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
可见原本的4元回归模型向后退了两次，最终稳定成了2元回归模型，与前面模型比较的结果一致。
《R语言实战》里提到了逐步回归法的局限：不是每个模型都评价了，不能保证选择的是“最佳”模型。比如上例中，从Murder ~ Population + Illiteracy + Income + Frost到Murder ~ Population + Illiteracy + Income再到Murder~Population+Illiteracy虽然AIC值确实在减少，但Murder ~ Population + Illiteracy + Frost没评价，如果遇到变量很多的情况下，逐步回归只沿一个方向回归，就有可能错过最优的回归方向。
library(leaps)
leaps <- regsubsets(Murder~Population+Illiteracy+Income+Frost,data=states,nbest=4)
plot(leaps,scale = "adjr2")
横坐标是变量，纵坐标是Adjusted R2，可见除截距项以外，只选定Population和Illiteracy这两个变量，可以使线性模型有最大的Adjusted R2。
全子集回归比逐步回归范围更广，模型优化效果更好，但是一旦变量数多了之后，全子集回归迭代的次数就很多，就会很慢。
事实上，变量的选择不是机械式地只看那几个统计指标，更主要的是根据数据的实际意义，从业务角度上来选择合适的变量。
线性模型变量的选择在《统计学习》后面的第6章还会继续讲到，到时继续综合讨论。
（3）交互项
交互项《统计学习》中花了一定篇幅来描写，但在《R语言实战》是在方差分析章节中讨论。添加变量间的交互项有时可以改善线性关系，提高Adjusted R2。针对数据的实际意义，如果两个基本上是独立的，也很难产生交互、产生协同效应的变量，那就不必考虑交互项；只有从业务角度分析，有可能产生协同效应的变量间才考虑交互项。
涉及到交互项有一个原则：如果交互项是显著的，那么即使变量不显著，也要放在回归模型中；若变量和交互项都不显著，则可以都不放。
（4）交叉验证
Andrew Ng的Machine Learning中就提到了，模型对旧数据拟合得好不一定就对新数据预测得好。因此一个数据集应当被分两训练集和测试集两部分（或者训练集、交叉验证集、测试集三部分），训练好的模型还要在新数据中测试性能。
所谓交叉验证，即将一定比例的数据挑选出来作为训练样本，另外的样本作保留样本，先在训练样本上获取回归方程，然后在保留样本上做预测。由于保留样本不涉及模型参数的选择，该样本可获得比新数据更为精确的估计。
在k 重交叉验证中，样本被分为k个子样本，轮流将k?1个子样本组合作为训练集，另外1个子样本作为保留集。这样会获得k 个预测方程，记录k 个保留样本的预测表现结果，然后求其平均值。
bootstrap包中的crossval()函数可以实现k重交叉验证。
shrinkage <- function(fit, k = 10) {
require(bootstrap)
# define functions
theta.fit <- function(x, y) {
lsfit(x, y)
}
theta.predict <- function(fit, x) {
cbind(1, x) %*% fit$coef
}
# matrix of predictors
x <- fit$model[, 2:ncol(fit$model)]
# vector of predicted values
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")
}
这个自定义的shrinkage()函数用来做k重交叉验证，比计算训练集和交叉验证集的R方差异。这个函数里涉及到一个概念：复相关系数。复相关系数实际上就是y和fitted(y)的简单相关系数。对于一元线性回归，R2就是简单相关系数的平方；对于多元线性回归，R2是复相关系数的平方。这个我没有成功地从公式上推导证明成立，就记下吧。这个方法用到了自助法的思想，这个在统计学习后面会细致讲到。
fit <- lm(Murder ~ Population + Income + Illiteracy +
Frost, data = states)
shrinkage(fit)
Original R-square = 0.5669502
10 Fold Cross-Validated R-square = 0.441954
Change = 0.1249963
可见这个4元回归模型在交叉验证集中的R2下降了0.12之多。若换成前面分析的2元回归模型——
fit2 <- lm(Murder ~ Population + Illiteracy , data = states)
shrinkage(fit2)
Original R-square = 0.5668327
10 Fold Cross-Validated R-square = 0.517304
Change = 0.04952868
这次R2下降只有约0.05。R2减少得越少，则预测得越准确。
5、模型应用
（1）预测
最重要的应用毫无疑问就是用建立的模型进行预测了。构建好模型后，可用predict()函数进行预测——
fit2 <- lm(Murder ~ Population + Illiteracy , data = states)
predict(fit2,
newdata = data.frame(Population=c(2000,3000),Illiteracy=c(1.7,2.2)),
interval = "confidence")
fit lwr upr
1 9.037174 8.004911 10.06944
2 11.301729 9.866851 12.73661
这里newdata提供了两个全新的点供模型来预测。还可以用interval指定返回置信区间（confidence）或者预测区间（prediction），这也反映了统计与机器学习的一个差异——可解释性。注意置信区间考虑的是平均值，而预测区间考虑的是单个观测值，所以预测区间永远比置信区间广，因此预测区间考虑了单个观测值的不可约误差；而平均值同时也把不可约误差给抵消掉了。
（2）相对重要性
有的时候需要解释模型中各个自变量对因变量的重要程度，简单处理可以直接看系数即可，《R语言实战》里自定义了一个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))
# correlations between original predictors and new orthogonal variables
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda^2
# regression coefficients of Y on orthogonal variables
beta <- solve(lambda) %*% rxy
rsquare <- colSums(beta^2)
rawwgt <- lambdasq %*% beta^2
import <- (rawwgt/rsquare) * 100
lbls <- names(fit$model[2:nvar])
rownames(import) <- lbls
colnames(import) <- "Weights"
# plot results
barplot(t(import), names.arg = lbls, ylab = "% of R-Square",
xlab = "Predictor Variables", main = "Relative Importance of Predictor Variables",
sub = paste("R-Square = ", round(rsquare, digits = 3)),
...)
return(import)
}
不要在意算法原理和代码逻辑这种细节，直接看结果：
fit <- lm(Murder ~ Population + Illiteracy + Income +
Frost, data = states)
relweights(fit, col = "lightgrey")
Weights
Population 14.723401
Illiteracy 59.000195
Income 5.488962
Frost 20.787442
