vlambda
2020-05-07

[R语言实战]第八章回归

第八章回归

回归分析是统计学的核心，通常指哪些用一个或多个预测变量来预测响应变量的方法

本章论述形式：

如何拟合和解释回归模型
回顾一系列鉴别模型潜在问题的方法，并学习解决它们
讨论变量选择问题
讨论一般性问题（模型在现实世界中表现如何）
预测变量重要性问题

8.1 回归的多面性

回归类型	用途
简单线性	用一个量化的解释变量预测一个量化的响应变量
多项式	用一个量化的解释变量预测一个量化的响应变量，模型的关系是n 阶多项式
多层	用拥有等级结构的数据预测一个响应变量（例如学校中教室里的学生）。也被称为`分层模型`、`嵌套模型`或`混合模型`
多元线性	用两个或多个量化的解释变量预测一个量化的响应变量
多变量	用一个或多个解释变量预测多个响应变量
Logistic	用一个或多个解释变量预测一个类别型响应变量
泊松	用一个或多个解释变量预测一个代表频数的响应变量
Cox 比例风险	用一个或多个解释变量预测一个事件（死亡、失败或旧病复发）发生的时间
时间序列	对误差项相关的时间序列数据建模
非线性	用一个或多个量化的解释变量预测一个量化的响应变量，不过模型是非线性的
非参数	用一个或多个量化的解释变量预测一个量化的响应变量，模型的形式源自数据形式，不事先设定
稳健	用一个或多个量化的解释变量预测一个量化的响应变量，能抵御强影响点的干扰

在本章内容中，重点是普通追小二乘回归法（OLS）包括简单线性回归、多项式回归和多元线性回归。

8.1.1 OLS回归的适用场景

OLS回归是通过预测变量的加权和来预测量化的因变量。主要困难有三个：发现问题，设计一个有用的可以测量的响应变量，以及收集合适的数据。

8.1.2 基础回顾

R函数如何拟合OLS回归模型、评价拟合优度，检验假设条件以及选择模型

##　８.２　ＯＬＳ回归

image-20200329144849858

OLS回归拟合模型的形式

其中，n伪观测的数目，k为预测变量的数目

等式中的变量	意义
Ŷi	第i次观测对应的因变量的预测值（具体来讲，它是在已知预测变量值的条件下，对Y分布估计的均值）
Xji	第i次观测对应的第j个预测变量值
β0	第i次观测对应的第j个预测变量值
	预测变量j的回归系数（斜率表示Xj改变一个单位所引起的Y的改变量）

我们的目标是通过减少响应变量的真实值与预测值的差值来获得模型参数（截距项和斜率）

为了能够恰当的解释OLS模型的参数，数据西要满足一下统计假设

正态性	对于固定的自变量值，因变量值成正态分布。
独立性	Yi值之间相互独立。
线性	因变量与自变量之间为线性相关
同方差性	因变量的方差不随自变量的水平不同而变化。也可称作不变方差。

8.2.1 用lm（）拟合回归模型

myfit <- lm(formula, data)#拟合最基本的函数

其中，fomula指的是要拟合的模型形式，data是一个数据框，表达式formula形式应同

Y ~ X1 + X2 + ... + Xk

~左侧为响应变量，右侧为预测变量

常用的符号

常用的符号

除了这些常用的符号，还有一些函数可以应用于lm()返回的对象来获得更多额外的信息

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

8.2.3 多项式回归

我们可以通过添加一个二项式来提高回归的预测精度

fit2 <- lm(weight ~ height + I(height^2), data=women)

以为^ 在R的表达式中存在特殊的含义，所以使用I()将括号的内同看作R的一个常规表达式

下面是一个拟合含二次项等式的结果

> fit2 <- lm(weight ~ height + I(height^2), data=women)> summary(fit2)
Call:lm(formula = weight ~ height + I(height^2), data = women)
Residuals: Min 1Q Median 3Q Max -0.50941 -0.29611 -0.00941 0.28615 0.59706 
Coefficients:#各项自变量的系数 Estimate Std. Error t value Pr(>|t|) (Intercept) 261.87818 25.19677 10.393 2.36e-07 ***height -7.34832 0.77769 -9.449 6.58e-07 ***I(height^2) 0.08306 0.00598 13.891 9.32e-09 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3841 on 12 degrees of freedomMultiple R-squared: 0.9995, Adjusted R-squared: 0.9994 F-statistic: 1.139e+04 on 2 and 12 DF, p-value: < 2.2e-16
> lm(weight ~ height + I(height^2), data=women)
Call:lm(formula = weight ~ height + I(height^2), data = women)
Coefficients:(Intercept) height I(height^2)  261.87818 -7.34832 0.08306

