Multiple imputation
山本 良平
2020/5/4
Introduction
データの欠落は非常に大きなトピック
ここでは、多重代入法を取り巻く計算上の問題について、表面にふれる
詳細については、Biostats学部のRod LittleのStatistical Analysis With Missing Dataコースをがおすすめ
http://blog.caspio.com/tech-tips/what-you-need-to-know-about-data-harvesting-and-how-to-prevent-it/
Types of missing data
欠落データは、以下の3つのカテゴリーに分類される
MCAR
完全にランダムに欠落しているデータは、どのオブザベーションが値を欠落しているかについての体系的なものを何も持たない
欠落度と観察された共変量または観察されていない共変量との間には関係がない
MAR
Missing At RandomはMCARよりも弱い
欠落はやはりランダムであるが完全に観察された変数に起因している
例えば、社会経済的地位の低い人々は、給与情報を提供する意欲が低い（しかし、我々は彼らのSESの地位を知っている）
重要なのは、Missingが観測されていない値によるものではないということ
MCARはMARを意味するが、その逆ではない
MNAR
データが無作為に欠落していない場合、欠落度は欠落データの値に依存する
打ち切りデータはこのカテゴリに該当する
例えば、体重が重い人は、体重を報告する可能性が低くなる
別の例として、ある反応を測定するデバイスは0.5以上の値しか測定できない
それ以下の値は欠落する
Determing which type of missingness
残念ながら、あなたのデータがどのカテゴリーに該当するかを決定する統計的な方法はない
MCARとMARは、データとその収集メカニズムに関するあなたの知識に基づいて行われる必要がある
データにはMNARを知らせることができるものは何もない
もし欠損が未観察値によるものであるならば，データには我々を助けるものは何もない
しかし，MNARはデフォルトの仮定であるべきである
洞察を得るために，共変量に基づいて欠落度を予測するロジスティック・モデルを適合させることができまる
もし何らかの関係を見つけたら、それはMARがMCARよりも現実的であることを暗示している。 しかし，それは代わりにMNARがないことを保証するものではない
同様に、もし関係が見つからない場合、それは次のような理由かもしれない
The missingness is MCAR.
The missingness is MNAR.
The missingness is really MAR, but you haven’t modeled it appropriately.
Methods to handle missingness
データがMCAR
- 完全ケース分析が有効
- 多重入力またはその他の入力方法が有効
データがMAR
- いくつかのComplete-case分析は、MCARよりも弱い仮定の下で有効
例: 線形回帰は、予測変数に条件付きで、ミッシングネスが応答から独立している場合、偏りがない
- 多重入力は有効である（偏りはあるが，その偏りは無視できる）
データがMNAR
- Missingnessを明示的にモデル化しなければならない; responseとMissingnessを共同でモデル化する
- いくつかの特定のケース（例：生存分析）では，MNARデータ（例：打ち切りデータ）は適切に扱われる
- 一般的には，このような状況を避けるために，可能な限りMARを仮定する
Multiple imputation
多重計算の一般的な手順は簡単
- m個の完全なケース・データ・セットを生成するために、いくつかの分布からランダムに引かれた値で欠落値をインプリメントする
- m個のデータセットのそれぞれで同じ分析を実行する
- 結果を何らかの方法でプールする
最初のステップについては、バリエーションがある
理解するのに最も単純なものの1つは、他のすべての変数（ミッシングネスを持つ変数と持たない変数）を使用して、ミッシングネスを持つ各変数を予測する一連の回帰モデルを含む このモデルに基づいて、我々は各欠落値の期待値の分布を得て、そこから標本化することができる
ステップ3で使用される正確な方法は、ステップ2からの分析に依存する
一般的に、プールされた推定統計量は、単にm個のデータセットのそれぞれからの統計量の平均であり、標準誤差は、複雑さが生じるかもしれない
mice in R
Rでは、Mouse(Multiple Imputation with Chained Equations)を使って、多重代入とその後の解析をおこなう
mice、with関数、pool関数をそれぞれ介して上記の3つのステップを実行する
mice
mice関数は連鎖方程式を介して入力を実行する
詳細にはあまり詳しく説明しないが、連鎖方程式は、ギブス・サンプラー（MCMCアプローチ）のバリエーションで、欠落値の推定値と変数の分布のパラメータの推定値（両方とも他の変数に条件付き）の間で反復処理を行う
連鎖方程式は、YtとYt-1の条件付き独立性のため、ほとんどのMCMCと比較して収束が速い
これは、推定値の “効率性”のクローズド・フォーム推定からの80年代と90年代のオリジナルの仕事から来ている
25%の欠損データと5つの入力で、効率性は95%以上
しかし、いくつかの最新の研究（http://sci2s.ugr.es/keel/pdf/specific/articulo/graham_olchowski_07.pdf） では、実際にはより多くの入力を行うことでより多くの利益を得ることができ、オリジナルの論文が発表されたときに問題となっていた計算要求は明らかにもはや懸念事項ではないという議論をしている
これらの著者は、m = 100のようなより大きな数を使うことを提案している
実際には、現在の標準は（私が知っているすべての分野で）m = 5のままである
懸念すべきもう1つの主要な議論、method関数は、各インピュートされた値がどのように計算されるかを制御する
いずれも連鎖方程式アプローチを採用していますが、詳細は大きく異なります。
あまり詳細に飛び込まずに、かなりの数があり、オリジナルのmiceの論文で十分に文書化されています。
我々の目的のために、デフォルトのpmm関数(予測平均マッチング)を使うことにします。
これは、ある意味でセミパラメトリックであるという利点があります（入力モデルが誤って指定された場合でも非線形関係が保持されます）。 Original mice paper：https://www.jstatsoft.org/article/view/v045i03
library(mice)## Warning: package 'mice' was built under R version 3.5.3##
## Attaching package: 'mice'## The following objects are masked from 'package:base':
##
## cbind, rbinddata(nhanes)
nhanes.imp <- mice(nhanes, m = 5, method = 'pmm')##
## 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
## 2 4 bmi hyp chl
## 2 5 bmi hyp chl
## 3 1 bmi hyp chl
## 3 2 bmi hyp chl
## 3 3 bmi hyp chl
## 3 4 bmi hyp chl
## 3 5 bmi hyp chl
## 4 1 bmi hyp chl
## 4 2 bmi hyp chl
## 4 3 bmi hyp chl
## 4 4 bmi hyp chl
## 4 5 bmi hyp chl
## 5 1 bmi hyp chl
## 5 2 bmi hyp chl
## 5 3 bmi hyp chl
## 5 4 bmi hyp chl
## 5 5 bmi hyp chl以下のような方法で補完を調べることができます。
# Show imputation 3
head(complete(nhanes.imp, 3))## age bmi hyp chl
## 1 1 25.5 1 238
## 2 2 22.7 1 187
## 3 1 28.7 1 187
## 4 3 24.9 2 206
## 5 1 20.4 1 113
## 6 3 24.9 2 184# Show just the imputed values. Columns are imputations, rows are observations
nhanes.imp$imp## $age
## [1] 1 2 3 4 5
## <0 rows> (or 0-length row.names)
##
## $bmi
## 1 2 3 4 5
## 1 25.5 35.3 25.5 30.1 27.2
## 3 30.1 30.1 28.7 30.1 30.1
## 4 21.7 27.2 24.9 27.4 20.4
## 6 22.5 21.7 24.9 27.4 25.5
## 10 27.4 26.3 22.7 27.4 22.5
## 11 33.2 28.7 22.0 30.1 27.2
## 12 28.7 27.4 27.4 27.5 22.5
## 16 22.7 20.4 27.2 33.2 30.1
## 21 22.5 33.2 22.5 33.2 28.7
##
## $hyp
## 1 2 3 4 5
## 1 1 1 1 1 1
## 4 1 2 2 2 2
## 6 2 2 2 1 1
## 10 2 1 1 1 1
## 11 1 1 1 1 1
## 12 2 1 2 1 2
## 16 1 1 1 1 1
## 21 1 1 1 1 1
##
## $chl
## 1 2 3 4 5
## 1 113 187 238 184 131
## 4 186 284 206 206 204
## 10 218 131 199 199 131
## 11 229 131 187 184 187
## 12 184 199 204 206 187
## 15 187 187 229 218 187
## 16 187 113 187 218 187
## 20 206 218 186 218 184
## 21 113 199 131 204 187
## 24 184 199 206 204 184# Extract the "tall" matrix which stacks the imputations
nhanes.comp <- complete(nhanes.imp, "long", include = TRUE)
# Now, the .imp variable identifies which imputation each belongs to.
# 0 (included because of `include = TRUE` above) is the un-imputed data.
table(nhanes.comp$.imp)##
## 0 1 2 3 4 5
## 25 25 25 25 25 25# Let's visualize the distribution of imputed and observed values.
# `cci` returns logical whether its input is complete at each observation.
nhanes.comp$bmi.NA <- cci(nhanes$bmi)
# Note that the above uses the recylcing properties of matrixes/data.frame:
# The `cci` call returns length 25; but because values are recylced to the total
# number of rows in nhanes.comp, it replicates 6 times.
head(nhanes.comp[, c("bmi", "bmi.NA")])## bmi bmi.NA
## 1 NA FALSE
## 2 22.7 TRUE
## 3 NA FALSE
## 4 NA FALSE
## 5 20.4 TRUE
## 6 NA FALSElibrary(ggplot2)
ggplot(nhanes.comp,
aes(x = .imp, y = bmi, color = bmi.NA)) +
geom_jitter(show.legend = FALSE,
width = .1)## Warning: Removed 9 rows containing missing values (geom_point).
デフォルトのメソッド(pmm)を使用しているので、極端な値を持たないことが保証されています。
他のメソッドではそうではないので、これはいいチェックになります。
With
withはRの一般的な関数で、attachを使ってコードを実行し、その後detachを実行します。
しかし、midsオブジェクト(mice()の出力)をwithに渡すと、入力を適切に処理する特別なmiceバージョンになります。
withは有効なRコマンドを第二引数として受け取り、それぞれのインputedデータセットに適用します。いくつかの例があります。
with(nhanes.imp, mean(bmi))## call :
## with.mids(data = nhanes.imp, expr = mean(bmi))
##
## call1 :
## mice(data = nhanes, m = 5, method = "pmm")
##
## nmis :
## age bmi hyp chl
## 0 9 8 10
##
## analyses :
## [[1]]
## [1] 26.372
##
## [[2]]
## [1] 27.012
##
## [[3]]
## [1] 26.032
##
## [[4]]
## [1] 27.656
##
## [[5]]
## [1] 26.368with(nhanes.imp, t.test(bmi ~ hyp))## call :
## with.mids(data = nhanes.imp, expr = t.test(bmi ~ hyp))
##
## call1 :
## mice(data = nhanes, m = 5, method = "pmm")
##
## nmis :
## age bmi hyp chl
## 0 9 8 10
##
## analyses :
## [[1]]
##
## Welch Two Sample t-test
##
## data: bmi by hyp
## t = -0.24546, df = 22.014, p-value = 0.8084
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.179585 2.506569
## sample estimates:
## mean in group 1 mean in group 2
## 26.27778 26.61429
##
##
## [[2]]
##
## Welch Two Sample t-test
##
## data: bmi by hyp
## t = 0.80703, df = 17.945, p-value = 0.4302
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.924624 4.324624
## sample estimates:
## mean in group 1 mean in group 2
## 27.3 26.1
##
##
## [[3]]
##
## Welch Two Sample t-test
##
## data: bmi by hyp
## t = -0.46461, df = 22.833, p-value = 0.6466
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.895945 1.834041
## sample estimates:
## mean in group 1 mean in group 2
## 25.88333 26.41429
##
##
## [[4]]
##
## Welch Two Sample t-test
##
## data: bmi by hyp
## t = 0.71665, df = 22.384, p-value = 0.481
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.503335 3.093335
## sample estimates:
## mean in group 1 mean in group 2
## 27.815 27.020
##
##
## [[5]]
##
## Welch Two Sample t-test
##
## data: bmi by hyp
## t = 1.0566, df = 11.24, p-value = 0.3129
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.797987 5.134829
## sample estimates:
## mean in group 1 mean in group 2
## 26.76842 25.10000with(nhanes.imp, lm(bmi ~ age + chl))## call :
## with.mids(data = nhanes.imp, expr = lm(bmi ~ age + chl))
##
## call1 :
## mice(data = nhanes, m = 5, method = "pmm")
##
## nmis :
## age bmi hyp chl
## 0 9 8 10
##
## analyses :
## [[1]]
##
## Call:
## lm(formula = bmi ~ age + chl)
##
## Coefficients:
## (Intercept) age chl
## 19.67355 -3.32501 0.06707
##
##
## [[2]]
##
## Call:
## lm(formula = bmi ~ age + chl)
##
## Coefficients:
## (Intercept) age chl
## 21.24074 -4.35818 0.07121
##
##
## [[3]]
##
## Call:
## lm(formula = bmi ~ age + chl)
##
## Coefficients:
## (Intercept) age chl
## 19.64958 -2.02677 0.05135
##
##
## [[4]]
##
## Call:
## lm(formula = bmi ~ age + chl)
##
## Coefficients:
## (Intercept) age chl
## 20.33830 -3.29865 0.06679
##
##
## [[5]]
##
## Call:
## lm(formula = bmi ~ age + chl)
##
## Coefficients:
## (Intercept) age chl
## 21.32152 -3.39089 0.05935別の方法として、nhanes.comp内の.impを反復処理して、それぞれの解析を行うこともできます。
一般的にはこの方法はお勧めできませんが、必要に応じてこれを行うことができます。
pool
最後に、lm, aov, lme (lme4から)のように、coef()とvcov()の両方のメソッドを持つオブジェクトを生成する任意のメソッド（見るためにmethod(coef)またはmethod(vcov)を実行します）に対して、プール関数は、ルービンのプールのルールに基づいて、個々の係数推定値を1つの “プールされた”推定値に結合することができます。
前述のように、ポイント推定値は厳密平均ですが、標準誤差はより複雑です。
また、線形モデルのプールされたR2推定値を得るための pool.r.squaredも有用です。
mod1 <- with(nhanes.imp, lm(bmi ~ age + chl))
pool(mod1)## Class: mipo m = 5
## term m estimate ubar b t dfcom
## 1 (Intercept) 5 20.44473864 8.550936059 6.601721e-01 9.343142564 22
## 2 age 5 -3.27990083 0.646572851 6.869321e-01 1.470891384 22
## 3 chl 5 0.06315382 0.000274842 6.183414e-05 0.000349043 22
## df riv lambda fmi
## 1 17.926991 0.09264559 0.08479015 0.1722571
## 2 5.237949 1.27490434 0.56042108 0.6671415
## 3 13.505505 0.26997678 0.21258403 0.3079966summary(pool(mod1),digits=3)## term estimate std.error statistic df p.value
## 1 (Intercept) 20.44473864 3.05665545 6.688598 17.926991 2.898417e-06
## 2 age -3.27990083 1.21280311 -2.704397 5.237949 4.056917e-02
## 3 chl 0.06315382 0.01868269 3.380338 13.505505 4.691733e-03pool.r.squared(mod1)## est lo 95 hi 95 fmi
## R^2 0.4990683 0.1295403 0.7773827 NaNConclusion
私がこのトピックを含めたのは、クックブックの手法に問題がある場合、複数のインピュートされたデータセット（この例ではnhanes.comp）を扱う際に混乱してしまうアナリストによく遭遇するからです。このような計算方法の基本的な理解があれば、そのような障害を乗り越えることができるのではないかと期待しています。
最後に、Stataの多重インputation手法は、より広い範囲のプールされた推定値をサポートしていますので、多重インputationが分析の要であるならば、チェックしてみることを検討してみてはいかがでしょうか。
Link: http://dept.stat.lsa.umich.edu/~jerrick/courses/stat701/notes/mi.html