（数据科学学习手札58）在R中处理有缺失值数据的高级方法

一、简介

　　在实际工作中，遇到数据中带有缺失值是非常常见的现象，简单粗暴的做法如直接删除包含缺失值的记录、删除缺失值比例过大的变量、用0填充缺失值等，但这些做法会很大程度上影响原始数据的分布或者浪费来之不易的数据信息，因此怎样妥当地处理缺失值是一个持续活跃的领域，贡献出众多巧妙的方法，在不浪费信息和不破坏原始数据分布上试图寻得一个平衡点，在R中用于处理缺失值的包有很多，本文将对最为广泛被使用的mice和VIM包中常用的功能进行介绍，以展现处理缺失值时的主要路径；

二、相关函数介绍

2.1 缺失值预览部分

　　在进行缺失值处理之前，首先应该对手头数据进行一个基础的预览：

　　1、matrixplot

　　效果类似matplotlib中的matshow，VIM包中的matrixplot将数据框或矩阵中数据的缺失及数值分布以色彩的形式展现出来，下面是利用matrixplot对R中自带的airquality数据集进行可视化的效果：

rm(list=ls())
library(mice)
library(VIM)

#以airquality数据为例
data(airquality)
data <- airquality
#利用矩阵热力图查看缺失情况，红色代表缺失
matrixplot(data)

　　红色部分即代表数据缺失值所在位置，通过这个方法，可以在最开始对数据整体的缺失情况有一个初步认识，如通过上图可以一眼看出变量Ozone缺失情况较为严重；

　　2、marginplot与marginmatrix

　　缺失值是否符合完全随机缺失是在对数据进行插补前要着重考虑的事情，VIM中的marginplot包可以同时分析两个变量交互的缺失关系，依然以airquality数据为例：

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

　　如上图所示，通过marginplot传入二维数据框，这里选择airquality中包含缺失值的前两列变量，其中左侧对应变量Solar.R的红色箱线图代表与Ozone缺失值对应的Solar.R未缺失数据的分布情况，蓝色箱线图代表与Ozone未缺失值对应的Solar.R未缺失数据的分布情况，下侧箱线图同理，当同一侧红蓝箱线图较为接近时可认为其对应考察的另一侧变量缺失情况比较贴近完全随机缺失，这种情况下可以放心大胆地进行之后的插补，否则就不能冒然进行插补；

　　与marginplot功能相似，marginmatrix在marginplot只能展现两个变量的基础上推广到多个变量两两之间，效果类似相关性矩阵图：

marginmatrix(data)

　　3、自编函数计算各个变量缺失比例

　　为了计算出每一列变量具体的缺失值比例，可以自编一个简单的函数来实现该功能：

> #查看数据集中每一列的缺失比例
> miss.prop <- function(x){sum(is.na(x))/length(x)}
> apply(data,2,miss.prop)
     Ozone    Solar.R       Wind       Temp      Month        Day 
0.24183007 0.04575163 0.00000000 0.00000000 0.00000000 0.00000000

　　通过自编的函数miss.prop，可以对每个变量中缺失值所占比例有个具体的了解；

2.2 mice函数

　　mice包中最核心的函数是mice()，其主要参数解释如下：

data: 传入待插补的数据框或矩阵，其中缺失值应表示为NA

m: 生成插补矩阵的个数，mice最开始基于gibbs采样从原始数据出发为每个缺失值生成初始值以供之后迭代使用，而m则控制具体要生成的完整初始数据框个数，在整个插补过程最后需要利用这m个矩阵融合出最终的插补结果，若m=1，则唯一的矩阵就是插补的结果；

method: 这个参数控制了传入数据框中每一个变量对应的插补方式，无缺失值的变量对应的为空字符串，带有缺失值的变量默认方法为"pmm"，即均值插补

