数据挖掘-使用R软件预测海藻数量

使用R软件预测海藻数量李强强2013.112/3/20231R软件R是一套完整的数据处理、计算和制图软件系统。其功能包括：数据存储和处理系统；数组运算工具（其向量、矩阵运算方面功能尤其强大）；完整连贯的统计分析工具；优秀的统计制图功能；简便而强大的编程语言：可操纵数据的输入和输出，可实现分支、循环，用户可自定义功能。R在语义上是函数设计语言。它允许在“语言上计算”。这使得它可以把表达式作为函数的输入参数，而这种做法对统计模拟和绘图非常有用。R是一个免费的自由软件。本案例使用的是R的3.0.1版。2/3/20232背景描述某些高浓度的有害藻类严重破坏着河流的生态环境，因此，能够监测并及早对海藻的繁殖进行预测对提高河流的质量是很有必要的。在约一年时间内，在不同的时间收集了多条不同河流的水样。每个水样测定了它们不同的化学性质和7种有害藻类的存在频率。还记录了如收集的季节、河流大小和水流速度。案例研究动机：1.化学监测相对人工检测价格便宜，且易于自动化。2.更好地了解藻类的频率和水样的某些化学性质以及其他特性(如季节、河流类型等)是如何相关的。2/3/202330海藻数据第一个数据集：包含200个水样。每条记录是一条河流在该年的同一季节的三个月内收集的水样的平均值。其中，每条记录由11个变量组成。其中3个名义变量：水样收集的季节、河流大小和河水速度。8个变量是水样的不同化学参数：最大pH值、最小含氧量、平均氯化物含量、平均硝酸盐含量、平均氨含量、平均正磷酸盐含量、平均磷酸盐含量和平均叶绿素含量。与之相关的是7种不同的有害藻类的频率数目。第二个数据集：140个不含7种藻类频率数目的测试集。2/3/202341数据加载1.点击文件菜单下的"改变工作目录"来设定当前工作目录。2.输入以下命令把文件中的数据读入：

algae<-read.table('Analysis.txt',

s=c('season',

'size',

'speed',

'mxPH','mnO2',

'Cl',

'NO3',

'NH4',

'oPO4','PO4','Chla',

'a1',

'a2',

'a3',

'a4',

'a5',

'a6',

'a7'),

na.strings=c('XXXXXXX'))3.点击文件菜单下的“保存工作空间”，输入文件名，退出，下次打开R后可通过拖拽的方式直接打开。2/3/202352数据摘要鉴于没有该问题领域足够的信息，首先了解一些数据的统计特性是一种较好的方式，它方便我们更好地理解问题。获取数据统计特性的方法之一是获取其描述性的统计摘要。命令如下：

summary(algae)对于名义变量，它给出了每个可能取值的频数。对于数值变量，它提供了均值、中位数、四分位数和极值等一系列统计信息。NA's表示缺失值的个数。通过观察这些值，我们可以了解到数据分布的偏度和分散情况。2/3/202363数据可视化(1)1.绘制变量mxPH的直方图的两种方式：

hist(algae$mxPH)

hist(algae$mxPH,prob=T) 区别在于前者给出的是频数，后者是区间的概率。2.绘制mxPH的Q-Q图： library(car) qqPlot(algae$mxPH,main='NormalQQplotofmaximumpH')

Q-Q图绘制变量值和正态分布的理论分位数的散点图。同时，它给出正态分布的95%置信区间的带状图(虚线)。main为设置图形的标题。2/3/202373数据可视化(2)3.绘制变量oPO4的箱图： boxplot(algae$oPO4,ylab="oPO4") rug(algae$oPO4,side=4) abline(h=mean(algae$oPO4,na.rm=T),lty=2) ylab为设置y轴标题；

rug函数绘制变量的实际值，side=4表示绘制在图的右侧(1在下方，2在左侧，3在上方)；

abline函数绘制水平线，mean表示均值，na.rm=T指计算时不考虑NA值，lty=2设置线型为虚线。2/3/202384数据清理 数据缺失的情形在实际问题中非常普遍。处理含有缺失值的数据时常用的几种策略：将含有缺失值的案例剔除。用中心趋势值来填补缺失值。根据变量之间的相关关系填补缺失值。根据案例之间的相似性填补缺失值。使用能够处理缺失值数据的工具(见下一节)。2/3/202394.1剔除缺失值(1)1.适用范围：含缺失值的记录在整个数据集中比例很小时。2.检查含缺失值的记录：

algae[!complete.cases(algae),

]3.剔除所有含缺失值的记录： algae<-na.omit(algae)4.找出每个记录中缺失值的个数： apply(algae,1,function(x)sum(is.na(x)))

