R语言实验七

实验七线性回归分析、回归诊断、稳健回归与逐步回归【实验类型】验证性【实验学时】2学时【实验目的】1、掌握回归分析的基本原理及线性回归的求解方法；2、掌握为什么要回归诊断以及回归诊断的内容与解决办法；3、掌握稳健回归、抗干扰回归的使用条件及逐步回归的基本思想、求解过程。【实验内容】1、一元线性与多元线性回归的计算；2、回归结果的诊断与残差分析；3、稳健回归、抗干扰回归与逐步回归分析。【实验方法或步骤】第一部分、课件例题：#1例7.1blood-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)#线性回归，其中公式线性回归，其中公式Y~1+X1+X2与Y~X1+X2(隐含常数项)等价；但公式Y~0+X1+X2或Y~X1+X2-1表示所得的回归方程中不含有常数项表示所得的回归方程中不含有常数项summary(lm.sol)#提取结果#2例7.2newdata-data.frame(X1=80,X2=40)predict(lm.sol,newdata,interval=prediction)#预测区间predict(lm.sol,newdata,interval=confidence)#置信区间#3例7.3（1）回归估计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)#提取结果#（2）计算预测值并绘图new-data.frame(x=seq(0.10,0.24,by=0.01))pp-predict(lm.sol,new,interval=prediction);pp#预测pc-predict(lm.sol,new,interval=confidence);pc#置信par(mai=c(0.8,0.8,0.2,0.2))matplot(new$x,cbind(pp,pc[,-1]),type=l,#对矩阵绘图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)legend(0.1,63,#添加图例c(Points,Fitted,Prediction,Confidence),pch=c(19,NA,NA,NA),lty=c(NA,1,5,2),col=c(orange,blue,red,brown))#4例7.4##调取数据集data(anscombe);anscombe##作回归并输出回归系数等值ff-y~x#将公式赋给变量fffor(iin1:4){ff[2:3]-lapply(paste(c(y,x),i,sep=),as.name)assign(paste(lm.,i,sep=),lmi-lm(ff,data=anscombe))}GetCoef-function(n)summary(get(n))$coeflapply(objects(pat=lm\\.[1-4]$),GetCoef)##绘图op-par(mfrow=c(2,2),mar=.1+c(4,4,1,1),oma=c(0,0,2,0))for(iin1:4){ff[2:3]-lapply(paste(c(y,x),i,sep=),as.name)plot(ff,data=anscombe,col=red,pch=21,bg=orange,cex=1.2,xlim=c(3,19),ylim=c(3,13))abline(get(paste(lm.,i,sep=)),col=blue)}mtext(Anscombe's4Regressiondatasets,outer=TRUE,cex=1.5)par(op)#5例7.5##(1)回归rt-read.table(F:/文档/大学课程/R语言/ch07/blood.dat,header=T)#读数据lm.sol-lm(Y~X,data=rt);lm.sol#线性回归summary(lm.sol)##(2)残差图pre-fitted.values(lm.sol);pre#提取回归值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)plot(pre,rst,xlab=FittedValues,ylab=StandardizedResiduals)##(3)对残差作回归rt$res-res;rt#原数据增加一列残差lm.res-lm(abs(res)~X,data=rt)#用残差绝对值与X作回归lm.ressummary(lm.res)##(4)用权重修正回归方程s-lm.res$coefficients[1]+lm.res$coefficients[2]*rt$X#计算残差的标准差，即s=β0+β1*xlm.weg-lm(Y~X,data=rt,weights=1/s^2)#用方差(标准差的平方)的倒数作为样本点的权重，这样做可以减少非齐性方差带来的影响lm.weg;summary(lm.weg)#6例7.6的求解：#输入数据，作回归分析fr-read.table(F:/文档/大学课程/R语言/ch07/buy_income.dat,header=T)X-fr$X;Y-fr$Ylm.sol-lm(Y~X);summary(lm.sol)#作残差图，确定B-C变换中的参数λlibrary(MASS)#加载MASS程序包op-par(mfrow=c(2,2),mar=.4+c(4,4,1,1),oma=c(0,0,2,0))plot(fitted(lm.sol),resid(lm.sol),cex=1.2,pch=21,col=red,bg=orange,xlab=FittedValue,ylab=Residuals)#残差图boxcox(lm.sol,lambda=seq(0,1,by=0.1))#确定参数图##Box-Cox变换后,作回归分析lambda-0.55;Ylam-(Y^lambda-1)/lambda#B-C变换lm.lam-lm(Ylam~X);summary(lm.lam)plot(fitted(lm.lam),resid(lm.lam),cex=1.2,pch=21,col=red,bg=orange,xlab=FittedValue,ylab=NewResiduals)#变换后的残差图beta0-lm.lam$coefficients[1]#回归系数beta1-lm.lam$coefficients[2]curve((1+lambda*(beta0+beta1*x))^(1/lambda),from=min(X),to=max(X),col=blue,lwd=2,xlab=X,ylab=Y)#回归曲线points(X,Y,pch=21,cex=1.2,col=red,bg=orange)mtext(Box-CoxTransformations,outer=TRUE,cex=1.5);par(op)#7例7.7的求解：#(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)回归诊断influence.measures(lm.sol)#回归诊断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)由回归诊断结果得到18号和19号样本点是强影响点由图可知：19号可能是异常值点，18号可能是强影响点。#(3)处理异常点、强影响点并进行回归修正n-length(intellect$x)weights-rep(1,n)#定义权重向量weights[18]-0.5#对18号样本点的权系数重新赋值为0.5其它的均为1,以减少18号的影响lm.correct-lm(y~1+x,data=intellect,subset=-19,weights=weights)#剔除19号样本点summary(lm.correct)#(4)绘制回归直线和散点图进行比较验证attach(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)#原回归直线text(x[c(19,18)],y[c(19,18)],label=c(19,18),adj=c(1.5,0.3))detach()abline(lm.correct,col=red,lwd=2,lty=5)#新回归直线legend(30,120,c(Points,Regression,CorrectReg),pch=c(19,NA,NA),lty=c(NA,1,5),col=c(orange,blue,red))#(5)进一步回归诊断以检验结果op-par(mfrow=c(2,2),mar=0.4+c(4,4,1,1),oma=c(0,0,2,0))plot(lm.correct,1:4)#绘制回归诊断图par(op)#8例7.8的求解：france-data.frame(x1=c(149.3,161.2,171.5,175.5,180.8,190.7,202.1,212.4,226.1,231.9,239.0),x2=c(4.2,4.1,3.1,3.1,1.1,2.2,2.1,5.6,5.0,5.1,0.7),x3=c(108.1,114.8,123.2,126.9,132.1,137.7,146.0,154.1,162.3,164.3,167.6),y=c(15.9,16.4,19.0,19.1,18.8,20.4,22.7,26.5,28.1,27.6,26.3))#建立数据框lm.sol-lm(y~1+x1+x2+x3,data=france)lm.sum-summary(lm.sol)coef(lm.sum)#提取模型系数#计算矩阵XTX特征值的最小值X-as.matrix(france[1:3]);Xmin(eigen(cor(X))$values)#cor()用来得到相关系数矩阵，eigen()计算特征值与特征向量#9例7.9：对例7.8的法国经济分析数据作岭回归分析。library(MASS)lm.rid-lm.ridge(y~1+x1+x