RでGLMMするときに使うglmer関数。要因の水準が3以上になると，どうしても多重比較したくなっちゃいますよね。けれど，デフォルトでは全水準の総当り比較はできない。そんなときにはlsmeans関数が便利です。ところが，glmer関数の出力結果をいれた変数を使ってlsmeans関数に渡してもなんかエラーメッセージが出てしまいました。

Error in format.default(nm[j], width = nchar(m[1, j]), just = "left") : 4 arguments passed to .Internal(nchar) which requires 3

こういうやつ。Rでエラーメッセージが出て，その原因がよくわかんなかったらもうそのエラーメッセージをそのままコピペしてググってしまうと。そうするとたいてい同じ悩みで困っている人の書き込みとかがヒットすることがしばしばあります。

するとこのサイトがヒット。

こんな感じ。

要因の順番は，多重比較したい要因が先にきます。つまり，上の例だと，condition条件の多重比較をpresentation条件ごとにやるという感じです。逆にすれば，presentation条件の多重比較をcondition条件ごとにやることになります。

もちろんコーディングを変えて全条件の組み合わせをつくればlsmeansとかしなくてもglmerだけで大丈夫ですし，例えば1, 2, 3と水準があったときに，1と2,3をまとめた水準を比較して違うか調べたいときとか（例えば英語母語話者が1で日本語母語話者が2, 韓国語母語話者が3で，母語話者群と学習者群が異なるか調べたいとか）には，やはり自分でコーディングを変えないといけません。カテゴリカル変数のコーディングの違いについてはこのサイトがわかりやすいです。

R Library: Contrast Coding Systems for categorical variables

しかしもはやこれでANOVAする必要が本当にないのではないかと…

そういえば，ついにdplyrの基礎がわかってなにこれめちゃ便利って感動したのでそのこともまた書きます。

なにをゆう　たむらゆう。

おしまい。

Nagoya.R #15の発表資料をslideshareにアップロードしていたのですが，どうもうまく表示されないので，speakerdeckにもアップロードしました。こちらだと正常に表示されるようです。

教えていただいたlangstat先生ありがとうございました。

なにをゆう　たむらゆう

おしまい。

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

なにをゆう　たむらゆう

おしまい。