函数apply()是元函数，可在某些条件下对对象应用其他函数。第二个参数“1”表示第一个参数algae中的对象的第一个维度，即行数据。第三个参数是临时函数，功能是计算对象x中NA的数量。(注：R中有TRUE=1,FALSE=0)2/3/2023104.1剔除缺失值(2)5.可根据4编写一个找出algae中含有给定数目缺失值的行。以下函数的功能是找出缺失值个数大于列数20%的行： library(DMwR) manyNAs(algae,0.2) 第二个参数如不指定，默认为0.2，下面的命令与上一条等价： manyNAs(algae)6.我们利用上面的函数来剔除缺失值较多的记录： algae<-algae[-manyNAs(algae),]

这里第二个参数的默认值为0.2。2/3/2023114.2用中心趋势值填补缺失值 代表中心趋势的值反映了变量分布的最常见值，这种方法也最自然、简便和快捷。对于接近正态的分布来说，均值是最佳选择；对偏态分布或有离群值的分布而言，中位数通常是更好的代表数据中心趋势的指标；对于名义变量，通常采用众数。可用以下函数完成填补所有缺失值： data(algae)

algae<-algae[-manyNAs(algae),

] algae<-centralImputation(algae) 上述方法特别适用于大数据集，但是这种方法可能导致较大的数据偏差，影响后期的数据分析。但使用复杂的无偏方法寻找最佳数据填补对大型数据集可能也不适用。2/3/2023124.3通过变量的相关性填补缺失值(1)

1.用以下命令获取变量之间的相关矩阵： data(algae) symnum(cor(algae[,4:18],use="complete.obs")) 其中，函数cor()产生相关值矩阵(忽略前3个名义变量)，use参数指计算相关值时忽略含有NA的记录。2.结果显示，有两个相关性较大的值：NH4和NO3之间，PO4和oPO4之间。 前者相关性不是特别明显(0.6~0.8)，考虑到只样本62和样本199含有过多的缺失值，若剔除它们，样本中NH4和NO3就不存在缺失值了。后者相关值很高(大于0.9)*，可用变量的相关性填补缺失值。*根据领域专家的解释，总的磷酸盐值包含正磷酸盐值。2/3/2023134.3通过变量的相关性填补缺失值(2)3.寻找PO4和oPO4之间的线性关系的方法： algae<-algae[-manyNAs(algae),] lm(formula=PO4~oPO4,data=algae) 我们得到线性模型：PO4=42.897+1.293*oPO4.4.剔除样本62和样本199后，仅样本28在PO4上有缺失值，我们用上面的线性关系来填补：

algae[28,"PO4"]<-42.897+1.293*algae[28,"oPO4"]

查看填补的记录： algae[28,]2/3/2023144.4通过案例的相关性填补缺失值1.度量相似性时，最常用的指标是欧式距离。我们可通过使用这种度量来寻找与任何含有缺失值的案例最相似的10个水样，并用它们填补缺失值。 方法一：简单计算这10个最近的案例的中位数并用中位数填补缺失值；若缺失值是名义变量则采用众数。 方法二：采用这些最相似数据的加权均值。这里用高斯核函数从距离获得权重。命令如下：

clean.algae<-knnImputation(algae,k=10)2.这种方法看起来更合理，但也可能存在不相关的变量扭曲相似性，甚至造成大型数据集的计算特别复杂等问题。因此，填补缺失值时，大多应根据分析领域的知识来确定。2/3/2023155获取预测模型用于预测海藻的两种模型：多元线性回归模型和回归树模型。线性回归不能使用有缺失值的数据集，而回归树模型可以很自然地处理带缺失值的数据。多元线性回归模型是最常用的统计数据分析方法，该方法给出了一个有关目标变量与一组解释变量关系的线性函数。由于多元线性回归模型中没有处理缺失值的方法，因此，我们可以做如下的数据预处理：data(algae)algae<-algae[-manyNAs(algae),]clean.algae<-knnImputation(algae,k=10)这里还是先移除缺失值较多的记录，然后根据训练集数据个案的相似性来填补缺失值。2/3/2023165.1线性回归模型(1) 建立用于预测海藻频率的线性回归模型：

lm.a1<-lm(a1~.,data=clean.algae[,1:12]) 函数lm()建立一个线性回归模型，其中，第一个参数给出了模型的函数形式。这个函数的形式是用数据中的其他所有变量来预测变量a1，第一个参数中的点“.”代表数据框中的所有除a1外的变量。如需要用变量mxPH和NH4来预测变量a1，就要定义模型为"a1~mxPH+NH4"。参数data是用来设定建模所用的数据集。2/3/2023175.1线性回归模型(2) 通过下面的代码，我们可以获取更多线性模型的信息：

