前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >【数据分析 R语言实战】学习笔记 第九章(中)多元回归分析 回归诊断

【数据分析 R语言实战】学习笔记 第九章(中)多元回归分析 回归诊断

作者头像
统计学家
发布2019-04-10 16:03:47
5K0
发布2019-04-10 16:03:47
举报
文章被收录于专栏:机器学习与统计学

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剔除后重新拟合多元线性回归方程

代码语言:javascript
复制
> lm.reg1=update(lm.reg,.~.-x1-x2)
代码语言:javascript
复制
> summary(lm.reg1)
代码语言:javascript
复制
Call:
代码语言:javascript
复制
lm(formula = y ~ x3 + x4 + x5 + x6, data = revenue)
代码语言:javascript
复制
Residuals:
代码语言:javascript
复制
    Min      1Q  Median      3Q     Max 
代码语言:javascript
复制
-325.62 -147.54   14.07  108.28  427.42 
代码语言:javascript
复制
Coefficients:
代码语言:javascript
复制
              Estimate Std. Error t value Pr(>|t|)    
代码语言:javascript
复制
(Intercept)  6.339e+04  2.346e+03  27.020 3.89e-14 ***
代码语言:javascript
复制
x3           6.584e-01  1.548e-02  42.523  < 2e-16 ***
代码语言:javascript
复制
x4          -5.438e-01  1.981e-02 -27.445 3.09e-14 ***
代码语言:javascript
复制
x5          -1.392e-01  1.918e-02  -7.256 2.80e-06 ***
代码语言:javascript
复制
x6          -1.803e-02  9.788e-03  -1.842   0.0854 .  
代码语言:javascript
复制
---
代码语言:javascript
复制
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
代码语言:javascript
复制
代码语言:javascript
复制
Residual standard error: 233.6 on 15 degrees of freedom
代码语言:javascript
复制
Multiple R-squared:  0.9999,   Adjusted R-squared:  0.9999 
代码语言:javascript
复制
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()的调用格式为:

代码语言:javascript
复制
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中的自由度。

对上一节的例子作逐步回归,每一步的分析都将在结果中显示:

代码语言:javascript
复制
> lm.step=step(lm.reg)
代码语言:javascript
复制
Start:  AIC=223.73
代码语言:javascript
复制
y ~ x1 + x2 + x3 + x4 + x5 + x6
代码语言:javascript
复制
代码语言:javascript
复制
       Df Sum of Sq      RSS    AIC
代码语言:javascript
复制
- x2    1     58668   775347 223.31
代码语言:javascript
复制
<none>                716679 223.73
代码语言:javascript
复制
- x1    1    101324   818003 224.38
代码语言:javascript
复制
- x6    1    183147   899826 226.28
代码语言:javascript
复制
- x5    1    811757  1528436 236.88
代码语言:javascript
复制
- x3    1  12515013 13231691 280.05
代码语言:javascript
复制
- x4    1  17048151 17764830 285.94
代码语言:javascript
复制
代码语言:javascript
复制
Step:  AIC=223.31
代码语言:javascript
复制
y ~ x1 + x3 + x4 + x5 + x6
代码语言:javascript
复制
代码语言:javascript
复制
       Df Sum of Sq      RSS    AIC
代码语言:javascript
复制
- x1    1     43428   818775 222.40
代码语言:javascript
复制
<none>                775347 223.31
代码语言:javascript
复制
- x6    1    218417   993764 226.27
代码语言:javascript
复制
- x5    1   1837987  2613334 245.61
代码语言:javascript
复制
- x4    1  22288387 23063734 289.16
代码语言:javascript
复制
- x3    1  97699918 98475265 318.19
代码语言:javascript
复制
代码语言:javascript
复制
Step:  AIC=222.4
代码语言:javascript
复制
y ~ x3 + x4 + x5 + x6
代码语言:javascript
复制
代码语言:javascript
复制
       Df Sum of Sq      RSS    AIC
代码语言:javascript
复制
<none>                818775 222.40
代码语言:javascript
复制
- x6    1    185123  1003898 224.47
代码语言:javascript
复制
- x5    1   2873929  3692704 250.52
代码语言:javascript
复制
- x4    1  41113643 41932417 299.12
代码语言:javascript
复制
- x3    1  98701814 99520589 316.40

在默认情况下,逐步回归的每一步过程都会打印出来,本例中逐步回归经历了三步,分别剔

除了不显著的自变量x,和x2 , AIC逐渐减小。最终,R会选择AIC最小的那个模型,即“最优”回归方程。

代码语言:javascript
复制
> summary(lm.step)
代码语言:javascript
复制
代码语言:javascript
复制
Call:
代码语言:javascript
复制
lm(formula = y ~ x3 + x4 + x5 + x6, data = revenue)
代码语言:javascript
复制
代码语言:javascript
复制
Residuals:
代码语言:javascript
复制
    Min      1Q  Median      3Q     Max 
