僕のTokyoRは終わっていない…
なぜなら復習をしていないからです。前回の日記の末尾に書いたものは何だったのか。
気を取り直して。
R for Data Scienceを読みましょう
TokyoRでも話がちょっと出ていましたが、R for Data Scienceの和訳が出るのはいつになることか。
いっそ、原著を買って読みましょう。私は3か月前に読み終えました^^
久々にページを開き、ほとんど忘れていることを確認済みです。
オライリーはもう直販しないんですっけ？
まぁ、amazonで購入できますので（私はamazonで購入しましたが、届くまでにやたら待ちました…）
また、こちらでも読むことが出来ますよ。
broomパッケージ
tidyverseには色々なパッケージが含まれています。tidyrとかreadrとかpurrrとか、ggplot2とかいろいろ。
でも、modelrやbroomって全然話題に上がりませんよね。
と思って検索したら…
そして
やはり情報が早い…。
しかしめげずに書かれていないことを探して書くのです。
broom::tidy
uryuさんの記事にもありますが、broom::tidyという関数があります。これは、lmやglmなどのclassをもったデータ（結果）をtidyにしてくれるというもの。戻り値は必ずdata.frameになるとのこと（ヘルプより）
いちいち結果をパースしなくてよいので便利ですね。t.testとか、wilcox.testの結果にも使えるようです。glanceとかaugmentなどもよさげです。
この辺りはリンク先やヘルプを見ると良いのです。追加で書くことは特にないのです。
で、色々見てる中でちょっと気づいたのですが、broom::tidyMCMCという関数があるのです。
これね。
犬四匹本もめでたく出版されましたし、MCMC系の結果のtidy化はこれから使いどころがあるかも…？
ということで、ちょっとこの関数を触ってみようと思います。
ちなみに、犬4匹本とはこちらの、
のことです。しかしその大きさ、厚さから、犬4匹本はポストに入らなかったために手元にありません…。
tidy化するには結果が必要なため、アヒル本のプログラムを使って作ってstanの結果を作ってみます。
ちなみに、アヒル本とはこちら
のことです。Rでstan/mcmc使うならコレ読んどけば間違いないといわれていますね。是非犬4匹本とともにお手元にどうぞ。
階層モデル以降は当時の自分のレベル的にちょっときつかったのですが、今なら行けそうかな…
アヒル本 model4-5のコードを使用します。
載せて良いのか確認とっていないので…プログラムの中身はサポートサイト（とアヒル本）をご参照ください…
今回はfitという変数に結果が入っているので、表示させてみます。
fit # Inference for Stan model: model4-5. # 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 # a -117.83 2.12 71.90 -262.76 -165.33 -116.11 -69.52 17.44 1148 1 # b 21.86 0.05 1.61 18.91 20.77 21.81 22.89 25.01 1123 1 # sigma 84.33 0.37 14.71 61.43 73.93 82.12 92.76 118.35 1609 1 # lp__ -93.57 0.04 1.23 -96.74 -94.17 -93.25 -92.66 -92.13 899 1 # Samples were drawn using NUTS(diag_e) at Tue Jul 25 22:54:39 2017. # 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).
tidyあるいはtidyMCMCの引数に与えると…
tidy(fit) # term estimate std.error # 1 a -117.83500 71.896853 # 2 b 21.86044 1.606462 # 3 sigma 84.32789 14.711962
一応、クラスを確認
class(tidy(fit)) # [1] "data.frame"
data.frameですね。
しかし、普通に表示させた結果と比べると物足りないですね。オプションをいじって色々表示させてみます
tidy(fit,estimate.method = "mean", conf.int = TRUE, droppars = NULL, rhat=TRUE, ess=TRUE, conf.level = 0.95, conf.method = "quantile") # term estimate std.error conf.low conf.high rhat ess # 1 a -117.83500 71.896853 -262.76378 17.44491 1.001971 1148 # 2 b 21.86044 1.606462 18.90667 25.01078 1.001802 1123 # 3 sigma 84.32789 14.711962 61.42625 118.34501 1.000426 1609 # 4 lp__ -93.56818 1.233635 -96.74103 -92.12595 1.004600 899
estimate.methodは"mean"か"median"が指定でき、conf.methodには"quantile"か"HPDinterval"が指定できます。その他、詳細はヘルプを参照下さい。
ただし、上記の結果に対してはHPDintervalは指定するとエラーになります。まだちゃんと理解できていませんが、実行したstanのプログラムがHPDintervalで出力していないからだと思いますが、まだちょっとわかっておりません。。。。
アヒル本著者のberoberoさんによれば、犬4匹本はHDIを使用している、ということなので、そちらで試せばよいのかなと思っています。まだ届いていないけど。 statmodeling.hatenablog.com
既知の情報かもしれませんが、気になったので一応まとめてみました。マサカリ下さいな。
おわりに
本当は、modelrとbroomにあるbootstrap()を調べて書こうと思ったのですが、もうちょっと調べたいのです…
Rは趣味でやっている関係で本業が忙しくなると色々止まります。今はちょっと余裕ある（本当か…？）ので、気になることを調べておきたいですね
追記 7/26
念のための追記。
broomはinstall.packages(tidyverse)で一緒にインストールされますが、library(tidyverse)ではロードされません。明示的にlibrary(broom)としてあげてください。
また、tidyverseに含まれているパッケージの確認は下記のようにすればできます(library(tidyverse)でロードしたあとはtidyverse::は不要です）
tidyverse::tidyverse_packages() # [1] "broom" "dplyr" "forcats" "ggplot2" "haven" "httr" "hms" "jsonlite" # [9] "lubridate" "magrittr" "modelr" "purrr" "readr" "readxl" "stringr" "tibble" #[17] "rvest" "tidyr" "xml2" "tidyverse"
さらに、下記で含まれているパッケージのバージョン、localとcranのものが比較できるようです。
tidyverse::tidyverse_deps() # # A tibble: 19 x 4 # package cran local behind # <chr> <chr> <chr> <lgl> # 1 broom 0.4.2 0.4.2 FALSE # 2 dplyr 0.7.2 0.7.2 FALSE # 3 forcats 0.2.0 0.2.0 FALSE # 4 ggplot2 2.2.1 2.2.1 FALSE # 5 haven 1.1.0 1.1.0 FALSE # 6 hms 0.3 0.3 FALSE # 7 httr 1.2.1 1.2.1 FALSE # 8 jsonlite 1.5 1.5 FALSE # 9 lubridate 1.6.0 1.6.0 FALSE # 10 magrittr 1.5 1.5 FALSE # 11 modelr 0.1.1 0.1.1 FALSE # 12 purrr 0.2.2.2 0.2.2.2 FALSE # 13 readr 1.1.1 1.1.1 FALSE # 14 readxl 1.0.0 1.0.0 FALSE # 15 rvest 0.3.2 0.3.2 FALSE # 16 stringr 1.2.0 1.2.0 FALSE # 17 tibble 1.3.3 1.3.3 FALSE # 18 tidyr 0.6.3 0.6.3 FALSE # 19 xml2 1.1.1 1.1.1 FALSE
あと、本文中でtidyを使ってstanの結果を見てみましたが、tidyMCMCと比較すると
identical(broom:::tidy.stanfit, broom:::tidyMCMC) # [1] TRUE
でしたので問題ないかなと。tidyはS3 generic functionなので、クラスがstanfitの場合はtidy.stanfitが呼ばれるのです。
class(fit) # [1] "stanfit" #attr(,"package") # [1] "rstan" typeof(fit) # [1] "S4" is(fit) # [1] "stanfit" str(fit) # 省略
また、
broom::augment(lm(mpg ~ wt, mtcars)) # 途中省略 # Warning messages: # 1: Deprecated: please use `purrr::possibly()` instead # うしろも省略
で、augmentはDeprecated、purrr::possiblyの使用を使えとのことでした。 glanceの方はといえば
broom::glance(lm(mpg ~ wt, mtcars)) # r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance #1 0.7528328 0.7445939 3.045882 91.37533 1.293959e-10 2 -80.01471 166.0294 170.4266 278.3219 # df.residual #1 30
で問題なく動いているので、一体何が気に入らなかったのかよくわかりません…purrr::possiblyについてはいずれ見てみようかなーと思います。
Enjoy!
