预测海藻数量

1、R语言实验设计题目：海藻数量数据预处理与数据探索性分析姓名：樊海涛班级：12-1学号：1207060106专业：统计学预测海藻数量基本任务：数据预处理、探索性数据分析和预测模型的构建问题描述与目标某些高浓度的有害海藻对河流生态环境的强大破坏是一个严重问题，她们不仅仅破坏 河流的生物，也破坏水质。能够检测并在早期对海藻的繁殖进行预测对提高河流质量是很有 必要的。针对这一问题的预测目标，在大约一年的时间里，在不同时间内收集了欧洲多条不同河流的水样。对于每个水样，测定了它们的不同化学性质以及 7种有害藻类的存在频率，在水样的收集过程中，也记录了一些其他特征，如收集的季节、河流的大小和水流的速度。本案

2、例研究的动机之一是化学检测价格便宜，并且易于自动化。而通过分析生物样品来识别水中的藻类要涉及显微镜检验，需要训练有素的工作人员， 因此即昂贵又缓慢，因此，构建一个可以基于化学性质来准确预测藻类的模型将有助于建立监测有害藻类的廉价的自 动化系统。另一个目的是更好地了解影响藻类频率的因素。也就是说，我们要了解藻类的频率和 水样的某些化学性质以及其他性质（季节、河流类型等）是如何相关的。数据说明来自于ERUDIT研究网络，并被用于 1999年的COIL国际数据分析竞赛。在 R软件中可 以直接获得。本案例有两个数据集，第一个数据集有200个水样。名为：algae,更精确的说，该数据集的每一条记录是同一

3、河流在该年的同一个季节的三个月内收集的水样的平均值。每条记录由11个变量构成。其中3个变量是名义变量，它们分别描述水样收集的季节、收集样品的河流大小和河水速度。余下的8个变量是所观测水样的不同化学参数，即最大PH值最小含氧量O2平均氯化物含量CI平均硝酸盐含量NO3-平均氨含量NH4+平均正磷酸盐含量 PO43-平均磷酸盐含量PO4平均叶绿素含量与这些参数相关的是7种不同有害藻类在相应水样中的频率数目。并未提供所观察藻类的名称的相关信息。第二个数据集由140个额外观测值构成。名为：，它们的基本结构和第一个数据集一样，但是它不包含7种藻类的频率数目。可视为测试集。本案例的主要目标是预测140个水

4、样中7种海藻的频率。在这种问题中，任务是建立预测模型，并预测在给定预测变量的取值时相应的目标变量的值。说明哪一个预测变量对目标变量有较大的影响。加载数据利用R软件载入DMwR添加包，里面有我们需要的名为algae的数据框。这个数据框里含有前面提到的200个观测值：library(DMwR)head(algae) season size speed mxPH mnO2 Cl NO3 NH4 oPO4winter small mediumspring small mediumautumn small mediumspring small mediumautumn small mediumwinte

5、r small highPO4 Chla a1 a2 a3 a4 a5 a6 a7函数head()将显示数据框的前 6行。数据框的每一行代表一个观测值。数据可视化和摘要为后面的数据鉴于开始我们对该领域一无所知，首先我们要了解一些数据的统计特性,为后面的数据处理与建模提供更多的信息。获取数据统计特性的一个方法是获取数据的描述性统计摘要。 summary(algae)seasonsizespeedmxPHautumn:40large :45high :84Min. :spring:53medium:84low :331st Qu.:summer:45small :71medium:83Median

6、 :winter:62Mean :3rd Qu.:Max. :NAs :1mnO2ClNO3Min.Min. :Min. :1st Qu.:1st Qu.:1st Qu.:Median :Median :Median :Mean :MeanMean3rd Qu.:3rd Qu.:3rd Qu.:Max. :Max. :Max. :NAs :2NAs:10NAs :2NH4oPO4PO4Min.Min.Min.1st Qu.:1st Qu.:1st Qu.:Median :Median :Median :MeanMeanMean3rd Qu.:3rd Qu.:3rd Qu.:Max. :Max.