接着画图

plot(women$height,women$weight, xlab="Height (in inches)", ylab="Weight (in lbs)",main = "forcast")lines(women$height,fitted(fit2))

[R语言实战]第八章回归

一般来说n次多项式生成一个n-1个弯曲的曲线？（不大懂这句话）

car包中的scatterplot函数可以方便地绘制二元关系图。

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

[R语言实战]第八章回归

身高与体重的散点图，直线式线性拟合，虚线为曲线平滑拟合，边界是箱线图

8.2.4 多元线性回归

下面将以state.x77为例子讲解多元线性回归

欲探究一个州的犯罪率和其他因素的关系，包括人口、文盲率、平均收入和结霜天数

lm()对于数据的处理需要数据框，因此我们先创建相关的数据框

> state=data.frame(state.x77[,c("Murder", "Population","Illiteracy", "Income", "Frost")])

在多元回归分析中，我们最好检查一下变量间的相关性。cor()提供了二变量之间的相关系数，car包中的scatterplot()会生成散点图矩阵

> statedf=data.frame(state.x77[,c("Murder", "Population","Illiteracy", "Income", "Frost")])> cor(statedf) Murder Population Illiteracy Income FrostMurder 1.0000000 0.3436428 0.7029752 -0.2300776 -0.5388834Population 0.3436428 1.0000000 0.1076224 0.2082276 -0.3321525Illiteracy 0.7029752 0.1076224 1.0000000 -0.4370752 -0.6719470Income -0.2300776 0.2082276 -0.4370752 1.0000000 0.2262822Frost -0.5388834 -0.3321525 -0.6719470 0.2262822 1.0000000> scatterplotMatrix(statedf)#使用默认参数，函数默认在非对角线区域绘制变量间的散点图，并添加平滑和线性拟合曲线。对角线区域绘制每个变量的密度图和轴须图。

现在使用lm()拟合多元线性回归模型。

> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost,+ data=statedf)> summary(fit)
Call:lm(formula = Murder ~ Population + Illiteracy + Income + Frost,  data = statedf)
Residuals: Min 1Q Median 3Q Max -4.7960 -1.6495 -0.0811 1.4815 7.6210 
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 ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.535 on 45 degrees of freedomMultiple R-squared: 0.567, Adjusted R-squared: 0.5285 F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08

通过分析可以认为，illiteracy在控制其他变量的情况下，每上升1%，谋杀率上升4.14%，同时在p<0.001的情况下显著不为零。

8.2.5 有交互项的多元线性回归

> fit <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)> summary(fit)
Call:lm(formula = mpg ~ hp + wt + hp:wt, data = mtcars)#这里设定了一个马力和车重的交互项
Residuals: Min 1Q Median 3Q Max -3.0632 -1.6491 -0.7362 1.4211 4.5513 
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 49.80842 3.60516 13.816 5.01e-14 ***hp -0.12010 0.02470 -4.863 4.04e-05 ***wt -8.21662 1.26971 -6.471 5.20e-07 ***hp:wt 0.02785 0.00742 3.753 0.000811 ***#显示出马力与车重的交互项是显著的---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.153 on 28 degrees of freedomMultiple R-squared: 0.8848, Adjusted R-squared: 0.8724 F-statistic: 71.66 on 3 and 28 DF, p-value: 2.981e-13

通过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,colors = rainbow(c(2,6,9)))

[R语言实战]第八章回归

通过调整wt的值来展示出不同wt情况下，hp和mpg的相关情况

8.3 回归诊断

上一节的内容中，我们并没有任何函数告知我们这样的模型是不是合适的。对于模型的信心依赖于它多大程度上满足，OLS模型统计假设。

我们现在使用confict()函数来判断一下8.2.4节中statesdf多元回归分析中的结果的适合程度

