R语言学习系列32回归分析.docx
- 文档编号:24581525
- 上传时间:2023-05-29
- 格式:DOCX
- 页数:28
- 大小:427.98KB
R语言学习系列32回归分析.docx
《R语言学习系列32回归分析.docx》由会员分享,可在线阅读,更多相关《R语言学习系列32回归分析.docx(28页珍藏版)》请在冰豆网上搜索。
R语言学习系列32回归分析
27.回归分析
回归分析是研究一个或多个变量(因变量)与另一些变量(自变量)之间关系的统计方法。
主要思想是用最小二乘法原理拟合因变量与自变量间的最佳回归模型(得到确定的表达式关系)。
其作用是对因变量做解释、控制、或预测。
回归与拟合的区别:
拟合侧重于调整曲线的参数,使得与数据相符;而回归重在研究两个变量或多个变量之间的关系。
它可以用拟合的手法来研究两个变量的关系,以及出现的误差。
回归分析的步骤:
(1)获取自变量和因变量的观测值;
(2)绘制散点图,并对异常数据做修正;
(3)写出带未知参数的回归方程;
(4)确定回归方程中参数值;
(5)假设检验,判断回归方程的拟合优度;
(6)进行解释、控制、或预测。
(一)一元线性回归
一、原理概述
1.一元线性回归模型:
Y=𝛽0+𝛽1X+ε
其中X是自变量,Y是因变量,𝛽0,𝛽1是待求的未知参数,𝛽0也称为截距;ε是随机误差项,也称为残差,通常要求ε满足:
①ε的均值为0;
②ε的方差为𝜎2;
③协方差COV(εi,εj)=0,当i≠j时。
即对所有的i≠j,εi与εj互不相关。
用最小二乘法原理,得到最佳拟合效果的
值:
,
2.模型检验
(1)拟合优度检验
计算R2,反映了自变量所能解释的方差占总方差的百分比,值越大说明模型拟合效果越好。
通常可以认为当R2大于0.9时,所得到的回归直线拟合得较好,而当R2小于0.5时,所得到的回归直线很难说明变量之间的依赖关系。
(2)回归方程参数的检验
回归方程反应了因变量Y随自变量X变化而变化的规律,若𝛽1=0,则Y不随X变化,此时回归方程无意义。
所以,要做如下假设检验:
H0:
𝛽1=0,H1:
𝛽1≠0;
①F检验
若𝛽1=0为真,则回归平方和RSS与残差平方和ESS/(N-2)都是𝜎2的无偏估计,因而采用F统计量:
来检验原假设β1=0是否为真。
②T检验
对H0:
𝛽1=0的T检验与F检验是等价的(t2=F)。
3.用回归方程做预测
得到回归方程
后,预测X=x0处的Y值
.
的预测区间为:
其中tα/2的自由度为N-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))
datas
ageheight
11876.1
21977.0
32078.1
42178.2
52278.8
62379.7
72479.9
82581.1
92681.2
102781.8
112882.8
122983.5
plot(datas)#绘制散点图
res.reg<-lm(height~age,datas)#做一元线性回归
summary(res.reg)#输出模型的汇总结果
Residuals:
Min1QMedian3QMax
-0.27238-0.24248-0.027620.160140.47238
Coefficients:
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)64.92830.5084127.71<2e-16***
age0.63500.021429.664.43e-11***
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
0.256on10degreesoffreedom
MultipleR-squared:
0.9888,AdjustedR-squared:
0.9876
F-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)age
64.9283220.634965
confint(res.reg,parm="age",level=0.95)#输出参数age的置信区间,若不指定parm将返回所有参数的置信区间
2.5%97.5%
age0.58727220.6826578
fitted(res.reg)#输出回归模型的预测值
123456789101112
76.3576976.9926677.6276278.2625978.8975579.5325280.1674880.8024581.4374182.0723882.7073483.34231
anova(res.reg)#输出模型的方差分析表
Response:
height
DfSumSqMeanSqFvaluePr(>F)
age157.65557.655879.994.428e-11***
Residuals100.6550.066
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
vcov(res.reg)#输出模型的协方差矩阵
(Intercept)age
.010*******
age-0.010766860.0004581642
residuals(res.reg)#输出模型的残差
123456789101112
-0.2576923080.0073426570.472377622-0.062587413-0.0975524480.167482517-0.2674825170.297552448-0.237412587-0.2723776220.0926573430.157692308
AIC(res.reg)#输出模型的AIC值
[1]5.161407
BIC(res.reg)#输出模型的BIC值
[1]6.616127
logLik(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.res
fitlwrupr
184.0203483.4683984.57228
284.6692184.0971185.24132
385.3180984.7236585.91254
485.9669785.3482586.58569
586.6158585.9711487.26056
(二)多元线性回归
一、基本原理
1.多元线性回归模型:
Y=𝛽0+𝛽1X1+…+𝛽NXN+ε
其中X1,…,XN是自变量,Y是因变量,𝛽0,𝛽1…,𝛽N是待求的未知参数,ε是随机误差项(残差),若记
多元线性回归模型可写为矩阵形式:
Y=Xβ+ε
通常要求:
矩阵X的秩为k+1(保证不出现共线性),且k ,其中I为N×N单位矩阵。 用最小二乘法原理,令残差平方和 最小,得到 为β的最佳线性无偏估计量(高斯-马尔可夫定理)。 2.𝜎2的估计和T检验 选取𝜎2的估计量: 则 假如t值的绝对值相当大,就可以在适当选定的置信水平上否定原假设,参数的1-α置信区间可由下式得出: 其中tα/2为与α%显著水平有关的t分布临界值。 3.R2和F检验 若因变量不具有0平均值,则必须对R2做如下改进: 随着模型中增添新的变量,R2的值必定会增大,为了去掉这种增大的干扰,还需要对R2进行修正(校正拟合优度对自由度的依赖关系): 做假设检验: H0: 𝛽1=…=𝛽N=0;H1: 𝛽1…,𝛽N至少有一个≠0; 使用F统计量做检验, 若F值较大,则否定原假设。 4.回归诊断 (1)残差图分析 残差图就是以残差 为纵坐标,某一个合适的自变量为横坐标的散点图。 回归模型中总是假定误差项是独立的正态分布随机变量,且均值为零和方差相等为𝜎2.如果模型适合于观察到的数据,那么残差作为误差的无偏估计,应基本反映误差的假设特征。 即残差图应该在零点附近对称地密布,越远离零点的地方就疏散(在形象上似有正态趋势),则认为模型与数据拟合得很好。 若残差图呈现如图(a)所示的形式,则认为建立的回归模型正确,更进一步再诊断“学生化残差”是否具有正态性: 图(b)表明数据有异常点,应处理掉它重新做回归分析(在SAS的REG回归过程步中用来度量异常点影响大小的统计量是COOKD统计量); 图(c)残差随x的增大而增大,图(d)残差随x的增大而先增后减,都属于异方差。 此时应该考虑在回归之前对数据y或x进行变换,实现方差稳定后再拟合回归模型。 原则上,当误差方差变化不太快时取变换 ;当误差方差变化较快时取变换logy或lny;当误差方差变化很快时取变换1/y;还有其他变换,如著名的Box-Cox幂变换 . 图(e)(f)表示选用回归模型是错误的。 (2)共线性 回归分析中很容易发生模型中两个或两个以上的自变量高度相关,从而引起最小二乘估计可能很不精确(称为共线性问题)。 在实际中最常见的问题是一些重要的自变量很可能由于在假设检验中t值不显著而被不恰当地剔除了。 共线性诊断问题就是要找出哪些变量间存在共线性关系。 (3)误差的独立性 回归分析之前,要检验误差的独立性。 若误差项不独立,那么回归模型的许多处理,包括误差项估计、假设检验等都将没有推导依据。 由于残差是误差的合理估计,因此检验统计量通常是建立在残差的基础上。 检验误差独立性的最常用方法,是对残差的一阶自相关性进行Durbin-Watson检验。 H0: 误差项是相互独立的;H1: 误差项是相关的 检验统计量: DW接近于0,表示残差中存在正自相关;如果DW接近于4,表示残差中存在负自相关;如果DW接近于2,表示残差独立性。 二、R语言实现 还是用lm()函数实现,不同是需要设置更复杂的formula格式: y~x1+x2——只考虑自变量的主效应(y=k1x1+k2x2),y~.表示全部自变量的主效应; y~x1+x2+x1: x2——考虑主效应和交互效应(y=k1x1+k2x2+k3x1x2); y~x1*x2——考虑全部主效应和交互效应的简写(效果同上); y~(x1+x2+x3)^2——考虑主效应以及至2阶以下的交互效应,相当于x1+x2+x3+x1: x2+x2: x3+x1: x3 y~x1%in%x2——x1含于x2,相当于x2+x2: x1 y~(x1+x2)^2-x1: x2——表示从(x1+x2)^2中去掉x1: x2 y~x1+I((x2+x3)^2)——使用I()函数,相当于用(x2+x3)^2计算出新变量h,然后y~x1+h function——在表达式中使用数学函数,例如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) lm.reg=lm(y~x1+x2+x3+x4+x5+x6,revenue) summary(lm.reg) Residuals: Min1QMedian3QMax -295.71-173.5226.5990.16370.01 Coefficients: EstimateStd.ErrortvaluePr(>|t|) (Intercept)6.046e+043.211e+0318.8298.12e-11*** x1-1.171e-018.638e-02-1.3560.19828 x23.427e-023.322e-021.0320.32107 x36.182e-014.103e-0215.0671.31e-09*** x4-5.152e-012.930e-02-17.5851.91e-10*** 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‘’1 Residualstandarderror: 234.8on13degreesoffreedom MultipleR-squared: 0.9999,AdjustedR-squared: 0.9999 F-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)表示添加一个新的变量 update(lm.reg,sqrt(.)~.)表示对因变量y开方,再重新拟合回归模型 lm.reg2<-update(lm.reg,.~.-x1-x2)#剔除自变量x1,x2 summary(lm.reg2) Residuals: Min1QMedian3QMax -325.62-147.5414.07108.28427.42 Coefficients: EstimateStd.ErrortvaluePr(>|t|) (Intercept)6.339e+042.346e+0327.0203.89e-14*** x36.584e-011.548e-0242.523<2e-16*** x4-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‘’1 Residualstandarderror: 233.6on15degreesoffreedom MultipleR-squared: 0.9999,AdjustedR-squared: 0.9999 F-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.42 Coefficients: EstimateStd.ErrortvaluePr(>|t|) (Intercept)6.339e+042.346e+0327.0203.89e-14*** x36.584e-011.548e-0242.523<2e-16*** x4-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‘’1 Residualstandarderror: 233.6on15degreesoffreedom MultipleR-squared: 0.9999,AdjustedR-squared: 0.9999 F-statistic: 3.476e+04on4and15DF,p-value: <2.2e-16 说明: 默认输出每步的结果(略),进行了3步回归,逐步剔除了自变量x1和x2,AIC值逐步减小,最终得到最优的模型。 drop1(lm.step) Singletermdeletions Model: y~x3+x4+x5+x6 DfSumofSqRSSAIC x319870181499520589316.40 x414111364341932417299.12 x5128739293692704250.52 x611851231003898224.47 lm.reg3<-lm(y~x3+x4+x5,revenue) summary(lm.reg3) Call: lm(formula=y~x3+x4+x5,data=revenue
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 语言 学习 系列 32 回归 分析
![提示](https://static.bdocx.com/images/bang_tan.gif)