7、 :Max. :NAs :2NAs :2NAs :2Chlaa1a2Min. :Min.Min. :1st Qu.:1st Qu.:1st Qu.:Median :MedianMedian :Mean :MeanMean :3rd Qu.:3rd Qu.:3rd Qu.:Max. :Max. :Max. :NAs :12a3a4a5Min. :Min.Min. :1st Qu.:1st Qu.:1st Qu.:Median :Median :Median :Mean :MeanMean3rd Qu.:3rd Qu.:3rd Qu.:Max. :Max. :Max. :a6a7Min. :Min

8、.1st Qu.:1st Qu.:Median :Median :Mean :Mean3rd Qu.:3rd Qu.:Max. :Max. :这个简单函数立即给出了数据的的统计特征概括， 对于名义变量，他给出了每个取值的变量的频数。例如，从结果中可知冬季采集的水样比其他季节更多，河流size为最大的有45个，河流的流速low比较少。对于数值型变量，R为我 们提供了四分之一位数，中位数，均值，四分之三位数，极值等一系列信息。这 些统计信息提供了变量值分布的初步信息，在变量有缺失值的情况下，字符串NA后面的数值即为缺失值的个数，通过中位数，均值，四分位数的信息，我们可以了解数据分布的偏度和分散情况

9、。且这些信息大多数都可以通过图形来表达出来。我们看一个例子。hist(algae$mxPH,probability = T)上面程序绘制变量 mxPH的直方图。如下图，变量mxPH的分布非常接近正 态分布，它的值大多数都集中在变量的均值附近，我们做 Q-Q来检验该变量是 否是正态分布，加载包car,使用里面的()绘制Q-Q图Histogram ofalgae$mxPHooaaeSmxPHlibrary(car)par(mfrow=c(1,2)hist(algae$mxPH,prob=T,xlab = ,ylim = 0:1,main = Histogram of maximum pH valu

10、e) lines(density(algae$mxPH,=T) rug(jitter(algae$mxPH) (algae$mxPH,main=Normal QQ plot of maximum pH)istogram of maximum pH Jormal QQ plot of maxi工Q.x后售co工Q.x后售co311 3 norm quantiles上图的左边是添加核密度曲线的直方图，而在 X轴附近绘制变量的实际值， 从而容易辨别离群点，例如我们可以发现有两个字显著低于其它所有值。这种数 据检查是很重要的，因为它能确定数据样本中可能出现的错误。右边的 Q-Q图， 他绘制了变量值与正

11、态分布的理论分位数的散点图，同时他给出正态分布的95%的置信区间的带状图，从右上图知，变量有几个小的值明显在95%置信区间之外， 它们不服从正态分布。卜面给出另一个数据检验的例子，用它来检验变量OP04: boxplot(algae$oPO4,ylab=Orthophosphate)rug(jitter(algae$oPO4),side = 2)abline(h=mean(algae$oPO4, = T),lty=2)寸 OOM 寸 OOM 0 力方 dsord。工to箱图上方小横线上面的小圆圈表示与其他值比较特别大的值，通常认为是离群值，这意味着箱线图给出大量的信息， 它不仅给出变量的中心趋

12、势，也给出了 变量的发散情况和离群值。上图中与X轴平行的直线，是变量的均值所在的位置，将均值线与中位数线 进行比较，就可以知道变量的多个离群值使得作为变量中心的均值产生了扭曲。变量oPO4的分布集中在较小的观测值范围，因此分布为正偏。大多数水样的oPO4值都比较低，也有几个水样的观测值较高，甚至特别高。有时，当有离群值时，需要确定哪些有离群值的观测，可以使用图形法，例 如绘制NH4的值，将会注意到一个特别大的值，我们可以使用一下方法识别特 大值相应的水样：plot(algae$NH4,xlab =)abline(h=mean(algae$NH4, = T),lty=1)abline(h=mea

13、n(algae$NH4, = T)+sd(algae$NH4, = T),lty=2)abline(h=median(algae$NH4, = T),lty=3)identify(algae$NH4)上面函数第一条绘制变量的所有值，调用函数 abline绘制三条有用的直线， 第一条为均值，第二条为均值加一个标准差，第三条为中位数，对于离群值的识 别尽管这三条不是必须的，但它们能提供有用的信息。sszS00L 0II150100200sszS00L 0II150100200从上图中知道，水样153的NH4值是一个极大异常值，20, 34,35,88,89,153等 水样中NH4值都异常高，可以把

14、这些水样选择出来，看看他们与别的水样中海 藻有什么种类有什么差异。下面给出几种其他类型的数据检验的例子， 这些例子应用R的添加包lattice , 做Trellis图形。假设需要研究海藻变量al的值得分布，然而，这里需要分布如何依赖于其 他变量，就需要新的变量，新的工具。条件绘图是依赖于某个特定因子的图形表示，因子是一个为一个取值为有限 集合的名义变量。例如，对于 size的不同取值，可以绘制变量 al的一组箱图。 每个箱图是对应于变量size的某个特定值的水样子集。通过这些箱线图可以研究 名义变量size如何影响变量a1值得分布。bwplot(sizea1,data=algae,ylab =

15、 River Size,xlab = Algal A1)smallfTtediurnsmallfTtediurn0204060SOAlgal A1从图知，在规模较小的河流中，海藻 al的频率越高。这种箱图的另外一个类型是分为箱图，它可以给出绘制变量的更多信息，R添加包Hmisc可以绘制分位箱图，下面绘制上面例子中al变量的条件分位箱图： library(Hmisc) bwplot(sizea1,data=algae,panel = ,probs=seq(.01,.49,by=.01),datadensity=T,ylab=RiverSize,xlab = Algal A1)23S323S3Al

16、gal Ai上图中的点代表代表不同大小的河流中海藻频数的均值，而图中的竖线分别代表变量的第一份位数，中位数和第三分位数。图中的小竖线代表数据的真实取 值，这些值分布信息由分位数图来体现。分位数箱图提供的信息要多于传统的箱 图。例如我们可以确认上面的观测结论， 小型河流有更高频率的海藻，但我们也 观察到小型河流的海藻频率的分布比其他类型的河流的海藻频率的分布分散。这种类型的条件绘图不局限于名义变量，也不局限于单个因子。只需把连续 变量离散化，也同样可以进行条件绘图，下面给出两因子条件绘图的例子。 考虑 变量a3在给定变量season和变量mnO2下的条件绘图，变量mnO2是一个连续 变量，代码如

17、下: m1 stripplot(seasona3|m1,data = algae!(algae$mnO2),)得到的图形如下，上面的代码第一行是调用函数()对连续变量mnO2离散化，把该变量转化为因子类型。参数 number设置需要的区间的个数。每个区间的观测 值的个数相等。第二个调用函数stripplot(),该函数是lattice包中的一个绘图函数，他根据 另一个变量season把变量的实际值绘制到不同的图形中。然后对变量mnO2的每个不同的区间绘制不同的图形。这些区间按照从左到右，从上到下的顺序来排列。左下方对应的是较小的mnO2值。IlliI0102DM4Clll11mlmlIwint

18、er口E0皿summeriid 皿ocooo ospfingom 0 O 000 O 08autumn匚cmoomlT1 winterrm coo oo0mUDCCMD 0000sumrrwon o ooo asptinyMMOOO OD口皿 O 8 口 0oocmci ooo &oautumn口口 o oIlliIIIII101020304033数据缺失处理在许多水样中，一些变量含有缺失值。这种情形在现实问题中非常普遍，这会导致一些不能处理缺失值的分析方法无法应用。当我们处理含有缺失值的数据时，可以运用以下几种最常见的策略：将含有缺失值得案例剔除根据变量之间的相关关系填补缺失值根据案例之间的

19、相似性填补缺失值使用能够处理缺失值数据的工具最后一种方式是最严格的，因为它限定可我们可以使用的工具。然而，当对这些可以处理缺失值的数据挖掘工具有信心时，它们可能是一个好的选择。在下面的章节中我们将通过举例来演示如何在 R中实现上述处理缺失值的方法。如果你决定试一试本书中的代码，你需要了解它们不是互相补充的，由于每节都用不同的方法来处理这些数据，这意味着当你选择其他方式来处理这些缺失时，需要重新读取含有所有缺失值的原始数据。重新读取数据的代码如下：library(DMwR)data(algae).将缺失部分剔除剔除含有缺失数据的记录非常容易实现，尤其是当这些记录所占的比例在可用数据集中非常小的时

20、候，这个选择就比较合理。在剔除某些变量中至少含有一个缺失数据的所有观测值时，最好先检查观测值，或者至少得到这些观测值的个数，例如：algae!(algae),nrow(algae!(algae),)1616个含缺失的样本函数()产生一个布尔值向量，该向量的元素个数与algae数据框中的行数相同，如果数据框的相应行中不含 NA值(即为一个完整的观测值)，函数返回值就是 TRUE前面提到过 的“ ！”运算符，它是取逻辑否，因此上述指令显示了含有缺失值得水样记录。为了从数据框中剔除16个样本，我们可以简单的输入:algae - (algae)即使我们决定不使用剔除所有包含缺失值记录的极端方法，我们也

21、可以剔除某些观测值。因为这些样本的缺失值太多，所以它们几乎是无用的样本，如果采用复杂的方法来填补缺失值，就会导致较大偏差。需要注意的是，如果执行了前面的命令， 你需要重新读取数据， 因为这个指令已经剔除了所有缺失数据，所以接下来的命令就没有意义！观察这个样本中的数据，我们可以看到第 62条和第199条记录中的11个解释变量有6个是缺失值。这种情 况下，最好是剔除它们：algae - algae-c(62,199),在有些问题中，由于大量记录中含有缺失值，用上面的观察方法来检查数据的缺失值是 不可行的，所以需要找到缺失值较多的样本所在的行。下面的代码可以找到海藻数据集中每行数据的缺失值个数：ap

22、ply(algae,1,function(x) sum(x)函数apply ()属于R中功能非常强大的一类函数。这类函数又称为元函数，它们可以在某些条件下对对象应用其他函数。对函数apply ()而言，它把一个函数应用到数据框的每一行。这个被应用的函数在apply ()函数的第三个参数中给出，对数据框的每一行都分别调用该函数。在这个案例中我们使用一个临时函数。它只有在调用apply ()函数时才存在。另外，函数 apply ()的第三个参数也可以是一个正常的函数的函数名。临时函数的功能是计算对象X中NA的数量。在R中逻辑值TRUE等于数值1,逻辑值FALSE等于数值0, 这意味着当加一个布尔值

23、向量时，得到向量中取值为TRUE的元素的个数。根据以上代码，可以编写一个程序找到algae中含有给定数目缺失值的行。可以如下应用该函数：data(algae)manyNAs(algae,1 62 199只有在前面操作中剔除了缺失值较多的几行数据后才需要调用函数data ()。函数manyNA是()的功能是找到缺失值个数大于列数20%的行。在第二个参数中可以设置一个精确的列数作为界限。因此，用下面的代码就无须知道含有缺失值较多的行的具体数量。algae - algae-manyNAs(algae),以上操作可以把样本中 62, 199号水样数据剔除。用高频率值来填补缺失值填补含有缺失值记录的另一

24、种方法是尝试找到这些缺失值最可能的值，同样，这里有多种策略可供选择，不同策略对逼近程度和算法复杂度的权衡不同，填补缺失数据最简单和快捷的方法是使用一些代表中心趋势的值，代表中心趋势的值反映了变量分布的最常见的值，因此中心趋势值是最自然的选择。有多个代表数据中心趋势的指标，例如平均值，中位数，众数等。最合适的选择由变量的分布决定。对于接近正态的分 布来说，所有的观测值都较好地聚集在平均值周围，平均值数就是最佳选择。然而，对于偏 态分布，或者离群值的变量来说， 选择平均值就不好。偏态分布的大部分值都聚集在变量分布的一侧，因此平均值不能作为最常见的代表。另一方面，离群值(极值)的存在会扭曲平 均值，

25、这就导致了平均值不具有代表性的问题。因此，在对变量分布进行检查之前选择平均值作为中心趋势的代表是不明智的，例如某些R的绘制工具。对偏态分布或者有离群值的分布而言，中位数是更好的代表数据中心趋势的指标。比如，样本algae48,中变量mxPH有缺失值。由于该变量分布近似正态分布我们可以选用平均值来填补这个“洞”，计算方法如下：algae48,mxPH - mean(algae$mxPH,=T)这里，函数mean ()计算数值向量的平均值，参数 =T使计算时忽略缺失数据。大多数时候采用一次填补一列中的所有缺失值而不是像上面那样一行一行地逐个填补。以变量Chla为例，这个变量在第12行上有缺失值。另

26、外，这也是平均值不能代表大多数变 量值的一种情况。事实上，Chla的分布偏向于较低的数值，并且它有几个极端值，这些都使 得平均值不能代表大多数的变量值。因此，我们使用中位数来填补这一类的缺失值：algae(algae$Chla),Chla options(digits = 1) cor(algae,4:18,use=)mxPH mnO2 Cl NO3 NH4 oPO4 PO4 Chla a1 mxPH mnO2 Cl NO3 NH4 oPO4 PO4 Chla a1 a2 a3 a4 a5 a6 a7a2 a3 a4 a5 a6 a7 mxPH mnO2 Cl NO3 NH4 oPO4 PO4

27、 Chla al a2 a3 a4 a5 a6 a7函数cor 0的功能是产生变量之间的相关值矩阵（因为前3个变量时名义变量，所以计算相关值时不考虑它们）。设定参数use= 时，R在计算相关值日忽略含有NA的记录。相关值在1 （或-1）周围表示相应的两个变量之间有强正（或负）线性相关关系。然后其他R函数可以得到变量间线性相关的近似函数形式，它可以让我们通过一个变量的值计算出另一个变量的值。函数cor 0的输出结果并不是很清晰，可用 symnum来改善结果的输出形式. symnum（cor（algae,4:18,use=）mP mO Cl NO NH o P Ch a1 a2 a3 a4 a5

28、a6 a7 mxPH 1 mnO2 1 Cl1 TOC o 1-5 h z NO31NH4, 1oPO4. .1PO4. .* 1Chlaal TOC o 1-5 h z a2 .1a31a4 . .1a51a6. . 1a71 attr(,legend)1 0 ，一 +，* 一 b1这种用符号表示相关值得方法更为清晰，特别是对于大的相关矩阵。在本案例的数据中，大多数变量之间是不相关的，然而，有两个例外：变量NH4和NO3之间，变量PO4和oPO4之间。后两个变量之间的相关值很高 (大于)。变量NH4和NO3之间 的相关性不是特别明显(为)，因此，根据它们来确定缺失数据是危险的。此外，因为样本

29、62和样本199有太多的变量含有缺失值，所以如果剔除它们，样本中的变量 NH4和NO3就没有缺失值了。至于变量PO4和oPO4,它们之间相关性可以帮助填补这两个变量的缺失值。为了达到这个目标，我们需要找到这两个变量之间的线性关系，方法如下： data(algae) algae lm(PO4oPO4,data = algae)Call:lm(formula = PO4 oPO4, data = algae)Coefficients:(Intercept)oPO4述该函数。线性模型是|P04 = 1 290Po4 + 42 9。如果这两个变量不是同时有缺失值，那 么可以用这个公式计算这些变量的缺失

30、值。在剔除样本62和样本199后，还剩下一个样本28在PO4上有缺失值，可以简单地使用上面的线性关系计算缺失值的填补值： algae28,PO4 algae28,PO41 48通过变量之间的相关关系求出水样28, PO4的缺失值填补为 48。对线性关系的研究使我们能够填充一些新的缺失值。然而还有几个观测值含有缺失值。可以试着探索案例数据中含有缺失值的变量和名义变量之间的关系，这可以通过R添加包lattice中的函数来绘制条件直方图来进行。绘图代码如下：histogram(mxPH|season,data=algae)summerwintersummerwinterSDnriilSDnriil6

31、7 B 910mxPH上图是绘制在不同的季节变量mxPH的直方图。mxPH直方图在不同季节没有显著的不同，下面对河流的大小size进行上面类似的分析，执行命令：histogram(mxPH|season,data=algae)B 73910laraemediumsnialle 7&9 ke t s & idmsPH那么从得到的直方图知，较小的河流有较小的mxPH值。对这种相关性研究，可以扩展到多个名义变量： stripplot(sizemxPH|speed,data=algae,jitter=T)上面的指令结果如下图，参数jitter=T说明丫的变量值进行小范围的随机排列，这些可以避免相同值之

32、间的互相重叠而导致导致失去具有某些特定值的观察值集中度的信息。这种类型的分析可以应用到其它含有缺失值的变量中，然而这种分析是一种繁琐的过 程，有大量的变量组合需要分析。适合名义变量较少的数据集的分析。通过探索案例之间的相似性来填补缺失值不同于探索数据集列之间的相关性，尝试使用行之间的相似性来填充缺失值。我们可以使用这种方法来填补剔除了那两个含有太多的NA值得样本外的其它缺失数据。如上所描述的方法假如两个水样是相似的，其中某些变量有缺失值，那么该缺失值很可能与另一个水样的值是相似的。为了使用这种直观的方法， 首先定义相似性的概念。 相似性经常由描述观察值的多元度量空间的变量所定义。在文献中有许多

33、度量相似性的指标，常用的是欧式距离。这个距离可以非正式的定义为任何两个案例的观察值之差的平方和，计算公式如下：下面要描述的方法是使用这种度量来寻找与任何含有缺失值的案例最相似的10个水样，并用它们来填充缺失值。我们考虑两种应用这些值的方法。第一种方法简单地计算这10个最相近的案例的中位数并用这个中位数来填充这些缺失值。第二种方法是采用这些相似数据的加权均值。权重的大小随着距待填补缺失值的个案的距离为d,则它的值在加权平均中的权重为：上面的方法可以通过添加包函数knnImputation()来实现，这个函数用一个欧式距离的变种来找到距离任何个案最近的k个邻居。在计算距离时都要对数据进行标准化。下

34、面介绍如何使用 knnImputation()函数 algae-algae-manyNAs(algae),algae-knnImputation(algae,k=10)algae如果使用中位数来填补缺失值，可以使用下面的代码algae algae 在运行以上的代码后，得到的数据框将不含缺失值。接下来，将建立一个用于预测海藻频率的线性回归模型: summaryCall:lm(formula = a1 ., data = , 1:12)Residuals:Min 1Q Median 3Q MaxCoefficients:Estimate Std. Error t value Pr(|t|) (In

35、tercept).seasonspring seasonsummer seasonwinter sizemedium sizesmall*speedlow speedmedium mxPH mnO2 Cl NO3*NH4oPO4PO4Chla Signif. codes: 0 *. 0 1Residual standard error: 18 on 182 degrees of freedomMultiple R-squared: ,Adjusted R-squared:F-statistic: on 15 and 182 DF, p-value:在解释函数summary。应用到线性模型对象所

36、得到的信息之前，先介绍 R处 理3个名义变量。当像上面一样进行模型构建时 R会生成一组铺助变量，即对每 个有k个水平的因子变量，R会生成k-1个铺助变量。这些铺助变量的值为 0或 者为1,当铺助变量为1表明该因子值出现，同时表明其他铺助变量的值为 0. 如果所有这k-1个铺助变量的取值都为0,则表明因子变量的取值为第 k个剩余 的值。对得到的线性模型对象应用函数 summary。，将给出所建模型的一些诊断信 息。首先是有关线性模型中数据拟合残差的。 残差应该是均值为0并且为这个态 分布。可以使用plot的命令来得到一系列的线性模型图，它们有助于了解模型的性 能，其中的一幅图形绘制拟合的目标变量

37、值和模型残差的散点图。R给出的另一图形是误差正态的Q-Q图，通过它可以检验是否符合应有的正态分布。Q 守 _Resi duals vs l-ittediriQ 守 _Resi duals vs l-ittedirifiI-10010203040fitted valuesIm（a1、.）norm quantites该模型解释的方差比例还是不太理想，大约为 32%,还可以拒绝目标变量不依赖预测变量的假设。检验某些系数的显著性，可能要考虑某些变量是否应该进 入模型中，下面将使用向后消元法。首先用函数anova()来精简线性,g型。当将 anova()应用到简单线性模型时， 这个函数提供一个模型拟合的

38、方差序贯分析， 也就是说，随着公式的项数的增加, 模型的残差平方和减少。对前面建立的模型进行方差分析。 anovaAnalysis of Variance TableResponse:a1Df Sum Sq Mean Sq F valuePr(F)season385size2 11401*speed23934*mxPH11329*mnO212287*Cl14304*NO313418*NH41404oPO414788*PO411406*Chla1377Residuals 182 56668 Signif. codes: 0 * *，.，1上面的结果表明season对减小模型拟合误差的贡献最小，下

39、面将它从模型中 剔除： summaryCall:lm(formulaa1size lm(formulaa1size + speed + mxPH + mnO2 + Cl + NO3 + NH4 +oPO4 + PO4 + Chla, data = , 1:12)Residuals: Min 1Q Median 3Q MaxCoefficients:Estimate Std. Error t value Pr(|t|)(Intercept).sizemedium sizesmall*speedlow speedmedium mxPH mnO2 Cl NO3*NH4 oPO4 PO4.Chla S

40、ignif. codes: 0 *. 0 1Residual standard error: on 185 degrees of freedomMultiple R-squared: ,Adjusted R-squared:F-statistic: on 12 and 185 DF, p-value:新模型的拟合指标提高到了 ,仍然不是太理想。下面使用anova()函数对两个模型进行比较正式的比较，但这次使用两个模型作为参数： anova,Analysis of Variance TableModel 1: a1 season + size + speed + mxPH + mnO2 + Cl

41、 + NO3 + NH4 + oPO4 + PO4 + ChlaModel 2: a1 size + speed + mxPH + mnO2 + Cl + NO3 + NH4 + oPO4 + PO4 +ChlaRSS Df Sum of Sq F Pr(F)1182 566682185 57116 -3上面的函数通过F检验对两个模型进行方差分析，据此评价两个模型是否有 显著差异。这种情况下，尽管误差平法和减少了-448,但比较的结果表示两者的差距并不显著。为了检验是否能移除更多的系数。为了简化向后消元过程， R有 一个函数来执行上面的所有过程。下面的代码对初始模型用向后消元法得到一个新的线性

42、模型。56668mnO2157363NH4157494PO4157567size 258526NO3 159008Step: AIC=a1 size + speed + mxPH + mnO2 + Cl + NO3 + NH4 + oPO4 + PO4 + ChlaDf Sum of Sq RSS AICspeed 257327-Chla157462-Cl157520-mnO2157572-mxPNH4157892-PO4157977-size259292-NO3159537-oPO4 157124Step: AIC=alsize +mxPH + mnO2 + Cla

43、lsize +mxPH + mnO2 + Cl+ NO3 + NH4 + oPO4 + PO4 +ChlaDf Sum of Sq RSS AICoPO4 157343Df Sum of Sq RSS AICoPO4 157343Chla 157550mnO2 157740Cl157799mxPH157810none57327NH4158047PO4158136size 259388NO3 159706Step: AIC=a1 size + mxPH + mnO2 + Cl + NO3 + NH4 + PO4 + ChlaDf Sum of Sq RSS AICChla 157551mnO2

44、157746Cl157814mxPHH4158047size 259393-PO4 163161Step: AIC=a1 size + mxPH + mnO2 + Cl+ NO3 + NH4 +Df Sum of Sq RSSmnO2 157986CH4158298mxPH 158384size 259768NO3160218PO4163860Step: AIC=AICa1 size + mxPH + Cl + NO3 +NH4 + PO4Df Sum of Sq RSSNH4 158517CxPH 158805si

45、ze 260464NO3160237PO4167084AICStep: AIC=a1 size + mxPH + Cl + NO3 +PO4-NO3 159713PO4Df Sum of Sq RSS AIC58517-mxPH159301-Cl159353-NO3160505-size261181-PO4167093可以通过下面的代码来获得最后模型的信息: summaryCall:lm(formula = al size + mxPH + Cl + NO3 + PO4, data =,1:12)Residuals:Min 1Q Median 3Q MaxCoefficients:Estima

46、te Std. Error t value Pr(|t|)(Intercept) sizemedium sizesmall mxPH ClNO3PO4*(Intercept) sizemedium sizesmall mxPH ClNO3PO4*Signif. codes: 0 *. 0 1Residual standard error:Multiple R-squared:,Residual standard error:Multiple R-squared:,Adjusted R-squared:这样的R这样的R2表明对海藻案例应用这个模型的解释的方差比例任然不是很乐观, 线性模型是不合适的回归树给出R的另一种回归模型，通过建立回归树来预测海藻a1出现的频率。由于这类模型可以处理缺