专栏首页机器学习与统计学【数据分析 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 ‘ ’ 1
Residual standard error: 233.6 on 15 degrees of freedom
Multiple 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.73
y ~ 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.94
Step:  AIC=223.31
y ~ 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.19
Step:  AIC=222.4
y ~ 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 ‘ ’ 1
Residual standard error: 233.6 on 15 degrees of freedom
Multiple 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 test
data:  y.res
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()

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

本文分享自微信公众号 - 机器学习与统计学(tjxj666)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2015-06-01

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 新恒结衣为什么是中国程序员共同的老婆?

    2001年,姐姐帮她报名参加了《nicola》杂志举办的模特比赛,13岁的新恒结衣意外胜出,开始进军娱乐圈。

    挖数
  • Python数据清洗实践

    “数据科学家们80%的精力消耗在查找、数据清理、数据组织上,只剩于20%时间用于数据分析等。”——IBM数据分析

    AI研习社
  • 经验分享 | 如何写好数据分析师简历?

    我们要确定怎么样简历是一份好数据分析师简历呢?那我们就要涉及到如何评价一个好数据分析师?一般来说,优秀的数据分析师有着很好的表达能力,能通过在二分钟对自己工作内...

    用户2769421
  • Apache Pulsar:灵活的可扩展的批流一体的系统架构

    导读:对于Apache Pulsar,一个经常被问的问题是:Apache Pulsar与现有的消息系统有什么根本的不同。我们之前在文章中介绍了Aache Pul...

    林一
  • 进阶指南 | 如何从数据分析师转型为数据科学家?

    如何从数据分析师华丽转型,成为一名数据科学家?好比“把大象装进冰箱”,成为“数据科学家”仅需简单三步:

    用户2769421
  • 数据化运营02:概念与趋势

    2019年,将针对数据化运营进行一系列的文章总结,期待能够形成一套科学和体系化的方法和指引。而过程中,随着思考和实践的深入,相关的方法论会有优化,甚至推翻重构,...

    用户1756920
  • 中手游《龙珠觉醒》热血重燃,全渠道部署腾讯云

    历时3年研发,中国手游集团(CMGE)超人气日本动漫IP授权大作《龙珠觉醒》烙印着三代龙珠粉的永恒记忆,于2月28日全平台首发上线,全渠道部署腾讯云。腾讯云满载...

    腾讯游戏云
  • 独家数据!《流浪地球》逆袭《新喜剧之王》的真正原因

    可以看到《流浪地球》的票房和排片量基本同步,都在稳步上涨,而《新喜剧之王》的排片量并没有票房下跌得那么厉害。

    挖数
  • 分析完花粥的30首歌词,女流氓名副其实

    花粥是谁?一个在年轻人中传唱度极高的民谣女歌手,93年生人,最早起源于豆瓣,曾好几首歌在豆瓣音乐排行榜上霸榜,而后与同出道于豆瓣的宋东野一起,在国内举办了多场巡...

    挖数
  • 数据分析师必掌握的统计学知识!

    概率是指的对于某一个特定事件的可能性的数值度量,且在0-1之间。我们抛一枚硬币,它有正面朝上和反面朝上两种结果,通常用样本空间S表示,S={正面,反面},而正...

    用户2769421

扫码关注云+社区

领取腾讯云代金券