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

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

《薛毅数学模型数学建模第八次作业多元分析实验.docx》由会员分享，可在线阅读，更多相关《薛毅数学模型数学建模第八次作业多元分析实验.docx（34页珍藏版）》请在冰点文库上搜索。

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

数学模型第八次作业多元分析实验

8.1实验目的与要求

●学会对数据进行线性回归分析、预测与回归诊断

●学会对数据进行方差分析和判别分析

●建立相应的统计模型，用R软件计算，并对计算结果进行分析和讨论

8.2基本实验

1.回归分析

为估计山上积雪融化后对下游灌溉的影响，在山上建立一个观测站，测量最大积雪深度X（米）与当年灌溉面积Y（公顷），测得连续10年的数据如表8.1所示。

（1）建立一元线性回归模型，求解，并验证系数、方程或相关系数是否通过检验；

（2）现测得今年的数据是X=7米，给出今年灌溉面积的预测值、预测区间和置信区间（α=0.05）；（3）将数据散点、回归预测值、回归预测区间和置信区间均匀化在一张图上，分析线性回归的拟合情况。

解：

（1）为了研究这些数据中所蕴含的规律性，根据10对数据利用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题得到的回归模型作回归诊断：

（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<-c（1907,1287,2700,2373,3260,3000,1947,2273,3113,2493）

+）

>lm.sol<-lm（y~1+x,data=intellect）

>summary（lm.sol）

结果同1题，调用influence.measures（）函数进行诊断：

>influence.measures（lm.sol）

得到回归诊断结果：

Influencemeasuresof

lm（formula=y~1+x,data=intellect）:

dfb.1_dfb.xdffitcov.rcook.dhatinf

1-0.34940.2706-0.44791.1650.099410.157

2-1.67001.5077-1.73200.8591.065030.413*

30.0232-0.0475-0.10531.4610.006270.126

4-0.02710.0056-0.08891.4230.004470.100

50.5656-0.6935-0.82081.4440.326490.349

6-0.04730.06630.09641.5940.005280.190

71.2525-1.05771.40850.4440.580590.229*

80.2329-0.15310.37861.1200.071170.120

9-0.18840.25360.34651.4740.064570.215

100.01330.00460.07301.4320.003020.100

做出回归诊断图并保存：

op<-par（mfrow=c（2,2）,mar=0.4+c（4,4,1,1）,oma=c（0,0,2,0））

plot（lm.sol,1:

4）

par（op）

savePlot（"diagnoses1",type="eps"）

得到回归诊断图：

由回归结果得知第2点和第7点为强影响点。

分析回归诊断图：

（1）第一张为残差图，由于残差的分布与回归值的分布相吻合，因此可以认为残差的方差满足齐性。

考虑点7所对应的样本为异常值。

第二张为正态QQ图，所有点基本在一条直线上，7号点较远一些，也就是说可以认为残差满足正态性。

（2）第三张图为标准差的平方根与预测值的散点图。

第7点的值较大在1.2以外，有可能是异常值点。

（3）第四章图给出了Cook距离值，从图上来看，2号点的Cook距离最大，这说明2号点可能是强影响点（高杠杆点）。

通过上述判断7号点是异常值点，将它在后面的计算中剔除。

2号点是强影响点，加权计算中减少它的影响。

>weights<-rep（1,n）;weights[2]<-c（0.5）

>lm.correct<-lm（y~1+x,data=intellect,subset=-7,

+weights=weights）

在程序中，subset=-7是去掉7号点。

weights<-rep（1,n）是将所有点的权赋为1，weights[2]<-0.5是再将18号点的权定义为0.5。

这样可以直观认为，2号点对回归方程的影响减少一半，再次进行回归分析：

>summary（lm.correct）

得到如下结果：

Call:

lm（formula=y~1+x,data=intellect,subset=-7,weights=weights）

WeightedResiduals:

Min1QMedian3QMax

-93.962-60.875-9.21336.037115.036

Coefficients:

EstimateStd.ErrortvaluePr（>|t|）

（Intercept）64.97118.670.5470.601

x373.7517.4021.4851.19e-07***

---

Signif.codes:

0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

Residualstandarderror:

71.08on7degreesoffreedom

MultipleR-squared:

0.9851,AdjustedR-squared:

0.9829

F-statistic:

461.6on1and7DF,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（x[c（7,2）],y[c（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","CorrectReg"）,

pch=c（19,NA,NA）,lty=c（NA,1,5）,

col=c（"orange","blue","red"））

savePlot（"diagnoses2",type="eps"）

图中实线是原始数据计算的结果，虚线是修正数据后的计算结果。

再次进行回归诊断：

>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）

>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<-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且比之前的值都低。

由此得到了最优回归方程：

4.方差分析І（单因素方差分析）

24只小鼠随机地分成3组，每组8只，每组喂养不同的饲料组，喂养一定时间后，测得小鼠肝中铁含量，结果如表8.3所示。

（1）试用单因素方差分析方法分析不同饲料的小鼠肝中铁含量是否有显著差异？

（2）如果有显著差异，哪些水平之间有显著差异？

（3）数据是否满足正态性和方差齐次性的要求，如果不满足请选择合适的分析方法，并给出你的最终结论。

解：

（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情况下小鼠肝中铁含量无显著