broom 0.3.7を使ってみる

R ggplot2 broom

R界隈の最新動向をストーキングすることに定評がある@u_ribo氏が、こんな情報をさらっと書いていました。情報網すごすぎる。

broomってそんなことできたっけ？？、と気になりすぎたのでbroomをちょっと触ってみます。

broomとは

たぶんこのVignetteを読むのが分かりやすいと思います。ここに書くコードもこのVigentteをパクり参考にしました。

broom: let's tidy up a bit

要は、クラスによってデータの構造が違うのを吸収して、統一的なフォーマットのきれい（tidy）なdata.frameをつくる、というものです。

tidy()

たとえば、以下のようなlmのオブジェクトがあるとします。

lmfit <- lm(mpg ~ wt, mtcars)
lmfit
#> 
#> Call:
#> lm(formula = mpg ~ wt, data = mtcars)
#> 
#> Coefficients:
#> (Intercept)           wt  
#>      37.285       -5.344

これを、色々モデルを変えてCoefficientsを比較したい、みたいな時はどうすればいいのでしょう。strで中身を見てみましょう。

str(lmfit)
#> List of 12
#>  $ coefficients : Named num [1:2] 37.29 -5.34
#>   ..- attr(*, "names")= chr [1:2] "(Intercept)" "wt"
#>  $ residuals    : Named num [1:32] -2.28 -0.92 -2.09 1.3 -0.2 ...
#>   ..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
#>  $ effects      : Named num [1:32] -113.65 -29.116 -1.661 1.631 0.111 ...
#>   ..- attr(*, "names")= chr [1:32] "(Intercept)" "wt" "" "" ...
...

この中から、目的の値を探し出すのはちょっと骨が折れます。しかもこれが、glmとか他のオブジェクトだとまた構造が違うわけなので探しなおしです。うむむ...

こういうめんどくささを避けるために、いちおうcoef(summary(lmfit))というイディオムは用意されています。

coef(summary(lmfit))
#>              Estimate Std. Error   t value     Pr(>|t|)
#> (Intercept) 37.285126   1.877627 19.857575 8.241799e-19
#> wt          -5.344472   0.559101 -9.559044 1.293959e-10

でも、これも完全ではありません。たとえばPr(>|t|)とかいうカラム名はクラスやモデルによって違うこともあるので、結果を組み合わせるときに苦戦します。

broomのtidy()は、このへんの違いを吸収してtidyなdata.frameをつくってくれる関数です。

tidy(lmfit) %>%
  knitr::kable(.)
term estimate std.error statistic p.value
(Intercept) 37.285126 1.877627 19.857575 0
wt -5.344472 0.559101 -9.559044 0

わあいきれいなテーブル！（小並感

もちろん、あらゆるクラスに対応しているわけではありません。対応しているクラスは、GithubレポジトリのREADMEのAvailable Tidiersという項目に表があります。

ということで、私はこのtidy()しか知らなかったんですが、augment()とglance()という関数がいつの間にか増えていました。

augment()

このaugment()が、ggplot2のfortify()に代わるもののようです。

ヘルプを見ると、

This generic originated in the ggplot2 package, where it was called "fortify."

ってばっちり書いてました。代わるもの、というか、ほんとにfortify()の延長線上にあるもののようです。

使い方はfortify()と同じです。（私あまり詳しくないので、間違ってたらツッコミをください…）

head(augment(lmfit))
#>           .rownames  mpg    wt  .fitted   .se.fit     .resid       .hat   .sigma      .cooksd  .std.resid
#> 1         Mazda RX4 21.0 2.620 23.28261 0.6335798 -2.2826106 0.04326896 3.067494 1.327407e-02 -0.76616765
#> 2     Mazda RX4 Wag 21.0 2.875 21.91977 0.5714319 -0.9197704 0.03519677 3.093068 1.723963e-03 -0.30743051
#> 3        Datsun 710 22.8 2.320 24.88595 0.7359177 -2.0859521 0.05837573 3.072127 1.543937e-02 -0.70575249
#> 4    Hornet 4 Drive 21.4 3.215 20.10265 0.5384424  1.2973499 0.03125017 3.088268 3.020558e-03  0.43275114
#> 5 Hornet Sportabout 18.7 3.440 18.90014 0.5526562 -0.2001440 0.03292182 3.097722 7.599578e-05 -0.06681879
#> 6           Valiant 18.1 3.460 18.79325 0.5552829 -0.6932545 0.03323551 3.095184 9.210650e-04 -0.23148309

head(fortify(lmfit))
#>                    mpg    wt       .hat   .sigma      .cooksd  .fitted     .resid   .stdresid
#> Mazda RX4         21.0 2.620 0.04326896 3.067494 1.327407e-02 23.28261 -2.2826106 -0.76616765
#> Mazda RX4 Wag     21.0 2.875 0.03519677 3.093068 1.723963e-03 21.91977 -0.9197704 -0.30743051
#> Datsun 710        22.8 2.320 0.05837573 3.072127 1.543937e-02 24.88595 -2.0859521 -0.70575249
#> Hornet 4 Drive    21.4 3.215 0.03125017 3.088268 3.020558e-03 20.10265  1.2973499  0.43275114
#> Hornet Sportabout 18.7 3.440 0.03292182 3.097722 7.599578e-05 18.90014 -0.2001440 -0.06681879
#> Valiant           18.1 3.460 0.03323551 3.095184 9.210650e-04 18.79325 -0.6932545 -0.23148309

fortify()の結果と比べてみると、違いは、rownameを使うのをやめて.rownamesというカラムになっていることと、.se.fitが増えていることくらいでしょうか。あまり大差はなさそうです。

対応しているクラスにはけっこう違いがあります。置き換わるのがいつごろになるかはなんとも言えなそうです。

methods("augment") %>%
  as.character %>%
  gsub("^augment\\.", "", .)
#>  [1] "coxph"         "data.frame"    "default"       "felm"          "kmeans"        "lm"            "lme"           "loess"         "merMod"       
#> [10] "mlm"           "nls"           "NULL"          "plm"           "rowwise_df"    "smooth.spline" "survreg"

methods("fortify") %>%
  as.character %>%
  gsub("^fortify\\.", "", .)
#>  [1] "cld"                      "confint.glht"             "data.frame"               "default"                  "glht"                    
#>  [6] "Line"                     "Lines"                    "lm"                       "map"                      "NULL"                    
#> [11] "Polygon"                  "Polygons"                 "SpatialLinesDataFrame"    "SpatialPolygons"          "SpatialPolygonsDataFrame"
#> [16] "summary.glht"

glance()

glance()は、R2 とか尤度とか、モデル全体に関わる統計量を1行のdata.frameにまとめる関数です。

glance(lmfit)
#>   r.squared adj.r.squared    sigma statistic      p.value df    logLik      AIC      BIC deviance df.residual
#> 1 0.7528328     0.7445939 3.045882  91.37533 1.293959e-10  2 -80.01471 166.0294 170.4266 278.3219          30

実例

kmeans法で、クラスタ数を変えて結果を見る、というのをbroomでやるとこうなる。というのが下のvignetteです。なんかすごく便利そうです。

Tidying k-means clustering

（この例ではいったんひとつのdata.frameにまとめてからggplot2でfacetしてますが、別々のdata.frameのまま別々にグラフをつくってgridExtra::grid.table()でまとめるほうがいいのでは？という気もします）

感想

いつのまにこんな進化してたんだbroom。すごい。