predictorMatrix: 因为mice中绝大部分方法是用拟合的方式以含缺失值变量之外的其他变量为自变量，缺失值为因变量构建回归或分类模型，以达到预测插补的目的，而参数predictorMatrix则用于控制在对每一个含缺失值变量的插补过程中作为自变量的有哪些其他变量，具体用法下文示例中会详细说明

maxit: 整数，用于控制每个数据框迭代插补的迭代次数，默认为5

seed: 随机数种子，控制随机数水平

在对缺失值插补过程中，非常重要的是为不同的变量选择对应的方法，即method中对应的输入，下表是每种算法对应的参数代号、适用数据类型和算法名称：

`方法代号`	适用数值类型	对应的具体算法名称
`pmm`	any	Predictive mean matching
`midastouch`	any	Weighted predictive mean matching
`sample`	any	Random sample from observed values
`cart`	any	Classification and regression trees
`rf`	any	Random forest imputations
`mean`	numeric	Unconditional mean imputation
`norm`	numeric	Bayesian linear regression
`norm.nob`	numeric	Linear regression ignoring model error
`norm.boot`	numeric	Linear regression using bootstrap
`norm.predict`	numeric	Linear regression, predicted values
`quadratic`	numeric	Imputation of quadratic terms
`ri`	numeric	Random indicator for nonignorable data
`logreg`	binary	Logistic regression
`logreg.boot`	binary	Logistic regression with bootstrap
`polr`	ordered	Proportional odds model
`polyreg`	unordered	Polytomous logistic regression
`lda`	unordered	Linear discriminant analysis
`2l.norm`	numeric	Level-1 normal heteroscedastic
`2l.lmer`	numeric	Level-1 normal homoscedastic, lmer
`2l.pan`	numeric	Level-1 normal homoscedastic, pan
`2l.bin`	binary	Level-1 logistic, glmer
`2lonly.mean`	numeric	Level-2 class mean
`2lonly.norm`	numeric	Level-2 class normal

　　在面对数据集具体情况时，对插补方法进行微调是很必要的步骤，在上面铺垫了这么多之后，下面在具体示例上进行演示，并引入其他的辅助函数；

2.3 利用mice进行缺失值插补——以airquality数据为例

　　因为前面对缺失值预览部分已经利用airquality进行演示，这里就不再赘述，直接进入正式插补部分，首先，我们将data传入mice函数，注意这里设置maxit为0以取得未开始迭代的初始模型参数：

#初始化插补模型,这里最大迭代次数选0是为了取得未开始插补的朴素模型参数
init <- mice(data, maxit = 0)

　　下面我们来看看取得的需要进行调整的重要参数的初始情况：

> #取得对于每一个变量的初始插补方法
> methods <- init$method
> methods
  Ozone Solar.R    Wind    Temp   Month     Day 
  "pmm"   "pmm"      ""      ""      ""      ""

　　可以看到对应缺失变量Ozone和Solar.R的插补拟合方法为pmm，下面我们把它们改成CART决策树回归：

#将变量Ozone的插补方法从pmm修改为norm
methods[c("Ozone")] <- 'cart'
#将变量Solar.R的插补方法从pmm修改为norm
methods[c("Solar.R")] <- 'cart'

　　接着我们来查看predictorMatrix参数：

> #取得对每一个变量进行拟合用到的变量矩阵，0代表不用到，1代表用到
> predM <- init$predictorMatrix
> predM
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

　　这里我们认为变量Month和Day是日期，与缺失变量无相关关系，因此将其在矩阵中对应位置修改为0使它们不参与拟合过程：

#调整参与拟合的变量
#这里认为日期对与其他变量无相关关系，因此令变量Month与变量Day不参与对其他变量的拟合插补过程
predM[, c("Month")] <- 0
predM[c("Month"),] <- 0
predM[, c("Day")] <- 0
predM[c("Day"),] <- 0

　　这样我们就完成了两个重要参数的初始化，下面我们进行正式的拟合插补：

#利用修改后的参数组合来进行拟合插补
imputed <- mice(data, method = methods, predictorMatrix = predM)

