第四课回归分析.docx

第四课回归分析.docx

《第四课回归分析.docx》由会员分享，可在线阅读，更多相关《第四课回归分析.docx（27页珍藏版）》请在冰豆网上搜索。

第四课回归分析

回归分析

一、回归分析：

是建立因变量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）###