第四课 回归分析.docx
《第四课 回归分析.docx》由会员分享,可在线阅读,更多相关《第四课 回归分析.docx(23页珍藏版)》请在冰豆网上搜索。
第四课回归分析
回归分析
一、回归分析:
是建立因变量Y与自变量X之间关系的模型,一元线性回归使用一个自变量X,多元线性回归使用超过一个自变量,如:
.
1线性回归模型
,
进行n次独立观测:
,i=1,2,…,n
称为残差向量
2回归系数显著性检验:
判断某一自变量
的系数
是否为0,如果统计量(
表示)的
值
(通常取0.05),认为
。
3回归方程显著性检验:
检查是否可用线性方程来处理数据。
即检验方程系数
是否全为0,如果统计量(
表示)的
值
,认为可以用线性方程来处理问题。
4相关性检验:
检验因变量
与自变量
相关性程度。
用相关系数的平方
来衡量,
的值接近1表示相关,接近0表示不相关。
二、线性回归模型的计算
1、lm()函数:
完成多元线性回归系数的估计、回归系数的检验、回归方程的检验等工作,返回值称为拟合结果的对象,本质上是一个具有类属性值lm的列表.
格式:
lm(formula,data,subset,weights,na.action,method=”qr”,model=TRUE,x=FALSE,y=FALSE,qr=TRUE,singular.ok=TRUE,contrasts=NULL,offset,…)
formula为模型公式,形如y~1+x1+x2的形式,表示常数项、
的系数和
的系数,去掉
公式中的1,其意义不变。
拟合成齐次线性模型:
y~0+x1+x2或y~-1+x1+x2的形式;
data为数据框,由样本数据构成;
subset为可选项,表示所使用的样本子集;
weights为可选向量,表示对应样本的权重;
na.action为函数,表示当数据中出现缺失数据(NA)的处理方法;
method为估计回归系数的计算方法,默认值为“qr”;
model、x、y、qr为逻辑变量,如果取TRUE,函数的返回值将给出模型的框架、模型矩阵、
响应变量,及QR分解;
singular.ok为逻辑变量,取FALSE表示奇异值拟合是错误的;
contrasts为可选列表;
offset为NULL,或者是数值向量。
…为附加参数。
2、summary()函数
lm()函数的返回值称为拟合结果的对象,是一个有类属性值lm的列表,有model,coefficients,residuals(残差)等成员,为了获得更多信息,summary()常与lm()函数一起使用,格式:
summary(object,correlation=FALSE,symbolic.cor=FALSE,…)
参数意义:
object为lm()函数生成的对象,correlation为逻辑变量,取TRUE表示给出估计参数的相关矩
阵,symbolic.cor为逻辑变量,取TRUE表示用符号形式给出估计参数的相关矩阵,此参数只有当correlation=TRUE时才有效。
例1年龄相等情况下,建立血压的收缩压
与体重
(kg),年龄
(岁数)有关,建立
与
,
的线性回归方程。
解
blood<-data.frame(
X1=c(76.0,91.5,85.5,82.5,79.0,80.5,74.5,79.0,85.0,76.5,82.0,95.0,92.5),
X2=c(50,20,20,30,30,50,60,50,40,55,40,40,20),
Y=c(120,141,124,126,117,125,123,125,132,123,132,155,147)
)
lm.sol<-lm(Y~1+X1+X2,data=blood)
summary(lm.sol)
运行结果:
Call:
lm(formula=Y~1+X1+X2,data=blood)###调出函数使用的模型
Residuals:
(残差)###列出残差
Min1QMedian3QMax
-4.0404-1.01830.46400.69084.3274
Coefficients:
##列出估计值、估计值的标准差、t统计量、对应t统计量的P值
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)-62.9633616.99976-3.7040.004083**
X12.136560.1753412.1852.53e-07***
X20.400220.083214.8100.000713***
---
极为显著***,高度显著**,显著*
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
2.854on10degreesoffreedom###残差的标准差
MultipleR-squared:
0.9461,AdjustedR-squared:
0.9354##相关系数的平方,修正相关系数的平方
F-statistic:
87.84on2and10DF,p-value:
4.531e-07###F统计量,F统计量对应的P值
回归方程与回归系数的检验都是显著的,回归方程为:
3预测区间与置信区间
如果回归效果显著,就可以利用回归方程进行预测,即对给定的回归自变量值,预测因变量可能的取值范围,是一个区间估计问题。
假设
=
,
回归方程的真实值(未知):
相应的估计值:
predict()函数计算
的估计值、预测区间和置信区间。
格式:
predict(object,newdata,se.fit=FALSE,scale=NULL,df=Inf,
interval=c(“none”,”confidence”,”prediction”),level=0.95,type=c(“response”,”terms”),terms=NULL,
na.action=na.pass,pred.var=res.var/weights,weights=1,…)
参数意义:
object为lm()函数得到的对象;
newdata为数据框,由预测点构成,如取默认值,将计算已知数据的回归值;
se.fit为逻辑变量,取TRUE表示输出预测值的标准差、自由度和残差尺度信息;
scale为计算标准差的尺度参数;
df为尺度参数的自由度;
interval为计算的区间类型,取“none”(默认值)表示不计算,取”confidence”表示计算置信区
间,取”prediction”表示计算预测区间;
level为置信水平,默认值为0.95,在计算置信区间和预测区间时用到;
type为预测类型,或者取“response”(默认值),或者取”terms”,当type=”terms”时,为全部
选择项;
terms为选择项,只能取NULL,1,2等;
na.action为函数,表示处理newdata中有缺失数据(NA)时的处理方法,默认预测为NA.pred.var
预测区间的方差;
weights为数值向量或单侧模型公式,用于预测时方差的权;
…附加参数。
例2继续上面的例1,设
,求
的估计值、预测区间和置信区间(置信水平为0.95)。
解
blood<-data.frame(
X1=c(76.0,91.5,85.5,82.5,79.0,80.5,74.5,79.0,85.0,76.5,82.0,95.0,92.5),
X2=c(50,20,20,30,30,50,60,50,40,55,40,40,20),
Y=c(120,141,124,126,117,125,123,125,132,123,132,155,147)
)
lm.sol<-lm(Y~1+X1+X2,data=blood)
newdata<-data.frame(X1=80,X2=40)
predict(lm.sol,newdata,interval="prediction")
predict(lm.sol,newdata,interval="confidence")
例分析合金强度
与合金中含碳量
之间的关系:
(它们之间的一组数据如下)
(1)完成一元线性回归的计算;
(2)计算自变量
在区间[0.10,0.23]内的回归方程的预测估计值、预测区间和置信区间(取
),并将数据点、预测估计曲线、预测区间曲线和置信区间曲线画在同一张图上。
解
##回归估计
x<-c(0.10,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.20,0.21,0.23)
y<-c(42.0,43.5,45.0,45.5,45.0,47.5,49.0,53.0,50.0,55.0,55.0,60.0)
lm.sol<-lm(y~1+x)
summary(lm.sol)
##计算预测值,并绘图
new<-data.frame(x=seq(0.10,0.24,by=0.01))##取x的范围值
pp<-predict(lm.sol,new,interval="prediction")
pc<-predict(lm.sol,new,interval="confidence")
par(mai=c(0.8,0.8,0.2,0.2))
matplot(new$x,cbind(pp,pc[,-1]),type="l",######绘出矩形各列pp和pc[,-1]共5列
xlab="X",ylab="Y",lty=c(1,5,5,2,2),
col=c("blue","red","red","brown","brown"),lwd=2)
points(x,y,cex=1.4,pch=21,col="red",bg="orange")###真实值x,y
legend(0.1,63,
c("Points","Fitted","Prediction","Confidence"),
###如果是字符串,用c(“”),而不用expr1数学表达式
pch=c(19,NA,NA,NA),lty=c(NA,1,5,2),
col=c("orange","blue","red","brown"))
savePlot("predict",type="eps")
三、回归诊断
为什么要做回归诊断?
因为有些数值例子得到的回归方程尽管通过了t检验和F检验,但做线性回归方程还是有问题的。
如Anscombe在1973年构造的数据集给出4组数据拟合简单的线性模型:
,得到的结果:
[[1]]
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)3.00009091.12474682.6673480.025734051
x10.50009090.11790554.2414550.002169629
[[2]]
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)3.0009091.12530242.6667580.025758941
x20.5000000.11796374.2385900.002178816
[[3]]
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)3.00245451.12448122.6700800.025619109
x30.49972730.11787774.2393720.002176305
[[4]]
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)3.00172731.12392112.6707630.025590425
x40.49990910.11781894.2430280.002164602
并且进一步检验4组数据的
、F值检验全合格,但绘出图分析4组数据做线性回归模型不是全合格:
只有第一组数据适合做线性回归模型。
1残差检验
检验模型误差是否满足正态性和方差齐性。
1)用残差图检验,残差为纵坐标,拟合值或数据观测序号或数据观测时间做横坐标,做出的散点图称为残差图。
(a)(b)(c)
正常情况异方差情况非线性情况
有异常值点的情况
检查1)有无有异常值点;2)画出残差的Q-Q散点图,若这些点位于一条直线上,则残差服从正态分布,否则不然。
注意:
上面图中只有(a)图是正常的,(b)图是喇叭口,回归值的大小与残差的波动有关,方差的假设有问题,(c)图表明线性模型不合适。
residuals()(或resid()函数)计算回归模型的残差,rstandard()计算回归模型的标准化残差,rstudent()计算回归模型的化学生化残差(删除第i个样本数据后得到的标准化残差)
格式:
residuals(object,...)resid(object,...)
rstandard(model,infl=lm.influence(model,do.coef=FALSE),
sd=sqrt(deviance(model)/df.residual(model)),...)
rstudent(model,infl=lm.influence(model,do.coef=FALSE),res=infl$wt.res,...)
参数object或model为lm生成的对象;infl为lm.influence返回值得到的影响结构;sd为模型的标准差;res为模型残差。
plot()函数绘出各种残差的散点图,格式:
plot(x,which=c(1:
3,5),...)),参数x由lm生成的对象,which是开关变量,其值1~6,默认子集(1,2,3,5),即绘出1、2、3和5号散点图,该函数可以绘出6张诊断图,第一张残差与预测值散点图,第二张残差的正态Q-Q图,第三张标准差的平方根与预测值的散点图等。
例20~60成年女性的血压与年龄之间的回归关系,并作残差分析
解
##
(1)回归
rt<-read.table("blood.dat",header=TRUE)
lm.sol<-lm(Y~X,data=rt);lm.sol
summary(lm.sol)
##
(2)残差图
pre<-fitted.values(lm.sol)###拟合值
res<-residuals(lm.sol)###残差
rst<-rstandard(lm.sol)####标准差
par(mai=c(0.9,0.9,0.2,0.1))
plot(pre,res,xlab="FittedValues",ylab="Residuals")
savePlot("resid-1",type="eps")
plot(pre,rst,xlab="FittedValues",ylab="StandardizedResiduals")
savePlot("resid-2",type="eps")
##(3)对残差作回归
rt<-read.table("blood.dat",header=TRUE)
rt$res<-res###########给rt数据框增加一列res
lm.res<-lm(abs(res)~X,data=rt);################利用残差的绝对值与自变量做回归
lm.res
summary(lm.res)
##(4)修正回归方程
s<-lm.res$coefficients[1]+lm.res$coefficients[2]*rt$X#####计算残差的标准差
lm.weg<-lm(Y~X,data=rt,weights=1/s^2);#####用标准差平方的倒数作权重
################权重的选取可减少非齐性方差带来的影响
lm.weg
summary(lm.weg)
2影响分析
影响分析就是检查是否有对估计有异常大数据(高杠杆点)情况。
影响分析的方法有:
DFFITS准则、Cook距离、COVRATIO准则、帽子值和帽子矩阵。
相关的R函数有:
influence.measures(),dffits(),dfbeta(),dfbetas(),cooks.distance(),covratio(),hatvalues(),hat();
格式:
influence.measures(model)
dffits(model,infl=,res=)
dfbeta(model,infl=lm.influence(model,do.coef=TRUE),...)
dfbetas(model,infl=lm.influence(model,do.coef=TRUE),...)
covratio(model,infl=lm.influence(model,do.coef=FALSE),res=weighted.residuals(model))
cooks.distance(model,infl=lm.influence(model,do.coef=FALSE),res=weighted.residuals(model),sd=sqrt(deviance(model)/df.residual(model)),hat=infl$hat,...)
hatvalues(model,infl=lm.influence(model,do.coef=FALSE),...)
hat(x,intercept=TRUE)
参数model为lm生成的对象,infl为lm.influence返回值得到的影响结构,res为模型残差,x为设计矩阵。
例通过智力测试数据建立智力随年龄变化的关系(在exam0606.R文件中)。
解
第一步:
运行,计算回归系数,并做回归系数与回归方程检验
##1.回归分析
intellect<-data.frame(
x=c(15,26,10,9,15,20,18,11,8,20,7,9,10,11,11,10,12,42,17,11,10),######年龄(月)
y=c(95,71,83,91,102,87,93,100,104,94,113,96,83,84,102,100,105,57,121,86,100)
##智力
)
lm.sol<-lm(y~1+x,data=intellect)
summary(lm.sol)
第二步做回归诊断
##2.回归诊断
op<-par(mfrow=c(2,2),mar=0.4+c(4,4,1,1),oma=c(0,0,2,0))
plot(lm.sol,1:
4)#####通过绘图分析有无异常值点
par(op)
savePlot("diagnoses1",type="eps")
influence.measures(lm.sol)#####通过函数计算分析有无异常值点,下面是运行结果:
dfb.1_dfb.xdffitcov.rcook.dhatinf
10.016640.003280.041271.1668.97e-040.0479
20.18862-0.33480-0.402521.1978.15e-020.1545
.....................
180.83112-1.11275-1.155782.9596.78e-010.6516*
190.143480.273170.853740.3962.23e-010.0531*
20-0.207610.10544-0.263851.0433.45e-020.0567
210.02791-0.016220.032981.1875.74e-040.0628
influence.measures()函数得到回归诊断结果共7列,1,2列(dfbetas)对应常数和变量x值,3列是DFFITS准则值,4列是COVRATIO准则值,5列是Cook距离,6列高杠杆值,7列是影响点记号。
(18,19号样本是强影响点记号)。
第三步处理强影响点
##3.回归修正
n<-length(intellect$x)
weights<-rep(1,n);weights[18]<-0.5
lm.correct<-lm(y~1+x,data=intellect,subset=-19,weights=weights)
summary(lm.correct)
第四步验证
##4.验证####在脚本里做
intellect<-data.frame(
x=c(15,26,10,9,15,20,18,11,8,20,7,9,10,11,11,10,12,42,17,11,10),######年龄(月)
y=c(95,71,83,91,102,87,93,100,104,94,113,96,83,84,102,100,105,57,121,86,100)
##智力
)
lm.sol<-lm(y~1+x,data=intellect)
attach(intellect)######把intellect调入内存,可操作它的元素
par(mai=c(0.8,0.8,0.2,0.2))
plot(x,y,cex=1.2,pch=21,col="red",bg="orange")########真值
abline(lm.sol,col="blue",lwd=2)###########修正前,lm.sol是拟合值
text(x[c(19,18)],y[c(19,18)],label=c("19","18"),adj=c(1.5,0.3))
###adj()调整符号坐标位置
detach()
abline(lm.correct,col="red",lwd=2,lty=5)##