> statedf<- as.data.frame(state.x77[,c("Murder", "Population",+ "Illiteracy", "Income", "Frost")])> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=statedf)> confint(fit) 2.5 % 97.5 %(Intercept) -6.552191e+00 9.0213182149Population 4.136397e-05 0.0004059867Illiteracy 2.381799e+00 5.9038743192Income -1.312611e-03 0.0014414600Frost -1.966781e-02 0.0208304170

结果表明，文盲率改变1%，谋杀率在95%的置信区间内[2.38,5,90]中变化

回归诊断技术提供了评价回归模型适用性的必要工具。我们接下来先探究R基础包中函数的标准方法，然后看看car包中改进了的新方法

8.3.1 标准方法

最简单的方法是对lm()函数的返回对象使用plot（）生成评价模型拟合情况的四幅图像。

fit <- lm(weight ~ height, data=women)par(mfrow=c(1,4),pin=c(1,1))plot(fit)

理解这些图形应该从OLS回归的统计假设下手

正态性

当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正态分布。“正态Q-Q图”（Normal Q-Q，第二个图）是在正态分布对应的值下，标准化残差的概率图。若满足正态假设，那么图上的应该落在呈45度角的直线上；若不是如此，那么就违反了正态性的假设。

####　独立性

你无法从这些图中分辨出因变量值是否相互独立，只能从收集的数据中来验证。上面的例子中，没有任何先验的理由去相信一位女性的体重会影响另外一位女性的体重。假若你发现数据是从一个家庭抽样得来的，那么可能必须要调整模型独立性的假设。

线性

若因变量与自变量线性相关，那么残差值与预测（拟合）值就没有任何系统关联。换句话说，除了白噪声，模型应该包含数据中所有的系统方差。在“残差图与拟合图”（Residuals vs Fitted，第一个图）中可以清楚地看到一个曲线关系，这暗示着你可能需要对回归模型加上一个二次项

同方差性

若满足不变方差假设，那么在“位置尺度图”（Scale-Location Graph，第三个）中，水平线周围的点应该随机分布。该图似乎满足此假设。

最后一幅残差与杠杆图提供了可能关注的单个观测点信息,可以从图像中鉴别出离群点,高杠杆点,和强影响点。

fit2 <- lm(weight ~ height + I(height^2), data=women)par(mfrow=c(2,2))plot(fit2)

这一组图显示出多项式回归的拟合效果比较理想，基本符合线性假设，残差正态性（第二张图，除了观测点13），和同方差性。从图四看，点15是强影响点。删除观测点13和15我们再进行拟合。

newfit <- lm(weight~ height + I(height^2), data=women[-c(13,15),])plot(newfit)

最后我们用一下这个方法来看states的问题

Rplot01

我们可以看到nevada这个点有比较明显的偏差，可以去除这个点进行模型构建。

8.2.3 改进的方法

car包提供了大量的函数来进行拟合和评估回归模型

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

正态性

qqPlot提供了更加精确的正态假设检验方法

qqPlot()函数提供了更为精确的正态假设检验方法，它画出了在n–p–1个自由度的t分布下的学生化残差

(我在Rtudio上没有试出来他的交互)

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

在图上我们可以很清楚地看到nevada的数据有着明显的偏倚，因此可以去检查nevada的数据

可视化误差还有其他的办法，比如使用residplot函数生成学生化残差柱状图（直方图），添加正态函数，核密度函数和轴须图。

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.误差的独立性

判断因变量或残差是否相互独立，最好是依据收集数据方式的先验知识。

car包里面一个叫做durbinWatsonTest()的函数能够检验误差的序列相关性

> durbinWatsonTest(fit) lag Autocorrelation D-W Statistic p-value 1 -0.2006929 2.317691 0.268 Alternative hypothesis: rho != 0#p值不显著则表明无自相关性，误差项相互独立。

滞后项lag1表明数据中的每个数据都是与其后一个数据中进行比较的。该检验适合时间独立的数据，对于非聚集型数据不适用。注意 durbinWatsonTest()使用自助法来导出p值

3.线性

通过成分残差图可以查看自变量与因变量是否成线性相关，也可以查看是否有不同于已设定线性模型的系统偏差，car包中的crPlots()函数可以完成这一工作

