数学建模作业 实验7多元分析实验.docx
《数学建模作业 实验7多元分析实验.docx》由会员分享，可在线阅读，更多相关《数学建模作业 实验7多元分析实验.docx（24页珍藏版）》请在冰点文库上搜索。
数学建模作业实验7多元分析实验
实验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="X",ylab="Y",cex=1.4,pch=19,col="red"）
得到如下图像：
分析图像，数据点大致落在一条直线附近，说明变量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:
Min1QMedian3QMax
-128.591-70.978-3.72749.263167.228
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）140.95125.111.1270.293
x364.1819.2618.9086.33e-08***
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
96.42on8degreesoffreedom
MultipleR-squared:
0.9781,AdjustedR-squared:
0.9754
F-statistic:
357.5on1and8DF,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="prediction",
+level=0.95）
得到如下结果：
fitlwrupr
12690.2272454.9712925.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="X",ylab="Y",cex=1.4,pch=19,col="red"）
>
>lm.sol<-lm（y~1+x）
>summary（lm.sol）
做出图像：
>abline（lm.sol,lwd=2,col="blue"）
>segments（x,fitted（lm.sol）,x,y,lwd=2,col="blue"）
标注图像：
>ex1<-expression（paste（"（",x[i],",",y[i],"）"））
>ex2<-expression（paste（"（",x[i],",",hat（y）[i],"）"））
>
>points（x[8],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）首先根据题意建立多元线性回归方程：
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:
Min1QMedian3QMax
-28.349-11.383-2.65912.09548.807
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）43.6500718.054422.4180.02984*
X11.785340.539773.3080.00518**
X2-0.083290.42037-0.1980.84579
X30.161020.111581.4430.17098
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
19.97on14degreesoffreedom
MultipleR-squared:
0.5493,AdjustedR-squared:
0.4527
F-statistic:
5.688on3and14DF,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
DfSumofSqRSSAIC
-X2115.75599.4109.32
5583.7111.27
-X31830.66414.4111.77
-X114363.49947.2119.66
Step:
AIC=109.32
Y~X1+X3
DfSumofSqRSSAIC
5599.4109.32
-X31833.26432.6109.82
-X115169.510768.9119.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:
Min1QMedian3QMax
-29.713-11.324-2.95311.28648.679
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）41.479413.88342.9880.00920**
X11.73740.46693.7210.00205**
X30.15480.10361.4940.15592
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
19.32on15degreesoffreedom
MultipleR-squared:
0.5481,AdjustedR-squared:
0.4878
F-statistic:
9.095on2and15DF,p-value:
0.002589
此时回归系数检验的水平已有显著提升，但X3项系数仍然不显著。
利用drop1（）函数计算：
>drop1（lm.step）
得到如下结果：
Singletermdeletions
Model:
Y~X1+X3
DfSumofSqRSSAIC
5599.4109.32
X115169.510768.9119.09
X31833.26432.6109.82
此时的结果说明，去掉X3项的时候，AIC值和残差平方值上升都是最小的，因此去掉X3项再次做线性回归：
>lm.opt<-lm（Y~X1,data=import）;
>summary（lm.opt）
得到结果如下：
Call:
lm（formula=Y~X1,data=import）
Residuals:
Min1QMedian3QMax
-31.486-8.282-1.6745.62359.337
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）59.25907.42007.9865.67e-07***
X11.84340.47893.8490.00142**
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
20.05on16degreesoffreedom
MultipleR-squared:
0.4808,AdjustedR-squared:
0.4484
F-statistic:
14.82on1and16DF,p-value:
0.001417
此时常数项的检测结果为极为显著，X1项系数为很显著。
方程式P-值为0.001417<0.05且比之前的值都低。
由此得到了最优回归方程：
3.方差分析I
解：
（1）首先提出假设H0不同饲料的小鼠肝中铁含量无显著差异，μ1=μ2=μ3；H1不同饲料的小鼠肝中铁含量有显著差异，μ1，μ2，μ3不全相等。
使用R软件求解，用数据框的格式输入数据，调用aov（）函数计算方差分析，编程如下：
>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.lm<-lm（X~A,data=mouse）
>anova（mouse.lm）
得到如下结果：
AnalysisofVarianceTable
Response:
X
DfSumSqMeanSqFvaluePr（>F）
A273.11836.5599.1040.001422**
Residuals2184.3294.016
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.从结果中看到P-值为0.001422<0.05，因此原假设是不成立的，拒绝H0，即不同饲料的小鼠肝中铁含量有显著差异。
（2）继续使用R软件来分析哪些水平之间有显著差异。
首先计算数据在各水平下的均值：
>attach（mouse）
>tapply（X,A,mean）
得到如下结果：
123
1.606253.107505.82375
可以看出不同饲料喂食下的小鼠肝中铁含量的均值已有明显差异。
再做多重t检测：
>pairwise.t.test（X,A）
得到如下结果：
PairwisecomparisonsusingttestswithpooledSD
data:
XandA
12
20.1489-
30.00120.0262
Pvalueadjustmentmethod:
holm
由计算结果得出结论，μ1与μ3、μ2与μ3是有显著差异的，而μ1与μ2没有显著差异。
即是说，喂食饲料A和喂食饲料B情况下小鼠肝中铁含量有显著差异；喂食饲料B和喂食饲料C情况下小鼠肝中铁含量有显著差异；喂食饲料A和喂食饲料B情况下小鼠肝中铁含量无显著差异。
进一步，使用plot（）函数画出线箱图并保存：
>plot（X~A,col=5:
7,
+main="Box-and-WhiskerPlotofMouseData"）
>detach（mouse）
>
>savePlot（"box_plot",type="eps"）
可以直观看到数据的水平及各因素之间的差异。
（3）根据题意，先编写程序，做Shapiro-Wilk正态性检验
>attach（mouse）
>tapply（X,A,shapiro.test）
得到如下结果：
$`1`
Shapiro-Wilknormalitytest
data:
X[[1L]]
W=0.8742,p-value=0.1656
$`2`
Shapiro-Wilknormalitytest
data:
X[[2L]]
W=0.8893,p-value=0.2306
$`3`
Shapiro-Wilknormalitytest
data:
X[[3L]]
W=0.985,p-value=0.9833
结果显示三组数据均数据满足正态性。
再用Bartlett函数做方差齐性检验：
>attach（mouse）
>bartlett.test（X,A）
得到如下结果：
Bartletttestofhomogeneityofvariances
data:
XandA
Bartlett'sK-squared=10.5677,df=2,p-value=0.005073
从结果中看到p-值为0.005073<0.05因此认为数据并不满足方差齐性。
对于只满足正态性，不满足齐性要求的数据，用函数oneway.test（）作方差分析：
>oneway.test（X~A,data=mouse）
得到方差的分析结果：
One-wayanalysisofmeans（notassumingequalvariances）
data:
XandA
F=10.3592,numdf=2.00,denomdf=10.51,p-value=0.003271
此时P-值较第一问计算时的结果有所增大，但是仍然满足p-值<0.05因此可以认为原假设是不成立的，拒绝H0，即不同饲料的小鼠肝中铁含量有显著差异。
4.方差分析II
解：
（1）设有A、B两个因素，因素A有3个水平A1、A2、A3，因素B有3个水平B1、B2、B3。
利用R软件来进行判断：
>tree.aov<-aov（Y~A+B+A:
B,data=tree）
>factory<-data.frame（
+Y=c（4.6,4.3,6.1,6.5,6.8,6.4,6.3,3.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）,
+A=gl（3,2,18）
+）
>factory.aov<-aov（Y~A+B+A:
B,data=factory）
>summary（factory.aov）
得到结果：
DfSumSqMeanSqFvaluePr（>F）
A25.4082.7046.1300.02090*
B27.8413.9218.8880.00740**
A:
B412.1923.0486.9100.00793**
Residuals93.9700.441
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1
结果显示有交互作用时在显著水平0.05下，因素A（反应温度）效应显著，而因素B（反应压力）和交互效应很显著。
（2）使用R软件求解各种反应温度下产量均值的估计：
>attach（factory）
>tapply（Y,A,mean）
得到如下结果：
123
4.6500004.5333335.750000
计算各种反应压力下产量均值的估计：
>tapply（Y,B,mean）
得到如下结果：
123
5.7833334.1666674.983333
计算同时考虑温度和压力下产量均值的估计：
>matrix（tapply（Y,A:
B,mean）,nr=3,nc=3,byrow=T,
+dimnames=list（levels（A）,levels（B）））
得到如下结果：
123
14.455.04.50
26.303.63.70
36.603.96.75
（3）从第一问结果因素A（反应温度）效应显著、而因素B（反应压力）和交互效应很显著来看第二问得到的数据即可得到答案。
各种反应温度下产量均值中3条件下最多；
各种反应压力下产量均值中1条件下最多；
交互效应下3、3条件的产量均值最多，且高于单独作用时的产量均值；
综合看来，选用3、3条件是最佳的，即采用80℃的反应温度
3公斤的反应压力时对生产最有利。
7.3加分实验
解：
根据题意，明确解题思路，要解决的问题是：
是否有理由认为某一厂家的产品比其他厂家的产品更“有营养”（高蛋白、低脂肪、高纤维、低糖等）？
也就是研究营养成分在不同厂家之间是否有显著性差异。
营养成分数据都是定量数据，因此可以采用方差分析的思想来解决这个问题。
为了数据表示的方便，我们将厂家A、B、C分别用数字1、2、3来表示。
由于数据量比较大，解答过程用SPSS软件进行计算，而没有选用R软件。
分析过程的显著性α统一设定为0.05.
解答过程：
1.先对数据做方差齐性检验，计算结果如下表所示：
方差齐性检验
Levene统计量
df1
df2
显著性
热量
6.665
2
40
.003
蛋白质
1.676
2
40
.200
脂肪
6.045
2
40
.005
钠
7.146
2
40
.002
纤维
2.428
2
40
.101
碳水化合物
.917
2
40
.408
糖
.729
2
40
.489
钾
3.266
2
40
.049
由上表可以看出，在0.05的显著性水平下，热量、脂肪、钠、钾三个变量没有通过方差齐性检验，其它都是方差齐性的。
因此对热量、脂肪、钠、钾三个变量做方差非齐性的方差分析，其余变量做方差齐性的方差分析模型。
2.方差分析
（1）方差齐性变量的方差分析结果：
ANOVA
平方和
df
均方
F
显著性
蛋白质
组间
.682
2
.341
.220
.804
组内
62.016
40
1.550
总数
62.698
42
纤维
组间
10.884
2
5.442
1.740
.189
组内
125.088
40
3.127
总数
135.972
42
碳水化合物
组间
130.318
2
65.159
4.131
.023
组内
630.868
40
15.772
总数
761.186
42
糖
组间
47.564
2
23.782
1.165
.322
组内
816.715
40
20.418
总数
864.279
42
从结果可以看出，在0.05的显著性水平下，三个厂商在碳水化合物上有显著性差异，其余变量没有显著性差异。
下面进一步进行两两比较分析，看不同厂商的差异程度，如下表所示：
多重比较
因变量
（I）厂商
（J）厂商
均值差（I-J）
标准误
显著性
95%置信区间
升级会员 