miceパッケージampute関数による欠測データの作成:(1)原理

はじめに

前回、混合モデルの説明に使用した欠測データの作成のため、Rのパッケージの一つ、miceパッケージ（参考資料1）に実装されている関数の「ampute」（参考資料2）を使用しました。

今回は、ampute関数の原理と使い方をもう少し詳しく見ていきたいと思います。

miceパッケージ

このパッケージは、欠測データの補完方法である「多重代入法(Multiple Imputation)」の一種である、Full Conditional Specification(FCS)法を実行するために用いられます。

元々、多重代入法の理論では、欠測が生じている変数も含めたすべての変数*1の同時確率分布を想定し、その上で欠測のない変数による条件付けのもとで欠測の生じた変数のサンプリングを行うこととされていました。
ところが具体的に同時確率分布を特定することが難しいケースもある*2ことから、補完したい変数一つ一つの「他のすべての変数で条件付けした場合の」確率分布を考えて、変数ごとに補完値を発生させるというFCS法が提案されました。FCS法の代表的なアルゴリズムが、パッケージの名称になっているMICE(Multiple Imputation by Chained Equation)です（参考資料3）。

miceパッケージは主に上記のFCS法によって欠測値を補完するために用いられますが、逆に「欠測を発生させる」という機能もあります。それが、以下で紹介するampute関数です。

一部の値を欠測にさせるampute関数

ampute関数は、データフレームに含まれる任意の変数のうち、一部の値を「NA」に置き換えます。つまり、人工的に欠測を生じさせます。

f:id:mstour:20210711195412j:plain

とりあえずデフォルトの設定通りに欠測を発生させるには、amputeの引数にデータフレームを入れます。返される結果の「amp」に欠測を発生させたデータフレームが格納されます。

# amputeの結果をresultに格納します。
# dataframeはデータフレーム名です。
result = ampute(dataframe)
# 作成された欠測データをmdに格納します。
md = result$amp

ですが、通常は欠測させたい変数、当該変数に欠測が生じる個体*3の割合、欠測パターン*4、欠測メカニズム（MCAR、MAR、MNAR）などを明確に指定したいことが多いため、細かな設定が必要となります。

参考資料4に沿って、ampute関数で用いられている欠測発生の原理を確認していきましょう。

データに欠測値がある場合の統計手法の性能を確認するため*5に、人為的に欠測を発生させたデータを用いてシミュレーションを行うことがあります。
欠測の起こる個体の割合についていくつかの異なるシナリオを考えたいような時には、狙った割合になるようにシミュレーションデータを作成できると嬉しいです。欠測の起こる個体の割合が10%・20%・30%それぞれの場合の解析結果への影響を見積もることで、予定している解析方法を使うことに問題はないかの判断材料が得られるかもしれません。

この時、欠測させる変数の数が1つだけであれば、それほど困難は発生しないでしょう。例えばロジスティック回帰を用いて欠測発生のモデルを作ることにすると、当該変数が欠測する確率が狙った値（例えば30%）になるようにパラメータを調節し、推定された各個体の欠測確率に基づいてデータを欠測させるかどうかを決めれば、おおむね30%の個体に欠測があるデータが得られるはずです。

しかし同様の方法で複数の変数を欠測させる場合、変数ごとの欠測割合をコントロールできても、「少なくとも一つの変数に欠測が生じた個体の割合」のように変数の組み合わせにおける欠測割合をコントロールするのは困難であることが知られています。そこで、複数の変数を同時に欠測させる「Multivariate amputation」というアプローチが考えられています。ampute関数は、この考え方に基づいています。

Multivariate amputation

欠測のある変数が1つだけならば、当然ですが欠測パターンは「当該変数に欠測が生じる」という1通りしかありません。しかし複数の変数に欠測が生じうる場合には、欠測パターンも複数存在します。
例えばA, B, Cという3変数があって、このうちBとCが欠測しうるような場合、欠測パターンは

