薛毅 数学模型 数学建模 第八次作业 多元分析实验.docx

数学模型 第八次作业 多元分析实验8.1实验目的与要求l 学会对数据进行线性回归分析、预测与回归诊断l 学会对数据进行方差分析和判别分析l 建立相应的统计模型，用R软件计算，并对计算结果进行分析和讨论8.2 基本实验1. 回归分析为估计山上积雪融化后对下游灌溉的影响，在山上建立一个观测站，测量最大积雪深度X（米）与当年灌溉面积Y（公顷），测得连续10年的数据如表8.1所示。（1）建立一元线性回归模型，求解，并验证系数、方程或相关系数是否通过检验；（2）现测得今年的数据是X=7米，给出今年灌溉面积的预测值、预测区间和置信区间(=0.05)；（3）将数据散点、回归预测值、回归预测区间和置信区间均匀化在一张图上，分析线性回归的拟合情况。解：(1) 为了研究这些数据中所蕴含的规律性，根据10对数据利用R软件作出散点图，编程如下： x y plot(x,y, xlab=X, ylab=Y, cex=1.4, pch=19, col=red)得到如下图像：分析图像，数据点大致落在一条直线附近，说明变量x和y之间大致可看作线性关系，假定有如下结构式：y=0+1x+其中0和1是未知常数，为回归系数，为其它随机因素对灌溉面积的影响，服从正态分布N(0,2)。利用R软件进行一元线性回归分析，并提取相应的计算结果： x y lm.sol 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 freedomMultiple R-squared: 0.9781, Adjusted R-squared: 0.9754 F-statistic: 357.5 on 1 and 8 DF, p-value: 6.33e-08Estimate项中给出了回归方程的系数估计，即0=140.95；1=364.18观查其中的评价参数易知对于0项的估计并不是很准确，不显著。但该方程总体通过了F统计数的检验，其p值为6.33e-08 new predict(lm.sol,new,+ interval=prediction,+ level=0.95)得到如下结果：fit lwr upr1 2690.227 2454.971 2925.484得到灌溉面积的预测值为2690.227、预测区间2454.971和置信区间(=0.05)为2925.484。(3) 利用R软件做出图像并保存，编程如下：先重复回归线性分析： x y plot(x,y, xlab=X, ylab=Y, cex=1.4, pch=19, col=red) lm.sol summary(lm.sol)做出图像： abline(lm.sol, lwd=2, col=blue) segments(x, fitted(lm.sol), x, y, lwd=2, col=blue)标注图像： ex1 ex2 points(x8, fitted(lm.sol)8, pch=19, cex=1.4, col=blue) text(c(5.7, 5.7), c(2400, 2100), labels = c(ex1, ex2)保存图像： savePlot(regression, type=eps)最终得到的图像如图所示：由图像可以直观看出此线性回归的拟合对于前4年的拟合误差比较大，误差最大的是第2年。对于后6年的拟合是比较吻合的。2. 回归诊断对1题得到的回归模型作回归诊断：（1）残差是否满足齐性、正态性的条件；（2）哪些点可能是异常值点；（3）如果有异常值点，则去掉异常值点，再作回归分析。解：使用1题得到的回归模型： 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 lm.sol summary(lm.sol)结果同1题，调用influence.measures()函数进行诊断： influence.measures(lm.sol)得到回归诊断结果：Influence measures of lm(formula = y 1 + x, data = intellect) : 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 做出回归诊断图并保存：op weights-rep(1, n); weights2 lm.correct-lm(y1+x, data=intellect, subset=-7,+ weights=weights)在程序中，subset=-7是去掉7号点。weights-rep(1,n)是将所有点的权赋为1，weights2 summary(lm.correct)得到如下结果：Call:lm(formula = y 1 + x, data = intellect, subset = -7, weights = weights)Weighted Residuals: Min 1Q Median 3Q Max -93.962 -60.875 -9.213 36.037 115.036 Coefficients: Estimate Std. Error t value Pr(|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 1 Residual 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验证如下，做出两次线性回归的对比图：attach(intellect)plot(x, y, cex=1.2, pch=21, col=red, bg=orange)abline(lm.sol, col=blue, lwd=2)text(xc(7, 2), yc(7, 2), label=c(7, 2), adj=c(1.5, 0.3)detach()abline(lm.correct, col=red, lwd=2, lty=5)legend(5, 3200, c(Points, Regression, Correct Reg), pch=c(19, NA, NA), lty=c(NA, 1,5), col=c(orange, blue, red)savePlot(diagnoses2, type=eps)图中实线是原始数据计算的结果，虚线是修正数据后的计算结果。再次进行回归诊断： op plot(lm.correct, 1:4) par(op) savePlot(diagnoses3, type=eps)所有结果具有改善。3. 回归分析和逐步回归研究同一地区土壤所含可给态磷（Y）的情况，得到18组数据如表8.2所示。表中X1为土壤内所含无机磷浓度，X2为土壤内溶于K2CO3溶液并受溴化物水解的有机磷，X3为土壤内溶于K2CO3溶液但不溶于溴化物水解的有机磷。（1） 建立多元线性回归方程模型，求解，并验证系数、方程或相关系数是否通过检验；（2）做逐步回归分析。解：(1)首先根据题意建立多元线性回归方程：Y=0+1X1+2X2+3X3+利用R软件进行求解，使用lm()函数，用函数summary()提取信息，写出R程序： import lm.sol 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 freedomMultiple 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.16102X3p-值为0.009227 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，则相应的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 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值和残差平方值上升都是最小的，因此去掉X3项再次做线性回归： lm.opt summary(lm.opt)得到结果如下：Call:lm(formula = Y X1, data = import)Residuals: Min 1Q Median 3Q Max -31.486 -8.282 -1.674 5.623 59.337 Coefficients: Estimate Std. Error t value Pr(|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 1 Residual 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此时常数项的检测结果为极为显著，X1项系数为很显著。方程式P-值为0.001417 mouse mouse.lm anova(mouse.lm)得到如下结果：Analysis of Variance TableResponse: X Df Sum Sq Mean Sq F value Pr(F) A 2 73.118 36.559 9.104 0.001422 *Residuals 21 84.329 4.016 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.从结果中看到P-值为0.001422 attach(mouse) tapply(X, A, mean)得到如下结果： 1 2 3 1.60625 3.10750 5.82375可以看出不同饲料喂食下的小鼠肝中铁含量的均值已有明显差异。再做多重t检测： 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由计算结果得出结论，1与3、2与3是有显著差异的，而1与2没有显著差异。即是说，喂食饲料A和喂食饲料B情况下小鼠肝中铁含量有显著差异；喂食饲料B和喂食饲料C情况下小鼠肝中铁含量有显著差异；喂食饲料A和喂食饲料B情况下小鼠肝中铁含量无显著差异。进一步，使用plot()函数画出线箱图并保存： plot(XA, col=5:7, + main=Box-and-Whisker Plot of Mouse Data) detach(mouse) savePlot(box_plot, type=eps)可以直观看到数据的水平及各因素之间的差异。(3) 根据题意，先编写程序，做Shapiro-Wilk正态性检验 attach(mouse) tapply(X,A,shapiro.test)得到如下结果：$1 Shapiro-Wilk normality testdata: X1L W = 0.8742, p-value = 0.1656$2 Shapiro-Wilk normality testdata: X2L W = 0.8893, p-value = 0.2306$3 Shapiro-Wilk normality testdata: X3L W = 0.985, p-value = 0.9833结果显示三组数据均数据满足正态性。再用Bartlett函数做方差齐性检验： attach(mouse) bartlett.test(X, A)得到如下结果：Bartlett test of homogeneity of variancesdata: X and A Bartletts K-squared = 10.5677, df = 2, p-value = 0.005073从结果中看到p-值为0.005073 oneway.test(XA, data=mouse)得到方差的分析结果：One-way analysis of means (not assuming equal variances)data: X and A F = 10.3592, num df = 2.00, denom df = 10.51, p-value = 0.003271此时P-值较第一问计算时的结果有所增大，但是仍然满足p-值 tree.aov factory factory.aov summary(factory.aov)得到结果：Df Sum Sq Mean Sq F value Pr(F) A 2 5.408 2.704 6.130 0.02090 * B 2 7.841 3.921 8.888 0.00740 *A:B 4 12.192 3.048 6.910 0.00793 *Residuals 9 3.970 0.441 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1结果显示有交互作用时在显著水平0.05下，因素A（反应温度）效应显著，而因素B（反应压力）和交互效应很显著。（2） 使用R软件求解各种反应温度下产量均值的估计： attach(factory) tapply(Y, A, mean)得到如下结果： 1 2 3 4.650000 4.533333 5.750000计算各种反应压力下产量均值的估计： tapply(Y, B, mean)得到如下结果：1 2 3 5.783333 4.166667 4.983333计算同时考虑温度和压力下产量均值的估计： matrix(tapply(Y, A:B, mean), nr=3, nc=3, byrow=T,+ dimnames=list(levels(A), levels(B)得到如下结果： 1 2 31 4.45 5.0 4.502 6.30 3.6 3.703 6.60 3.9 6.75(3) 从第一问结果因素A（反应温度）效应显著、而因素B（反应压力）和交互效应很显著来看第二问得到的数据即可得到答案。各种反应温度下产量均值中3条件下最多；各种反应压力下产量均值中1条件下最多；交互效应下3、3条件的产量均值最多，且高于单独作用时的产量均值；综合看来，选用3、3条件是最佳的，即采用80的反应温度3公斤的反应压力时对生产最有利。8.3 加分实验（早餐方便粥数据分析）某超市销售三家厂商43种品牌的早餐方便粥，表8.8列出了这些方便粥的具体数据。你的任务是分析这些数据，是否有理由认为某一厂家的产品比其他厂家的产品更“有营养”（高蛋白、低脂肪、高纤维、低糖等）。解：根据题意，明确解题思路，要解决的问题是：是否有理由认为某一厂家的产品比其他厂家的产品更“有营养”（高蛋白、低脂肪、高纤维、低糖等）？也就是研究营养成分在不同厂家之间是否有显著性差异。营养成分数据都是定量数据，因此可以采用方差分析的思想来解决这个问题。为了数据表示的方便，我们将厂家A、B、C分别用数字1、2、3来表示。由于数据量比较大，解答过程用SPSS软件进行计算，而没有选用R软件。分析过程的显著性统一设定为0.05.解答过程：1. 先对数据做方差齐性检验，计算结果如下表所示：方差齐性检验Levene 统计量df1df2显著性热量6.665240.003蛋白质1.676240.200脂肪6.045240.005钠7.146240.002纤维2.428240.101碳水化合物.917240.408糖.729240.489钾3.266240.049由上表可以看出，在0.05的显著性水平下，热量、脂肪、钠、钾三个变量没有通过方差齐性检验，其它都是方差齐性的。因此对热量、脂肪、钠、钾三个变量做方差非齐性的方差分析，其余变量做方差齐性的方差分析模型。2. 方差分析（1）方差齐性变量的方差分析结果：ANOVA平方和df均方F显著性蛋白质组间.6822.341.220.804组内62.016401.550总数62.69842纤维组间10.88425.4421.740.189组内125.088403.127总数135.97242碳水化合物组间130.318265.1594.131.023组内630.8684015.772总数761.18642糖组间47.564223.7821.165.322组内816.7154020.418总数864.27942从结果可以看出，在0.05的显著性水平下，三个厂商在碳水化合物上有显著性差异，其余变量没有显著性差异。下面进一步进行两两比较分析，看不同厂商的差异程度，如下表所示：多重比较因变量(I) 厂商(J) 厂商均值差 (I-J)标准误显著性95% 置信区间下限上限碳水化合物LSD12-.66181.3101.616-3.3101.98634.5882*1.8858.020.7778.40021.66181.3101.616-1.9863.31035.2500*1.8486.0071.5148.98631-4.5882*1.8858.020-8.400-.7772-5.2500*1.8486.007-8.986-1.514*. 均值差的显著性水平为 0.05。从上表看出：碳水化合物，厂商1和厂商3有显著性差异，厂商2与厂商3有显著性差异，厂商1和厂商2没有显著性差异。其均值图为：（2）方差非齐性变量的方差分析结果：针对热量、脂肪、钠、钾三个变量。主体间效应的检验源因变量III 型平方和df均方FSig.校正模型热量2237.510a21118.7553.476.041脂肪4.035b22.0173.517.039钾5163.382c22581.691.579.565钠50024.129d225012.0654.686.015截距热量352416.2741352416.2741094.961.000脂肪37.604137.60465.563.000钾200603.0501200603.05044.979.000钠862181.1181862181.118161.520.000厂商热量2237.51021118.7553.476.041脂肪4.03522.0173.517.039钾5163.38222581.691.579.565钠50024.129225012.0654.686.015误差热量12874.11840321.853脂肪22.94240.574钾178397.083404459.927钠213516.569405337.914总计热量515800.00043脂肪68.00043钾490000.00043钠1663950.00043校正的总计热量15111.62842脂肪26.97742钾183560.46542钠263540.69842a. R 方 = .148（调整 R 方 = .105）b. R 方 = .150（调整 R 方 = .107）c. R 方 = .028（调整 R 方 = -.020）d. R 方 = .190（调整 R 方 = .149）从上表可以看出，在0.05的显著性水平下，不同厂商的热量、脂肪