R语言学习系列32

27.回归分析回归分析是研究一个或多个变量（因变量）与另一些变量（自变量）之间关系的统计方法。主要思想是用最小二乘法原理拟合因变量与自变量间的最佳回归模型（得到确定的表达式关系）。其作用是对因变量做解释、控制、或预测。回归与拟合的区别：拟合侧重于调整曲线的参数，使得与数据相符；而回归重在研究两个变量或多个变量之间的关系。它可以用拟合的手法来研究两个变量的关系，以及出现的误差。回归分析的步骤：（1）获取自变量和因变量的观测值；（2）绘制散点图，并对异常数据做修正；（3）写出带未知参数的回归方程；（4）确定回归方程中参数值；（5）假设检验，判断回归方程的拟合优度；（6）进行解释、控制、或预测。（一）一元线性回归一、原理概述一元线性回归模型：Y=o+1X+£其中X是自变量，Y是因变量，0,1是待求的未知参数，0也称为截距；£是随机误差项，也称为残差，通常要求£满足:£的均值为0；£的方差为2；③协方差COV(£,£.)=0，当i的时。即对所有的i却£与£互不1j1j相关。用最小二乘法原理，得到最佳拟合效果的B,B值：o1E3-x)(y-y)p=^-^=!，B=y-Bx1厂o1乙(X一X)2ii=1模型检验拟合优度检验计算R2,反映了自变量所能解释的方差占总方差的百分比，值越大说明模型拟合效果越好。通常可以认为当R2大于0.9时，所得到的回归直线拟合得较好，而当R2小于0.5时，所得到的回归直线很难说明变量之间的依赖关系。回归方程参数的检验回归方程反应了因变量Y随自变量X变化而变化的规律，若广0，则Y不随X变化，此时回归方程无意义。所以，要做如下假设检验：H0:广0,H1：1#0；①F检验若广0为真，则回归平方和RSS与残差平方和ESS/(N-2)都是2的无偏估计，因而采用F统计量：来检验原假设禹=0是否为真。②T检验对H0：1=0的T检验与F检验是等价的(t2=F)。用回归方程做预测得到回归方程Y=B+BX后，预测X=x处的Y值y=B+BX.0i000ioy0的预测区间为：其中t,，的自由度为N-2.a/2二、R语言实现使用lm()函数实现，基本格式为：lm(formula,data,subset,weights,na.action,

