今までデータをあれこれいじったり記述統計出したりするのはExcelで色々やってからRにいれて分析という流れでした。dplyr使ったら便利そう…というのはなんとなくわかっていたんですが，自分がよく扱うようなデータセットの構造に対してどうやって使えば今までやってたことが楽になるのかいまいちイメージがわかずにずっと放置していました。最近色々あってようやく基礎中の基礎だけを覚えたので以下メモ書きです。

まずはサンプルのデータを作ります

#事前・事後・遅延事後でCAFのデータを取ってるというデザイン
pre.fluency <-rnorm(50,mean=10,sd=2)
pre.accuracy<-rnorm(50,mean=7,sd=2)
pre.complexity<-rnorm(50,mean=4,sd=2)
post.fluency<-rnorm(50,mean=14,sd=4)
post.accuracy<-rnorm(50,mean=10,sd=3)
post.complexity<-rnorm(50,mean=7,sd=2)
delayed.fluency<-rnorm(50,mean=12,sd=2)
delayed.accuracy<-rnorm(50,mean=8,sd=1)
delayed.complexity<-rnorm(50,mean=6,sd=1)

#それぞれの列を横にくっつける
dat<-cbind(pre.fluency,pre.accuracy,pre.complexity,post.fluency,post.accuracy,post.complexity,delayed.fluency,delayed.accuracy,delayed.complexity)
dat<-as.data.frame(dat) #データフレーム型に変換
dat$subject<-rep(1:50) #実験参与者のID列をつける

---

とりあえずこんな感じでいまデータを持ってるとする

head(dat)

##   pre.fluency pre.accuracy pre.complexity post.fluency post.accuracy
## 1   11.916630     7.480500       4.644798     8.431656     13.070707
## 2   10.719751     7.090714       6.566487    13.652658     11.462576
## 3   11.980852     6.182161       5.731289    16.580042     11.287683
## 4   11.109400     7.256950       5.897070    16.966183      9.497555
## 5   10.488391     5.449349       2.673285    15.778588      9.994575
## 6    9.076144     8.621538       2.122296     9.465845      7.577611
##   post.complexity delayed.fluency delayed.accuracy delayed.complexity
## 1        8.639772       12.961946         7.740613           7.042510
## 2        6.397390        9.921947         7.520258           6.951231
## 3        6.917590       10.478401         9.378301           4.972710
## 4        5.897965        6.738074         8.983224           5.045168
## 5        6.214526       10.694636         8.510670           5.184044
## 6        7.644466        8.666421         8.863990           4.719461
##   subject
## 1       1
## 2       2
## 3       3
## 4       4
## 5       5
## 6       6

まずは縦横変換

参考までに以前tidyrの紹介を書いたブログ→変数が多い研究デザイン系の縦横変換(tidyr)

library(tidyr)
dat%>%
  gather(key=variable,value=score,-subject)%>%
    separate(col=variable,into=c('test','measure'))->dat2
head(dat2) #test, measure, scoreという列ができてるか確認

##   subject test measure     score
## 1       1  pre fluency 11.916630
## 2       2  pre fluency 10.719751
## 3       3  pre fluency 11.980852
## 4       4  pre fluency 11.109400
## 5       5  pre fluency 10.488391
## 6       6  pre fluency  9.076144

---

悩み

でもさーこうやってlong型にしちゃうとさー例えば事前のときの流暢さの指標の平均値とか出せないじゃん？wide型になってたら全部横にならんでるから列ごとに平均求めればいいだけじゃん？

long型になってても層ごとの記述統計出せちゃう

そう，dplyrならね

library(dplyr)
dat2%>%
    group_by(test)%>%
    group_by(measure,add=T)%>%
  summarise_each(funs(mean,sd,min,max),-subject) #summarise(mean=mean(score),sd=sd(score),min=min(score),max=max(score))でも同じ。もともとsummarise_eachは複数列に同じ関数を適用するときに威力を発揮するものだけど

## Source: local data frame [9 x 6]
## Groups: test [?]
## 
##      test    measure      mean        sd        min       max
##     (chr)      (chr)     (dbl)     (dbl)      (dbl)     (dbl)
## 1 delayed   accuracy  8.050833 1.0231177  5.5291358 10.250133
## 2 delayed complexity  6.043974 0.8633653  4.1740878  7.908483
## 3 delayed    fluency 11.966818 2.2327496  6.7231359 15.657528
## 4    post   accuracy  9.529745 2.4594148  2.7470626 16.876272
## 5    post complexity  7.211917 1.9710329  2.3591702 11.759346
## 6    post    fluency 14.624789 3.8821630  7.3521391 21.161635
## 7     pre   accuracy  6.941050 1.7114477  3.3338811 10.798454
## 8     pre complexity  4.326063 2.1271468 -0.7051241 10.391882
## 9     pre    fluency 10.239550 2.0043330  5.8044643 14.556087

※このデータの場合，test,measure,score,subjectの4つしかなく，グルーピング変数で2つ使っていて，-subjectするので必然的にscoreだけの計算になるけど，もっと列数多くて記述統計出したいのは1つだけの時は列を指定してあげる（下の例参照）

---

追記

• もしも人ごとに1試行が1行になっているようなlong型の場合は，最初に人ごとにgroup_byすればOK

#こんな感じ
#dat2%>%
#  group_by(subject)%>%
#  group_by(test)%>%
#  group_by(measure,add=T)%>%
#  summarise_each(funs(mean,sd,min,max),score)

• もし値に欠損があるときに関数が複数あるとちょっとめんどうだけど…

#こんな感じ
#dat2%>%
#  group_by(subject)%>%
#  group_by(test)%>%
#  group_by(measure,add=T)%>%
#  summarise_each(funs(mean(.,na.rm=T),sd(.,na.rm=T),min(.,na.rm=T),max(.,na.rm=T)),score)

---

ハッ，なに言ってるのこの人そんなのdplyrの威力をぜ～んぜんっ活かせてないんだからね！っていうのもあるかと思うのでもっと勉強したいと思いますがとりあえずtidyrからのdplyrを外国語教育研究のありそうなデータセットの形でやってみるとこんな感じになりますよというお話でした。

参考サイト

dplyrを使いこなす！基礎編

R dplyr, tidyr でのグルーピング/集約/変換処理まとめ

なにをゆう　たむらゆう。

おしまい。

20160516追記

よくよく考えてみると，そもそもCAFって指標ごとに分析するからaccuracy, fluency, complexityの3列で待ってる形じゃないと分析できなくない？

そんなときはgather->separateのあとにspreadしてあげる

dat2%>%
  spread(measure,score) ->dat3 #measureごとにscoreの列を作る
head(dat3)

##   subject    test  accuracy complexity   fluency
## 1       1 delayed  7.740613   7.042510 12.961946
## 2       1    post 13.070707   8.639772  8.431656
## 3       1     pre  7.480500   4.644798 11.916630
## 4       2 delayed  7.520258   6.951231  9.921947
## 5       2    post 11.462576   6.397390 13.652658
## 6       2     pre  7.090714   6.566487 10.719751

あとは上のやり方と同じ

dat3%>%
  group_by(test)%>%
  summarise_each(funs(mean,sd,min,max),-subject)

## Source: local data frame [3 x 13]
## 
##      test accuracy_mean complexity_mean fluency_mean accuracy_sd
##     (chr)         (dbl)           (dbl)        (dbl)       (dbl)
## 1 delayed      8.050833        6.043974     11.96682    1.023118
## 2    post      9.529745        7.211917     14.62479    2.459415
## 3     pre      6.941050        4.326063     10.23955    1.711448
## Variables not shown: complexity_sd (dbl), fluency_sd (dbl), accuracy_min
##   (dbl), complexity_min (dbl), fluency_min (dbl), accuracy_max (dbl),
##   complexity_max (dbl), fluency_max (dbl)

列数が多いので表示されていないけど，出力はデータフレームなので，出力結果を何かの変数にいれてあげれば良い

---

ちなみに，実をいうとwide型の段階で列ごとにapplyしたほうが(ry

library(psych)
apply(dat,2,describe)

## $pre.fluency
##   vars  n  mean sd median trimmed  mad min   max range skew kurtosis   se
## 1    1 50 10.24  2  10.45   10.24 2.21 5.8 14.56  8.75    0     -0.6 0.28
## 
## $pre.accuracy
##   vars  n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## 1    1 50 6.94 1.71   7.15    6.93 1.79 3.33 10.8  7.46 0.06    -0.43 0.24
## 
## $pre.complexity
##   vars  n mean   sd median trimmed  mad   min   max range skew kurtosis
## 1    1 50 4.33 2.13   4.09    4.33 2.07 -0.71 10.39  11.1 0.09     0.33
##    se
## 1 0.3
## 
## $post.fluency
##   vars  n  mean   sd median trimmed  mad  min   max range  skew kurtosis
## 1    1 50 14.62 3.88  14.78   14.67 5.19 7.35 21.16 13.81 -0.09    -1.27
##     se
## 1 0.55
## 
## $post.accuracy
##   vars  n mean   sd median trimmed  mad  min   max range skew kurtosis
## 1    1 50 9.53 2.46   9.48    9.44 2.69 2.75 16.88 14.13 0.24     0.67
##     se
## 1 0.35
## 
## $post.complexity
##   vars  n mean   sd median trimmed  mad  min   max range  skew kurtosis
## 1    1 50 7.21 1.97   6.93    7.22 2.18 2.36 11.76   9.4 -0.07    -0.31
##     se
## 1 0.28
## 
## $delayed.fluency
##   vars  n  mean   sd median trimmed  mad  min   max range  skew kurtosis
## 1    1 50 11.97 2.23  12.42   12.09 2.51 6.72 15.66  8.93 -0.42    -0.52
##     se
## 1 0.32
## 
## $delayed.accuracy
##   vars  n mean   sd median trimmed  mad  min   max range  skew kurtosis
## 1    1 50 8.05 1.02      8    8.06 1.06 5.53 10.25  4.72 -0.15     -0.2
##     se
## 1 0.14
## 
## $delayed.complexity
##   vars  n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## 1    1 50 6.04 0.86   5.97    6.05 1.01 4.17 7.91  3.73 0.02    -0.65 0.12
## 
## $subject
##   vars  n mean    sd median trimmed   mad min max range skew kurtosis   se
## 1    1 50 25.5 14.58   25.5    25.5 18.53   1  50    49    0    -1.27 2.06

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

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

補足

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パッケージの方が使い勝手が良さそうですね（ｱﾚﾚ

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

今日はこのへんで。

なにをゆう　たむらゆう

おしまい。

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

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

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

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

まずはパッケージの準備

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

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

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

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

dat %>%

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

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

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

dat2 %>%

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

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

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

おおおおー！！！！！

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

 

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

なにをゆう　たむらゆう

おしまい。

 

もう2ヶ月くらい前の話ですが，Nagoya.R #12で発表しました。

今更って感じなのですが，ちょっと今お手伝いでLMEをやっていて，自分でも理解があやふやな点がぼろぼろと出てきたのでメモ的に。

僕は発表中にlmer関数と，モデル比較に使えるstep関数の使い方を説明したのですが，どうやらstep関数はlmer関数で出力したもののみに対応している模様。というか，正規分布を仮定したモデルにしか適用できないみたいです。なので，lmerでfamily指定を使ってポアソン分布(poisson) や二項分布(binomial)を指定した場合には出力がされません。また，lmer関数でfamily指定すると警告メッセージでglmer関数を使うようにと言われます。なので，正規分布以外でやるときはglmer関数を使うほうがいいかもしれません。回帰の式の入力はほぼ一緒です。

ただし，glmer関数は結果の出力の解釈が実はちょっと難しくて，一発だけじゃ全水準の多重比較までみれないんですよね。なので，ダミー変数にいれたものの順番を入れ替えて計算を回していくようなのですが，それに使う引数がstartってやつっぽいのですよね。そんないちいちダミー変数入れ替えるなんてめちゃめんどくさいわけで，これで引数指定して一番最初にいれる水準を指定できるようになっているみたいなのです。ただしちょっと使い方がまだよく理解できていなくてですね…

実際に自分の研究でちゃんと使えるようになるにはここは避けて通れないわけなのですが，そっちばっかりに手を回しているわけにもいかず，RのヘルプやマニュアルとにらめっこしてはGoogle検索して…とかやってたらなんか1日終わっているみたいな幸せなのか不幸せなのかわからない昨日今日です。

いろいろ終わってません

なにをゆう　たむらゆう

おしまい。