> library(car)> crPlots(fit)

如果发现图形存在非线性，可能表明你对于预测变量的函数形式建模不够充分，需要添加一些其他的曲线成分，比如说多项式项或者对一个或者多个变量进行变换,或者使用其他回归变体而不是线性回归。

4.同方差性

car包提供了两个有用的函数，可以判断误差方差是否恒定。ncvTest()函数生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。

如果检验显著，则说明存在异方差性（误差方差不恒定）

读了好几次都没懂

spreadLevelPlot()函数创建了一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值和拟合值的关系。

> library(car)> ncvTest(fit)Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 1.746514, Df = 1, p = 0.1863

8.3.3 线性模型假设的综合检验

终于重启真的有点莫名的小开心

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

> library(gvlma)> gvmodel <- gvlma(fit)> summary(gvmodel)ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONSUSING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:Level of Significance= 0.05Call:gvlma(x=fit)Value p-value DecisionGlobal Stat 2.773 0.597 Assumptions acceptable.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.

从输出项（Global Stat中的文字栏）我们可以看到数据满足OLS回归模型所有的统计假设（p=0.597）。若Decision下的文字表明违反了假设条件（比如p<0.05），你可以使用前几节讨论的方法来判断哪些假设没有被满足

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)Population Illiteracy Income Frost1.2 2.2 1.3 2.1> sqrt(vif(fit)) > 2 Population Illiteracy Income FrostFALSE FALSE FALSE FALSE

8.4 异常观测值

一个全面的回归分析要覆盖对异常值的分析，包括离群点、高杠杆值点和强影响点

8.4.1 离群点

离群点是指那些模型预测效果不佳的观测点。它们通常有很大的、或正或负的残差,一种鉴别离群点的方法是Q-Q图，落在置信区间带外的点就可以被认为是离群点。

car包提供了一种离群点的统计检验方法，outlierTest()可以求得最大化标准残差的绝对值bonferroni调整后的p值

> library(car)> outlierTest(fit)rstudent unadjusted p-value Bonferonni pNevada 3.5 0.00095 0.048

你可以看到Nevada被判定为离群点（p=0.048）。注意，该函数只是根据单个最大（或正或负）残差值的显著性来判断是否有离群点。若不显著，则说明数据集中没有离群点；若显著，则你必须删除该离群点，然后再检验是否还有其他离群点存在。

8.4.2 高杠杆值点

高杠杆值观测点，即与其他预测变量有关的离群点。换句话说，它们是由许多异常的预测变量值组合起来的，与响应变量值没有关系。

高杠杆值的观测点可通过帽子统计量（hat statistic）判断。对于一个给定的数据集，帽子均值为p/n，其中p是模型估计的参数数目（包含截距项），n是样本量。一般来说，若观测点的帽子值大于帽子均值的2或3倍，就可以认定为高杠杆值点。

statedf=data.frame(state.x77[,c("Murder", "Population","Illiteracy", "Income", "Frost")])fit <- lm(Murder ~ Population + Illiteracy + Income + Frost,data=statedf)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)

交互过程，点击图上的点，全部选择完之后点击esc

水平线标注的即帽子均值2倍和3倍的位置。定位函数（locator function）能以交互模式绘图：单击感兴趣的点，然后进行标注，停止交互时，用户可按Esc键退出，或从图形下拉菜单中选择Stop，或直接右击图形。

高杠杆值点可能是强影响点，也可能不是，这要看它们是否是离群点。

8.4.3 强影响点

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

例如，若移除模型的一个观测点时模型会发生巨大的改变，那么你就需要检测一下数据中是否存在强影响点了.

检测影响点的两种方法：

Cook距离或者称为D统计量

一般来说，Cook’s D值大于4/(n–k–1)，则表明它是强影响点m,其中n为样本量大小，k是预测变量数目

state=data.frame(state.x77[,c("Murder", "Population","Illiteracy", "Income", "Frost")])states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy", "Income", "Frost")])cutoff <- 4/(nrow(states)-length(fit$coefficients)-2)#创建D值plot(fit, which=4, cook.levels=cutoff)abline(h=cutoff, lty=2, col="red")#绘制D值