method="qr",...)其中，formula为要拟合的回归模型的形式，一元线性回归的格式为:y〜x，y表示因变量，x表示自变量，若不想包含截距项，使用y〜x-1；data为数据框或列表；subset选取部分子集；weights取NULL时表示最小二乘法拟合，若取值为权重向量，则用加权最小二乘法；na.action设定是否忽略缺失值；method指定拟合的方法，目前只支持“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))datasageheightTOC\o"1-5"\h\z1876.11977.02078.12178.22278.82379.72479.92581.12681.22781.82882.82983.5plot(datas)#绘制散点图res.reg<-lm(height~age,datas)#做一元线性回归summary(res.reg)Residuals:summary(res.reg)Residuals:#输出模型的汇总结果Min1QMedian3QMax-0.27238-0.24248-0.027620.160140.47238Coefficients:tvaluePr(>|t|)tvaluePr(>|t|)127.71<2e-16***29.664.43e-11***(Intercept)64.92830.5084age0.63500.0214Signif.codes:0‘***’0.001‘心’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:0.256on10degreesoffreedomMultipleR-squared:0.9888,AdjustedR-squared:0.9876F-statistic:880on1and10DF,p-value:4.428e-11说明：输出了残差信息Residuals；回归系数估计值、标准误、t统计量值、p值，可得到回归方程:height=64.9283+0.6350*age回归系数p值(<2e-16，4.43e-11)很小，非常显著的手0***也表示显著程度非常显著。拟合优度R2=0.9888>0.5,表示拟合程度很好。F统计量=880，p值=4.428e-11远小于0.05,表示整个回归模型显著，适合估计height这一因变量。coefficients(res.reg)#返回模型的回归系数估计值(Intercept)age64.9283220.634965confint(res.reg,parm="age”,level=0.95)#输出参数age的置信区间，若不指定parm将返回所有参数的置信区间2.5%97.5%age0.58727220.6826578fitted(res.reg)#输出回归模型的预测值12345678910111276.3576976.9926677.6276278.2625978.8975579.5325280.1674880.8024581.4374182.0723882.7073483.34231anova(res.reg)#输出模型的方差分析表Response:heightDfSumSqMeanSqFvaluePr(>F)age157.65557.655879.994.428e-11***Residuals100.6550.066—Signif.codes:0‘***’0.001‘心’0.01‘*’0.05‘.’0.1‘’1vcov(res.reg)#输出模型的协方差矩阵(Intercept)age(Intercept)0.-0.00age-0.010766860.0004581642residuals(res.reg)#输出模型的残差123456789101112-0.80.0073426570.2-0.0-0.00.7-0.70.8-0.7-0.20.00.8AlC(res.reg)#输出模型的AIC值5.161407BlC(res.reg)#输出模型的BIC值6.616127logLik(res.reg)#输出模型的对数似然值'logLik.'0.4192965(df=3)abline(res.reg)#给散点图加上一条回归线par(mfrow=c(2,2))plot(res.reg)#绘制回归诊断图说明：分别是残差与拟合值图，二者越无关联越好，若有明显的曲线关系，则说明需要对线性回归模型加上高次项；残差的Q-Q图，看是否服从正态分布；标准化残差与拟合值图，也叫位置-尺度图，纵坐标是标准化残差的平方根，残差越大，点的位置越高，用来判断模型残差是否等方差，若满足则水平线周围的点应随机分布；残差与杠杆图，虚线表示Cooks距离(每个数据点对回归线的影响力)等高线，从中可以鉴别出离群点(第3个点较大，表示删除该数据点，回归系数将有实质上的改变，为异常值点)、高杠杆点、强影响点。datas<-datas[-3,]#删除第3个样本点，重新做一元线性回归res.reg2<-lm(height~age,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.resfitlwrupr84.0203483.4683984.5722884.6692184.0971185.2413285.3180984.7236585.9125485.9669785.3482586.5856986.6158585.9711487.26056多元线性回归一、基本原理1.多元线性回归模型：Y=0+1Xi+..・+/产其中X1,…,XN是自变量，Y是因变量，0，］...,N是待求的未知参数，£是随机误差项（残差），若记多元线性回归模型可写为矩阵形式：Y=Xp+s通常要求：矩阵X的秩为k+1（保证不出现共线性），且k<N;£为正态分布，E（£）=0和E（££’）二21，其中I为NXN单位矩阵。用最小二乘法原理，令残差平方和最小，得到为P的最佳线性无偏估计量（高斯一马尔可夫定理）。2.2的估计和T检验选取2的估计量：则假如t值的绝对值相当大，就可以在适当选定的置信水平上否定原假设，参数的1-a置信区间可由下式得出：其中t〃为与a%显著水平有关的t分布临界值。a/23.R2和F检验若因变量不具有0平均值，则必须对R2做如下改进：随着模型中增添新的变量，R2的值必定会增大，为了去掉这种增大的干扰，还需要对R2进行修正（校正拟合优度对自由度的依赖关系）：做假设检验：H0：1二••=N=0；H1：1...，N至少有一个NO；使用F统计量做检验，若F值较大，则否定原假设。4.回归诊断（1）残差图分析残差图就是以残差£=y-y为纵坐标，某一个合适的自变量为横坐标的散点图。回归模型中总是假定误差项是独立的正态分布随机变量，且均值为零和方差相等为2.如果模型适合于观察到的数据，那么残差作为误差的无偏估计，应基本反映误差的假设特征。即残差图应该在零点附近对称地密布，越远离零点的地方就疏散（在形象上似有正态趋势），则认为模型与数据拟合得很好。若残差图呈现如图（a）所示的形式，则认为建立的回归模型正确，更进一步再诊断“学生化残差”是否具有正态性：图（b）表明数据有异常点，应处理掉它重新做回归分析（在SAS的REG回归过程步中用来度量异常点影响大小的统计量是COOKD统计量）；图（c）残差随x的增大而增大，图（d）残差随x的增大而先增后减，都属于异方差。此时应该考虑在回归之前对数据y或x进行变换，实现方差稳定后再拟合回归模型。原则上，当误差方差变化不太快时取变换、厅；当误差方差变化较快时取变换logy或lny；当误差方差变化很快时取变换1/y；还有其他变换，如著名的Box-Cox幕变换"-1.力图(e)(f)表示选用回归模型是错误的。共线性回归分析中很容易发生模型中两个或两个以上的自变量高度相关，从而引起最小二乘估计可能很不精确(称为共线性问题)。在实际中最常见的问题是一些重要的自变量很可能由于在假设检验中t值不显著而被不恰当地剔除了。共线性诊断问题就是要找出哪些变量间存在共线性关系。误差的独立性回归分析之前，要检验误差的独立性。若误差项不独立，那么回归模型的许多处理，包括误差项估计、假设检验等都将没有推导依据。由于残差是误差的合理估计，因此检验统计量通常是建立在残差的基础上。检验误差独立性的最常用方法，是对残差的一阶自相关性进行Durbin-Watson检验。H0:误差项是相互独立的；H1:误差项是相关的检验统计量：DW接近于0，表示残差中存在正自相关；如果DW接近于4,表示残差中存在负自相关；如果DW接近于2,表示残差独立性。二、R语言实现还是用函数实现，不同是需要设置更复杂的formula格式：y〜x1+x2只考虑自变量的主效应(y=k1x1+k2x2)，y〜.表示全部自变量的主效应；y〜x1+x2+x1:x2考虑主效应和交互效应(y=kix1+k2x2+k3x1x2)；y〜x1*x2——考虑全部主效应和交互效应的简写(效果同上)；y〜(x1+x2+x3)A2考虑主效应以及至2阶以下的交互效应，相当于x1+x2+x3+x1:x2+x2:x3+x1:x3y〜x1%in%x2x1含于x2，相当于x2+x2:x1y~(x1+x2)A2-x1:x2表示从(x1+x2)A2中去掉x1:x2y~x1+I((x2+x3)A2)使用I()函数，相当于用(x2+x3)A2计算出新变量h，然后y〜x1+hfunction在表达式中使用数学函数，例如log(y)〜x1+x2三、实例例2现有1990〜2009年财政收入的数据revenue.txt：各变量分别表示：y：财政收入(亿元)x1：第一产业国内生产总值(亿元)x2：第二产业国内生产总值(亿元)x3：第三产业国内生产总值(亿元)x4：人口数(万人)x5：社会消费品零售总额(亿元)x6：受灾面积(万公顷)做多元线性回归分析。setwd("E：/办公资料/R语言/R语言学习系列/codes")revenue=read.table("revenue.txt",header=TRUE)Im.reg=lm(y~x1+x2+x3+x4+x5+x6,revenue)summary(lm.reg)Residuals:Min1QMedian3QMax-295.71-173.5226.5990.16370.01Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)6.046e+043.211e+0318.8298.12e-11'加HZ?x1-1.171e-018.638e-02-1.3560.19828x23.427e-023.322e-021.0320.32107x36.182e-014.103e-0215.0671.31e-0944*x4-5.152e-012.930e-02-17.5851.91e-1044*x5-1.104e-012.878e-02-3.8370.00206**x6-1.864e-021.023e-02-1.8230.09143.—Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:234.8on13degreesoffreedomMultipleR-squared:0.9999,AdjustedR-squared:0.9999F-statistic:2.294e+04on6and13DF,p-value:<2.2e-16说明：拟合优度R2=0.9999,效果非常好。但是多元回归时，自变量个数越多，拟合优度必然越好，所以还要看检验效果和回归系数是否显著。结果解释、回归方程、回归预测与前文类似(略)。结合显著性代码可看出：x1和x2不显著，x6只在0.1显著水平下显著，故应考虑剔除x1和x2.R语言中提供了update()函数，用来在原模型的基础上进行修正，还可以对变量进行运算，其基本格式为：update(object,formula.,...,evaluate=TRUE)其中，object为前面拟合好的原模型对象；formula指定模型的格式，原模型不变的部分用“.”表示，只写出需要修正的地方即可，例如update(lm.reg,.〜.+x7)表示添加一个新的变量叩date(lm.reg,sqrt(.)〜.)表示对因变量y开方，再重新拟合回归模型lm.reg2<-update(lm.reg,.~.-x1-x2)#剔除自变量x1,x2summary(lm.reg2)Residuals:Min1QMedian3QMax-325.62-147.5414.07108.28427.42Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)6.339e+042.346e+0327.0203.89e-14***x36.584e-011.548e-0242.523<2e-16***TOC\o"1-5"\h\zx4-5.438e-011.981e-02-27.4453.09e-14***x5-1.392e-011.918e-02-7.2562.80e-06***x6-1.803e-029.788e-03-1.8420.0854.—Signif.codes:0‘心*’0.001‘心’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:233.6on15degreesoffreedomMultipleR-squared:0.9999,AdjustedR-squared:0.9999F-statistic:3.476e+04on4and15DF,p-value:<2.2e-16逐步回归多元线性回归模型中，并不是所有的自变量都与因变量有显著关系，有时有些自变量的作用可以忽略。这就需要考虑怎样从所有可能有关的自变量中挑选出对因变量有显著影响的部分自变量。逐步回归的基本思想是，将变量一个一个地引入或剔出，引入或剔出变量的条件是“偏相关系数”经检验是显著的，同时每引入或剔出一个变量后，对已选入模型的变量要进行逐个检验，将不显著变量剔除或将显著的变量引入，这样保证最后选入的所有自变量都是显著的。逐步回归每一步只有一个变量引入或从当前的回归模型中剔除，当没有回归因子能够引入或剔出模型时，该过程停止。R语言中，用step()函数进行逐步回归，以AIC信息准则作为选入和剔除变量的判别条件。AIC是日本统计学家赤池弘次，在熵概念的基础上建立的：AIC=2(p+1)-2ln(L)其中，p为回归模型的自变量个数，L是似然函数。注：AIC值越小越被优先选入。基本格式：step(object,direction=,steps=,k=2,...)其中，object是线性模型或广义线性模型的返回结果；direction确定逐步回归的方法，默认“both”综合向前向后法，“backward”向后法(先把全部自变量加入模型，若无统计学意义则剔出模型)，“forward”向前法(先将部分自变量加入模型，再逐个添加其它自变量，若有统计学意义则选入模型)；steps表示回归的最大步数，默认1000；k默认=2,输出为AIC值，=log(n)有时输出BIC或SBC值。另外，有时还需要借助使用drop1(object)和add1(object)函数，其中object为逐步回归的返回结果，判断剔除或选入一个自变量，AIC值的变化情况，以筛选选入模型的自变量。lm.step<-step(lm.reg)summary(lm.step)Call:lm(formula=y~x3+x4+x5+x6,data=revenue)Residuals:Min1QMedian3QMax-325.62-147.5414.07108.28427.42Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)6.339e+042.346e+0327.0203.89e-14***x36.584e-011.548e-0242.523<2e-16***TOC\o"1-5"\h\zx4-5.438e-011.981e-02-27.4453.09e-14***x5-1.392e-011.918e-02-7.2562.80e-06***x6-1.803e-029.788e-03-1.8420.0854.—Signif.codes:0‘***’0.001‘心’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:233.6on15degreesoffreedomMultipleR-squared:0.9999,AdjustedR-squared:0.9999F-statistic:3.476e+04on4and15DF,p-value:<2.2e-16

