数据缺失的机制分为三类

假设数据缺失是完全随机缺失,但是如果数据缺失太多,也可能造成很大的问题。一般来说,数据的缺失量
小于数据总量的5%是可以接受的,那么如果某个特征或样本的缺失数据量超过5%,那么可以考虑是否需要留
下该特征,我们可以写一个简单的函数来检查某一特征或者某一样本的缺失量是否超过5%

使用 mice 中的 airquality 来进行介绍

library(mice)
## Loading required package: lattice
## 
## Attaching package: 'mice'
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
data <- airquality

查看数据的缺失情况

pMiss <- function(x) {sum(is.na(x))/length(x)*100}
apply(data,2,pMiss)
##     Ozone   Solar.R      Wind      Temp     Month       Day 
## 24.183007  4.575163  0.000000  0.000000  0.000000  0.000000
apply(data,1,pMiss)
##   [1]  0.00000  0.00000  0.00000  0.00000 33.33333 16.66667  0.00000
##   [8]  0.00000  0.00000 16.66667 16.66667  0.00000  0.00000  0.00000
##  [15]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
##  [22]  0.00000  0.00000  0.00000 16.66667 16.66667 33.33333  0.00000
##  [29]  0.00000  0.00000  0.00000 16.66667 16.66667 16.66667 16.66667
##  [36] 16.66667 16.66667  0.00000 16.66667  0.00000  0.00000 16.66667
##  [43] 16.66667  0.00000 16.66667 16.66667  0.00000  0.00000  0.00000
##  [50]  0.00000  0.00000 16.66667 16.66667 16.66667 16.66667 16.66667
##  [57] 16.66667 16.66667 16.66667 16.66667 16.66667  0.00000  0.00000
##  [64]  0.00000 16.66667  0.00000  0.00000  0.00000  0.00000  0.00000
##  [71]  0.00000 16.66667  0.00000  0.00000 16.66667  0.00000  0.00000
##  [78]  0.00000  0.00000  0.00000  0.00000  0.00000 16.66667 16.66667
##  [85]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
##  [92]  0.00000  0.00000  0.00000  0.00000 16.66667 16.66667 16.66667
##  [99]  0.00000  0.00000  0.00000 16.66667 16.66667  0.00000  0.00000
## [106]  0.00000 16.66667  0.00000  0.00000  0.00000  0.00000  0.00000
## [113]  0.00000  0.00000 16.66667  0.00000  0.00000  0.00000 16.66667
## [120]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [127]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [134]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [141]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [148]  0.00000  0.00000 16.66667  0.00000  0.00000  0.00000

summary the dataset

summary(data)
##      Ozone           Solar.R           Wind             Temp      
##  Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00  
##  1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:72.00  
##  Median : 31.50   Median :205.0   Median : 9.700   Median :79.00  
##  Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88  
##  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                                       
##      Month            Day      
##  Min.   :5.000   Min.   : 1.0  
##  1st Qu.:6.000   1st Qu.: 8.0  
##  Median :7.000   Median :16.0  
##  Mean   :6.993   Mean   :15.8  
##  3rd Qu.:8.000   3rd Qu.:23.0  
##  Max.   :9.000   Max.   :31.0  
## 

Ozone变量缺失最多,有37个缺失

利用mice包查看数据缺失的模式

md.pattern(data)

##     Wind Temp Month Day Solar.R Ozone   
## 111    1    1     1   1       1     1  0
## 35     1    1     1   1       1     0  1
## 5      1    1     1   1       0     1  1
## 2      1    1     1   1       0     0  2
##        0    0     0   0       7    37 44

从上面的输出结果可以得到,有111个样本是完整的,35个样本缺失Ozone特征值,Solar.R变量有7个缺失值
Ozone有37个缺失值

使用VIM包也可以对数据模式进行可视化

## Loading required package: colorspace
## Loading required package: grid
## Loading required package: data.table
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
## VIM is ready to use. 
##  Since version 4.0.0 the GUI is in its own package VIMGUI.
## 
##           Please use the package to use the new (and old) GUI.
## Suggestions and bug-reports can be submitted at: https://github.com/alexkowa/VIM/issues
## 
## Attaching package: 'VIM'
## The following object is masked from 'package:datasets':
## 
##     sleep

