こんにちは、すきとほる疫学徒です。
R Tipsでは、私自身の備忘録も兼ねて、おすすめのRパッケージを紹介します。
80％くらいは備忘録が目的なので、殴り書きになりますことをお許しください。
多重代入法とは？
miceは、Rで多重代入法を使用するためのパッケージです。
多重代入とは、シンプルに言うと「欠損している変数を、欠損していない他の変数で予測してやり、埋めちゃおうぜ」という欠損の対処法です。
欠損に対しては、
・欠損を持つ対象を全て削除するComplete case analyses
・欠損に平均値を代入する方法
・多重代入法
などさまざまな選択肢があり得ます。
しかしながら、欠損を削除、平均値などの一意の値を代入するなどの方法は、バイアスに繋がることが分かっており、医学研究においては使用が推奨されていません。
そんな中で人気を博してきたのが今回紹介する多重代入法です1。
今回は、多重代入法の説明ではなく、それをRで行うためのパッケージの解説ですので、多重代入法自体には深入りしません。
さっくり説明しますと、
多重代入法をイラストにするとこんな感じ。
なぜ”代入”かと言うと、欠損していないデータから、欠損しているデータを予測して代入するからです。
また、なぜ”多重”かと言うと、欠損による推定の不安定性を結果に反映させるため、図のように複数データセットを作成し、それぞれで代入・解析を行い、最終的に結果を統合するからです。
私の身近で研究をしている人の中で、なぜか完全に欠損を除外したComplete case analysesをやる方々が少なくありませんでした。
Complete case analysesは選択バイアスをもたらすリスクがあり、またサンプルサイズを減らすことで解析精度が低下しますので、私は基本的には用いるべきではないと思っています。
ですので、「どうしてComplete case analysesをするんですか？」と聞いてみたのです。
すると、
「他にどんな方法があるか知らない」、
「多重代入法というのがあるのは知ってるけれど、どうすれば良いかわからない」
という声がけっこう聞こえてきました。
そんなことから、今回はRで多重代入法を行うパッケージであるmiceを紹介することにしたのでした。
ちなみに、多重代入法自体の解説はこちらのBMJの論文に丁寧に書かれているので、併せててどうぞ。
https://www.bmj.com/content/338/bmj.b2393
また、欠損に対してなんでもかんでもmiceを使えば良いと言うわけではなく、まず行うべきは欠損の特徴を把握することです。
その方法はこちらの記事で紹介していますので、後ほどご覧ください。
mice()の使い方
まず、パッケージをインストールし、読み込みます。
install.packages("mice")
library(mice)
そうしたら、次は欠損の代入を行います。
df_miced <- mice(df,
m = 20,
maxit = 50,
method = "pmm",
seed = 1234)
dfで代入対象のデータを指定、
mは代入の結果として作成するデータセットの数（私はいつも20）、
maxitは代入の反復回数（私はいつも50）、
methodは代入方法の指定、
そしてseedは代入のための乱数設定です。
なぜseedを設定する必要があるかというと、この乱数を決めておかないと、解析を回すたびにランダムに乱数が決まってしまい、毎回代入結果が変わってしまうからです。
すると、解析の再現可能性を確保することができませんね。
miceでは、pmm（Predictive mean matching）以外にも様々な代入方法を選択可能ですが、今回の記事の趣旨は多重代入法自体の説明ではないため、そこは割愛します。
余談ですが、mice()のデフォルト設定では、dfに含まれる全ての変数が多重代入法に使用されます。
状況によっては、「dfのうち、この変数はmice()には使いたくない」ということがあるかもしれません。
そんな時は、subset()を使って、以下のように指定してあげればOKです。
df_miced <- mice(subset(df,
select = c(選びたい変数),
m = 20,
maxit = 50,
method = "pmm",
seed = 1234))
さて、これで欠損が代入された複数のデータセットができました。
次は、それぞれのデータセットごとに解析を行いましょう。
今回は、重回帰分析を行うことにします。
miced_lm <- with(data = df_miced,
lm(outcome ~ variable1 + variable2 + variable3))
with()パッケージを用いますが、冒頭で分析対象とするデータセットを指定している以外は、lm()で通常の重回帰分析の式を書いているだけですね。
余談ですが、mice()は適応可能な分析パッケージと、適応不可能なパッケージがあります。
lm()やglm()などのメジャーなパッケージは問題なく使用できます。
しかし、混合効果モデルのためのlmer()や、一般化推定方程式のためのgee()などをwith()の中に使ってしまうと、残念ながらエラーとなります。
幸いなことに、lmer()に関してはbroom.mixed()を使うことで、適応可能になります。
install.packages("broom.mixed")
library(broom.mixed)
医療大規模データベース研究をしていると、地域や病院のようなデータの階層性を考慮せねばならない状況が多発するので、broom.mixed()によってmice()とlmer()が共存できることを知った時には、開発者への感謝がとまりませんでした。
mice()と共存していないパッケージに関しては、どうやら代入されたそれぞれのデータセットに対して分析を行い、最終的には手計算（もしくは自作関数）で結果の統合を行わねばならないようです（mice()はルービンズルールと呼ばれる方法で結果の統合を行っていますが、それを手計算でやるということですね）。
「そんなはずはないだろ！」と現実を受け入れきれず、長らくネットの海を漂ってきました。
しかし、「これだ！」といった回答に出会うことができず、またTwitterの猛者たちも結局は自作関数で対処しているとのことだったので、ここは諦めて自力でやるしかないようです。
さて、with()によって代入された各データセットごとに解析を行いました。
次はいよいよ結果の統合です。
miced_pool <- pool(miced_lm)
summary(miced_pool)
このように、pool()に入れ込んであげるだけですね。
あとはsummary()を使えば、通常lm()をかけた時と同じような結果が表示されます。
mice()を使用する際の注意点
投入する変数を必要最低限にする
mice()は解析負荷が重く、医療大規模データベース研究のようにサンプルサイズが100万を超えることもあるようなデータでは、極めて不安定になります。
まず、時間が数時間から数日かかりますし、すぐにR自体がフリーズします。
3日間かけて解析し、朝起きていざパソコンを見たらフリーズしていたなんてことになったら、もう立ち直れません。
ですので、少しでも解析負荷を減らすための工夫が必要になります。
多重代入に使用しない変数はそもそもデータセットから除外するか、それか先ほど紹介したsubset()から除外するようにしましょう。
少ないデータセット、反復回数でmice()が回るか検証する
mice()、with()、pool()と回し、最後の最後で「やべ、必要な変数忘れた」とか、「あ、変数のカテゴリーミスった」とかケアレスミスが発覚したら、mice()を回す前の自分を呪い殺したくなってしまいます。。。
ですので、まずは「解析さいごまでちゃんと回るかな？」トライアルとして、少なめのデータセット、反復回数でmice()を回し、pool()までもっていきましょう。
私はいつも適当に、
m = 3
maxit = 3
とかにして、トライアルをやっています。
また、トライアルではmice()後に作成された複数データセットから、一つを抽出して実際に欲しい形で代入が行われているかを確認した方が良いでしょう。
complete()パッケージを使えば、代入済みの複数データセットから任意のデータセットを取り出すことができます。
miced_comp <- complete(df_miced, action = 1)
“action =”によって何番目のデータセットを抽出するかを指定します。
その後、miced_compは通常のデータセットと同様に扱うことができますので、
「本当に欠損が埋まっているか」
「おかしな挙動をしている変数がないか」
などと確認してやると良いでしょう。
多重代入したデータセットを保存しよう
先ほどもお伝えした通り、mice()は解析負荷が重く、場合によっては一つのmice()だけで数日間かかることもあります。
サブグループ解析、感度分析を含めてトータルの解析回数が10回になるなんてことも医療大規模データベース研究では珍しくありませんが、となると1回で数日かかるmice()を10回も回さねばならなくなります。
別の記事で「Rでは基本は加工したデータそのものは保存せず、加工するための設計図を保存する」とお伝えしましたが、今回は例外にあたります。
mice()を終えた後は、何よりもまず先に多重代入したデータセット自体を保存しましょう。
保存方法は、
save(data, file = "パス/名前.Rdata")
です。
fileには保存先をしていするため、まず保存先のパスを書き、その後にデータ名を書きます。
この時、Rのデータセットのデータ形式である”.Rdata”を忘れずに付けるようにしましょう。
mice()が終わった直後に保存できるように、mice()を回す際にはセットでsave()も回すようにした方が良いです。
なお、保存したデータを読み込むには、
load(file = "パス/名前.Rdata")
でOKです。
終わりに
いかがでしたでしょうか？
「多重代入法ってなんだか難しい！」と思って敬遠されていた方もいらっしゃるかもしれませんが、コードだけ見れば実装はさほど難しくないと思えてきませんか？
ぜひぜひ、お手元のデータで試してみてください。
なお、繰り返しますが本記事は多重代入法をRで実装する方法をお伝えすることを目的としており、多重代入法自体の解説を行なっておりません。
ですので、もし「多重代入法そのものがまだわからん」という方がいらっしゃいましたら、ぜひともそちらも資料を読み込んでみてくださいませ。
すきとほる疫学徒からのお願い
本ブログは、読者の方が自由に記事の金額を決められるPay What You Want方式を採用しています。
「勉強になった！」、「次も読みたい！」と本ブログに価値を感じてくださった場合は、以下のボタンをクリックし、ご自身が感じた価値に見合うだけの寄付を頂戴できますと幸いです。
もちろん価値を感じなかった方、また学生さんなど金銭的に厳しい状況にある方からのご寄付は不要です。
引き続き情報発信していく活力になりますので、ぜひお気持ちに反しない範囲でご寄付をお願い致します！
