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

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

MENU

【R】estimatr::lm_robustで回帰分析のテーブル表をhuxregで出す。固定効果付きで。

以前、RでStata的な回帰分析を手軽にやるにはestimatr::lm_robustが良さそう、という記事を書きました。

keita43a.hatenablog.com


実際に自分でもlm_robustを結構使っていたのですが、このlm_robustにもちょっと不便な点がいくつかあって、その一つが回帰分析の表なんですよね。自動で出力できない。

いや、できるんですよ。xtableやtexregなんかと使うとできるんですが、やっぱりstargazerが一番きれいかつ、細かいカスタマイズが簡単で気に入ってるんです。しかし、stargazer側がアップデートされないのでできない。estimatrパッケージは悪くないんです。

estimatr側としては、stargazer用に、starprepというlmオブジェクトからロバストな標準誤差を計算してstargazer側に渡す関数を作るなど涙ぐましい努力をしています。この使い方も上述の記事に書きました。

しかし、固定効果の変数が多かったりして、lmオブジェクトにfactor()で多くのダミー変数を入れたモデルで推定したりするとstarprepで標準誤差が計算できなかったりします。

なんとか簡単に、楽にlm_robustの結果を回帰分析テーブルにしたい。そこで見つけたのがhuxtableパッケージです。

回帰分析用関数huxreg()

この部分はhuxregの説明で、基本的に次の説明書に基づいています。
Regression Tables with huxreg

huxtableはRのテーブル表作成パッケージです。結構便利なパッケージで、latex, html, MS Word, MS Excel, MS PowerPoint, Markdownなどにエクスポートもできるし、Rのコンソール上できれいなテーブルを表示できるので作業中も見やすい。tidyverseとも相性が良いので、テーブルの操作もdplyr使ってできます。

huxreg()を使ってみる

シンプルに、lmオブジェクトなどを引数に入れることで、テーブルを出してくれます。

# ggplot2のダイヤモンドのデータを使う
data(diamonds, package = 'ggplot2')


# 任意の回帰分析を走らせる。
lm1 <- lm(price ~ carat + depth, diamonds)
lm2 <- lm(price ~ depth + factor(color, ordered = FALSE), diamonds)
lm3 <- lm(log(price) ~ carat + depth, diamonds)
> huxreg(lm1, lm2, lm3)
──────────────────────────────────────────────────────────────────────────
                            (1)               (2)              (3)        
                    ──────────────────────────────────────────────────────
  (Intercept)            4045.333 ***      6491.466 ***        7.313 ***  
                         (286.205)         (730.537)          (0.074)     
  carat                  7765.141 ***                          1.971 ***  
                          (14.009)                            (0.004)     
  depth                  -102.165 ***       -53.835 ***       -0.018 ***  
                           (4.635)          (11.815)          (0.001)     
  factor(color,                             -95.142                       
  ordered = FALSE)E                                                       
                                            (62.037)                      
  factor(color,                             554.742 ***                   
  ordered = FALSE)F                                                       
                                            (62.374)                      
  factor(color,                             832.357 ***                   
  ordered = FALSE)G                                                       
                                            (60.338)                      
  factor(color,                            1324.183 ***                   
  ordered = FALSE)H                                                       
                                            (64.296)                      
  factor(color,                            1929.902 ***                   
  ordered = FALSE)I                                                       
                                            (71.561)                      
  factor(color,                            2164.044 ***                   
  ordered = FALSE)J                                                       
                                            (88.144)                      
                    ──────────────────────────────────────────────────────
  N                     53940             53940            53940          
  R2                        0.851             0.032            0.847      
  logLik              -472488.441       -522908.139       -26617.649      
  AIC                  944984.882       1045834.277        53243.298      
──────────────────────────────────────────────────────────────────────────
  *** p < 0.001; ** p < 0.01; * p < 0.05.                                 

Column names: names, model1, model2, model3

このテーブルがそのままコンソールに出てくれるので、見やすい。
作業中にモデル同士を比較したいときなどに間違いが減りそうです。

簡単なカスタマイズ

もちろんいろんなカスタマイズができます。すべて紹介するのは多すぎるので、元のドキュメントを見てほしいのですが、例えばモデル名を変えるときは次のようにする。

huxreg('Price' = lm1, 'Log price' = lm3)

モデル名が表示されている。

──────────────────────────────────────────────────
                     Price          Log price     
              ────────────────────────────────────
  (Intercept)      4045.333 ***        7.313 ***  
                   (286.205)          (0.074)     
  carat            7765.141 ***        1.971 ***  
                    (14.009)          (0.004)     
  depth            -102.165 ***       -0.018 ***  
                     (4.635)          (0.001)     
              ────────────────────────────────────
  N               53940            53940          
  R2                  0.851            0.847      
  logLik        -472488.441       -26617.649      
  AIC            944984.882        53243.298      
