R による整列ランク変換を使った分散分析
整列ランク変換（Aligned Rank Transform：ART）という言葉を聞いたことありますか？ 私はなかったです。ART（えーあーるーてぃー。「あーと」ではない）について、Tokyo.R 第61回で LT をしてきました。補足も兼ねて ART についてまとめます。
きっかけ
先日、自分の専門分野の先行研究を調べていて論文（Heidenreich, et al. 2015）を読んでいたら次のような記述に出会いました。
“As our data were not normally distributed (D(243)=2.193, p<0.01), we employed the aligned rank transform (ART) procedure for hypothesis testing. For non-normal distributed data, the ART procedure is more robust and powerful than the traditional analysis of variance (ANOVA)”
超訳すると「データが正規分布に従わなかったので、aligned rank transform (ART) を使って検定やった。非正規分布のデータなら、ARTを使ったほうが普通の分散分析よりロバストでパワフルだ」。
ART について、ネット上で検索してもほとんど日本語の情報がありません。かろうじて「整列ランク変換」という訳語が見つかるくらいです。そこで少し調べてみました。
ART法はタイプIエラーをコントロールできるエクセレントな方法だ。研究者がやりがちな古典的分散分析と比べてパワーがある（Mansouri, et al. 2004）
ART はどんな分布のデータにも使えるわけじゃないけど、ロバストだし一般的なノンパラメトリック検定より望ましい性質を持っている（Higgins, et al. 1990）
ART 検定はパワフルでロバストで使いやすい（Mansouri 1998）
交互作用を調べるとき、いくつかの仮定が成り立たない場合、パラメトリックな分析は必ずしもよいわけではない。代替の策として、調整ランク変換検定（adjusted rank transform test）がある。それはパラメトリックな手法とノンパラメトリックな手法の中間みたいなものだ（Leys & Schumann 2010）
つまり、分散分析したい、特に交互作用に注目したい時にデータが正規分布に従わない場合、従来の分散分析は（厳密には）使えないので、整列ランク変換を施してから分散分析すればいいよ、ということだと思います。
R パッケージがないか探したところ、ARTool というものがちゃんとありました。CRAN に登録されています。使い方もカンタン。以下、ART を使った分散分析の方法です。例題として、性別と世代の不倫の許容度に対する違い・交互作用を見たいとします。応答変数は「結婚していても場合によっては不倫を許せる？（１：許せない～４：許せる）」です。4件法でスミマセン。
０．データの分布の確認
ggplot(rensyu, aes(furin)) + geom_histogram(fill = "red", alpha = 0.7, binwidth = 0.5)
あきらかに正規分布には従わないっぽい。もちろん、正規性の検定はすべきですが、ここでは飛ばします。
１．ARTool を使えるようにする。
install.packages("ARTool") # 1回のみでOK。 library(ARTool)
２．データを整列ランク変換する。
henkan = art(furin ~ nendai*sex, data = rensyu) summary(henkan)
art() で変換しますが、式はよくあるタイプです。式は「art(furin ~ nendai + sex + nendai:sex, data = rensyu)」を短縮した書き方にしています。どちらでもOK。lm() のときと同じです。
変換後、summary(henkan) の出力結果に注意します。すべての数値が０になっていることが望ましいです。今回は残念ながらそうならなかったので、警告が出ています。しかし、練習なので先に進めます。
## Warning in summary.art(henkan): F values of ANOVAs on aligned responses not ## of interest are not all ~0. ART may not be appropriate. ## Aligned Rank Transform of Factorial Model ## ## Call: ## art(formula = furin ~ nendai * sex, data = rensyu) ## ## Column sums of aligned responses (should all be ~0): ## nendai sex nendai:sex ## 0 0 0 ## ## F values of ANOVAs on aligned responses not of interest (should all be ~0): ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.0000000 0.0000000 0.0000000 0.0144000 0.0007531 0.0854000
３．分散分析を実行する。
anova(henkan) ## Analysis of Variance of Aligned Rank Transformed Data ## ## Table Type: Anova Table (Type III tests) ## Model: No Repeated Measures (lm) ## Response: art(furin) ## ## Df Df.res F value Pr(>F) Sum Sq Sum Sq.res ## 1 nendai 4 1261 3.7856 0.0045704 ** 1935456 161178877 ## 2 sex 1 1261 100.7097 < 2.22e-16 *** 10944212 137033949 ## 3 nendai:sex 4 1261 4.6317 0.0010275 ** 2369093 161247964 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
「anova(変換したデータ)」でOK。実際には、anova.art() が動いて、デフォルトでタイプ３の分散分析をやってくれます。世代と性別の主効果、交互作用は有意になりました（なお、効果量 は世代がほぼ0、性別が0.02、交互作用が0.01でほとんどインパクトはありません）。いちおう、交互作用のグラフをプロットします。
library(ggplot2)
ggplot(rensyu) +
aes(x = nendai, color = sex, group = sex, y = furin) +
stat_summary(fun.y = mean, geom = "point", size = 2) +
stat_summary(fun.y = mean, geom = "line", size = 1) +
xlab("世代") + ylab("不倫の許容度") +
scale_colour_discrete(name ="性別")
伝統的な分散分析と結果を比べてみる。
aaa = lm(furin ~ nendai*sex, data = rensyu) library(car) Anova(aaa, type="III") ## Anova Table (Type III tests) ## ## Response: furin ## Sum Sq Df F value Pr(>F) ## (Intercept) 224.03 1 535.4585 < 2e-16 *** ## nendai 2.39 4 1.4269 0.22278 ## sex 2.22 1 5.3061 0.02141 * ## nendai:sex 3.02 4 1.8051 0.12541 ## Residuals 527.59 1261 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
同じくタイプ３でやってみましたが、こちらでは世代の主効果と交互作用が有意になりませんでした。いい感じに対比できましたが、データが怪しいので、結果はあくまで参考用として受け取ってください。
まとめ
整列ランク変換を利用した分散分析では交互作用が有意になった。
ただし、今回のデータはちゃんと整列ランク変換できてはいないのでなんとも言えない。
今後、データが正規分布に従わない場合、使ったほうがよい手法かもしれない。
さいごに
多重比較もできます。繰り返しのデータも分析可能のようです。詳しくはパッケージ作者の Github ページを参考にしてください。今後、必要あれば、使っていきたい分析法です。
参考文献
ARTool パッケージ
https://github.com/mjskay/ARToolFeys, J. (2016). Nonparametric Tests for the Interaction in Two-way Factorial Designs Using R. The R Journal, 8(1), 367-378.
Heidenreich, S., Wittkowski, K., Handrich, M., & Falk, T. (2015). The dark side of customer co-creation: exploring the consequences of failed co-created services. Journal of the Academy of Marketing Science, 43(3), 279-296.
Higgins, J. J., Blair, R. C., & Tashtoush, S. (1990). The aligned rank transform procedure.
Leys, C., & Schumann, S. (2010). A nonparametric method to analyze interactions: The adjusted rank transform test. Journal of Experimental Social Psychology, 46(4), 684-688.
Mansouri, H. (1998). Multifactor analysis of variance based on the aligned rank transform technique. Computational statistics & data analysis, 29(2), 177-189.
Mansouri, H., Paige, R. L., & Surles, J. G. (2004). Aligned rank transform techniques for analysis of variance and multiple comparisons. Communications in Statistics-Theory and Methods, 33(9), 2217-2232.
Wobbrock, J.O., Findlater, L., Gergle, D. and Higgins, J.J. (2011). The Aligned Rank Transform for nonparametric factorial analyses using only ANOVA procedures. Proceedings of the ACM Conference on Human Factors in Computing Systems (CHI ‘11). Vancouver, British Columbia (May 7-12, 2011). New York: ACM Press, 143-146.
目を通してないけど、重要そうな文献
Salter, K. C., & Fawcett, R. F. (1993). The ART test of interaction: a robust and powerful rank test of interaction in factorial models. Communications in Statistics-Simulation and Computation, 22(1), 137-153.
Sawilowsky, S. S. (1990). Nonparametric tests of interaction in experimental design. Review of Educational Research, 60(1), 91-126.
RMarkdown で参考文献を作るときにちょっとつまずいた話
Rを使うドキュメントはできるだけ RStudio内 で完結するように（コピペが面倒なのと，スライドも配布するPDF資料もRStudioでできるので）環境を整えているところですが，参考文献も自動で挿入できるようにしたいと思いました。やってみて，つまずいたところをメモします。
RMarkdown チートシートを見ると，2ページ目に「Citations and Bibliographies」というセクションがあります。YAML ヘッダーに簡単な追記をすればできるようなので，やって見たところ，次のように怒られました。
pandoc: Filter pandoc-citeproc not found
エラーメッセージをそのまま検索すると，こちらの答えがありました。
どうやら，
brew install pandoc pandoc-citeproc
でいいそうです（Mac で homebrew が使える場合）。確かに，pandoc-citeproc をインストールした後は，YAMLヘッダーに参考文献のリスト（.bib）を用意すれば，うまく参考文献が表示されるようになります。
今回，日本語文献だったので，順番がうまくいきませんでしたが，論文ではないので，まぁOKかなと思います。記述したYAMLヘッダーはこちら。
--- title: "Rの操作の復習と簡単な分散分析" author: Hikaru GOTO date: 2017年4月15日（土） output: rmdshower::shower_presentation: self_contained: false theme: ribbon ratio: 4x3 bibliography: mylib.bib csl: journal-of-business-research.csl ---
なお，.bib ファイルの作成はこちらを参考にしました。.csl は zotero のリポジトリからダウンロードしています。JBR は私が比較的よく見る論文誌です。
ちなみに，RMarkdown で参考文献を作るときの公式ガイドはこちらです。RMarkdown で引用するとき方法と参考文献一覧を表示させるは公式ガイドを参照してください。
P.S.1
RMarkdown チートシートは日本語訳もあります。超便利！ 感謝。
https://www.rstudio.com/wp-content/uploads/2016/11/Rmarkdown-cheatsheet-2.0_ja.pdf
P.S.2
zotero のスタイルを使った場合は，ネットに接続してないと knit が失敗します。
ggdendro と PAC 分析を使ったサービス学会の発表
2017年3月28日（火）、広島県情報プラザというところで開催されたサービス学会第５回国内大会に参加してきました。仕事以外で広島に滞在したのは小学校の修学旅行ぶりです。
http://ja.serviceology.org/events/domestic2017.html
今回は研究をスタートし始めたばかりのネタで、調査・分析が進んでおらず、いろいろと参加者と相談したいという理由で、ポスター発表にしました。
PAC 分析というものを用いたのですが、この分析法ではデンドログラムを使います。せっかくなので、ggdendro を使ってみようと直前までポスターをいじっていました。PAC分析の手法についてはこちらを参照してください（ネット上にもたくさん情報はあります）。
- 作者: 内藤哲雄
- 出版社/メーカー: ナカニシヤ出版
- 発売日: 2003/04
- メディア: 単行本
- クリック: 10回
- この商品を含むブログ (3件) を見る
R で普通にクラスター分析するとこんな感じ。
ポスターを作るときはこれをパワーポイントに貼り付けて、回転させてトリミングして整えていきますが、ちょっと色合いがさみしい感じ。私のポスターは茶系で地味なので、ちょっと色がほしかったのです（ggplot のあのグレーが好き）。
というわけで、ggdendro を使うことにしました。何もしないと ggdendro は背景色が白です。
ここから ggplot2 に持って行き、縦軸と横軸を消していきます。例によって、kazutan 先生の資料を参考にしました。
https://rpubs.com/kazutan/ggdendro_test
できあがりはこちら。ggplot っぽくなりました。
1行目で ggplot に渡すデータを作り、5行目でデンドログラムを横倒しにして、6～8行目で縦軸と横軸を削除しています。あとはこの図をパワーポイントに持っていき、情報を追加したら完成。
いつも直前にホテルでかんづめして作ってます。
ちなみに、A０サイズで7000円ほどかかりました。高いけど、直前まで出力を受け付けてくれたキンコーズ広島本通店さん、ありがとう。
R Markdown スライドとPDF 文章（日本語含む）作成メモ
月1回のゼミ前にRを使った統計分析の勉強会をすることになりました。せっかくなので（これまで何度か挫折した）R Markdown でやろうと思い立ち，今回は成功したので資料づくり（スライドとPDF）のメモを書いておきます。自分が忘れないように。
1）R Markdown でスライドを作ろうと思い，以前，本家のページを読んだ。
rmarkdown.rstudio.com
2）ただ，やはり日本語の情報のほうがストレスが少ないので，こちらを熟読した。 わかりやすくまとまっていて個人的大感謝祭。
https://kazutan.github.io/kazutanR/Rmd_intro.html
3）YAML はコピペがよい。複雑なものは作らないので，達人たちの力を借りる。「revealjs::revealjs_presentation:」の最後のコロンに気づかず，入力を忘れていて，「どーしてもエラーが出る」と1時間くらい悩んでしまった反省から。
4）作ったスライドをPDFにする方法はこちらにある。PDFにすれば Speakerdeck とかにも簡単にアップできて便利。
https://kazutan.github.io/SappoRoR6/rmd_slide.html#/pdf
5）作った資料（日本語含む）をPDF文書にするときの YAML ヘッダーの書き方はこちらにある。xelatex の YAML ヘッダーを使う場合は，IPAフォントを入れておかないとエラーが出た。
RStudio で Rmd ファイルから日本語 PDF を作成する方法 | miyazakikenji
以上です。
1〜2年前は「エラー出た。解決策がない（日本語で）」と困っていたのですが，ここ最近，有志の方が情報を出してくれて，本当に便利になりました。これでまた1つ，パワーポイントに結果の図とか数値を貼り付ける作業から解放されそうです。
旅人と RStudio
このエントリーは RStudio Advent Calendar 2016 の23日目のエントリーです。ちなみ，本日は Osaka.stan 第2回に参加してきます。
さて，私は仕事柄，出張が多く，2016年は12月末までに224泊する予定です。だいたい1日ごとに移動しています。たぶん，日本でいちばんビジネスホテルで RStudio を使っている人間だと思います（なんの自慢にもならないですけどね！）。
せっかく出張が多いので，全国の名所でRStudioを使っている様子の写真を撮ってアップしようと企画したのですが，なんだか恥ずかしいのでやめました。
というわけで，本エントリーでは R や RStudio を付き合い始めてだいたい2年。日頃，とても勉強させていただいているRコミュニティーへの感謝も込めて，自分なりの学びや使い方を振り返ってみたいと思います。対象は例によって，統計分析やRを独学で学ばないといけなくて，ある程度時間よりもお金を少しかけられる人（つまり社会人）向けです。
パイプ演算子の出し方
RStudio のショートカットで %>% を出すには Mac の場合，シフト＋コマンド＋M です。でも，どーしてもこのショートカットが苦手なんです。いちいちキーボードを見ないと打てないし，ホームポジションから離れてしまう……。
で，現在は，Alfred Powerpac のスニペット機能を使っています。textexpander でもいいのでしょうが，月額課金よりは買い切り型でいろいろできる Alfred Powerpac を選びました。35英ポンド。
↓こんな感じでセミコロン2つで %>% が出るように設定しています。
他にも，よく読み込むパッケージを登録していたりします*1。
R と RStudio に慣れるには
身近に詳しい人がいればいいのでしょうが，私の場合はいません。コツコツ学べる才能があればいいのでしょうが，私にはそんなにありません。やらざるを得ない環境というか，追い詰められないとできないタイプです。ある程度コストをかけると「やらねば！」という気持ちになるので，そうしています。私がけっこう勉強になったなぁと思うものを2つ挙げます。
1．DataCamp
無料でいくつか受講できます。例えば「Introduction to R」。最初に動画で簡単な説明があって，RStudio っぽい画面で操作をしていきます。英語ですが，だいたい雰囲気はつかめます。私は昨年，3ヶ月ほど課金していくつか受講しました。
1ヶ月30ドルくらい。得るものと比べたら安いくらいです。3ヶ月やれば，だいたい「これできるようになりたいな〜」が受講しきれます。難しいのはやらないし。でも最近，またちゃんとやり直そうと思っています。
2．『RStudioではじめるRプログラミング入門』
簡単な例からスタートして，少しずつ難しくなっていきます。私が Rの 本で最後まで演習をしながら読み通せたのはこの本でした。「プログラミングとかよくわかんない」という私のようなものにはちょうどよく，勉強になりました。ただ，まったく未経験者にはいきなりこの本は難しいと思うので『はじめてのR』をざっと目を通してからのほうが無難です。
- 作者: Garrett Grolemund,大橋真也,長尾高弘
- 出版社/メーカー: オライリージャパン
- 発売日: 2015/03/25
- メディア: 大型本
- この商品を含むブログを見る
なお，いちいち本なんか読んでられないという方には近い内容のものがオンライン講座になって発売されています。こちらです↓
およそ160ドルですが，私の場合，たまたま何かのセールのときに99ドルでした。
英語はちょっと厳しいという場合は，Udemy がいいです。これもセールだったので登録しました（見てない）。いまもセールで1800円のようです。お得ですね。
上記のように，課題的なものがあるとか，ステップバイステップのものが私にはあっているようでした。
もし RStudio に飽きてしまったら…
最近はいろいろと選択肢があります。R のアドカレ15日目のように Atom内で動かせる Hydrogen もありますし，Jupyter notebook とか，nteract とか，Sublimetext 内で動かしたりとか。
最近はもっぱら Beaker Notebook にしています。Jupyter と違って，初めから等幅フォントで分析結果が見やすかったし，なんだか雰囲気がいい。でも，ちょっと起動が遅い。今日，Osaka.stan で Beaker Notebook を使っているグレーのスーツを着ている人がいたら，それはたぶん私です。
それでも RStudio を手放せない理由
実は個人的に RStudio でいちばん気に入っている機能はパッケージ管理です。新しく使うパッケージのインストールも更新も楽チン。しばらく放置していて，たくさん更新があると，一気にアップデートができるので気持ちがいい。かゆい所に手が届くというか，そういうのがRStudio のいいところだと私は思います。
以上，RStudio Advent Calendar 23日目でした。
*1:.Rprofile に書いておけ，という意見もあるでしょうが，毎回これらのパッケージを使うわけではないので，必要な時に読み込んでいます。
ミーハー.R
本エントリーは R Advent Calendar 2016 11日目のものです。
分析ネタではありません。日曜日なので箸休め的なエントリーにしようと思います。
さて，所属している社会人大学院のゼミに2ヶ月ほど前，新しい方が入ってきました。月1回のペースでゼミが始まる1時間ほど前に統計勉強会を有志と行っているのですが，初めてのその方が何気な〜く「Rって敷居高いですよね〜」と発したとき，私は頭が真っ白になってしまいました。「ええっ！？ そんなことはない！！」と。
その後もモンモンとして，今に至ります。ひょっとして自分はRに浸りきっていて，もはや抜け出せない状態になってしまっているのかもしない。それは危険です（なんとなくそう思う）。せっかくの機会なので，その気持ちを整理したいと思います。
私が R に触れたキッカケ
ちゃんと統計分析を学んだのは専門職大学院に行ってからですが，そこで初めてアルファベット4文字の超有名分析ソフトに触れました。学籍があれば1年間1000円でのフル機能が使えたのです。課題ごとにレポートを提出していたのですが，そのソフトから出力されたグラフをワードにぺたっと貼ったときに，強烈にダサいと感じてなえたことがあります。
（いま冷静に見たら，そうでもないかな）
きっとカッコよくする方法があるのだと思います。「いや，研究は中身が大事だ。外見は二の次」という意見もわかるのですが，でも，なんかもうダメ。ぜひカッコよくしてほしいです･･････。
とはいえ，修論を書くときに手法として量的分析を選んだ*1ので，しばらく使っていました。修論ではSEMを使うことにしたのですが，追加料金なしでAMOSも使えたため，続けて使っていたのです。
そこでもう1つ衝撃的なものを見ました。こちらです。
トラックですよ。トラック。大人のアソビゴゴロってやつですかね。ユニークかもしれないけど，まったくカッコよくない。
当時，心の余裕がなかった私はこれを理解できず，さらには分析がうまく進まないため（40個くらい変数があって，それをチマチマ動かして分析するのはたいへん），何か方法があるのではとRに走りました。結局，Rでもうまく行かなくて最終的にはMplusを使って修論を書き上げたのですが，これがRに触れるキッカケでした。
大胆な仮説
ここまで書いて気づきました。ひょっとして，ミーハーなだけなのでは，と。そこで大胆な仮説を立てみます。「Rユーザーはミーハーである」。
確かめるべく，5名に「こういう画面ってカッコよくないですか？」と投げかけたところ，1名だけが即答でイエスと答えてくれました（某Y社の偉い方）。こんな画面をカチャカチャいじっているのはカッコいいと思うんですよ。
イメージ図
でも残念ながら，この仮説は採択されないようです。きっとこれはここで「イエス」と答えたら自分がミーハーであることを公言したような物ですから，みんな素直に答えてくれなかったんだと主言っています。だから，私の心の中では有意差があります。今後も調査を継続していきたいと思います。
これから tidyverse に慣れるなら
最後に少しだけまともな情報提供を。「Rの情報は基本，英語」というのが敷居が高いと思われるもっとも大きな理由の1つだと思います。おそらく，Rにそれほど抵抗がない人は英語がそれなりにできる（少なくとも読める）のではないでしょうか。アドベントカレンダーを見てる方なら，きっとそう。
というわけで，12月23日に “R for Data Science”という本が発売されます（国外で）。もとの原稿（？）が公開されていて，ちょくちょく読んでいるのですが，初級者（自称を除く）から脱却できそうレベルな人にちょうどよいと思います。purrrのやさしい解説があります。いっしょに勉強していきましょう！
*1:「質的分析より先に量的分析に慣れておくべき。その逆はたいへんだから」という先生たちの教えより。
rstanで分析をするのに必要なR力（りょく）
Stan Advent Calendar 2016 の 4日目です。非マニア枠です。
本日のエントリーは Stan を普通に使える人、情報科学・データ科学な人を対象にしていません。少し R に慣れてきて、ベイズ統計モデリングにも興味を持ち始めた人向け。ずばりターゲットは次のようなレベルの人です。
いざ R を使って分析するときには使い方を忘れていて検索しながら使う。
dplyr の使い方に慣れておらず、めんどくさくなって結果エクセルを使ってしまう。
つまり私。そんな中途半端なレベルな人（私）向けに『Stan と R でベイズ統計モデリング』（以下，アヒル本とします*1）を読んでいるときにつまずくであろうところを書いておきます。アヒル本を読む・理解するには、適度に R と ggplot2 を使えるスキルが必要だからです。
さて，R を使った統計分析の本をいくつか読んで、簡単な分析をできるようになったという人が Stan を使ったベイズ統計モデリングをするときにつまずくのは次のようなところだと思います。
- 結果が出るまで待ち時間が長い。
- リスト形式でデータを渡すのを戸惑う。
- 結果が入ったオブジェクトを表示させたらズラーッと出てきてしまい、見たい結果を見るために必死に上へスクロールする。
- タイプミスする。
- 結果からほしいデータを取り出せない。
1. 待ち時間が長い
私の場合，簡単なモデルでも30秒くらいかかります。
lm()と違って，パッと結果が帰ってこないので，学習のリズムが狂いますよね。これは仕方ない。オススメは BGM をかけることです。ガンダム UC で流れていた MOBILE SUIT です。無事ちゃんと結果が出るか，ドキドキ感が増します。
2．リスト形式にとまどう
真の R の初級者はあまりリスト形式に触れることはありません。だいたいデータフレーム形式で間に合うので。lm(y ~ x, data=d)な感じで。しかし，Stan を使う場合，
> data = list(N=row(d), A=d$A, Score=d$Score/200, M=d$M)
というように分析に使うデータを別に作らなければなりません。「リスト形式ってなんだよ、めんどくさい」となるんですよ*2。
リスト形式とは、Todo リストみたいなものだと思ってください*3。Todo リストには、いろんな種類、内容を書けます。R のリスト形式も同じ。いろいろなデータをヒトマトメにできるのです。上の list だと N には1つの数値しか入っていませんが、A や Score にはたくさんの数値が入っています。そういうのをまとめて1つのものにするのがリスト。
Stan を使うには専用のデータを作ってあげる必要があって、それはリスト形式なんだと納得しましょう。慣れるしかありません。
3．結果の表示が多過ぎ！
（RStudioの場合）「サンプリングもうまくいった。よし，結果見てみよう」
> fit
どばっと結果が出ます。そしてパラメータを確認したいので上へ必死にスクロール。
こういうのを10回くらいやっていて「ひょっとしたらいい方法あるんじゃないか」と思いました。ありました。たとえば、b1・b2 というパラメータを見たい場合は次のようにします*4。
> print(fit, par="b") Inference for Stan model: model5-6b. 4 chains, each with iter=2000; warmup=1000; thin=1; post-warmup draws per chain=1000, total post-warmup draws=4000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat b[1] 3.57 0 0.10 3.39 3.51 3.57 3.63 3.76 976 1.01 b[2] 0.26 0 0.04 0.19 0.24 0.26 0.29 0.34 1820 1.00 b[3] 0.30 0 0.15 0.01 0.20 0.30 0.40 0.59 967 1.01 Samples were drawn using NUTS(diag_e) at Sat Dec 3 10:54:04 2016. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).
事後対数尤度も見たいというなら、
> print(fit, par=c("b","lp__"))
便利ですね。Rhat を小数点とか細かく見てみたいなという場合は例えばprint(fit, par="b", digits=3)とすればOKです。詳しい人は「当たり前」と思ってしまうレベルですが、真の初級者にはこういう tips を知っておくと，途中でめんどくさくなってやめるリスクが減ります。
4．タイプミス
私の場合、1つの Stan ファイルを実行するのに4～5回くらいエラーが出ます。本当にささいなミスで。英単語のつづりが間違っていたとか、閉じカッコが足りなかったとか、引用符がなかったとか*5。これの解決法はわかりません。誰か教えてほしいです。print を pritn とやっちゃうくらいなので。
5．ほしいデータを取り出せない。
アヒル本にあるようにrstan::extract(fit)で結果を取り出せます。取り出せますが，そのデータをいじるのには慣れが必要です。「fit の中身ってどうなっているんだろうと」と思ってstr()を使ってもよくわからない。なんか表示結果に圧倒されてしまう･･･。
> ms = rstan::extract(fit) > str(ms) List of 4 $ b : num [1:4000, 1:3] 3.48 3.57 3.56 3.52 3.62 ... ..- attr(*, "dimnames")=List of 2 .. ..$ iterations: NULL .. ..$ : NULL $ lambda: num [1:4000, 1:50] 3.64 3.69 3.67 3.66 3.71 ... ..- attr(*, "dimnames")=List of 2 .. ..$ iterations: NULL .. ..$ : NULL $ m_pred: num [1:4000, 1:50] 34 40 54 32 43 32 32 40 36 40 ... ..- attr(*, "dimnames")=List of 2 .. ..$ iterations: NULL .. ..$ : NULL $ lp__ : num [1:4000(1d)] 6897 6896 6898 6898 6897 ... ..- attr(*, "dimnames")=List of 1 .. ..$ iterations: NULL
数字とアルファベットがたくさんあって，イラッときますね。でも落ち着いて。普通にfit = stan(～～～)としたときは合計4000個のMCMCサンプルが入っています。$b のところを見ましょう。推定したパラメータ b はb1, b2, b3の3つある（「1：3」のところ）。4000行3列にまとまっている状態です。なので，b1 のヒストグラムはhist(ms$b1[,1])で書けます。最大事後(MAP)推定値も max(ms$b1[,1]) で出ます。
ggplot を使うときはちょっとやっかいです。そもそもms$b[,1]はデータフレーム形式ではないから、そのまま ggplot には使えません。いったんデータフレーム形式に変換します。data.frame()のカッコの中が大事。データフレームを作ったときの列名をつけておきます。
> b1 = data.frame(b1 = ms$b[,1]) > str(b1) 'data.frame': 4000 obs. of 1 variable: $ b1: num 3.48 3.57 3.56 3.52 3.62 ... > head(b1) b1 1 3.484667 2 3.572416 3 3.559215 4 3.521891 5 3.615658 6 3.657821
b1という列名がついた4000行1列のデータフレームになりました。ggplot でヒストグラムを作るなら
> ggplot(data=b1) + geom_histogram(aes(x=b1)) + theme_bw()
というように、いつもどおりの方法でできます。本エントリーの内容がわかっていれば，アヒル本のコードも読み解けるはず。コツコツ取り組んでいくと、Stan への理解が深まるだけでなく、R の操作や作図能力もぐっと上がると思います。
最後に
周りに統計に詳しくてすぐに質問できる人がいない社会科学系の社会人学生の方はいっしょにがんばっていきましょう！ そして、積極的に読書会や勉強会に参加しましょう！ 意外と年齢層高いから大丈夫。
Osaka.Stan は12月23日のようです（正式な募集はまだ）。
Hommachi.Stats は12月10日です。
