用R语言的MICE包对缺失数据进行多重插补（一）-- 缺失数据分析

在进行数据分析时，缺失数据是一个令人头痛的问题。数据缺失的原因五花八门，修补这些缺失数据的方法也是因情况而异。插补法（Imputation）是对一种对缺失数据进行调整的方法。该方法用多个可能的值来填充一个缺失的数据项，从而反映了缺失数据的不确定性。本例以R语言的MICE包为例，说明如何使用多重插补方法对缺失值进行估计。原文地址参见https://datascienceplus.com/imputing-missing-data-with-r-mice-package/，译文稍稍做了些改动。

如果缺失的数据量相对于数据集的大小而言非常小，那么丢掉少量具有缺失特征的样本可能是我们的最佳选择，这样我们就可以进行无偏分析。但是，在某些情况下，丢弃可用的数据点会减少原有数据中包含的信息。因此，我们需要对缺失数据进行填补和修复。
在某些情况下，我们可以用平均值代替缺失值来快速填充缺失数据，但是这类简单的方法通常会在数据中引入偏差。例如，在应用均值替换时，数据的平均值虽然保持不变，但是方差会减小，这可能不是我们所想要的。R语言中的MICE包可以帮助您使用合理的数据值来估算缺失值， 这些合理的值是根据缺失数据点所在的数据列中已有数据的分布估算出来的。
在这篇文章中，我们将使用空气质量数据集（R语言中提供的airquality数据包）来估算缺失值。首先，我将从数据集中删除一些数据点。

library(mice)  
data <- airquality
data[4:10,3] <- rep(NA,7)
data[1:5,4] <- NA

对分类变量而言，我们通常不对这类变量进行缺失值插补。虽然有人经常用观察值来替换缺失的分类变量，然而，这种做法的合理性需要进一步验证。 在airquality数据集中，我们没有分类变量。我们可以使用summary（）函数查看数据集的均值、分布及缺失值等概要信息。

data <- data[-c(5,6)]
summary(data)
     Ozone           Solar.R           Wind             Temp      
 Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :57.00  
 1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:73.00  
 Median : 31.50   Median :205.0   Median : 9.700   Median :79.00  
 Mean   : 42.13   Mean   :185.9   Mean   : 9.806   Mean   :78.28  
 3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00  
 Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00  
 NA's   :37       NA's   :7       NA's   :7        NA's   :5   

从上面的数据概要表中，我们可以看到Ozone（臭氧层）变量缺失最多，一共有37个样本缺失Ozone数据。下面，我们来看看缺失数据的模式。

数据缺失的机制
Little和Ruth(1987)把数据缺失的机制分为三类：

1. 完全随机缺失(missing completely at random, MCAR)
   所缺失的数据是完全随机的，缺失发生的概率既与已观察到的数据无关，也与未观察到的数据无关。这是一种比较理想的情况。
2. 随机缺失(missing at random, MAR)
   数据的缺失不是完全随机的。缺失数据发生的概率与所观察到的变量是有关的，而与未观察到的数据的特征是无关的。这是一个比较严重的问题，在这种情况下，我们需要进一步检查数据收集过程，并尝试了解数据为什么丢失。 例如，如果在一项问卷调查中，大多数人没有回答某个问题，他们为什么这么做，是问题不清楚吗？
3. 不可忽略的缺失(non-ignorable missing ,NIM)，亦称为非随机缺失(not missing at
   random, NMAR)，也有研究者将其称为MNAR(missing not at random)。
   缺失数据不仅依赖于其它变量，又依赖于变量本身，这种缺失即为不可忽略的缺失。

假设数据缺失是完全随机缺失（MCAR）。但是，如果数据缺失太多，也可能造成很大的问题。 一般来说，数据的缺失量小于数据总量的5％是可以接受的。 如果某个特征或样本的缺失数据量超过5％，那么您可考虑是否需要留下该特征或样本。 为此，我们可写一个简单的R函数检查某一特征（列）和样本（行）的缺失量是否超过5%，其代码如下所示：

pMiss <- function(x){sum(is.na(x))/length(x)*100}
apply(data,2,pMiss)
apply(data,1,pMiss)
运行结果如下：
    Ozone   Solar.R      Wind      Temp 
24.183007  4.575163  4.575163  3.267974 
  [1]  25  25  25  50 100  50  25  25  25  50  25   0   0   0   0
 [16]   0   0   0   0   0   0   0   0   0  25  25  50   0   0   0
 [31]   0  25  25  25  25  25  25   0  25   0   0  25  25   0  25

从上表我们看到，Ozone这列数据大约缺失了25％的数据点，因此我们可能会考虑在接下来的数据分析中剔除该数据列，如果要保留该数据列，我们可能需要收集更多数据。而其他三个变量的缺失值比例均低于5％，所以我们可以保留它们。就样本而言，由于每个样本有四个变量，如果某个样本缺失一个变量的值，那么该样本缺失的数据量为25％；如果缺失两个变量的值，那么该样本缺失的数据量为50％，这时，如果可能的话，我们应该剔除该样本。

利用R语言的Mice包来查看数据的缺失模式

R语言的Mice包提供了一个函数md.pattern（），该函数可以帮助我们更好地了解缺少数据的模式。

md.pattern(data)
    Temp Solar.R Wind Ozone   
104    1       1    1     1  0
 34    1       1    1     0  1
  4    1       0    1     1  1
  3    1       1    0     1  1
  3    0       1    1     1  1
  1    1       0    1     0  2
  1    1       1    0     0  2
  1    1       0    0     1  2
  1    0       1    0     1  2
  1    0       0    0     0  4
       5       7    7    37 56

上述函数的输出结果告诉我们，有104个样本是完整的，34个样本仅在Ozone数据列缺失一个数据列，有4个样本在Solar.R数据列缺失数据。
使用VIM包可以帮助我们对数据的缺失模式进行可视化，代码及可视化的结果如下所示：

library(VIM)
aggr_plot <- aggr(data, col=c('navyblue','red'), 
                  numbers=TRUE, 
                  sortVars=TRUE, 
                  labels=names(data), 
                  cex.axis=.7, 
                  gap=3, 
                  ylab=c("数据缺失模式直方图","模式"))

从上图我们可以直观地看到，近70％的样本没有缺失任何信息，大约22％的样本缺失Ozone数据列。 另一个有用的可视化方法是下面显示的盒状图：

marginplot(data[c(1,2)])

上图只显示了两个变量（Ozone和Solar.R）的缺失数据项，及这两个数据项之间的交汇图。在这张图上，左边的红色矩形代表了缺失Ozone数据项，但是Solar.R数据项没有缺失的样本分布，而左边的蓝色矩形则显示了剩余的数据点的分布。类似的，该图下面的红色矩形表示缺失Solar.R数据项但是没有缺失Ozone数据项的样本点。如果我们的缺失数据缺失服从完全随机分布（MCAR），那么我们在图上的红色矩形和蓝色矩形应该看起来很相似。

在对缺失数据分析后，我们就可以对它们进行插补了，详细步骤请见下一篇。