《统计分析软件：使用R与Python》 课件 第7章 R 语言回归分析

7.1一元回归模型7.2多元回归模型7.3回归模型的拟合优度7.4回归模型诊断7.5模型选择7.6模型的预测第7章R语言回归分析R语言回归分析一元回归模型lm函数，R语言公式回归模型诊断残差分析、异常点检测、共线性检测（条件数、VIF、FG检验）模型选择模型的预测多元回归模型系数置信区间、标准化回归模型的拟合优度需要建立数学模型,使得能够根据自变量的数

值预测因变量的大小,或者解释因变量的变化.血压与年龄刹车距离与车速薪金与资历、教育程度、工作岗位认识程度的限制客观事物的复杂人们关心的(因)变量受另外几个(自)变量的关联性(非因果性)

影响,并且存在众多随机因素.无法分析对象内在的因果关系回归分析(RegressionAnalysis)的具体步骤收集一组包含因变量和自变量的数据；选定因变量与自变量之间的模型，利用数据按照最小二乘准则计算模型中的系数；利用统计分析方法对不同的模型进行比较，找出与数据拟合得最好的模型；检验得到的模型是否适合于这组数据；利用模型对因变量作出预测或解释.

7.1

一元回归模型

R语言模型formula

I就表示其中所有的运算符都是普通意义上的算术运算符思考题

例：一元回归分析示例物理学家James．D．Forbes试图通过水的沸点来估计海拔高度，他知道通过气压计测得的大气压可用于得到海拔高度，气压越低，高度越高，他测量了17个地方水的沸点（℉）及大气压数据，并且对数据作了简单的处理，得到了较为明确的数学关系，所提数据如下：X<-matrix(c(194.5,20.79,1.3179,131.79,194.3,20.79,1.3179,131.79,197.9,22.40,1.3502,135.02,198.4,22.67,1.3555,135.55,199.4,23.15,1.3646,136.46,199.9,23.35,1.3683,136.83,200.9,23.89,1.3782,137.82,201.1,23.99,1.3800,138.00,201.4,24.02,1.3806,138.06,201.3,24.01,1.3805,138.05,203.6,25.14,1.4004,140.04,204.6,26.57,1.4244,142.44,209.5,28.49,1.4547,145.47,208.6,27.76,1.4434,144.34,210.7,29.04,1.4630,146.30,211.9,29.88,1.4754,147.54,212.2,30.06,1.4780,147.80),ncol=4,byrow=T,dimnames=list(1:17,c("F","h","log","log100")))forbes=as.data.frame(X)

相关数据如下：首先，画自变量和因变量之间的散点图，进行数据探索分析。从散点图上发现X和Y的排列基本是在一条直线附近，那么我们可以假设X和Y的关系是线性的。R语言实现有三种方式可以实现最小二乘法的简单线性回归，假设数据框为forbeslm(forbes$log100~forbes$F)lm(log100~F,data=forbes)attach(forbes)

;lm(log100~F);#处理完后用detach解除绑定>lm.sol<-lm(log100~F,data=forbes)>summary(lm.sol)通过P值(就是Pr那一列)来查看对应的解释变量x的显著性，通过将p值与0.05进行比较，若该值小于0.05，就可以说该变量与被解释变量存在显著的相关性。MultipleR-squared和AdjustedR-squared这两个值，又被称为”拟合优度“和”修正的拟合优度“，是指回归方程对样本的拟合程度，这里我们可以看到，修正的拟合优度为0.99946，表示拟合程度超过五成，这个值越高越好。F-statistic，也称为F检验，常用于判断方程整体的显著性实验，其p值<2.2e-16，显然小于0.05，方程在P=0.05的水平上是通过显著性检验的。回归方程

可以将得到的直线方程添加到散点图上。提取模型信息lm()

的返回值是一个模型拟合结果对象；技术上就是

"lm"

的一个结果列表类。关于拟合模型的信息可以用能调用对象类

"lm"

的泛型函数显示，提取，图示等等。这包括add1coefeffectskappapredictresidualsaliasdeviancefamilylabelsprintstepanovadrop1formulaplotprojsummary