A:欠測なし、B:欠測あり、C:欠測なし
A:欠測なし、B:欠測なし、C:欠測あり
A:欠測なし、B:欠測あり、C:欠測あり

の3通りが考えられます（すべての変数が「欠測なし」というパターンは除く）*6。
このように欠測しうる変数が複数ある場合、個々の変数ごとに欠測確率のモデル*7を設定するアプローチでは、複数の変数が同時に欠測する割合を思うようにコントロールするのが難しくなります。
そこで、Multivariate amputationでは、「個々の変数」ではなく「起こりうる欠測パターン」ごとに欠測モデルを設定する方針を取ります。

Multivariate amputationは、下記の模式図のように

データセットを「起こりうる欠測パターン」ごとに分割し、
パターンごとの部分データセットにて欠測を人工的に発生させ、
最後に各部分データセットを結合する、

という手順をとります。

f:id:mstour:20210711195433j:plain — Schouten et al.(2018)をもとに筆者作成

図の左から順に、手順を見ていきましょう。
なお、欠測を生じさせる元のデータセットには欠測がない（完全データ）ものとします。

1. 起こりうる欠測パターンごとにデータセットを分割する

まず、用意したデータセットを、起こりうる欠測パターン*8ごとに分割します。

例として、変数A, B, Cが存在するデータセットにて、以下の3つの欠測パターンを生じさせたい場合を考えましょう。

A:欠測なし、B:欠測あり、C:欠測なし　 (1)
A:欠測なし、B:欠測なし、C:欠測あり　 (2)
A:欠測なし、B:欠測あり、C:欠測あり　(3)

はじめに、欠測の起こる個体において(1)〜(3)の欠測パターンがそれぞれ何%を占めるかを決めます。この割合をampute関数では「freq」といい、各パターンのfreqの値に応じてデータセットを分割します（すべての変数が欠測なしというパターンは含みません）。
例えば欠測パターン(1)が40%、(2)が40%、(3)が20%を占めるとすると、(1)のパターンには元のデータセットの40%のレコードからなる部分データセットが割り当てられます（2と3も同様）。

2. それぞれの変数が欠測確率へ与える影響の大きさを設定する

続いて、欠測パターンごと（1.で分割した部分データセットごと）に、各変数が欠測へ影響する度合いを設定します。

欠測データ解析の理論では、ある変数が欠測する確率は何らかの変数の値に応じて決まると考えます。例えばMARは、欠測する確率は他の何らかの変数に影響されると想定しており、MNARは欠測した変数自身の値も欠測確率に影響する可能性を考慮に入れます。
そこで、各個体 $i$ の欠測確率が「欠測に関与する変数の値 $y$ 」と「その変数の影響度 $w$ 」を掛けたものからなる以下の「Weighted sum score」によって決まるものと考えます。

$\displaystyle wss_i = w_1 y_{1i} + w_2 y_{2i} + \cdots + w_m y_{mi}$

上式の影響度（重み） $w_1, \cdots, w_m$ の設定に基づき、ロジスティックモデルを通して各個体の欠測確率が算出されます*9。次のステップにて、実際に欠測を発生させます。

なお、欠測させたい変数が $y_{1i}$ だとすると、

すべての重みが0・・・ $y_{1i}$ の欠測メカニズムはMCAR
重み $w_1$ が0で、他の重みの少なくとも1つは0でない・・・ $y_{1i}$ の欠測メカニズムはMAR
重み $w_1$ が0ではない・・・ $y_{1i}$ の欠測メカニズムはMNAR

という状況に対応します。

各変数の重みは、ampute関数では「weights」により指定します。

3. 指定した欠測の発生割合に基づいて欠測を発生させる

ampute関数では、何らかの変数が欠測する個体の割合を「prop」で指定しておきます。
Weighted sum scoreの算出後、欠測パターンに応じたそれぞれの部分データセットに対して、propで指定した割合の個体に欠測を発生させます。
どの個体に欠測を発生させるかは、Weighted sum scoreの値からロジスティックモデルによって算出した欠測確率に基づき決定されます*10。

