实验7 多元分析实验 1. 回归分析 解:(1) 根据题意,对数据利用R软件作出散点图 > x<-c(5.1,3.5,7.1,6.2,8.8,7.8,4.5,5.6,8.0,6.4)
> y<-c(1907,1287,2700,2373,3260,3000,1947,2273,3113,2493) > plot(x,y, xlab=\得到如下图像:
Y1500200025003000456X789
分析图像,数据点大致落在一条直线附近,说明变量x和y之间大致可看作线性关系,假定有如下结构式: y=β0+β1x+ε
其中β0和β1是未知常数,为回归系数,ε为其它随机因素对灌溉面积的影响,ε服从正态分布N(0,σ2)。
利用R软件进行一元线性回归分析,并提取相应的计算结果: > x<-c(5.1,3.5,7.1,6.2,8.8,7.8,4.5,5.6,8.0,6.4)
> y<-c(1907,1287,2700,2373,3260,3000,1947,2273,3113,2493) > lm.sol<-lm(y ~ 1+x) > summary(lm.sol) 得到如下结果: Call:
lm(formula = y ~ 1 + x)
Residuals:
Min 1Q Median 3Q Max
-128.591 -70.978 -3.727 49.263 167.228
Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) 140.95 125.11 1.127 0.293 x 364.18 19.26 18.908 6.33e-08 *** ---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 96.42 on 8 degrees of freedom
Multiple R-squared: 0.9781, Adjusted R-squared: 0.9754 F-statistic: 357.5 on 1 and 8 DF, p-value: 6.33e-08
Estimate项中给出了回归方程的系数估计,即β0=140.95;β1=364.18 观查其中的评价参数易知对于β0项的估计并不是很准确,不显著。 但该方程总体通过了F统计数的检验,其p值为6.33e-08<0.05 由此得到的回归方程为:Y=140.95+364.18X
(2)若现测得今年的数据是X=7米,则有X=X0=7,置信水平为0.95,此时利用R软件求解,编程如下:
> new<-data.frame(x=7) > predict(lm.sol,new,
+ interval=\+ level=0.95) 得到如下结果:
fit lwr upr
1 2690.227 2454.971 2925.484
得到灌溉面积的预测值为2690.227、预测区间2454.971和置信区间(α=0.05)为2925.484。
(3)利用R软件做出图像并保存,编程如下: 先重复回归线性分析:
> x<-c(5.1,3.5,7.1,6.2,8.8,7.8,4.5,5.6,8.0,6.4)
> y<-c(1907,1287,2700,2373,3260,3000,1947,2273,3113,2493) > plot(x,y, xlab=\>
> lm.sol<-lm(y ~ 1+x) > summary(lm.sol) 做出图像:
> abline(lm.sol, lwd=2, col=\
> segments(x, fitted(lm.sol), x, y, lwd=2, col=\标注图像:
> ex1<-expression(paste(\> ex2<-expression(paste(\>
> points(x[8], fitted(lm.sol)[8], pch=19, cex=1.4, col=\
> text(c(5.7, 5.7), c(2400, 2100), labels = c(ex1, ex2)) 保存图像:
> savePlot(\
最终得到的图像如图所示:
25003000(xi,yi)^)(xi,yiY15002000456X789
由图像可以直观看出此线性回归的拟合对于前4年的拟合误差比较大,误差最大的是第2年。对于后6年的拟合是比较吻合的。
2. 回归分析和逐步回归 解:
(1)首先根据题意建立多元线性回归方程:
Y=β0+β1X1+β2X2+β3X3+ε
利用R软件进行求解,使用lm()函数,用函数summary()提取信息,写出R程序:
> import<-data.frame(
+ X1=c(0.4,0.4,3.1,0.6,4.7,1.7,9.4,10.1,11.6,12.6,10.9,23.1,23.1,21.6,23.1,1.9,26.8,29.9), + X2=c(52,23,19,34,24,65,44,31,29,58,37,46,50,44,56,36,58,51),
+ X3=c(158,163,37,157,59,123,46,117,173,112,111,114,134,73,168,143,202,124), + Y= c(64,60,71,61,54,77,81,93,93,51,76,96,77,93,95,54,168,99) + )
> lm.sol<-lm(Y~X1+X2+X3, data=import) > summary(lm.sol) 得到如下结果: Call:
lm(formula = Y ~ X1 + X2 + X3, data = import)
Residuals:
Min 1Q Median 3Q Max -28.349 -11.383 -2.659 12.095 48.807
Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) 43.65007 18.05442 2.418 0.02984 * X1 1.78534 0.53977 3.308 0.00518 ** X2 -0.08329 0.42037 -0.198 0.84579 X3 0.16102 0.11158 1.443 0.17098 ---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 19.97 on 14 degrees of freedom
Multiple R-squared: 0.5493, Adjusted R-squared: 0.4527 F-statistic: 5.688 on 3 and 14 DF, p-value: 0.009227
所以得到回归方程为:
Y=43.65007 +1.78534X1 -0.08329X2+0.16102X3
p-值为0.009227<0.05方程本身是通过检测的,各项系数的检验结果为:
常数项显著;X1项系数很显著;X2项系数不显著;X3项系数不显著。有两项系数没有通过检验,总体来说拟合并不理想。(2) 利用R软件进行逐步回归:
> lm.step<-step(lm.sol)
得到如下结果: Start: AIC=111.27 Y ~ X1 + X2 + X3
Df Sum of Sq RSS AIC - X2 1 15.7 5599.4 109.32
Step: AIC=109.32 Y ~ X1 + X3
Df Sum of Sq RSS AIC
从程序的运行结果可以看到,用全部变量作回归方程时,AIC值为111.27。如果去
掉变量X2,则相应的AIC值为109.32;如果去掉变量X3则相应的AIC值为111.77;如果去掉变量X1则相应的AIC值为119.66。软件去掉X2项,进入下一轮运算,给出结果:
> summary(lm.step) 得到运算结果: Call:
lm(formula = Y ~ X1 + X3, data = import)
Residuals:
Min 1Q Median 3Q Max -29.713 -11.324 -2.953 11.286 48.679
Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) 41.4794 13.8834 2.988 0.00920 ** X1 1.7374 0.4669 3.721 0.00205 ** X3 0.1548 0.1036 1.494 0.15592 ---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 19.32 on 15 degrees of freedom
Multiple R-squared: 0.5481, Adjusted R-squared: 0.4878 F-statistic: 9.095 on 2 and 15 DF, p-value: 0.002589
此时回归系数检验的水平已有显著提升,但X3项系数仍然不显著。 利用drop1()函数计算: > drop1(lm.step)
得到如下结果: Single term deletions
Model: Y ~ X1 + X3
Df Sum of Sq RSS AIC
此时的结果说明,去掉X3项的时候,AIC值和残差平方值上升都是最小的,因此去掉X3项再次做线性回归:
> lm.opt<-lm(Y~X1,data=import); > summary(lm.opt) 得到结果如下: Call:
lm(formula = Y ~ X1, data = import)
Residuals: