数学建模作业 实验7多元分析实验.docx

1、数学建模作业 实验7多元分析实验实验7 多元分析实验1. 回归分析解：（1） 根据题意，对数据利用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

2、 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

3、Multiple 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.97

4、1 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

5、)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）首先根据题意建立多元线性回归方程：Y=0+1X1+2X2+3X3+利用R软件进行求解，使用lm()函数，用函数summary()提取信息，写出R程序： import lm.sol summa

6、ry(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 X

7、3 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.

8、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从程序的运行结果

9、可以看到，用全部变量作回归方程时，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

10、 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

11、 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(

12、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.0

13、5 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

14、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 an

15、d 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_pl

16、ot, 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,

17、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 assumin

18、g 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

19、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

20、.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条件的产量均值最多，且高于单独作用时的产量均值；

21、综合看来，选用3、3条件是最佳的，即采用80的反应温度3公斤的反应压力时对生产最有利。7.3 加分实验解：根据题意，明确解题思路，要解决的问题是：是否有理由认为某一厂家的产品比其他厂家的产品更“有营养”（高蛋白、低脂肪、高纤维、低糖等）？也就是研究营养成分在不同厂家之间是否有显著性差异。营养成分数据都是定量数据，因此可以采用方差分析的思想来解决这个问题。为了数据表示的方便，我们将厂家A、B、C分别用数字1、2、3来表示。由于数据量比较大，解答过程用SPSS软件进行计算，而没有选用R软件。分析过程的显著性统一设定为0.05.解答过程：1. 先对数据做方差齐性检验，计算结果如下表所示：方差齐性检验

22、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.016

23、401.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% 置信区间