最後に、欠測を発生させた部分データセットを統合することで、元のデータセットに欠測を生じさせた新たなデータセットが完成します。

以上のステップによって、複数の変数を対象とする場合にも意図した割合の欠測を容易に発生させることができます。
ではいよいよRでの実践を・・・といきたいところですが、長くなってしまったので、続きは次の記事に回すことにします。

おわりに

欠測を人為的に発生させるための方法である「miceパッケージ・ampute関数」について、まずはどのような原理に基づいているかを紹介しました。
自分で欠測メカニズムのモデルを（例えばロジスティック回帰モデルで）定義するとなると、各変数が欠測確率へ与える影響を表すパラメータの値をうまく調節する必要があります。ampute関数ではそのような複雑な作業をせずとも、全データのうちどの程度に欠測を発生させたいかを直接指定すればよいため*11、とても使い勝手が良いと考えています。

次回は実際にRを動かしてより具体的な話をしていきたいと思います。

今回は以下の資料を参考にしました。
[1] Package ‘mice’ Version 3.13.0.
https://cran.r-project.org/web/packages/mice/mice.pdf
[2] ampute: Generate missing data for simulation purposes.
https://www.rdocumentation.org/packages/mice/versions/3.13.0/topics/ampute
[3] 高井啓二・星野崇宏・野間久史(2016), "欠測データの統計科学", 岩波書店.
[4] Schouten, R. M., Lugtig, P., & Vink, G. (2018). Generating missing values for simulation purposes: a multivariate amputation procedure. Journal of Statistical Computation and Simulation, 88(15), 2909-2930.
https://www.tandfonline.com/doi/full/10.1080/00949655.2018.1491577
[5] Multivariate amputation using ampute.
https://mran.microsoft.com/snapshot/2017-02-20/web/packages/mice/vignettes/ampute.html

*1:ここで「すべて」というのは、「欠測を補完するプロセスにおいて登場するすべての変数」を意味します。例えば個人の識別番号のように、補完の対象とならない・あるいは欠測を説明する役には立たないような変数は含まれません。

*2:すべての変数に正規分布を仮定した場合には同時確率分布も正規分布であることがわかりますが、例えば連続型分布と離散型分布が混ざっているような場合には具体的な同時確率分布を定義することは難しそうです。

*3:以降、解析の単位は「個体」とします。またデータセットには1個体につき1レコードのみ格納されているこものとします。

*4:欠測する変数の組み合わせをまとめたものです。例えば変数A, B, Cが同時に欠測しうる場合、(A, B, C)の組が欠測パターンの1つの要素となります。

*5:例えば、MARではない欠測が起こった場合に、欠測がMARであると仮定する混合効果モデルで解析するとどのくらい結果に悪影響があるか、などを知りたいものとします。

*6:ここでは可能性として考えられるパターンを挙げていますが、人工的に欠測データを作る場合に必ずしもすべてのパターンを含めないといけないわけではありません。また現実のデータでも、必ずしも3通りすべての欠測パターンが存在しているとは限りません。

*7:MCAR以外の欠測メカニズムを考えるには、欠測する個体の割合だけでなく、他の変数が欠測確率に与える影響の大きさや、他の変数と欠測確率とを結びつける関数の形などを特定する必要があります。

*8:人工的に発生させたい欠測パターンのことです。このセクションの最初に挙げた例だと、「Bのみ欠測」「Cのみ欠測」「BとCの両方が欠測」という3パターンに対応する部分データセットがそれぞれ作成されます。

*9:ampute関数ではロジスティックモデルしか実装されていないようですが、モデルの概形（右上がりか右下がりかなど）は「type」にて指定できます。

*10:欠測する個体の割合が意図したものになるように、propで指定した値に応じてロジスティックモデルの位置は調整されます。

*11:とはいえ、その他にもいくつかの設定は必要です。