北京工业大学-数学建模8-多元分析实验201312.doc

数学建模作业8多元分析实验一、 基本实验1.回归分析为估计山上积雪融化后对下游灌溉的影响，在山上建立一个观测站，测量最大积雪深度X(米)与当年灌溉面积Y(公顷)，测得连续10年的数据如表8.1所示。表8.1 10年中最大积雪深度与当年灌溉面积的数据XYXY15.1190767.8300023.5128774.5194737.1270085.6227346.2237398.0311358.83260106.42493（1）建立一元线性回归模型，求解，并验证系数、方程或相关系数是否通过检验；（2）观测得今年的数据是X=7米，给出今年灌溉面积的预测值、预测区间和置信区间（=0.05）；（3）将数据散点、回归预测值、回归的预测区间和置信区间均画在一张图上，分析线性回归的拟合情况。解：（1）利用R软件中的lm()求回归参数并作相应的检验。写出相应的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|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 1Residual standard error: 96.42 on 8 degrees of freedomMultiple R-squared: 0.9781, Adjusted R-squared: 0.9754 F-statistic: 357.5 on 1 and 8 DF, p-value: 6.33e-08系数检验：常数项截距Intercept 的p值0.05,接受0假设，说明不显著。x的系数的p值0.05,拒绝0假设，x的系数极为显著，有*标示。方程检验：F-statistic的p值为6.33e-080.05,拒绝0假设，通过回归方程的检验。相关系数检验：R2=0.9781,调整的R2=0.9754, 说明因变量与自变量是相关的。除了常数项截距，回归方程通过了回归参数的检验与回归方程的检验，得到的回归方程为：（2）利用predict()函数进行预测，编写R程序：new predict(lm.sol, new, interval=confidence) fit lwr upr1 2690.227 2613.35 2767.105有分析结果可知：今年灌溉渠的预测值为2690.227公顷；预测区间是2613.35, 2767.105 (单位：公顷)；置信区间（=0.05）为2454.971, 2925.484 (单位：公顷)。（3）参考exam0806.R,编写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)new - data.frame(x = seq(3.5, 8.8, by=0.1)pp-predict(lm.sol, new, interval=prediction)pc-predict(lm.sol, new, interval=confidence)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, green, green), lwd=2)points(x,y, cex=1.4, pch=21, col=red, bg=orange)legend(4, 3500,c(Points, Fitted, Prediction, Confidence), pch=c(19, NA, NA, NA), lty=c(NA, 1,5,2), col=c(orange, blue, red, brown)savePlot(predict, type=eps)运行结果如图8.1：图中显示散点、回归预测值、回归的预测区间和置信区间。散点紧密的散布在回归直线两侧，只有一点落在置信区间外侧。散点分布与回归直线趋势一致，总体效果还可以。图8.1 数据的散点、回归预测值、预测区间与置信区间2.回归诊断对1题得到的回归模型做回归诊断：（1）残差是否满足齐性、正态性的条件；（2）那些点可能是异常值点；（3）如果有异常值点，则去掉异常值点，再作回归分析。解：利用函数plot()绘制各种残差的散点图，编写程序如下：调用influence.measures()函数并作回归诊断图，R程序如下：influence.measures(lm.sol)op - par(mfrow=c(2,2), mar=0.4+c(2,2,1,1),oma= c(2,2,0,0) plot(lm.sol, 1:4) par(op)Influence measures of lm(formula = y 1 + x) : dfb.1_ dfb.x dffit cov.r cook.d hat inf1 -0.3494 0.2706 -0.4479 1.165 0.09941 0.157 2 -1.6700 1.5077 -1.7320 0.859 1.06503 0.413 *3 0.0232 -0.0475 -0.1053 1.461 0.00627 0.126 4 -0.0271 0.0056 -0.0889 1.423 0.00447 0.100 5 0.5656 -0.6935 -0.8208 1.444 0.32649 0.349 6 -0.0473 0.0663 0.0964 1.594 0.00528 0.190 7 1.2525 -1.0577 1.4085 0.444 0.58059 0.229 *8 0.2329 -0.1531 0.3786 1.120 0.07117 0.120 9 -0.1884 0.2536 0.3465 1.474 0.06457 0.215 10 0.0133 0.0046 0.0730 1.432 0.00302 0.100图8.2 回归诊断图先进性回归诊断结果分析，得到的回归诊断结果共7列，第1,2列为dfbetas值（对应于常数和变量x）;第3列为DFFITS准则值；第4列为COVRATIO准则值；第5列为Cook距离；第6列为帽子值（高杠杆值）；第7列为影像店记号。有诊断结果可知2号和7号点是强影响点。下面分析回归诊断图8.2：第一张图是残差图，散点基本在联系周围，可以认为残差的方差满足齐性。第二张图是正态QQ图，除7号点外，基本都在一条直线上，也就是说除了7号点外，残差满足正态性。第三张图是标准差的平方根与预测值的散点图，7号值最大接近1.4，可能是异常值点。第四张图是Cook距离值，从图上来看2号距离最大，说明2号点可能是强影响点（高杠杆点）。综上：（1）残差满足齐性和正态性条件。（2）7号点可能是异常值点，2号点是强影响点（3）将2号点和7号点去掉，再作回归分析。处理强影响点：将7号点去掉 ，2号点赋值0.5的权重，加权计算减少它的影响。编写R程序如下：intellect-data.frame(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(y1+x, data=intellect)summary(lm.sol)n-length(intellect$x)weights-rep(1,n); weights2-0.5lm.correct|t|) (Intercept) 64.97 118.67 0.547 0.601 x 373.75 17.40 21.485 1.19e-07 *-Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 71.08 on 7 degrees of freedomMultiple R-squared: 0.9851, Adjusted R-squared: 0.9829 F-statistic: 461.6 on 1 and 7 DF, p-value: 1.192e-07与修正之前相比，R2与调整的R2有所提高。回归直线为：检验：修正后是否有效，看一下回归诊断的结果，编写R程序：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.3 修正后的回归诊断图，从四张图中可见残差齐次性、正态性，标准差的平方根与预测值，以及Cook距离等结果均有所改善。图8.3 修正后的回归诊断图3. 回归分析和逐步回归研究同一地区土壤所含可给态磷（Y）的情况，得到18组数据如表8.2所示。表中X1为土壤内所含无机磷浓度，X2为土壤内溶于K2CO3溶液并受溴化物水解的有机磷，X3为土壤内溶于K2CO3溶液但不溶于溴化物水解的有机磷。表8.2 某地区土壤所含可给态磷的情况X1X2X3YX1X2X3Y10.452158641012.6581125120.423163601110.9371117633.11937711223.1461149640.634157611323.1501347754.72459541421.644739361.765123771523.1561689579.4444681161.93614354810.131117931726.858202168911.629173931829.95112499（1）建立多元线性回归方程模型，求解，并验证系数、方程或相关系数是否通过检验；（2）作逐步回归分析。解：（1）输入数据，做多元线性回归：cement-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|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 1Residual standard error: 19.97 on 14 degrees of freedomMultiple R-squared: 0.5493, Adjusted R-squared: 0.4527 F-statistic: 5.688 on 3 and 14 DF, p-value: 0.009227由运行结果可知：常数项一般相关，X1系数高度相关，但是X2,X3系数的P值0.05,不相关，没有通过检验。方程F检验P值0.05,通过检验。相关系数R2,及调整R2在0.5左右，相关性一般。（2）下面用函数step()做逐步回归：lm.step-step(lm.sol)运行结果为：Start: AIC=111.27Y X1 + X2 + X3 Df Sum of Sq RSS AIC- X2 1 15.7 5599.4 109.32 5583.7 111.27- X3 1 830.6 6414.4 111.77- X1 1 4363.4 9947.2 119.66Step: AIC=109.32Y X1 + X3 Df Sum of Sq RSS AIC 5599.4 109.32- X3 1 833.2 6432.6 109.82- X1 1 5169.5 10768.9 119.09从程序的运行结果可见，当用全部变量作回归方程时，AIC值为111.27，接下来的数据表显示，如果去掉X2,则相应的ACI值为109.32，此时ACI值达到最小，因此R软件自动去掉变量X2,进行下一轮计算。R软件认为此时得到“最优”回归方程。下面我们分析一下回归结果，用函数summary()提取相关信息，看看是否是最优的回归方程:summary(lm.step)Call:lm(formula = Y X1 + X3, data = cement)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 1Residual standard error: 19.32 on 15 degrees of freedomMultiple 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 deletionsModel:Y X1 + X3 Df Sum of Sq RSS AIC 5599.4 109.32X1 1 5169.5 10768.9 119.09X3 1 833.2 6432.6 109.82从运算结果可见如果去掉X3,则AIC值从109,32增加到109.82，增加的最少。去掉X3,残差的平方和上升的也最少，拟合越好的方程，残差的平方和应越小，从这两项指标，应该再去掉变量X3.lm.opt|t|) (Intercept) 59.2590 7.4200 7.986 5.67e-07 *X1 1.8434 0.4789 3.849 0.00142 * -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 20.05 on 16 degrees of freedomMultiple R-squared: 0.4808, Adjusted R-squared: 0.4484 F-statistic: 14.82 on 1 and 16 DF, p-value: 0.001417最终结果还是满意的，尽管R2还是0.5左右。但是系数检验的显著性水平已经很好了。最后得到的“最优”回归方程为：4. 方差分析I(单因素方差分析)24只小鼠随机地分成3组，每组8只，每组喂养不同的饲料组，喂养一定时间后，测得小鼠肝中含铁量，结果如表8.3所示。表8.3 不同饲料组小鼠肝脏中铁含量（单位：ug/g）饲料12345678A1.001.011.131.141.702.012.232.63B0.961.231.541.962.943.685.596.96C2.073.724.504.906.006.848.2310.33（1）试用单因素方差分析方法分析不同饲料的小鼠肝中的铁含量是否有显著差异？（2）如果有显著差异，那些水平之间有显著差异？（3）数据是否满足正态性和方差齐性的要求，如果不满足请选择合适的分析方法，并给出你的最终结论。解：（1）利用方差分析表的方法进行分析：设小白鼠喂养的饲料为因素，三种不同的饲料为三个水平，喂养饲料后的小白鼠肝脏中铁含量视为来自三个正态分布总体（i=1,2,3）的样本观测值。问题归结为：。参考exam0813.R编写R程序：mouse-data.frame( X=c(1.00, 1.01, 1.13, 1.14, 1.70, 2.01, 2.23, 2.63, 0.96, 1.23, 1.54, 1.96, 2.94, 3.68, 5.59, 6.96, 2.07, 3.72, 4.50, 4.90, 6.00, 6.84, 8.23, 10.33), A=factor(rep(1:3, c(8,8,8)mouse.lmF) A 2 7 3.118 36.559 9.104 0.001422 *Residuals 21 84.329 4.016 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1由计算结果显示，P值0.05.因此拒绝原假设，即认为不同饲料的小鼠肝中铁含量有显著差异。（2）在（1）中F检验的结论是拒绝H0,下面进一步检验：首先计算各个因子间的均值，再用pairwise.t.test()作多重t检验。参考exam0814.R:求数据在各水平下的均值，编写R程序为：attach(mouse); tapply(X, A, mean)运行结果为： 1 2 3 1.60625 3.10750 5.82375再作多重t检验，用调整方法用缺省值，即Holm方法：R程序为：pairwise.t.test(X,A)运行结果为： Pairwise comparisons using t tests with pooled SD data: X and A 1 2 2 0.1489 - 3 0.0012 0.0262P value adjustment method: holm从计算结果可知，u1与u2,u1与u3均有显著差异，而u2与u3没有显著差异。即喂养A饲料与其他两种饲料小鼠肝中铁含量有显著差异，而喂养B饲料和C饲料小鼠肝中铁含量差异不显著。 （3）对模型的正态性和齐性进行检验：假设：H0: 三组数据满足正态分布。对数据做Shapiro-Wilk正态性检验：编写R程序：tapply(X,A, shapiro.test)运行结果为：$1 Shapiro-Wilk normality testdata: X1LW = 0.8742, p-value = 0.1656$2 Shapiro-Wilk normality testdata: X2LW = 0.8893, p-value = 0.2306$3 Shapiro-Wilk normality testdata: X3LW = 0.985, p-value = 0.9833有计算结果可知，三组数据的P值均大于0.05，接受原假设，即数据是满足正态分布的。下面做Bartlett方差齐性检验：假设: 编写R程序如下：bartlett.test(X,A)运行结果为： Bartlett test of homogeneity of variancesdata: X and ABartletts K-squared = 10.5677, df = 2, p-value = 0.005073由于P值=0.0050730.05,故拒绝原假设。说明数据不满足齐次性。当数据满足正态性，而不满足齐次性的时候，可以用函数oneway.test()作方差分析。假设：喂养饲料后的小白鼠肝脏中铁含量视为来自三个正态分布总体（i=1,2,3）的样本观测值。编写R程序：oneway.test(XA, data=mouse)运行结果为： One-way analysis of means (not assuming equal variances)data: X and AF = 10.3592, num df = 2.00, denom df = 10.51, p-value = 0.003271有分析结果P值0.05,拒绝原假设，即认为不同饲料的小鼠肝中铁含量有显著差异。与（1）的结论一致。该题最终结论：三组数据满足正态性，但不满足方差齐性。因此用oneway.test()作方差分析。通过分析喂养不同饲料的小鼠肝中铁含量有显著差异。其中，喂养A饲料与其他两种饲料小鼠肝中铁含量有显著差异，而喂养B饲料和C饲料小鼠肝中铁含量差异不显著。4. 方差分析II(双因素方差分析)为了提高化工厂的产品质量，需要寻求最优反应温度与反应压力的配合，为此选择如下水平：A: 反应温度（oC） 60 70 80B: 反应压力（公斤） 2 2.5 3在每个AiBj条件下做两次实验，其产量如表8.4所示。表8.4 实验数据A1A2A3B14.64.36.16.56.86.4B26.36.73.43.84.03.8B34.74.33.93.56.57.0（1）对数据作方差分析（考虑有交互作用的情况）；（2）计算各种反应温度下产量均值的估计，各种反应压力下产量均值的估计，以及同时考虑温度和压力下产量均值的估计；（3）通过（1）与（2）计算结果来说明，在今后的生产中，我们将如何选择生产的反应温度和反应压力，使得这些条件对生产最有利（注意，一定要说明你的理由）。解：（1）假定，各xijk相互独立，数据可以分解成：判断A,B及交互效应的影响是否显著等价于检验下列假设：输入数据，用aov()函数求解，用summary()函数列出方差分析信息。编写R程序如下：fanying-data.frame ( Y=c(4.6, 4.3, 6.1, 6.5, 6.8, 6.4, 6.3, 6.7, 3.4, 3.8, 4.0, 3.8, 4.7, 4.3, 3.9, 3.5, 6.5, 7.0), B=gl(3, 6, 18, labels= paste(B, 1:3, sep=), A=gl(3, 2, 18, labe