　　随着程序运行完，我们需要的结果便呼之欲出，但在取得最终插补结果前，为了严谨起见，需要对模型的统计学意义进行分析，下面以Ozone为例：

　　1、查看模型中Ozone对应的拟合公式：

> #查看Ozone主导的拟合公式
> imputed$formulas['Ozone']
$Ozone
Ozone ~ 0 + Solar.R + Wind + Temp
<environment: 0x000000002424d410>

　　可以看到，Ozone对应的公式与前面predictorMatrix参数中经过修改的保持一致；

　　2、基于上述公式为合成出的m=5个数据框分别进行拟合：

> #把上面的公式填入下面的lm()内
> fit <- with(imputed, lm(Ozone ~ Solar.R + Wind + Temp))
> 
> #查看fit中对应每一个插补数据框的回归显著性结果
> fit
call :
with.mids(data = imputed, expr = lm(Ozone ~ Solar.R + Wind + 
    Temp))

call1 :
mice(data = data, method = methods, predictorMatrix = predM)

nmis :
  Ozone Solar.R    Wind    Temp   Month     Day 
     37       7       0       0       0       0 

analyses :
[[1]]

Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)

Coefficients:
(Intercept)      Solar.R         Wind         Temp  
  -53.36269      0.05791     -3.37688      1.51550  


[[2]]

Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)

Coefficients:
(Intercept)      Solar.R         Wind         Temp  
  -55.70273      0.06189     -3.25456      1.50816  


[[3]]

Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)

Coefficients:
(Intercept)      Solar.R         Wind         Temp  
  -55.45601      0.05659     -2.90911      1.47244  


[[4]]

Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)

Coefficients:
(Intercept)      Solar.R         Wind         Temp  
  -70.00005      0.07008     -2.56784      1.59856  


[[5]]

Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)

Coefficients:
(Intercept)      Solar.R         Wind         Temp  
  -72.08381      0.05252     -2.71741      1.67060

　　3、将上述5个拟合模型融合并查看融合后的显著性水平

> #计算fit中模型的显著性
> pooled <- pool(fit)
> summary(pooled)
                estimate   std.error statistic       df      p.value
(Intercept) -61.32105668 21.84851057 -2.806647 53.60143 6.964478e-03
Solar.R       0.05979673  0.02144304  2.788632 90.71532 6.449107e-03
Wind         -2.96516062  0.67509037 -4.392243 29.06277 1.361726e-04
Temp          1.55305117  0.23531131  6.599985 78.14524 4.468489e-09

　　可以看到从截距项，到每一个变量的p值都远远小于0.05，至少在0.05显著性水平下每个参数都具有统计学意义；

　　4、对5个合成出的数据框在缺失值位置进行融合，这里需要用到新的函数complete，其主要有下面三个参数：

data: 前面mice函数输出的结果

action: 当只希望从合成出的m个数据框中取得某个单独的数据框时，可以设置action参数，如action=3便代表取得m个数据框中的第3个

mild: 逻辑型变量，当为TRUE时，会输出包含全部m个合成数据框的列表

　　获悉上列参数意义后，若只想抽取某个数据框如第3个：

result <- complete(imputed, action = 3)
matrixplot(result)

　　可以看到，取回第3个数据框，每个缺失值都已被补全，若希望得到5个合成数据框的融合结果，则需要自编函数：

#取得所有合成数据框组成的列表
complete(imputed, mild = T)
all.imputed <- complete(imputed, action = 'all')
#对得到的m个合成数据框缺失值部分求平均
avgDF <- function(data,all.imputed){
  baseDF <- all.imputed[[1]][is.na(data)]
  for(child in 2:length(all.imputed)){
    baseDF <- baseDF + all.imputed[[child]][is.na(data)]
  }
  data[is.na(data)] <- baseDF/length(all.imputed)
  
  return(data)
}

#得到最终平均数据框
result <- avgDF(data,all.imputed)
matrixplot(result)