使用R软件预测海藻数量李强强2013.113/10/20201:37AM2R软件•R是一套完整的数据处理、计算和制图软件系统。其功能包括:数据存储和处理系统;数组运算工具(其向量、矩阵运算方面功能尤其强大);完整连贯的统计分析工具;优秀的统计制图功能;简便而强大的编程语言:可操纵数据的输入和输出,可实现分支、循环,用户可自定义功能。•R在语义上是函数设计语言。它允许在“语言上计算”。这使得它可以把表达式作为函数的输入参数,而这种做法对统计模拟和绘图非常有用。•R是一个免费的自由软件。本案例使用的是R的3.0.1版。3/10/20201:37AM3背景描述•某些高浓度的有害藻类严重破坏着河流的生态环境,因此,能够监测并及早对海藻的繁殖进行预测对提高河流的质量是很有必要的。•在约一年时间内,在不同的时间收集了多条不同河流的水样。每个水样测定了它们不同的化学性质和7种有害藻类的存在频率。还记录了如收集的季节、河流大小和水流速度。•案例研究动机:1.化学监测相对人工检测价格便宜,且易于自动化。2.更好地了解藻类的频率和水样的某些化学性质以及其他特性(如季节、河流类型等)是如何相关的。3/10/20201:37AM40海藻数据•第一个数据集:包含200个水样。每条记录是一条河流在该年的同一季节的三个月内收集的水样的平均值。•其中,每条记录由11个变量组成。其中3个名义变量:水样收集的季节、河流大小和河水速度。8个变量是水样的不同化学参数:最大pH值、最小含氧量、平均氯化物含量、平均硝酸盐含量、平均氨含量、平均正磷酸盐含量、平均磷酸盐含量和平均叶绿素含量。与之相关的是7种不同的有害藻类的频率数目。•第二个数据集:140个不含7种藻类频率数目的测试集。3/10/20201:37AM51数据加载1.点击文件菜单下的改变工作目录来设定当前工作目录。2.输入以下命令把文件中的数据读入:algae-read.table('Analysis.txt',col.names=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后可通过拖拽的方式直接打开。3/10/20201:37AM62数据摘要•鉴于没有该问题领域足够的信息,首先了解一些数据的统计特性是一种较好的方式,它方便我们更好地理解问题。•获取数据统计特性的方法之一是获取其描述性的统计摘要。命令如下:summary(algae)•对于名义变量,它给出了每个可能取值的频数。•对于数值变量,它提供了均值、中位数、四分位数和极值等一系列统计信息。NA's表示缺失值的个数。•通过观察这些值,我们可以了解到数据分布的偏度和分散情况。3/10/20201:37AM73数据可视化(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为设置图形的标题。3/10/20201:37AM83数据可视化(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设置线型为虚线。3/10/20201:37AM94数据清理数据缺失的情形在实际问题中非常普遍。处理含有缺失值的数据时常用的几种策略:–将含有缺失值的案例剔除。–用中心趋势值来填补缺失值。–根据变量之间的相关关系填补缺失值。–根据案例之间的相似性填补缺失值。–使用能够处理缺失值数据的工具(见下一节)。3/10/20201:37AM104.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)3/10/20201:37AM114.1剔除缺失值(2)5.可根据4编写一个找出algae中含有给定数目缺失值的行。以下函数的功能是找出缺失值个数大于列数20%的行:library(DMwR)manyNAs(algae,0.2)第二个参数如不指定,默认为0.2,下面的命令与上一条等价:manyNAs(algae)6.我们利用上面的函数来剔除缺失值较多的记录:algae-algae[-manyNAs(algae),]这里第二个参数的默认值为0.2。3/10/20201:37AM124.2用中心趋势值填补缺失值代表中心趋势的值反映了变量分布的最常见值,这种方法也最自然、简便和快捷。对于接近正态的分布来说,均值是最佳选择;对偏态分布或有离群值的分布而言,中位数通常是更好的代表数据中心趋势的指标;对于名义变量,通常采用众数。可用以下函数完成填补所有缺失值:data(algae)algae-algae[-manyNAs(algae),]algae-centralImputation(algae)上述方法特别适用于大数据集,但是这种方法可能导致较大的数据偏差,影响后期的数据分析。但使用复杂的无偏方法寻找最佳数据填补对大型数据集可能也不适用。3/10/20201:37AM134.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)*,可用变量的相关性填补缺失值。*根据领域专家的解释,总的磷酸盐值包含正磷酸盐值。3/10/20201:37AM144.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,]3/10/20201:37AM154.4通过案例的相关性填补缺失值1.度量相似性时,最常用的指标是欧式距离。我们可通过使用这种度量来寻找与任何含有缺失值的案例最相似的10个水样,并用它们填补缺失值。方法一:简单计算这10个最近的案例的中位数并用中位数填补缺失值;若缺失值是名义变量则采用众数。方法二:采用这些最相似数据的加权均值。这里用高斯核函数从距离获得权重。命令如下:clean.algae-knnImputation(algae,k=10)2.这种方法看起来更合理,但也可能存在不相关的变量扭曲相似性,甚至造成大型数据集的计算特别复杂等问题。因此,填补缺失值时,大多应根据分析领域的知识来确定。3/10/20201:37AM165获取预测模型用于预测海藻的两种模型:多元线性回归模型和回归树模型。线性回归不能使用有缺失值的数据集,而回归树模型可以很自然地处理带缺失值的数据。•多元线性回归模型是最常用的统计数据分析方法,该方法给出了一个有关目标变量与一组解释变量关系的线性函数。由于多元线性回归模型中没有处理缺失值的方法,因此,我们可以做如下的数据预处理:data(algae)algae-algae[-manyNAs(algae),]clean.algae-knnImputation(algae,k=10)•这里还是先移除缺失值较多的记录,然后根据训练集数据个案的相似性来填补缺失值。3/10/20201:37AM175.1线性回归模型(1)建立用于预测海藻频率的线性回归模型:lm.a1-lm(a1~.,data=clean.algae[,1:12])函数lm()建立一个线性回归模型,其中,第一个参数给出了模型的函数形式。这个函数的形式是用数据中的其他所有变量来预测变量a1,第一个参数中的点“.”代表数据框中的所有除a1外的变量。如需要用变量mxPH和NH4来预测变量a1,就要定义模型为a1~mxPH+NH4。参数data是用来设定建模所用的数据集。3/10/20201:37AM185.1线性回归模型(2)通过下面的代码,我们可以获取更多线性模型的信息:summary(lm.a1)首先,给出数据拟合的残差(residuals)。残差应该是均值为0且为正态分布的。其次,对于每个多元线性回归方程的系数(变量),给出它的估计值和标准误差,并使用t检验来验证系数为0的假设检验。再者,给出模型与数据的吻合度,即模型所能解释的数据变差的比例。R-squared越接近于1说明模型拟合得越好,越小则代表模型拟合得越差。最后,还可以检验任何解释变量与目标变量的依赖关系。3/10/20201:37AM195.2回归树模型(1)•因模型解释的方差比例太低,才约32%,故实际验证结果表明:对海藻案例应用线性模型是不合适的。用线性思维去考虑非线性问题,得不到理想的结果。•我们考虑使用回归树预测。建立回归树:library(rpart)data(algae)algae-algae[-manyNAs(algae),]rt.a1-rpart(a1~.,data=algae[,1:12])•函数的形式是用数据中其他所有变量来预测a1,data是用来设定建模所用的数据集。3/10/20201:37AM205.2回归树模型(2)回归树rt.a1的图形表示的两种方法:plot(rt.a1)text(rt.a1)或prettyTree(rt.a1)建立回归树通常分两步。最初,生成一棵较大的树,然后通过统计估计删除底部的一些结点来对树进行修剪。这样是为了防止过度拟合。3/10/20201:37AM216模型评价和选择使用已有的训练数据获得模型的性能指标是不可靠的,因为这些计算是有偏的。实际上,有的模型可以很容易获得训练数据的零误差预测。然而,这一优秀性能很难推广到目标变量值未知的新样本上。这种现象我们通常称为过度拟合训练数据。我们需要一个模型,使它在未知数据上有可靠的预测性能。k折交叉验证是获得模型性能可靠估计的一种常用方法,它适用于小数据集。3/10/20201:37AM22k折交叉验证方法K折交叉验证:初始采样分割成K个子样