Rで回帰分析のモデルを並べる方法

回帰分析のモデル比較をする表をExcelで作らずRでキレイに整形できるという話を聞いたので少し調べてみた。Rのstargazerパッケージ、texregパッケージできれいに整形できるらしい。いずれも、TeX、HTML、テキストの3種類に出力できる能力がある。

下準備として回帰分析の結果を作る。AERパッケージに入っているCPS1985というデータを使う。

library(AER) #パッケージの呼び出し
head(CPS1985) #CPS1985データの行頭のみ表示
?CPS1985 #データについて説明の表示

このデータを使い4つの回帰分析を作り、それぞれlm0~lm3に格納する。

lm0 <- lm(wage ~ age, data=CPS1985) #年齢が賃金を決めるモデル
lm1 <- update(lm0, ~. + gender) # 性別をモデルに加える
lm2 <- update(lm1, ~. + education) # 学歴をモデルに加える
lm3 <- update(lm2, ~. + experience) # 仕事の経験年数をモデルに加える

memiscパッケージ

テキストでまとめる方法。 memiscパッケージの mtable( ) 関数を使う。

library(memisc)
mtable(lm0, lm1, lm2, lm3)
==============================================================
                     lm0        lm1        lm2        lm3     
--------------------------------------------------------------
  (Intercept)      6.167***   6.929***  -4.843***  -1.957     
                  (0.723)    (0.720)    (1.244)    (6.835)    
  age              0.078***   0.085***   0.113***  -0.367     
                  (0.019)    (0.018)    (0.017)    (1.120)    
  gender: female             -2.275***  -2.335***  -2.344***  
                             (0.430)    (0.388)    (0.389)    
  education                              0.827***   1.307     
                                        (0.075)    (1.120)    
  experience                                        0.481     
                                                   (1.121)    
--------------------------------------------------------------
  R-squared            0.0        0.1        0.3        0.3   
  adj. R-squared       0.0        0.1        0.2        0.2   
  sigma                5.1        4.9        4.5        4.5   
  F                   17.2       23.0       59.9       44.9   
  p                    0.0        0.0        0.0        0.0   
  Log-likelihood   -1622.8    -1609.1    -1553.4    -1553.3   
  Deviance         13635.9    12954.1    10514.6    10510.9   
  AIC               3251.6     3226.2     3116.8     3118.6   
  BIC               3264.5     3243.4     3138.2     3144.3   
  N                  534        534        534        534     
==============================================================

きれいにまとめられる。 memiscパッケージは結果をいろいろな形式で出力なできるそうなので、またエントリーを改めたい。

stargazerパッケージ

stargazerパッケージをインストー

install.packages("stargazer")
library(stargazer)

上の回帰分析をまとめると次のようになる。

stargazer(lm0, lm1, lm2, lm3) # TeX形式
stargazer(lm0, lm1, lm2, lm3, type = "html") #HTML形式

デフォルトにもう少し手を入れていこう。

stargazer(lm0, lm1, lm2, lm3, type = "html", 
          title            = "表1 回帰分析のモデル比較", # タイトルを入れる
          covariate.labels = c("年齢", "性別-女性", "教育年数", 
                               "就労経験", "定数"), #独立変数の名前をいれる
          dep.var.labels   = "賃金"
          )

はてなブログは透明の罫線が表示される困った使用なので、画像として張った。下部にHTML出力のリンクをつけたので、生の出力結果がみたい場合はリンク先を見てほしい。

f:id:iDES:20170625155435p:plain HTML出力

texregパッケージ

次に試すのはtexregパッケージ。デフォルトは次のコマンド。

library(texreg)
htmlreg(list(lm0, lm1, lm2, lm3))

タイトルや変数名などを少し加えてみる。

htmlreg(
  list(lm0, lm1, lm2, lm3),
  caption.above = TRUE,
  caption = "表1 回帰分析のモデル比較",
  custom.coef.names = c("定数","年齢", "性別-女性", "教育年数", 
                        "就労経験")
  )

HTML出力

HTMLの出力はtexregパッケージの方がキレイに出力できるようだ。