答应大家今天要写一写关于多重插补的文章。查阅了很多资料，发现了一篇非常棒的英文文献，翻译给大家。
附上原文链接：
文章名称：mice: Multivariate Imputation by Chained Equations in R
文章链接：https://www.jstatsoft.org/article/view/v045i03
多重插补原理图
加载mice包 library('mice')
今天所用的数据包含了四个变量：age(年龄组)，bmi(体重指数)，hyp(血压状况) ，anchl(胆固醇水平)。数据作为数据框存储，缺失值表示为NA。
> nhanes
age bmi hyp chl
1 1 NA NA NA
2 2 22.7 1 187
3 1 NA 1 187
4 3 NA NA NA
5 1 20.4 1 113
6 3 NA NA 184
7 1 22.5 1 118
···
1
检查缺失值
缺失值的数量可以按以下方式计算和可视化：
> md.pattern(data)
age hyp bmi chl
13 1 1 1 1 0
3 1 1 1 0 1
1 1 1 0 1 1
1 1 0 0 1 2
7 1 0 0 0 3
0 8 9 10 27
13行是完整的（25行中）。有一行只缺少bmi，有七行只知道age。缺失值总数等于(7×3)+(1×2)+(3×1)+(1×1)=27。大多数缺失值(10)发生在chl中。
另一种研究模式的方法是计算每个模式对所有变量的观测数。一对变量完全可以有四个缺失模式：两个变量都被观察到(模式rr)，第一个变量被观察到，第二个变量缺失(模式rm)，第一个变量丢失，第二个变量被观察到(模式mr)，并且两者都缺失(模式mm)。我们可以使用md.pairs()函数来计算所有变量对在每个模式中的频率，如:
> md.pairs(nhanes)
$`rr`
age bmi hyp chl
age 25 16 17 15
bmi 16 16 16 13
hyp 17 16 17 14
chl 15 13 14 15
$rm
age bmi hyp chl
age 0 9 8 10
bmi 0 0 0 3
hyp 0 1 0 3
chl 0 2 1 0
$mr
age bmi hyp chl
age 0 0 0 0
bmi 9 0 1 2
hyp 8 0 0 1
chl 10 3 3 0
$mm
age bmi hyp chl
age 0 0 0 0
bmi 0 9 8 7
hyp 0 8 8 7
chl 0 7 7 10
因此，对(bmi，chl)有13对完全观察到，3对bmi未观察到，2对bmi缺失但有hyp观察，7对bmi和hyp均缺失。请注意，这些数字加起来等于总样本大小。
2
多重插补
多重插补可以通过调用mice()来完成，如下所示：
> imp <- mice(nhanes,seed = 23109)
iter imp variable
1 1 bmi hyp chl
1 2 bmi hyp chl
1 3 bmi hyp chl
1 4 bmi hyp chl
1 5 bmi hyp chl
2 1 bmi hyp chl
2 2 bmi hyp chl
2 3 bmi hyp chl
···
其中，多重插补的数据集存储在mids类的对象imp中。检查结果是什么样子
> print(imp)
Class: mids
Number of multiple imputations: 5
Imputation methods:
age bmi hyp chl
"" "pmm" "pmm" "pmm"
PredictorMatrix:
age bmi hyp chl
age 0 1 1 1
bmi 1 0 1 1
hyp 1 1 0 1
chl 1 1 1 0
插补值是根据默认的方法生成的，这就是对于数值数据，预测平均匹配(pmm)。默认的多个插补值数量等于m=5。
3
诊断检验
多重插补的一个重要步骤是评估插补值是否可信。插补值应该是如果没有丢失就可以得到的值。插补值应该接近数据。显然不可能的数据值(例如负数、怀孕的父亲)不应出现在插补值的数据中。插补值应尊重变量之间的关系，并反映其“真实”值的适当不确定性。对插补数据的诊断检查提供了一种检查插补值合理性的方法。bmi的插补值存储为
> imp$imp$bmi
1 2 3 4 5
1 28.7 27.2 22.5 27.2 28.7
3 30.1 30.1 30.1 22.0 28.7
4 22.7 27.2 20.4 22.7 20.4
6 24.9 25.5 22.5 21.7 21.7
10 30.1 29.6 27.4 25.5 25.5
11 35.3 26.3 22.0 27.2 22.5
12 27.5 26.3 26.3 24.9 22.5
16 29.6 30.1 27.4 30.1 25.5
21 33.2 30.1 35.3 22.0 22.7
每一行的第一个数对应于bmi中缺失的数据。这些列就是等待的插补值。完整的数据集组合了观察值和插补值。完整的数据集可以被观察
> complete(imp)
age bmi hyp chl
1 1 28.7 1 199
2 2 22.7 1 187
3 1 30.1 1 187
4 3 22.7 2 204
5 1 20.4 1 113
6 3 24.9 2 184
7 1 22.5 1 118
8 1 30.1 1 187
···
complete()函数从imp对象中提取五个输入数据集，作为一个包含125条记录的长矩阵(行堆叠)。nhances中缺失的条目现在已经由第一列(五个)插补值来填充。第二个完整的数据集可以通过complete(imp，2)获得。对于观察到的数据，它与第一个完整的数据集相同，但插补数据是不同的。
检查原始数据和插补后数据的分布通常是有用的。这样做的一种方法是在mice中使用函数stripplot()
> stripplot(imp,pch=19,cex=1.2,alpha=.3)
图显示这四个变量的分布情况。观察值到蓝色点，插补值是红色点。age面板只包含蓝色点，因为age已经是完整数据。此外，请注意，红色点相当好地跟随蓝色点，包括分布中的间隙，例如，对于chl。
下图中每个插补数据集的chl和bmi散点图是由xyplot()函数绘制：
> xyplot(imp,bmi ~ chl / .imp,pch=20,cex=1.2)
插补值是用红色绘制的。蓝色的点在不同的面板上是相同的，但是红色的点是不同的。红点与蓝色数据的形状大致相同，这表明，如果没有遗漏，它们可能是合理的插补值。红点之间的差异代表了我们对真实(但未知)值的不确定性。
在完全随机缺失下，观测数据和插补数据的单变量分布是一致的。在随机缺失下，它们可以是不同的，无论是在位置上还是在传播上，但它们的多元分布假设是相同的。有许多其他方法可以查看已完成的数据。
4
分析插补数据
假设完整数据分析是age和bmi的线性回归。为此，我们可以使用函数with.mids()，这是一个包装器函数，它将完整的数据模型应用于每个已插补的数据集：
> fit <- with(imp,lm(chl ~ age + bmi))
fit对象具有mira类，包含五个完整的数据分析结果。这些可以按以下方式汇总
> print(pool(fit))
Class: mipo m = 5
estimate ubar b t dfcom df riv
(Intercept) 5.958468 3575.717813 1649.43830 5555.043768 22 9.216954 0.5535465
age 29.725360 82.553435 115.72029 221.417779 22 4.331856 1.6821147
bmi 5.138790 3.686724 0.91985 4.790544 22 12.907767 0.2994040
lambda fmi
(Intercept) 0.3563115 0.4616878
age 0.6271599 0.7288640
bmi 0.2304164 0.3271721
经多次插补后，我们发现bmi有显著性影响。列fmi包含Rubin(1987)中定义的缺失信息的部分，列是lambda可归因于缺失数据(λ=(b/b+m)/t)的总方差的比例。合并的结果会受到模拟误差的影响，因此取决于mice()函数的seed参数。为了最小化模拟误差，我们可以使用更多的插补，例如m=50。这样做很容易
> imp50 <- mice(nhanes,m=50,seed = 23109)
> fit <- with(imp50,lm(chl ~ age + bmi))
> print(pool(fit))
Class: mipo m = 50
estimate ubar b t dfcom df riv
(Intercept) -7.276001 3193.70850 1200.0871230 4417.797368 22 14.30395 0.3832813
age 32.075691 78.93751 52.6214971 132.611433 22 11.58145 0.6799547
bmi 5.470019 3.33579 0.9781109 4.333463 22 15.32201 0.2990815
lambda fmi
(Intercept) 0.2770813 0.3606366
age 0.4047458 0.4863912
bmi 0.2302254 0.3142527
我们发现，age和bmi实际上都有显著的影响。这是很好的结果。
提取插补结果，以第二个插补数据为例吧
data2 <- complete(imp50,2)
> data2
age bmi hyp chl
1 1 28.7 1 187
2 2 22.7 1 187
3 1 22.0 1 187
4 3 22.5 1 186
5 1 20.4 1 113
6 3 25.5 2 184
7 1 22.5 1 118
8 1 30.1 1 187
9 2 22.0 1 238
10 2 20.4 1 238
11 1 28.7 1 199
12 2 20.4 2 199
13 3 21.7 1 206
14 2 28.7 2 204
15 1 29.6 1 229
16 1 27.2 1 187
17 3 27.2 2 284
18 2 26.3 2 199
19 1 35.3 1 218
20 3 25.5 2 186
21 1 29.6 2 218
22 1 33.2 1 229
23 1 27.5 1 131
24 3 24.9 1 186
25 2 27.4 1 186
data2便是完整的数据啦。
通过翻译这篇文章，在下终于见识到了大神长什么样。真实啊，这文章写的太帅了！
希望各位看得愉快。
下期再见。
你可能还想看
等你很久啦，长按加入古同社区
