【R】新しい回帰分析表のパッケージ {modelsummary}
Rで計量経済や統計分析やる時に、結果をきちんと理解しながらモデル作りたいし、できたモデルの結果を書き出すのも間違いなくやりたいですよね。
新しい回帰分析表のパッケージを発見したので、ざっと試してみました。
Rの回帰分析のパッケージ
Rでの回帰分析の表を書き出すのには {stargazer}というパッケージがあって、かなり柔軟にかつ自動的にhtmlやLatexに書き出しができるのですが、作者の都合か、更新が止まってしまってます。更新が止まってしまうと、新しい分析パッケージの結果とかを受け付けられなくなるんですよね。特に計量経済における回帰分析で役立つ{estimatr}の結果を受け付けないのがつらいです。
他のパッケージで書き出すことも可能です。例えば{huxtable}パッケージのhuxreg関数は、かんたんに書き出せるし、Rのコンソールに見やすい結果を出力してくれるので良いです。しかし、まだ発展途上ということで、{stargazer}ほどの柔軟性はないのと、Latexで出力する際に表の位置が真ん中にいかないという問題がありまして、まだ実用的ではないかなという個人的な印象です。
以前なんとかestimatrの出力結果を固定効果付きでhuxregで書き出すという記事も書きました。
今のところは作業中はhuxregで結果見ながら作業して、最終的に出力するときはいろんなワークアラウンドを駆使してstargazerできれいな表を書き出すという形をとっています。
{modelsummary}パッケージ
そんな時に、以下のTweetを見つけました。
Before we get going, let's install version 0.3.0: install.packages("modelsummary") Done? Great! You now have the power to create PDF/LaTeX tables like this one: pic.twitter.com/M1e4bF7cQD
— Vincent Arel-Bundock (@VincentAB) May 27, 2020
新しい統計モデルのサマリー表ですって？？
しかも、スレッドの次のTweetを見てみると、なかなかきれいな表が出力されています。
パッケージの説明はこちら。
さっそくインストール
{modelsummary}がCRANに上がっているというツイートだったので、CRANからダウンロードしましょう。
dplyrの最新版に依存してるみたいなので、tidyverseも一緒に更新しておきましょう。
install.package("modelsummary") install.package("tidyverse")
シンプルに回帰分析の出力
とりあえずmtcarsデータ使って、回帰分析を行って、結果を出力してみます。
結果はmsummary()関数を使います。
reg <- lm(mpg ~ cyl+disp, mtcars) msummary(reg)
この関数を使うと、コンソールではなくViewerにきれいなテーブルが出力されます。
２つ以上のモデルを同時に表示させる場合には、summary()のようにオブジェクトをそのまま並べるのではなく、リストにします。
リストの要素指定時にモデル名をつけることができて、これが列名になるようです。
regs <- list() regs[['Model A']] <- lm(mpg ~ cyl+disp, mtcars) regs[['Model B']] <- lm(mpg ~ cyl+disp+hp, mtcars) msummary(regs)
出力フォーマット
コンソールに表示
オプションでmarkdownと指定すれば、huxregのようにコンソールに表示することも可能です。
msummary(regs, 'markdown')
| |Model A |Model B | |:-----------|:-------|:-------| |(Intercept) |34.661 |34.185 | | |(2.547) |(2.591) | |cyl |-1.587 |-1.227 | | |(0.712) |(0.797) | |disp |-0.021 |-0.019 | | |(0.010) |(0.010) | |hp | |-0.015 | | | |(0.015) | |Num.Obs. |32 |32 | |R2 |0.760 |0.768 | |Adj.R2 |0.743 |0.743 | |AIC |167.1 |168.0 | |BIC |173.0 |175.3 | |Log.Lik. |-79.573 |-79.009 |
出力フォーマットは以下のものに対応しているようです。
| フォーマット | オプション | |
|---|---|---|
| HTML | 'html' | |
| LaTeX | 'latex' | |
| Text / Markdown / ASCII | 'markdown' | |
| RTF (Microsoft Word互換) | 'ファイル名.rtf' | |
| Images: JPG and PNG | 'ファイル名.png' | |
| gt テーブルオブジェクト(R) | 'gt |
指定もかんたんで、例えばLatexコードを出力したいときは
msummary(regs, 'latex')
と打つだけでOK.
ファイルに保存したい場合は
msummary(regs, 'result.tex') # Tex ファイル msummary(regs, 'result.png') # PNGファイル
情報の編集
テーブル出力パッケージとして簡単にきれいな表を出力できることはわかったが、問題はどの程度柔軟に情報が編集できるかである。
見た目
まずテーブルの見た目については、modelsummaryはkableExtraやgtの柔軟性を応用して編集できるようになっているようだ。
私はgtはまだ使ったことがないが、テーブルを画像出力する時に力を発揮するらしい。
Rmarkdown/knitrで出力するケースではknitr::kableのエクステンションであるkableExtraが良い。
HTML用にもLatex用にも編集がしやすく、きれいなテーブルを出力できる。
htmlについてはgtパッケージが、latexについてはhtmlパッケージがデフォルトになっているらしいが、次のようなオプションで変更することも可能らしい。
例えばhtmlのテーブルをgtパッケージの関数で修飾してみると以下のようになる。
msummary(regs, output = 'html')%>% gt::tab_style(style = cell_fill(color = "lightcyan"), # 行に色を付ける locations = cells_body(rows = 3)) %>% gt::tab_style(style = cell_text(weight = "bold"), #文字を太くする locations = cells_body(rows = 3))
回帰分析表の要素
回帰分析表の要素については基本的な編集機能が実装されているようだ。
標準誤差・t統計量・p値・信頼区間
回帰分析表の係数の推定値の下に表示される統計量。経済学では標準誤差を示すのが一般的だが、他の統計量もstatisticオプションを変更することで示すことができる。
msummary(regs, statistic = 'std.error') # 標準誤差 (デフォルト) msummary(regs, statistic = 'p.value') # p値 msummary(regs, statistic = 'statistic') # t統計量 msummary(regs, statistic = 'conf.int', conf_level = .99) # 99% 信頼区間
表のタイトルとメモ
表のタイトルもかんたん。titleオプションで指定します。
メモも簡単にnotesオプションで追加できるし、複数のメモをリストで表示できる。
msummary(regs,
title = "表１：mtcarsの結果", # タイトル
notes = list('メモ1. 南無妙法蓮華経', # メモ
'メモ２. アーメン'))
変数名の変更
変数名の変更は、coef_mapで行う。
予め変更前後の変数名の対応表をベクトルにしておくと便利そうだ。
新しく入るほうが右に来るらしい。(dplyrのrenameと逆なので注意）
# 変数名対応のベクトル
var_nam = c("mpg" = "Mile per gallon", "(Intercept)" = "Constant","cyl" = "Cylinder", "disp" = "Displacement",
"hp" = "Horse Power")
msummary(regs, coef_map = var_nam)
モデルフィットなどの統計量
msummary()ではデフォルトでR二乗などのモデルフィットの統計量が表示されている。
不必要な統計量はgof_omitオプションで非表示にできる。
msummary(regs, gof_omit = "R2|AIC|BIC")
統計量のリストはgof_mapで見ることができる。
> modelsummary::gof_map # A tibble: 14 x 4 raw clean fmt omit <chr> <chr> <chr> <lgl> 1 nobs Num.Obs. %.0f FALSE 2 r.squared R2 %.3f FALSE 3 adj.r.squared Adj.R2 %.3f FALSE 4 AIC AIC %.1f FALSE 5 BIC BIC %.1f FALSE 6 logLik Log.Lik. %.3f FALSE 7 deviance Deviance %.2f TRUE 8 df.residual DF Resid %.0f TRUE 9 df.null DF Null %.0f TRUE 10 sigma Sigma %.3f TRUE 11 statistic Statistics %.3f TRUE 12 p.value p %.3f TRUE 13 df DF %.0f TRUE 14 null.deviance Deviance Null %.2f TRUE
P値の星
stargazerと違って、msummaryではp値の星はデフォルトではつかない。（最近は星つけないジャーナルも多い）
starオプションをTRUEにすると、デフォルト（10%, 5%, 1%)の基準で星をつけてくれる。
もちろんカスタムも可能。"*"を変更することで好きな文字にすることもできる。
msummary(regs, stars = TRUE) # デフォルト, `* p < 0.1, ** p < 0.05, *** p < 0.01` msummary(regs, stars = c("*" = .05, "**" = .01, "***" = .001)) # カスタム
小数点以下の調整
小数点以下の調整はfmtオプションで行う。sprintfのフォーマットをとるので、
例えば小数点第２位までを表示したければ以下のように書く。
msummary(regs, fmt = '%.2f') # 小数点以下第2位まで
huxregと違って、observationの数はきちんと整数になっているのがいいですね。
ちなみに、科学的表記法にする場合はfをeに変える。
> # 例 > sprintf("%.3f", 0.5342) [1] "0.534" > sprintf("%.3e", 0.5342) [1] "5.342e-01
行の追加
例えば固定効果を入れたかどうか(Yes/No)などを書いた行を追加したい場合、stargazerではそれ専用のオプションが合ったが、
msummaryではhuxregのようにカスタムされた行を追加することで対応するらしい。
しかし、huxregより追加は容易である。
add_rowsオプションで追加する。表と同じ行数の文字ベクトルを作成して、listにして追加する。
(2020年９月19日追記) Ver. 0.4.1 より add_rows はデータフレームしか受け付けなくなったので、文字ベクトルではなくデータフレームで情報を与える必要があります。
##row1 <- c('固定効果1', 'No', 'No') # Ver 0.4.1.以前 ##row2 <- c('固定効果2', "No","No") # Ver 0.4.1.以前 # 0.4.1.よりデータフレームで与える # tibble::tribbleで行ごとにデータを与えてtibbleを作れる。 rows <- tribble(~term, ~`Model A`, ~`Model B`, 'FE 1', 'No', 'No', 'FE 2', "No","No") # Sep 20, added. Now the add_rows accepts only data.frame attr(rows, 'position') <- c(15,16) # position attributeでどこに挿入するか指定する。15-16行目に追加 ## msummary(regs, add_rows = list(row1, row2))# Ver 0.4.1.以前 msummary(regs, add_rows = rows)
また、上の表な手作業ではなく、モデルをカスタマイズして固定効果を表示する方法も提案されていますが、別記事でまとめたいと思います。
で、estimatrは使えるのか？
大変使いやすそうだし、broom::tidyの出力に基づいているので柔軟そうだ。
ではestimatr::lm_robustは使えるのか？
結論から言うと、使えます!
regs_rob <- list() regs_rob[["Model A rob"]] <- lm_robust(mpg ~ cyl+disp, se_type = "stata", fixed_effects = ~ vs + am,mtcars) regs_rob[["Model B rob"]] <- lm_robust(mpg ~ cyl+disp+hp, se_type = "stata", fixed_effects = ~ vs + am,mtcars) msummary(regs_rob)
自動で標準誤差のタイプまで表示してくれていますね。
まとめ
- {modelsummary}はシンプルで使いやすい回帰分析表出力パッケージ
- broom::tidyの結果を受け付けるので、最近開発された統計モデルパッケージ等にも対応
- gtやkableExtraを使ってテーブルが修飾できるし、様々なフォーマットへの出力に対応 (html, latex, RTF, text, 画像ファイル)
- 多くのオプションでテーブルの内容もカスタマイズ可能
まだ実際に論文書くのには使ってはいませんが、よさそうなので使ってみます！
使ったコードはここにおいてます。