──────────────────────────────────────────────────
  *** p < 0.001; ** p < 0.01; * p < 0.05.         

Column names: names, Price, Log price

その他、引数でカスタムできるもの。

  • 統計的有意を示すスター(*)の基準の変更や非表示 (star)
#少し緩い基準
huxreg(lm1, lm3, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01)) 
# スターを表示しない
huxreg(lm1, lm3, stars = NULL) 
  • 標準誤差が示されているかっこ内の内容を変える(p値や信頼区間など) (error_format)
# 標準誤差の代わりにp値を表示
huxreg(lm1, lm3, error_format = '({p.value})')
# 標準誤差の代わりに信頼区間を表示
huxreg(lm1, lm3, ci_level = .99, error_format = '{conf.low} to {conf.high}')
  • 表の下にメモ(ノート)の追加を行う。(note)
# {stars}を入れると有意を示すスターの基準を表示する。
huxreg(lm1, lm3, note = 'かっこ内は標準誤差を表示. {stars}.')
  • 小数点以下の桁数を調整する。(number_format)
huxreg(lm1, lm3, number_format = 2)

モデルごとの統計量の表示

デフォルトではサンプルサイズ、R二乗値、対数尤度、AICが表示されるが、これを調整できる。表示できる統計量は、broomパッケージのglanceという関数を対象オブジェクトに適用したときに表示されるものだ。

たとえば

> broom::glance(lm1)
# A tibble: 1 x 11
  r.squared adj.r.squared sigma statistic p.value    df  logLik    AIC    BIC
      <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>   <dbl>  <dbl>  <dbl>
1     0.851         0.851 1542.   153635.       0     3 -4.72e5 9.45e5 9.45e5
# … with 2 more variables: deviance <dbl>, df.residual <int>

というふうに、各統計量が一覧で表示される。statisticはF統計量である。

たとえばF統計量と、F検定のp値を対数尤度とAICの代わりに表示したいときは、次のようにする。

> huxreg(lm1, lm3, statistics = c('# observations' = 'nobs', 'R squared' = 'r.squared', 'F statistic' = 'statistic','P value' = 'p.value'))
────────────────────────────────────────────────────
                        (1)              (2)        
                 ───────────────────────────────────
  (Intercept)        4045.333 ***        7.313 ***  
                     (286.205)          (0.074)     
  carat              7765.141 ***        1.971 ***  
                      (14.009)          (0.004)     
  depth              -102.165 ***       -0.018 ***  
                       (4.635)          (0.001)     
                 ───────────────────────────────────
  # observations    53940            53940          
  R squared             0.851            0.847      
  F statistic      153634.765       149771.327      
  P value               0.000            0.000      
────────────────────────────────────────────────────
  *** p < 0.001; ** p < 0.01; * p < 0.05.           

Column names: names, model1, model2

たとえばこれを、lm_robustでやると、標準誤差のタイプも表示される。

lm4 <- lm_robust(price ~ carat + depth,diamonds)
lm5 <- lm_robust(price ~ carat + depth,se_type = "stata",diamonds)
>broom::glance(lm5)
  r.squared adj.r.squared statistic p.value df.residual     N se_type
1 0.8506755     0.8506699  48782.02       0       53937 53940     HC1

HC1とHC2が数字扱いで、小数点以下が表示されてしまっているのが微妙。。。

> huxreg(lm4, lm5, statistics = c('# observations' = 'nobs', 'R squared' = 'r.squared', 'SE Type' = 'se_type'))
──────────────────────────────────────────────────
                        (1)             (2)       
                 ─────────────────────────────────
  (Intercept)       4045.333 ***    4045.333 ***  
                    (369.246)       (369.176)     
  carat             7765.141 ***    7765.141 ***  
                     (25.109)        (25.105)     
  depth             -102.165 ***    -102.165 ***  
                      (5.947)         (5.946)     
                 ─────────────────────────────────
  # observations   53940           53940          
  R squared            0.851           0.851      
  SE Type            HC2.000         HC1.000      
──────────────────────────────────────────────────
  *** p < 0.001; ** p < 0.01; * p < 0.05.         

Column names: names, model1, model2

latexなどにエクスポート

huxregで作られるオブジェクトはhuxtableオブジェクトなのだが、これを付属の関数でlatexや様々なフォーマットに変換できる。

> regtab1 = huxreg('Price' = lm1, 'Log price' = lm3)
> print_latex(regtab1)

# ここにLatexのコードが出力される。ここでは省略

> # HTML
> print_html(regtab1)

