読者です 読者をやめる 読者になる 読者になる

Knowledge As Practice

サービス経営の研究をしている社会人大学院生の統計分析勉強メモです。それ以外のネタは Medium に書いています。

rstanで分析をするのに必要なR力(りょく)

Stan と R でベイズ統計モデリング ベイズデータ分析 統計分析 ノンプログラマーのためのR入門

Stan Advent Calendar 2016 の 4日目です。非マニア枠です。

 
本日のエントリーは Stan を普通に使える人、情報科学・データ科学な人を対象にしていません。少し R に慣れてきて、ベイズ統計モデリングにも興味を持ち始めた人向け。ずばりターゲットは次のようなレベルの人です。

 

  • いざ R を使って分析するときには使い方を忘れていて検索しながら使う。

  • dplyr の使い方に慣れておらず、めんどくさくなって結果エクセルを使ってしまう。

 
つまり私。そんな中途半端なレベルな人(私)向けに『Stan と R でベイズ統計モデリング』(以下,アヒル本とします*1)を読んでいるときにつまずくであろうところを書いておきます。アヒル本を読む・理解するには、適度に R と ggplot2 を使えるスキルが必要だからです。

 
さて,R を使った統計分析の本をいくつか読んで、簡単な分析をできるようになったという人が Stan を使ったベイズ統計モデリングをするときにつまずくのは次のようなところだと思います。

  1. 結果が出るまで待ち時間が長い。
  2. リスト形式でデータを渡すのを戸惑う。
  3. 結果が入ったオブジェクトを表示させたらズラーッと出てきてしまい、見たい結果を見るために必死に上へスクロールする。
  4. タイプミスする。
  5. 結果からほしいデータを取り出せない。

 

1. 待ち時間が長い

私の場合,簡単なモデルでも30秒くらいかかります。 f:id:hikaru1122:20161203224300p:plain

 
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

 
f:id:hikaru1122:20161203174841g:plain
どばっと結果が出ます。そしてパラメータを確認したいので上へ必死にスクロール。

 
こういうのを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日です。

*1:ほんとは Rstan 本なんでしょうけども,表紙がアヒルというのが衝撃だったので。

*2:これを解消するのに、glmmStan や rstanarm や brms を使えば大丈夫ですけど・・・。

*3:たしか、DataCamp でこんな説明が出てきた気がします。

*4:ひょっとしたら本の中にそういう記述があったかもしれない。

*5:RStudioなどは補完してくれますが、それでもつい消しちゃって足りないときもあるんです。

5章まで読了『StanとRでベイズ統計モデリング』

ベイズデータ分析 ノンプログラマーのためのR入門 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!)

ggplot2: Elegant Graphics for Data Analysis (Use R!)

 
R・ggplot2 にそんなに慣れていない人が詰まりそうな部分がわかってきたので(自分のことなので),12月4日分の Stan Advent Calender で書く予定です。

Stan Advent Calendar 2016 - Qiita

*1:R Analytic Flowとかベイズの話とか発表者が不在でも発表ができるとか。

*2:修行が足りず,こういう本は電子書籍だと読めない派なのです・・・。

4章まで読了。『StanとRでベイズ統計モデリング』

統計分析 ベイズデータ分析 ノンプログラマーのための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章から先は少しずつやります。

*1:来年1月か2月までに読めたらいいなぁ,というゆっくりペースで行きます。JAGSを使った分析結果の論文を年内に書くのが優先なので。

*2:でも,正しいlp__の読み方は何だろう? 勝手に「えるぴーちょめちょめ」と読んでます。

統計的因果推論勉強会 第5回を開催しました!

統計的因果推論 統計分析 ノンプログラマーのためのR入門

経営学系統計分析エンドユーザーのため統計的因果推論勉強会第5回を行いました。最近週休3日制を打ち出した某有名大企業の新しいオフィス内での開催で,とても楽しかったです。参加者は7名(クローズド)。

speakerdeck.com

 
今回は宮川本の第5章(の特にバックドア基準)と星野本第4章「共変量選択と無視できない欠測」が勉強範囲でした。目的は2つ。1つは強く無視できる割り当て条件とバックドア基準が似たようなものであることを理解すること。もう1つはRでの実習。Rnotebookを使うと便利でした。配布する資料もカンタンに作れるので。

試用したものはこちらにあります。

github.com

 
Rでの実習は傾向スコアマッチングを星野本の巻末付録に載っているもので写経(使用したデータはいろいろねつ造したフィクションのものです),IPW法を岩波データサイエンス第3巻のサポートページに載っているもので写経しました。

 
c統計量(このcは大文字のほうがいいのか,小文字なのかわかりません・・・)については,私たちはお手軽に rms パッケージの lrm() でやったほうがラクです。「C」と記載されているところがc統計量です。

 
f:id:hikaru1122:20161029205730p:plain

 
星野本第4章までで,論文を読む・書くツールはだいたい出そろっているように思います。しかし,いろいろ謎は残ります。最大の謎は推定した因果効果の大きさ・小ささをどう判断したらよいか,です。効果量(effect size)のようにある程度の基準があれば便利なのですが・・・。

 
さて,最新の『行動計量学』には一般化傾向スコアを用いた論文が出ています。今度はそれもみんなで読んでみたいですね。

