R语言实战-topic6回归分析

1、 Topic6 回归 什么是回归 回归其实是一个广义的概念，通指那些用一个或多个预测变量（也称自变量或解释变 量）来预测响应变量（也称因变量、效标变量或结果变量）的方法。通常，回归分析 可以用来挑选与响应变量相关的解释变量，可以描述两者的关系，也可以生成一个等式，通过解释变量来预测响应变量。 表8-1回归分析的各种变体 回归类型 用途 简单线性 当项式 多层 用一个量化的解释变量预测一个量化的响应变量 用一个乐化的解群变展侦测一个星化的响应变晶.模型的关系是月阶军项式 用拥有等级结构的数据预测一个响应变量I例如学校中教室里的学生L也被称为分层模型.被套模型或混合模型 名元线件篓变后Logist

2、ic泊松 Cox比例风险 时间序列非线性非搴数稳健 用两个或方个量化的解择变量预测一个量化的响应变量 用一个或多个解释变景预测多个响应变量 用一个成名个解释变量前测一个类别型响应变后 用一个或多个解释变量预测一个代表频数的响应变量 用一个或多个解释变后预测一个事件（死亡、失败或用病夏发）发电的时间 对溟差项相关的时间序列数据建模 用一个或2个最化的解释变成预测一个量化的响应变盘,不过模型地等线性的 用一个或多个最化的解样变最预测一个量化的响应变型.模型的形式源H数据形式,不事先设定用个或多个量化的解释变量预测一个量化的响应垩量，能抵御强影响点的干扰 我们的重点是普通最小二乘（OLS）回归法，包

3、括简单线性回归、多项式回归和多元线性回归。 OLS回归是现今最常见的统计分析方法，其他回归模型（Logistic回归和泊松回归）将在第13章介绍 OLS回归的适用情境 OLS回归是通过预测变量的加权和来预测量化的因变量，其中权重是通过数据估计而 得的参数 示例 一名工程师想找出跟桥梁退化有关的最重要的因素，比如使用年限、交通流量、桥梁 设计、建造材料和建造方法、建造质量以及天气情况，并确定它们之间的数学关系。 他从一个有代表性的桥梁样本中收集了这些变量的相关数据，然后使用OLS回归对数 据进行建模。解决以下几个方面的问题： 在众多变量中判断哪些对预测桥梁退化是有用的，得到它们的相对重要性，从而

4、关 注重要的变量。 根据回归所得的等式预测新的桥梁的退化情况 （预测变量的值已知， 但是桥梁退化程度未知） ，找出那些可能会有麻烦的桥梁。 利用对异常桥梁的分析，获得一些意外的信息。比如他发现某些桥梁的退化速度比预测的更快或更慢，那么研究这些“离群点”可能会有重大的发现，能够帮助理解桥梁退化的机制。 OLS回归模型的形式： 彳=A+力西+自乙，=1打 n为观测的数目，k为预测变量的数目，Y?第i次观测对应的因变量的预测值，X/第i次观测对应的第j个预测变量值，用截距项，片预测变量j的回归系数 目标是使残差平方和最小： -寓”匹十十A*J谪 用lm（）拟合回归模型 myfit-lm（formul

5、a,data） 表达式（formula）形式如下： Y。XI+X2+Xk 表8-3对拟合线性模型非常有用的其他函数 函数用途 surranary( 展示拟台模型的详细结果 coefficients() 列出拟合模型的模型零数（截环项和斜率） contint() 提供模型参数的置信区间（默认95%） fitted() 列出拟合模型的捌漫值 rfisidviaLs) 列出拟合模鳖的残差值 anova) 生成一个拟合校地的疔差分析表，或者比较两个或更多拟合模型的方差分析去 vcov() 列出模型参数的办力差矩阵 AICO 输出赤池信息统订量 plot() 生成评价拟合模型的诊斯图 prodiCT(

6、用拟合模型对新的数据柒顼测响应变量值 当回归模型包含一个因变量和一个自变量时，我们称为简单线性回归。当只有一个预测变量，但同时包含变量的哥（比如，X,X2,X3）时，我们称为多项式回归。当有不止一个预测变量时，则称为多元线性回归 简单线性回归 通过身高来预测体重(women数据集) fit|t|) (Intercept)-87.516675.93694-14.741.71e-09* height3.450000.0911437.851.09e-14* Signif.codes:0*0.001*0.01*0.05.0.11 Residualstandarderror:1.525on13degre

7、esoffreedomMultipleR-squared:0.991,AdjustedR-squared:0.9903 F-statistic:1433on1and13DF,p-value:1.091e-14 women$weight 1115117120123126129132135139142146150154159164 fitted(fit) 12345678 112.5833116.0333119.4833122.9333126.3833129.8333133.2833136.7333 9101112131415 140.1833143.6333147.0833150.5333153