summary(lm.a1) 首先，给出数据拟合的残差(residuals)。残差应该是均值为0且为正态分布的。 其次，对于每个多元线性回归方程的系数(变量)，给出它的估计值和标准误差，并使用t检验来验证系数为0的假设检验。 再者，给出模型与数据的吻合度，即模型所能解释的数据变差的比例。R-squared越接近于1说明模型拟合得越好，越小则代表模型拟合得越差。 最后，还可以检验任何解释变量与目标变量的依赖关系。

2/3/2023185.2回归树模型(1)因模型解释的方差比例太低，才约32%，故实际验证结果表明：对海藻案例应用线性模型是不合适的。用线性思维去考虑非线性问题，得不到理想的结果。我们考虑使用回归树预测。建立回归树： library(rpart) data(algae) algae<-algae[-manyNAs(algae),] rt.a1<-rpart(a1~.,data=algae[,1:12])函数的形式是用数据中其他所有变量来预测a1，data是用来设定建模所用的数据集。2/3/2023195.2回归树模型(2) 回归树rt.a1的图形表示的两种方法： plot(rt.a1) text(rt.a1) 或

prettyTree(rt.a1) 建立回归树通常分两步。最初，生成一棵较大的树，然后通过统计估计删除底部的一些结点来对树进行修剪。这样是为了防止过度拟合。2/3/2023206模型评价和选择 使用已有的训练数据获得模型的性能指标是不可靠的，因为这些计算是有偏的。实际上，有的模型可以很容易获得训练数据的零误差预测。然而，这一优秀性能很难推广到目标变量值未知的新样本上。这种现象我们通常称为过度拟合训练数据。我们需要一个模型，使它在未知数据上有可靠的预测性能。 k折交叉验证是获得模型性能可靠估计的一种常用方法，它适用于小数据集。2/3/202321k折交叉验证方法

K折交叉验证：初始采样分割成K个子样本，一个单独的子样本被保留作为验证模型的数据，其他K-1个样本用来训练。交叉验证重复K次，每个子样本验证一次，平均K次的结果或者使用其它结合方式，最终得到一个单一估测。常见的选择是k=10。有时还会重复进行多次K折交叉验证以获得更加可靠的估计。 总之，在做一项预测任务时要做出以下决策： 1)为预测任务选择模型； 2)选择比较模型性能的评估指标； 3)选择获取评估指标的可靠估计的实验方法。2/3/202322进行模型比较的函数 在R中，有函数experimentalComparison()可用来进行模型之间的选择和比较。它有三个参数：1)用于比较的数据集；2)需要比较的可选模型；3)实验过程中的系数。此函数适用于任何模型和任何数据，从这个意义上说，它是一个泛型函数。 使用者提供一组实现待比较的模型的函数，其中每一个函数应该对训练集和测试集实现一个完整的“训练+测试+评估”周期。在评估过程的每一次迭代中，调用这些函数。这些函数返回一个向量，其元素为交叉验证中用户需要的性能评估指标量。2/3/202323两个目标模型的函数cv.rpart<-function(form,train,test,...){ m<-rpartXse(form,train,...) p<-predict(m,test) mse<-mean((p-resp(form,test))^2) c(nmse=mse/mean((mean(resp(form,train))-resp(form,test))^2))}cv.lm<-function(form,train,test,...){ m<-lm(form,train,...) p<-predict(m,test) p<-ifelse(p<0,0,p) mse<-mean((p-resp(form,test))^2) c(nmse=mse/mean((mean(resp(form,train))-resp(form,test))^2))}2/3/202324函数的说明11.我们用NMSE(标准化后的平均绝对误差)作为回归树模型和线性回归模型的性能评估指标。2.这些函数的前三个参数是公式、训练数据和测试数据。3.特殊参数"..."可用在任意的R函数中，它允许一个特定函数具有可变的参数，它用来给实际模型传递所需要的额外参数。4.函数resp()用于根据公式获得数据集的目标变量值。2/3/202325模型的交叉验证比较

定义好用于模型学习和测试的函数后，我们可按下列代码进行模型的交叉验证比较：

res<-experimentalComparison(

c(dataset(a1~.,clean.algae[,1:12],'a1')),

c(variants('cv.lm'),

variants('cv.rpart',se=c(0,0.5,1))),

cvSettings(3,10,1234)

)

参数一的形式为dataset(<formula>,<dataframe>,<label>);

参数二包含可选的模型方法并通过variants()指定，该函数第一个参数为目标模型的函数名，第二个作为可选参数；

参数三设定交叉验证试验的参数(repetition,fold,seed)。2/3/202326查看比较结果1.查看比较结果的摘要：

summary(res) 从结果中可知，apart.v1有最优的平均NMSE值。2.验证结果也可转化成可视化图形：