# ここにhtmlのコードが出力される。ここでは省略

Latexを使う人が注意しないといけないのは、huxtableが出力するlatexコードはいくつかの追加パッケージが必要なので、そのまま使うとエラーが出る。必要なパッケージは次の関数で確認できる。出力は現時点での私のPCでの出力結果。

> report_latex_dependencies()

\usepackage{array}
\usepackage{caption}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{colortbl}
\usepackage{multirow}
\usepackage{hhline}
\usepackage{calc}
\usepackage{tabularx}
\usepackage{threeparttable}
\usepackage{wrapfig}
% These are LaTeX packages. You can install them using your LaTex management software,
% or by running `huxtable::install_latex_dependencies()` from within R.
% Other packages may be required if you use non-standard tabulars (e.g. tabulary).

さらに便利なのは、そのままファイル出力ができる関数である。
作ったテーブルを一つのPDF, MSワード、エクセルなどなどに出力する簡単な関数が用意されている。*1
以下のquick_pdfのpdfをdocx,html, xlsx,pptx,rtf,latexに変えるとそれぞれ対応するファイルに出力される。

quick_pdf(regtab1, file = "regression_table_1.pdf")

固定効果の表示は?

使ってみると、非常に使いやすいパッケージである。
しかし、問題はstargazerに比べてhuxregの引数の柔軟性が少ないことだろう。

例えば、固定効果モデルを推定したら、どの固定効果を含めたか表示したい時がある。
しかし、lm_robustオブジェクトをナイーブに含めると以下のように表示される。

#固定効果なしのモデル
lm4 <- lm_robust(price ~ carat + depth,diamonds)

#固定効果を入れたモデル
lm6 <- lm_robust(price ~ carat + depth, fixed_effects = color,diamonds)

> huxreg("No FEs"=lm4,"Color FE" = lm6)
───────────────────────────────────────────────
                   No FEs         Color FE     
              ─────────────────────────────────
  (Intercept)    4045.333 ***                  
                 (369.246)                     
  carat          7765.141 ***    8070.639 ***  
                  (25.109)        (26.386)     
  depth          -102.165 ***     -89.761 ***  
                   (5.947)         (5.780)     
              ─────────────────────────────────
  N             53940           53940          
  R2                0.851           0.865      
───────────────────────────────────────────────
  *** p < 0.001; ** p < 0.01; * p <            
  0.05.                                        

Column names: names, No FEs, Color FE
Warning message:
In huxreg(`No FEs` = lm4, `Color FE` = lm6) :
  Unrecognized statistics: logLik, AIC
Try setting `statistics` explicitly in the call to `huxreg()`

固定効果付きのモデルの方は、切片は表示されていないが、特に固定効果の変数であるcolorについては何も表示されない。今の所、これを表示するような引数はhuxregには見当たらない。

自力で行を追加する

しかし、huxregが含まれるパッケージhuxtableには柔軟な関数が用意されていて、huxtableオブジェクトに自由に行や列を追加できる。そこで、まずhux()関数で、追加したい行を作ってから、add_rows()で追加する方法を考えてみた。

# 追加したい行の作成。set_bottom_borderで該当セルの下に線を引いている。
FE_row = hux("Color FE","No","Yes") %>% set_bottom_border(1,2:3,TRUE)

# 上の図で、表の7行目(depthの標準誤差が表示されている行)の下に追加したいので、after = 7とする。

add_rows(regtab2,FE_row, after = 7)

このようにすると、思った通りの表ができた。

───────────────────────────────────────────────
                   No FEs         Color FE     
              ─────────────────────────────────
  (Intercept)    4045.333 ***                  
                 (369.246)                     
  carat          7765.141 ***    8070.639 ***  
                  (25.109)        (26.386)     
  depth          -102.165 ***     -89.761 ***  
                   (5.947)         (5.780)     
              ─────────────────────────────────
  Color FE      No              Yes            
              ─────────────────────────────────
  N             53940           53940          
  R2                0.851           0.865      
───────────────────────────────────────────────
  *** p < 0.001; ** p < 0.01; * p <            
  0.05.                                        

Column names: names, No FEs, Color FE

いちいち自分で作らないといけないのが若干面倒だが、わりと簡単に追加できることがわかった。
この方法なら、固定効果に限らず、自由に行を追加できるのでカスタマイズ性は高いですね。

huxtableパッケージ自体は非常に協力かつ柔軟で、テーブルの装飾なども簡単にできるし、latexのエクスポートも柔軟で、さらにdplyrパッケージと互換性があるのが強いので、おすすめパッケージである。他の表作成パッケージとの比較もされているので参照されたい。

*1:MSワードやパワポならflextable, エクセルならopenexlsxパッケージが別途必要