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日です。
5章まで読了『StanとRでベイズ統計モデリング』
昨日はJapan.Rに初参加。4〜5時間ずっとスクリーンを見ていたので眼精疲労がスゴい。8割方よくわからなかった（領域が違うので）けれど，おもしろい収穫もありました*1。
さて，『StanとRでベイズ統計モデリング』を5章まで読み進めました。私は出張が多いので，移動時間にちょくちょく読むスタイル。実際のコードは新幹線の中や在来線特急の中で動かしています。1つStanコードを動かすのに4〜5回は修正します。大文字にすべきところを小文字にしていたり，スペルミスだったり，セミコロンが抜けていたり･･･。根気よく慣れるしかないのでしょう。
とはいえ，4章・5章でキツいのは Stan そのものの理解というより，R・ggplot2 の理解のほうでした。だから1つ章を終わらせるのに，4週間ほどかかった･･･。特に ggplot2 はそんなに慣れていないから大変。geom_ribbon, geom_pointrange とかなんとか使ったことがないので，その理解に時間を費やしました。もう勢いで “ggplot2: Elegant Graphics for Data Analysis” も注文したくらい（まだ届いてないけど）*2。
ggplot2: Elegant Graphics for Data Analysis (Use R!)
- 作者: Hadley Wickham
- 出版社/メーカー: Springer
- 発売日: 2016/06/16
- メディア: ペーパーバック
- この商品を含むブログを見る
R・ggplot2 にそんなに慣れていない人が詰まりそうな部分がわかってきたので（自分のことなので），12月4日分の Stan Advent Calender で書く予定です。
4章まで読了。『StanとRでベイズ統計モデリング』
当初1〜8章まででいいかなと目次を見てつぶやいたら，著者の方から「10章も読んだほういい」というアドバイスをいただきました(感謝です。読むモチベーションもあがります)。ということで，1〜8章＋10章を読むことを目標にします*1。自分の復習として感想を記録しておきます。
環境構築
今年からMacを使うようになり，Stanを使える環境はまだ作っていませんでした。次のページのとおりにやれば，問題なく完了しました。
https://github.com/stan-dev/rstan/wiki/Installing-RStan-on-Mac-or-Linux
※「If using g++ version 4.9 or higher〜〜」のところは不要だった。
読む前に
これまで1人でコソコソ学んできて（Kobe.Rを中心に他のRの勉強会には参加したことがあります），久保先生の『データ解析のための統計モデリング入門（みどり本）』と豊田先生の『基礎からのベイズ統計学』をなんとか読み通しています。なので，「はじめに」に書いてある「前提とする知識」はなんとかクリアしていると思います。
1〜3章は在来線特急の中で読める。
今週の出張のお供として読んでいます。1〜3章までは実際にコードを書いてみるところはほぼないので，電車の中で揺られても読めます。読んでみた第1印象は「とってもていねい！」ということ。用語が1つ1つ説明されています。「thinning って何？」「lp__って何？」とこれまで疑問に思っていたのですが，解決することができました*2。
3章2節には統計モデリングの手順があります。こういうことって，たぶん身近に詳しい人がいれば体得できるんでしょうが，私のような独学のものには助かります。
辛うじて突破の4章
Stanコードの書き方がていねいに解説されています。「そういう順番（やルール）で書けばいいのか」と目からウロコでした。私はプログラマでもデータサイエンティストでもなく，コードを書くときのルールをよく知らないので，こういう内容は助かります。
この章から実習があります。ついつい先を読んでしまいますが，ただ通読しているだけでは，わかった気になるので，実際に打ち込んで行きました。N とすべきところを n としていたり，ささいなことでエラーになるのを試行錯誤しながら，取り組んでいくことになるので，経験値が増えそうです。あと，model4-4.stan を読んでいて「このN_newって何？ 突然出てきた！」と戸惑いましたが，サンプルコードを読んで解決できました。
52ページの解説は自分のR力のなさで，ちょっとぼんやりとした理解になりました（apply が苦手）。練習問題も（3）まではできたけど，R力のなさで（4）ができなかった…。残念。答えを見て，こんな簡単なものを解けない自分に自己嫌悪。
5章から先は少しずつやります。
