前置き
記述統計について書いた前回に引き続き, 回帰分析の結果をスマートに表にまとめる方法を複数のパッケージを比較しつつ紹介する. 前回に引き続き, xtable, latex, stagazer 関数の他, 新たに texreg についても紹介したい.
経済学系の論文では, 回帰分析の結果は
| 変数1 | 1.23*** |
| (0.012) | |
| 変数2 | 2.34 |
| (1.23) | |
| N | 100 |
のように係数の推定値の下に括弧書きで標準誤差を表記し, 有意な係数にはアスタリスクを添えるというフォーマットが多い (有意水準は大抵は *: p<.1, **: p < .05, ***: p < .01). stargazer パッケージは R-statistics blog の記事で, texreg は R-bloggers の記事で取り上げられているが, この例ではいずれも右寄せにしてるせいで有意水準を表すアスタリスクがついている数値がずれていて, 見栄えが悪い……*1
texreg は
Description texreg converts coefficients, standard errors, significance stars, and goodness-of-fit statis- tics of statistical models into LaTeX tables or HTML tables/MS Word documents or to nicely formatted screen output for the R console for easy model comparison. A list of several models can be combined in a single table. The output is highly customiz- able. New model types can be easily implemented.
とのことなので記述統計表こそ対応していないが, 回帰分析の結果表を LaTeX 以外にも HTML 及び MS ワードに対応した表を出力してくれるらしい*2. ただし, texreg という名前が示すように, 回帰分析の結果だけしか対応していないようだ. 一方, stargazer は, 前回言及したように記述統計も独力で出力してくれるが,
Description Produces LaTeX code and ASCII text for well-formatted tables that hold regression analysis results from several models side-by-side, as well as summary statistics.
なので, 原則 LaTeX 形式しか出力出来ない.
対応パッケージ比較
では次に, それぞれカバーしている手法はどれくらい多いだろうか. 計量経済学でよく用いられる推定手法は, いわゆる最小二乗法による回帰分析の他に, 一般化線形モデルの一種である変量効果モデル・固定効果モデル, あるいはロジスティック回帰やプロビットモデル, さらにそれを拡張した多項ロジット (multinominal logit), 順序ロジット (ordered logit), 入れ子ロジット (nested logit) などの離散選択モデルと呼ばれる手法がある. これらは統計推理の応用なので, 自然科学や機械学習の分野でも用いられているようだが, 操作変数法と呼ばれる独特の, 純粋な統計学からしたら異端に思われるような手法もよく用いられる (GMM, 2SLS, 3SLSなど). もちろん, AR, MA, ARMA といった時系列分析の手法もよく用いられる. また, 最近ではいわゆるサバイバル分析も用いられる.
よって, これらの手法を R で全てカバーしようとした場合, かなりの数のパッケージが必要となる*3.
stargazer 関数の description には
The stargazer command produces LaTeX code and ASCII text for well-formatted tables that hold regression analysis results from several models side-by-side. It can also output summary statis- tics and data frame content. stargazer supports model objects from aftreg (eha), betareg (betareg), binaryChoice (sampleSelection), bj (rms), brglm (brglm), coeftest (lmtest), coxph (survival), coxreg (eha), clm (ordinal), clogit (survival), cph (rms), dynlm (dynlm), ergm (ergm), errorsarlm (spdep), gam (mgcv), gee (gee), glm (stats), Glm (rms), glmer (lme4), gls (nlme), Gls (rms), gmm (gmm), heckit (sampleSelection), hurdle (pscl), ivreg (AER), lagarlm (spdep), lm (stats), lmer (lme4), lmrob (robustbase), lrm (rms), maBina (erer), mclogit (mclogit), mlogit (mlogit), mlreg (eha), multinom (nnet), nlmer (lme4), ols (rms), phreg (eha), plm (plm), pmg (plm), polr (MASS), psm (rms), rem.dyad (relevent), rlm (MASS), rq (quantreg), Rq (rms), probit (sampleSelection), selection (sampleSelection), svyglm (survey), survreg (survival), tobit (AER), weibreg (eha), zeroinfl (pscl), as well as from the implementation of these in zelig. It also supports the following zelig models: "relogit", "cloglog.net", "gamma.net", "probit.net" and "logit.net".
とあるため, ひと通りの手法はカバーしてると言えるだろう (膨大な量なので全てのパッケージを使用した場合のレビューはしない.). texreg についても, 概ね対応具合は同じようだ.
構文
stargazer
前回書き忘れたが, type="text" とすると latex コードでなくプレーンテキストで表示される. 「とりあえず手っ取り早く R 上に見やすい結果表を出したい」というときに使うといいだろう.
前回の記述統計に続いて, 次の引数が有効である.
dep.var.labels引数は被説明変数につけるラベルで, ヘッダーに出力される. 質的変数を尺度変換し被説明変数とするときなどは候補が複数あることが考えられるので, そういう場合に用いられると考えられる.omit.statは表示しない統計量. 詳しくはマニュアルの stargazer_stat_code_list を参照.- 有意水準については,
star.charで記号を,star.cutoffsで有意水準のカットオフポイントを決められる.ただし 1から3つまで. ciは係数の下の括弧内に標準誤差ではなく信頼区間を表すかどうか. ベクトル型なので変数ごとに指定できる. 信頼区間のレベルはci.levelで設定可能 (デフォルトは95%). 信頼区間を表示した場合の出力サンプルは R-statistics blog の記事 の Table 5 などを参考に.
次のサンプルコードはマニュアルに掲載されているものを一部改変したもの. 事務労働者に対する職場環境についてのアンケート結果である attitude データセットを使用している. 詳細はヘルプ参照.
library(stargazer) ## 2通りのOLSモデル (変数選択が異なる) linear.1 <- lm(rating ~ complaints + privileges + learning + raises + critical, data=attitude) linear.2 <- lm(rating ~ complaints + privileges + learning, data=attitude) ## rating を 70 を境に分割し, プロビットモデルを適用した場合. attitude$high.rating <- (attitude$rating > 70) probit.model <- glm(high.rating ~ learning + critical + advance, data=attitude, family = binomial(link = "probit")) ## stargazer で結果表示 stargazer(linear.1, linear.2, probit.model, title="stargazer による回帰分析の結果")
align=F のままなので, 数値の位置がずれて見栄えが悪い. また, F統計量と残差標準誤差は蛇足に感じるし, 表が横に間延びしてしまっている. そこで, これらを省略し, かつ align=T でやり直してみる.
stargazer(linear.1, linear.2, probit.model, title="stargazer による回帰分析の結果 (修正版)", omit.stat=c("f","ser"), align=T)
回帰分析の結果の場合は, align=T にしてもサンプルサイズ (Observations) の行が自動で中央寄せにされる. これはありがたい.
ただ, まだ表が間延びしている. LaTeXコードを確認すると, 変数ごとに1行ぶんのスペースがあることが分かった. これは no.space=T とすることで省ける*4. 行ごとに挿入するスペースの幅は指定できないようだ. 個人的には係数と対応する標準誤差の列を狭め, それ以外を少し空けるのがベストだと思うのだが.
stargazer(linear.1, linear.2, probit.model, title="stargazer による回帰分析の結果 (最終版)", omit.stat=c("f","ser"), align=T, no.space=T)
また, 表に外部から持ってきた数値を掲載することもできる. coef, se, t, p 引数でそれぞれ係数, 標準誤差, t-値, p-値の箇所に別の値を掲載することができる. さらに, apply.coef, apply.se, ... で数値でなく関数を与えることができるようだ.
その他, style="all" とすると t-値や p-値など利用可能な統計量を全て出力する. よって, type="text" と併用することで, 分析中に結果の早見表ができる.
stargazer(linear.1, linear.2, probit.model, out="", style="all", type="text")
その他にも対応する各パッケージの手法に応じた引数が存在する.
texreg
texreg 関数も stargazer と同様に, 複数のモデルの結果を表記できる. ただし, 2つ以上のモデルの場合は, リスト形式でしか受け付けない (stargazer は単に,で並べる方法でも, リストでも良い). なお, texでなくhtml形式で出力したい場合は, htmlreg 関数を用いる. ただし cssを細かく設定することは出来ない. R のコンソールに表示させたい場合は screenreg を用いる. この3つの関数のシンタックスはほぼ同じ.
l引数に, 表に掲載したいモデルをリストの形式で与える.file引数に出力先のファイル名を与える.starsで有意水準のレベルを設定. いくつでもよい.symbol有意水準を表す記号を設定できる. htmlreg の場合はstar.symbolを使う.digits小数の表示桁数.dcolumndcolumn を使用するかどうか.FALSEにすると数値は中央揃えになる.leading.zerostargazer のinitial.zeroと同じ.ci.forceは stargazer のciと同じ. 信頼区間のレベルはci.force.levelcustom.coef.namesで説明変数の名称を変更.custom.model.namesもある. ただし, 説明説明と同じ数のベクトルを与え, 変更しないものはNAとする.reorder.coefで各説明変数の行の並び替え. 同様にreorder.gofで当てはまりの指標の並び替えも可能.omit.coefは表示を省略する説明変数. stargazer のomit同様, 正規表現で評価される.
これらは, custom.*.names, omit.*, reorder.* の順に実行される. よって, omit.* は改名後の名称でマッチさせる必要がある.
決定係数や尤度などの当てはまりの指標の表示を省略したいときは, extract 関数を用いる. とりあえずは screenreg で結果を表示してみる.
library(texreg) # stargazer と同じモデルを利用 models <-list(linear.1, linear.2, probit.model) screenreg(l=models)
とすると, 当てはまりの指標として, 決定係数, 自由度修正済み決定係数, AIC, BIC, 対数尤度, deviance が表示される. もう少し簡潔に, OLSは自由度修正済み決定係数, プロビットは対数尤度のみにしてみる (出力結果省略).
linear.3 <- extract(model=linear.1, include.rsquared=F) linear.4 <- extract(model=linear.2, include.rsquared=F) probit.2 <- extract(model=probit.model, include.aic=F,include.bic=F,include.deviance=F) screenreg(l=list(linear.3, linear.4, probit.2), reorder.gof=c(2,1,3))
しかし, はっきり言って面倒臭いので出力したものをテキストエディタで削ってもいい...
override.coef,override.se, ... は数値に別の任意の数値を上書きする場合に使用する. R-bloggers の記事では通常の最小二乗法に対し clustering-robust な標準誤差を表示させるために override.se を使用していた. stargazer の apply のように関数を与えることは出来ない.centercenter 環境で囲むかどうか (texreg, htmlreg のみ)booktabsbooktabs を使用するかどうか.tabletable 環境で囲うかどうかsideways横倒しにするかどうか.table=Tでないと意味がない.use.packages必要なパッケージを読み込む記述も出力するかどうか. デフォルトではTRUE.float.posはフロート位置の設定. stargazer のfloat.envと同じ.caption\caption{} の設定.caption.above\caption{} を表の上側に持っていくか. デフォルトでは下側.label\label{} の設定.
texreg(l=models,file="texreg.tex", caption="texregを用いた場合", digits=3, booktabs=T, dcolumn=T, center=T, use.packages=F, caption.above=T, custom.model.names=c("OLS 1", "OLS 2", "Probit")) # html への出力 # エンコードを指定できないため, 日本語を使用すると文字化けする可能性があるので注意. htmlreg(l=models,file="htmlreg.html", caption="htmlregを用いた場合", digits=3, booktabs=T, dcolumn=T, center=T, use.packages=F, html.tag=T, head.tag=T, body.tag=T)
texreg, htmlreg それぞれの結果は以下の通り.
| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 11.011 | 11.258 | -7.476* |
| (11.704) | (7.318) | (3.570) | |
| complaints | 0.692*** | 0.682*** | |
| (0.149) | (0.129) | ||
| privileges | -0.104 | -0.103 | |
| (0.135) | (0.129) | ||
| learning | 0.249 | 0.238 | 0.164** |
| (0.160) | (0.139) | (0.053) | |
| raises | -0.033 | ||
| (0.202) | |||
| critical | 0.015 | -0.001 | |
| (0.147) | (0.044) | ||
| advance | -0.062 | ||
| (0.042) | |||
| R2 | 0.715 | 0.715 | |
| Adj. R2 | 0.656 | 0.682 | |
| Num. obs. | 30 | 30 | 30 |
| AIC | 26.175 | ||
| BIC | 31.780 | ||
| Log Likelihood | -9.087 | ||
| Deviance | 18.175 | ||
| ***p < 0.001, **p < 0.01, *p < 0.05 | |||
よって, texreg も stargazer とほぼ同等のオプションを備えている. しかし, LaTeX 形式では dcolumn を使用した場合, サンプルサイズの位置調整がなされない. これはサンプルサイズの桁が大きくなると見栄えが悪くなってくる.
また, texreg には plotreg 関数なるものも存在する. これは係数の (点) 推定値と信頼区間を視覚化する. ただし係数のスケールが違いすぎるとよくわからないので微妙……
xtable, latex を使う場合
lm, glm オブジェクトを合成してマトリックスを自作するところから始まるが, まず, lm, glm オブジェクト等には標準誤差が保存されていない. そこで, summary 関数によって summary.lm オブジェクトや summary.glm オブジェクトを得る必要がある.
linear.1.sum <- summary(linear.1) linear.2.sum <- summary(linear.2) probit.model.sum <- summary(probit.model) # coefficient 要素に推定値, 標準誤差, t-値/z-値, p-値のマトリックスがある. linear.1.sum$coefficients probit.model.sum$coefficients
よって, 最も簡単な方法は, これらのマトリックスを直接 xtable や latex に渡すことである. 当てはまりの指標やサンプルサイズも併記した方が良いため, 実際にはもう少し複雑になる. なにより stargazer や texreg でやってくれる処理を手動でやる意味もないので, これは特殊なパッケージや独自の手法を使う時に限られるだろう.
前後編のまとめ
- stargazer は手軽に体裁の整った表を出力してくれる.
- texreg も stargazer に次いで簡単なシンタックスで表を出力してくれる.
- html や MSワード/LibreOffice にも対応できる.
screenreg,plotregを使えば R での分析中に結果早見表を作成できる.- 一方で, 記述統計の出力は出来ない.
screenregは標準誤差と信頼区間を同時に表示することは出来ない.- dcolumn を使用してもサンプルサイズの表示を位置調整してくれない.
- 標準誤差と信頼区間を同時に表示することは出来ない.
latexは LaTeX 形式の出力について最も豊富なオプションを持つ.format.dfを使えば R での分析中に結果を見やすくまとめられる.- 一方で, 表の元となるマトリックスは自ら作成する必要があるため, 手間がかかる.
- 回帰分析の結果を見やすくまとめるのが面倒.
xtableは最も簡単な操作で表を出力してくれる- 一方で, 最低限のオプションしかない
latex同様, 表は別途作成しなければならない.
参照
- ill-identified diary, Rでの分析結果をLaTeXに出力するパッケージの比較 (前編)
- R-statistics blog, Tailor Your Tables with Stargazer: New Features for LaTeX and Text Output
- R-bloggers, texreg: A package for beautiful and easily customizable LaTeX regression tables from R
- William Revelle, Package ‘psych’ (pdf)
- David B. Dahl, package 'xtable' (pdf)
- Frank E Harrell Jr, Package ‘Hmisc’ (pdf)
- Marek Hlavac, Package ‘stargazer’ (pdf)
- Philip Leifeld, Package 'texreg' (pdf)
- RjpWiki, LaTeX
