《R语言学习系列32-回归分析(共26页).docx》由会员分享，可在线阅读，更多相关《R语言学习系列32-回归分析(共26页).docx（26页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上27. 回归分析回归分析是研究一个或多个变量（因变量）与另一些变量（自变量）之间关系的统计方法。主要思想是用最小二乘法原理拟合因变量与自变量间的最佳回归模型（得到确定的表达式关系）。其作用是对因变量做解释、控制、或预测。回归与拟合的区别：拟合侧重于调整曲线的参数，使得与数据相符；而回归重在研究两个变量或多个变量之间的关系。它可以用拟合的手法来研究两个变量的关系，以及出现的误差。 回归分析的步骤： （1）获取自变量和因变量的观测值；（2）绘制散点图，并对异常数据做修正；（3）写出带未知参数的回归方程；（4）确定回归方程中参数值；（5）假设检验，判断回归方程的拟合优度；（
2、6）进行解释、控制、或预测。（一）一元线性回归一、原理概述 1. 一元线性回归模型： Y=𝛽0+𝛽1X+其中 X是自变量，Y是因变量，𝛽0,𝛽1是待求的未知参数，𝛽0也称为截距；是随机误差项，也称为残差，通常要求满足： 的均值为0； 的方差为𝜎2； 协方差COV(i, j)=0，当ij时。即对所有的ij, i与j 互不相关。用最小二乘法原理，得到最佳拟合效果的值：， 2.模型检验(1) 拟合优度检验计算R2，反映了自变量所能解释的方差占总方差的百分比，值越大说明模型拟合效果越好。通常可以认为当R2大
3、于0.9时，所得到的回归直线拟合得较好，而当R2小于0.5时，所得到的回归直线很难说明变量之间的依赖关系。(2) 回归方程参数的检验回归方程反应了因变量Y随自变量X变化而变化的规律，若𝛽1=0，则Y不随X变化，此时回归方程无意义。所以，要做如下假设检验：H0: 𝛽1=0, H1: 𝛽10； F检验 若𝛽1=0为真，则回归平方和RSS与残差平方和ESS/(N-2)都是𝜎2的无偏估计，因而采用F统计量：来检验原假设1=0是否为真。 T检验对H0: 𝛽1=0的T检验与F检验是等价的（t2=F）。3. 用回
4、归方程做预测得到回归方程后，预测X=x0处的Y值.的预测区间为：其中t/2的自由度为N-2. 二、R语言实现使用lm()函数实现，基本格式为：lm(formula, data, subset, weights, na.action,method=qr, .)其中，formula为要拟合的回归模型的形式，一元线性回归的格式为：yx，y表示因变量，x表示自变量，若不想包含截距项，使用yx-1；data为数据框或列表；subset选取部分子集；weights取NULL时表示最小二乘法拟合，若取值为权重向量，则用加权最小二乘法；na.action设定是否忽略缺失值；method指定拟合的方法，目前只支
5、持“qr”（QR分解），method=“model.frame”返回模型框架。三、实例例1 现有埃及卡拉马村庄每月记录儿童身高的数据，做一元线性回归。datas-data.frame(age=18:29,height=c(76.1,77,78.1,78.2,78.8,79.7,79.9,81.1,81.2,81.8,82.8,83.5)datas age height1 18 76.12 19 77.03 20 78.14 21 78.25 22 78.86 23 79.77 24 79.98 25 81.19 26 81.210 27 81.811 28 82.812 29 83.5plot
6、(datas) #绘制散点图res.reg|t|) (Intercept) 64.9283 0.5084 127.71 2e-16 *age 0.6350 0.0214 29.66 4.43e-11 *-Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 0.256 on 10 degrees of freedomMultiple R-squared: 0.9888,Adjusted R-squared: 0.9876 F-statistic: 880 on 1 and 10 DF, p-value: 4.
7、428e-11说明：输出了残差信息Residuals；回归系数估计值、标准误、t统计量值、p值，可得到回归方程：height=64.9283+0.6350*age回归系数p值（0.5, 表示拟合程度很好。F统计量=880，p值=4.428e-11远小于0.05，表示整个回归模型显著，适合估计height这一因变量。coefficients(res.reg) #返回模型的回归系数估计值(Intercept) age 64. 0.confint(res.reg,parm=age,level=0.95) #输出参数age的置信区间，若不指定parm将返回所有参数的置信区间 2.5 % 97.5 %a
8、ge 0. 0.fitted(res.reg) #输出回归模型的预测值 1 2 3 4 5 6 7 8 9 10 11 1276.35769 76.99266 77.62762 78.26259 78.89755 79.53252 80.16748 80.80245 81.43741 82.07238 82.70734 83.34231anova(res.reg) #输出模型的方差分析表Response: height Df Sum Sq Mean Sq F value Pr(F) age 1 57.655 57.655 879.99 4.428e-11 *Residuals 10 0.655
9、 0.066 -Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1vcov(res.reg) #输出模型的协方差矩阵 (Intercept) age(Intercept) 0. -0.age -0. 0.residuals(res.reg) #输出模型的残差 1 2 3 4 5 6 7 8 9 10 11 12-0. 0. 0. -0. -0. 0. -0. 0. -0. -0. 0. 0.AIC(res.reg) #输出模型的AIC值1 5.BIC(res.reg) #输出模型的BIC值1 6.logLik(res.reg) #输出模型的对数似然值lo
10、g Lik. 0. (df=3)abline(res.reg) #给散点图加上一条回归线par(mfrow=c(2,2)plot(res.reg) #绘制回归诊断图说明：分别是残差与拟合值图，二者越无关联越好，若有明显的曲线关系，则说明需要对线性回归模型加上高次项；残差的Q-Q图，看是否服从正态分布；标准化残差与拟合值图，也叫位置-尺度图，纵坐标是标准化残差的平方根，残差越大，点的位置越高，用来判断模型残差是否等方差，若满足则水平线周围的点应随机分布；残差与杠杆图，虚线表示Cooks距离（每个数据点对回归线的影响力）等高线，从中可以鉴别出离群点（第3个点较大，表示删除该数据点，回归系数将有实质
11、上的改变，为异常值点）、高杠杆点、强影响点。datas-datas-3, #删除第3个样本点，重新做一元线性回归res.reg2-lm(heightage,datas)summary(res.reg2) 新的回归方程为：height=64.5540+0.6489*age，拟合优度R2=0.993，拟合效果变得更好。 #用回归模型预测ages-data.frame(age=30:34)pre.res-predict(res.reg2,ages,interval=prediction,level=0.95) #注意predict函数的第1个参数必须是回归模型的自变量数据构成的数据框或列表pre.r
12、es fit lwr upr1 84.02034 83.46839 84.572282 84.66921 84.09711 85.241323 85.31809 84.72365 85.912544 85.96697 85.34825 86.585695 86.61585 85.97114 87.26056（二）多元线性回归一、基本原理 1. 多元线性回归模型：Y=𝛽0+𝛽1X1+ 𝛽NXN+其中 X1, , XN是自变量，Y是因变量，𝛽0, 𝛽1, 𝛽N是待求的未知参数，是随机误差项（残差），若记
13、多元线性回归模型可写为矩阵形式：Y=X+通常要求：矩阵X的秩为k+1（保证不出现共线性）, 且k|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 *
14、 x6 -1.864e-02 1.023e-02 -1.823 0.09143 . -Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 234.8 on 13 degrees of freedomMultiple R-squared: 0.9999,Adjusted R-squared: 0.9999 F-statistic: 2.294e+04 on 6 and 13 DF, p-value: 2.2e-16 说明：拟合优度R2=0.9999, 效果非常好。但是多元回归时，自变量个数越多，拟合优度必然越
15、好，所以还要看检验效果和回归系数是否显著。结果解释、回归方程、回归预测与前文类似（略）。结合显著性代码可看出：x1和x2不显著，x6只在0.1显著水平下显著，故应考虑剔除x1和x2.R语言中提供了update()函数，用来在原模型的基础上进行修正，还可以对变量进行运算，其基本格式为：update(object, formula., ., evaluate=TRUE)其中，object为前面拟合好的原模型对象；formula指定模型的格式，原模型不变的部分用“.”表示，只写出需要修正的地方即可，例如update(lm.reg, .+x7) 表示添加一个新的变量update(lm.reg, sqr
16、t(.).) 表示对因变量y开方，再重新拟合回归模型lm.reg2|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
17、 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（三）逐步回归多元线性回归模型中，并不是所有的自变量都与因变量有显著关系，有时有些自变量的作用可以忽略。这就需要考虑怎样从所有可能有关的自变量中挑选出对因变量有显著影响的部分自变量。逐步回归的基本思想是，将变量一个一个地引入或剔出，引入或剔出变量的条件是“偏相关系数”经检验是显著的，同时
18、每引入或剔出一个变量后，对已选入模型的变量要进行逐个检验，将不显著变量剔除或将显著的变量引入，这样保证最后选入的所有自变量都是显著的。逐步回归每一步只有一个变量引入或从当前的回归模型中剔除，当没有回归因子能够引入或剔出模型时，该过程停止。R语言中，用step()函数进行逐步回归，以AIC信息准则作为选入和剔除变量的判别条件。AIC是日本统计学家赤池弘次，在熵概念的基础上建立的：AIC=2(p+1)-2ln(L)其中，p为回归模型的自变量个数，L是似然函数。注：AIC值越小越被优先选入。基本格式：step(object, direction=,steps=, k=2, .)其中，object是线
19、性模型或广义线性模型的返回结果；direction确定逐步回归的方法，默认“both”综合向前向后法，“backward”向后法（先把全部自变量加入模型，若无统计学意义则剔出模型），“forward”向前法（先将部分自变量加入模型，再逐个添加其它自变量，若有统计学意义则选入模型）；steps表示回归的最大步数，默认1000；k默认=2, 输出为AIC值，= log(n)有时输出BIC或SBC值。 另外，有时还需要借助使用drop1(object)和add1(object)函数，其中object为逐步回归的返回结果，判断剔除或选入一个自变量，AIC值的变化情况，以筛选选入模型的自变量。lm.st
20、ep|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 de
21、grees 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说明：默认输出每步的结果（略），进行了3步回归，逐步剔除了自变量x1和x2，AIC值逐步减小，最终得到最优的模型。drop1(lm.step)Single term deletionsModel:y x3 + x4 + x5 + x6 Df Sum of Sq RSS AIC 222.40x3 1 316.40x4 1 299.12x5 1 250
22、.52x6 1 224.47lm.reg3|t|) (Intercept) 6.284e+04 2.494e+03 25.191 2.66e-14 *x3 6.614e-01 1.651e-02 40.066 2e-16 *x4 -5.467e-01 2.118e-02 -25.813 1.81e-14 *x5 -1.412e-01 2.053e-02 -6.877 3.72e-06 *-Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 250.5 on 16 degrees of freedomMult
23、iple R-squared: 0.9999,Adjusted R-squared: 0.9998 F-statistic: 4.032e+04 on 3 and 16 DF, p-value: 2.2e-16 说明：使用drop1()函数考察分别剔除每个自变量，AIC值变化的情况，可以看出不剔除x6与剔除x6，AIC值只从222.40变大到224.47，相对其它自变量变化很小。所以，可以考虑剔除掉x6，重新做多元线性回归。（四）回归诊断回归分析之后，还需要从残差的随机性、强影响分析、共线性方面进行诊断。一、残差诊断(1) 残差y.res-lm.reg3$residuals #回归模型的残差y
24、.fit0.05, 接受原假设，即残差服从正态分布。(2) 标准化残差残差与数据的数量级有关，除以标准误差后得到标准化残差。理想的标准化残差服从N(0,1).rs-rstandard(lm.reg3) #得到标准化残差plot(rsy.fit,main=标准残差图)shapiro.test(rs) #标准化残差的正态性检验Shapiro-Wilk normality testdata: rsW = 0.97766, p-value = 0.9004(3) 学生化残差为了回避标准化残差的方差齐性假设，使用学生化残差。rst-rstudent(lm.reg3)plot(rsy.fit,main=学
25、生化残差图)shapiro.test(rst)Shapiro-Wilk normality testdata: rstW = 0.97463, p-value = 0.848(4) 残差自相关性的Durbin-Watson检验使用car包中的函数：durbinWatsonTest(model,alternative=c(two.side, positive, negative)H0：序列不存在自相关性library(car)durbinWatsonTest(lm.reg3) lag Autocorrelation D-W Statistic p-value 1 -0. 2.42579 0.77
26、 Alternative hypothesis: rho != 0二、强影响分析对参数估计或预测值有异常影响的数据，称为强影响数据。回归模型应当具有一定的稳定性，若个别一两组数据对估计有异常大的影响，剔除后将得到与原来差异很大的回归方程，从而有理由怀疑原回归方程是否真正描述了变量间的客观存在的关系。1. 反映这种强影响的统计量有4种及函数：Leveragehatvalues(model)DEFITSdffits(model)Cooks距离cooks.distance(model) COVRATIOcovratio(model)另外，influence.measures(model)函数，可以汇
27、总上述4种统计量，判断强影响点。influence.measures(lm.reg3) Influence measures of lm(formula = y x3 + x4 + x5, data = revenue) : dfb.1_ dfb.x3 dfb.x4 dfb.x5 dffit cov.r cook.d hat inf18 0. -3.04124 -0. 2. -3.40945 0.810 2.14e+00 0.6347 *19 0. -0.09558 -0. 0. 1.61704 0.515 5.04e-01 0.3127 *20 -1. 1.56506 1. -1. -3.3
28、3452 1.426 2.25e+00 0.6996 * 说明：判断出第18, 19, 20个样本是强影响点。2. Bonferroni离群点检验使用car包中的函数outlierTest(model)library(car)outlierTest(lm.reg3)No Studentized residuals with Bonferonni p 0.05Largest |rstudent|: rstudent unadjusted p-value Bonferonni p18 -2. 0.02064 0.4128注：去掉强影响点，重新做多元线性回归（略）。三、共线性诊断回归分析中很容易发生
29、模型中两个或两个以上的自变量高度相关，从而引起最小二乘估计可能很不精确（称为共线性问题）。在实际中最常见的问题是一些重要的自变量很可能由于在假设检验中t值不显著而被不恰当地剔除了。共线性诊断问题就是要找出哪些变量间存在共线性关系。 1. 模型条件数检验使用函数kappa(z, exact=FALSE, )，其中，z为矩阵XTX，或lm、glm的返回对象；exact设置是否精确计算。 一般认为：当K100时不存在多重共线性；当100K1000时存在较强的多重共线性；当K1000时存在严重的多重共线性。x-scale(revenue,3:8) #取出自变量数据，做标准化xx=crossprod(x
30、) #求xx即矩阵的叉积kappa(xx)1 6132.142 2. 方差膨胀因子（VIF）检验使用car包中的函数vif(model)，该函数还能判断哪些自变量间存在共线性。一般认为：当vif10时不存在多重共线性；当10vif100时，存在较强的多重共线性；当vif100时存在严重的多重共线性。lm.reg0), 则接近奇异的程度就会比小很多，从而消除了多重共线性。考虑到变量的量纲，应先对数据进行标准化，仍记为X，称为的岭回归估计，k称为岭参数。显然岭回归估计的值比最小二乘估计值稳定，当k=0时，岭回归就是普通的最小二乘估计。由于岭参数k不是唯一确定的，故岭回归估计实际上是回归参数的一个估
31、计族（关于k的函数），该函数曲线称为岭迹。根据岭迹图可以选择合适的k值，称为岭迹法。其一般原则是： (1) 各回归系数的岭估计基本稳定；(2) 最小二乘估计的回归系数符合不合理时，岭估计参数的符合变的合理；(3) 回归系数每月不合乎实际意义的绝对值；(4) 残差平方和增加不太多。岭回归用MASS包中的函数lm.ridge()实现，基本格式为：lm.ridge(formula, data, subset, na.action, lambda=0,.)其中，formula为回归公式的形式，data为数据框，subset设置其子集；lambda为岭参数的标量或矢量。例3 商品销售量y，考虑4个自变量：可支配收入x1, 商品价格指数x2, 商品的社会保有量x3, 其它消费品平均价格指数x4