plot(res)3.可通过以下代码查看模型所对应的参数： getVariant("cv.rpart.v1",res)2/3/202327多个预测任务同时进行DSs

<-

sapply(

names(clean.algae)[12:18], function(x,names.attrs){ f<-as.formula(paste(x,"~.")) dataset(f,clean.algae[,c(names.attrs,x)],x) }, names(clean.algae)[1:11])res.all<-experimentalComparison(

DSs, c(variants('cv.lm'),

variants('cv.rpart',se=c(0,0.5,1)) ), cvSettings(5,10,1234))2/3/202328函数的说明21.该代码首先创建用于比较7个预测任务的数据集向量。对每一个预测问题需构建一个公式，该公式由一个字符串构成，它是数据集中相应的需要预测的目标变量和符号"~."连接而成的。然后，该字符串通过函数as.formula()转换为一个R公式。2.这次采用重复5次10折交叉验证以提高统计结果的显著性。3.本条指令运行的时间稍长(处理器：Intel(R)Core(TM)i5-3210MCPU@2.50GHz，2501Mhz，2个内核，4个逻辑处理器，约需1分钟)。2/3/202329模型对不同海藻的结果1.所有海藻交叉验证结果的可视化：

plot(res.all) 图中显示有几个结果很差，即其NMSE值明显大于1。这意味着测试结果比简单采用目标变量的均值这一基准模型还要差！2.查看每个问题对应的最优模型的代码：

bestScores(res.all)

其结果说明，只有海藻1的预测结果尚可。3.考虑用组合法进行模型构建。通过产生大量可选模型并把其进行组合，这样得到的模型可以克服单个模型的局限性。2/3/202330随机森林1.随机森林作为组合模型的代表，它由大量的树模型(回归树或分类树)构成。每个树是完全生长的(没有事后剪枝)，在树生长的每一步，最好的结点分割方法将从变量集合的一个随机子集中选取。回归任务的预测采用组合中预测结果的平均值。2.以下代码是包含三个版本的随机森林模型的交叉验证，在组合中每个模型有不同数目的数。2/3/202331随机森林的交叉验证library(randomForest)cv.rf<-function(form,train,test,...){ m<-randomForest(form,train,...) p<-predict(m,test) mse<-mean((p-resp(form,test))^2) c(nmse=mse/mean((mean(resp(form,train))-resp(form,test))^2))}res.all<-experimentalComparison(

DSs, c(variants('cv.lm'),

variants('cv.rpart',se=c(0,0.5,1)),

variants('cv.rf',ntree=c(200,500,700))), cvSettings(5,10,1234))2/3/202332模型应用的结果1.采用下面的函数来证实组合方法的优势：

bestScores(res.all) 对于某些问题，随机森林给出很好的结果。但像海藻7，结果还不能令人满意。2.还不知道这些最佳模型和其它模型之间的区别是否显著。即，采用别的随机数据我们能得到相似结果的可能性有多大？ 我们用Wilcoxon(威尔科克森)检验来判断显著性。上述结果表明，对海藻1、2、4和6，模型"cv.rf.v3"最好，检验代码为：

compAnalysis(res.all,against='cv.rf.v3',

datasets=c('a1','a2','a4','a6'))2/3/202333显著性分析1.结果中的“sig.X”列是我们需要的信息。没有任何标识符则意味着相应的模型和"cv.rf.v3"模型之间有显著差异的可能性低于95%。加号意味着相应模型的平均性能估计指标显著高于模型"cv.rf.v3"。由于好的模型对应较低的NMSE值，所以该模型的性能比模型"cv.rf.v3"差。减号的含义相反。2.从结果可知，不同版本之间的随机森林模型的差异在统计上通常不显著。与其他模型相比，在大部分情况下，随机森林具有显著的优势。3.参数against和datasets取不同的值，可对在其他海藻上有最优性能的模型进行类似分析。2/3/2023347预测海藻频率 本案例的目的：预测140个水样的7个海藻的频率值。 我们已经采用交叉验证的过程给出了最佳的预测模型，下面应该应用所有可得的训练数据来构建模型，并将得到的模型应用到测试数据集。为了避免回归树采用它自身的缺失值的处理方法，我们采用了k近邻法填补数据框clean.algae的NA值。随机森林本身没有处理缺失值的方法，我们把数据框clean.algae作为它的训练集数据。1.为每种藻类选择最优的预测模型：

bestModelsNames<-sapply(bestScores(res.all),

function(x)x['nmse','system'])

learners<-c(rf='randomForest',rpart='rpartXse')

funcs<-learners[sapply(strsplit(bestModelsNames,'\\.'),

function(x)x[2])]2/3/202335parSetts<-lapply(bestModelsNames, function(x)getVariant(x,res.all)@pars)