其中一些常用的泛型函数可以简洁描述如下。anova(object_1,

object_2)比较一个子模型和原模型，并且产生方差分析表。coef(object)提取回归系数(矩阵)。全称：coefficients(object).deviance(object)残差平方和，如有权重可加权。formula(object)提取模型信息。plot(object)产生4个图，显式残差，拟合值和一诊断图。print(object)简要打印一个对象的内容。常常隐式使用。residuals(object)提取残差(矩阵)，有权重时可加权，省略方式：

resid(object)。step(object)通过增加或者减少模型中的项并且保留层次来选择合适的模型。在逐步搜索过程中，aiC(akaike信息规范)值最大的模型将会被返回。summary(object)显示较详细的模型拟合结果。7.2

多元回归模型

多元回归分析示例以R语言MASS包中的cement数据为例cement数据包括美国Portland市某种水泥产生的热量数据，共13组。水泥凝固时放出的热量y（响应变量）和水泥混合物中的四种化学成分x1,x2,x3,x4.用全部的变量来解释热量y，对应的R代码和结果:可以发现，F检验的p值很小，说明回归方程是显著的，但在各个系数的t检验中，p值均大于0.05，说明每个系数都是不显著的，回归方程可能包括不重要的变量或存在共线性等问题，这样的模型不能直接应用。回归系数的置信区间

其中object参数为拟合的回归模型，parm参数表示计算置信区间的变量名称标准化的回归系数标准化回归系数，它是在对自变量和因变量同时进行标准化处理后所得到的回归系数，数据经过标准化处理后消除了量纲、数量级等差异的影响，使得不同变量之间具有可比性，因此可以用标准化回归系数来比较不同自变量对因变量的作用大小。我们用scale()函数给出标准化回归系数的回归模型代码R标准化scale函数

set.seed(1);x<-runif(10)#Manuallyscaling(x-mean(x))/sd(x)ss=scale(x)identical(attr(ss,"scaled:center"),mean(x))identical(attr(ss,"scaled:scale"),sd(x))[1]-0.8717643-0.52873940.11708951.1960620-1.07712101.1644732[1]TRUE[1]TRUE未标准化回归系数体现的是自变量变化对因变量的绝对作用大小标准化回归系数反映的是不同自变量对因变量的相对作用大小，可以显示出不同自变量对因变量影响的重要性。如果用标准化回归系数构建方程，得到的结论是有偏的，因为此时自变量和因变量的数据都发生了转化，成为了标准化数据，因此标准化回归系数不能直接用于回归方程预测。7.3

回归模型的拟合优度

7.4

回归模型诊断

拟合好一个回归模型并不是直接就能够使用了，而是先要对回归分析中的假设（高斯-马尔可夫条件）以及数据进行检验与分析。检验回归分析中的假设是否合理。对数据的诊断7.4.1残差的分类

R语言中可利用rstandard()调用回归模型的标准化残差。

7.4.2线性假设诊断判断线性假设，可以对lm()函数返回的模型对象使用plot()函数。该函数可生成评价模型拟合情况的多个图形，默认返回四个图形。可使用plot()函数中的残差-拟合值（Residualvsfitted）图判断回归模型中线性关系是否成立。plot(lm.sol,which=1)残差与预测值的变化之间不存在明显的规律，可以认为基本满足线性关系。ResidualsvsFitted7.4.3

残差分析和异常点检测高斯-马尔可夫条件中涉及随机误差项的假设较多，而误差是无法被观测的，因此无法直接利用误差检验模型假设是否成立。我们可以通过判断残差是否符合误差项的某些性质，达到检验模型是否符合相应假设的目的。在R中，我们可以使用plot()函数中的正态Q-Q图（NormalQQ-plot）判断Forbes数据判断残差是否服从正态分布。plot(lm.sol,which=2)除开12号样本点外，其他数据点分布趋于一条直线,说明残差近似服从正态分布。链接：不同情形下的NormalQQPlots在R中，我们还可以使用plot()函数中的尺度-位置（Scale-Location）图显示的残差分布情况，判断同方差性：若满足不变方差假设，那么在位置尺度图中，水平线周围的点应该随机分布。plot(lm.sol,which=3)大部分数据点呈现出随机的分布，但部分数据（如12号点）的标准化残差平方根过大，说明原数据集中不满足同方差的假设。

