R语言实验七.docx

精选优质文档-倾情为你奉上实验七 线性回归分析、回归诊断、稳健回归与逐步回归【实验类型】验证性【实验学时】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) # 线性回归，其中公式线性回归，其中公式 Y1+X1+X2 与 YX1+X2( 隐含常数项 ) 等价；但公式Y0+X1+X2 或 YX1+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(i in 1:4) ff2: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(i in 1:4) ff2: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's 4 Regression data sets", outer = TRUE, cex=1.5)par(op)#5 例7.5# (1) 回归rt<-read.table("F:/文档/大学课程/R语言/ch07/blood.dat", header=T) # 读数据lm.sol<-lm(YX, 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="Fitted Values", ylab="Residuals")plot(pre, rst, xlab="Fitted Values", ylab="Standardized Residuals")# (3) 对残差作回归rt$res<-res;rt # 原数据增加一列残差lm.res<-lm(abs(res)X, data=rt) # 用残差绝对值与 X 作回归lm.ressummary(lm.res)# (4) 用权重修正回归方程s<-lm.res$coefficients1+lm.res$coefficients2*rt$X # 计算残差的标准差，即 s= 0 + 1 *xlm.weg<-lm(YX, data=rt, weights=1/s2) # 用方差 ( 标准差的平方 ) 的倒数作为样本点的权重，这样做可以减少非齐性方差带来的影响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(YX); 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="Fitted Value", ylab="Residuals") # 残差图boxcox(lm.sol, lambda=seq(0, 1, by=0.1) # 确定参数图#Box-Cox 变换后 , 作回归分析lambda<-0.55; Ylam<-(Ylambda-1)/lambda #B-C 变换lm.lam<-lm(YlamX); summary(lm.lam)plot(fitted(lm.lam), resid(lm.lam), cex=1.2, pch=21, col="red", bg="orange", xlab="Fitted Value", ylab="New Residuals") # 变换后的残差图beta0<-lm.lam$coefficients1 # 回归系数beta1<-lm.lam$coefficients2curve(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-Cox Transformations", 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(y1+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) # 定义权重向量weights18<-0.5 #对18号样本点的权系数重新赋值为0.5其它的均为1,以减少 18 号的影响lm.correct<-lm(y1+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(xc(19, 18), yc(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", "Correct Reg"), 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(y1+x1+x2+x3, data=france)lm.sum<-summary(lm.sol)coef(lm.sum) # 提取模型系数# 计算矩阵 X T X 特征值的最小值X<-as.matrix(france1:3);Xmin(eigen(cor(X)$values) #cor()用来得到相关系数矩阵，eigen()计算特征值与特征向量#9 例 7.9 ：对例 7.8 的法国经济分析数据作岭回归分析。library(MASS)lm.rid<-lm.ridge(y1+x1+x2+x3, data=france, lambda=seq(0, 0.2, length=50)par(mai=c(0.9, 0.9, 0.2, 0.2)plot(lm.rid) # 绘制岭迹图# 选择岭参数，作岭回归分析lm.ridge(y1+x1+x2+x3, data=france, lambda=0.04)#10 例 7.10 的求解：rt<-data.frame(x=c(20,19.6,19.6,19.4,18.4,19, 19,18.3,18.2,18.6,19.2,18.2,18.7,18.5,18,17.4,16.5,17.2,17.3, 17.8,17.3,18.4,16.9), y=c(1.0,1.2,1.1,1.4,2.3,1.7,1.7,2.4,2.1, 2.1,1.2,2.3,1.9,2.4,2.6,2.9,4.0,3.3,3.0,3.4,2.9,1.9,3.9) # 数据lm.ls<-lm(y1+x, data=rt); summary(lm.ls) # 线性回归par(mai=c(0.9, 0.9, 0.3, 0.2) # 绘图参数plot(rt, cex=1.2, col=2, pch=19, xlim=c(10, 20), ylim=c(1,4)abline(lm.ls, lwd=2, col=1) # 绘散点和回归图points(rt$x15, rt$y15, cex=1.2, col=4, lwd=2) # 标 15 号点text (rt$x15, rt$y15, labels="15", adj=c(-0.5, 0)rt$x15<-10 # 给 15 号点重新赋值points(rt$x15, rt$y15, pch=21, cex=1.2, col=4, lwd=2, bg="red") # 标记新的 15 号点text (rt$x15, rt$y15, labels="15", adj=c(-0.5, 0)lm.ls<-lm(y1+x, data=rt); summary(lm.ls) # 再线性回归abline(lm.ls, lwd=2, col=4, lty=2) # 再绘回归图#11 例 7.11 ：用最小一乘估计对例 7.10 中的数据 ( 原始数据和带有错误的数据 ) 作最小一乘估计。lm.ls<-lm(y1+x, data=rt) # 线性回归( 最小二乘)Q<-function(beta, data) sum(abs(data,2-beta1-beta2*data,1) #绝对值函数绝对值函数|yi-(0+1*xi)| par(mai=c(0.9, 0.9, 0.3, 0.2)plot(rt, cex=1.2, col=2, pch=19, xlim=c(10, 20), ylim=c(1,4) # 散点图z<-optim(lm.ls$coefficients, Q, data=rt); z$par #最小一乘估计( 原数据)abline(z$par, lwd=2, col=1, lty=1) # 绘回归直线图points(rt$x15, rt$y15, cex=1.2, col=4, lwd=2) #标记15 号点text (rt$x15, rt$y15, labels="15", adj=c(-0.5, 0)rt$x15<-10 # 给15 号点重新赋值points(rt$x15, rt$y15, pch=21, cex=1.2, col=4, lwd=2, bg="red")text (rt$x15, rt$y15, labels="15", adj=c(-0.5, 0)z<-optim(lm.ls$coefficients, Q, data=rt); z$par #最小一乘估计( 新数据)abline(z$par, lwd=2, col=4, lty=5) #再绘回归直线图#12 例 7.12 ：用 rlm() 对例 7.10 中带有错误的数据作 M 估计。lm.ls<-lm(y1+x, data=rt) # 线性回归( 最小二乘)par(mai=c(0.9, 0.9, 0.3, 0.2) # 图形参数plot(rt, cex=1.2, col=2, pch=19, xlim=c(10, 20), ylim=c(1,4)abline(lm.ls, lwd=2, col=1, lty=1) # 画原始数据最小二乘估计图形points(rt$x15, rt$y15, cex=1.2, col=4, lwd=2) # 标记15 号点text (rt$x15, rt$y15, labels="15", adj=c(-0.5, 0)rt$x15<-10 # 给15 号点重新赋值（错误数据）points(rt$x15, rt$y15, pch=21, cex=1.2, col=4, lwd=2, bg="red")text (rt$x15, rt$y15, labels="15", adj=c(-0.5, 0)library(MASS) # 加载包lm.ro1<-rlm(y1+x, data=rt, method = "M", maxit=50) # M 估计lm.ro1; abline(lm.ro1, lwd=3, col=4, lty=4) # 画M 估计直线lm.ro2<-rlm(y1+x, data=rt, method = "MM", maxit=50) # MM 估计lm.ro2; abline(lm.ro2, lwd=3, col=2, lty=5) # 画MM#13 例 7.13 ：用 lqs() 对例 7.10 中带错误的数据作抗干扰回归。lm.ols<-lm(yx, data=rt); lm.ols # 线性回归( 最小二乘)par(mai=c(0.9, 0.9, 0.3, 0.2)plot(rt, cex=1.2, col=2, pch=19, xlim=c(10, 20), ylim=c(1,4)points(rt$x15, rt$y15, cex=1.2, col=4, lwd=2) # 标记15 号点text (rt$x15, rt$y15, labels="15", adj=c(-0.5, 0)rt$x15<-10 # 改动数据points(rt$x15, rt$y15, pch=21, cex=1.2, col=4, lwd=2, bg="red")text (rt$x15, rt$y15, labels="15", adj=c(-0.5, 0)library(MASS) # 加载包lqs.lts<-lqs(yx, data=rt, method="lts"); lqs.lts # lqs 回归( 多种方法)lqs.lqs<-lqs(yx, data=rt, method="lqs"); lqs.lqslqs.lms<-lqs(yx, data=rt, method="lms"); lqs.lmslqs.S<-lqs(yx, data=rt, method="S"); lqs.Sabline(lm.ols, lwd=2, col=1) # 画线abline(lqs.lts, lwd=2, lty=2, col="red")abline(lqs.lqs, lwd=2, lty=3, col="blue")abline(lqs.lms, lwd=2, lty=4, col="brown")abline(lqs.S, lwd=2, lty=5, col="green")#14 例 7.14 的求解：cement<-data.frame( X1=c( 7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10), X2=c(26, 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68), X3=c( 6, 15, 8, 8, 6, 9, 17, 22, 18, 4, 23, 9, 8), X4=c(60, 52, 20, 47, 33, 22, 6, 44, 22, 26, 34, 12, 12), Y =c(78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5, 93.1,115.9, 83.8, 113.3, 109.4)lm.sol<-lm(Y X1+X2+X3+X4, data=cement) #线性回归summary(lm.sol)# 逐步回归 ( 默认自动增减变量 )lm.step<-step(lm.sol)summary(lm.step)# 人工方法去掉X 4lm.opt<-lm(Y X1+X2, data=cement)summary(lm.opt)#例 7.14 的另一种求解方法：# 从增加变量的角度逐步回归lm0<-lm(Y1, data=cement) # 从 Y 中不含自变量开始lm.ste<-step(lm0, scope = X1+X2+X3+X4, k=4) # 回归时的搜索范围为全部变量 X 1 , X 2 , X 3 , X 4#15 例 7.15 的求解：# 输入数据，由散点图 ( 略 ) 知，可进行二次多项式回归x <- c(37.0, 37.5, 38.0, 38.5, 39.0, 39.5, 40.0, 40.5, 41.0, 41.5, 42.0, 42.5, 43.0)y <- c(3.40, 3.00, 3.00, 3.27, 2.10, 1.83, 1.53, 1.70, 1.80, 1.90, 2.35, 2.54, 2.90) # 输入数据summary(lm.sol<-lm(y1+x+I(x2) # 二次变二元回归# 预测x 0 处的值predict(lm.sol, data.frame(x=38.75), interval="p")# 计算预测区间和置信区间 , 并绘图new <- data.frame(x = seq(37, 43, by=0.1)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(40, 4.5, c("Points", "Fitted", "Prediction", "Confidence"), pch=c(19, NA, NA, NA), lty=c(NA, 1,5,2), col=c("orange", "blue", "red", "brown") # 添加图例#16 例 7.16 的求解：cars.lo <- loess(dist speed, cars, degree = 1, control =loess.control(surface = "direct"); summary(cars.lo)#degree=1 表示局部线性 , control 为让预测值向两端延伸newdata<-data.frame(speed = seq(2, 28, 1)cars.pr <- predict(cars.lo, newdata, se = TRUE )# 取 se=T, 表示除得到预测值 $fit 外 , 还可得到预测的标准差 $se.fitpar(mai=c(0.9, 0.9, 0.5, 0.1)plot(cars, cex=1.2, pch=19, col="magenta", main="Stopping Distance versus Speed", xlab="Speed (mph)", ylab="Distance (ft)", xlim=c(2,28), ylim=c(0, 120)lines(newdata$speed, cars.pr$fit, lty=1, col="blue", lwd=2)# 绘制约 95%(2 原则 ) 的置信区间曲线lines(newdata$speed, cars.pr$fit - 2*cars.pr$se.fit, lty=2, col="red", lwd=2)lines(newdata$speed, cars.pr$fit + 2*cars.pr$se.fit, lty=2, col="red", lwd=2)#17 例 7.17 的求解：dfr<-data.frame(x = c(129.0, 140.0, 108.5, 88.0, 185.5, 195.0, 105.5, 157.5, 107.5, 77.0, 81.0, 162.0, 162.0, 117.5), y = c(7.5, 141.5, 28.0, 147.0, 22.5, 137.5, 85.5, -6.5, -81.5, 3.0, 56.5, -66.5, 84.0, -38.5), z = -c(1.22, 2.44, 1.83, 2.44, 1.83, 2.44, 2.44, 2.74, 2.74, 2.44, 2.44, 2.74, 1.22, 2.74)dfr.loess <- loess(zx+y, dfr, degree=2, span=0.75) # 二次summary(dfr.loess)dfr.mar <- list(x=seq(min(dfr$x), max(dfr$x), 5), y=seq(min(dfr$y), max(dfr$y), 5)dfr.lo <- predict(dfr.loess, expand.grid(dfr.mar)#expand.grid() 表示将向量生成网络形式的数据框，得到的dfr.lo 为一个 length(x) × length(y) 的矩阵par(mai = c(0.9,0.9,0.3,0.2)contour(dfr.mar$x, dfr.mar$y, dfr.lo, xlab = "X", ylab = "Y", levels = seq(-1.1, -2.9, -0.1), cex = 0.7) # 等值线图#18 例 7.18 的求解：# 非线性回归分析cl<-read.table("F:/文档/大学课程/R语言/ch07/time_cl.dat", header=T) # 输数据nls.sol<-nls(ya+(0.49-a)*exp(-b*(x-8), data=cl, start=list(a=0.1, b=0.01) # 非线性回归summary(nls.sol) # 提取回归结果# 绘制散点图和回归预测曲线newdata<-data.frame(x=seq(8, 42, by=1)nls.pre<-predict(nls.sol, newdata)par(mai=c(0.9, 0.9, 0.2, 0.1)plot(cl, pch=19, col=2, xlab="X", ylab="Y")lines(newdata$x, nls.pre, lwd=2, col=4)#19 # 输入数据 , 作 Logistic 回归x <- 0 : 5 #6 种电流强度y <- c(0, 9, 21, 47, 60, 63) # 对应 6 种电流的响应次数Y <- cbind(y, 70-y) # 构造 2 列矩阵 , 表示成功和失败的次数glm.sol <- glm(Yx, family=binomial(link=logit) # 使用二项分布族 , 连接函数为 logit( 默认 )summary(glm.sol)# 预测（电流为 3.5mA 时，有响应的牛的概率）predict(glm.sol, data.frame(x=3.5), type = "response")# 控制（若有 50% 的牛有响应，电流强度为多少）c<- -glm.sol$coefficients1/glm.sol$coefficients2c # 当 P=0.5 时 , ln(P/(1-P)=0, 从而 x = - 0 / 1# 绘制响应的比例与回归曲线d<- seq(0, 5, length=100)names<-c("logit", "probit", "cauchit") # 不同连接函数col <- c("black", "blue", "brown")lty<-c(1, 2, 4)par(mai=c(0.9,0.9, 0.2,0.1)plot(x, y/70, pch=19, cex=1.2, col=2, xlab="X", ylab="proibability", ylim=c(0,.95)for (i in 1:3) glm.sol <- glm(Yx, family=binomial(link=namesi) pre<-predict(glm.sol, data.frame(x=d), type = "response") lines(d, pre, lwd=2, col=coli, lty=ltyi) # 绘制不同曲线legend(0.2, 0.9, legend=c("points", names), col=c("red", col), pch=c(19, NA, NA, NA), lty=c(NA, lty), lwd=2)#20 例 7.20 的求解：life<-read.table("F:/文档/大学课程/R语言/ch07/leukemia_life.dat",header = T)glm.sol<-glm(YX1+X2+X3, family=binomial, data=life)summary(g