答應(yīng)大家今天要寫(xiě)一寫(xiě)關(guān)于多重插補(bǔ)的文章。查閱了很多資料纵寝，發(fā)現(xiàn)了一篇非常棒的英文文獻(xiàn)绸狐，翻譯給大家谣妻。
附上原文鏈接：
文章名稱(chēng)：mice: Multivariate Imputation by Chained Equations in R
文章鏈接：https://www.jstatsoft.org/article/view/v045i03
多重插補(bǔ)原理圖
加載mice包 library('mice')
今天所用的數(shù)據(jù)包含了四個(gè)變量：age(年齡組)兽叮，bmi(體重指數(shù))芬骄，hyp(血壓狀況) ，anchl(膽固醇水平)鹦聪。數(shù)據(jù)作為數(shù)據(jù)框存儲(chǔ)账阻，缺失值表示為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
檢查缺失值
缺失值的數(shù)量可以按以下方式計(jì)算和可視化：
> 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。缺失值總數(shù)等于(7×3)+(1×2)+(3×1)+(1×1)=27。大多數(shù)缺失值(10)發(fā)生在chl中蒲牧。
另一種研究模式的方法是計(jì)算每個(gè)模式對(duì)所有變量的觀(guān)測(cè)數(shù)撇贺。一對(duì)變量完全可以有四個(gè)缺失模式：兩個(gè)變量都被觀(guān)察到(模式rr)，第一個(gè)變量被觀(guān)察到造成，第二個(gè)變量缺失(模式rm)显熏，第一個(gè)變量丟失雄嚣，第二個(gè)變量被觀(guān)察到(模式mr)晒屎，并且兩者都缺失(模式mm)。我們可以使用md.pairs()函數(shù)來(lái)計(jì)算所有變量對(duì)在每個(gè)模式中的頻率缓升，如:
> 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
因此鼓鲁，對(duì)(bmi，chl)有13對(duì)完全觀(guān)察到港谊，3對(duì)bmi未觀(guān)察到骇吭，2對(duì)bmi缺失但有hyp觀(guān)察，7對(duì)bmi和hyp均缺失歧寺。請(qǐng)注意燥狰，這些數(shù)字加起來(lái)等于總樣本大小。
2
多重插補(bǔ)
多重插補(bǔ)可以通過(guò)調(diào)用mice()來(lái)完成斜筐，如下所示：
> 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
···
其中龙致，多重插補(bǔ)的數(shù)據(jù)集存儲(chǔ)在mids類(lèi)的對(duì)象imp中。檢查結(jié)果是什么樣子
> 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
插補(bǔ)值是根據(jù)默認(rèn)的方法生成的顷链，這就是對(duì)于數(shù)值數(shù)據(jù)目代，預(yù)測(cè)平均匹配(pmm)。默認(rèn)的多個(gè)插補(bǔ)值數(shù)量等于m=5嗤练。
3
診斷檢驗(yàn)
多重插補(bǔ)的一個(gè)重要步驟是評(píng)估插補(bǔ)值是否可信榛了。插補(bǔ)值應(yīng)該是如果沒(méi)有丟失就可以得到的值。插補(bǔ)值應(yīng)該接近數(shù)據(jù)煞抬。顯然不可能的數(shù)據(jù)值(例如負(fù)數(shù)霜大、懷孕的父親)不應(yīng)出現(xiàn)在插補(bǔ)值的數(shù)據(jù)中。插補(bǔ)值應(yīng)尊重變量之間的關(guān)系革答，并反映其“真實(shí)”值的適當(dāng)不確定性战坤。對(duì)插補(bǔ)數(shù)據(jù)的診斷檢查提供了一種檢查插補(bǔ)值合理性的方法。bmi的插補(bǔ)值存儲(chǔ)為
> 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
每一行的第一個(gè)數(shù)對(duì)應(yīng)于bmi中缺失的數(shù)據(jù)蝗碎。這些列就是等待的插補(bǔ)值湖笨。完整的數(shù)據(jù)集組合了觀(guān)察值和插補(bǔ)值。完整的數(shù)據(jù)集可以被觀(guān)察
> 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()函數(shù)從imp對(duì)象中提取五個(gè)輸入數(shù)據(jù)集蹦骑，作為一個(gè)包含125條記錄的長(zhǎng)矩陣(行堆疊)慈省。nhances中缺失的條目現(xiàn)在已經(jīng)由第一列(五個(gè))插補(bǔ)值來(lái)填充。第二個(gè)完整的數(shù)據(jù)集可以通過(guò)complete(imp，2)獲得边败。對(duì)于觀(guān)察到的數(shù)據(jù)袱衷，它與第一個(gè)完整的數(shù)據(jù)集相同，但插補(bǔ)數(shù)據(jù)是不同的笑窜。
檢查原始數(shù)據(jù)和插補(bǔ)后數(shù)據(jù)的分布通常是有用的致燥。這樣做的一種方法是在mice中使用函數(shù)stripplot()
> stripplot(imp,pch=19,cex=1.2,alpha=.3)
圖顯示這四個(gè)變量的分布情況。觀(guān)察值到藍(lán)色點(diǎn)排截，插補(bǔ)值是紅色點(diǎn)嫌蚤。age面板只包含藍(lán)色點(diǎn)，因?yàn)?strong>age已經(jīng)是完整數(shù)據(jù)断傲。此外脱吱，請(qǐng)注意，紅色點(diǎn)相當(dāng)好地跟隨藍(lán)色點(diǎn)认罩，包括分布中的間隙箱蝠，例如，對(duì)于chl垦垂。
下圖中每個(gè)插補(bǔ)數(shù)據(jù)集的chl和bmi散點(diǎn)圖是由xyplot()函數(shù)繪制：
> xyplot(imp,bmi ~ chl / .imp,pch=20,cex=1.2)
插補(bǔ)值是用紅色繪制的宦搬。藍(lán)色的點(diǎn)在不同的面板上是相同的，但是紅色的點(diǎn)是不同的劫拗。紅點(diǎn)與藍(lán)色數(shù)據(jù)的形狀大致相同间校，這表明，如果沒(méi)有遺漏杨幼，它們可能是合理的插補(bǔ)值撇簿。紅點(diǎn)之間的差異代表了我們對(duì)真實(shí)(但未知)值的不確定性。
在完全隨機(jī)缺失下差购，觀(guān)測(cè)數(shù)據(jù)和插補(bǔ)數(shù)據(jù)的單變量分布是一致的四瘫。在隨機(jī)缺失下，它們可以是不同的欲逃，無(wú)論是在位置上還是在傳播上找蜜，但它們的多元分布假設(shè)是相同的。有許多其他方法可以查看已完成的數(shù)據(jù)稳析。
4
分析插補(bǔ)數(shù)據(jù)
假設(shè)完整數(shù)據(jù)分析是age和bmi的線(xiàn)性回歸洗做。為此，我們可以使用函數(shù)with.mids()彰居，這是一個(gè)包裝器函數(shù)诚纸，它將完整的數(shù)據(jù)模型應(yīng)用于每個(gè)已插補(bǔ)的數(shù)據(jù)集：
> fit <- with(imp,lm(chl ~ age + bmi))
fit對(duì)象具有mira類(lèi)，包含五個(gè)完整的數(shù)據(jù)分析結(jié)果陈惰。這些可以按以下方式匯總
> 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
經(jīng)多次插補(bǔ)后畦徘，我們發(fā)現(xiàn)bmi有顯著性影響。列fmi包含Rubin(1987)中定義的缺失信息的部分，列是lambda可歸因于缺失數(shù)據(jù)(λ=(b/b+m)/t)的總方差的比例井辆。合并的結(jié)果會(huì)受到模擬誤差的影響关筒，因此取決于mice()函數(shù)的seed參數(shù)。為了最小化模擬誤差杯缺，我們可以使用更多的插補(bǔ)蒸播，例如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
我們發(fā)現(xiàn)萍肆，age和bmi實(shí)際上都有顯著的影響袍榆。這是很好的結(jié)果。
提取插補(bǔ)結(jié)果匾鸥，以第二個(gè)插補(bǔ)數(shù)據(jù)為例吧
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便是完整的數(shù)據(jù)啦蜡塌。
通過(guò)翻譯這篇文章碉纳，在下終于見(jiàn)識(shí)到了大神長(zhǎng)什么樣勿负。真實(shí)啊，這文章寫(xiě)的太帥了劳曹！
希望各位看得愉快奴愉。
下期再見(jiàn)。
你可能還想看
等你很久啦锭硼，長(zhǎng)按加入古同社區(qū)
