Stata分析機能例題を利用してもっと知ろう
多重代入～Multiple Inputation～

多重代入を行う理由

架空の喫煙と心臓発作の関係を調査した症例対照研究データを例に考えます。

変数attackとsmokesの関係をロジットモデルで確認します。

データセットには154の観測値がありますが、欠損値を含むため推定に使用されたのは132のみであることがわかります。
misstable summarizeで欠損値を確認します。

bmiが欠損している観測値が22個あることがわかります。
上記のロジットモデルではこれらを完全に排除して推定しましたが、多重代入を使えば、これらの情報を分析に取り込むことができます。

多重代入はシミュレーションを使用してデータの欠損に対して柔軟に対応する、統計技術です。
以下のステップで行われます

• 代入 ： 選択した代入モデルを使用して、M個の完全な（欠損値の無い）データセットを作成します。
• 推定 ： 代入したデータセットごと（m=1,2,…,M）に、本来行いたい推定を行います。
• プーリング ： M個の分析から得られたそれぞれの結果を1つの結果に統合します。

多重代入の特徴は、

1. 代入データを使用して既存の推定行う
2. 代入と推定が独立している

の2点です。

欠損データを扱う手法で一般的なものが、完全ケース分析またはリストワイズ削除と単一代入ですが、これらには、多重代入と比較すると、いくつか点で問題があります。
リストワイズ削除では、欠損値を含む観測値を削除するため、それらの情報は失われます。
欠損値が多数あれば、分析結果の有効性は多重代入のものに比べ低くなります。
また、残された観測値が母集団とかけ離れていた場合、バイアスのかかったパラメータが導き出されてしまいます。

先のロジットモデルでは、リストワイズ削除を使用しています。
観測値を削除した標本では、ageは統計的に有意ではありませんが、多重代入によって全ての観測値を使用した分析では、ageの有意性が明らかになります。

リストワイズ削除とは異なり、単一代入は欠損値を削除せず、既知の値として扱います。
これでは、分散が過小評価され、精度と信頼区間が過大になり、結果として有意性検定が楽観的になります。
これに対して多重代入では、複数回代入を行うことで、サンプリングの分散を保証しています。
代入のステップが分析と独立している点において、データ収集・代入ではより正確な代入を行うことが、分析の段階では様々な分析を行うことができます。

上の心臓発作の例では、bmiの欠損値を、線形回帰を使用して代入し、ロジットモデルで推定を行います。
推定では欠損を考慮して調整を行う必要がありません。

多重代入の理論的背景

多重代入はベイズ統計の技術を利用し、既知のサンプルから発生させたランダム化分布による点推定値と信頼区間の統計的な妥当性を保証するためにRubin（1996）の定義を使用します。
多重代入では、代入データを作成する代入モデルと分析モデルの2つを設定しなければなりません。
分析モデルでは、興味のあるパラメータ$ Q $の完全データ推定値 $ \hat{Q} $、および$ \hat{Q} $に関するサンプリングの妥当性$ U $を推定します。
プーリング段階では、完全データで推定された個々の$ (\hat{Q}, U) $を$ (\hat{Q}_{MI}, T) $に統合して反復代入を行います。

代入モデルと分析モデルがどちらもベイジアンである場合、ベイジアンモデルによる欠損データの事後分布の反復代入を行います。
統合された推定値$ \hat{Q}_{MI} $とサンプリング分散$ T=W+B $は、事後平均とQの分散に近似します。
$ W $は代入データセット内の分散（完全データの推定分散$ U $の平均）、$ B $は代入データセット間の分散（反復代入間の推定分散 $ \hat{Q}_{MI} $ ）です。
事後平均と分散が事後分散を正しく表していれば、統合された推測結果は、ベイズ統計の観点からも、頻度統計の観点からも、有効であると判断できます。

反復代入による推測の統計的妥当性は、

1. 適切な代入によって多重代入が行われているか
2. 分析段階の$ (\hat{Q}, U) $に基づく完全データ推論が適切であるか

に依ります。
$ (\hat{Q}, U) $に基づく完全データ推論は$ \hat{Q}\sim N \{Q, Var(\hat{Q})\} $であり、$ U $がサンプルの分布$ Var(\hat{Q}) $の一致推定量である場合に有効です。
多重代入のランダム化妥当性は代入を無限に行うことで保障されますが、実際の代入回数は、$ M $の値を決定しなければなりません。
ここで、再度$ Q $を次のように定義します。

\[ T_{M}^{-1/2}(Q-\hat{Q}_{M})\sim t_{\nu M} \]

$ \hat{Q}_M $は$ Q $の完全データ推定値$ M $の平均、$ T_{M}=W+(1+1/M)B $、また$ t_{\nu M} $は代入回数と欠損している情報の割合に依存する、自由度$ \nu_{M} $のStudentの$ t $分布です。

Mはどの程度の大きさが良いか

