#####第八章：回归#####
#回归分析都是统计学的核心。
#它其实是一个广义的概念，通指那些用一个或多个自变量（也称预测变量或解释变量）来预测因变量(也称响应变量、效标变量或结果变量的方法。
#通常，回归分析可以用来挑选与因变量相关的自变量，可以描述两者的关系，也可以生成一个方程，通过自变量来预测因变量。
#由于回归分析在现代统计学中非常重要，本章将对其进行一些深度讲解。
#首先，我们将看一看如何拟合和解释回归模型，其次，回顾一系列识别模型潜在问题的方法，并学习如何解决它们。
#然后，我们将探究变量选择问题。对于所有可用的自变量，如何确定哪些变量包含在最终的模型中？
#再次，我们将讨论一般性问题。模型在现实世界中的表现到底如何？
#最后，我们看看相对重要性问题。模型所有的自变量中，哪个最重要，哪个第二重要，哪个最无关紧要？
####8.1 回归的多面性####
#表8-1记录了回归分析的各种变体
#简单线性：用一个连续型自变量预测一个连续型因变量
#多项式：用一个连续型自变量预测一个连续型因变量，模型的关系是n阶多项式
#多元线性：用两个或多个自变量预测一个连续型因变量
#多层：用拥有层级结构的数据预测一个因变量（例如学校中教室里的学生）。也被称为分层模型、嵌套模型或混合模型
#多变量：用一个或多个自变量预测多个因变量
#Logistic：用一个或多个自变量预测一个分类因变量
#泊松：用一个或多个自变量预测一个代表频数的因变量
#Cox比例风险：用一个或多个自变量预测一个事件（死亡、失败或旧病复发）发生的时间
#时间序列：对误差项相关的时间序列数据建模
#非线性：用一个或多个自变量预测一个连续型因变量，不过模型是非线性的
#非参数：用一个或多个自变量预测一个连续型因变量，模型的形式源自数据形式，不事先设定
#稳健：用一个或多个自变量预测一个连续型因变量，能抵御强影响点的干扰
###本章中，我们的重点是普通最小二乘（OLS）回归法，包括简单线性回归、多项式回归和多元线性回归。
###第13章将介绍其他回归模型，包括Logistic回归和泊松回归。
#8.1.1 OLS回归的适用场景
#OLS回归是通过自变量的加权和来预测连续型因变量，其中权重是通过数据估计而得的参数。
####8.2 OLS回归####
#以下是chatgpt对OLS回顾的定义，有兴趣可以浅读以下（John注):
##OLS回归，即最小二乘法（Ordinary Least Squares）回归，是一种广泛用于线性回归分析的统计方法。OLS回归的目标是通过最小化预测值和实际值之间的残差平方和来确定模型参数，即截距项和斜率。在OLS回归中，通常假定数据满足某些统计假设，包括正态性、独立性、线性关系和同方差性。当数据违反这些假设时，统计显著性检验的结果和置信区间的准确性可能会受到影响。
#根据预测变量的数量，OLS回归可以分为几种类型。当回归模型只包含一个因变量和一个自变量时，称为简单线性回归。当只有一个预测变量，但同时包含变量的幂（如X²，X³）时，称为多项式回归。当有多个预测变量时，则称为多元线性回归。
#总的来说，OLS回归是一种强大的统计工具，可以用于预测和解释因变量与自变量之间的关系。然而，它也有一些限制和假设，需要在应用时予以注意。
##为了恰当解释OLS回归模型，数据必须满足：
##正态性、独立性、线性和同方差性（方差齐性）
#####8.2.1 用函数lm()拟合回归模型
#拟合线性模型最基本的函数是lm()，其格式为：
myfit <- lm(formula, data)
#formula是要拟合的模型公式，
#data是一个数据框，包含用于拟合的数据。
#拟合公式形式如下
Y ~ X1 + X2 + ... +Xk
#~左边是因变量，右边为各个自变量
#自变量用+分割
#表8-2中记录常用的符号：
# ~ 分隔符号，左边为因变量，右边为自变量。例如，要通过x、z和w预测y，代码为y ~ x + z + W
# + 分隔自变量
# : 表示自变量的交互项。例如，要通过x、z及x与z的交互项预测y，代码为y~x+z+x:z
# * 表示所有可能的交互项的简洁方式。代码y~x*z*w可展开为y~x+z+w+x:z+x:w+z:w+x:z:w
# ^ 表示交互项达到某个次数.代码y ~（x+z+w）^2可展开为y~x+z+w+x:z+x:w+z:w
# . 表示数据框内除因变量外的所有变量的占位符。例如，若一个数据框包含变量x、y、z和w.代码y~.可展开为y~x+z+w
# - 减号，表示从方程中移除某个变量。例如，y~(x+z+w)^2-x:w可展开为y~x+z+w+x:z+z:w
# -1 删除截距项。例如，表达式y~x-1拟合y在x上的回归，并强制直线通过原点
# I() 从算术的角度来解释括号中的元素。例如，y~x+（z+w)^2将展开为y~x+z+w+ z:w。而代码y~x+I((z + w)^2)将展开为y~x+h，h是一个由z和w的平方和创建的新变量
# function 可以在表达式中用的数学函数。例如，log(y)~x+z+w表示通过x、z和w来预测log(y)
###除了1m(），表8-3还列出了其他一些对做简单或多元回归分析有用的函数。
#拟合模型后，将这些函数应用于1m（）返回的对象，可以得到更多额外的模型信息。
#summary():展示拟合模型的详细结果
#coefficients():列出拟合模型的模型参数（截距项和斜率）
#confint():提供模型参数的置信区间（默认95%）
#fitted():列出拟合模型的预测值
#residuals():列出拟合模型的残差值
#anova():生成一个拟合模型的方差分析表，或者比较两个或更多拟合模型的方差分析表
#vcov():列出模型参数的协方差矩阵
#AIC():输出赤池信息量准则
#plot():生成评价拟合模型的诊断图
#predict():用拟合模型对新的数据集预测因变量值
#####8.2.2 简单线性回归
#让我们通过一个回归示例来熟悉表8-3中的函数。
#基础安装中的数据集women提供了15个年龄在30~39岁间女性的身高和体重信息。
#我们想通过身高来预测体重，通过获得一个方程，我们可以分辨出那些过重或过轻的个体。
#代码清单8-1提供了分析过程，图8-1展示了结果图形。
fit <- lm(weight ~ height, data = women) #
summary(fit) #结果中标注Estimate即为截距和斜率
women$weight #列出实际的体重值
fitted(fit) #列出预测的体重值
residuals(fit) #列出残差值（实际-预测）
plot(women$height, women$weight,
xlab = "Height(in inches)",
ylab = "Weight(in pounds)") #绘制体重~身高的点图）
abline(fit) #添加拟合曲线
#按照summary中的结果，可以得出预测方程为Weight=-87.52+3.45×Height
#方程拟合解释率99.1%
#####8.2.3 多项式回归
#图8-1可以看出，可以添加一个二次项（即X²）来提高回归的预测精度。
#二次项拟合方程如下：
fit2 <- lm(weight ~ height + I(height ^ 2), data = women)
#I(height^2)表示向预测方程添加一个身高的平方项。
#I函数将括号的内容看作R的一个常规表达式。
#因为^（参见表8-2）符号在表达式中有特殊的含义，会调用并不需要的东西，所以此处必须要用这个函数。
#代码清单8-2如下：
fit2 <- lm(weight ~ height + I(height ^ 2), data = women)
summary(fit2)
plot(women$height, women$weight,
xlab = "Height(in inches)",
ylab = "Weight(in pounds)")
lines(women$height, fitted(fit2))
#根据summary函数中Estimate结果，预测方程为Weight=261.87-7.3×Heigth+0.083×Height²
#方程拟合解释率99.9%
#多项式回归仍是线性的，即使使用更为复杂的运算，预测方程仍是线性的
#详见书第174页
#n次多项式生成一个带有n-1个弯曲的曲线，
#拟合三次多项式可用
fit3 <- lm(weight ~ height + I(height ^ 2) + I(height ^ 3), data = women)
#虽然可用，但更高的项显得没有必要
#####8.2.4 多元线性回归
#当自变量不止一个时，简单线性回归就变成了多元线性回归，分析也稍微复杂些。
#以R自带的state.x77数据集为例，我们想探究一个州的犯罪率和其他因素的关系，包括人口、文盲率、收入和结霜天数（温度在冰点以下的平均天数）。
#lm()函数需要使用数据框，但state.x77是矩阵，所以需要进行转换。
states <- as.data.frame(state.x77[, c("Murder", "Population", "Illiteracy", "Income", "Frost")])
#多元回归分析中，第一步是检查自变量间的相关性。
#cor()函数提供了二变量之间的相关系数
#car包中的scatterplotMartix()函数则会生成散点图矩阵
install.packages("car")
states <- as.data.frame(state.x77[, c("Murder", "Population", "Illiteracy", "Income", "Frost")])
cor(states)
library(car)
scatterplotMatrix(states, smooth = F, main = "Scatter Plot Matrix")
#代码清单8-4
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 带交互项的多元线性方程
#代码清单8-5:
fit <- lm(mpg ~ hp + wt + hp:wt, data = mtcars)
summary(fit)
#根据结果可以得出预测方程为
#mpg = 49.81 - 0.12hp - 8.2wt + 0.027×hp×wt
#通过effects包中的effect()函数可以用图形展示交互项的结果
install.packages("effects")
#其格式为
plot(effect(term, mod, , xlevels), multiline = TRUE)
#其中term即模型要绘制的项，mod为通过lm()拟合的模型，xlevels是一个列表，指定变量要设定的常量值
#multiline=True选项表示添加相应的直线，
#lines选项指定每条线的类型，1为实线，2为虚线；3为点线
library(effects)
plot(effect("hp:wt", fit,,list(wt = c(2.2, 3.2, 4.2))),
lines = c(1,2,3), multiline = TRUE)
####8.3 回归模型诊断####
#在8.2节中，我们用lm（）函数来拟合OLS回归模型，通过summary（)函数获取模型参数和相关统计量。
#但是，没有任何输出告诉我们模型是否合适，我们对模型参数推断的信心依赖于它在多大程度上满足OLS回归模型统计假设。
#数据的无规律性或者错误设定了自变量与因变量的关系，这将致使我们的模型产生巨大的偏差。
#一方面，我们可能得出某个自变量与因变量无关的结论，但事实上它们是相关的；
#另一方面，情况可能恰好相反。当我们的模型应用到真实世界中时，预测效果可能很差，误差显著。
#我们通过confint（）函数的输出来看看8.2.4节中states多元回归的问题。
states <- as.data.frame(state.x77[, c("Murder", "Population", "Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
confint(fit)
#其结果显示了95%的置信区间
#####8.3.1 标准方法
#R基础安装中提供了大量检验回归分析中统计假设的方法。
#最常见的方法就是对1m()函数返回的对象使用p1ot（）函数，生成评估模型拟合情况的4幅图形。
#下面是简单线性回归的例子：
fit <- lm(weight ~ height, data = women)
par(mfrow = c(2,2)) #将图调整为2×2
plot(fit) #直接输出回归方程的图出现的是正态性、线性、独立性和同方差性检验的四张图
par(mfrow = c(1,1))
#生成的图形见图8-5。par(mfrow=c(2，2))将plot()函数绘制的4幅图形组合在一个大的2×2的图中。
#第2个par()函数又恢复为显示单个图形。
#正态性：Q-Q图，符合正态性应该是斜向上45°的斜线
#独立性：从四幅图无法判断是否独立，只能从数据中得到
#线性：看残差-拟合值图，如果是线性的则为随机分布，如果为曲线则需要添加更多的项，远离0点则说明有异常值，扇形或不均匀分散则代表方差不齐
#同方差性：如果方差齐，位置-尺位图则为随机分布
#最后一个残差-杠杆图可以观测到离群值、高杠杆点和强影响点。
#以上诊断图在R中有更好的方法，后续会讲到。
#####8.3.2 改进的方法
#car包中提供了大量的函数，增强拟合和评估回归模型的能力。
#表8-4
#qqPlot()：分位数比较图
#durbinWatsonTest()： 对误差自相关性做Durbin-Watson检验
#crPlots()： 成分与残差图
#ncvTest()： 对非恒定的误差方差做得分检验
#spreadLevelPlot()： 分散水平图
#outlierTest()： Bonferroni 离群点检验
#avPlots()： 添加的变量图形
#influencePlot()： 回归影响图
#vif()： 方差膨胀因子
##1、正态性
#与基础包中的 plot（)函数相比，gqplot（）函数提供了更为精确的正态假设检验方法，
#它绘制出了在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 = list(method = "identilfy"),
simulate = TRUE, main = "Q-Q Plot")
#qqplot()函数生成的概率图见图8-6。
#id=list(method="identify")选项能够交绘图——待图形绘制后，用鼠标单击图形内的点，将会标注函数中labels选项的设定值。
#“Esc”键，或者单击图形右上角的Finish按钮，都将关闭这种交互模式。(我无法实现交互，不知道为什么，John注)
#此处，我已经识别出了Nevada异常点。
#当simulate=TRUE 时，95%的置信区间将会用参数自助法（自助法请参见第12章）生成。
#这里还必须注意以下Nevada异常点
states["Nevada",] #显示Nevada异常点的值
fitted(fit)["Nevada"] #显示预测值
residuals(fit)["Nevada"] #显示残差值
rstudent(fit)["Nevada"] #计算学生化值
##2、误差的独立性
#判断因变量值（或残差）是否相互独立，最好的方法是依据收集数据方式的先验知识。
#car包提供了一个可做Durbin-Watson检验的函数，能够检测误差的序列相关性。
#在多元回归中，使用下面的代码可以做Durbin-Watson检验：
library(car)
durbinWatsonTest(fit)
#结果中p值大于0.05说明没有相关性，误差项之间独立。
#滞后项（lag=1）表明数据集中每个数据都是与其后一个数据进行比较的。
#该检验适用于时间独立的数据，对于非聚集型的数据并不适用。
#注意，durbinWatsonTest（）函数使用自助法（参见第12章）来导出p值。
#如果添加了选项simulate=TRUE，则每次运行测试时获得的结果都将略有不同。
##线性
#通过成分残差图(component plus residual plot)也称偏残差图(partial residual plot）
#可以看看因变量与自变量之间是否呈非线性关系，
#也可以看看是否有不同于已设定线性模型的系统偏差，图形可用car包中的crPlots（）函数绘制。
library(car)
crPlots(fit)
#若图形存在非线性，则说明我们可能对自变量的函数形式建模不够充分那么就需要添加一些曲线成分，
#比如多项式项，或对一个或多个变量进行变换（如用log(x)代替x），或用其他回归变体形式而不是线性回归。
#从图8-7中可以看出，成分残差图证实了我们的线性假设，线性模型形式对该数据集看似是合适的。
##4、同方差性
#car包还提供了两个有用的函数，可以判断误差方差是否恒定不变。
#ncvTest（）函数生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。若检验显著，则说明存在异方差性（误差方差不恒定）。
#spreadLevelPlot（）函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系。
#函数应用如代码清单8-6所示。
library(car)
ncvTest(fit) #p>0.05说明同方差假设
spreadLevelPlot(fit) #点在水平最佳拟合曲线周围随机分布代表同方差性
#spreadLevelPlot(fit)结果建议Suggested power transformation: 1.209626
#说明要将Y值进行幂次转换，即Y^1.209626更符合曲线，如果图形为非水平，例如结果为0.5，则需要进行Y%0.5转换
#如果结果接近1，则不需要转换
#####8.3.3 多重共线性
#多重共线性是指线性回归模型中的解释变量（或自变量）之间存在线性相关关系。
#当这种相关关系存在时，模型的估计可能会失真或难以准确估计。
#这通常是由于经济数据的限制或模型设计不当导致的。
#完全共线性的情况并不多见，一般出现的是在一定程度上的共线性，即近似共线性。
#为了检测和处理多重共线性，可以使用一些统计方法，如计算变量间的相关系数、方差膨胀因子（Variance Inflation Factor, VIF）等。
#一般情况下vif＞10就表明存在多重共线性问题
#代码清单8-7如下：
library(car)
vif(fit)
vif(fit) > 10 # problem?
####8.4 异常观测值####
#一个全面的回归分析要覆盖对异常观测值的分析，包括离群点、高杠杆值点和强影响点。
#这些数据点需要更深入的研究，因为它们在一定程度上与其他观测点不同，可能对结果产生较大的负面影响。
###8.4.1离群点
#离群点是指那些模型预测效果不佳的观测点。
#它们通常有很大的、或正或负的残差。
#正的残差说明模型对因变量的预测值过小，负的残差则说明模型对因变量的预测值过大。
#落在置信区间带外的点即可被认为是离群点。
#另外一个粗糙的判断准则：标准化残差值大于2或者小于-2的点可能是离群点，需要我们特别关注。
#car包提供了一种离群点的统计方法：
#outlierTest（）函数可以求得标准化残差绝对值Bonferroni调整后的p值
library(car)
outlierTest(fit)
#，该函数只是根据单个最大(或正或负)残差值的显著性来判断是否有离群点。
#若不显著，则说明数据集中没有离群点；若显著，则必须删除该离群点，然后再检验是否还有其他离群点存在。
###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倍的位置。
#定位函数（locatorfunction）能以交互模式绘图：单击感兴趣的点，然后进行标注，停止交互时，用户可按“Ese”键退出，或单击图形右上角的Finish按钮。
#高杠杆值点可能是强影响点，也可能不是，这要看它们是否是离群点。
###8.4.3强影响点
#强影响点，即对模型参数估计值的影响有些比例失衡的点。
#例如，若移除模型的一个观测点时模型会发生巨大的改变，那么我们就需要检测一下数据中是否存在强影响点了。
#有两种方法可以检测强影响点：Cook距离，或称D统计量，以及变量添加图(added variable plot)。
#一般来说，Cook'sD值大于 4/(n-k-1)，则表明它是强影响点，其中n为样本量，k是自变量数目。
#可通过如下代码绘制Cook'sD图（图8-10)：
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'sD图有助于识别强影响点，但是并不提供关于这些点如何影响模型的信息。
#变量添加图弥补了这个缺陷。对于一个因变量和k个自变量，我们可以如下图创建k个变量添加图。
#所谓变量添加图，即对于每个自变量Xk，绘制Xk，在其他k-1个自变量上回归的残差值相对于因变量在其他k-1个自变量上回归的残差值的关系图。
#car包中的avPlots（）函数可提供变量添加图：
library(car)
avPlot(fit, ask = FALSE, id = list(method = "idnetify"))
#结果如图8-11所示。图形一次生成一个，用户可以通过单击点来判断强影响点。
#按下“Esc'键，或单击图形右上角的 Finish 按钮，便可移动到下一个图形。我已在左下图中识别出 Alaska为强影响点。
#图中的直线表示相应自变量的实际回归系数。
#我们可以想象删除某些强影响点后直线的变，以此来估计它的影响效果。
#当然，利用car包中的influencePlot（）函数，我们还可以将离群点、杠杆值和强影响点的信息整合到一幅图形中：
library(car)
influencePlot(fit, id = "noteworthy → TRUE", main = "Influence Plot",
sub = "Circle size is proportional to Cook`s distance")
#Nevada和Rhode Island是离群点，Califonia和Hawaii有高杠杆值，Nevada和Alaska为强影响点。
#如果将id=TRUE替换为id=1ist(method="identify")，则我们通过单击鼠标来交互式识别观测点（按“Esc”键或者单击Finish按钮退出）。
#图8-12 影响图。
#纵坐标超过+2 或小于-2 的州可被认为是离群点，水平轴超过 0.2 或0.3的州有高杠杆值（自变量值的异常组合）。
#圆圈大小与影响成比例，圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影响点
####8.5 改进措施####
###8.5.1删除观测点
#删除离群点通常可以提高数据集对于正态假设的拟合度，而强影响点会干扰结果，通常也会被删除。
#删除最大的离群点或者强影响点后，模型需要重新拟合。
#若离群点或强影响点仍然存在，重复以上过程直至获得比较满意的拟合。
#对删除观测点需持谨慎态度。
###8.5.2变量变换
#当模型不符合正态性、线性或者同方差性假设时，一个或多个变量的变换通常可以改善或调整模型效果。
#若了出现较大的偏斜，那么进行log变换是很有用的。
#常见的变量变换见书第191页。
#当模型违反正态假设时，通常可以对因变量尝试某种变换。car包中的powerTransform()函数通过λ的最大似然估计来正态化变量Xλ。
#由此产生的转换称为 Box-Cox变换。在代码清单8-8中，此变换用于数据集states。
library(car)
summary(powerTransform(states$Murder))
#结果表明可以用Murder的0.6次方来正态化转换Murder
#由于0.6很接近0.5，我们可以尝试用平方根变换来提高模型正态性的符合程度。
#当违反了线性假设时，对自变量进行变换常常会比较有用。
#car包中的boxTidwell()函数通过获得自变量幂数的最大似然估计来改善线性关系。
#下面的例子用州的人口和文盲率来预测谋杀率，对模型进行了Box-Tidwell变换：
library(car)
boxTidwell(Murder ~ Population + Illiteracy, data = states)
#结果提示Population的0.86939次方和Illiteracy的1.35812次方能大大改善线性关系。
#但依据Pr(>|t|)值又提示不需要进行变换
#最后，因变量变换还能改善异方差性（误差方差非恒定）。
#在代码清单8-8中，我们可以看到car包中spreadLevel Plot()函数提供的幂次变换提高了同方差性，
#不过，states例子满足了同方差性假设，不需要进行变量变换。
###8.5.3增删变量
#改变模型的变量将会影响模型的拟合度。
#有时，添加一个重要变量可以解决我们已经讨论过的许多问题，删除一个冗余变量也能达到同样的效果。
#删除变量在处理多重共线性时是一种非常重要的方法。
#如果我们仅仅是做预测，那么多重共线性并不构成问题，但是如果还要对每个自变量进行解释，那么就必须解决这个问题。
#最常见的方法就是删除某个存在多重共线性的变量（某个变量vif>10）。
#另外一个可用的方法便是岭回归——多元回归的变体——专门用来处理多重共线性问题。
###8.5.4尝试其他方法
#正如刚才提到的，处理多重共线性的一种方法是拟合一种不同类型的模型（本例中是岭回归）。
#其实，如果存在离群点和/或强影响点，可以使用稳健回归模型替代OLS回归。
#如果违背了正态性假设，可以使用非参数回归模型。
#如果存在显著的非线性，可以尝试非线性回归模型。
#如果违背了误差独立性假设，还能用那些专门研究误差结构的模型，比如时间序列模型或者多层次回归模型。
#最后，我们还能转向广泛应用的广义线性模型，它能适用于许多OLS 回归假设不成立的情况。
#在第13章中，我们将会介绍其中一些方法。
#至于什么时候需要提高OLS回归拟合度，什么时候需要换一种方法，这些判断是很复杂的，需要依靠你对学科知识的理解，判断出哪个模型能够提供最佳结果。
####8.6 选择“最佳”回归模型####
###8.6.1模型比较
#用基础安装中的anova（）函数可以比较两个嵌套模型的拟合优度。
#所谓嵌套模型，即它的回归方程的项完全包含在另一个模型中。
#在states的多元回归模型中，我们发现Income和Frost 的回归系数不显著，此时我们可以检验不含这两个变量的模型与包含这两项的模型的预测效果是否一样好（见代码清单8-9）。
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)
#此处，模型1嵌套在模型2中。
#anova（）函数同时还对是否应该在Population和Illiteracy之外还要添加Income和Frost到线性模型中进行了检验。
#由于检验不显著（p=0.994），我们可以得出结论：不需要将这两个变量添加到线性模型中，可以将它们从模型中删除。
#AIC(Akaike Information Criterion， 赤池信息量准则）也可以用来比较模型，它考虑了模型的统计拟合度以及用来拟合的参数数目。
#AIC值较小的模型要优先选择，它说明模型用较少的参数获得了足够的拟合度。
#该准则可用AIC（）函数实现（见代码清单8-10）。
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
fit2 <- lm(Murder ~ Population +Illiteracy, data = states)
AIC(fit2, fit1)
#此处AIC值表明没有Income和Frost的模型更佳
###8.6.2变量选择
#从大量候选变量中选择最终的自变量有以下两种流行的方法：
#逐步回归法（stepwise method）和全子集回归（all-subsets regression）。
#1.逐步回归
#逐步回归中，模型会一次添加或者删除一个变量，直到达到某个判停准则为止。
#例如，向前逐步回归（forward stepwise regression)每次添加一个自变量到模型中，直到添加变量不会使模型有所改进为止。
#向后逐步回归（backward stepwise regression）从模型包含所有自变量开始，一次删除一个变量，直到会降低模型质量为止。
#而向前向后逐步回归（stepwise stepwise regression,通常称作逐步回归，以避免听起来太冗长），结合了向前逐步回归和向后逐步回归的方法：变量每次进入一个，但是每一步中，变量都会被重新评估，对模型没有贡献的变量将会被删除；自变量可能会被添加、删除好几次，直到获得最优模型为止。
#逐步回归法的实现依据增删变量的准则不同而不同。
#R基础包中的step（)函数可以实现逐步回归模型（向前、向后和向前向后），依据的是AIC准则。
#代码清单8-11中，我们用向后逐步回归来处理多元回归问题。
states <- as.data.frame(state.x77[, c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
step(fit, direction = "backward")
#开始时模型包含4个（全部）自变量，然后每一步中，AIC列提供了删除一个行中变量后模型的AIC值，
#<none>中的AIC值表示没有变量被删除时模型的AIC。
#第1步，Frost被删除。AIC从97.75降低到95.75；
#第2步，Income 被删除，AIC继续下降，成为93.76。
#然后再删除变量将会增加AIC，因此终止选择过程。
#逐步回归法其实存在争议，虽然它可能会找到一个好的模型，但是不能保证模型就是最佳模型，
#因为不是每一个可能的模型都被评估了。
#2.全子集回归
#顾名思义，全子集回归是指所有可能的模型都会被检验。
#分析员可以选择展示所有可能的结果，也可以展示n个不同子集大小（一个、两个或多个自变量）的最佳模型。
#例如，若nbest=2，那么先展示两个最佳的单自变量模型，然后展示两个最佳的双自变量模型，以此类推，直到包含所有的自变量。
#全子集回归可用leaps包中的regsubsets（）函数实现。我们可以通过R平方、调整R平方或Mallows Cp统计量等准则来选择最佳模型。
#正如我们所看到的，R平方的含义是自变量解释因变量的程度；
#调整R平方与之类似，但考虑了模型的参数数目。
#R平方总会随着自变量数目的增加而增加。当与样本量相比，自变量数目很大时，容易导致过拟合。
#R平方很可能会丢失数据的偶然变异信息，而调整R平方则提供了更为真实的总体R平方估计。
#在代码清单8-12中，我们对states数据集进行了全子集回归。
#leaps包在一单个图形中展示结果，但是我发现许多人对此图感到困惑。
#以下代码用表格的形式展示了部分结果，我相信这些结果更易于理解。
install.packages("leaps")
#安装的时候找不到镜像，用下面的代码直接指定镜像原后再次安装（John注）
#options("repos" = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
library(leaps)
states <- as.data.frame(state.x77[, c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
leaps <- regsubsets(Murder ~ Population +Illiteracy + Income +
Frost, data = states, nbest = 4)
subsTable <- function(obj, scale) {
x <- summary(leaps)
m <- cbind(round(x[[scale]],3), x$which[, -1])
colnames(m)[1] <- scale
m[order(m[, 1]),]}
subsTable(leaps, scale = "adjr2")
#表中的每一行代表一个模型。
#第1列是模型中自变量的数目，
#第2列是描述每个模型拟合度的度量值（本例中为调整R平方），每一行按该度量值进行排序。（注意：可以用其他比例替代adir2，有关其他比例的选项请参阅?regsubsets）。
#行中的1和0表示模型添加或删除了哪些变量。
#例如，只有一个自变量Income的模型的调整R平方为0.033，具有自变量 Population、Illiteracy和Income的模型的调整 R平方为 0.539。
#相比之下，只有自变量 Pooulation和Illiteracy 的模型的调整 R平方为 0.548。
#从本例中可以看到自变量更少的模型实际拥有更大的调整R平方（对非调整R平方不适用）。此表格表明具有两个自变量的模型是最佳模型。
#大部分情况中，全子集回归要优于逐步回归，因为它考虑了更多模型。
#但是，当有大量自变量时，全子集回归会很慢。
#一般来说，变量自动选择应该被看作是模型选择的一种辅助方法，而不是直接方法。
#拟合效果佳而没有意义的模型对我们毫无帮助，对学科背景知识的理解才能最终指引我们获得理想的模型。
####8.7 深层次分析####
#下面介绍评估模型泛化能力和自变量相对重要性的方法。
###8.7.1交叉验证
#从定义来看，回归方法本就是用来从一堆数据中获取最优模型参数的。
#对于 OLS 回归，通过使得预测误差（残差）平方和最小和对因变量的解释度（R平方）最大，可获得模型参数。
#因为方程只是最优化已给出的数据，所以在新数据集上表现并不一定好。
#在本章讨论了一个例子，运动生理学家想通过个体锻炼的时长和强度、年龄、性别与BMI来预测消耗的卡路里数。
#如果用 OLS 回归方程来拟合该数据，那么我们获得的仅仅是对一个特定的观测值集合最大化 R平方的模型参数。
#但是，研究员想用该方程预测一般个体消耗的卡路里数，而不是原始研究中的卡路里数。
#我们知道该方程对于新观测值表现并不一定好，但是预测的损失会是多少呢？我们可能并不知道。
#通过交叉验证法，我们便可以评估回归方程的泛化能力。
#所谓交叉验证，是指将一定比例的数据挑选出来作为训练样本，另外的样本作保留样本，先在训练样本上获取回归方程，然后在保留样本上做预测。
#由于保留样本不涉及模型参数的选择。该样本可获得比新数据更为精确的估计。
#在k重交叉验证中，样本被分为k个子样本，轮流将k-1个子样本组合作为训练集，另外1个子样本作为保留集。
#这样会获得k个预测方程，记录k个保留样本的预测表现结果，然后求其平均值。
#（当n是观测值总数目，且k等于n时，该方法又称作刀切法，jackknifing。）
#bootstrap包中的crossval(）函数可以实现k重交叉验证。
#在代码清单8-13中，shrinkage()函数对模型的R平方统计量做了k重交叉验证。
install.packages("bootstrap")
library(bootstrap)
shrinkage <- function(fit, k=10, seed = 1){
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]
set.seed(seed)
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("Oringinal R-square =", r2, "\n")
cat(k, "Fold Cross-Validated R-square =", r2cv, "\n")
}
#以上函数定义了shrinkage（）函数,创建了一个包含自变量和预测值的矩阵，
#可获得初始R平方和残差的标准误，以及交叉验证的R平方和残差标准误
#以下对states数据集中的所有变量进行回归，在用shrinkage函数做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.356）。
#（注意，观测值被随机分配到个群组中，因此用随机数种子可让结果可重现。）
#通过选择有更好泛化能力的模型，还可以用交叉验证来挑选变量。
#例如，含两个自变量(Population和Illiteracy)的模型，比全变量模型R平方减少得更少：
###8.7.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))
zfit <- lm(Murder ~ Population + Income +Illiteracy + Frost, data = zstates)
coef(zfit)
#此处可以看到，当人口、收入和结霜天数不变时，文盲率每增加一个标准差将会使谋杀率增加0.68个标准差。
#根据标准回归系数，我们可认为Illiteracy是最重要的自变量，而Frost是最不重要的。
#还有许多其他方法可定量分析相对重要性。
#比如，可以将相对重要性看作每个自变量（本身或与其他自变量组合）对R平方的贡献。
#相对权重（relative weight）是一种比较有前景的新方法，它是对所有可能子模型添加一个自变量引起的R平方平均增加量的一个近似值。
#代码清单8-14提供了一个生成相对权重的函数。
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-Sguare=", round(rsquare, digits = 3)),
...)
return(import)
}
###代码清单8-15中，将函数relweiahts()应用到states 数据集，根据人口、收人和结霜天数预测谋杀率。
states <- as.data.frame(state.x77[, c("Murder", "Population", "Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Income +Illiteracy + Frost, data = states)
relweights(fit, col = "blue")
#以上这段代码我无法运行，不知道为什么（John注）
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
