R语言学习系列32回归分析报告.docx

上传人:b****6 文档编号:7354382 上传时间:2023-01-23 格式:DOCX 页数:23 大小:426.51KB
下载 相关 举报
R语言学习系列32回归分析报告.docx_第1页
第1页 / 共23页
R语言学习系列32回归分析报告.docx_第2页
第2页 / 共23页
R语言学习系列32回归分析报告.docx_第3页
第3页 / 共23页
R语言学习系列32回归分析报告.docx_第4页
第4页 / 共23页
R语言学习系列32回归分析报告.docx_第5页
第5页 / 共23页
点击查看更多>>
下载资源
资源描述

R语言学习系列32回归分析报告.docx

《R语言学习系列32回归分析报告.docx》由会员分享,可在线阅读,更多相关《R语言学习系列32回归分析报告.docx(23页珍藏版)》请在冰豆网上搜索。

R语言学习系列32回归分析报告.docx

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互不相关。

用最小二乘法原理,得到最优拟合效果的

值:

(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=“〞返回模型框架。

三、实例

例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

 

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

说明:

输出了残差信息Residuals;

回归系数估计值、标准误、t统计量值、p值,可得到回归方程:

height=+*age

回归系数p值〔<2e-16,4.43e-11〕很小,非常显著的≠0;***也表示显著程度非常显著。

拟合优度R2=0.9888>0.5,表示拟合程度很好。

F统计量=880,p值=,表示整个回归模型显著,适合估计height这一因变量。

coefficients(res.reg)#返回模型的回归系数估计值

(Intercept)age

confint(res.reg,parm="age",level=0.95)#输出参数age的置信区间,假设不指定parm将返回所有参数的置信区间

2.5%97.5%

fitted(res.reg)#输出回归模型的预测值

123456789101112

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

 

residuals(res.reg)#输出模型的残差

123456789101112

AIC(res.reg)#输出模型的AIC值

BIC(res.reg)#输出模型的BIC值

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+489*age,拟合优度R2=,拟合效果变得更好。

#用回归模型预测

ages<-data.frame(age=30:

34)

pre.res<-predict(res.reg2,ages,interval="prediction",level=0.95)#注意predict函数的第1个参数必须是回归模型的自变量数据构成的数据框或列表

fitlwrupr

 

〔二〕多元线性回归

一、根本原理

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

说明:

拟合优度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

〔三〕逐步回归

多元线性回归模型中,并不是所有的自变量都与因变量有显著关系,有时有些自变量的作用可以忽略。

这就需要考虑怎样从所有可能有关的自变量中挑选出对因变量有显著影响的局部自变量。

逐步回归的根本思想是,将变量一个一个地引入或剔出,引入或剔出变量的条件是“偏相关系数〞经检验是显著的,同时每引入或剔出一个变量后,对已选入模型的变量要进展逐个检验,将不显著变量剔除或将显著的变量引入,这样保证最后选入的所有自变量都是显著的。

逐步回归每一步只有一个变量引入或从当前的回归模型中剔除,当没有回归因子能够引入或剔出模型时,该过程停止。

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

说明:

默认输出每步的结果〔略〕,进展了3步回归,逐步剔除了自变量x1和x2,AIC值逐步减小,最终得到最优的模型。

drop1(lm.step)

Singletermdeletions

Model:

y~x3+x4+x5+x6

DfSumofSqRSSAIC

 

3<-lm(y~x3+x4+x5,revenue)

3)

Call:

lm(formula=y~x3+x4+x5,data=revenue)

Residuals:

Min1QMedian3QMax

-336.34-186.821.5289.46437.84

Coefficients:

EstimateStd.ErrortvaluePr(>|t|)

(Intercept)6.284e+042.494e+0325.1912.66e-14***

x36.614e-011.651e-0240.066<2e-16***

x4-5.467e-012.118e-02-25.8131.81e-14***

x5-1.412e-012.053e-02-6.8773.72e-06***

---

Signif.codes:

0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

Residualstandarderror:

250.5on16degreesoffreedom

MultipleR-squared:

0.9999,AdjustedR-squared:

0.9998

说明:

使用drop1()函数考察分别剔除每个自变量,AIC值变化的情况,可以看出不剔除x6与剔除x6,AIC值只从222.40变大到224.47,相对其它自变量变化很小。

所以,可以考虑剔除掉x6,重新做多元线性回归。

〔四〕回归诊断

回归分析之后,还需要从残差的随机性、强影响分析、共线性方面进展诊断。

一、残差诊断

(1)残差

y.res<-lm.reg3$residuals#回归模型的残差

y.fit<-predict(lm.reg3)#回归模型的预测值

plot(y.res~y.fit,main="残差图")#绘制残差图,以预测值作为横坐标

说明:

从图形看,残差分布比拟均匀,大致满足随机性。

shapiro.test(y.res)#残差的正态性检验

Shapiro-Wilknormalitytest

 

说明:

p值=0.2622>0.05,承受原假设,即残差服从正态分布。

(2)标准化残差

残差与数据的数量级有关,除以标准误差后得到标准化残差。

理想的标准化残差服从N(0,1).

rs<-rstandard(lm.reg3)#得到标准化残差

plot(rs~y.fit,main="标准残差图")

shapiro.test(rs)#标准化残差的正态性检验

Shapiro-Wilknormalitytest

data:

rs

(3)学生化残差

为了回避标准化残差的方差齐性假设,使用学生化残差。

rst<-rstudent(

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 幼儿教育 > 少儿英语

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1