Rubin（1987, 114）では、無限回の$ M $と比較して、欠損している情報の割合が50％までであれば、代入回数が2回であっても有限なMとの漸近的な相対有効性(RE)はおよそ90％であるとしています。
多くの文献では、$M=5$ （欠損割合が50％で、REが95％）であれば有効な推測を行うのに十分であるとしています。
一般的には、欠損している情報量だけではなく、分析モデルとデータ数を考量して、代入回数を決定する必要があります。
$ M $が50以上であれば、安定した結果を得られるとしている研究もあります。
代入計算を行う際の実効性とサンプリング誤差の低減を考慮して、Stataでは少なくとも20回の代入を行うことをお勧めします。

欠損データに関する仮定

多重代入では、データが欠損するメカニズムについて仮定を設けていませんが、ほとんどの代入方法では欠損のメカニズムが無視可能な（ignorable）ものであることを前提としています。
データが欠損するメカニズムが他のデータに依らない、missing completely at random（MCAR）であれば、欠損値を削除しても推定結果の一致性と有効性は保たれます。
データが欠損するメカニズムは他のデータに依存するmissing at random（MAR）の場合、リストワイズ削除を行うと、推定結果にバイアスが生じます。
欠損のメカニズムがMARで、モデルパラメータと欠損パターンのパラメータが異なる場合、欠損パターンを考慮する必要はなく、代入モデルを単純化できます。
実際の分析では、欠損メカニズムが無視可能であるかの判別は難しく、上記の仮定を用いるには慎重な判断が必要です。

データの欠損パターン

$ N \times p $のデータ行列$ Y=(Y_{1}, Y_{2}, \dots , Y_{p}) $で、変数は$ p $個で、観測値は$ N $個を考えます。行列$ R_{1} $が$ Y $の欠損パターンを表しているとします。

\[ R_{1}= \begin{pmatrix} 1 & 1 & 1 \\ 0 & 0 & 1 \\ 0 & 1 & 1 \\ 0 & 1 & 1 \end{pmatrix} \]

$ Y_{j} $が観測値$ i $において観測されれば$ R_{1}^{ij} $は1、観測出来なければ0です。$ R_{1} $の1列目と3列目を入れ替えると、$ Y $の欠損パターンはモノトーン（単調）ということができます。

\[ R_{1}= \begin{pmatrix} 1 & 1 & 1 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 0 \end{pmatrix} \]

行を上のように並べ替えると、欠損値の単調性がさらに明らかになります。
これに対して、モノトーンではないパターンは次ようになります。$ R_{2} $の最初の2列では、片方の列の欠損値がもう一方を説明するようなパターンが見られません。

\[ R_{2}= \begin{pmatrix} 1 & 1 & 1 \\ 1 & 1 & 0 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{pmatrix} \]

このモノトーンパターンは代入手順を大幅に単純化します。欠損パターンがモノトーンでの多変量代入は、独立した単変量代入の繰り返しに置き換えることができ、代入モデルを柔軟に作成できます。

適切な代入方法

多重代入は次の場合、妥当性を担保できます。

1. 多重代入の推定量$ \hat{Q}_{MI} $が漸近的に平均$ \hat{Q} $の正規分布と分散共分散推定$ B $に従う。
2. 代入内の分散推定$ W $は、$ Var(\hat{Q}_{MI}) $の分散共分散推定$ U $の一致推定である。

実際には、これらの基準で妥当性を判断することは難しく、Rubin（1987, 125-127）は代入モデルと欠損メカニズムのベイジアンモデルの事後分布を利用して代入することを提案しています。
予測平均マッチングと連鎖方程式を除き、Stataで利用できる代入手法は、欠損データのベイズ統計による事後推定分布を使用します。

代入モデルは、欠損データのメカニズムに関連する予測因子を含む必要があります。
例えば、分析モデルで2つの変数が相関していることがわかり、代入モデルでそれらのいずれかを削除してしまうと、相関の推定値に0に向かうバイアスが生じます。
また、分析モデルのアウトカム変数を代入モデルで使用しないことも推定値にバイアスを生じさせます。
サーベイデータでは、加重やストラータ、クラスタIDなどの構造変数も代入モデルに含めなければなりません。
一般的には、分析モデルで使用する全ての変数は、代入モデルにも含めるべきです。

不適切な代入モデルがどの程度の影響を与えるかは、観測データに対応する代入データの比率によります。
代入した値の割合が小さければ、推定結果に大きな影響を与えることはありません。

多重代入データの分析

代入データが準備できましたら、それぞれのデータセットで推定を行い、結果をRubinのcombination ruleに基づいて統合します。
推定段階では、様々な推定方法を使用できますが、それぞれの推定方法は本来、結果を統合すること想定していません。
このため、尤度比検定や適合度検定などの推定後機能は多重代入ではそのまま使用できません。

Stataでの多重代入

