一篇文章写清楚一个问题，关注我，自学python！

哈哈，写了好几篇R了，应该改成关注我自学R了，R和python我都会写，今天给大家写写如何用R做多重插补。

多重插补介绍

多重插补（MI）是一种基于重复模拟的处理缺失值的方法。在面对复杂的缺失值问题时，MI是最常用的方法，它将从一个包含缺失值的数据集中生成一组完整的数据集（通常是3到10个）。在每个模拟数据集中，缺失数据将用蒙特卡洛方法来填补。此时，标准的统计方法便可应用到每个模拟的数据集上，通过组合输出结果给出估计的结果，以及引入缺失值时的置信区间。R语言中可利用Amelia.mice和mi包来执行这些操作。今天写一写mice包中提供的方法。

实例操练

本文使用airquality这个数据集（R中自带的，哈哈）

首先我制造一点缺失数据，代码如下：

data <- airquality
data[4:10,3] <- rep(NA,7)
data[1:5,4] <- NA

然后把数据集中的分类变量去除，并看一看数据分布情况，代码如下：

data <- data[-c(5,6)]
summary(data)

可以看到数据集中有4个连续性变量的分布和每个变量的缺失值个数。

缺失数据分类

缺失数据的分类大概有3中（详细请看我之前的文章）：

• 非随机缺失
• 随机缺失
• 完全随机缺失

基本上，我们最期望自己的缺失数据是完全随机缺失，因为这种缺失最好处理，但是也不能缺多，最多不要超过5%。如果一个样本的特征或者一个变量缺失超过了5%我们就会考虑将这个样本或者特征剔除。所以接下来要做的就是看一看在我们的数据集中样本或者变量的缺失比例，代码如下

pMiss <- function(x){sum(is.na(x))/length(x)*100}
apply(data,2,pMiss)
apply(data,1,pMiss)

我们可以看到Ozone这个变量的缺失比例都快到25%了，这个时候就要考虑丢掉它或者补数据。再看样本，因为我们每一个样本就4个变量，如果一个样本的变量缺失达到2肯定得丢掉。

使用mice查看数据缺失类型

mice包有一个非常好用的查看缺失的方法md.pattern()

library(mice)
md.pattern(data)

从上面的表格和图形，我们很容易看出数据集中有104个完整样本，其中Ozone这个变量缺失最多。

当然啊，上面的方法并不直接告诉我们数据集缺失数据的比例，看缺失比例，代码如下：

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"))

上面的输出就可以很清晰的让我们知道，数据缺失比例是多少，很容易看出，我们的数据集中，有70%是完整数据。

插补缺失

多重插补的代码如下：

tempData <- mice(data,m=5,maxit=50,meth='pmm',seed=500)
summary(tempData)

上面的代码会输出插补的迭代过程和上图，给大家解释一下各个参数的意思：

• m=5 代表插补的数据集个数，默认为5
• meth='pmm' 代表插补方法，PMM为默认方式：预测均值匹配，你也可以输入methods(mice)看一看R提供的别的插补方法。

怎么看插补完成的数据呢？很简单，比如我现在想看Ozone这个变量的完整数据，只需要以下代码：

tempData$imp$Ozone

现在不是已经插补完成了嘛，此时就需要获取完整数据集，代码很简单：

completedData <- complete(tempData,1)

上面的代码就是取我们插补的第一个数据集，我相信，你肯定知道取第二个数据集的代码怎么写了。

对比原始数据和插补数据

我们可以画个图看看插补数据和原始数据的分布情况，首先是点图，代码如下

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

在上面这个图中，我们需要看的是紫红色的点也就是插补的点和蓝色的点也就是原始数据的分布是不是一样的，上面的图大概可以看出来，插补的数据和原始数据的分布还是差不多的，说明插补的还好。

还有另一个图：密度图也可以看看

densityplot(tempData)

在上图中，紫红色的都是插补数据集的分布情况，蓝色的是真实数据的分布情况，但是好像这个图也不能分清楚哪个插补的数据集最好，不过也没关系。最终我们都要合并的。

合并

假设我们要对数据做的分析是线性回归，你可能会问我到底该选哪一个数据集做呢？哈，mice直接可以都做好直接给我输出最终合并的结果。

代码如下：

modelFit1 <- with(tempData,lm(Temp~ Ozone+Solar.R+Wind))
summary(pool(modelFit1))

在上面的代码中，modelFit1中就包含了用插补数据集做出来的模型结果，然后用pool方法可以将将所有的模型结果合并在一起，从输出的结果看，只有Ozone这个变量是显著的。

小结

今天给大家写了多重插补的做法，感谢大家耐心看完。发表这些东西的主要目的就是督促自己，希望大家关注评论指出不足，一起进步。内容我都会写的很细，用到的数据集也会在原文中给出链接，你只要按照文章中的代码自己也可以做出一样的结果，一个目的就是零基础也能懂，因为自己就是什么基础没有从零学Python的，加油。

（站外链接发不了，请关注后私信回复“数据链接”获取本头条号所有使用数据）

往期精彩：

R数据分析：如何用R做验证性因子分析及画图，实例操练

R数据分析：结构方程模型画图以及模型比较，实例操练

python数据分析：数据缺失机制及利用missingno库对缺失可视化