回归分析概述

回归分析的应用

我们主要的困难有三个：发现有趣的问题， 设计一个有用的、可以测量的响应变量，以及收集合适的数据。

回归分析的分类

回归分析有多重分类，用的比较多的是logistic回归和Lasso回归。

普通最小二乘（ OLS）回归法

• 普通最小二乘（ OLS）回归法，包括简单线性回归、多项式回归和多元线性回归。 OLS回归是现今最常见的统计分析方法。OLS回归是通过预测变量的加权和来预测量化的因变量，其中权重是通过数据估计而得的参数。
• 利用OLS法通过一系列的预测变量来预测响应变量（OLS 回归说是在预测变量上“回归”响应变量——其名也因此而来）。
• OLS法的公式及含义：

OLS法的公式及含义

• OLS回归主要是通过减少响应变量的真实值与预测值的差值来获得模型参数（截距项和斜率），即使得残差平方和最小。

回归的R实现

简单线性回归

• 用一个量化的预测变量来预测一个量化的响应变量。
• 常用线性回归模型进行拟合。

多元线性回归

• 当预测变量不止一个时，简单线性回归就变成了多元线性回归。
• 多元回归分析中，第一步最好检查一下变量间的相关性。
  cor()函数提供了二变量之间的相关系数， car包中scatterplotMatrix()函数则会生成散点图矩阵。两者均是描述了两两预测变量之间的相关关系。
• lm(y~x1+x2+x3函数)可进行多元线性回归分析。此种多元线性回归并没有考虑交互项。

> states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])   ### 建立待研究变量的子数据框
> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)  ##进行多元线性拟合
> summary(fit)  ##输出拟合的详细结果

• 含交互项的多元线性回归。若两个预测变量的交互项显著，说明响应变量与其中一个预测变量的关系依赖于另外一个预测变量的水平。

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

通过effects包中的effect()函数，来用图形展示交互项的结果。例如如车重均值3.2各加减1，即取wt2.2,3.2,4.2来观察马力（hp）对车mpg的影响（mtcars数据）。

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)

• 在模型拟合好后，进行推断之前，必须对方法中暗含的统计假设进行检验。

回归诊断

• 使用lm()函数来拟合OLS回归模型，通过summary()函数获取模型参数和相关统计量。但是，没有任何输出告诉你模型是否合适，你对模型参数推断的信心依赖于它在多大程度上满足OLS模型统计假设。
• 在根据模型进行统计推断要建立在你的数据满足统计假设的前提之上。回归诊断技术向你提供了评价回归模型适用性的必要工具，它能帮助发现并纠正问题。
• 首先，对拟合后的模型进行OLS回归统计假设的验证。

标准方法

fit <- lm(weight ~ height, data=women)  ##对响应变量及解释变量进行拟合
par(mfrow=c(2,2))
plot(fit)  ##plot()函数可画出四副图形

OLS回归统计假设判定.png

“正态Q-Q图”（ Normal Q-Q，右上）是在正态分布对应的值下，标准化残差的概率图。当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正态分布。若满足正态假设，那么图上的点应该落在呈45度角的直线上；否则就违反了正态性的假设。
“残差图与拟合图”（ Residuals vs Fitted，左上）：若因变量与自变量线性相关，那么残差值与预测（拟合）值就没有任何系统关联。在“残差图与拟合图”中可以清楚地看到一个曲线关系，这暗示着你可能需要对回归模型加上一个二次项。反映线性假设。
“位置尺度图”（ Scale-Location Graph，左下）：若满足不变方差假设，那么在“位置尺度图”（ Scale-Location Graph，左下）中，水平线周围的点应该随机分布。该图似乎满足此假设。反映同方差性假设。
“残差与杠杆图”（ Residuals vs Leverage，右下）：提供了你可能关注的单个观测点的信息。从图形可以鉴别出离群点、高杠杆值点和强影响点。

残差与杠杆图的解释.png

• 根据拟合后的情况必要时进行二次拟合，分析。直至符合OLS回归假设。必要时可删除数据（谨慎！！！）来满足模型拟合更好。