Stataの多重代入では、代入から結果のプーリングまで、複数の変数を扱う様々な代入方法を利用できます。
分析とプーリング段階では、mi estimateコマンドで推定した結果を統合します。
一般的な推定方法で係数を統合します。
解析段階では、一般的な推定のみならず、代入データの管理や診断を行うことができます。
ここでは、心臓発作のデータの例を使用して、推定を行います。
目的は、

1. bmiの欠損値を、線形回帰で代入
2. 代入したデータを、mi estimateを利用してロジスティック回帰で分析

の2つです。
これらの2つのステップを行う前に、データをmiで使用できるように準備します。
最初に、データをmiデータであることを宣言します。

ここでは、メモリを効率的に利用できるようにmlong形式を使用します。
代入を行うmi imputeのために、まず、代入変数を登録します。
一般には、mi registerにimputed, passiveまたはregularを使用して、分析に関連するすべての変数を登録することをお勧めします。

シミュレーションエラーを抑えるため、add(20)で20個の代入を作成します。
また、再現性のため、rseed()オプションで乱数シード値を設定します。

22個の欠損しているbmiの値が代入されました。
代入したデータと元データを比べて確認するために、元データ（m=0）と最初の代入データ（m=1）および最後の代入データ（m=20）に対して、mi xeqでsummarizeを実行します。

代入データは理想的な形に見えますので、mi estimateプリフィックスコマンドでロジットモデルにフィットします。

リストワイズ削除を使用した、上記のlogitモデルと比較し、欠損によって検出できなくなっていたageの有意性を確認できました。

まとめ

多重代入はシミュレーションベースの手法です。
個々の欠損値を真の値に近づけるのではなく、妥当性のある統計推測を行うことを目的としています。
多重代入は、

1. 欠損データのメカニズムに対して、代入手法が適切であるとき
2. 完全データの分析が適切であるとき

に、有効な推測が可能です。

欠損データの割合が小さい場合は、代入回数が小さくても（5から20）問題はありません。
この割合が高かったり、またはデータ構造によっては、最大で100回（またはそれ以上）の代入を用意する必要があります。

代入回数が小さい場合、多重代入の分布はStudentの$ t $（または$ F $）に従いします。
残差の自由度は$ M $と欠損情報の比率に依存します。
代入回数が大きい場合、分布は正規分布(または$ \chi ^{2} $)に近づきます。

代入モデルで使用する変数が分析モデルより制限されている場合、多重代入の前提が満たされず、有効な推測ができません。
一方、分析モデルが制限されていれば、推測は有効ではありますが、保守的になります。
分析の前提が異なる場合、その結果はバイアスしてしまいます。
多重代入の枠組みでは、尤度や偏差など特定の指標を正しく評価できません。
このような推定後機能は、多重代入の結果には直接適用できません。

代入データのスタイル

miデータには、flongsep、flong、mlong、およびwideの4種類があります。

• flong
  m=0,m=1,…,m=Mの代入データが個別の.dtaデータセットとなります。
  m=0がpat.dtaという名前で保存されれば、m=1は_1_pat.dta、m=2は_2_pat.dtaとなります。
  flongsepはfull long and separateの略です。
• flong
  m=0,m=1,…,m=Mの代入データを、_N=N+M×Nの単独のデータセットにします。
  NはM=0の観測数です。flongはfull longの略です。
• mlong
  m=0,m=1,…,m=Mの代入データをN=N+M×nの単独のデータセットにします。
  nはm=0における欠損値を含む観測値の数です。
• wide
  m=0,m=1,…,m=Mの代入データを_N=Nの単独のデータセットにします。
  多重代入に使用する全ての変数にM個の変数が追加されます。
  変数bpがm=0で欠損値を含んでいる場合、m=1では変数_1_bp、m=2では変数_2_bpを作成します。

まとめ

miデータは4つのフォーマットのいずれかで保存します。
これらの詳細は、上記をご覧ください。
miデータにはm=1,2,…,$ M $の番号を持つ$ M$ 個の代入が含まれ、m=0は欠損値を含むオリジナルのデータです。
miデータの変数は、imputed、passiveまたはregularにいずれかとして登録されます。これら以外はunregisteredとなります。
unregisterの変数はほぼregularの変数と同じように扱われます。

　　　　　

regularの変数は通常、欠損値を含みません。
欠損値がある場合は、m>0において代入が行われません。

　　　　　

imputed変数はm=0において欠損値を含みます。
m>0において代入されます。

　　　　　

passive変数はimputed、regular、および他のpassive変数の組み合わせです。
imputed変数がm=0において、.より大きい値、.a、.bから.zを含む場合、これらの値は、m>0においても欠損値のままです。

参考文献

Rubin, D. B. 1987. Multiple Imputation for Nonresponse in Surveys. New York: Wiley.

Rubin, D. B. 1996. Multiple imputation after 18+ years. Journal of the American Statistical Association 91: 473-489.