代码语言:javascript
复制
-325.62 -147.54   14.07  108.28  427.42 
代码语言:javascript
复制
代码语言:javascript
复制
Coefficients:
代码语言:javascript
复制
              Estimate Std. Error t value Pr(>|t|)    
代码语言:javascript
复制
(Intercept)  6.339e+04  2.346e+03  27.020 3.89e-14 ***
代码语言:javascript
复制
x3           6.584e-01  1.548e-02  42.523  < 2e-16 ***
代码语言:javascript
复制
x4          -5.438e-01  1.981e-02 -27.445 3.09e-14 ***
代码语言:javascript
复制
x5          -1.392e-01  1.918e-02  -7.256 2.80e-06 ***
代码语言:javascript
复制
x6          -1.803e-02  9.788e-03  -1.842   0.0854 .  
代码语言:javascript
复制
---
代码语言:javascript
复制
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
代码语言:javascript
复制
代码语言:javascript
复制
Residual standard error: 233.6 on 15 degrees of freedom
代码语言:javascript
复制
Multiple R-squared:  0.9999,   Adjusted R-squared:  0.9999 
代码语言:javascript
复制
F-statistic: 3.476e+04 on 4 and 15 DF,  p-value: < 2.2e-16

逐步回归直接得到根据AIC选择的最优模型,回归方程所有的检验都是显著的,得到的方程为。

9.3回归诊断及R实现

回归分析完成后,我们仅从显著性检验的角度了解回归效果,但模型的其他特性还有待商榷,例如异常值、共线性等问题,所以我们应该立即进行回归诊断。本节主要包括三个部分的内容:残差诊断、影响分析以及多重共线性诊断。

9.3.1残差诊断

通过分析残差我们可以发现数据是否存在异常值。异常值有两种:一种是“真的”,指由于模型的缺陷、数据违背统计假设、特殊个案等因素形成的异常值;还有一种“假的”的异常值,是由于失误造成的,比如数据录入错误、计算错误、测量错误等。残差也分为几类:普通残差、标准化残差、学生化残差等。

(1)普通残差

利用最小二乘法计算回归模型时,假设中对残差的要求是满足独立性和方差齐性的。所以提取模型残差后,我们要通过画图和检验作残差诊断。残差图是以数据序号或v的拟合值为横坐标、残差为纵坐标的散点图,常见的图形形状:

提取财政收入案例中的模型残差,并绘制残差图,查看残差分布情况

代码语言:javascript
复制
> y.res=lm.reg$residual
代码语言:javascript
复制
> y.fit=predict(lm.reg)
代码语言:javascript
复制
> plot(y.res~y.fit,main="残差图")

从图中可以看出,残差的分布较为均匀,直观上大致满足模型的假设条件,再利用Shapiro-Wilk检验判断正态性。

代码语言:javascript
复制
> shapiro.test(y.res)
代码语言:javascript
复制
代码语言:javascript
复制
        Shapiro-Wilk normality test
代码语言:javascript
复制
代码语言:javascript
复制
data:  y.res
代码语言:javascript
复制
W = 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()

代码语言:javascript
复制
hatvalues(model, infl = lm.influence(model, do.coef = FALSE), ...)
代码语言:javascript
复制
hat(x, intercept = TRUE)

其中的参数设置与残差函数类似,model是回归分析Im()返回的对象;x为设计矩阵。

代码语言:javascript
复制
> hii=hatvalues(lm.step)
代码语言:javascript
复制
> p=4;n=20
代码语言:javascript
复制
> hii>2*(p+1)/n
代码语言:javascript
复制
    1     2     3     4     5     6     7     8     9    10    11    12 
代码语言:javascript
复制
 TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
代码语言:javascript
复制
   13    14    15    16    17    18    19    20 
代码语言:javascript
复制
FALSE FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE 

在财政收入的回归诊断中,有三个样木点的Leverage高于2(p+1)ln(判断结果为TRUE ),说明可能是异常值点。

( 2 ) DFFITS统计量

代码语言:javascript
复制
> dff=dffits(lm.step)
代码语言:javascript
复制
> dff>2*sqrt((p+1)/n)
代码语言:javascript
复制
    1     2     3     4     5     6     7     8     9    10    11    12 
代码语言:javascript
复制
FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
代码语言:javascript
复制
   13    14    15    16    17    18    19    20 
代码语言:javascript
复制
FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE 

根据DFFITS准则,发现第19个样本点的影响较大。

(3) Cook's距离

代码语言:javascript
复制
> cook=cooks.distance(lm.step)
代码语言:javascript
复制
> cook>4/n
代码语言:javascript
复制
    1     2     3     4     5     6     7     8     9    10    11    12 
代码语言:javascript
复制
 TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
代码语言:javascript
复制
   13    14    15    16    17    18    19    20 
代码语言:javascript
复制
FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE 

根据Cook's距离判断的异常值点分别是第18,19, 20个样本,与Leverage和DFFITS统计量判断的强影响点一致。

( 4 ) COVRATIO准则

R中用函数covratio()计算COVRATIO统计量。

