R语言实战（四）回归

本文对应《R语言实战》第8章：回归

回归是一个广义的概念，通指那些用一个或多个预测变量（也称自变量或解释变量）来预测响应变量（也称因变量、效标变量或结果变量）的方法。通常，回归分析可以用来挑选与相应变量相关的解释变量，可以描述两者的关系，也可以生成一个等式，通过解释变量来预测响应变量。

回归分析的各种变体

回归类型	用途
简单线性	用一个量化的解释变量预测一个量化的响应变量
多项式	用一个量化的解释变量预测一个量化的响应变量，模型的关系是n阶多项式
多元线性	用两个或多个量化的解释变量预测一个量化的响应变量
多变量	用一个或多个解释变量预测多个响应变量
Logistic	用一个或多个解释变量预测一个类别型响应变量
泊松	用一个或多个解释变量预测一个代表频数的响应变量
Cox比例风险	用一个或多个解释变量预测一个事件（死亡、失败或旧病复发）发生的时间
时间序列	对误差项相关的时间序列数据建模
非线性	用一个或多个量化的解释变量预测一个量化的响应变量，不过模型是非线性的
非参数	用一个或多个量化的解释变量预测一个量化的响应变量，模型的形式源自数据形式，不事先设定
稳健	用一个或多个量化的解释变量预测一个量化的响应变量，能抵御强影响点的干扰

本章主要介绍OLS（最小二乘回归法）

统计假设：

正态性：对于固定自变量值，因变量值成正态分布

独立性：因变量值之间相互独立

线性：因变量与自变量之间线性相关

同方差性：因变量的方差不随自变量的水平不同而变化

基本格式：

myfit <- lm(formula, data)

　　

表达式中常用符号 

~	分隔符号，左边为响应变量，右边为解释变量
+	分隔预测变量
:	表示预测变量的交互项 e.g. x:z
*	表示所有可能交互项的简洁方式 e.g. x * z = x + z + x:z
^	表示交互项达到某个次数
.	表示除因变量外的所有变量
-	表示从formula中移除某变量
-1	删除截距项（强制回归直线通过原点）
I()	从算术的角度来解释括号中的元素（即在I中的表达式符号为算术符号）
function	可以在表达式中用的数学函数

拟合模型的其他常用函数

summary()	展示拟合模型的详细结果
coefficients()	列出拟合模型的模型参数
confint()	提供模型参数的置信区间（默认95%）
fitted()	列出拟合模型的预测值
residuals()	列出拟合模型的残差值
anova()	生成一个拟合模型的方差分析表，或者比较两个或更多拟合模型的方差分析表
vcov()	列出模型参数的协方差矩阵
AIC()	输出赤池信息统计量
plot()	生成评价拟合模型的诊断图
predict()	用拟合模型对新的数据集预测响应变量值

代码示例

#简单线性回归
fit <- lm(weight ~ height, data = women)
summary(fit)

#多项式回归
fit2 <- lm(weight ~ height + I(height ^ 2), data = women)
summary(fit2)
fit3 <- lm(weight ~ height + I(height ^ 2) + I(height ^ 3), data = women)
summary(fit3)

#多元线性回归
fit4 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
summary(fit4)

#有交互项的多元线性回归
fit5 <- lm(mpg ~ hp + wt + hp:wt)
summary(fit5)

　　

回归诊断：

标准方法：plot(fit)

四张图：Residuals vs Fitted(残差图与拟合图), Normal Q-Q(正态Q-Q图), Scale-Location(位置尺度图), Residuals vs Leverage(残差与杠杆图)

诊断统计假设方式：

正态性：Normal Q-Q图上的点在45度直线上

独立性：无法从图中读出，但可以收集的数据中验证

线性：Residuals vs Fitted随机分布（若有明显趋势，则暗示回归模型有相应的趋势分布）

同方差性：Scale-Location水平线周围的点随机分布

Residuals vs Leverage图可以鉴别离群点、高杠杆值点和强影响点（由于该图可读性差，通过其他方式完成这个工作）

改进的方法：car包

qqplot()	分位数比较图
durbinWatsonTest()	对误差自相关性做Durbin-Watson检验
crPlots()	成分与残差图
ncvTest()	对非恒定的误差方差做得分检验
spreadLevelPlot()	分散水平检验
outlierTest()	Bonferroni离群点检验
avPlots()	添加的变量图形
inluencePlot()	回归影响图
scatterplot()	增强的散点图
scatterplotMatrix()	增强的散点图矩阵
vif()	方差膨胀因子

正态性：

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

#id.method = “identify”选项能够交互式绘图，鼠标单击图形内的点，将会标注函数中label选项的设定值
#simulate = TRUE时，95%的置信区间将会用参数自助法
library(car)
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
qqPlot(fit, labels = row.names(states), id.method = “identify”,
simulate = TRUE, main = “Q-Q Plot”)

　　

误差的独立性： 

最好的方法是依据收集数据方式的先验知识。car包中Durbin-Watson检验函数也可以检测误差的序列相关性（p值不显著说明无自相关性）。

