ちょっとわけあって，欠損値の処理について勉強するソルジャー業務機会がありました。そこで，多重代入法（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 = 語彙テストの得点

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

dat1 <-dat[,-1]

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

install.packages(“VIM”)

library(VIM)

#欠損情報を入手

m <-aggr(dat1)

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

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

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

 

#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)

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

dat1.out$covMatrices

#補完データを書き出し

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

#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.out, ols.out$coef)

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

}

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

b.out

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

se.out

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

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

combined.results

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

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

補足

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

install.packages(“mice”)

library(mice)

imp_data<-mice(dat1,m = 5)

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

summary(pool(fit))

まぁまぁ結果は近い？

なんかmiceパッケージの方が使い勝手が良さそうですね（ｱﾚﾚ

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

今日はこのへんで。

なにをゆう　たむらゆう

おしまい。

先日ついにtidyrパッケージで縦横変換をやってみた。今までずっと覚えなきゃなと思ってて結局億劫でやってなかったんだけれど手伝いで仕方なく。

データハンドリングで，横に広がる（列数の多い）wide-formatのdataを分析に便利なlong-formatに変換するというのは，分析をやるまえのデータ整形の過程で必ず通らなければならない道。Excel上で地道にやるのもあるにはあるだろうけど，Rのtidyrを使えばすぐできるというお話。これ系の話はググればネット上にゴロゴロ記事が転がっている。なので，今回は具体的に外国語教育研究系のデータ分析ではどんなときに役に立つかなーという話です。特に無駄に変数が多い時なんかは結構役に立つと思います。

sample data

適当にこんな感じでデータが入っているとします（数値は乱数）。事前・事後・遅延事後のテストがあって，それぞれで流暢さ（fluency），正確さ（accuracy），複雑さ（complexity）の指標があるみたいな。1列目には被験者番号があって，それぞれの列名は”test.measure”という感じでドット区切りにしてあります。で，これを，「事前・事後・遅延事後」の3つのカテゴリカル変数からなる”test”列と，3つの指標（CAF）の”measure”列の2列に分解しましょうという話し。

R上に上のようなデータを読み込みます。

read.table関数で読み込んだ場合（※MacなのでちょっとWindowsと違います）

まずはパッケージの準備

install.packages(“tidyr”) #パッケージのインストール

library(tidyr) #tidyrパッケージの呼び出し

datに入っているデータフレームをまずは，変数名の列（variable）とその値（value）の列に集約します。key=変数名の列，value=値の列という引数指定です。

#%>%はパイプ演算子。以下の関数で扱うデータフレームを指定するということです

dat %>%

#-subjectでsubjectの列を除外。変換したdataをdat2に入れる

gather(key = variable,value = value,-subject) -> dat2

long formatのデータに変換されました

これで変数名を1列に集約できました。ただし，この列にはpre, post, delayedというテスト実施時期の情報と，fluency, accuracy, complexityという指標の情報が混ざってしまっています。これを，2列に分けます。

dat2 %>%

#colで分けたい列名，intoで分けた後の列名を指定。sepでセパレータを指定しますが，デフォルトはあらゆる非アルファベットになっているので，ドットなら指定しなくても大丈夫

separate(col = variable, into = c(‘test’, ‘measure’)) -> dat 3

さて，dat3の中身を確認してみると…

long formatのデータに変換されました

おおおおー！！！！！

というわけで，横に変数がたくさんある場合もこのtidyrのgatherとseparateを使えば簡単にlong型に変換できそうです。

参考サイト：http://meme.biology.tohoku.ac.jp/students/iwasaki/rstats/tidyr.html

なにをゆう　たむらゆう

おしまい。

ストレスが溜まっているのかなんなのか帰り道にセブン-イレブンでチーズとワインを買って飲むとかいう珍しいことをしているという言い訳を最初に書いておこうと思います。それから，自戒もたっぷり込めていることも申し添えておきたいと思います。

データ分析の方法の相談についてです。色んな所で起こっていることだと思うし色んな人がきっと同じようなことを考えているだろうと思います。一言でいうと，データを取る前に相談してみましょうということです。データを取ったあとに，「このデータってどうやって分析すればいいですか？」っていう質問は結構つらいところがあったりします（恐れ多くもなぜだかこういう質問をされる立場に私もなってしまいました）。なぜなら，本当にその質問者の方が調べたいこと，明らかにしたいことが，そのデータ収集方法で明らかにすることができるのかちょっと怪しかったりすることがあるからです。

そういう状況にならないために重要なことは，研究課題を見つけて，それを具体的な実験計画なりに落としこむ際に，その先のことまでちょっと考えてからデータ収集をしてみることだと思います。つまり，

こういう方法を使うと，こういうデータを収集することができる。それをこういう分析にかければ，こういう結果が得られるはず。もしこういう結果が得られた場合，それはこういうことを意味していると解釈できるだろう。

みたいなことを考えてから実験を組んで，データを取っていただきたいなということです。そうすれば，データを取ってからどうしたらいいかという質問が来ることはなくなるからです。上述の事でイマイチ自分がうまく説明できないとなったら，その時点で相談をすればよいと思います。そうすれば，実はその方法で得られたデータで，質問者の方の主張をするのはちょっと違うので，こういうデータ収集をしないといけませんよ，なんてアドバイスができたりするかもしれないからです。これって仮説検証型の研究でも探索的な研究でも同じだと思います。たとえ探索的な研究だったとしても，得られた結果に対してなんらかの解釈を加える必要はあるはずですから。研究課題がwh疑問文で始まる形であったとしても，結果に対しては，きっとこういう理由でこうだったのだろうということを言えなければいけないわけですし。

データをとって，分析をしてみて，これであってるのか不安だから相談してみようということもあるのかもしれませんが，それって別にデータ取る前にもできることなんですよね。

こういう研究課題をもっていて，これを明らかにするためにこういう分析をしてこういう結果が得られたらこういう解釈ができると思うのですが，これで合ってますか？

ってデータ取る前に相談すればすむと思うのです。なのにどうして，「データを取ったあとにこれってどうやって分析すればいいんですか？」ってなってしまうのでしょう。

これはあくまで持論なのですが，研究って，研究のデザインとか実験計画ですべて決まると私は思っています。面白いなと思う研究はやっぱりデザインがしっかりしていて，分析の結果から得られたことの解釈がとってもシンプル（注1)なんです（そういう研究が良い研究だと思う私からすると，discussionを長々と書いている論文やそれが良しとされる風潮，また「discussionが短い」みたいなコメントってなんだかなと思ったりします）。それって，デザインを練る段階で色んなことを考えているからなんですよね。研究の計画を立てる段階で，データの分析方法まで明確にイメージできない場合は，データ収集に入るべきではないと私は考えています。私の場合は基本的に，ある研究で行われることを考えることは，R上で分析を行う際のデータセットの形をイメージすることとほぼ同義です。変数のデータフレームがイメージできていて，それをどうやって分析すれば何をどうしていることになるのか，こういう結果になったらそれはどう解釈するのか，また別の結果がでたらそれをどう解釈するのかまで考えられて初めて研究のデータ収集ができるのではないでしょうか。

こういった考えが当たり前になれば「統計屋さん」（私は違いますけど）とか言われる人たちもずいぶん楽になるのではないかなと思っています。

なにをゆう　たむらゆう

おしまい。

*注1: だからといって，t検定一発で明らかになることってそんなにはないとは思います。少なくとも今なら昔はt検定一発だったこともGLMMするでしょう。例えば正答率の平均値をt検定というのはもはや時代遅れになりつつあるわけですし。