通过图形可以判断Alaska、Hawaii和Nevada是强影响点。若删除这些点，将会导致回归模型截距项和斜率发生显著变化。注意，虽然该图对搜寻强影响点很有用，但我逐渐发现以1为分割点比4/(n–k–1)更具一般性。若设定D=1为判别标准，则数据集中没有点看起来像是强影响点。

上图虽然有助于鉴别强影响点，但是却不提供如何影响模型的信息。

变量添加图

变量添加图弥补了这个缺陷。对于一个响应变量和k个预测变量，你可以如下图创建k个变量添加图。所谓变量添加图，即对于每个预测变量Xk，绘制Xk在其他k–1个预测变量上回归的残差值相对于响应变量在其他k–1个预测变量上回归的残差值的关系图。car包中的avPlots()函数可提供变量添加图：

library(car)statedf=data.frame(state.x77[,c("Murder", "Population","Illiteracy", "Income", "Frost")])fit <- lm(Murder ~ Population + Illiteracy + Income + Frost,data=statedf)avPlots(fit, ask=FALSE, id.method="identify")

(现在这个包已经可以做到自动鉴定强影响点了)

图中的直线表示相应预测变量的实际回归系数。你可以想象删除某些强影响点后直线的改变，以此来估计它的影响效果。

例如，来看左下角的图（“Murder | others” vs. “Income | others”），若删除点Alaska，直线将往负向移动。事实上，删除Alaska，Income的回归系数将会从0.000 06变为–0.000 85。

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

rm(list=ls())library(car)statedf=data.frame(state.x77[,c("Murder", "Population","Illiteracy", "Income", "Frost")])fit <- lm(Murder ~ Population + Illiteracy + Income + Frost,data=statedf)influencePlot(fit, id.method="identify", main="Influence Plot", sub="Circle size is proportional to Cook's distance")

影响图。纵坐标超过+2或小于–2的州可被认为是离群点，水平轴超过0.2或0.3的州有高杠杆值（通常为预测值的组合）。圆圈大小与影响成比例，圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影响点

不知道为什么同样的代码标不出他的效果

上图反映出反映出Nevada和Rhode Island是离群点，New York、California、Hawaii和Washington有高杠杆值，Nevada、Alaska和Hawaii为强影响点。

8.5 改进措施

8.5.1 删除观测点

删除离群点通常可以提高数据集对于正态假设的拟合度，而强影响点会干扰结果，通常也会被删除。删除最大的离群点或者强影响点后，模型需要重新拟合。若离群点或强影响点仍然存在，重复以上过程直至获得比较满意的拟合。

8.5.2 变量变换

当模型不符合正态性、线性或者同方差性假设时，一个或多个变量的变换通常可以改善或调整模型效果。变换多用Yλ 替代Y，λ的常见值和解释见表8-5。若Y是比例数，通常使用logit变换[ln(Y/1–Y)]。

当模型违反正态假设的时候，通常可以对相应变量尝试某种变化来进行拟合。

car包中的powerTransform()函数通过λ的最大似然估计来正态化变量Xλ

> library(car)> summary(powerTransform(states$Murder))bcPower Transformation to NormalityEst.Power Std.Err. Wald Lower Bound Wald Upper Boundstates$Murder 0.6 0.26 0.088 1.1Likelihood ratio tests about transformation parametersLRT df pvalLR test, lambda=(0) 5.7 1 0.017LR test, lambda=(1) 2.1 1 0.145

结果表明，你可以用Murder0.6 来正态化变量Murder。由于0.6很接近0.5，你可以尝试用平方根变换来提高模型正态性的符合程度。但在本例中，λ=1的假设也无法拒绝（p=0.145），因此没有强有力的证据表明本例需要变量变换，这与图8-9的Q-Q图结果一致。

当违反了线性假设时，对预测变量进行变换常常会比较有用。car包中的boxTidwell()函数通过获得预测变量幂数的最大似然估计来改善线性关系。下面的例子用州的人口和文盲率来预测谋杀率，对模型进行了Box-Tidwell变换：

> library(car)> boxTidwell(Murder~Population+Illiteracy,data=states)Score Statistic p-value MLE of lambdaPopulation -0.32 0.75 0.87Illiteracy 0.62 0.54 1.36