## 
##  Variables sorted by number of missings: 
##  Variable      Count
##     Ozone 0.24183007
##   Solar.R 0.04575163
##      Wind 0.00000000
##      Temp 0.00000000
##     Month 0.00000000
##       Day 0.00000000

从图片中可以看出70%的样本没有缺失信息,23%的样本缺失Ozone特征

另一种有用的可视化方法

这个函数一次只能显示两个变量的缺失情况,左边红色boxplot代表了缺失Ozone数据项,但是Solar.R数据项
没有缺失的样本分布,而左边的蓝色boxplot则显示了剩余的数据点的分布。类似的,下面的两个图形也是这
么进行描述的
如果我们的数据服从完全随机分布的话,那么红色 boxplot 和蓝色 boxplot应该很相似

Imputing the missing data

tempData <- mice(data,
                 m = 5,  # 产生查补数据集的个数,默认是5个
                 maxit = 50,
                 method = "pmm", # 计算插值的方式,这里是用的预测均值
                 seed = 2019,
                 printFlag = F)
summary(tempData)
## Class: mids
## Number of multiple imputations:  5 
## Imputation methods:
##   Ozone Solar.R    Wind    Temp   Month     Day 
##   "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

查看Ozone的插补情况

tempData$imp$Ozone
##       1   2   3   4   5
## 5    19   8  19  19  14
## 10   30  21  30  23  23
## 25   19  14   6   8  18
## 26    4  37  37  19  13
## 27   24  18  32  23  24
## 32   52  45  45  35  32
## 33   18  16  23  16   7
## 34    1  37  32  18  32
## 35   49  63  20  89  16
## 36   89  35  59 110  49
## 37   30  16  12  23  30
## 39   66 135  61  77  85
## 42   78  78  85  76  78
## 43   76  76  78  73  85
## 45   45  28  18  32  23
## 46   32  46  28  63  52
## 52   39  28  71  35  20
## 53   61  40 168  78  59
## 54   64  29  78  23  47
## 55  108  47  64  48  59
## 56   35  31  46  28  16
## 57   37  29  63  20  29
## 58   65  12  23  44  12
## 59   16   7  39  65  20
## 60   30  21   9   9  23
## 61   82  35  85  73  59
## 65   39  31  23  45  36
## 72   47  28  59  63  45
## 75   89  39  52  37  20
## 83   71  23  35  59  16
## 84   59  63  44  35  29
## 102 135 135  73  77  96
## 103  28  59  36  20  44
## 107  16  11  10  24   7
## 115  20  41  11  30  14
## 119  97  50  79  78 110
## 150  30  11   9  14   9

结果展示了每一个观测样本Ozone变量的插补值

查看每一个变量的插值方法

tempData$method
##   Ozone Solar.R    Wind    Temp   Month     Day 
##   "pmm"   "pmm"      ""      ""      ""      ""

现在,我们就可以对原始数据进行最后操作

completedata <- complete(tempData,1)    # 1 表示使用5个插值数据集的第一个

查看插补前后数据集的分布情况

首先,画一个 Ozone和其他变量间的散点图

xyplot(tempData,Ozone~Wind+Temp+Solar.R,pch = 18,cex = 1)

上图中,品红色的点(imputed)与蓝色的点(observed)几乎匹配,预示者插补方法很可靠

另一个有用的图是密度曲线图

densityplot(tempData)

密度图中,每一个插补数据集用品红色表示,实际观测值用蓝色表示

还有一个图用于展示预测的插值在变量中的分布情况

stripplot(tempData,pch = 20, cex = 1.2)

Pooling

假设在后续的步骤中,你想对数据进行一个线性回归分析,那么应该选取那个插补数据集呢?MICE 包 提供了很方便的方法,先用每一个插补的数据集进行模型建立,最后将几个模型综合评估,得到最优的
结果

with-pool

modelfit <- with(tempData,lm(Temp~Ozone+Solar.R+Wind))
summary(pool(modelfit))
##                estimate   std.error statistic        df      p.value
## (Intercept) 72.39012292 2.934926474 24.665055  93.36193 0.000000e+00
## Ozone        0.16548481 0.028554047  5.795494  24.27635 5.406154e-06
## Solar.R      0.01054416 0.008871597  1.188530  17.18176 2.507933e-01
## Wind        -0.34603937 0.208260501 -1.661570 115.30217 9.931570e-02