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章から先は少しずつやります。
統計的因果推論勉強会 第5回を開催しました！
経営学系統計分析エンドユーザーのため統計的因果推論勉強会第5回を行いました。最近週休3日制を打ち出した某有名大企業の新しいオフィス内での開催で，とても楽しかったです。参加者は7名（クローズド）。
今回は宮川本の第5章（の特にバックドア基準）と星野本第4章「共変量選択と無視できない欠測」が勉強範囲でした。目的は2つ。1つは強く無視できる割り当て条件とバックドア基準が似たようなものであることを理解すること。もう1つはRでの実習。Rnotebookを使うと便利でした。配布する資料もカンタンに作れるので。
試用したものはこちらにあります。
Rでの実習は傾向スコアマッチングを星野本の巻末付録に載っているもので写経（使用したデータはいろいろねつ造したフィクションのものです），IPW法を岩波データサイエンス第3巻のサポートページに載っているもので写経しました。
c統計量（このcは大文字のほうがいいのか，小文字なのかわかりません･･･）については，私たちはお手軽に rms パッケージの lrm() でやったほうがラクです。「C」と記載されているところがc統計量です。
星野本第4章までで，論文を読む・書くツールはだいたい出そろっているように思います。しかし，いろいろ謎は残ります。最大の謎は推定した因果効果の大きさ・小ささをどう判断したらよいか，です。効果量（effect size）のようにある程度の基準があれば便利なのですが･･･。
さて，最新の『行動計量学』には一般化傾向スコアを用いた論文が出ています。今度はそれもみんなで読んでみたいですね。
勉強会ツールとしての R Notebooks
RStudio をプロジェクタで映すと，ソースエディタが小さく見にくいです。必要なところだけを集中して，見やすくすることってできないだろうか，という動機から，もう1つの選択肢である R Notebooks（あーるのーとぶっく*1）を使えるようにしてみました。
RStudio Preview edition をインストールしてからの動作。
1. 左上の＋ボタンから R Notebook を選ぶ。
2. 「田」のアイコンをクリックして Zoom Source をクリック。そうすると，R Notebook だけが表示されるようになって，プロジェクタに映すとき便利になる。
3. メニューの View から Zoon In をクリックして，望みの大きさにする。元の大きさに戻すときは Acutual Size をクリックすればOK。
4. また元の4ペインに戻す場合は，「田」から Show All panes をクリック。
勉強会のときはこれで乗り切れそう。
Jupyter Notebook と R Notebooks の使い分け
Jupyter Notebook はブラウザ内で動作するので，勉強会中にこのウェブページ見せたいとかよくネット検索するなら Jupyter が便利だと思います。R Notebooks などは次の補足を参考にしてください。
補足
R Notebooks がどんなものかを使う前に知りたいならこちら。
わかりやすい kazutan 先生のスライド
R Notebooks を使う準備
1．RStudio Preview edition をインストール。 ダウンロード先
Desktop Version で OK。
2．基本 R Notebooks が使える状態になっている（2016年10月6日現在）。
3．でも，knit するときに「rmarkdown パッケージがないよ」と叱られるので，普通にインストールする。ただ特に分析結果を出力するつもりがないなら，インストールしなくてもいい思う。
*1:本当は複数形だけど，日本語に複数形はないので。
R を使うための Jupyter notebook インストールメモ
【2016年9月27日追記あり（いちばん下）】
--
土曜と日曜（2016年9月24日・25日），石川・品川の人たちと合宿に参加しました。2日目に3時間ほどのRの実習をする機会に恵まれ，つたないながらファシリテーター的な役割を務めました。ゼミは多国籍チーム（4カ国）なので，スライドは英語で，しゃべりは日本語中心*1。
RStudio を使ってやってみたところ，起動画面がぱっと見で複雑そうな印象を与えるかもしれない，とふと思いました（実際，起動した瞬間に参加者から「あ，なんか難しそう」という空気が広がったような。）。
もう1つの可能性として，jupyter notebook（じゅぴたーのーとぶっく）があるかなと思いインストール。注意点と手順を2ステップで半年後・1年後の自分のためにメモしておきます。mac 向け。
【ステップ1】
↓の1〜6．を実行。5．のとき anaconda のバージョンは最新のものにした。
【ステップ2】
↓に従って，jupyter notebook をインストール。
無事，jupyter notebook で R が使えるようになりました。でも，初心者向け R 実習に使うには厳しいです。理由は2つ。インストールが難しい（っていうか，面倒。時間もかかる）のと，コードの補完機能がないから（あれは初心者には便利のようです）。
とは言え，画面はスッキリして威圧感ない。なんか出力結果に罫線があるし，見やすい。個人的には好き。
そしてここまで書いて気付いた。 R notebookっていう選択肢があるじゃないか。また今度の機会にしよう･･･。
●参考
【追記1】
ターミナルで毎回 jupyter notebook と打つのは面倒なので，note というエイリアスを設定した。エイリアスの設定方法は次のページがわかりやすい。
【追記2】
コードの補完機能がない，と上で書きましたがウソでした。あります。コードを書いているときに，tab で出ます。
【追記3】
等幅フォントにするにはちょっと設定が必要。以下、mac の場合をメモ。
１．ホームディレクトリのライブラリフォルダに移動。
２．隠しフォルダ .jupyter を表示（コマンド＋シフト＋．）
３．その中に custom フォルダを作る。
４．さらにその中に custome.css を保存する。
CSSファイルの中身は↓。
.CodeMirror pre, .output pre { font-family: Monaco, monospace; }
デキる人はコマンドライン一発なのだろうが、ムリなのでアナログ（？）でやった。
参考ページは↓
http://rakuishi.com/archives/jupyter-font-family/
研究法の向き不向き
昨日から「次世代知識科学特論」というすごい名前のクラスが始まった。JAIST は人間力なんとかとか，大きなタイトルが好きみたい（私は嫌いじゃない）。どういう内容なのかは興味ある人はシラバスで検索してください。
この授業は，集中して質的な研究法も量的な研究法にも慣れましょう，という意図があるクラスのように思う。4日目と5日目は機械学習とかデータマイニングの講座になる。本来はドクターコース用だけど，今年からマスターの人たちにも開放されたらしい。そのため夜18時半〜22時までなのに，受講者がけっこう多い。
今日から2日間は文化人類学の先生によるエスノグラフィの演習。観察しながらメモを取るというのが難しい。途中であれこれと妄想してしまう。基本，量的なアプローチが好きだから，バランスを取るためにもエスノグラフィみたいなものもできるようになりたいと思う。でも，妄想が止まらない。きっとこの人の人生は〜とか，考えてしまう。
自分にそういう傾向があるということを知っただけでもよかったかも。ひょっとして，こういうアプローチには向いてないのかな･･･。まあ，まだ始まったばかりだし，明日もエスノグラフィの演習が続くので，少しずつ慣れていこう。在学中に1本，エスノグラフィで論文を書いてみたい･･･。
★ ★ ★
本日の進捗＝1643文字