改进的方法

正态性检验

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

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")

误差的独立性

• 判断因变量值（或残差）是否相互独立，最好的方法是依据收集数据方式的先验知识。
• car包提供了一个可做Durbin-Watson检验的函数，能够检测误差的序列相关性。该检验适用于时间独立的数据，对于非聚集型的数据并不适用。

> durbinWatsonTest(fit)
lag Autocorrelation D-W Statistic p-value
1     -0.201          2.32        0.282     ##p值不显著（ p=0.282）说明无自相关性，误差项之间独立
Alternative hypothesis: rho != 0  

线性

• 通过成分残差图（ component plus residual plot）也称偏残差图（ partial residual plot），你可以看看因变量与自变量之间是否呈非线性关系。
• 使用car包中的crPlots()函数绘制图像。(若图形存在非线性，则说明你可能对预测变量的函数形式建模不够充分，那么就需要添加一些曲线成分，比如多项式项，或对一个或多个变量进行变换（如用log(X)代替X），或用其他回归变体形式而不是线性回归。)

> library(car)
> crPlots(fit)

成分残差图.png

同方差性

• ncvTest()函数生成一个记分检验，若检验显著，则说明存在异方差性。
• spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图（水平分布图），展示标准化残差绝对值与拟合值的关系。

> library(car)
> ncvTest(fit)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare=1.7 Df=1 p=0.19   ##p无显著性说明同方差性
> spreadLevelPlot(fit)  ##用来绘制尺度位置图（分布水平图）
Suggested power transformation: 1.2

分布水平图.png

线性模型假设的综合验证

gvlma包中的gvlma()函数能对线性模型假设进行综合验证，同时还能做偏斜度、峰度和异方差性的评价。它给模型假设提供了一个单独的综合检验（通过/不通过）综合验证。

> library(gvlma)
> gvmodel <- gvlma(fit)
> summary(gvmodel)
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance= 0.05
Call:
gvlma(x=fit)
           Value   p-value    Decision
Global Stat 2.773 0.597 Assumptions acceptable.  #p>0.05线性模型假设全部接受
Skewness 1.537 0.215 Assumptions acceptable.
Kurtosis 0.638 0.425 Assumptions acceptable.
Link Function 0.115 0.734 Assumptions acceptable.
Heteroscedasticity 0.482 0.487 Assumptions acceptable.  #p>0.05同方差性假设接受

多重共线性

• 多重共线性（ multicollinearity）对于解释多元回归的结果非常重要。它会导致模型参数的置信区间过大，使单个系数解释起来很困难。
• 多重共线性可用统计量VIF（ Variance Inflation Factor， 方差膨胀因子）进行检测。
• car包中的vif()函数提供VIF值。一般原则下， vif 平方根>2就表明存在多重共线性问题。
• 代码示例：

>library(car)
>vif(fit)   #计算VIF值
# Population Illiteracy Income Frost
      #1.2      2.2        1.3  2.1
>sqrt(vif(fit)) > 2 # 判定是否大于2，即是否存在方差共线性
# Population Illiteracy Income Frost
# FALSE        FALSE     FALSE FALSE

• 多重共线性一般处理是删除变量。多元回归的变体岭回归专门用来处理多重共线性问题。

异常观测点的判定

• 利用car包中的influencePlot()函数绘制影响图，将离群点、杠杆值和强影响点的信息整合到一幅图形中：

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

影响图.png

违反OLS假设的处理

• 删除观测
• 变量转换
  car包中的powerTransform()函数通过λ的最大似然估计来正态化变量Xλ；

> library(car)
> summary(powerTransform(states$Murder))
bcPower Transformation to Normality
           Est.Power Std.Err. Wald Lower Bound Wald Upper Bound
states$Murder 0.6      0.26            0.088            1.1 ##Murder的0.6次方来改善正态化
Likelihood ratio tests about transformation parameters
                   LRT  df pval
LR test, lambda=(0) 5.7 1 0.017
LR test, lambda=(1) 2.1 1 0.145   ##Murder的lambda=1次方不能拒绝，因此Murder可不用变量转换

