[Statistics][Stata]Multiple Imputation 多重代入法について

Missing value 欠損値について

欠損値を考えるときに、様々な法則を考慮しないといけない。アンケート調査で、たまたまある項目を答え忘れてしまったとき、これはmissing completely at random (MCAR)という。欠損値がランダムに生じていない場合、例えば「あなたは恋人がいますか」という質問は、恋人がいない人にとっては答えたくない質問で、恋人がいない人は有意にこの質問に答えていないかもしれない。このような欠損値はmissing not at random (MNAR)である。この場合、欠損しているかどうかが欠損値そのものに依存していることをもって、not at randomというのである。missing at random (MAR)は、欠損しているかどうかが欠損値そのものに依存していないものである。女性の多くが体重を自己申告しておらず、体重が多いか少ないかには依存していない場合、この欠損はMARだと言える。

この分類はRubin (1976)による。

欠損値に対する対応

最もわかりやすいのが欠損値がある症例をすべて除外して解析する完全ケース分析 (complete-case analysis)である。その他欠損値をダミー変数（\(9\)や\(99\)などとすることが多い）として補完 (imputation)する場合や、平均値あるいは回帰式で補完する場合がある。ただし、こういったimputationはデータの歪みを生じるためあまり推奨されていない。

Multiple Imputation 多重代入法とは

上の弱点を補う方法として推奨されているのがMultiple Imputation 多重代入法である。イメージとしては，オリジナルのデータセットから欠損値を予測したデータを複数作成し、個々の補完データから関連を推定、推定結果の統合という流れになる。

Multiple Imputationを用いる場合としては、MCARでも欠損値が多い場合や、MARの場合となるが、臨床的にMARなのかMNARなのかを正確に識別するのは困難であり、一般的なツールとして多重代入法が用いられることが多い。完全ケース分析でも多重代入法でも同じ結論になるということをもって、論文の内的妥当性を保証する、という使い方がなされることが多い。

Multiple Imputationの実例

わかりにくいので、実際にデータセットを用いて多重代入法を行う。ここではStataで行うが、SPSSでもRでもSASでも可能である。SPSSではMissing Valuesを購入する。

Stata付属の心筋梗塞データセットを用いる。

use http://www.stata-press.com/data/r16/mheart0.dta

データは以下のようになる。

BMIに欠損値が含まれているので、.になっている箇所がある。

補完する変数が一つの場合

上のデータセットが使える。まずはMIすることを宣言。データセットの形式としては、以下の\(4\)つがある。

わかりにくいのですべてやってみて説明するが、まずは流れを掴む。

mi set wide/mlong/flong/flongset
mi register imputed imputer_vars
mi register passive passive_vars
mi register regular regular_vars

imputedは欠損を補完する変数で、この場合はBMIになる。passiveは欠損があり、imputeされる変数からの変形で補完される。regularは欠損の有無に関わらず、補完が行われない変数になる。

mi register imputed bmi

登録が済んだら、いよいよ実際にimputationを行う。

mi impute chained (regress) reg_var (logit) logit_var (mlogit) mlogit_var (ologit) ologit_var = regular_var, add(n) rseed(m)

regressには連続変数、logitには二値変数、mlogitにはカテゴリー変数、ologitには順序変数を投入する。regular varは上で宣言した補完が行われない変数すべてである。addの後の整数\(n\)で何セットデータセットを生成するか、rseedの後の整数\(m\)で結果を再現するためseedを指定する。この\(m\)は適当な数字で良い。下で実際に行ってみよう。

mi set wide

mi set wide

これを行うと、_mi_missという変数が生成される。この時点ではすべて\(0\)である。

mi register imputed bmi

このコードを実行することで、bmiの中で欠損しているものの_mi_missが\(1\)になる。

mi impute chained (regress) bmi = attack smokes age female hsgrad marstatus alcohol hightar, add(5) rseed(2021)

これでMIは終了である。データセットを見ると、wideの場合補完されたデータが横につながる様子がわかる。

mi set flong

ここも上と同じように進める。

mi set flong
mi register imputed bmi
mi impute chained (regress) bmi = attack smokes age female hsgrad marstatus alcohol hightar, add(5) rseed(2021)

今度は

mi set flong

を行った時点で_mi_miss, _mi_m, _mi_idという\(3\)つの変数が生成される。_mi_missというのが欠損値の有無で\(0\)か\(1\)をとる変数で、_mi_idは各個人のID、_mi_mはadd(n)でデータセットを\(n\)回生成したとき、その何番目にあたるデータセットなのかという番号である。実際にregisterやimputeを行う前は_mi_missと_mi_mは\(0\)である。

mi register imputed bmi

を行うと、欠損値がある場合_mi_missが\(1\)になる。

mi impute chained (regress) bmi = attack smokes age female hsgrad marstatus alcohol hightar, add(5) rseed(2021)

を行うと、縦データが生成される。_mi_mが\(0\)のときは欠損値があるままで、それ以外の場合は欠損値が補完されていることが分かる。

mi set mlong

ここは結果だけ示す。

mi set mlong
mi register imputed bmi
mi impute chained (regress) bmi = attack smokes age female hsgrad marstatus alcohol hightar, add(5) rseed(2021)

これはやっていることはflongとほとんど同じなのだが、flongでは欠損値のないデータもずらずら縦に並んでいくのに対して、mlongでは欠損値のあるものだけを縦に並べる。

