多重插补(MI)是一种基于重复模拟的处理缺失值的方法。在面对复杂的缺失值问题时,MI是最常用的方法,它将从一个包含缺失值的数据集中生成一组完整的数据集(通常是3到10个)。在每个模拟数据集中,缺失数据将用蒙特卡洛方法来填补。此时,标准的统计方法便可应用到每个模拟的数据集上,通过组合输出结果给出估计的结果,以及引入 缺失值时的置信区间。R语言中可利用Amelia.mice和mi包来执行这些操作。我们重点学习mice包中提供的方法。
通过mice包应用多重插补的步骤如下:
mice()函数首先从一个包含缺失数据的数据框开始,然后返回一个或多个(默认为5个)完整数据集的对象。每个完整数据集都是通过对原始数据框中的缺失数据进行插补而生成的。由于插补有随机的成分,因此每个完整的数据集都略有不同。然后,with()函数可依次对每个完整数据集应用统计模型(如线性模型或广义线性模型)。最后,pool()函数将这些单独的分析结果整合为一组结果。最终模型的标准误差和p值都将准确地反映出由于缺失值和多重插补而产生的不确定性。
基于mice包的分析通常符合以下几个分析过程:
library(mice)
imp <- mice(mydata,m)
fit <- with(imp, analysis)
pooled <- pool(fit)
summary(pooled)
mydata是一个包含缺失值的矩阵或数据框
imp是一个包含m个插补数据集的列表对象,同时还含有完成插补过程的信息。默认地,m=5;
analysis是一个表达式对象,用来设定应用与m个插补数据集的统计分析方法。方法包括做线性回归模型的lm()函数、做广义线性模型的glm()函数、做广义可加模型的gam(),以及做负二型模型的nbrm()函数。表达式在函数的括号中,~左边是响应变量,右边是预测变量(用+符号分隔开)。
fit是一个包含m个单独统计分析结果的列表对象。
pooled是一个包含这m个统计分析平均结果的列表对象。
现将多重插补法应用到sleep数据集上,我们设定随机种子为1234,应用代码如下:
> library(mice)
> data(sleep, package = "VIM")
> imp <- mice(sleep, seed = 1234)
iter imp variable
1 1 NonD Dream Sleep Span Gest
1 2 NonD Dream Sleep Span Gest
1 3 NonD Dream Sleep Span Gest
1 4 NonD Dream Sleep Span Gest
1 5 NonD Dream Sleep Span Gest
2 1 NonD Dream Sleep Span Gest
2 2 NonD Dream Sleep Span Gest
2 3 NonD Dream Sleep Span Gest
2 4 NonD Dream Sleep Span Gest
2 5 NonD Dream Sleep Span Gest
3 1 NonD Dream Sleep Span Gest
3 2 NonD Dream Sleep Span Gest
3 3 NonD Dream Sleep Span Gest
3 4 NonD Dream Sleep Span Gest
3 5 NonD Dream Sleep Span Gest
4 1 NonD Dream Sleep Span Gest
4 2 NonD Dream Sleep Span Gest
4 3 NonD Dream Sleep Span Gest
4 4 NonD Dream Sleep Span Gest
4 5 NonD Dream Sleep Span Gest
5 1 NonD Dream Sleep Span Gest
5 2 NonD Dream Sleep Span Gest
5 3 NonD Dream Sleep Span Gest
5 4 NonD Dream Sleep Span Gest
5 5 NonD Dream Sleep Span Gest
> fit <- with(imp,lm(Dream~Span + Gest))
> pooled <- pool(fit)
> summary(pooled)
est se t
(Intercept) 2.546199168 0.254689696 9.997260
Span -0.004548904 0.012039106 -0.377844
Gest -0.003916211 0.001468788 -2.666287
df Pr(>|t|) lo 95
(Intercept) 52.12563 1.021405e-13 2.035156222
Span 51.94538 7.070861e-01 -0.028707741
Gest 55.55683 1.002562e-02 -0.006859066
hi 95 nmis fmi lambda
(Intercept) 3.0572421151 NA 0.08710301 0.05273554
Span 0.0196099340 4 0.08860195 0.05417409
Gest -0.0009733567 4 0.05442170 0.02098354
此处,我们可以看到Span的回归系数不显著(p=0.08),Gest的系数在p<0.01的水平下很显著。若将这些结果与利用完整数据分析法所得的结果进行对比,我们会发现背离的结论相同。当控制寿命不变时,妊娠期与做梦时长有一个显著的、负相关的关系。完整数据分析法基于42个有完整数据的动物,而此处的分析法基于整个数据集中全部62个动物的数据。另外,fmi栏也展示了确实信息(即由于引入缺失数据而引起的编译所占整体不确定性的比例)。
我们可以通过检查分析过程所创建的对象来获取更多的差不信息,例如,来看看imp对象的汇总信息:
> imp
Multiply imputed data set
Call:
mice(data = sleep, seed = 1234)
Number of multiple imputations: 5
Missing cells per column:
BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger
0 0 14 12 4 4 4 0 0 0
Imputation methods:
BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger
"" "" "pmm" "pmm" "pmm" "pmm" "pmm" "" "" ""
VisitSequence:
NonD Dream Sleep Span Gest
3 4 5 6 7
PredictorMatrix:
BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger
BodyWgt 0 0 0 0 0 0 0 0 0 0
BrainWgt 0 0 0 0 0 0 0 0 0 0
NonD 1 1 0 1 1 1 1 1 1 1
Dream 1 1 1 0 1 1 1 1 1 1
Sleep 1 1 1 1 0 1 1 1 1 1
Span 1 1 1 1 1 0 1 1 1 1
Gest 1 1 1 1 1 1 0 1 1 1
Pred 0 0 0 0 0 0 0 0 0 0
Exp 0 0 0 0 0 0 0 0 0 0
Danger 0 0 0 0 0 0 0 0 0 0
Random generator seed value: 1234
从输出结果可以看到,五个数据集同时被创建,预测均值(pmm)匹配法被用来处理每个含缺失数据的变量。BodyWgt、BrainWgt、Pred、Exp和Danger没有进行插补(“”),因为它们并没有缺失数据。VisitSequence从左到右展示了插补的变量,从NonD开始,至Gest结束。最后,预测变量矩阵(PredictorMatrix)展示了进行插补过程的含有缺失数据的变量,它们利用了数据集中其他变量的信息。(在矩阵中,行代表插补变量,列代表为插补提供信息的变量,1和0分别表示使用和未使用)。
通过提取imp对象的子成分,可以观测到实际的插补值。如:
> imp$imp$Dream
1 2 3 4 5
1 1.0 0.5 0.5 0.5 0.3
3 2.6 2.1 1.5 1.8 1.3
4 3.4 3.1 3.4 1.2 3.4
14 0.3 0.5 0.5 0.3 1.2
24 1.8 1.3 3.6 0.9 5.6
26 2.3 3.1 2.0 2.6 2.1
30 1.2 0.3 3.4 2.6 2.3
31 3.4 0.5 0.6 1.0 0.5
47 0.5 1.5 1.5 2.2 3.4
53 0.3 0.5 0.5 0.5 0.6
55 0.5 0.9 2.6 2.7 2.4
62 1.0 2.1 0.5 3.9 3.6
展示了在Dream变量上有缺失值的12个动物的5次插补值。检查该矩阵可以帮助我们判断插补值是否合理。若睡眠时长出现了负值,差不将会停止。
利用complete()函数可以观察m个插补数据集中的任意一个,格式为:
complete(imp,action=#)
其中,#指定m个完整数据集中的一个来展示,例如:
> dataset3 <- complete(imp,action = 3)
> dataset3
BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger
1 6654.000 5712.00 2.1 0.5 3.3 38.6 645.0 3 5 3
2 1.000 6.60 6.3 2.0 8.3 4.5 42.0 3 1 3
3 3.385 44.50 11.0 1.5 12.5 14.0 60.0 1 1 1
4 0.920 5.70 13.2 3.4 16.5 2.0 25.0 5 2 3
5 2547.000 4603.00 2.1 1.8 3.9 69.0 624.0 3 5 4
6 10.550 179.50 9.1 0.7 9.8 27.0 180.0 4 4 4
7 0.023 0.30 15.8 3.9 19.7 19.0 35.0 1 1 1
8 160.000 169.00 5.2 1.0 6.2 30.4 392.0 4 5 4
9 3.300 25.60 10.9 3.6 14.5 28.0 63.0 1 2 1
10 52.160 440.00 8.3 1.4 9.7 50.0 230.0 1 1 1
11 0.425 6.40 11.0 1.5 12.5 7.0 112.0 5 4 4
12 465.000 423.00 3.2 0.7 3.9 30.0 281.0 5 5 5
13 0.550 2.40 7.6 2.7 10.3 18.0 46.0 2 1 2
14 187.100 419.00 3.3 0.5 3.1 40.0 365.0 5 5 5
15 0.075 1.20 6.3 2.1 8.4 3.5 42.0 1 1 1
16 3.000 25.00 8.6 0.0 8.6 50.0 28.0 2 2 2
17 0.785 3.50 6.6 4.1 10.7 6.0 42.0 2 2 2
18 0.200 5.00 9.5 1.2 10.7 10.4 120.0 2 2 2
19 1.410 17.50 4.8 1.3 6.1 34.0 28.0 1 2 1
20 60.000 81.00 12.0 6.1 18.1 7.0 21.0 1 1 1
21 529.000 680.00 11.0 0.3 10.8 28.0 400.0 5 5 5
22 27.660 115.00 3.3 0.5 3.8 20.0 148.0 5 5 5
23 0.120 1.00 11.0 3.4 14.4 3.9 16.0 3 1 2
24 207.000 406.00 8.2 3.6 12.0 39.3 252.0 1 4 1
25 85.000 325.00 4.7 1.5 6.2 41.0 310.0 1 3 1
26 36.330 119.50 11.0 2.0 13.0 16.2 63.0 1 1 1
27 0.101 4.00 10.4 3.4 13.8 9.0 28.0 5 1 3
28 1.040 5.50 7.4 0.8 8.2 7.6 68.0 5 3 4
29 521.000 655.00 2.1 0.8 2.9 46.0 336.0 5 5 5
30 100.000 157.00 7.4 3.4 10.8 22.4 100.0 1 1 1
31 35.000 56.00 3.3 0.6 3.8 16.3 33.0 3 5 4
32 0.005 0.14 7.7 1.4 9.1 2.6 21.5 5 2 4
33 0.010 0.25 17.9 2.0 19.9 24.0 50.0 1 1 1
34 62.000 1320.00 6.1 1.9 8.0 100.0 267.0 1 1 1
35 0.122 3.00 8.2 2.4 10.6 12.7 30.0 2 1 1
36 1.350 8.10 8.4 2.8 11.2 4.5 45.0 3 1 3
37 0.023 0.40 11.9 1.3 13.2 3.2 19.0 4 1 3
38 0.048 0.33 10.8 2.0 12.8 2.0 30.0 4 1 3
39 1.700 6.30 13.8 5.6 19.4 5.0 12.0 2 1 1
40 3.500 10.80 14.3 3.1 17.4 6.5 120.0 2 1 1
41 250.000 490.00 7.7 1.0 8.4 23.6 440.0 5 5 5
42 0.480 15.50 15.2 1.8 17.0 12.0 140.0 2 2 2
43 10.000 115.00 10.0 0.9 10.9 20.2 170.0 4 4 4
44 1.620 11.40 11.9 1.8 13.7 13.0 17.0 2 1 2
45 192.000 180.00 6.5 1.9 8.4 27.0 115.0 4 4 4
46 2.500 12.10 7.5 0.9 8.4 18.0 31.0 5 5 5
47 4.288 39.20 11.0 1.5 12.5 13.7 63.0 2 2 2
48 0.280 1.90 10.6 2.6 13.2 4.7 21.0 3 1 3
49 4.235 50.40 7.4 2.4 9.8 9.8 52.0 1 1 1
50 6.800 179.00 8.4 1.2 9.6 29.0 164.0 2 3 2
51 0.750 12.30 5.7 0.9 6.6 7.0 225.0 2 2 2
52 3.600 21.00 4.9 0.5 5.4 6.0 225.0 3 2 3
53 14.830 98.20 2.1 0.5 2.6 17.0 150.0 5 5 5
54 55.500 175.00 3.2 0.6 3.8 20.0 151.0 5 5 5
55 1.400 12.50 8.4 2.6 11.0 12.7 90.0 2 2 2
56 0.060 1.00 8.1 2.2 10.3 3.5 100.0 3 1 2
57 0.900 2.60 11.0 2.3 13.3 4.5 60.0 2 1 2
58 2.000 12.30 4.9 0.5 5.4 7.5 200.0 3 1 3
59 0.104 2.50 13.2 2.6 15.8 2.3 46.0 3 2 2
60 4.190 58.00 9.7 0.6 10.3 24.0 210.0 4 3 4
61 3.500 3.90 12.8 6.6 19.4 3.0 14.0 2 1 1
62 4.050 17.00 17.9 0.5 19.4 13.0 38.0 3 1 1
上述代码展示了多重插补过程中创建的第三个完整数据集。
处理缺失值的其他方法
R语言还支持其他一些处理缺失值的方法,例如:
最后,还有两种仍在使用但已过时的方法:成对删除和简单插补
处理含缺失值的数据集时,成对删除常作为行删除的备选方法使用。对于成对删除,观测只是当它含缺失数据的变量涉及某个特定分析时才会被删除,考虑如下代码:
> cor(sleep,use = "pairwise.complete.obs")
BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger
BodyWgt 1.00000000 0.93416384 -0.3759462 -0.1093833 -0.3071859 0.30245056 0.6511022 0.05949472 0.3382737 0.13358123
BrainWgt 0.93416384 1.00000000 -0.3692177 -0.1051388 -0.3581020 0.50925268 0.7472425 0.03385548 0.3678004 0.14587888
NonD -0.37594625 -0.36921766 1.0000000 0.5142539 0.9627147 -0.38443179 -0.5947028 -0.31818462 -0.5437566 -0.48385220
Dream -0.10938331 -0.10513879 0.5142539 1.0000000 0.7270870 -0.29574535 -0.4508987 -0.44747050 -0.5372245 -0.57933653
Sleep -0.30718591 -0.35810203 0.9627147 0.7270870 1.0000000 -0.41020239 -0.6313262 -0.39583497 -0.6422845 -0.58774241
Span 0.30245056 0.50925268 -0.3844318 -0.2957453 -0.4102024 1.00000000 0.6148488 -0.10254416 0.3603522 0.06177846
Gest 0.65110218 0.74724248 -0.5947028 -0.4508987 -0.6313262 0.61484879 1.0000000 0.20050426 0.6382790 0.37861701
Pred 0.05949472 0.03385548 -0.3181846 -0.4474705 -0.3958350 -0.10254416 0.2005043 1.00000000 0.6182460 0.91604245
Exp 0.33827367 0.36780037 -0.5437566 -0.5372245 -0.6422845 0.36035221 0.6382790 0.61824597 1.0000000 0.78720311
Danger 0.13358123 0.14587888 -0.4838522 -0.5793365 -0.5877424 0.06177846 0.3786170 0.91604245 0.7872031 1.00000000
在这个例子中,任何两个变量的相关系数都只利用了仅这两变量的可用观测(忽略其他变量)。比如BodyWgt和BrainWgt基于62个动物的数据,而BodyWgt和NonD基于42个动物的数据,Dream和NonDream基于46个动物的数据。
虽然成对删除似乎利用了所有可用数据,但实际上每次计算都只用了不同的数据子集。这将会导致一些扭曲的、难以解释的结果,因此不建议使用这个方法。
简单插补,即用某个值(如均值、中位数或众数)来替换变量中的缺失值。若使用均值替换,Dream变量中的缺失值可用1.97(Dream的均值)来替换,NonD中的缺失值可用8.67(NonD的均值)来替换。这些替换时非随机的,这意味着不会引入随机误差。
简单插补的一个优点是,解决“缺失值问题”时不会减少分析过程中可用的样本量。虽然简单插补用法简单,但对于非MCAR的数据会产生有偏的结果。若缺失数据的数目非常大,那么简单插补很可能会低估标准差‘曲解变量间的相关性,并会产生不正确的统计检验p值。与成对删除一样,在解决缺失数据问题时也应尽量避免使用该方法。