8、.9833157.4333160.8833 residuals(fit) 123456 2.416666670.966666670.516666670.06666667-0.38333333-0.83333333 789101112 -1.28333333-1.73333333-1.18333333-1.63333333-1.08333333-0.53333333 131415 0.016666671.566666673.11666667 plot(women$height,women$weight, main=WomenAge30-39, xlab=Height(ininches),ylab

9、=Weight(inpounds)abline(fit) womenAge30-39 Migrilininch?l 通过输出结果，可以得到预测等式： 而福二-87.52+3都川例亚 多项式回归 通 过 添 加 一 个 二 次 项 ( 即X2)来 提 高 回 归 的 预 测 精 度fit2|t|) (Intercept)261.8781825.1967710.3932.36e-07* height-7.348320.77769-9.4496.58e-07* I(heightA2)0.083060.0059813.8919.32e-09* Signif.codes:0*0.001*0.01*0.0

10、5.0.11 Residualstandarderror:0.3841on12degreesoffreedom MultipleR-squared:0.9995,AdjustedR-squared:0.9994 F-statistic:1.139e+04on2and12DF,p-value:2.2e-16 plot(women$height,women$weight, main=WomenAge30-39, xlab=Height(ininches),ylab=Weight(inlbs)lines(women$height,fitted(fit2) 新的预测等式为： Weight=26L8-7

11、35,Height+0.083丈Hei鲍线性模型与非线性模型 多项式等式仍可认为是线性回归模型，因为等式仍是预测变量的加权和形式(本例中是身高和身高的平方)。即使这样的模型： =A*iog%+乂疝羽 仍可认为是线性模型(参数项是线性的)，能用这样的表达式进行拟合： YJlog(XJ+sin(XJ 相反，下面的例子才能算是真正的非线性模型： 这种非线性模型可用nls()函数进行拟合 一般来说，n次多项式生成一个nT个弯曲的曲线。拟合三次多项式，可用： fit3-Lm(weigliCheight+I(heightA2)+1(helghtA3),data=women)虽然更高次的多项式也可用，但我发

12、现使用比三次更高的项几乎没有必要。 用scatterplot绘制二元关系图 library(car) scatterplot(weightheight,data=women, spread=FALSE,smoother.args=list(lty=2),pch=19, main=WomenAge30-39”,xlab=Height(inches),ylab=Weight(lbs.) WorrtanAge304S 多元线性回归 当预测变量不止一个时，简单线性回归就变成了多元线性回归 多元回归分析中，第一步最好检查一下变量间的相关性。cor()函数提供了二变量之间 的相关系数，car包中scatt

13、erplotMatrix()函数则会生成散点图矩阵 states-as.data.frame(state.x77,c(Murder,Population, Illiteracy,Income,Frost)cor(states) MurderPopulationIlliteracyIncomeFrost Murder1.00000000.34364280.7029752-0.2300776-0.5388834 Population0.34364281.00000000.10762240.2082276-0.3321525 Illiteracy0.70297520.10762241.0000000

14、-0.4370752-0.6719470 Income-0.23007760.2082276-0.43707521.00000000.2262822 Frost-0.5388834-0.3321525-0.67194700.22628221.0000000 library(car) scatterplotMatrix(states,spread=FALSE,smoother.args=list(lty=2),main=ScatterPlotMatrix) ScatterPlotMaINx wa15aMi scatterplotMatrix()函数默认在非对角线区域绘制变量间的散点图，并添加平滑

15、和线性拟合曲线。对角线区域绘制每个变量的密度图和轴须图。 使用lm()函数拟合多元线性回归模型 states-as.data.frame(state.x771c(Murder,Population, Illiteracy,Income,Frost)fit|t|)(Intercept)1.235e+003.866e+000.3190.7510Population2.237e-049.052e-052.4710.0173* Illiteracy4.143e+008.744e-014.7382.19e-05* Income6.442e-056.837e-040.0940.9253 Frost5.81

16、3e-041.005e-020.0580.9541 Signif.codes:0*0.001*0.01*0.05.0.11 Residualstandarderror:2.535on45degreesoffreedom MultipleR-squared:0.567,AdjustedR-squared:0.5285 F-statistic:14.73on4and45DF,p-value:9.133e-08 有交互项的多元线性回归 以mtcars数据框中的汽车数据为例，若你 对汽车重量和马力感兴趣，可以把它们作为预测变量，并包含交互项来拟合回归模型 fit|t|)(Intercept)49.80

17、8423.6051613.8165.01e-14*hp-0.120100.02470-4.8634.04e-05* wt-8.216621.26971-6.4715.20e-07* hp:wt0.027850.007423.7530.000811* Signif.codes:0*0.001*0.01*0.05.0.11 Residualstandarderror:2.153on28degreesoffreedom MultipleR-squared:0.8848,AdjustedR-squared:0.8724 F-statistic:71.66on3and28DF,p-value:2.981

18、e-13 通过effects包中的effect()函数，你可以用图形展示交互项的结果。格式为：library(effects) plot(effect(hp:wt,fit,list(wt=c(2.2,3.2,4.2),multiline=TRUE) nWcrplcK 勇MlSW酎211 从图中可以很清晰地看出， 随着车重的增加， 马力与每加仑汽油行驶英里数的关系减弱了。 当wt=4.2时，直线几乎是水平的，表明随着hp的增加，mpg不会发生改变。 回归诊断 通过confint()函数的输出来看看8.2.4节中states多元回归的问题 states-as.data.frame(state.x7

19、7,c(Murder,Population, Illiteracy,Income,Frost) fit-lm(MurderPopulation+Illiteracy+Income+Frost,data=states)confint(fit) 2.5%97.5% (Intercept)-6.552191e+009.0213182149 Population4.136397e-050.0004059867 Illiteracy2.381799e+005.9038743192 Income-1.312611e-030.0014414600 Frost-1.966781e-020.0208304170

20、 结果表明，文盲率改变1%,谋杀率就在95%的置信区间2.38,5.90中变化。另外，因为Frost。的置信区间包含0,所以可以得出结论：当其他变量不变时，温度的改变与谋杀率无关。 检验回归分析中统计假设的方法简单线性回归 fit-lm(weightheight,data=women)par(mfrow=c(2,2) plot(fit) 二次拟合诊断图 newfit-lm(weightheight+I(heightA2),data=women)par(opar) par(mfrow=c(2,2) plot(newfit)FittedvaivRS 1201341MMQ网 ResiduatevsF

21、led EEIrn苜y PiBEiylQuartiles rJir1.J, 5CAteLocafi rl一/npr9FU7nBFUrn.v-J 000(ILOE0.10DJSQ3):期 R檄mzisvsleversje HCMIY刖fgLev*w 最后，我们再应用这个基本的方法，来看看states的多元回归问题。 opar-par(no.readonly=TRUE) states-as.data.frame(state.x77,c(Murder,Population, Illiteracy,Income,Frost)fit-lm(MurderPopulation+Illiteracy+Inco

22、me+Frost,data=states)par(mfrow=c(2,2)plot(fit)par(opar) 线性模型假设的综合验证 gvlma能对线性模型假设进行综合验证，同时还能做偏斜度、峰度和异方差性的评价library(gvlma) gvmodel2#problem? PopulationIlliteracyIncomeFrost1.2452822.1658481.3458222.082547 sqrt(vif(fit)2#problem? PopulationIlliteracyIncomeFrostFALSEFALSEFALSEFALSE 离群点检测 outlierTest(fi

23、t) rstudentunadjustedp-valueBonferonnip Nevada3.5429290.000950880.047544 高杠杆值点 高杠杆值观测点，即与其他预测变量有关的离群点。它们是由许多异常的预测变量值组合起来的，与响应变量值没有关系。 高杠杆值的观测点可通过帽子统计量(hatstatistic)判断。对于一个给定的数据集，帽子均值为p/n,其中p是模型估计的参数数目(包含截距项)，n是样本量。一般来说,若观测点的帽子值大于帽子均值的2或3倍，就可以认定为高杠杆值点 hat.plot-function(fit) p-length(coefficients(fit)

24、 n-length(fitted(fit) plot(hatvalues(fit),main=IndexPlotofHatValues) abline(h=c(2,3)*p/n,col=red,lty=2) identify(1:n,hatvalues(fit),names(hatvalues(fit) ) hat.plot(fit) 强影响点 强影响点，即对模型参数估计值影响有些比例失衡的点。例如，若移除模型的一个观测点时模型会发生巨大的改变。 有两种方法可以检测强影响点：Cook距离，或称D统计量，以及变量添加图(added variableplot)o一般来说，CookSD值大于4/(n

25、*T),则表明它是强影响点，其中n为样本量大小，k是预测变量数目。 cutoff.)11h,i,.11.I.llIII iirir VW3CJM 口(HEE 卬 EFAwderlitFaqr*lnE_ane-*PiCBl| 利用car包中的influencePlot()函数，还可以将离群点、杠杆值和强影响点的信息整合到一幅图形中： influencePlot(fit,id.method=identify,main=InfluencePlot, sub=CirclesizeisproportialtoCooksDistance) inniufrnceRui Ciltss3曲片，幅pEHioCgM

26、OslRrte 改进措施 “如果通过回归诊断发现了问题，那么能做些什么呢？”有四种方法可以处理违背回归假设的问题： 1 .删除观测点 2 .变量变换 3 .增删变量 4 .使用其他回归方法 选择“最佳”的回归模型 模型比较 用基础安装中的anova()函数可以比较两个嵌套模型的拟合优度。所谓嵌套模型，即它 的一些项完全包含在另一个模型中。在states的多元回归模型中，我们发现Income和 Frost的回归系数不显著， 此时你可以检验不含这两个变量的模型与包含这两项的模型预测效果是否一样好。 用anova()函数比较 states-as.data.frame(state.x77,c(Murder,Population, Illiteracy,Income,Frost) fit2-lm(MurderPopulation+Illiteracy,data=states) fit1F) 47289.25 45289.1720.0785050.00610.9939 用AIC来比较模型 fitl-lm(MurderPopulation+Illiteracy+Income+Frost,data=states) fit2-lm(MurderPopulation+Illiteracy,data=states) AIC(fit1,fit2) dfAIC