最终得到最优的模型。说明：默认输出每步的结果（略），进行了3步回归，逐步剔除最终得到最优的模型。SingletermdeletionsModel:y~x3+x4+x5+x6DfSumofSqRSSAIC<none>818775222.40x31316.40x41299.12x5128739293692704250.52x611851231003898224.47dropl(lm.step)了自变量x1和x2，AIC值逐步减小，lm.reg3<-lm(y~x3+x4+x5,revenue)summary(lm.reg3)Call:lm(formula=y~Residuals:Min1Q-336.34-186.82Coefficients:Estimate(Intercept)6.284e+04x36.614e-01x4-5.467e-01x5-1.412e-01x3+x4Median1.52+x5,data=revenue)3QMax89.46437.84Std.Error2.494e+031.651e-022.118e-022.053e-02tvalue25.19140.066-25.813-6.877Pr(>|t|)2.66e-14<2e-161.81e-143.72e-06Signif.codes:0‘***’0.001‘心’0.01‘*’0.05‘.’0.1‘Residualstandarderror:250.5on16degreesoffreedomMultipleR-squared:0.9999,AdjustedR-squared:0.9998F-statistic:4.032e+04on3and16DF,p-value:<2.2e-16说明：使用drop1()函数考察分别剔除每个自变量，AIC值变化的情况，可以看出不剔除x6与剔除x6,AIC值只从222.40变大到224.47，相对其它自变量变化很小。所以，可以考虑剔除掉x6，重新做多元线性回归。（四）回归诊断回归分析之后，还需要从残差的随机性、强影响分析、共线性方面进行诊断。一、残差诊断残差y.res<-lm.reg3$residuals#回归模型的残差y.fit<-predict(lm.reg3)#回归模型的预测值plot(y.res~y.fit,main="残差图”)#绘制残差图，以预测值作为横坐标说明：从图形看，残差分布比较均匀，大致满足随机性。shapiro.test(y.res)#残差的正态性检验Shapiro-Wilknormalitytestdata:y.resW=0.94206,p-value=0.2622说明：p值=0.2622>0.05,接受原假设，即残差服从正态分布。标准化残差残差与数据的数量级有关，除以标准误差后得到标准化残差。理想的标准化残差服从N(0,1).rs<-rstandard(lm.reg3)#得到标准化残差plot(rs~y.fit,main="标准残差图”)shapiro.test(rs)#标准化残差的正态性检验Shapiro-Wilknormalitytestdata:rsW=0.97766,p-value=0.9004学生化残差为了回避标准化残差的方差齐性假设，使用学生化残差。rst<-rstudent(lm.reg3)plot(rs~y.fit,main="学生化残差图”)shapiro.test(rst)Shapiro-Wilknormalitytestdata:rstW=0.97463,p-value=0.848⑷残差自相关性的Durbin-Watson检验使用car包中的函数：durbinwatsonTest(model,alternative=c("two.side","positive","negative"))H0：序列不存在自相关性library(car)durbinWatsonTest(lm.reg3)lagAutocorrelationd-wStatisticp-value2.425790.77rho2.425790.77rho!=0Alternativehypothesis:二、强影响分析对参数估计或预测值有异常影响的数据，称为强影响数据。回归模型应当具有一定的稳定性，若个别一两组数据对估计有异常大的影响，剔除后将得到与原来差异很大的回归方程，从而有理由怀疑原回归方程是否真正描述了变量间的客观存在的关系。1.反映这种强影响的统计量有4种及函数：Leveragehatvalues(model)DEFITS——dffits(model)Cook’s距离cooks.distance(model)COVRATIOcovratio(model)另外，influence.measures(model)函数，可以汇总上述4种统计量，判断强影响点。influence.measures(lm.reg3)Influencemeasuresoflm(formula=y~x3+x4+x5,data=revenue):dfb.1_dfb.x3dfb.x4dfb.x5dffitcov.rcook.dhatinfTOC\o"1-5"\h\z0.344152-3.04124-0.4612242.916617-3.409450.8102.14e+000.6347*0.679128-0.09558-0.7071590.3090341.617040.5155.04e-010.3127*-1.7022341.565061.816261-1.982696-3.334521.4262.25e+000.6996*说明：判断出第18,19,20个样本是强影响点。2.Bonferroni离群点检验使用car包中的函数outlierTest(model)library(car)outlierTest(lm.reg3)NoStudentizedresidualswithBonferonnip<0.05Largest|rstudent|:rstudentunadjustedp-valueBonferonnip18-2.5866470.020640.4128注：去掉强影响点，重新做多元线性回归(略)。三、共线性诊断回归分析中很容易发生模型中两个或两个以上的自变量高度相关，从而引起最小二乘估计可能很不精确(称为共线性问题)。在实际中最常见的问题是一些重要的自变量很可能由于在假设检验中t值不显著而被不恰当地剔除了。共线性诊断问题就是要找出哪些变量间存在共线性关系。模型条件数检验使用函数kappa(z,exact=FALSE,…)，其中，z为矩阵XTX，或lm、glm的返回对象；exact设置是否精确计算。一般认为：当K<100时不存在多重共线性；当100WK<1000时存在较强的多重共线性；当KN1000时存在严重的多重共线性。x<-scale(revenue[,3:8])#取出自变量数据，做标准化xx=crossprod(x)#求x’x即矩阵的叉积kappa(xx)[1]6132.142方差膨胀因子(VIF)检验使用car包中的函数vif(model)，该函数还能判断哪些自变量间存在共线性。般认为：当vif<10时不存在多重共线性；当10Wvif<100时，存在较强的多重共线性；当vifN100时存在严重的多重共线性。Im.reg<-lm(y~x1+x2+x3+x4+x5+x6,revenue)vif(lm.reg)x1x2x3x4x5x6196.993779777.7625221014.24830610.484018342.6054901.278766cor(revenue$x2,revenue$x3)#x2和x3的vif值最大，考察二者的相关性[1]0.9977899可见，x2和x3存在严重的共线性，应该考虑剔除其中的一个。岭回归多元线性回归分析中，我们会在众多变量中选择对因变量显著性影响大的那些自变量。但常常会遇到一个问题：在某些情况下，增加或剔除一个自变量后，回归系数变化很大甚至改变符号。主要原因就是变量之间存在多重共线性。岭回归分析是一种专用于共线性数据分析的有偏估计回归方法，实质上是一种改良的最小二乘估计法，通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价，获得回归系数更为符合实际、更可靠的回归方法。基本原理：当自变量间存在多重共线性时，有|X『X周0，考虑加上一个正常数矩阵kI,(k>0),则xtx+kI接近奇异的程度就会比xtx小很多，从而消除了多重共线性。考虑到变量的量纲，应先对数据进行标准