勉強会ツールとしての R Notebooks

統計分析 研究に役立つツール ノンプログラマーのためのR入門

RStudio をプロジェクタで映すと,ソースエディタが小さく見にくいです。必要なところだけを集中して,見やすくすることってできないだろうか,という動機から,もう1つの選択肢である R Notebooks(あーるのーとぶっく*1)を使えるようにしてみました。

RStudio Preview edition をインストールしてからの動作。

1. 左上の+ボタンから R Notebook を選ぶ。 f:id:hikaru1122:20161007013327p:plain  

2. 「田」のアイコンをクリックして Zoom Source をクリック。そうすると,R Notebook だけが表示されるようになって,プロジェクタに映すとき便利になる。 f:id:hikaru1122:20161007013345p:plain  

3. メニューの View から Zoon In をクリックして,望みの大きさにする。元の大きさに戻すときは Acutual Size をクリックすればOK。 f:id:hikaru1122:20161007013412p:plain

 
4. また元の4ペインに戻す場合は,「田」から Show All panes をクリック。 f:id:hikaru1122:20161007013424p:plain

 
勉強会のときはこれで乗り切れそう。

 

Jupyter Notebook と R Notebooks の使い分け

 
Jupyter Notebook はブラウザ内で動作するので,勉強会中にこのウェブページ見せたいとかよくネット検索するなら Jupyter が便利だと思います。R Notebooks などは次の補足を参考にしてください。

 

補足

 

R Notebooks がどんなものかを使う前に知りたいならこちら。

www.youtube.com

 

わかりやすい kazutan 先生のスライド

RStudioで R Notebook機能を試してみる

 

R Notebooks を使う準備

1.RStudio Preview edition をインストール。 ダウンロード先

www.rstudio.com

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 のバージョンは最新のものにした。

qiita.com

 

【ステップ2】

↓に従って,jupyter notebook をインストール。

Installation · IRkernel

 
無事,jupyter notebook で R が使えるようになりました。でも,初心者向け R 実習に使うには厳しいです。理由は2つ。インストールが難しい(っていうか,面倒。時間もかかる)のと,コードの補完機能がないから(あれは初心者には便利のようです)。

 
とは言え,画面はスッキリして威圧感ない。なんか出力結果に罫線があるし,見やすい。個人的には好き。

f:id:hikaru1122:20160926022201p:plain

 
そしてここまで書いて気付いた。 R notebookっていう選択肢があるじゃないか。また今度の機会にしよう・・・。

 
●参考

speakerdeck.com

 

【追記1】

ターミナルで毎回 jupyter notebook と打つのは面倒なので,note というエイリアスを設定した。エイリアスの設定方法は次のページがわかりやすい。

qiita.com

【追記2】

コードの補完機能がない,と上で書きましたがウソでした。あります。コードを書いているときに,tab で出ます。

*1:第一声は英語にしたけど,英語と日本語を交えながらやっていたら,だんだんシン・ゴジラ石原さとみっぽくなってきた感じがしてやめた。

研究法の向き不向き

その他 社会人博士

昨日から「次世代知識科学特論」というすごい名前のクラスが始まった。JAIST人間力なんとかとか,大きなタイトルが好きみたい(私は嫌いじゃない)。どういう内容なのかは興味ある人はシラバスで検索してください。

www.jaist.ac.jp

 
この授業は,集中して質的な研究法も量的な研究法にも慣れましょう,という意図があるクラスのように思う。4日目と5日目は機械学習とかデータマイニングの講座になる。本来はドクターコース用だけど,今年からマスターの人たちにも開放されたらしい。そのため夜18時半〜22時までなのに,受講者がけっこう多い。

 
今日から2日間は文化人類学の先生によるエスノグラフィの演習。観察しながらメモを取るというのが難しい。途中であれこれと妄想してしまう。基本,量的なアプローチが好きだから,バランスを取るためにもエスノグラフィみたいなものもできるようになりたいと思う。でも,妄想が止まらない。きっとこの人の人生は〜とか,考えてしまう。

 
自分にそういう傾向があるということを知っただけでもよかったかも。ひょっとして,こういうアプローチには向いてないのかな・・・。まあ,まだ始まったばかりだし,明日もエスノグラフィの演習が続くので,少しずつ慣れていこう。在学中に1本,エスノグラフィで論文を書いてみたい・・・。

 
★ ★ ★
本日の進捗=1643文字

クリエイティブ・コモンズ・ライセンス
この 作品 は クリエイティブ・コモンズ 表示 - 継承 4.0 国際 ライセンスの下に提供されています。