图中12号样本点的标准化残差较大，属于离群点，但该点杠杆值和Cook距离不大，故不是高杠杆点。我们可以用plot()函数中的Cook距离图进一步观察各点的Cook距离plot(lm.sol,which=4)一般来说，小样本的Cook距离大于1，大样本的Cook距离大于4/n算作高杠杆点。对于离群点的处理，可以去掉12号样本点再拟合可以看出，去掉12号样本点后回归方程系数没有太大的变化，但是系数的标准差和残差标准差有很大的变化，平均减少约3倍左右，R平方的值也有提高。7.4.4多重共线性检测多重共线性是指在多元回归模型中，自变量（解释变量）之间存在一定的线性相关关系而使模型估计失真或估计精度不高。在Gauss-Markov条件满足的条件下，多重共线性不改变估计量的无偏性，但各共线变量参数的最小二乘法估计值的精度很低，容易使结果变得不显著，进而无法正确判断各自变量对因变量的影响。R中常见验证变量之间是否存在共线性的方法：条件数法、方差膨胀因子(VIF)、Farrar-Glauber方法

多重共线性问题#ExactCollinearitygen_exact_collin_data=function(num_samples=100){x1=rnorm(n=num_samples,mean=80,sd=10)x2=rnorm(n=num_samples,mean=70,sd=5)x3=2*x1+4*x2+3y=3+x1+x2+rnorm(n=num_samples,mean=0,sd=1)data.frame(y,x1,x2,x3)}set.seed(1234)exact_collin_data=gen_exact_collin_data()当我们尝试使用所有预测变量在R中拟合回归模型时会发生什么？>head(exact_collin_data)yx1x2x31143.487267.9293472.07262427.14922154.097582.7742967.62641439.05423164.359990.8444170.32997466.00874127.731456.5430267.48761386.03655153.472984.2912565.87001435.06256159.656085.0605670.83495456.4609>exact_collin_fit=lm(y~x1+x2+x3,data=exact_collin_data)>summary(exact_collin_fit)Call:lm(formula=y~x1+x2+x3,data=exact_collin_data)Residuals:Min1QMedian3QMax-3.1933-0.67340.10600.58232.7650Coefficients:(1notdefinedbecauseofsingularities)EstimateStd.ErrortvaluePr(>|t|)(Intercept)1.2829981.5371400.8350.406x11.0080510.009633104.649<2e-16***x21.0176640.01874754.284<2e-16***x3NANANANA---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:0.9624on97degreesoffreedomMultipleR-squared:0.9929, AdjustedR-squared:0.9928F-statistic:6809on2and97DF,p-value:<2.2e-16

>X=cbind(1,as.matrix(exact_collin_data[,-1]))#inverse>solve(t(X)%*%X)Errorinsolve.default(t(X)%*%X):systemiscomputationallysingular:reciprocalconditionnumber=1.44535e-18fit1=lm(y~x1+x2,data=exact_collin_data)fit2=lm(y~x1+x3,data=exact_collin_data)fit3=lm(y~x2+x3,data=exact_collin_data)all.equal(fitted(fit1),fitted(fit2))##[1]TRUE我们看到三个模型的拟合值完全相同。WHY？这是因为x3包含来自x1和x2的所有信息的结果。只要模型中包含x1或x2之一，x3就可以用于从未包含的变量中恢复信息。虽然拟合值都相同，但模型的估计系数却大不相同。x2的符号在fit2、fit3两个模型中甚至是相反的！>coef(fit1)(Intercept)x1x21.2829981.0080511.017664>coef(fit2)(Intercept)x1x30.51975030.49921950.2544160>coef(fit3)(Intercept)x2x3-0.2290789-0.99843900.5040257

