データ分析メモと北欧生活

旧Untitled Note. データ分析、計量経済・統計とR、水産管理、英語勉強、海外生活などについて備忘録や自分の勉強のOutputの場所として

MENU

【R】 broomパッケージ入門: Rのモデル出力を整然化(tidy)する

入門記事というより、俺が入門(おさらい)した内容をまとめています。

分析結果を自動的に出力するために modelsummary::msummary() をもう少し深く勉強しているうちに、broomパッケージの内容をもうちょっときちんと理解せんといかんな、と思っておさらいすることにしました。

keita43a.hatenablog.com

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
  • r.squared : R二乗値
  • adj.r.squared : 調整済みR二乗値
  • sigma : 残差平方和
  • statistic : F統計量
  • p.value : F検定のp値
  • df : 係数の自由度
  • logLik : 対数尤度
  • AIC : 赤池情報量基準
  • BIC : ベイジアン情報量基準
  • deviance : 逸脱度
  • df.residual : 残差の自由度

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など。