【数据分析 R语言实战】学习笔记 第九章（中）多元回归分析 回归诊断
【数据分析 R语言实战】学习笔记 第九章（中）多元回归分析 回归诊断
9.2多元线性回归
多元线性回归分析同样由函数lm()完成，但参数formula的表达式应表示为多元形式
例：
财政收入相关统计数据
多元回归模型为：
>revenue=read.table("D:/ProgramFiles/RStudio/revenue.txt",header=T) #读取数据 >lm.reg=lm(y~x1+x2+x3+x4+x5+x6,data=revenue) > summary(lm.reg) Call: lm(formula = y ~ x1 +x2 + x3 + x4 + x5 + x6, data = revenue) Residuals: Min 1Q Median 3Q Max -295.71 -173.52 26.59 90.16 370.01 Coefficients: Estimate Std. Error t valuePr(>|t|) (Intercept) 6.046e+04 3.211e+03 18.829 8.12e-11 *** x1 -1.171e-01 8.638e-02 -1.356 0.19828 x2 3.427e-02 3.322e-02 1.032 0.32107 x3 6.182e-01 4.103e-02 15.067 1.31e-09 *** x4 -5.152e-01 2.930e-02 -17.585 1.91e-10 *** x5 -1.104e-01 2.878e-02 -3.837 0.00206 ** x6 -1.864e-02 1.023e-02 -1.823 0.09143 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1 Residual standarderror: 234.8 on 13 degrees of freedom MultipleR-squared: 0.9999, Adjusted R-squared: 0.9999 F-statistic: 2.294e+04on 6 and 13 DF, p-value: < 2.2e-16
计算结果显示，回归模型的拟合优度0.9999，说明模型的拟合效果较好，但在多元情况下的自变量个数越多，拟合优度会越高，还要看检验的结果:回归方程的F检验一下分显著(p值很小，回归系数x1, x2不显著，x6仅在0.1的显著性水平下显著。可以写出回归方程为
Y=60460-0.1171X1+0.03427X2+0.6182 X3-0.5152X4-0.11 104X5-0.01864X6
与一元线性模型相同，使用函数predict()可以对以后年份的人口增长率作点预测和区间预测。
提取线性拟合模型信息的函数：
在上面的拟合结果中，我们发现自变量x1, x2并不显著，说明第一、二产业国内生产总值对财政收入的解释意义并不显著，应当从模型中剔除，最简单的方式是重写拟合模型
lm.reg=lm(y~x3+x4+x5+x6,data=revenue)
R中的函数update()是专门用于修正模型的函数，在原模型的基础上，不仅可以添加或删除
某些项得到新的模型，还可以对变量进行运算，如对因变量取对数、开方等。其调用格式为
update (object,formula.,…,evaluate=TRUE)
object表示已经拟合好的模型对象，例如存储lm(),glm()的拟合结果;formula指定模型的表达式，原模型中不变的部分用点“.”表示，只写出需要修正的地方即可，例如
Im.reg1=update(lm.reg，~.+x4)表示添加一个新的变量。
lm.reg2=update(lm.reg,sqrt(.)~.)表示对因变量Y作开方运算后再拟合回归模型。
将上例中的x3剔除后重新拟合多元线性回归方程
> lm.reg1=update(lm.reg,.~.-x1-x2)> summary(lm.reg1)Call:lm(formula = y ~ x3 + x4 + x5 + x6, data = revenue)Residuals: Min 1Q Median 3Q Max -325.62 -147.54 14.07 108.28 427.42 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.339e+04 2.346e+03 27.020 3.89e-14 ***x3 6.584e-01 1.548e-02 42.523 < 2e-16 ***x4 -5.438e-01 1.981e-02 -27.445 3.09e-14 ***x5 -1.392e-01 1.918e-02 -7.256 2.80e-06 ***x6 -1.803e-02 9.788e-03 -1.842 0.0854 . ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 233.6 on 15 degrees of freedomMultiple R-squared: 0.9999, Adjusted R-squared: 0.9999 F-statistic: 3.476e+04 on 4 and 15 DF, p-value: < 2.2e-16去除x1 , x2后的方程仍然十分显著，剩余的自变量系数均比较显著，这时拟合的回归方程为Y=63390+0.6584X3-0.5438X4-0.1392X5-0.01803X6
9.2.4逐步回归
在实际分析中，我们使用多元线，性模型描述变量之间的关系时，无法事先了解哪些变量之间的关系显著，就会考虑很多的潜在自变量。若用上一节的方法一一剔除变量，建模过程将变得非常烦琐，所以一般采用逐步回归法。逐步回归建模时，按偏相关系数的大小次序(即变量对y影响程度)将自变量逐个引入方程，对引入的每个自变量的偏相关系数进行统计检验，效应显著的自变量留在回归方程内，如此循此继续遴选下一个自变量。
R中进行逐步回归的函数是step()，以AIC信息准则作为添加或删除变量的判别方法。AIC准则由日本统计学家赤池弘次创立，建立在嫡的概念基础上，一般情况AIC表示为AIC=2(P+1)-2ln(L)
其中，P是回归模型中自变量的个数，L是似然函数。AIC用于寻找解释性最好且包含最少自
由参数的模型，所以选择变量时优先考虑的模型应是AIC值最小的那一个。step()的调用格式为:
step(object, scope, scale = 0,direction = c("both", "backward", "forward"),trace = 1, keep = NULL, steps = 1000, k = 2, ...)其中，object是线性模型或广义线性模型的分析结果;scope确定逐步搜索的范围，是一个公
式或包含upper. lower的列表;direction确定逐步回归的方法，默认值”both”指”一切子集回归法”，"backward”指后退法，"forward”指前进法。trace若是正值，则逐步回归分析的过程将
打印出来:keep是一个过滤器的功能，通常keep选择对象元素的一个子集并返回;steps表示回归的最大步数;k指AIC中的自由度。
对上一节的例子作逐步回归，每一步的分析都将在结果中显示：
> lm.step=step(lm.reg)Start: AIC=223.73y ~ x1 + x2 + x3 + x4 + x5 + x6 Df Sum of Sq RSS AIC- x2 1 58668 775347 223.31<none> 716679 223.73- x1 1 101324 818003 224.38- x6 1 183147 899826 226.28- x5 1 811757 1528436 236.88- x3 1 12515013 13231691 280.05- x4 1 17048151 17764830 285.94Step: AIC=223.31y ~ x1 + x3 + x4 + x5 + x6 Df Sum of Sq RSS AIC- x1 1 43428 818775 222.40<none> 775347 223.31- x6 1 218417 993764 226.27- x5 1 1837987 2613334 245.61- x4 1 22288387 23063734 289.16- x3 1 97699918 98475265 318.19Step: AIC=222.4y ~ x3 + x4 + x5 + x6 Df Sum of Sq RSS AIC<none> 818775 222.40- x6 1 185123 1003898 224.47- x5 1 2873929 3692704 250.52- x4 1 41113643 41932417 299.12- x3 1 98701814 99520589 316.40在默认情况下，逐步回归的每一步过程都会打印出来，本例中逐步回归经历了三步，分别剔
除了不显著的自变量x，和x2 , AIC逐渐减小。最终，R会选择AIC最小的那个模型，即“最优”回归方程。
> summary(lm.step)Call:lm(formula = y ~ x3 + x4 + x5 + x6, data = revenue)Residuals: Min 1Q Median 3Q Max -325.62 -147.54 14.07 108.28 427.42 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.339e+04 2.346e+03 27.020 3.89e-14 ***x3 6.584e-01 1.548e-02 42.523 < 2e-16 ***x4 -5.438e-01 1.981e-02 -27.445 3.09e-14 ***x5 -1.392e-01 1.918e-02 -7.256 2.80e-06 ***x6 -1.803e-02 9.788e-03 -1.842 0.0854 . ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 233.6 on 15 degrees of freedomMultiple R-squared: 0.9999, Adjusted R-squared: 0.9999 F-statistic: 3.476e+04 on 4 and 15 DF, p-value: < 2.2e-16逐步回归直接得到根据AIC选择的最优模型，回归方程所有的检验都是显著的，得到的方程为。
9.3回归诊断及R实现
回归分析完成后，我们仅从显著性检验的角度了解回归效果，但模型的其他特性还有待商榷，例如异常值、共线性等问题，所以我们应该立即进行回归诊断。本节主要包括三个部分的内容:残差诊断、影响分析以及多重共线性诊断。
9.3.1残差诊断
通过分析残差我们可以发现数据是否存在异常值。异常值有两种:一种是“真的”，指由于模型的缺陷、数据违背统计假设、特殊个案等因素形成的异常值;还有一种“假的”的异常值，是由于失误造成的，比如数据录入错误、计算错误、测量错误等。残差也分为几类:普通残差、标准化残差、学生化残差等。
(1)普通残差
利用最小二乘法计算回归模型时，假设中对残差的要求是满足独立性和方差齐性的。所以提取模型残差后，我们要通过画图和检验作残差诊断。残差图是以数据序号或v的拟合值为横坐标、残差为纵坐标的散点图，常见的图形形状：
提取财政收入案例中的模型残差，并绘制残差图，查看残差分布情况
> y.res=lm.reg$residual> y.fit=predict(lm.reg)> plot(y.res~y.fit,main="残差图")
从图中可以看出，残差的分布较为均匀，直观上大致满足模型的假设条件，再利用Shapiro-Wilk检验判断正态性。
> shapiro.test(y.res) Shapiro-Wilk normality testdata: y.resW = 0.96214, p-value = 0.5873从结果可知，p值=0.5873远远大于显著性水平0.05,故不能拒绝原假设，说明数据服从正态分布。
(2)标准化残差
R中用函数rstandard()提取标准化残差，调用格式为
rstandard(modelinfl=lm.influence(model,do.coef=FALSE)sd=sqrt(deviance(model)/df .residual (model))
(3)学生化残差
R中用rstudent()计算标准化残差，调用格式与rstandard()类似，参数res指定模型残差。
rstudent(model，infl=lm.influence(model, do.coef=FALSE)，res=infl$wt.res,…)
9.3.2影响分析
回归诊断要研究的另一个重要问题，是对参数估计或预测值有异常影响的数据，称为强影响
数据。回归模型应当具有一定的稳定性，如果个别一两组数据对估计有异常大的影响，当我们剔除这些数据之后，将得到与原来差异很大的经验回归方程，从而我们将有理由怀疑原回归方程是否真正描述了变量之间的客观存在的相依关系。因此我们要考察每组数据对参数估计的影响大小，这部分内容在}口i归诊断中统称为影响分析。
(1)Leverage
R中计算Leveragel狗函数为hatvalues()和hat()
hatvalues(model, infl = lm.influence(model, do.coef = FALSE), ...)hat(x, intercept = TRUE)其中的参数设置与残差函数类似，model是回归分析Im()返回的对象;x为设计矩阵。
> hii=hatvalues(lm.step)> p=4;n=20> hii>2*(p+1)/n 1 2 3 4 5 6 7 8 9 10 11 12 TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 13 14 15 16 17 18 19 20 FALSE FALSE FALSE FALSE FALSE TRUE FALSE TRUE 在财政收入的回归诊断中，有三个样木点的Leverage高于2(p+1)ln(判断结果为TRUE )，说明可能是异常值点。
( 2 ) DFFITS统计量
> dff=dffits(lm.step)> dff>2*sqrt((p+1)/n) 1 2 3 4 5 6 7 8 9 10 11 12 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 13 14 15 16 17 18 19 20 FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE 根据DFFITS准则，发现第19个样本点的影响较大。
(3) Cook's距离
> cook=cooks.distance(lm.step)> cook>4/n 1 2 3 4 5 6 7 8 9 10 11 12 TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 13 14 15 16 17 18 19 20 FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE 根据Cook's距离判断的异常值点分别是第18,19, 20个样本，与Leverage和DFFITS统计量判断的强影响点一致。
( 4 ) COVRATIO准则
R中用函数covratio()计算COVRATIO统计量。
> covratio(lm.step) 1 2 3 4 5 6 7 0.7418824 1.2893321 1.6064006 0.3699462 1.6073969 0.7233305 1.6635792 8 9 10 11 12 13 14 1.6561510 0.9874913 1.5804153 1.4751229 1.4619963 1.5685943 1.7414143 15 16 17 18 19 20 1.7245615 1.7228882 1.6753871 1.3233294 0.4506737 2.2775795 > options(digits=3)> influence.measures(lm.step)9.3.3多重共线性诊断
多重共线性是指线性回归模型中的解释变量之间由J二存在线性关系或近似线性关系，而使模型难以估计准确，这种现象在经济数据中尤为普遍。一般来说，产生多重共线性的原因主要有三个方面:
经济变量相关的共同趋势。例如经济繁荣期，各基木经济变量都趋于增长。
滞后变量的引入。例如消费受当期收入和前期收入的影响。
样本资料的限制。完全符合理论模型要求的样本数据很难收集，实际样木或多或少都会存在某种程度的共线性。
(1)特征根分析
R中计算矩阵特征根和特征向量的函数是eigen()，调用格式为
eigen(x, symmetric, only.values = FALSE, EISPACK = FALSE)对财政收入案例中的6个变量样本组成的设计矩阵计算特征值:
> options(digits=3)> xx=cor(revenue[3:8])> eigen(xx)$values[1] 4.980418 0.838539 0.162066 0.012847 0.005649 0.000481$vectors [,1] [,2] [,3] [,4] [,5] [,6][1,] -0.445 -0.0735 0.00746 0.8536 0.023853 -0.25950[2,] -0.444 -0.0582 0.26445 -0.0788 0.618251 0.58376[3,] -0.443 -0.0729 0.30977 -0.4538 0.137132 -0.69122[4,] -0.414 -0.1640 -0.86909 -0.2125 0.000528 0.03283[5,] -0.443 -0.0834 0.26783 -0.1168 -0.773369 0.33617[6,] 0.201 -0.9757 0.08372 0.0200 0.017073 0.00212最小的特征根为0.000481，十分接近0,说明6个自变量间一定存在多重共线性，根据对应
的特征向量[,6]还可以写出共线性关系。正是山f这种共线性的存在，逐步回归中系统自动剔除了变量x1、x2
(2)条件数
R软件提供了计算矩阵条件数的函数kappa()，其调用格式为
kappa(z, exact = FALSE, norm = NULL, method = c("qr", "direct"), ...)z为计算的矩阵:exact表示逻辑值，若为TRUE表示精确计算条件数，默认为近似计算;method
指定使用的方法，默认为QR分解。
> kappa(xx)[1] 6132在财政收入的例子中，包含所有变量样木数据的设计矩阵条件数是6132>1000，故认为多重
共线性十分严重。
(3)方差扩大因子(VIF)
> lm.reg=lm(y~x1+x2+x3+x4+x5+x6,data=revenue)> library(DAAG)> vif(lm.reg,digits=3) x1 x2 x3 x4 x5 x6 197.00 778.00 1010.00 10.50 343.00 1.28 计算结果显示，除了X6以外所有变量的方差扩大因子均大于10，说明模型中存在很强的多重共线性。使用方差扩大因子VlF判断共线性还有一个好处是，它可以初步判断l哪些变量之间存在共线性。例如，变量X2和X3方差扩大因子最大，初步判断它们之间可能存在很强的相关性，计算它们的相关系数:
> attach(revenue)> cor(x2,x3)[1] 0.998可知，X2和X3的线性相关系数高达0.998,因此可以肯定它们之间存在严重的共线性，所以上面作逐步回归的过程中就剔除了变量X2
腾讯云开发者
