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

1、简介

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

 

2、相关函数介绍

2.1  缺失值预览部分

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

　　一、matrixplot

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

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

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

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

　　二、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)

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

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

> #查看数据集中每一列的缺失比例
> 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为例：

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

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

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

　　二、基于上述公式为合成出的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  

　　三、将上述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显著性水平下每一个参数都具备统计学意义；

　　四、对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)

 

 

　　以上就是本文的所有内容，若有错误之处望斧正。

 参考资料：

https://stefvanbuuren.name/Winnipeg/Lectures/Winnipeg.pdf

https://www.rdocumentation.org/packages/mice/versions/3.5.0/topics/mice