线性：

成分残差图（component plus residual plot，或称偏残差图，partial residual plot）

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

同方差性：

ncvTest()函数生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。（若检验显著，则说明存在异方差性）

spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系。另外会输出一个建议幂次变换（suggested power transformation），含义是：经过p次幂变换，非恒定的方差将会平稳。

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

使用gvlma包中的函数gvlma()

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

　　

输出项可以看到满足OLS回归模型的所有统计假设，若p<0.05可认为违反假设条件 

多重共线性：

两个预测变量之间存在高度的相关性，造成模型参数的置信区间过大，使单个系数解释起来很困难。（简单来说，回归系数测量的是当其他预测变量不变时，某个预测变量与响应变量之间的关系，而有两个预测变量之间相关时，就变成了假定某变量不变，测量响应变量与该变量之间的关系）

多重共线性可用统计量VIF（variance inflation factor, 方差膨胀因子）进行检测。VIF的平方根表示变量回归参数的置信区间能膨胀为与模型无关的预测变量的程度。car包中的vif()函数可以计算。一般原则下，就表明存在多重共线性问题。

library(car)
vif(fit)
sqrt(vif(fit)) > 2
#返回各变量及布尔值，TRUE代表存在问题，FALSE代表不存在问题

异常观测值：离群点、高杠杆值、强影响点

离群点：指那些模型预测效果不佳的观测点（很大的残差）。鉴别：Q-Q图中落在置信区间外的点，或标准化残差值大于2或小于-2的点。

car包中的函数outlierTest()可以检测到单个离群点（p值显著），要检验其他离群点，需要先将第一步的离群点删除重新检测。

高杠杆值点：与其他预测变量有关的离群点，即许多异常的预测变量值组合起来的，与响应变量值没有关系。

鉴别：通过帽子统计量（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)    

　　

强影响点：即对模型参数估计值影响有些比例失衡的点。 

检测：Cook距离（或称D统计量）或者变量添加图（added variable plot）。一般来说，Cook’s D值大于4/(n-k-1),则表明它是强影响点，其中n为样本量大小，k是预测变量数目。（注：如果以D=1为判别标准可能更具有一般性）

#绘制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”)

　　

Cook’s D图有助于鉴别强影响点，但是并不能提供关于这些点如何影响模型的信息。变量添加图弥补了这个缺陷。 

library(car)
avPlots(fit, ask = FALSE, onepage = TRUE, id.method = “identify”)

　　

将离群点、高杠杆值和强影响点信息整合到一幅图形中： 

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的有高杠杆值，圆圈大小与影响成比例，圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影响点。

改进措施：针对违背回归假设的问题

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

删除观测点：删除离群点或强影响点。谨慎。若是数据记录错误、没有遵守规程等的离群点可以删除，否则发掘为何该观测点不同于其他点也是研究的方向之一。

变量变换：模型不符合正态性、线性或者同方差性假设时，一个或多个变量变换可以改善效果。模型违反正态假设时，通常对响应变量尝试变换；违反线性假设时，对预测变量变换有效，同时也可以改善异方差性。需要说明的是，变换后的变量需要有意义，否则很难解释。

增删变量：删除变量可以处理多重共线性。但是如果仅是预测并不需要处理，处理多重共线性最重要的是对预测变量进行解释。

选择“最佳”的回归模型

#检验嵌套模型的拟合优度
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
fit2 <- lm(Murder ~ Population + Illiteracy, data = states)
anova(fit2, fit1)

　　

 如果检验不显著（p > 0.05），则使用较少变量的线性模型。

使用AIC（Akaike Information Criterion，赤池信息准则）比较模型：AIC值越小的模型越优先选择。

变量选择：逐步回归法（stepwise method）；全子集回归（all-subsets regression）

逐步回归法：向前逐步回归（forward stepwise）向后逐步回归（backward stepwise）向前向后逐步回归（stepwise stepwise）

library(MASS)
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
stepAIC(fit1, direction = “backward”)

　　

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

全子集回归法：所有可能的模型都会被检验，可用leaps包中的regsubsets()函数实现。一般来说，变量的自动选择应该被看作是对模型选择的一种辅助方法，而不是直接方法。拟合效果佳而没有意义的模型对你毫无帮助，主题背景知识的理解才能最终指引你获得理想的模型。

交叉验证：将一定比例的数据挑选出来作为训练样本，另外的样本作保留样本，先在训练样本上获得回归方程，然后在保留样本上做预测。

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

相对重要性：比较标准化的回归系数，表示当其他预测变量不变时，该预测变量一个标准差的变化可引起的响应变量的预期变化（以标准差单位度量）。在进入回归分析之前，可用scale()函数将数据标准化为均值为0、标准差为1的数据集，这样用R回归即可获得标准化的回归系数。（注意scale()函数返回的是一个矩阵，而lm()函数要求一个数据框）

zstates <- as.data.frame(scale(states))
zfit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data = zstates)
coef(zfit)

　　

另外可以根据相对权重判断相对重要性。比较各个回归系数，最大的就是最重要的。