niszetの日記

10年目エンジニアが最近勉強したことを忘れないようにメモをする

(R調べもの) broom::tidyとかtidyMCMCとか (7/26 追記)

僕のTokyoRは終わっていない…

なぜなら復習をしていないからです。前回の日記の末尾に書いたものは何だったのか。

気を取り直して。

R for Data Scienceを読みましょう

TokyoRでも話がちょっと出ていましたが、R for Data Scienceの和訳が出るのはいつになることか。
いっそ、原著を買って読みましょう。私は3か月前に読み終えました^^
久々にページを開き、ほとんど忘れていることを確認済みです。

オライリーはもう直販しないんですっけ?
まぁ、amazonで購入できますので(私はamazonで購入しましたが、届くまでにやたら待ちました…)

shop.oreilly.com

また、こちらでも読むことが出来ますよ。

r4ds.had.co.nz

broomパッケージ

tidyverseには色々なパッケージが含まれています。tidyrとかreadrとかpurrrとか、ggplot2とかいろいろ。
でも、modelrやbroomって全然話題に上がりませんよね。

と思って検索したら…

notchained.hatenablog.com

そして

uribo.hatenablog.com

やはり情報が早い…。
しかしめげずに書かれていないことを探して書くのです。

broom::tidy

uryuさんの記事にもありますが、broom::tidyという関数があります。これは、lmやglmなどのclassをもったデータ(結果)をtidyにしてくれるというもの。戻り値は必ずdata.frameになるとのこと(ヘルプより)
いちいち結果をパースしなくてよいので便利ですね。t.testとか、wilcox.testの結果にも使えるようです。glanceとかaugmentなどもよさげです。
この辺りはリンク先やヘルプを見ると良いのです。追加で書くことは特にないのです。

で、色々見てる中でちょっと気づいたのですが、broom::tidyMCMCという関数があるのです。
これね。

www.rdocumentation.org

犬四匹本もめでたく出版されましたし、MCMC系の結果のtidy化はこれから使いどころがあるかも…?
ということで、ちょっとこの関数を触ってみようと思います。

ちなみに、犬4匹本とはこちらの、

www.kyoritsu-pub.co.jp

のことです。しかしその大きさ、厚さから、犬4匹本はポストに入らなかったために手元にありません…。

tidy化するには結果が必要なため、アヒル本のプログラムを使って作ってstanの結果を作ってみます。

ちなみに、アヒル本とはこちら

www.kyoritsu-pub.co.jp

のことです。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!