R语言实战（18）—处理缺失数据的高级方法
引言：上一章我们学习了一系列用于二分类的机器学习方法，包括逻辑回归分类方法、传统决策树、条件推断树、集成性的随机森林以及支持向量机。这一期我们就来学习如何处理缺失数据吧。
后台回复“R语言实战”即可获取二维码加入R语言实战学习讨论群。
本章中，我们将学习处理缺失数据的传统方法和现代方法，主要使用 VIM 和 mice 包。数据来源：VIM 包提供的哺乳动物睡眠数据sleep，该数据研究了62种哺乳动物的睡眠变量（因变量）、生态学变量（自变量）和体质变量间的关系（自变量）。
install.packages(c("VIM","mice"))18.1 处理缺失值的步骤
(1) 识别缺失数据；
(2) 检查导致数据缺失的原因；
(3) 删除包含缺失值的实例或用合理的数值代替（插补）缺失值。
图18-1 处理不完整数据的方法，以及R中相关的包和函数
要完整介绍处理缺失数据的方法，用一本书的篇幅才能做到。本章，我们只是学习探究缺失值模式的方法，并重点介绍三种最流行的处理不完整数据的方法（推理法、行删除法和多重插补法）。
18.2 识别缺失值
背景知识：
NA （不可得）代表缺失值，
NaN （不是一个数）代表不可能值。
符号 Inf 和 Inf 分别代表正无穷和负无穷。
函数 is.na() 、 is.nan() 和is.infinite() 可分别用来识别缺失值、不可能值和无穷值。
表181 is.na() 、 is.nan() 和 is.infinite() 函数的返回值示例
识别具体的缺失值 is.na() 、 is.nan()
矩阵或数据框中没有缺失值的行 函数 complete.cases()+ sum() 和 mean() 函数
举例：
# 例子1y <- c(1, 2, 3, NA)is.na(y)c(FALSE, FALSE,FALSE, TRUE)# 例子2# 加载数据集data(sleep, package="VIM")# 列出没有缺失值的行sleep[complete.cases(sleep),]# 列出有一个或多个缺失值的行sleep[!complete.cases(sleep),]#逻辑值 TRUE 和 FALSE 分别等价于数值1和0，用sum和mean函数简化结果。> sum(is.na(sleep$Dream))#12个缺失值[1] 12> mean(is.na(sleep$Dream))#19%的实例在此变量上有缺失值[1] 0.19> mean(!complete.cases(sleep))#数据集中32%的实例包含一个或多个缺失值[1] 0.32
notes:
第一， complete.cases() 函数仅将 NA 和 NaN 识别为缺失值，无穷值（ Inf 和 Inf ）被当作有效值。第二，必须使用与本章中类似的缺失值函数来识别R数据对象中的缺失值。像 myvar == NA 这样的逻辑比较无法实现。
18.3 探索缺失值模式
18.3.1 列表显示缺失值
mice 包中的 md.pattern() 函数可生成一个以矩阵或数据框形式展示缺失值模式的表格.
library(mice)data(sleep, package="VIM")md.pattern(sleep)BodyWgt BrainWgt Pred Exp Danger Sleep Span Gest Dream NonD42 1 1 1 1 1 1 1 1 1 1 011 1 10 0 22 1 1 1 1 1 0 1 1 1 0 21 1 1 1 1 1 1 0 0 1 1 21 1 1 0 1 1 0 0 31 1 1 1 1 0 1 0 0 30 0 0 0 0 4 4 4 1 2 14 38
解读：0表示变量的列中有缺失值，1则表示没有缺失值。第一行表述了“无缺失值”的模式（所有元素都为1）。第二行表述了“除了 Span 之外无缺失值”的模式。第一列表示各缺失值模式的实例个数，最后一列表示各模式中有缺失值的变量的个数。此处可以看到，有42个实例没有缺失值，仅2个实例缺失了 Span 。9个实例同时缺失了 NonD 和 Dream的值。数据集包含了总共(42×0)+(2×1)+…+(1×3)=38个缺失值。最后一行给出了每个变量中缺失值的数目。
18.3.2 图形探究缺失数据
VIM包很多类似的函数，本节我们主要学习三种 aggr() 、matrixplot() 和 scattMiss()
aggr()
library("VIM")aggr(sleep, prop=FALSE, numbers=TRUE)aggr(sleep, prop=TRUE, numbers=TRUE)
图18-2 aggr() 生成的 sleep 数据集的缺失值模式图形
matrixplot()
matrixplot(sleep)图18-3 sleep 数据集按实例（行）展示真实值和缺失值的矩阵图。矩阵按 BodyWgt重排。
marginplot() 函数可生成一幅散点图，在图形边界展示两个变量的缺失值信息。
marginplot(sleep[c("Gest","Dream")], pch=c(20),col=c("darkgray", "red", "blue"))
图18-4 做梦时长与妊娠期时长的散点图，边界展示了缺失数据的信息
scattMiss()
18.3.3 用相关性探索缺失值
用指示变量（1表示缺失，0表示存在）替代数据集中的缺失数据，生成更的矩阵有时被称作影子矩阵，通过该方法可以简化探索不同缺失变量之间的关系。
> x <- as.data.frame(abs(is.na(sleep)))#将缺失值转化为为1和0并赋予给新的矩阵中head(sleep, n=5)#观察原数据的前几行BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger1 6654.000 5712.0 NA NA 3.3 38.6 645 3 5 32 1.000 6.6 6.3 2.0 8.3 4.5 42 3 1 33 3.385 44.5 NA NA 12.5 14.0 60 1 1 1NA 25 5 2 35 4head(x, n=5)#转换后的数据集BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger1 0 0 1 1 0 0 0 0 0 00000#提取含（但不全部是）缺失值的变量y <- x[which(apply(x,2,sum)>0)]#列出这些指示变量间的相关系数cor(y)Sleep Span Gest0.486 0.015 -0.142Dream 0.907 1.000 0.204 0.038 -0.129Sleep 0.486 0.204 1.000 -0.069 -0.069Span 0.015 0.038 -0.069 1.000 0.198Gest -0.142 -0.129 -0.069 0.198 1.000# Dream 和 NonD 常常一起缺失（r=0.91）。相对可能性较小的是 Sleep 和 NonD 一起缺失=0.49），以及 Sleep 和 Dream （r=0.20）#含缺失值变量与其他可观测变量间的关系cor(sleep, y, use="pairwise.complete.obs")NonD Dream Sleep Span GestBodyWgt 0.227 0.223 0.0017 -0.058 -0.054BrainWgt 0.179 0.163 0.0079 -0.079 -0.073NA NA NA -0.043 -0.046NA -0.189 0.117 0.228Sleep -0.080 -0.080 NA 0.096 0.040-0.065Gest 0.202 0.051 0.1597 -0.175 NA0.048 -0.068 0.2025 0.023 -0.201Exp 0.245 0.127 0.2608 -0.193 -0.1930.2089 -0.067 -0.204Warning message:In cor(sleep, y, use = "pairwise.complete.obs") :the standard deviation is zero
18.4 理解缺失数据的来由和影响
识别缺失数据的数目、分布和模式有两个目的：(1) 分析生成缺失数据的潜在机制；(2) 评价缺失数据对回答实质性问题的影响，以便判断哪种统计方法最适合用来分析你的数据。
例如我们想知道：
缺失数据的比例多大？
缺失数据是否集中在少数几个变量上，抑或广泛存在？
缺失是随机产生的吗？
缺失数据间的相关性或与可观测数据间的相关性，是否可以表明产生缺失值的机制？等等。
如果是不太重要的不太重要的变量上，可以删除，然后再进行正常的数据分析。如果有一小部分数据（如小于10%）随机分布在整个数据集中（MCAR），那么我们可以分析数据完整的实例。如果可以假定数据是MCAR或者MAR，可以应用多重插补法来获得有效的结论。如果数据是NMAR，则需要借助专门的方法。
18.5 理性处理不完整数据方法一
当数据存在冗余信息或有外部信息可用时，推理法可用来恢复缺失值。
推理方法会根据变量间的数学或者逻辑关系来填补或恢复缺失值。
举例：
1、在 sleep 数据集中，变量 Sleep 是 Dream 和 NonD 变量的和。若知道了它们中的任意两个变量，你便可以推导出第三个。
2、我们需要了解各代群体（依据出生年代区分，如沉默的一代、婴儿潮一代、婴儿潮后期一代、无名一代、千禧一代）在工作与生活间的平衡差异。调查对象都被问及了他们的出生日期和年龄，如果出生日期缺失，你便可以根据他们的年龄和其完成调查时的日期来填补他们的出生年份（以及他们所属的年代群体），这样便可使调查问卷完整。
3、推理研究法常常需要创造性和想法，同时还需要许多数据处理技巧，而且数据的恢复可能是准确的（如睡眠的例子）或者近似的（性别的例子）。下一节我们将探究一种通过删除观测来创建完整数据集的方法。
18.6 完整实例分析（行删除）方法二
当数据是MCAR，后续样本量的减少对统计检验效力不会造成很严重的影响时，行删除法非常有用。
行删除法假定数据MCAR（即完整的观测只是全数据集的一个随机子样本）的前提下应用的。
2个主要的函数：na.omit 函数和 complete.cases()函数
# mydata 中所有包含缺失数据的行都被删除，把结果存储到newdata 中newdata <- mydata[complete.cases(mydata),]newdata <- na.omit(mydata)#例子：用行删除法处理数据后再计算相关系数，探索睡眠研究中变量间的关系options(digits=1)cor(na.omit(sleep))cor(sleep,use="complete.obs")#可得出一样的数据BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp DangerBodyWgt 1.00 0.96 -0.4 -0.07 -0.3 0.47 0.71 0.10 0.4 0.260.3 0.15-0.53Dream -0.07 -0.07 0.5 1.00 0.7 -0.27 -0.41 -0.40 -0.5 -0.57Sleep -0.34 -0.34 1.0 0.72 1.0 -0.38 -0.61 -0.40 -0.6 -0.60Span 0.47 0.63 -0.4 -0.27 -0.4 1.00 0.65 -0.17 0.3 0.01Gest 0.71 0.73 -0.6 -0.41 -0.6 0.65 1.00 0.09 0.6 0.31Pred 0.10 -0.02 -0.4 -0.40 -0.4 -0.17 0.09 1.00 0.6 0.930.41 0.32 -0.6 -0.50 -0.6 0.32 0.57 0.63 1.0 0.790.26 0.15 -0.5 -0.57 -0.6 0.01 0.31 0.93 0.8 1.00#研究寿命和妊娠期对睡眠中做梦时长的影响，可应用行删除法的线性回归fit <- lm(Dream ~ Span + Gest, data=na.omit(sleep))summary(fit)Call:= Dream ~ Span + Gest, data = na.omit(sleep))Residuals:Min 1Q Median 3Q Max-0.915 -0.221 0.382 4.183CoefficientsEstimate Std. Error t value Pr(>|t|)2.480122 0.298476 8.31 3.7e-10 ***0.013130 -0.04 0.9710.002081 -2.11 0.041 *---codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 1 on 39 degrees of freedomMultiple R-squared: 0.167, Adjusted R-squared: 0.125: 3.92 on 2 and 39 DF, p-value: 0.0282#解读：动物妊娠期越短，做梦时长越长（控制寿命不变）；而控制妊娠期不变时，寿命与做梦时长不相关。整个分析基于有完整数据的42个实例。如果 data=na.omit(sleep) 被 data=sleep替换，m() 将使用有限的行删除法定义。只有用函数拟合的、含缺失值的变量（本例是 Dream 、Span 和 Gest ）对应的实例才会被删除，这时数据分析将基于44个实例。
行删除法假定数据MCAR（即完整的观测只是全数据集的一个随机子样本）。此例中，我们假定42种动物是62种动物的一个随机子样本。如果违反了MCAR假设，回归参数的结果将是有偏的，行删除法由于减少了样本数量，统计效率会下降，比如此例中就减少了32%的样本量。接下来，我们将探讨一种能够利用整个数据集的方法（可以囊括那些含缺失值的观测）。
18.7 多重插补方法三
当你认为数据是MCAR或MAR，并且缺失数据问题非常复杂时，多重插补将是一个非常实用方法。
多重插补（MI）是一种基于重复模拟的处理缺失值的方法。
本章主要介绍了 mice 包提供的多重插补法（MI）。
mice 包的插补分析过程：
library(mice)
imp < mice(data, m)
fit < with(imp, analysis)
pooled < pool(fit)
summary(pooled)
其中，
data 是一个包含缺失值的矩阵或数据框。
imp 是一个包含m个插补数据集的列表对象，同时还含有完成插补过程的信息。默认m为5。
analysis 是一个表达式对象，用来设定应用于m个插补数据集的统计分析方法。方法包括做线性回归模型的 lm() 函数、做广义线性模型的 glm() 函数、做广义可加模型的gam() ，以及做负二项模型的 nbrm() 函数。表达式在函数的括号中， ~ 的左边是响应变量，右边是预测变量（用 + 符号分隔开）。
fit 是一个包含m个单独统计分析结果的列表对象。
pooled 是一个包含这m个统计分析平均结果的列表对象。
将多重插补法应用到 sleep 数据集，如下所示
library(mice)data(sleep, package="VIM")imp <- mice(sleep, seed=1234)[...output deleted to save space...]fit <- with(imp, lm(Dream ~ Span + Gest))pooled <- pool(fit)summary(pooled)est se t df Pr(>|t|) lo 95(Intercept) 2.58858 0.27552 9.395 52.1 8.34e-13 2.03576Span -0.00276 0.01295 -0.213 52.9 8.32e-01 -0.02874Gest -0.00421 0.00157 -2.671 55.6 9.91e-03 -0.00736hi 95 nmis fmi(Intercept) 3.14141 NA 0.0870Span 0.02322 4 0.0806Gest -0.00105 4 0.0537
我们进一步可以通过检查分析过程所创建的对象来获取更多的插补信息，如 imp 对象的汇总信息：
impMultiply imputed data setCall:= sleep, seed = 1234)Number of multiple imputations: 5Missing cells per column:BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred0 0 14 12 4 4 4 0Exp Danger0 0Imputation methods:BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred"" "pmm" "pmm" "pmm" "pmm" "pmm" ""Exp Danger""VisitSequence:NonD Dream Sleep Span Gest3 4 5 6 7PredictorMatrix:BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger0 0 0 0 0 001111 1 1Gest 1 1 1 1 1 1 0 1 1 10 0 0 0Exp 0 0 0 0 0 0 0 0 0 00 0 0Random generator seed value: 1234
再而，提取 imp 对象的子成分，可以观测到实际的插补值：
imp$imp$Dream1 2 3 4 51 0.5 0.5 0.5 0.5 0.03 2.3 2.4 1.9 1.5 2.44 1.2 1.3 5.6 2.3 1.314 0.6 1.0 0.0 0.3 0.524 1.2 1.0 5.6 1.0 6.626 1.9 6.6 0.9 2.2 2.030 1.0 1.2 2.6 2.3 1.431 5.6 0.5 1.2 0.5 1.447 0.7 0.6 1.4 1.8 3.653 0.7 0.5 0.7 0.5 0.555 0.5 2.4 0.7 2.6 2.662 1.9 1.4 3.6 5.6 6.6
利用 complete() 函数可以观察m个插补数据集中的任意一个。格式为：complete(imp, action=#)，其中 # 指定m个完整数据集中的一个来展示，比如：
# 展示了多重插补过程中创建的第三个完整数据集。
dataset3 <- complete(imp, action=3)dataset3BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger645 3 5 342 3 1 31344deleted to save space...]
如果你对缺失值的多重插补法有进一步兴趣，可以进一步参考该书中的相关内容及资源。
18.8 处理缺失值的其他方法方法四
最后，还有两种仍在使用中的缺失值处理方法，但它们已经过时，都应被舍弃，分别是成对删除（pairwise deletion）和简单插补（simple imputation）。
18.8.1 成对删除
对于成对删除，很少使用，观测只是当它含缺失数据的变量涉及某个特定分析时才会被删除。
cor(sleep, use="pairwise.complete.obs")BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger0.30 0.7 0.06 0.3 0.130.7 0.03 0.4 0.15NonD -0.38 -0.37 1.0 0.5 1.0 -0.38 -0.6 -0.32 -0.5 -0.48Dream -0.11 -0.11 0.5 1.0 0.7 -0.30 -0.5 -0.45 -0.5 -0.581.0 -0.41 -0.6 -0.40 -0.6 -0.59Span 0.30 0.51 -0.4 -0.3 -0.4 1.00 0.6 -0.10 0.4 0.06Gest 0.65 0.75 -0.6 -0.5 -0.6 0.61 1.0 0.20 0.6 0.38Pred 0.06 0.03 -0.3 -0.4 -0.4 -0.10 0.2 1.00 0.6 0.92-0.5 0.5 -0.6 0.36 0.6 0.62 1.0 0.79Danger 0.13 0.15 -0.5 -0.6 -0.6 0.06 0.4 0.92 0.8 1.00
此例中，任何两个变量的相关系数都只利用了仅这两变量的可用观测（忽略其他变量）。比如BodyWgt 和 BrainWgt 基于62种（所有变量下的动物数）动物的数据，而 BodyWgt 和 NonD 基于42种动物的数据， Dream 和 NonDream 则基于46种动物的数据。虽然成对删除似乎利用了所有可用数据，但实际上每次计算都只用了不同的数据子集。这将会导致一些扭曲的、难以解释的结果，所以我建议不要使用该方法。
18.8.2 简单（非随机）插补
简单插补，即用某个值（如均值、中位数或众数）来替换变量中的缺失值。若使用均值替换，NonD 中的缺失值可用8.67来替换（两个值分别是Dream 和 NonD 的均值）。注意这些替换是非随机的，这意味着不会引入随机误差（与多重插补不同）。但是对于非MCAR的数据会产生有偏的结果。若缺失数据的数目非常大，那么简单插补很可能会低估标准差、曲解变量间的相关性，并会生成不正确的统计检验的p值。
18.9 小结
在本章中，我们学习了一些鉴别缺失值和探究缺失值模式的方法。学习了产生缺失值的机制，以及分析它们对后续可能产生的影响。同时回顾了三种流行的缺失值处理方法：推理法、行删除法和多重插补。下一章，我们将探究高级作图方法，使用 ggplot2 包创建交互式多元图形。
后台回复“R语言实战”即可获取二维码加入R语言实战学习讨论群。