回帰分析のモデル比較をする表をExcelで作らずRでキレイに整形できるという話を聞いたので少し調べてみた。Rのstargazerパッケージ、texregパッケージできれいに整形できるらしい。いずれも、TeX、HTML、テキストの３種類に出力できる能力がある。
下準備として回帰分析の結果を作る。AERパッケージに入っているCPS1985というデータを使う。
library(AER) #パッケージの呼び出し head(CPS1985) #CPS1985データの行頭のみ表示 ?CPS1985 #データについて説明の表示
このデータを使い４つの回帰分析を作り、それぞれ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 = "表１ 回帰分析のモデル比較", # タイトルを入れる
covariate.labels = c("年齢", "性別－女性", "教育年数",
"就労経験", "定数"), #独立変数の名前をいれる
dep.var.labels = "賃金"
)
はてなブログは透明の罫線が表示される困った使用なので、画像として張った。下部にHTML出力のリンクをつけたので、生の出力結果がみたい場合はリンク先を見てほしい。
texregパッケージ
次に試すのはtexregパッケージ。デフォルトは次のコマンド。
library(texreg) htmlreg(list(lm0, lm1, lm2, lm3))
タイトルや変数名などを少し加えてみる。
htmlreg(
list(lm0, lm1, lm2, lm3),
caption.above = TRUE,
caption = "表１ 回帰分析のモデル比較",
custom.coef.names = c("定数","年齢", "性別－女性", "教育年数",
"就労経験")
)
HTMLの出力はtexregパッケージの方がキレイに出力できるようだ。
追記: textregパッケージについて新しいエントリをしました。(2020/01/09)