car包中的boxTidwell()函数通过获得预测变量幂数的最大似然估计来改善线性关系。

> library(car)
> boxTidwell(Murder~Population+Illiteracy,data=states)
          Score Statistic p-value MLE of lambda
Population    -0.32         0.75       0.87
Illiteracy     0.62         0.54       1.36   ##两者p>0.05均拒绝lamda转换，可不用变量转换

• 删除变量
• 改用其他方法：
  如果存在离群点和/或强影响点，可以使用稳健回归模型替代OLS回归。如果违背了正态性假设，可以使用非参数回归模型。如果存在显著的非线性，能尝试非线性回归模型等等。

选择“最佳”的回归模型

模型比较（两个模型）

• 用基础安装中的anova()函数可以比较两个嵌套模型的拟合优度。所谓嵌套模型，即它的一些项完全包含在另一个模型中。

> states <- as.data.frame(state.x77[,c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
> fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost,
data=states)  ##模型1，包括 Income + Frost两个自变量
> fit2 <- lm(Murder ~ Population + Illiteracy, data=states)
> anova(fit2, fit1) ##模型2
Analysis of Variance Table
Model 1: Murder ~ Population + Illiteracy
Model 2: Murder ~ Population + Illiteracy + Income + Frost
Res.Df RSS Df  Sum of Sq F Pr(>F)
1 47 289.246
2 45 289.167 2 0.079 0.0061 0.994  ##p=0.994>0.05提示检验不显著，可以得出结论：不需要将这两个变量添加到线性模型中，可以将它们从模型中删除

• AIC（ Akaike Information Criterion，赤池信息准则）也可以用来比较模型，它考虑了模型的统计拟合度以及用来拟合的参数数目。 AIC值较小的模型要优先选择，它说明模型用较少的参数获得了足够的拟合度。此准则可用AIC()函数实现：

fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost,
data=states)
fit2 <- lm(Murder ~ Population + Illiteracy, data=states)
AIC(fit1,fit2)
    df   AIC
fit1 6 241.6429
fit2 4 237.6565  ##模型2的AIC值较小，考虑应用，即剔除后两个变量

模型变量的选择

包括：逐步回归法（ stepwise method）和全子集回归（ all-subsets regression）

• 逐步回归法（ stepwise method）
  又包括向前逐步回归（ 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")   ##向后逐步回归法

• 全子集回归法（ all-subsets regression）
  全子集回归是指所有可能的模型都会被检验。全子集回归可用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, nbest=4)
plot(leaps, scale="adjr2")

深层次分析

交叉验证

对于OLS回归，通过使得预测误差（残差）平方和最小和对响应变量的解释度（ R平方）最大，可获得模型参数。由于等式只是最优化已给出的数据，所以在新数据集（验证数据集）上表现并不一定好。
所谓交叉验证，即将一定比例的数据挑选出来作为训练样本，另外的样本作保留样本，先在训练样本上获取回归方程，然后在保留样本上做预测。由于保留样本不涉及模型参数的选择，该样本可获得比新数据更为精确的估计。

• bootstrap 包 中 的 crossval() 函 数 可 以 实 现 k 重 交 叉 验 证
• 构建shrinkage()函数对模型的R平方统计量做了k重交叉验证

shrinkage <- function(fit, k=10){  ##shrinkage()函数，验证次数默认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) ##crossval() 函 数 可 以 实 现 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")
}

10次交叉验证结果.png

图中显示基于初始样本的R平方（ 0.567）过于乐观了。对新数据更好的方差解释率估计是交叉验证后的R平方（ 0.448）。

变量的相对重要性

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

预测变量的标准化回归系数.png

• 从图中看出，illiteracy标准化的回归系数最大，我们可认为Illiteracy是最重要的预测变量，而Frost是最不重要的。

注：本文仅供个人学习之用，使用时请注明文章引用自R in Action: Data Analysis and Graphics with R, Second Edition by Robert I. Kabacoff。