在(0,10)之间,没有多重共线性.[10,100)有较强多重共线性.大于100，有严重的多重共线性.例：考虑一个有6个回归自变量的线性回归问题，原始数据列在表中，这里共有12组数据，除第一组外，自变量x1,x2,…,x6的其余11组数据满足线性关系：x1+x2+x3+x4=10，试用求矩阵条件数的方法，分析出自变量间存在多重共线性。例：R语言多重共线性collinear<-data.frame(y=c(10.006,9.737,15.087,8.422,8.625,16.289,5.958,9.313,12.96,5.541,8.756,10.937),x1=rep(c(8,0,2,0),c(3,3,3,3)),x2=rep(c(1,0,7,0),c(3,3,3,3)),x3=rep(c(1,9,0),c(3,3,6)),x4=rep(c(1,0,1,10),c(1,2,6,3)),x5=c(0.541,0.13,2.116,-2.397,-0.046,0.365,1.996,0.228,1.38,-0.798,0.257,0.44),x6=c(-0.099,0.07,0.115,0.252,0.017,1.504,-0.865,-0.055,0.502,-0.399,0.101,0.432))先计算出对应矩阵的特征值：可以利用定义计算条件数：e$values[1]/e$values[6][1]2195.908#条件数大于100有严重的多重共线性。

e$vectors[,6][1]0.4476797190.4211402800.5416891240.5733718720.0060521270.002166594也可直接用kappa()函数计算条件数:kappa(xx,exact=TRUE)##[1]2196方差膨胀因子(VIF)

Farrar-Glauber方法

用mctest包的omcdiag函数进行检验发现标准化行列式的值为0.001，这是非常小的。卡方检验统计量的计算值为56.951，它非常显著，因此暗示模型中存在多重共线性。这促使我们对多重共线性的位置进行下一步Farrar-Glauber检验（F检验）mctest::imcdiag(lmcoliner)VIF、TOL和Wi列分别为方差膨胀因子、容差和F检验的结果。F检验表明:变量“x1”或“x2”或“x3”或“x4”可能是多重共线性的根本原因。正如预期的那样，“x1”和“x2,x3,x4”之间高度偏相关，具有统计显著性。Farrar-Glauber检验指出这几个变量是所有多重共线性问题的根本原因。我们去掉“x2,x3,x4”，只保留”x1”,对应的VIF指标如下：VIF的值较小，说明剩余变量中不存在共线性问题

模型诊断-car包实用函数qqPlot()#分位数比较图durbinWastonTest()#误差的独立性Durbin-Waston检验crPlots()#成分和残差图,偏残差图，用于检测因变量与每一个自变量之间的线性关系是否显著outlierTest()#Bonferroni离群点检验avPlots()#添加的变量图形，对于每个预测变量Xk，绘制Xk在其他k-1个预测变量上回归的残差值influencePlot()#回归影响图,将离群点、杠杆值和强影响点的信息整合到一幅图形中ncvTest()#同方差性检测vif()#检验多重共线性7.5模型选择当预测模型中有许多预测变量时，模型选择方法允许自动选择预测变量的最佳组合以构建最佳的回归模型。一种很自然的模型选择策略是穷举所有可能的预测变量组合，然后从中选择最佳的回归模型。这种方法称为最佳子集回归或全子集回归。另外一个是逐步回归，即按顺序比较包含不同个数预测变量组合的线性回归模型。7.5.1最佳子集回归最佳子集回归可用leaps包中的regsubsets()函数实现参数nvmax的值表示模型中包含的最大预测变量的个数。最好的2变量模型只包含x1和x2（y~x1+x2）最好的3变量模型则是（y~x1+x2+x4）这些最佳模型哪个最好调整后的R2告诉我们最好的模型是第3个模型即具有3个预测变量的模型(y~x1+x2+x4)。但是，使用BIC和Cp标准，则应该选择具有2个变量的模型（y~x1+x2）。可根据实际情况，选择“最佳”的模型。Allmodelsarewrong,butsomeareuseful7.5.2逐步回归(stepwiseregression)逐步回归中，模型会一次添加或者删除一个变量，直到达到某个判停准则为止。向前逐步回归-forwardstepwise-每次添加一个预测变量到模型中，直到添加变量不会使模型有所改进为止。向后逐步回归-backwardstepwise-从模型包含所有预测变量开始，一次删除一个变量直到降低模型质量为止。向前向后逐步回归（stepwisestepwise，通常称作逐步回归），结合了向前逐步回归和向后逐步回归的方法，变量每次进入一个，但是每一步中，变量都会被重新评价，对模型没有贡献的变量将会被删除，预测变量可能会被添加、删除好几次，直到获得最优模型为止。逐步回归