つまり、_mi_mが\(0\)のときは最初のデータで、_mi_mが\(0\)でないときは、欠損値のあるデータだけ縦に並べていくということになる（_mi_idをみると、欠損値のある症例だけ掲載されていることが分かる）。flongsepはあまり使わないので割愛。

MIした結果で推測

このデータセットを用いてそのままregeressやlogitをしてはいけない。MIをしたデータセットで推測を行う場合は、

mi estimate:

の後に通常の回帰分析やロジスティック回帰を行う。生存時間分析では

mi stset:

である。新しい変数を作る場合は

mi xeq:

の後にgenerateやreplaceを行う。また、

mi xeq n:regress ..., noisily

などとすることで、n番目のinputationデータの解析を行うことができる。noisilyを指定することで、強制的に出力できるので、エラーが出ている場合など役に立つ。

解析を行ってみると、

mi estimate: logit attack smokes age bmi female

などとなる。mi estimateを付けずに行う解析と比べてみると良い。

ロジスティック回帰で、オッズ比を表示したいときは、

mi est, or: logit outcome var1 var2...

とする。

補完する変数が2つ以上の場合

mi set mlong
mi register imputed var1 var2
mi imputed chained (logit) var1 (regress) var2 = regular var, add(n) rseed(m) force augment

で良い。forceをオプションで用いると、imputeされる値が欠落している場合でもimputeを続行する。augumentオプションが指定された場合、完全予測はカテゴリーデータの代入の過程でロジスティック、順序ロジスティック、あるいは多項ロジスティック代入法を用いて行われる。つまり、二値変数ですべてが\(0\)か\(1\)になるという完全予測を防いでくれる。

エラーについて

“Omitted terms vary error following mi estimate:”のエラーが出るときの多くが、変数の登録忘れである。つまり、

mi est: logit var1 var2 var3

などとしたときにvar3の登録をしていないとエラーとなる。imputeでもpassiveでもregularでもいいので、解析に用いる変数は登録しなくてはいけない。

または、complete解析で、説明変数のreferenceとなるグループでのアウトカム発生件数が\(0\)でも解析できない。どういうことかというと、そもそもMIする前の解析でロジスティック回帰分析が成立しないならば、いくらMIをしても無意味だということである。例えば、心疾患がアウトカム、肥満が説明変数だとする。次のようなテーブルの場合、そもそもロジスティック回帰分析ができない。

この状態でいくらMIを行っても、そもそも発生件数が\(0\)をreferenceとしてロジスティック回帰分析が成立せず、また発生件数が\(0\)件のものをいくらimputationしても\(0\)件のままなので、

mi est: logit...

としても解析は不可能である。このような場合、カテゴリー数を減らすか、あるいは

mi imputed mvn var1 var2 ...  = regular_var1 regular_var2 ..., add(n)

とうようにmvnを指定する。MIした結果を用いてそのままロジスティック回帰分析も行うことが可能であるが、結果の解釈には注意が必要である。

また、少ないNに対して多数の変数を指定すると、ロジスティック回帰の場合perfect predictionになってエラーがでる。その場合は、説明変数のi.を削除して連続変数として解析するか、変数自体のカテゴリーの仕方を変更すると良い。

縦データに対してMultiple Imputationする場合

実はこれが本題。パネルデータなどで縦データの解析を行うことも多いが、縦データの場合欠損値があるタイムポイントに対して、他の時点のデータも重要なデータになりうるので、これを利用しないともったいない。つまり、個人内での属性を考慮しなくてはならない。下のUCLAのweb siteがすべてなのだが、わかりやすく日本語にしてみる。

 How can I perform multiple imputation on longitudinal data using ICE? | Stata FAQ

stats.oarc.ucla.edu

まずはデータの読み込みから。

use "https://stats.idre.ucla.edu/stat/stata/faq/mi_longi.dta", replace

データ構造は以下の様になっている。

どれだけ欠損値があるかは、コマンドnmissingを用いると便利。インストールしていない場合はsscでインストールする。

ssc install nmissing
nmissing

まずはMIしますよと宣言。ここではmlongで良い。

mi set mlong

misstableコマンドで欠損値についてのTableを作成することもできる（これは便利）。

mi misstable summarize female private ses read math

また、

mi misstable patterns var

で欠損値のパターンもTable形式で出力される。

データのreshape

mi misstableで欠損値がどのようなものか大まかな把握が出来たら、縦データを横データに変換する。

mi reshape wide read math, i(id) j(time)

このコマンドを行うと、_mi_idや_mi_miss, _mi_mが生成される。重要なのは、time変数が消失し、reshapeで指定したreadやmathに\(0, 1, 2\)とつくことである。

後は通常のMIと同様に、imputeする変数を登録する。

mi register imputed private ses read1 read2 read3 math1 math2 math3

登録が済んだら実際のimputation。

mi impute chained (logit) private (ologit) ses (regress) read1 read2 read3 math1 math2 math3 = female, add(5) rseed(2021)

Imputationできたら横データから再度縦データに戻す。

mi reshape long read math, i(id) j(time)

これからは通常のMIと同じ。

mi estimate: xtreg read math time female, i(id)

関連リンク

時間の扱い
Difference in Difference 差の差分法
Tableをコマンドひとつで作成する方法「table1_mc」について