【R】 broomパッケージ入門: Rのモデル出力を整然化(tidy)する
入門記事というより、俺が入門(おさらい)した内容をまとめています。
分析結果を自動的に出力するために modelsummary::msummary() をもう少し深く勉強しているうちに、broomパッケージの内容をもうちょっときちんと理解せんといかんな、と思っておさらいすることにしました。
broom パッケージとは?
Rにおける様々な統計分析モデルで分析した結果って、各々の関数のオブジェクトとして出力されるのですが、
それを一律にtidyなオブジェクトとして出力しよう、というアイデアのパッケージです。
Rで行った分析をそのまま見るだけなら、コンソールに出力して見るだけでいいのですが、
レポートに含めたり、また分析された結果を加工したり別のモデルに利用する時に、出力されたデータがtidyでないと煩雑なプロセスを踏んで加工する必要が出てきます。そういった手間を省き、整然な出力を行うことができるのがbroomパッケージです。
整然データ (tidy data)とは?
R界の神、Hadley Wickam氏によって提案された一律の構造に則ったデータ形式のことです。
基本的に各行(→)が一つ一つの観察値であり、各列(↓)が一つの変数である形を持っています。
Rのオブジェクトは主にリストの形で、いろいろな形のオブジェクトを一つにまとめたような形になっていますが、関数によって形式や含む要素がバラバラのため、扱いにくいため、これをtidy dataのルールに則って出力するパッケージがbroomということになります。
broomの3つの関数
broomには3つの異なった整然化を行うメソッドがあります。
1. tidy : 統計モデルの結果を整然データとして出力する。
2. augment : 統計モデルに入力する元データに、出力結果を付加する。(例:予測値など)
3. glance : 統計モデルの結果に関する一行のサマリーを出力する。
一つ一つ見ていきます。
1. tidy: 係数などの整然化
例としてmtcarsのデータを使って線形回帰を行います。
出力された結果の
library(broom) lmfit <- lm(mpg ~ wt, mtcars)
出力されたオブジェクトはlmというオブジェクトです。
> class(lmfit) [1] "lm"
一般的には summary()関数を使って結果をコンソールに出力します。
> summary(lmfit) Call: lm(formula = mpg ~ wt, data = mtcars) Residuals: Min 1Q Median 3Q Max -4.5432 -2.3647 -0.1252 1.4096 6.8727 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 37.2851 1.8776 19.858 < 2e-16 *** wt -5.3445 0.5591 -9.559 1.29e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.046 on 30 degrees of freedom Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446 F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10
もちろんコンソールで見る限りは見やすい結果なのですが、この結果をさらに加工したり、係数だけ取り出して図を作ったりするのに手間がかかります。
しかし、tidy()を使えば一発で整然化されたデータフレームで出力されます。
> tidy(lmfit) # A tibble: 2 x 5 term estimate std.error statistic p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) 37.3 1.88 19.9 8.24e-19 2 wt -5.34 0.559 -9.56 1.29e-10
上では、各行が各係数の結果、各列は係数についての推定値、標準誤差、t-統計量、p値がそれぞれ入っています。
これは各行が観察値、各列が変数、という整然データの概念に一致しています。
2. augment: 予測値などの整然化
lmによって推定した結果を用いて予測値を追加したいとします。元データの各行(すなわち各観察値)に対して予測値を追加したい場合は以下のコマンドで簡単に追加できます。
> augment(lmfit) # A tibble: 32 x 10 .rownames mpg wt .fitted .se.fit .resid .hat .sigma .cooksd .std.resid <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 Mazda RX4 21 2.62 23.3 0.634 -2.28 0.0433 3.07 0.0133 -0.766 2 Mazda RX4 Wag 21 2.88 21.9 0.571 -0.920 0.0352 3.09 0.00172 -0.307 3 Datsun 710 22.8 2.32 24.9 0.736 -2.09 0.0584 3.07 0.0154 -0.706 4 Hornet 4 Drive 21.4 3.22 20.1 0.538 1.30 0.0313 3.09 0.00302 0.433 5 Hornet Sportabout 18.7 3.44 18.9 0.553 -0.200 0.0329 3.10 0.0000760 -0.0668 6 Valiant 18.1 3.46 18.8 0.555 -0.693 0.0332 3.10 0.000921 -0.231 7 Duster 360 14.3 3.57 18.2 0.573 -3.91 0.0354 3.01 0.0313 -1.31 8 Merc 240D 24.4 3.19 20.2 0.539 4.16 0.0313 3.00 0.0311 1.39 9 Merc 230 22.8 3.15 20.5 0.540 2.35 0.0314 3.07 0.00996 0.784 10 Merc 280 19.2 3.44 18.9 0.553 0.300 0.0329 3.10 0.000171 0.100 # … with 22 more rows
推定に使用したmpgとwtという変数に加えて、以下の変数が追加されている。追加された変数は . (ドット)から始まる。
- 予測値(.fitted),
- 予測値の標準誤差(.se.fit),
- 残差 (.resid),
- 標準化残渣 (.std.resid),
- 射影行列の対角(.hat)
- 各行が除かれて推定された場合の残差標準誤差の推定値 (.sigma)
- 各行が除かれて推定された場合の予測値の変化(クック距離) (.cooksd)
3. glance: モデル適合度などの整然化
glance()関数を使うと、モデルの出力のうち、モデル適合値など一つの値で代表される統計量を一行のデータフレームで示す。
> glance(lmfit) # A tibble: 1 x 11 r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> 1 0.753 0.745 3.05 91.4 1.29e-10 2 -80.0 166. 170. 278. 30
broomの使用例
etimatrの結果
Rに最初から入っているパッケージだけではなく、様々なパッケージが利用できる。
estimatrの結果もtidyにできるが、lm()とは少し出力される変数が異なる。
> library(estimatr) > # mtcarsを使った推定 > lmrob <- lm_robust(mpg ~ wt, mtcars) > tidy(lmrob) term estimate std.error statistic p.value conf.low conf.high df outcome 1 (Intercept) 37.285126 2.2689318 16.43290 1.514499e-16 32.651349 41.918903 30 mpg 2 wt -5.344472 0.6832764 -7.82183 9.959194e-09 -6.739908 -3.949035 30 mpg
dplyrのパイプの中で使う。
これは効果検証入門の2.2.7からの抜粋ですが、
cor_visit_treatment <- lm( data = biased_data, formula = treatment ~ visit + channel + recency + history) %>% tidy()
のように、lmで推定したものをパイプ(%>%)でtidy関数に受け渡してしまうと、そのまま整然化されたデータフレームが取得できます。
まとめて回帰分析を実行する。
これも効果検証入門から抜粋です。(2.2.3)
以下のコードは事前にmodelsというデータフレームに複数の推定式を入れています。
二行目で、各モデル(formula)に対して、lmという関数をbiased_dataというデータを使って実行するという繰り返し(map)を行い、
3行目で、得られた各結果(model)に対して、tidyという関数をそれぞれ実行するという命令になっています。
df_models <- models %>% mutate(model = map(.x = formula, .f = lm, data = biased_data)) %>% mutate(lm_result = map(.x = model, .f = tidy))
エレガントですね。
詳しいコードについては効果検証入門を参照してください。
コードはこちらにあります。
使えるモデル・パッケージ
broomを使えるモデルやパッケージはどんどん増えていて、全てを挙げるのは難しいですが、計量経済学関連だと以下のようなパッケージが使えるのがありがたいです。(ご存じだと思いますがパッケージ::関数という表記です。)
- estimatr: ロバスト標準誤差、操作変数、固定効果など基本的な計量経済学の推定ができる。estimatr::lm_robust, iv_robustなど。
- lfe : 複数の固定効果モデルを含めて推定できるlfe::felmなど。
- fixest : 複数の固定効果モデルを速く推定するfixest::feolsなど
- mfx : ロジット・プロビットなどの限界効果を計算するmfx::mfx, logitmfx, negbinmfx, poissonmfx, probitmfxなど
- spatialreg : 空間回帰分析のできるspatialreg::sarlmなど
- gmm : GMM(一般化モーメント法)の推定ができるgmm::gmmなど。
参考文献
効果検証入門
tidyverse ブログ
https://www.tidyverse.org/blog/2020/07/broom-0-7-0/#new-tidier-methods
Introduction to broom
https://cran.r-project.org/web/packages/broom/vignettes/broom.html