预测海藻数量(R语言)

1. 问题描述与目标
我希望通过建立预测模型预测河流中有害海藻的数量。本案例的目的是更好地了解影响藻类频率的因素。也就是说，我们要了解藻类的频率和水样的某些化学性质以及其他特征(如季节、河流类型等)是如何相关的。

2. 数据说明
有两个数据集，第一个数据集有200个水样。该数据集的每一条记录是同一条河流在该年的同一个季节的三个月内收集的水样的平均值。
每条记录由11个变量构成。其中3个变量是名义变量，它们分别描述水样收集的季节、收集河流的大小和河水速度。余下的8个变量是所观察水样的不同化学参数，即最大pH值、最小含氧量（O2）、平均氯化物含量(cl)、平均硝酸盐含量(NO3)、平均氨含量(NH4)、平均正磷酸含量(PO4)、平均磷酸盐含量(PO4)、平均叶绿素含量。与这些参数相关的是7种不同有害藻类在相应水样中的频率数目。并未提供所观察藻类的名称的有关信息。

第二个数据集由140个额外观察值构成。它们的基本结构和第一个数据集一样，但是它不包含7种藻类的频率数目。本案例的主要目的是预测140个水样中7种藻类的频率。

在这种问题中，任务是建立预测模型，并预测在给定预测变量的取值时相应的目标变量的值。预测模型也可能会说明哪一个预测变量对目标变量有较大的影响，即模型可能提供影响目标变量因素的一个综合描述。

3. 数据加载到R
两种方法：1)利用提供的R添加包；2)下载相应数据的文本文件，然后把文件读入到R中

第一种方法：安装R添加包，该添加包提供了案例的数据集和函数，代码如下：install.packages(‘DMwR’) 只要载入了R添加包，就直接有了一个名为algae的数据框。这个数据框含有前面提到的200个观测值:

>library(DMwR)
>head(algae)

数据框可以看做一种有列名称的矩阵或者表格。函数head()将显示数据框的前6行。

第二种方法：先从网站http://www.dcc.fc.up.pt/~ltorgo/DataMiningWithR/下载三个文本文件。“Analysis.txt”中包含200个水样，而”Eval.txt”中包含140个测试样本。“Sols.txt”包含140个测试样本的藻类频率。将文本文件读入到R中的代码如下：

