实用回归分析实验 (R软件的使用)
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
式
class(X1)
scan命令的用法：
③scan的第三种用法
• x<-scan() • 输入数据后回车
4、调用R 自带的数据
• data() • women • 可以对此数据做统计分析，一元回归
二、线性回归模型（用命令计算参数、就参数
的假设检验等）
• 1、做回归系数的估计，回归系数的显著性 检验、回归系数的标准误差、F检验、决定 系数、调整的决定系数用lm命令。
2、误差的独立性检验（p107）
• 目的检验误差是否是独立的 • durbinWatsonTest##D—W检验 • ###杜宾和沃特森于1951年提出的适用于小
样本的一种检验方法。DW检验只能使用用 检验随机扰动项具有一阶自回归形式的序 列相关性的检验。这种检验是建立计量经 济学模型中常用方法。
and)),col=4)
• library(car)
• spreadLevelPlot(lm.baoxian)
• 画出标准化残差绝对值与拟合值之间的关 系
2、误差的独立性检验
• durbinWatsonTest(lm.baoxian)##误差的独 立性（结果看p值较大，暂时接受原假设， 可认为误差独立）
2、qq图检验 plot(lm.sol,2)
Q-Q图的解释
7、残差分析中几个概念
• 标准化残差（内学生化残差）
学生化残差SRE
• y.rst<-rstandard(lm.baoxian)###标准化残 差
教材120删除残差
8、残差分析图及解释
• y.fit<-predict(lm.baoxian)##yi的预测值
• left<-A[,1]-A[,2]*qt(1-alpha/2, df)
• right<-A[,1]+A[,2]*qt(1-alpha/2, df)
• rowname<-dimnames(A)[[1]]
• colname<-c("Estimate", "Left", "Right")
• matrix(c(A[,1], left, right), ncol=3,
图分析：
• 红线是用数据和核估计法拟合的残差的概 率密度函数图，蓝色的线是理论的正态分 布的概率密度函数图。从上图中可看成出 数据是双峰的，红色线和蓝色线差别很大 可以认为残差不服从正态分布。
残差的直方图hist
• hist
自己编写命令：
• y.res<-resid(lm.baoxian) • 标准化残差直方图和概率密度图（当然也可以不标准
ቤተ መጻሕፍቲ ባይዱ
残差的直方图及其和正态分布概率密度函数 的比较
• residplot<-function(fit,nbreaks=10){ • z<-rstudent(fit) • hist(z,breaks=nbreaks,freq=FALSE,xlab="Studentized Residual", • main="Distribution of Errors") • rug(jitter(z),col="brown") • curve(dnorm(x,mean=mean(z),sd=sd(z)),add=TRUE,col="blue",lwd=
3、线性crPlots(lm.baoxian)
4、同方差性（四种方法）
• 1.plot(lm.baoxian)
library(car) spreadLevelPlot(lm.baoxian) • 画出标准化残差绝对值与拟合值之间的关系,从图形中可知点在水平的
最佳拟合曲线周围呈水平随机分布，基本满足同方差的假定，代码结 果建议幂次变换。
• y<c(26.2,17.8,31.3,23.1,27.5,36.0,14.1,22.3, 19.6,31.3,24.0,17.3,43.2,36.4,26.1)
• baoxian<-data.frame(x=x,y=y)
• lm.baoxian<-lm(y~x,data=baoxian) • summary(lm.baoxian) • qqPlot(lm.baoxian)
6、残差分析
• ①提取普通残差 • names(lm.baoxian) • residuals(lm.baoxian) • ##用命令直接显示残差 • 或者 • y.res<-resid(lm.baoxian)
②残差的正态性检验（两种方法画图和统计推断）
• 1、y.res<-resid(lm.baoxian) • shapiro.test(y.res)
•
dimnames = list(rowname, colname ))
•}
回归系数的置信区间另外一种求法(利用命令confint)
5、新值x=4.2的预测值及其置信区 间
• > new<-data.frame(x=4.2) • > y.fit<-
predict(lm.baoxian,new,interval="predictio n",level=0.95) • > y.fit • fit lwr upr • 1 30.93912 25.71222 36.16601
2) • lines(density(z)$x,density(z)$y,col="red",lwd=2,lty=2) • legend("topright",legend= • c("Normal Curve","Kernel Density
Curve"),lty=1:2,col=c("blue","red"),cex=.7) •} • residplot(lm.baoxian)
c(26.2,17.8,31.3,23.1,27.5,36.0,14.1,22.3,19.6,31.3,24.0,17.3,43.2,3 6.4,26.1) • length(x)##检验输入的x数据长度是否正确 • length(y) )##检验输入的x数据长度是否正确 • baoxian<-data.frame(x=x1,y=y1)##建立以x,y为列名，x1,y1向量的值 ###为数据数据框 • baoxian #显示前面建立的数据框，检查数据的正确性
• plot(y.res~y.fit) • y.rst<-rstandard(lm.baoxian)###标准化残差 • plot(y.rst~y.fit)###画标准化残差图 • plot(lm.baoxian) • 残差ＱＱ图＃＃检验残差是否符合线性回归模型的假定 • plot(lm.baoxian,2) • plot(y.res~x)##横坐标为x的残差图 • plot(lm.baoxian,1)#普通残差与拟合值的残差的图 • plot(lm.baoxian,3)#标准化残差的开方与拟合值的残差图 • plot(lm.baoxian,4)#画cook统计量 • par(mfrow=c(2,2)) • plot(lm.baoxian)#画四个图
化） • y.restand<-rstudent(lm.baoxian) • hist(y.restand,freq=FALSE) • lines(density(y.restand),col=5) • x<-seq(-2,2,0.1) • lines(x,dnorm(x,mean=mean(y.restand),sd=sd(y.rest
实用回归分析实验课
R软件的使用
R软件与线性回归模型相关的函数
• 输入数据： • 数据格式要求是数据框，下面以 P17 example2.1为例子讲述 一、当要处理的数据较少时，可直接在R的窗口输入数据 1、x1<-c(3.4,1.8,4.6,2.3,3.1,5.5,0.7,3.0,2.6,4.3,2.1,1.1,6.1,4.8,3.8) • y1<-
3、对回归方程做方差分析 anova(lm.baoxian)
4、回归系数的置信区间
(p是自变量的个数)
• ### 求线性模型系数的区间估计
• beta.int<-function(fm,alpha=0.05){
• A<-summary(fm)$coefficients
• df<-fm$df.residual
• durbinWatsonTest(lm.baoxian)
3、线性
• 通过成分残差图也称偏残差图可以看因变 量与自变量之间是否成线性关系，也可以 看看是否有不同于已设定的线性模型的系 统偏差。crPlots
4、同方差性检验
• 1.残差图plot(lm.baoxian,1) • 2.spreadLevelPlot(lm.baoxian)##加载包
• length(x)
• length(y) • plot(x,y)##画散点图 • cor(x,y)##计算相关系数 • cor(x,y)^2##计算决定系数或者用lm命令做的结果中找 • baoxian<-data.frame(x=x,y=y)##用数据x,y做一元回归
• summary(baoxian)
• lm.baoxian<-lm(y~x,data=baoxian)
• summary(lm.baoxian)
• ##汇总模型的信息，可以得到残差，参数 估计，参数估 计的标准差以及t检验，F 检 验
2、在散点图上加上拟合回归直线
• plot(x,y,main="散点图和拟合直线") • abline(lm(y~x))##在散点图上加上拟合回归直线
ncvTest(lm.baoxian)#检验同方差性
5、线性模型的综合检验
• gvlma包中gvlma命令
• Pena和Slate(2006)年编写，能对线性模型 假设进行综合验证，同时还能做偏倚度、 峰度和异方差性的检验。即它给模型提供 一个单独的综合检验（通过或者不通过）
• x<-c(3.4,1.8,4.6,2.3,3.1,5.5,0.7,3.0,2.6,4.3,2.1,1.1,6.1,4.8,3.8)
• y<c(26.2,17.8,31.3,23.1,27.5,36.0,14.1,22.3,19.6,31.3,24.0,17.3,43. 2,36.4,26.1)
残差分析：
• 从上图中看残差的散点图均值在0，再看红 色的线条，说明残差的散点图有一定的规 律，不满足误差等方差的假定，也不满足 线性回归的假设。红线表明残差中还有些 规律，我们的一元线性回归模型没有找出 规律，即模型不合适。
Q-Q图解释
• 从图中看实际数据的分位点和标准正态分 布的理论分位点基本在一条直线上。有几 个点不在直线上也正常，数据的随机性
三、回归诊断(要安装加载包“car”)
• 1、误差的正态性检验 • 基础包内的 plot(lm.baoxian,2)##检验残差是否符合正态
分布 qqPlot(lm.baoxian)##car包中命令 residplot画残差的直方图和概率密度图与理论
的正态分布进行比较。 K-S检验 ks.test shapiro.test ##正态性检验
• （若不满足你将会看到一个非水平的曲线。）
x<-x^(0.3817295 ) lm.baoxian1<-lm(y~x) spreadLevelPlot(lm.baoxian1)
斯皮尔曼等级相关系数
• y.res<-residuals(lm.baoxian) cor.test(y.res,x,method="spearman")
②用scan命令可以把数据读入后转为矩阵
X <- matrix(scan(“example22.data”, 0),nrow=20,ncol=3, byrow=TRUE)##读入数据 后变成一个20*3的矩阵，按行读入
colnames(X)<-c(“year”,“y”,“x”) ##改变矩阵的列名 X##显示X中的数据 class(X)##数据结构为矩阵 X1<-as.data.frame(X)##数据结构转换为数据框形
###car功能类似于基础包中 plot(lm.baoxian,1) • 2.cor.test##选择方法是Spearman相关系数 • 是残差和x的斯皮尔曼等级相关系数检验 • 3.ncvTest(lm.baoxian)
• library(car)
• x<c(3.4,1.8,4.6,2.3,3.1,5.5,0.7,3.0,2.6,4.3,2.1 ,1.1,6.1,4.8,3.8)
class(X1)
scan命令的用法：
③scan的第三种用法
• x<-scan() • 输入数据后回车
4、调用R 自带的数据
• data() • women • 可以对此数据做统计分析，一元回归
二、线性回归模型（用命令计算参数、就参数
的假设检验等）
• 1、做回归系数的估计，回归系数的显著性 检验、回归系数的标准误差、F检验、决定 系数、调整的决定系数用lm命令。
2、误差的独立性检验（p107）
• 目的检验误差是否是独立的 • durbinWatsonTest##D—W检验 • ###杜宾和沃特森于1951年提出的适用于小
样本的一种检验方法。DW检验只能使用用 检验随机扰动项具有一阶自回归形式的序 列相关性的检验。这种检验是建立计量经 济学模型中常用方法。
and)),col=4)
• library(car)
• spreadLevelPlot(lm.baoxian)
• 画出标准化残差绝对值与拟合值之间的关 系
2、误差的独立性检验
• durbinWatsonTest(lm.baoxian)##误差的独 立性（结果看p值较大，暂时接受原假设， 可认为误差独立）
2、qq图检验 plot(lm.sol,2)
Q-Q图的解释
7、残差分析中几个概念
• 标准化残差（内学生化残差）
学生化残差SRE
• y.rst<-rstandard(lm.baoxian)###标准化残 差
教材120删除残差
8、残差分析图及解释
• y.fit<-predict(lm.baoxian)##yi的预测值
• left<-A[,1]-A[,2]*qt(1-alpha/2, df)
• right<-A[,1]+A[,2]*qt(1-alpha/2, df)
• rowname<-dimnames(A)[[1]]
• colname<-c("Estimate", "Left", "Right")
• matrix(c(A[,1], left, right), ncol=3,
图分析：
• 红线是用数据和核估计法拟合的残差的概 率密度函数图，蓝色的线是理论的正态分 布的概率密度函数图。从上图中可看成出 数据是双峰的，红色线和蓝色线差别很大 可以认为残差不服从正态分布。
残差的直方图hist
• hist
自己编写命令：
• y.res<-resid(lm.baoxian) • 标准化残差直方图和概率密度图（当然也可以不标准
ቤተ መጻሕፍቲ ባይዱ
残差的直方图及其和正态分布概率密度函数 的比较
• residplot<-function(fit,nbreaks=10){ • z<-rstudent(fit) • hist(z,breaks=nbreaks,freq=FALSE,xlab="Studentized Residual", • main="Distribution of Errors") • rug(jitter(z),col="brown") • curve(dnorm(x,mean=mean(z),sd=sd(z)),add=TRUE,col="blue",lwd=
3、线性crPlots(lm.baoxian)
4、同方差性（四种方法）
• 1.plot(lm.baoxian)
library(car) spreadLevelPlot(lm.baoxian) • 画出标准化残差绝对值与拟合值之间的关系,从图形中可知点在水平的
最佳拟合曲线周围呈水平随机分布，基本满足同方差的假定，代码结 果建议幂次变换。
• y<c(26.2,17.8,31.3,23.1,27.5,36.0,14.1,22.3, 19.6,31.3,24.0,17.3,43.2,36.4,26.1)
• baoxian<-data.frame(x=x,y=y)
• lm.baoxian<-lm(y~x,data=baoxian) • summary(lm.baoxian) • qqPlot(lm.baoxian)
6、残差分析
• ①提取普通残差 • names(lm.baoxian) • residuals(lm.baoxian) • ##用命令直接显示残差 • 或者 • y.res<-resid(lm.baoxian)
②残差的正态性检验（两种方法画图和统计推断）
• 1、y.res<-resid(lm.baoxian) • shapiro.test(y.res)
•
dimnames = list(rowname, colname ))
•}
回归系数的置信区间另外一种求法(利用命令confint)
5、新值x=4.2的预测值及其置信区 间
• > new<-data.frame(x=4.2) • > y.fit<-
predict(lm.baoxian,new,interval="predictio n",level=0.95) • > y.fit • fit lwr upr • 1 30.93912 25.71222 36.16601
2) • lines(density(z)$x,density(z)$y,col="red",lwd=2,lty=2) • legend("topright",legend= • c("Normal Curve","Kernel Density
Curve"),lty=1:2,col=c("blue","red"),cex=.7) •} • residplot(lm.baoxian)
c(26.2,17.8,31.3,23.1,27.5,36.0,14.1,22.3,19.6,31.3,24.0,17.3,43.2,3 6.4,26.1) • length(x)##检验输入的x数据长度是否正确 • length(y) )##检验输入的x数据长度是否正确 • baoxian<-data.frame(x=x1,y=y1)##建立以x,y为列名，x1,y1向量的值 ###为数据数据框 • baoxian #显示前面建立的数据框，检查数据的正确性
• plot(y.res~y.fit) • y.rst<-rstandard(lm.baoxian)###标准化残差 • plot(y.rst~y.fit)###画标准化残差图 • plot(lm.baoxian) • 残差ＱＱ图＃＃检验残差是否符合线性回归模型的假定 • plot(lm.baoxian,2) • plot(y.res~x)##横坐标为x的残差图 • plot(lm.baoxian,1)#普通残差与拟合值的残差的图 • plot(lm.baoxian,3)#标准化残差的开方与拟合值的残差图 • plot(lm.baoxian,4)#画cook统计量 • par(mfrow=c(2,2)) • plot(lm.baoxian)#画四个图
化） • y.restand<-rstudent(lm.baoxian) • hist(y.restand,freq=FALSE) • lines(density(y.restand),col=5) • x<-seq(-2,2,0.1) • lines(x,dnorm(x,mean=mean(y.restand),sd=sd(y.rest
实用回归分析实验课
R软件的使用
R软件与线性回归模型相关的函数
• 输入数据： • 数据格式要求是数据框，下面以 P17 example2.1为例子讲述 一、当要处理的数据较少时，可直接在R的窗口输入数据 1、x1<-c(3.4,1.8,4.6,2.3,3.1,5.5,0.7,3.0,2.6,4.3,2.1,1.1,6.1,4.8,3.8) • y1<-
3、对回归方程做方差分析 anova(lm.baoxian)
4、回归系数的置信区间
(p是自变量的个数)
• ### 求线性模型系数的区间估计
• beta.int<-function(fm,alpha=0.05){
• A<-summary(fm)$coefficients
• df<-fm$df.residual
• durbinWatsonTest(lm.baoxian)
3、线性
• 通过成分残差图也称偏残差图可以看因变 量与自变量之间是否成线性关系，也可以 看看是否有不同于已设定的线性模型的系 统偏差。crPlots
4、同方差性检验
• 1.残差图plot(lm.baoxian,1) • 2.spreadLevelPlot(lm.baoxian)##加载包
• length(x)
• length(y) • plot(x,y)##画散点图 • cor(x,y)##计算相关系数 • cor(x,y)^2##计算决定系数或者用lm命令做的结果中找 • baoxian<-data.frame(x=x,y=y)##用数据x,y做一元回归
• summary(baoxian)
• lm.baoxian<-lm(y~x,data=baoxian)
• summary(lm.baoxian)
• ##汇总模型的信息，可以得到残差，参数 估计，参数估 计的标准差以及t检验，F 检 验
2、在散点图上加上拟合回归直线
• plot(x,y,main="散点图和拟合直线") • abline(lm(y~x))##在散点图上加上拟合回归直线
ncvTest(lm.baoxian)#检验同方差性
5、线性模型的综合检验
• gvlma包中gvlma命令
• Pena和Slate(2006)年编写，能对线性模型 假设进行综合验证，同时还能做偏倚度、 峰度和异方差性的检验。即它给模型提供 一个单独的综合检验（通过或者不通过）
• x<-c(3.4,1.8,4.6,2.3,3.1,5.5,0.7,3.0,2.6,4.3,2.1,1.1,6.1,4.8,3.8)
• y<c(26.2,17.8,31.3,23.1,27.5,36.0,14.1,22.3,19.6,31.3,24.0,17.3,43. 2,36.4,26.1)
残差分析：
• 从上图中看残差的散点图均值在0，再看红 色的线条，说明残差的散点图有一定的规 律，不满足误差等方差的假定，也不满足 线性回归的假设。红线表明残差中还有些 规律，我们的一元线性回归模型没有找出 规律，即模型不合适。
Q-Q图解释
• 从图中看实际数据的分位点和标准正态分 布的理论分位点基本在一条直线上。有几 个点不在直线上也正常，数据的随机性
三、回归诊断(要安装加载包“car”)
• 1、误差的正态性检验 • 基础包内的 plot(lm.baoxian,2)##检验残差是否符合正态
分布 qqPlot(lm.baoxian)##car包中命令 residplot画残差的直方图和概率密度图与理论
的正态分布进行比较。 K-S检验 ks.test shapiro.test ##正态性检验
• （若不满足你将会看到一个非水平的曲线。）
x<-x^(0.3817295 ) lm.baoxian1<-lm(y~x) spreadLevelPlot(lm.baoxian1)
斯皮尔曼等级相关系数
• y.res<-residuals(lm.baoxian) cor.test(y.res,x,method="spearman")
②用scan命令可以把数据读入后转为矩阵
X <- matrix(scan(“example22.data”, 0),nrow=20,ncol=3, byrow=TRUE)##读入数据 后变成一个20*3的矩阵，按行读入
colnames(X)<-c(“year”,“y”,“x”) ##改变矩阵的列名 X##显示X中的数据 class(X)##数据结构为矩阵 X1<-as.data.frame(X)##数据结构转换为数据框形
###car功能类似于基础包中 plot(lm.baoxian,1) • 2.cor.test##选择方法是Spearman相关系数 • 是残差和x的斯皮尔曼等级相关系数检验 • 3.ncvTest(lm.baoxian)
• library(car)
• x<c(3.4,1.8,4.6,2.3,3.1,5.5,0.7,3.0,2.6,4.3,2.1 ,1.1,6.1,4.8,3.8)