代码语言:javascript
复制
> covratio(lm.step)
代码语言:javascript
复制
        1         2         3         4         5         6         7 
代码语言:javascript
复制
0.7418824 1.2893321 1.6064006 0.3699462 1.6073969 0.7233305 1.6635792 
代码语言:javascript
复制
        8         9        10        11        12        13        14 
代码语言:javascript
复制
1.6561510 0.9874913 1.5804153 1.4751229 1.4619963 1.5685943 1.7414143 
代码语言:javascript
复制
       15        16        17        18        19        20 
代码语言:javascript
复制
1.7245615 1.7228882 1.6753871 1.3233294 0.4506737 2.2775795 
代码语言:javascript
复制
> options(digits=3)
代码语言:javascript
复制
> influence.measures(lm.step)
代码语言:javascript
复制

9.3.3多重共线性诊断

多重共线性是指线性回归模型中的解释变量之间由J二存在线性关系或近似线性关系,而使模型难以估计准确,这种现象在经济数据中尤为普遍。一般来说,产生多重共线性的原因主要有三个方面:

经济变量相关的共同趋势。例如经济繁荣期,各基木经济变量都趋于增长。

滞后变量的引入。例如消费受当期收入和前期收入的影响。

样本资料的限制。完全符合理论模型要求的样本数据很难收集,实际样木或多或少都会存在某种程度的共线性。

(1)特征根分析

R中计算矩阵特征根和特征向量的函数是eigen(),调用格式为

代码语言:javascript
复制
eigen(x, symmetric, only.values = FALSE, EISPACK = FALSE)

对财政收入案例中的6个变量样本组成的设计矩阵计算特征值:

代码语言:javascript
复制
> options(digits=3)
代码语言:javascript
复制
> xx=cor(revenue[3:8])
代码语言:javascript
复制
> eigen(xx)
代码语言:javascript
复制
$values
代码语言:javascript
复制
[1] 4.980418 0.838539 0.162066 0.012847 0.005649 0.000481
代码语言:javascript
复制
代码语言:javascript
复制
$vectors
代码语言:javascript
复制
       [,1]    [,2]     [,3]    [,4]      [,5]     [,6]
代码语言:javascript
复制
[1,] -0.445 -0.0735  0.00746  0.8536  0.023853 -0.25950
代码语言:javascript
复制
[2,] -0.444 -0.0582  0.26445 -0.0788  0.618251  0.58376
代码语言:javascript
复制
[3,] -0.443 -0.0729  0.30977 -0.4538  0.137132 -0.69122
代码语言:javascript
复制
[4,] -0.414 -0.1640 -0.86909 -0.2125  0.000528  0.03283
代码语言:javascript
复制
[5,] -0.443 -0.0834  0.26783 -0.1168 -0.773369  0.33617
代码语言:javascript
复制
[6,]  0.201 -0.9757  0.08372  0.0200  0.017073  0.00212

最小的特征根为0.000481,十分接近0,说明6个自变量间一定存在多重共线性,根据对应

的特征向量[,6]还可以写出共线性关系。正是山f这种共线性的存在,逐步回归中系统自动剔除了变量x1、x2

(2)条件数

R软件提供了计算矩阵条件数的函数kappa(),其调用格式为

代码语言:javascript
复制
kappa(z, exact = FALSE,
代码语言:javascript
复制
      norm = NULL, method = c("qr", "direct"), ...)

z为计算的矩阵:exact表示逻辑值,若为TRUE表示精确计算条件数,默认为近似计算;method

指定使用的方法,默认为QR分解。

代码语言:javascript
复制
> kappa(xx)
代码语言:javascript
复制
[1] 6132

在财政收入的例子中,包含所有变量样木数据的设计矩阵条件数是6132>1000,故认为多重

共线性十分严重。

(3)方差扩大因子(VIF)

代码语言:javascript
复制
> lm.reg=lm(y~x1+x2+x3+x4+x5+x6,data=revenue)
代码语言:javascript
复制
> library(DAAG)
代码语言:javascript
复制
> vif(lm.reg,digits=3)
代码语言:javascript
复制
     x1      x2      x3      x4      x5      x6 
代码语言:javascript
复制
 197.00  778.00 1010.00   10.50  343.00    1.28 

计算结果显示,除了X6以外所有变量的方差扩大因子均大于10,说明模型中存在很强的多重共线性。使用方差扩大因子VlF判断共线性还有一个好处是,它可以初步判断l哪些变量之间存在共线性。例如,变量X2和X3方差扩大因子最大,初步判断它们之间可能存在很强的相关性,计算它们的相关系数:

代码语言:javascript
复制
> attach(revenue)
代码语言:javascript
复制
> cor(x2,x3)
代码语言:javascript
复制
[1] 0.998

可知,X2和X3的线性相关系数高达0.998,因此可以肯定它们之间存在严重的共线性,所以上面作逐步回归的过程中就剔除了变量X2

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2015-06-01,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 机器学习与统计学 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档