在分析数据集时，常常会碰到一些缺失值，如果缺失值的数量相对总体来说非常小，那么直接删除缺失值就是一种可行的方法。但某些情况下，直接删除缺失值可能会损失一些有用信息，此时就需要寻找方法来补全缺失值。今天小编给大家介绍一个用来处理缺失值的 R 包——MICE，本文为译文，原文链接[1]及参考文章[2]见文末。
小编在原文的基础上找到了一种确定最佳插补集的方法，文章有点长，但是干货满满，希望大家耐心阅读呀～
mice 简介
mice包帮助我们用可信的数据值来填补缺失值，这些可信的数据值是根据原始数据分布特征得到的。该包为多元缺失数据创建多个输入(替换值)，其中每个不完全变量由一个单独的模型输入。MICE 算法支持输入的数据类型有：连续的、二值的、无序分类和有序分类数据。
数据处理
本文，我们将使用 R 自带的一个空气质量数据集airquality来估算缺失的值。为了介绍 mice 包的用法，先从数据集中删除一些数据点，制造一个缺失数据集。
data <- airquality data[4:10,3] <- rep(NA,7) data[1:5,4] <- NA
使用summary()查看数据特征，可以看出 Ozone (臭氧)是缺失数据最多的变量。下面我们将深入挖掘缺失数据的特征。
> summary(data) Ozone Solar.R Wind Temp Month Day Min. : 1.00 Min. : 7.0 Min. : 1.700 Min. :57.00 Min. :5.000 Min. : 1.0 1st Qu.: 18.00 1st Qu.:115.8 1st Qu.: 7.400 1st Qu.:73.00 1st Qu.:6.000 1st Qu.: 8.0 Median : 31.50 Median :205.0 Median : 9.700 Median :79.00 Median :7.000 Median :16.0 Mean : 42.13 Mean :185.9 Mean : 9.806 Mean :78.28 Mean :6.993 Mean :15.8 3rd Qu.: 63.25 3rd Qu.:258.8 3rd Qu.:11.500 3rd Qu.:85.00 3rd Qu.:8.000 3rd Qu.:23.0 Max. :168.00 Max. :334.0 Max. :20.700 Max. :97.00 Max. :9.000 Max. :31.0 NA's :37 NA's :7 NA's :7 NA's :5
确定缺失数据的类型
缺失数据有两大类别：
- MCAR : missing completely at random，数据是完全随机缺失的。这是缺失数据的理想场景。
- MNAR : missing not at random ，数据不是随机缺失的。这种情况非常严重，此时需要检查数据收集过程并试图找出造成数据缺失的环节。
假设数据是 MCAR ，那么缺失值过多也可能是个问题。对于大型数据集，通常安全的最大阈值为总阈值的 5% 。 如果某个样本（或特征）的缺失数据量超过5%，可以考虑删除该样本（或特征）。因此，我们建立一个简单的函数pMiss()检查是否有超过 5% 缺失值的特征（列）和样本（行）:
pMiss <- function(x){round(sum(is.na(x))/length(x),3)}
将该函数分别应用在行和列上：
> apply(data, 1, pMiss) [1] 0.167 0.167 0.167 0.333 0.667 0.333 0.167 0.167 0.167 0.333 0.167 0.000 0.000 0.000 0.000 [16] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.167 0.167 0.333 0.000 0.000 0.000 [31] 0.000 0.167 0.167 0.167 0.167 0.167 0.167 0.000 0.167 0.000 0.000 0.167 0.167 0.000 0.167 [46] 0.167 0.000 0.000 0.000 0.000 0.000 0.167 0.167 0.167 0.167 0.167 0.167 0.167 0.167 0.167 [61] 0.167 0.000 0.000 0.000 0.167 0.000 0.000 0.000 0.000 0.000 0.000 0.167 0.000 0.000 0.167 [76] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.167 0.167 0.000 0.000 0.000 0.000 0.000 0.000 [91] 0.000 0.000 0.000 0.000 0.000 0.167 0.167 0.167 0.000 0.000 0.000 0.167 0.167 0.000 0.000 [106] 0.000 0.167 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.167 0.000 0.000 0.000 0.167 0.000 [121] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 [136] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.167 [151] 0.000 0.000 0.000
就样本（行）而言（不考虑 Month 和 Day 特征，因为没有缺失值），只要缺失一个特征，每个样本就会丢失25%的数据。如果样本缺少2个或更多特征(>50%)，应尽可能删除该样本。
> apply(data, 2, pMiss) Ozone Solar.R Wind Temp Month Day 0.242 0.046 0.046 0.033 0.000 0.000
可以看到 Ozone （臭氧）缺失了近 25% 的数据，因此我们可以考虑删除它，或者收集更多的测量数据。其他变量低于 5% 的阈值可以保留。
使用 mice 包寻找缺失数据的特征
mice 包提供了一个很好的函数md.pattern() 来寻找缺失值的特征。
> library(mice) > md.pattern(data) Month Day Temp Solar.R Wind Ozone 104 1 1 1 1 1 1 0 34 1 1 1 1 1 0 1 3 1 1 1 1 0 1 1 1 1 1 1 1 0 0 2 4 1 1 1 0 1 1 1 1 1 1 1 0 1 0 2 1 1 1 1 0 0 1 2 3 1 1 0 1 1 1 1 1 1 1 0 1 0 1 2 1 1 1 0 0 0 0 4 0 0 5 7 7 37 56
输出结果表明有 104 个样本没有缺失值，有 34 个样本只缺失了 Ozone 的数据，以此类推。
- 使用
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("Histogram of missing data","Pattern"))
这幅图帮助我们了解到，大约 68% 的样本没有丢失任何信息，22% 的样本丢失了 Ozone 值，剩下的样本也存在其他缺失的形式。
- 使用边际箱线图
marginplot()可视化
marginplot(data[c(1,2)])
虽然一次只能绘制两个变量，但是也能获得一些有用信息。左边的红箱显示了缺失 Ozone 的 Solar.R 的分布，蓝箱表示剩余数据点的分布。底部的红箱显示了缺失 Solar.R 的 Ozone 分布。如果我们对数据是 MCAR 类型的假设是正确的，那么红箱和蓝箱将非常相似。
填补缺失值
这里用到的是 mice()函数，所需的主要参数如下：
data：包含缺失值的数据框或矩阵。缺失值被编码为 NA。
m：多重插补法的数量，默认为 5。
method：指定数据中每一列的输入方法。1）数值型数据适用 pmm；2）二分类数据适用 logreg；3）无序多类别数据适用 ployreg；4）有序多分类变量适用 polr。默认方法为 pmm 。
maxit：迭代次数，默认为 5 。
seed：一个整数，由set.seed()产生，作为偏移随机数生成器。默认情况下，不使用随机数生成器。
针对本文所用的数据集进行插补:
> tempData <- mice(data,m=5,maxit=50,meth='pmm',seed=500) > summary(tempData) Class: mids Number of multiple imputations: 5 Imputation methods: Ozone Solar.R Wind Temp Month Day "pmm" "pmm" "pmm" "pmm" "" "" PredictorMatrix: Ozone Solar.R Wind Temp Month Day Ozone 0 1 1 1 1 1 Solar.R 1 0 1 1 1 1 Wind 1 1 0 1 1 1 Temp 1 1 1 0 1 1 Month 1 1 1 1 0 1 Day 1 1 1 1 1 0
结果中的 PredictorMatrix 是预测变量矩阵，行代表插补变量，列代表为插补提供信息的变量，1和0表示使用和未使用。以第一行为例，Ozone存在缺失值，并利用了其他五个变量的信息来进行数值插补。
如果需要查看某个变量的插补值，可以使用下面的语句：
> tempData$imp$Ozone 1 2 3 4 5 5 21 36 20 45 41 10 20 7 7 4 18 25 14 14 12 18 18 26 8 32 37 32 1 27 12 13 11 32 20 32 32 20 63 59 30 ......
右边第一列的数字 5 表示 Ozone 变量的第 5 个观测值，这一行的其他五个数值表示每次插值后的结果。
完成插补后，接下来可以使用complete()函数返回完整的数据集，action的参数值表示选择第几次的插补值来填补原始数据集。
completedData <- complete(tempData,action = 5)
查看原始数据和插补值的分布情况
这里使用密度图来查看，蓝色线代表原始数据，红色线代表每一次插补得到的数据。
densityplot(tempData)
选择合适的插补值
从上面的密度图可以发现多重插补后的数据拟合度有好有坏，如果需要利用数据建模，那么必然要选择拟合效果最好的一个插补值。目前，小编还未发现 mice 包能直接提取最佳插补数据集的方法，只能间接通过图片比较，如果同学有更好的方法欢迎在评论区留言，让小编和其他同学可以一起学习～
下面的语句可以看到每一次插补所得到的数据集，蓝色是原始数据集，第一幅红色图为一重插补，以此类推。
densityplot(tempData,~ Ozone + Solar.R + Wind + Temp | .imp)
上述图形，没有将插补后的数据与原始数据比较，可以采用下面语句先在左边图形（该图形包含了全部插值结果）中确定拟合度最好的线，然后通过修改imp的值，直到在右侧图形中找到那条线。
# 以 Ozone 为例 densityplot(tempData,~ Ozone | .imp == c(1,2))
确定好最佳插补重数后，利用上文提到的complete（）函数就可以得到最终数据集啦。
参考资料
[1]
原文链接: https://datascienceplus.com/imputing-missing-data-with-r-mice-package/
[2]
参考文章: https://zhuanlan.zhihu.com/p/36436260