结果显示，使用变换Population0.87 和Illiteracy1.36 能够大大改善线性关系。但是对Population（p=0.75）和Illiteracy（p=0.54）的计分检验又表明变量并不需要变换。

响应变量变换还能改善异方差性（误差方差非恒定）。在代码清单8-7中，你可以看到car包中spreadLevelPlot()函数提供的幂次变换应用，不过，states例子满足了方差不变性，不需要进行变量变换。

8.5.3 增删变量

删除变量在处理多重共线性时是一种非常重要的方法。如果你仅仅是做预测，那么多重共线性并不构成问题，但是如果还要对每个预测变量进行解释，那么就必须解决这个问题。最常见的方法就是删除某个存在多重共线性的变量（某个变量vif >2）.

8.6 选择“最佳”的回归模型

8.6.1 模型比较

用基础安装中的anova()函数可以比较两个嵌套模型的拟合优度。

举个例子：在states的多元回归模型中，我们发现Income和Frost的回归系数不显著，此时你可以检验不含这两个变量的模型与包含这两项的模型预测效果是否一样好

> 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)Analysis of Variance TableModel 1: Murder ~ Population + IlliteracyModel 2: Murder ~ Population + Illiteracy + Income + FrostRes.Df RSS Df Sum of Sq F Pr(>F)1 47 289.2462 45 289.167 2 0.079 0.0061 0.994

在这里进行了是不是应该添加income和forst到线性模型中进行了校验，检验结果不显著，所以可以得出结论，不需要。

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 AICfit1 6 241.6429fit2 4 237.6565

此处AIC值表明没有Income和Frost的模型更佳。注意，ANOVA需要嵌套模型，而AIC方法不需要。

8.6.2 变量选择

从大量候选变量中选择最终的预测变量有以下两种流行的方法

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")#调整向前向后Start: AIC=97.75Murder ~ Population + Illiteracy + Income + FrostDf Sum of Sq RSS AIC- Frost 1 0.02 289.19 95.75- Income 1 0.06 289.22 95.76<none> 289.17 97.75- Population 1 39.24 328.41 102.11- Illiteracy 1 144.26 433.43 115.99Step: AIC=95.75Murder ~ Population + Illiteracy + IncomeDf Sum of Sq RSS AIC- Income 1 0.06 289.25 93.76<none> 289.19 95.75- Population 1 43.66 332.85 100.78- Illiteracy 1 236.20 525.38 123.61Step: AIC=93.76Murder ~ Population + IlliteracyDf Sum of Sq RSS AIC<none> 289.25 93.76- Population 1 48.52 337.76 99.52- Illiteracy 1 299.65 588.89 127.31Call:lm(formula=Murder ~ Population + Illiteracy, data=states)Coefficients:(Intercept) Population Illiteracy1.6515497 0.0002242 4.0807366

开始时模型包含4个（全部）预测变量，然后每一步中，AIC列提供了删除一个行中变量后模型的AIC值，/ 中的AIC值表示没有变量被删除时模型的AIC。第一步，Frost被删除，AIC从97.75降低到95.75；第二步，Income被删除，AIC继续下降，成为93.76。然后再删除变量将会增加AIC，因此终止选择过程。

2. 全子集回归

全子集回归是指所有可能的模型都会被检验,可以选择展示所有可能的结果，也可以展示n个不同子集大小（一个、两个或多个预测变量）的最佳模型（通过参数调整）

全子集回归可用leaps包中的regsubsets()函数实现。你能通过R平方、调整R平方或Mallows Cp统计量等准则来选择“最佳”模型。

R平方含义是预测变量解释响应变量的程度；调整R平方与之类似，但考虑了模型的参数数目。

R平方总会随着变量数目的增加而增加。当与样本量相比，预测变量数目很大时，容易导致过拟合。R平方很可能会丢失数据的偶然变异信息，而调整R平方则提供了更为真实的R平方估计。

Mallows Cp统计量也用来作为逐步回归的判停规则

广泛研究表明，对于一个好的模型，它的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")
library(car)subsets(leaps, statistic="cp",main="Cp Plot for All Subsets Regression")abline(1,1,lty=2,col="red")

第一张图

第一行中（图底部开始），可以看到含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）是最佳模型。

image-20200425174623061

这张图是可以交互的,点个地方出角标