>algae <- read.table('Analysis.txt',
          header=F,
          dec='.',
          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'))

参数header = F表示要读的文件的第一行不包括变量名。dec=’.’指出数值使用’.’字符分隔小数位。这两个参数设置可以省略，因为这里使用了它们的默认值。

4. 数据缺失
当我们处理含有缺失值的数据时，可以有几种处理策略：
1.将含有缺失值的案例剔除
2.根据变量之间的相关关系填补缺失值
3.根据案例之间的相似性填补缺失值
4.使用能够处理缺失值数据的工具

4.1 将缺失部分剔除
剔除含有缺失数据的记录非常容易实现，尤其当这些记录所占的比例在可用数据集中非常小的时候，这个选择就比较合理。

在剔除某些变量中至少含有一个缺失数据的所有观测值时，最好先检查观测值，或者至少得到这些观测值的个数，例如：

>algae[!complete.cases(algae),]

函数complete.cases()产生一个布尔值向量，该向量的元素个数与algae数据框中的行数相同，如果数据库的相应行中不含NA值(即为一个完整的观测值)，函数返回的值就是TURE。

当然我们也可以用如下代码统计不完整数据的个数：

>nrow(algae[!complete.cases(algae),])
[1] 16

除了统计个数之外，我们还可以删除这16条不完整数据的样本:

>algae <- na.omit(algae)

几个处理NA数据很有用的函数：
na.fail 返回不包含缺失值记录的数据
na.omit 返回剔除了缺失值记录的数据
na.pass 返回没有作任何处理的数据

即使我们不适用剔除所有包含缺失值记录的极端方法，我们也可以剔除某些缺失值太多的样本，这些样本几乎是无用样本，如果采用复杂的方法来填补缺失值，就会导致较大偏差。我们观察这16个缺失值记录的数据，可以看到62条和199条记录中的11个解释变量有6个是缺失值。

    season   size  speed mxPH mnO2    Cl   NO3 NH4    oPO4     PO4  Chla   a1   a2  a3   a4  a5  a6  a7

62  summer  small medium 6.40   NA    NA    NA  NA      NA  14.000    NA 19.4  0.0 0.0  2.0 0.0 3.9 1.7

199 winter  large medium 8.00  7.6    NA    NA  NA      NA      NA    NA  0.0 12.5 3.7  1.0 0.0 0.0 4.9

所以我们可以剔除它们：

>algae <-algae[-c(62,199),]

在有些问题中，由于大量记录中含有缺失值，用上面的观察方法来检查数据的缺失值是不可行的，我们需要找出缺失值较多的样本所在的行。下面的代码可以帮助我们找出海藻数据集中每行数据的缺失值的个数：

apply(algae,1,function(x) sum(is.na(x)))

apply()函数，它可以把任何其他函数应用到一个多维对象的各个维度上。这个被应用的函数在apply()函数的第三个参数中给出，对数据框的每一行都分别调用该函数。在这里，临时函数的功能是计算对象x中NA的数量。在R中逻辑值TURE等于数值1，逻辑值FALSE等于数值0，这意味着当加一个布尔值向量时，得到向量中TURE的元素的个数。

根据以上代码，可以找出algae中含有给定数目缺失值的行：

>data(algae)
>manyNAs(algae,0.2)

函数manyNAs()的功能是找出缺失值个数大于列数20%的行。函数第二个参数的默认值为0.2。

用下面的代码就无须知道含有缺失值较多的行的具体数量：

>algae <- algae[-manyNAs(algae),]

4.2 用最高频率值来填补缺失值
填补缺失数据最简便和便捷的方法是使用一些代表中心趋势的值。代表中心趋势的值反映了变量分布的最常见值。有多个代表数据中心趋势的指标，例如平均值、中位数、众数等。最适合的选择由变量的分布决定。
对于接近正态的分布来说，所有的观测值都较好地聚集在平均值周围，平均值数就是最佳选择。然而，对于偏态分布，或者有离群值的变量来说，选择平均值就不好。因此，在对变量分布进行检查之前选择平均值作为中心趋势的代表是不明智的。

1. 使用平均值来填补缺失值，代码如下：

>algae[48,'mxPH'] <- mean(algae$mxPH,na.rm=T)

对样本[48]中变量mxPHd的缺失值用平均值填补，因为该变量分布近似正态分布。这里函数mean()计算数值向量的平均值，参数na.rm = T使计算时忽略缺失数据NA.

1. 大多数时候采用一次填补一列中的所有缺失值而不是像上面那样一行一行地逐个填补。以Chla为例，这个变量在第12行上有缺失值。并且Chla的分布偏向于较低的数值，并且它有几个极端值，这些使得平均值不能代表大多数的变量值。因此我们使用中位数来填补这一类的缺失值：

>algae[is.na(algae$Chla),'Chla'] <- median(algae$Chla,na.rm=T)

添加包中提供的函数centralImputation()可以用数据的中心趋势值来填补数据集的所有缺失值。对于数值型变量，该函数用中位数；对于名义变量，它采用众数，代码如下：

>data(algae)
>algae <- algae[-manyNAs(algae),]
>algae <- centralImputation(algae)

第二行代码先去掉数据缺失严重的样本62和199
第三行用centalImputaion()填补缺失值

由于缺失值的存在会导致某些方法不能使用，所有使用上面的方法填补缺失值通常也认为不是很好的方法。虽然上述的简单方法速度快，特别适用于大数据集，但是它可能导致较大的数据偏差，影响后期的数据分析工作。然而，使用无偏方法来寻找最佳数据填补值复杂，对于大型数据挖掘问题可能并不适用。

4.3 用过变量的相关关系来填补缺失值
另一种获取缺失值较少偏差估计值的方法是探寻变量之间的相关关系。如下命令可以得到变量间的相关值：

>cor(algae[,4:18],use="complete.obs")

函数cor()的功能是产生变量之间的相关值矩阵(前3个变量是名义变量，所以计算相关值时不考虑它们)。设定参数use = “complete.obs”时，R在计算相关值时忽略含有NA的记录。相关值在1(或-1)周围表示相应的两个变量之间有强正(或负)线性相关关系。然后利用其他R函数可以得到变量间线性相关的近似函数形式，它可以让我们通过一个已知变量的值计算出另一个未知变量的值。

            mxPH        mnO2          Cl         NO3         NH4         oPO4         PO4        Chla          a1
mxPH  1.00000000 -0.10269374  0.14709539 -0.17213024 -0.15429757  0.090229085  0.10132957  0.43182377 -0.16262986

函数cor()的输出结果并不是很清晰，但可以通过symnum()来改善结果的输出形式。

>symnum(cor(algae[,4:18],use="complete.obs"))

我们发现：变量NH4和NO3之间，变量PO4和oPO4之间的相关性比较强。

开始查找变量之间的相关性之前，我们要剔除样本62和样本199，因为它们有太多变量含有缺失值。样本中的变量NH4和NO3就没有缺失值了。至于变量PO4和oPO4,它们之间的相关性可以帮助填补这两个变量的缺失值。

我们需要找到这两个变量之间的线性相关关系，代码如下：

>data(algae)
>algae <- algae[-manyNAs(algae),]
>lm(PO4 ~ oPO4,data=algae)

Call:
lm(formula = PO4 ~ oPO4, data = algae)

Coefficients:
(Intercept)         oPO4  
     42.897        1.293  

函数lm()可以用来获取形如$Y=\beta_0+\beta_1x_1+...+\beta_nx_n$ 的线性模型。如果这两个变量不是同时有缺失值，那么可以通过这个公式计算这些变量的缺失值。线性模型是：$PO4 = 42.897+1.293*oPO4$。

样本28在变量PO4上有缺失值，可以使用上面的线性关系计算缺失值：

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

假设变量有多个缺失值时，我们可以通过构造函数来解决问题：