全子集回归可用leaps包的regsubsets()函数实现。R自带的step函数实现了最基本的逐步回归模型。MASS包中的stepAIC()函数可以实现逐步回归模型（向前、向后和向前向后），依据的是精确AIC准则。Call:lm(formula=y~x1+x2+x3+x4,data=cement)Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)62.405470.07100.8910.3991

x11.55110.74482.0830.0708.x20.51020.72380.7050.5009x30.10190.75470.1350.8959x4-0.14410.7091-0.2030.8441---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘

’1Residualstandarderror:2.446on8degreesoffreedomMultipleR-squared:0.9824,AdjustedR-squared:0.9736F-statistic:111.5on4and8DF,p-value:4.756e-07

#解释变量不显著考虑所有可用预测变量的模型作为初始模型Step处理和其他软件不同，R软件根据AIC准则，自动选择合适的解释变量.>summary(lm.step)Call:lm(formula=y~x1+x2+x4,data=cement)Residuals:Min1QMedian3QMax-3.0919-1.80160.25621.28183.8982Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)71.648314.14245.0660.000675***x11.45190.117012.4105.78e-07***x20.41610.18562.2420.051687.x4-0.23650.1733-1.3650.205395---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:2.309on9degreesoffreedomMultipleR-squared:0.9823, AdjustedR-squared:0.9764F-statistic:166.8on3and9DF,p-value:3.323e-08drop1(lm.step)与改进：>drop1(lm.step)SingletermdeletionsModel:y~x1+x2+x4DfSumofSqRSSAIC<none>47.9724.974x11820.91868.8860.629x2126.7974.7628.742x419.9357.9025.420删除变量会改变根据AIC最小准则得到的step结果，根据drop1分析，删除x4使AIC值上升最小，为尽量遵守AIC准则的同时，提高解释变量的显著性，删除x4是“无奈之举”。drop1(lm.step)与改进：newsol=lm(y~x1+x2,data=cement)summary(newsol)Call:lm(formula=y~x1+x2,data=cement)Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)52.577352.2861723.005.46e-10***x11.468310.1213012.112.69e-07***x20.662250.0458514.445.03e-08***Residualstandarderror:2.406on10degreesoffreedomMultipleR-squared:0.9787,AdjustedR-squared:0.9744F-statistic:229.5on2and10DF,p-value:4.407e-09

比较：step()处理后的残差和可决系数：Residualstandarderror:2.309AdjustedR-squared:0.9764在MASS包中，stepAIC()函数实现逐步回归模型stepAIC(object,scope,scale=0,direction=c("both",

"backward","forward"),trace=1,keep=NULL,steps=1000,use.start=FALSE,k=2,...)“backward”表示从给定模型中依次删除预测变量，它产生一系列复杂性递减的模型，直到获得最佳模型；“forward”表示使用数据集中的所有可用变量，将预测变量按顺序添加到给定模型中，它产生一系列复杂性增加的模型，直到获得最佳模型；“both”（默认）表示向前-向后搜索，在每一步，决定是包含还是排除预测变量，与之前的模式不同，先前被排除/包括的预测变量可以稍后被包括/排除。stepAIC()函数默认考虑AIC作为度量准则

第一步返回的结果是”Step：AIC=29.77”，第二步移除了x3，因为它给出了最低的BIC,此时BIC变为“Step:AIC=27.23”。最后一步移除了x4,BIC值进一步减少，改进为: