タグ別アーカイブ: R

RStudio 1.0.143

MacのRStudioを新しいversion(1.0.143)にアップデートしました。

何が変わっているのかなと思ったら,Rmdでコードを実行した時に,もはやコードの実行結果がコンソール上に表示されなくなっていてびっくり。これはこれで便利なのかもしれないけど,慣れていなさすぎて逆に戸惑っています…

もちろん,コンソール上に直接打ち込んでいたり,R Script上で打ち込んでいれば結果はコンソールに出るのですが,RmdファイルのRチャンクをつかっている場合,そのチャンクの下に結果が表示されるようになっていて,大きいデータはGUIでポチポチするようになっています。

別ウィンドウに出せるようにもなっていて,Notebook Outputと呼ばれているようです。

多分出力の結果によってこのアウトプットのされ方も色々変わってくるのでしょうね。細かいことはまだよくわかってませんが,ちょっと使ってびっくりしたのはこの仕様でした。

個人的にはRStudioは好きで今まで使ってきていたのですが,これはなんかちょっと違うような…

なにをゆう たむらゆう

おしまい。

※追記

プロットも右下のペインではなくRチャンクの下に表示されるようになっているようです。

[R] データフレームの中の英文の語数を数える

以前,エクセルで英文の中の特定の単語までの語数を数えるという記事を書きました。今回はそれよりもっと単純で,Rのデータフレームの中の1部に英文が入っていて,その英文の語数を数えてデータフレームに付け足したいなっていう作業です。

%e3%83%96%e3%83%ad%e3%82%b0%e8%a8%98%e4%ba%8b%e7%94%a8_20170218

こんな感じでロング型のデータがdatという変数にあるとします。このデータフレームのなかの,sentenceという列には英文が入っていますが,それぞれ長さが違います。文が長いほうが反応時間も長くなっちゃいますから,そのために英文に含まれる語数を数えてそれも分析に使いましょうということです。僕がやっているような反応時間を扱う実験ではよくありがちな手続きです。本来なら刺激文作るときにそれもちゃんと数えておいて(エクセルで語数数えるのは簡単です。前掲の記事の一部がそれです),結果が出力される際にちゃんと語数の情報も含めておけばいいんですけど,まぁ今回はそれ忘れちゃったんでRの中でやってしまおうという事後的な対処法です。

さてさて前置きが長くなりました。Rでの文字列処理の基本は「分ける」・「数える」だと思います。僕はそんなにRで文字列処理しないんですが,今回の場合で言うと,データフレームの要素に含まれる英文をスペースで区切り,区切ったときにいくつかに分けられた要素の数を数えることで語数をカウントするという手続きになります。

今回使う関数は次の4つです。

  • (as.character関数)
  • strsplit関数
  • length関数
  • sapply関数

順番に説明していきます。strsplit関数は

strsplit(分けたいもの,区切り文字)

という感じで使います。今回の場合,分けたいのはdatの中のsentenceという列(dat$sentence)で,スペースで区切りますから,

strsplit(dat$sentence, ” “)

とすれば,

%e3%83%96%e3%83%ad%e3%82%b0%e8%a8%98%e4%ba%8b%e7%94%a8_20170218_2

こんな感じで単語ごとに区切ってくれます。この出力は1文ずつのリスト形式になっていて,そのリストの要素が英文を構成する単語になっています。もしも,このstrsplit関数がうまくいかなければ,

is.character(dat$sentence)

として,データフレーム中の英文が文字列形式で格納されているかどうか確かめてみてください。因子型(factorial)になっていると,strsplit関数は適用されません。これを変えるのが,as.character関数です。

strsplit(as.character(dat$sentence))

または,

dat$sentence<-as.character(dat$sentence)
strsplit(dat$sentence, ” “)

としてみてください。

さて,これでひとまず「区切る」作業はできましたが,まだ「数える」にはいたっていませんよね。数えるためによく使われるのが要素の数を数えるlength関数です。ただし,例えば

length(strsplit(as.character(dat$sentence)))

としてもうまくいきません。これでは出力されたリストの要素を数えるので,元のデータフレームの行数が表示されるはずです。ここでは,「リストのそれぞれの要素に対してlength関数を適用したい」わけなので,もうひと工夫必要です。そこで使うのがapply系列のsapply関数です。この関数は,

sapply(リスト,関数)

という形で引数を与えてあげると,リストの各要素に関数を適用してくれます。これと同じことができるのがlapply関数です。lapply関数を使うと,結果の出力もリスト形式で返されます。ただ,ここでは語数が並んだベクトル形式のものを手に入れて,それをもとのデータフレーム(dat)に追加したいわけなので,sapply関数のほうがベターです。lapply関数にすると,そのあとリストをベクトルに直すunlist関数を使う必要があって二度手間になってしまうので。というわけで,

sapply(strsplit(as.character(dat$sentence), ” “), length)

こうしてあげれば「スペースで区切る-> 語数を数える -> ベクトル形式で結果をゲット」という作業が完了します。あとは,

dat$num.words<-sapply(strsplit(as.character(dat$sentence), ” “), length)

としてあげれば,datに新しくnum.wordsという列が追加され,そこにはsentenceの列に入っている英文の語数が格納されることになります。めでたしめでたし。といきたいところなのですが,最後におまけでもうひとつ。

上の書き方だと,いくつもの関数が入れ子状態になっていて,一見して作業がわかりにくいですよね?そこでパイプ演算子を使って同じことをもう少しわかりやすく書き直してあげます。

dat$sentence%>% #datの中のsentenceという列を取ってくる
as.character%>% #文字列型に変換
strsplit(. , ” “)%>% #スペースで区切る
sapply(., length) -> dat$num.words #リストの要素ごとに語数を数えてdatというデータフレームのnum.wordsという列にいれる

こうやって書くと,作業の過程が見えてわかりやすいですよね。ちなみにこのパイプ演算子はパッケージ依存です。dplyrやtidyrを使うときに使うやつですね。たぶんこれらのパッケージを読み込んでいればパイプ演算子も使えるようになっているはずです。もともとはmagrittrというパッケージに入っているのがパイプ演算子が%>%なんですよね(たしか)?

というわけでlapplyの出力をunlistしてできたできたと思っていたらsapply関数の存在に気づいてなんだよ二度手間じゃんかよって思ったのでブログ記事にしました。

なにをゆう たむらゆう。

おしまい。

 

 

 

[R] 欠損値を1つでも含む行を抽出

ものすごい初歩的なことだったのに,調べてもなかなかヒットせずに格闘したのでメモ。

RでNAを処理する方法はたくさんあって,ある特定の列にNAが含まれる行だけ引っ張ってくるなどは結構簡単だったりする

subset(dat, is.na(dat$A))

とかでいい。ただし,データフレームにある複数の変数の中で1つでもNAが入っているかどうかとなるとちょっとどうしたらいいんだこれは,となる。

na.omit(dat)としてNAが含まれている行を全削除すればいいのでは?となるかもしれないが,実は今扱っているデータがlong型のうえにかなり複雑なデータフレームになっていて,同じidが複数行にまたがってしまっているという状態。実験のデザインがかなり複雑なのでこれはこれで仕方がない。

つまり,複数行のうち1行でもNAがあるidは,NAが含まれない行でも削除しないといけない。そのためには,NAが含まれる行のidを特定する必要があるわけだ。そこで,欠損値を含む行を抽出し,その中からid列だけを取ってくるという手続きを取る必要がでてきた。こうすれば,あとはdplyrのfilter関数なり,subset関数なりでそのidが含まれる行をすべて抜いて新しいデータセットを作ればよい。

それで,この,NA行の特定をベクトル形式で返してくれる関数が,complete.casesという関数。is.na関数にデータフレームを渡すと,その結果は行列(matrix)形式になる。これはこれで,どの行のどの列にNAがあるかを確かめるのにいいのだけれど,こちらが欲しいのは,「何行目」にNAが含まれるのかという情報だけ。そこで,complete.cases関数を以下のような感じで使う。

subset(dat, complete.cases(dat)==F)

あるいは,

dat%>%
filter(complete.cases(dat)==F)%>%
select(id)

上の方法だと,結果はデータフレームだけど,下の方法だとid列を引っ張ってくるところまでやるので,これをなにか別の変数(ここではmissIDにしてみる)にいれてあげて,

subset(dat,id != missID)

とかやれば,id列でmissIDに該当しない行だけを抽出してくれる。下書き状態で放置されていたからとりあえず最後まで書いたけど,自分がなんのデータに対してこの手続をしていたのかはもう思い出せないw

データ分析業務ばかりやってるのは楽しいからいいんですけどね。

なにをゆう たむらゆう。

おしまい。

※追記(2017.03.11)

subset(dat,id != missID)

とやってもうまくいかないことがあるようです,というか多分うまくいきません(自分でやってみてダメだったので)。多分ですが原因の1つはmissIDがデータフレーム形式になっている場合。ベクトルになおしてください。もう1つの原因は,subsetなりdplyrのfilterなりで複数の条件を設定する場合に,参照先がベクトル形式の場合には通常の論理演算記号が機能しないことです。多分こちらのほうが大問題。時間が経ってから書いたからかこのことをすっかり忘れていましたすみません。正しくは%in%を使って,

subset(dat,id !ID %in% missID)

とします。この場合,データフレームの中の参照先の列名の前に!をつけます。これをつけないと「当てはまるもの」を取り出すことになります。

dplyrを使ってみる:long型から記述統計

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

// <![CDATA[
window.onload = function() {
var imgs = document.getElementsByTagName('img'), i, img;
for (i = 0; i

 

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

#事前・事後・遅延事後で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

glmerからのlsmeans

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

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

Error in R lsmeans() function: Error in format.default(nm[j], width = nchar(m[1, j]), just = “left”)

この掲示板の書き込みを見てみると,どうやらRのバージョンが古い(R3.2.0以前)とこのエラーメッセージが表示されるらしいです。なんか根本的な解決にはなってないような気もしますが…

とはいえ私もR3.2.0ユーザーだったので,これを気に新しいRをインストールすることに。R.3.3.0が最新らしいです。最新のに変えると今まで使ってたパッケージに不具合が出たりすることがあって不安でしたが,実際にRのバージョンを新しくすればglmerのあとに多重比較ができました。

fit1<-glmer(RT ~presentation*condition+ (1|subject) + (1|item), family = Gamma (link= “identity”), data=dat,glmerControl (optimizer=”bobyqa”)

lsmeans(fit1, pairwise~condition|presentation)

こんな感じ。

要因の順番は,多重比較したい要因が先にきます。つまり,上の例だと,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の基礎がわかってなにこれめちゃ便利って感動したのでそのこともまた書きます。

なにをゆう たむらゆう。

おしまい。