Rで多重代入法:Ameliaパッケージ

ちょっとわけあって,欠損値の処理について勉強するソルジャー業務機会がありました。そこで,多重代入法(MI: Multiple Imputation)という方法をRで実行する方法を少しかじったのでメモ代わりに残しておきます。

ちなみに,欠損値の分析をどうするかという話は全部すっ飛ばしますのでそのあたりは下記リンクなどをご参照ください。

欠損値があるデータの分析

DARM勉強会第3回 (missing data analysis)

多重代入法に関してはこのあたりの資料をどうぞ。

多重代入法の書き方 公開用

様々な多重代入法アルゴリズムの比較(リンク先直PDF)

ざっくりと多重代入法はなにをやるかというと,手法によってアルゴリズムは異なりますが,欠損値を推定して補完したデータをたくさん作って,そのデータを元にして行った分析(e.g., 回帰分析なり分散分析なり)の結果得られたパラメタの推定値を統合して,欠損値がない場合の真のパラメタ値を推定しましょうみたいなことです。

Rには色々パッケージがあって,miパッケージ,normパッケージ,Ameliaパッケージ,miceパッケージ等々たくさんパッケージが用意されています。私はAmeliaパッケージを使ってみたので,その報告です。Ameliaパッケージは,bootstrapped EMアルゴリズムという手法を用いて欠損値補完を行っています。以下スクリプト例です。

#まずはデータセットの下準備。こんな仮想データだとします。

#id =参与者ID, reading = 読解テストの得点,grammar = 文法テストの得点,vocab = 語彙テストの得点

Screen Shot 2016-01-01 at 15.14.42

 #id列はいらないのでid列を抜いたデータセットを作る

dat1 <-dat[,-1]

Screen Shot 2016-01-01 at 15.15.31

#欠損を可視化するにはVIMというパッケージが便利

install.packages(“VIM”)

library(VIM)

#欠損情報を入手

m <-aggr(dat1)

#うまくいくとこんな感じで欠損情報を可視化してくれます

左側が欠損の割合で,右側が欠損のパターン。赤い所が欠損です。

左側が欠損の割合で,右側が欠損のパターン。赤い所が欠損です。

#数値で確認する場合はmの中身を見る

 Screen Shot 2016-01-01 at 15.17.48

#Multiple Imputation by a bootstrappted EM algorithm (Amelia Package)

#パッケージのインストール

install.packages(“Amelia”)

library(Amelia)

#amelia()関数をまずは使います

#引数は以下のとおり

#x = data.frame

#m = 何個のデータセットを作るかの指定

#dat1.outという変数に,多重代入した結果を入れます

dat1.out <-amelia(x = dat1,m = 5)

#中身はsummaryで確認します

summary(dat1.out)

Screen Shot 2016-01-01 at 15.19.32

#補完データの分散共分散行列をみる

dat1.out$covMatrices

Screen Shot 2016-01-01 at 15.20.12

#補完データを書き出し

#separate =Fと指定すると,データが1つのファイルに書きだされます

#separate = Tにすると,データセット1つにつき1ファイルで書きだされます

#dat1.ameria.csv, dat1.ameria1.csv, … dat1.ameria5.csvみたいな感じで番号をつけてくれます

#orig.data=Tで,オリジナルのデータを出力する際に含めるかどうかの指定ができます

write.amelia(obj = dat1.out,“dat1.amechan”,separate=T,orig.data = F)

#AmeliaView()を使うと,GUIで操作できます(重いです)

AmeliaView()

#多重代入した結果得られた補完データをまとめる作業をします

#例として重回帰分析をやってみます

#独立変数 = grammar, vocab

#従属変数 = reading

m =5

b.out <-NULL

se.out <-NULL

for(i in 1:m){

  ols.out <-lm(reading ~ grammar+vocab,data=dat1.out$imputations[[i]])

  b.out <-rbind(b.outols.out$coef)

  se.out <-rbind(se.out, coef(summary(ols.out))[,2])

}

#b.outにはbetaが5つ入っています

b.out

Screen Shot 2016-01-01 at 16.26.42

#se.outには標準誤差(SE)が入っています

se.out

Screen Shot 2016-01-01 at 16.26.54

#mi.meld()で,5つのbetaとSEをまとめます

combined.results<-mi.meld(q = b.out, se = se.out)

combined.results

Screen Shot 2016-01-01 at 16.28.26

SEがちょっと広めですかね

というわけで,まとめる段階が手動で関数書いていますが,それ以外は割りと簡単にできます。ファイルの書き出しなんかもできますし。

補足

miceパッケージを使うともっと簡単に多重代入->分析->統合ができるようです。

install.packages(“mice”)

library(mice)

imp_data<-mice(dat1,m = 5)

fit <-with(imp_data,lm(reading ~ grammar + vocab))

summary(pool(fit))

まぁまぁ結果は近い?

まぁまぁ結果は近い?

参考:Tokyo r #37 Rubin’s Rule

なんかmiceパッケージの方が使い勝手が良さそうですね(アレレ

ちょっとまた時間があったらもう少し勉強してみたいと思います。

今日はこのへんで。

なにをゆう たむらゆう

おしまい。

広告

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中