对于不同子集大小，基于Mallows Cp统计量的四个最佳模型。越好的模型离截距项和斜率均为1的直线越近。图形表明，你可以选择这几个模型，其余可能的模型都可以不予考虑：含Population和Illiteracy的双变量模型；含Population、Illiteracy和Frost的三变量模型，或Population、Illiteracy和Income的三变量模型（它们在图形上重叠了，不易分辨）；含Population、Illiteracy、Income和Frost的四变量模型。

大部分情况中，全子集回归要优于逐步回归，因为考虑了更多模型。但是，当有大量预测变量时，全子集回归会很慢。

8.7 深层次分析

8.7.1 交叉验证

从定义来看，回归方法本就是用来从一堆数据中获取最优模型参数。对于OLS回归，通过使得预测误差（残差）平方和最小和对响应变量的解释度（R平方）最大，可获得模型参数。由于等式只是最优化已给出的数据，所以在新数据集上表现并不一定好

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

在k重交叉验证中，样本被分为k个子样本，轮流将k–1个子样本组合作为训练集，另外1个子样本作为保留集。这样会获得k个预测方程，记录k个保留样本的预测表现结果，然后求其平均值。（当n是观测总数目，且k为n时，该方法又称作刀切法，jackknifing。）

bootstrap 包中的crossval() 函数可以实现k 重交叉验证。在代码清单8-15 中，shrinkage()函数对模型的R平方统计量做了k重交叉验证。

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)^2r2cv <- cor(y, results$cv.fit)^2cat("Original R-square =", r2, "\n")cat(k, "Fold Cross-Validated R-square =", r2cv, "\n")cat("Change =", r2-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)Original R-square=0.56710 Fold Cross-Validated R-square=0.4481Change=0.1188

可以看到，基于初始样本的R平方（0.567）过于乐观了。对新数据更好的方差解释率估计是交叉验证后的R平方（0.448）.

通过选择有更好泛化能力的模型，还可以用交叉验证来挑选变量。例如，含两个预测变量（Population和Illiteracy）的模型，比全变量模型R平方减少得更少（0.03 vs. 0.12）

> fit2 <- lm(Murder ~ Population + Illiteracy,data=states)> shrinkage(fit2)Original R-square=0.566832710 Fold Cross-Validated R-square=0.5346871Change=0.03214554

8.7.2 相对重要性

用于评价模型中每一个变量的重要性。

最简单的莫过于比较标准化的回归系数，它表示当其他预测变量不变时，该预测变量一个标准差的变化可引起的响应变量的预期变化（以标准差单位度量）。

> 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)(Intercept) Population Income Illiteracy Frost-9.406e-17 2.705e-01 1.072e-02 6.840e-01 8.185e-03

此处可以看到，当其他因素不变时，文盲率一个标准差的变化将增加0.68个标准差的谋杀率。根据标准化的回归系数，我们可认为Illiteracy是最重要的预测变量，而Frost是最不重要的

相对权重（relative weight）是一种比较有前景的新方法，它是对所有可能子模型添加一个预测变量引起的R平方平均增加量的一个近似值。

vlambda博客 学习文章列表

[R语言实战]第八章 回归

第八章回归

8.1 回归的多面性

8.1.1 OLS回归的适用场景

8.1.2 基础回顾

8.2.1 用lm（）拟合回归模型

8.2.3 多项式回归

8.2.4 多元线性回归

8.2.5 有交互项的多元线性回归

8.3 回归诊断

8.3.1 标准方法

正态性

线性

同方差性

8.2.3 改进的方法

正态性

2.误差的独立性

3.线性

4.同方差性

8.3.3 线性模型假设的综合检验

8.3.4 多重共线性

8.4 异常观测值

8.4.1 离群点

8.4.2 高杠杆值点

8.4.3 强影响点

Cook距离或者称为D统计量

变量添加图

8.5 改进措施

8.5.1 删除观测点

8.5.2 变量变换

8.5.3 增删变量

8.6 选择“最佳”的回归模型

8.6.1 模型比较

8.6.2 变量选择

1.逐步回归

2. 全子集回归

8.7 深层次分析

8.7.1 交叉验证

8.7.2 相对重要性

标签:

推荐阅读

相关文章

vlambda博客
学习文章列表

[R语言实战]第八章回归