Rで計量経済の回帰分析やるならestimatrパッケージが良さそう。

計量経済学の分析をRでやろう！とした時、計量経済ならではのモデルや検定などがあるんですが、ちょっと前はRだとめんどくさくて、やっぱstataだよね～となってた部分があると思うんですが、Rの隆盛に伴っていろいろパッケージが出てくるようになった気がします。

基本的な計量経済学の教科書をRでやろう！という試みがこれやこれで行われているのですが、何しろRの進化は早い。やっと使えるようになったと思ったパッケージが置き換えられてたりとかザラにありますね。基本的な計量経済モデルを推定するのにはAERやsandwichm, plmなどのパッケージをそれぞれ駆使して推定するものだったのですが、包括的で便利そうなパッケージを見つけたのでメモしておきます。

{estimatr} パッケージ
不均一分散での回帰分析
係数仮説検定
- 追記(Sep.2.2019)
固定効果モデル
クラスター分散
操作変数
ggplot2の中で使う
回帰分析の結果のテーブル
参考文献

{estimatr} パッケージ

declaredesign.org

estimatrって名前がいかにもRパッケージっぽいが、ググるとだいたいestimatorに直されてしまう、ちょっと微妙なググラビリティ。

基本的な計量経済学って回帰分析なわけですが、パラメータの同定を重視する計量経済学では推定値がバイアスを持っていないか非常に注意して使うわけです。
一番基本的な線形モデルの推定方法は最小二乗法(OLS)だけど、バイアスのない推定量を得るにはたくさんの仮定を置かないといけない。仮定が満たされない場合はそれに応じていろいろ推定方法が考案されています。

たとえば、均一分散の仮定が満たされていない場合は不均一分散のもとで標準誤差を推定したり、パネルデータで自己相関の問題を解決するために固定効果をいれてみたり、外生性の仮定が満たされていない場合は操作変数を用いるなど。

これらをわりと簡単なコードでやるのに、RだとAERパッケージとかいろいろな関数が用意されていたのですが、estimatr パッケージが一番手軽で包括的、かつtidyverseシリーズとの連携も良さそう。

あと、かなりstataでの結果と同じものにすることを意識していて、標準誤差のタイプを "stata" にすると Stata でやった推定と同じ結果が出てきます。

Rで回帰分析やるならlm()ですが、問題は標準誤差をいじりにくいこと。従来なら coeftest() で vcov を指定するなどで計算できますが、
Stataでお尻に vce(robust) をつけるだけ、みたいな感じで手軽に不均一分散での標準誤差を推定したい場合は lm_robust() がよい。

下の例ではestmatrのチートシートに基づいて、計量経済学の定番教科書であるWooldridgeのデータで回帰分析を行っている。ちなみにWooldridgeは"Introductory Econometrics"の方。"Econometric Analysis of Cross Section and Panel Data" の方ではないです。アメリカでは前者をBaby Wooldridge, 後者をPapa Wooldridgeといって呼び分けたりしてます。
データがこちらでダウンロードも可能。

不均一分散での回帰分析

lm_robust()だけで推定すると、分散共分散行列はデフォルトでHC2で推定される。
引数のse_typeで指定すると、推定方法を変えることができる。
stataと同じにしたい場合はse_type = "stata"でできるようです。

以下は、gpa3のデータでWooldridgeのChapter 8の内容をlmとcoeftestを使った従来の方法と、estimatrパッケージの lm_robust() を使った方法でやっています。結果は同じ。

library(tidyverse)
library(estimatr)
library(lmtest)
library(car)

# lm + coeftest:  "old" way 

reg_lm = lm(cumgpa ~ sat + hsperc + tothrs + female + black + white,
                    data = gpa3, subset = (spring == 1))

coeftest(reg_lm, vcov=hccm) # hccm in car library

# estimatr::lm_robust : "new" way

reg_estm = lm_robust(cumgpa ~ sat + hsperc + tothrs+female + black + white,
                     data = gpa3 %>% filter(spring == 1),
                     se_type = "HC3")

summary(reg_estm)

もとの推定がHC3でやっているので、HC3 にしてます。

不均一分散のみが問題の場合は、系列相関(Autocorrelation)はないという仮定なので、分散共分散行列は対角成分のみ正(それぞれの分散)でそれ以外がゼロになる。
n個のデータから推定した分散共分散行列を $\hat{\Omega}$ とすると、各対角成分が $\omega_1,...,\omega_n$ となる。均一分散の仮定では $\omega_i$ はすべて $\hat{\sigma}^2$ であるが、不均一分散の場合はそれぞれモデル方法によって異なる。そのモデルの方法がそれぞれHC0からHC4まで名付けられている。sandwichパッケージではHC4まで実装されているが estimatr はHC3までしか実装されていないようだ。

se_type	内容
classical	均一分散。lm()と同じ結果。
HC0	$\omega_i = \hat{u}_i^2$ . White(1980)によるモデル。漸近的に一致。
stata	Stataのvce(robust)と同じ結果。
HC1	Stataと同じ結果。 $\omega_i = \frac{n}{n-k}\hat{u}_i^2$ .
HC2	$\omega_i = \frac{1}{1-h_i}\hat{u}_i^2$
HC3	$\omega_i = \frac{1}{(1-h_i)^2}\hat{u}_i^2$
HC4	$\omega_i = \frac{1}{(1-h_i)^{\min\{4, h_{i}/\bar{h} \}}}\hat{u}_i^2$

$h_i$ はハット行列Hの対角成分でXは右辺の説明変数の行列。 $H = X(X'X)^{-1}X'$
HC1からHC3はMacKinnon and White (1985)で不均一分散下での推定で、サンプルサイズが小さいときにできるだけバイアスが少ない方法として提案されていて、HC3が最もパフォーマンスがよいと後の論文で結論づけられているようです。lm_robustはHC2がデフォルトです。ここらへんはsandwichパッケージのvignettesが詳しいが、estimatrでもそれぞれの推定方法に関する理論的な数学をまとめている。

係数仮説検定

Joint null hypothesis testって日本語でなんていうんだっけ？回帰分析で線形の仮説を検定するために、restrictiveなモデルとそうでないモデルを推定してF testを行って比較しますが、それも簡単にできます。 lh_robust() という関数を使います。

例として、上で推定したGPAの式のうち、黒人(black)と白人(white)のダミー変数の係数が両方０である、という仮説を検定してみます。

# Joint test
# black = 0 and white = 0

lh1 = lh_robust(cumgpa ~ sat + hsperc + tothrs + female + black + white,
          data = gpa3 %>% filter(spring == 1),
          linear_hypothesis = c("black = 0","white = 0"),
          se_type = "HC3")

summary(lh1$lh)

結果は以下のようになります。この標準誤差で検定すると、どちらも０である、という帰無仮説を棄却できないという結果になります。

> summary(lh1$lh)
Linear hypothesis test

Hypothesis:
black = 0
white = 0

Model 1: restricted model
Model 2: cumgpa ~ sat + hsperc + tothrs + female + black + white

  Res.Df Df  Chisq Pr(>Chisq)
1    361                     
2    359  2 1.3449     0.5104

追記(Sep.2.2019)

コメントにいただきましたが、version 0.18.0では、上記の結果が再現できず、summary(lh1$lh)がblack=0とwhite=0それぞれのt検定のみを返すようです。
（この記事はversion 0.16.0を使って書きました。）

マニュアルを見るかぎり変更されたわけではないのでバグだと思うのですが、現状はcarパッケージのlinearHypothesisを代用するしかないみたいです。*1
linearHypothesisという関数はlmオブジェクトの他にlm_robustオブジェクトも受け付けるので、lm_robustで推定した場合わざわざlmでやり直さなくても大丈夫です。
lh_robustオブジェクトはlm_robustを含んでいるので、いかのように$でlm_robustを指定することで、joint hypothesis testを実行できます。

> car::linearHypothesis(lh1$lm_robust, c("black = 0", "white = 0"))
Linear hypothesis test

Hypothesis:
black = 0
white = 0

Model 1: restricted model
Model 2: cumgpa ~ sat + hsperc + tothrs + female + black + white

  Res.Df Df  Chisq Pr(>Chisq)
1    361                     
2    359  2 1.3449     0.5104

固定効果モデル

固定効果モデル(Fixed Effect)はパネルデータで観察されない個人や時間の特徴が従属変数に影響を与える場合は、それを捕捉しないと推定誤差に系列相関が生まれてしまう（各時点で同じ個人の効果が効き続けるので）という問題が起こりうるときに使う方法。誤差項を $\varepsilon_{it} = a_{i} + c_{y} + u_{it}$ に分解して、時点間や個人間での相関を $a_{i}$ や $c_{t}$ に求めて推定します。各時点で使ってる人は皆さんご存知だと思うので詳しい説明は省きますが、Rでこれを使う場合、一番簡単なのはlm()の中で、個人や時点などをダミー変数として扱うこと。

fe_lm = lm(lwage ~ married + union + factor(year) * educ + factor(nr),
data = wagepan)

summary(fe_lm)

結果の頭の部分は以下。この方法だとSummaryしたときに特に興味ない固定効果自体も表示されてしまうので、結果が長くなってしまって非常に見にくい。

> summary(fe_lm)

Call:
lm(formula = lwage ~ married + union + factor(year) * educ + 
    factor(nr), data = wagepan)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1521 -0.1256  0.0109  0.1608  1.4834 

Coefficients: (1 not defined because of singularities)
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)            1.4114905  0.4283573   3.295 0.000993 ***
married                0.0548205  0.0184126   2.977 0.002926 ** 
union                  0.0829785  0.0194461   4.267 2.03e-05 ***
factor(year)1981      -0.0224159  0.1458885  -0.154 0.877893    
factor(year)1982      -0.0057613  0.1458558  -0.039 0.968494

plmパッケージを使うと以下のようになる。これは Using R for Introductory Econometricsのコードと同じ。

wagepan = foreign::read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/wagepan.dta")
wagepan.p = pdata.frame(wagepan, index = c("nr","year"))

pdim(wagepan.p)

fe_plm = plm(lwage ~ married + union + factor(year) * educ, data = wagepan.p, model = "within")
summary(fe_plm)

これをlm_robustでやると、fixed_effectsという引数を指定してやることができる。このほうが推定が早いらしい。
なぜか factor(year)*educ の形で指定するとエラーになるので式を分けて書いた。

fe_estm = lm_robust(lwage ~ married + union + educ + factor(year):educ,
          fixed_effects = ~ nr + year,
          se_type = "classical",
          data = wagepan)

結果は

> summary(fe_estm)
1 coefficient  not defined because the design matrix is rank deficient


Call:
lm_robust(formula = lwage ~ married + union + educ + factor(year):educ, 
    data = wagepan, fixed_effects = ~nr + year, se_type = "classical")

Standard error type:  classical 

Coefficients: (1 not defined because the design matrix is rank deficient)
                      Estimate Std. Error t value  Pr(>|t|)   CI Lower CI Upper   DF
married                0.05482    0.01841  2.9773 2.926e-03  0.0187210  0.09092 3799
union                  0.08298    0.01945  4.2671 2.029e-05  0.0448527  0.12110 3799
educ                        NA         NA      NA        NA         NA       NA   NA
educ:factor(year)1981  0.01159    0.01226  0.9448 3.448e-01 -0.0124562  0.03563 3799

lmと同じ係数が推定されている。educがドロップするのはどうしてだっけ？

ただ実際に使ってみると、fixed_effectsで指定すると永遠に終わらない推定が普通にlmの中でfactorを使ってやると、わりとすぐ終わったことがあった。固定効果に使うカテゴリーのかぶり具合で計算がおかしくなることがあるみたいだが理由はまだわかっていない。パッケージの説明には”Do not pass multiple-fixed effects with intersecting groups”とあるので、計算上おかしくなるのだと思う。たくさん固定効果を入れる場合はこの引数を使わずに、factor()で入れていくほうがよさそう。

クラスター分散

Wooldridgeの14章にクラスター分散の方法が載っているのだけど、例がなぜかFirst Differenceのモデルで、lm_robustだとFDは推定できないので、代わりにここに載っている例を使うことにした。

従来の方法だとplmで推定した後、coeftestでvcovHCでクラスターを指定して計算します。

p.df <- read.table("http://www.kellogg.northwestern.edu/faculty/petersen/htm/papers/se/test_data.txt")
names(p.df) <- c("firmid", "year", "x", "y")
p.df <- pdata.frame(p.df, index = c("firmid", "year"), drop.index = F, row.names = T)
head(p.df)

pm1 <- plm(y ~ x, data = p.df, model = "pooling")
coeftest(pm1, vcov = vcovHC(pm1, type = "HC0", cluster = "group", adjust = T))

結果

> coeftest(pm1, vcov = vcovHC(pm1, type = "HC0", cluster = "group", adjust = T))

t test of coefficients:

            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.029680   0.066939  0.4434   0.6575    
x           1.034833   0.050540 20.4755   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

これもlm_robustなら簡単な指定で推定できる。clustersでクラスタリングするグループを指定する。

reg_c_estm = lm_robust(y ~ x, data = p.df,
          clusters = firmid, se_type = "CR0")

lm_robustでclustersを指定すると、デフォルトではse_typeがCR2になるが、上の結果と同じものにするために、CR0で指定した。
結果は以下の通り

> lm_robust(y ~ x, data = p.df,
+           clusters = firmid, se_type = "CR0")
              Estimate Std. Error    t value     Pr(>|t|)   CI Lower  CI Upper  DF
(Intercept) 0.02967972 0.06693896  0.4433848 6.576796e-01 -0.1018372 0.1611967 499
x           1.03483344 0.05054005 20.4755131 4.361347e-68  0.9355359 1.1341310 499

se_typeはCR0, stata, CR2が指定できる。

se_type	内容
CR0	通常のクラスタリング。クラスタ内での相関を認める推定。
stata	Stataと同じ結果。 $\frac{N-1}{N-K}\frac{S}{S-1}$ という項をCR0にかけて、有限サンプル化でバイアスがかからないように補正している
CR2	HC2のアイデアをクラスタリングに延長したもの。 Pustejovsky and Tipton (2018)

詳しくはestimatrの数学ノート参照

操作変数

計量経済学の中で、特に因果効果の同定・推定に使われるツールが操作変数（手段変数, Instrumental Variable, IV）だろう。
従属変数Yと説明変数Xの間に内生性がある場合に、Xと相関のある外生変数Zを使ってXのうち外生的な部分だけで推定するというものである。

詳しい説明は他の文献に譲るとして、WooldridgeのChapter 15のコードは以下のようになる。従来は,RでIVを使うときはAERパッケージのivregを使うのが主流だ。以下の練習データは Card(1995)の有名な論文がもとで、教育年数の賃金に対する効果を推定する際に、教育年数は内生変数である（教育→賃金は成り立つが、教育が、観察されないが賃金に影響を与える変数（先天的な能力など）と相関しているため、推定値にバイアスがかかる）ので、外生的（と仮定される）「大学への近さ」を操作変数とし、大学への近さと大学へ行く年数に相関があることを利用した推定を行っている。

式の指定は y ~ x1 + x2という式で、x1が内生変数、zが操作変数だとすると、 y ~ x1 + x2 | z + x2 というように、| で分けて第1ステージの推定式を書く。Stataと違って第一ステージの式を全て書かないと行けない。

library(AER)

card <- foreign::read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/card.dta")

ols = lm(log(wage) ~ educ + exper + I(exper^2) + black + smsa + south + 
           smsa66 + reg662 + reg663 + reg664 + reg665 + 
           reg666 + reg667 + reg668 + reg669, data = card)

# IV
iv = ivreg(log(wage) ~ educ + exper + I(exper^2) + black + smsa + south + 
           smsa66 + reg662 + reg663 + reg664 + reg665 + 
           reg666 + reg667 + reg668 + reg669 |
             nearc4 + exper + I(exper^2) + black + smsa + south + 
             smsa66 + reg662 + reg663 + reg664 + reg665 + 
             reg666 + reg667 + reg668 + reg669, data = card)

summary(ols)
summary(iv)

これをestimatrで書き換える。 iv_robust() という関数を使う。
式の指定はivregと同様である。

iv_estm = iv_robust(log(wage) ~ educ + exper + I(exper^2) + black + smsa + south + 
                      smsa66 + reg662 + reg663 + reg664 + reg665 + 
                      reg666 + reg667 + reg668 + reg669 |
                      nearc4 + exper + I(exper^2) + black + smsa + south + 
                      smsa66 + reg662 + reg663 + reg664 + reg665 + 
                      reg666 + reg667 + reg668 + reg669,
                    se_type = "classical",data = card,
                    diagnostics = TRUE)

summary(iv_estm)

educの係数はOLSだと0.0746933だが、IVだと ivreg でも iv_robust でも0.131504と推定されている。
OLSとIVの推定値が異なるかどうか統計的に検定するのがHausman test (ハウスマン検定？ホースマン検定？）だが、
これもiv_robustでは diagnostics という引数をTRUEにしてやると、IVに関連する検定を一通り行って報告してくれる。

Diagnostics:
                 numdf dendf  value  p.value    
Weak instruments     1  2994 13.256 0.000276 ***
Wu-Hausman           1  2993  1.168 0.279973    
Overidentifying      0    NA     NA       NA

操作変数の数が内生変数より多いときはOveridentifying testで操作変数が外生かどうか検定できますが、この場合は内生変数１に対して操作変数が１なので検定はできてません。

ggplot2の中で使う

ggplotで描く推定曲線の信頼区間を robust で推定したい場合、ggplot内で stat_function の method で指定可能。

library(ggplot2)

ggplot(mtcars, aes(mpg,hp)) +
  geom_point() +
  stat_smooth(method = "lm",fill = "blue") +
  stat_smooth(method = "lm_robust",fill = "orange") + 
  theme_bw()

f:id:keita43a:20190416042513p:plain — ggplot2 + lm_robust

lmで推定した標準誤差が青い影で、lm_robustで推定したものがオレンジ。少し異なっているのがわかる。

回帰分析の結果のテーブル

texやrmarkdownでそのままhtmlやtexに書き出したりすることが多くなってきているので、新しい分析パッケージの問題はそのオブジェクトが既存のテーブル書き出す系パッケージに受け入れられるかどうか。今の所texregパッケージはestimatrのオブジェクトをそのまま書き出せるようです。xtableとhuxtableはdata.frameを受けるので、推定したオブジェクトを %>% でtidyに通してやると書き出せるようです。

問題はStargazer。きれいな回帰分析のテーブルを書き出すのに便利なパッケージなのですが、しばらく更新が止まっているようです。
そこで、estimatrの結果をstargazerで使うときは、lmで推定したモデルを受け渡し、stargazer 内で改めて標準誤差を指定し直す、というワークアラウンドが提案されています。そのためにstarprepという関数が用意されている。

declaredesign.org

下の例では上で推定した reg_lmを受け渡した後、stargazerの関数内で標準誤差を指定する引数 se に starprep のオブジェクトを受け渡すことで、標準誤差をStataと同じロバスト標準誤差(HC1)にする設定にしている。

library(stargazer)

stargazer(reg_lm,se = starprep(reg_lm, se_type = "stata"),type = "html")


	Dependent variable:

	cumgpa

sat	0.001^***
	(0.0002)

hsperc	-0.009^***
	(0.001)

tothrs	0.003^***
	(0.001)

female	0.303^***
	(0.059)

black	-0.128
	(0.123)

white	-0.059
	(0.115)

Constant	1.470^***
	(0.224)


Observations	366
R²	0.401
Adjusted R²	0.391
Residual Std. Error	0.469 (df = 359)
F Statistic	39.982^*** (df = 6; 359)

Note:	^p<0.1; ^p<0.05; ^**p<0.01

以上が計量経済での回帰分析に便利そうなパッケージ estimatrの紹介でした。

関連して以下の記事もどうぞ。
keita43a.hatenablog.com

keita43a.hatenablog.com

Rでミクロ計量経済学やるのに詳しい教科書が出たようです。（洋書）

Learning Microeconometrics with R (Chapman & Hall/CRC The R Series)

作者:Adams, Christopher P.
発売日: 2020/12/03
メディア: ハードカバー

参考文献

Card, David (1995) “Using geographic variation in college proximity to estimate the return to schooling.” Aspects of Labour Economics: Essays in Honour of John Vanderkamp. University of Toronto Press

MacKinnon JG, White H (1985). “Some Heteroskedasticity-Consistent Covariance Matrix
Estimators with Improved Finite Sample Properties.” Journal of Econometrics, 29, 305–
325. doi:10.1016/0304-4076(85)90158-7.

White H (1980). “A Heteroskedasticity-Consistent Covariance Matrix and a Direct Test for
Heteroskedasticity.” Econometrica, 48, 817–838.

*1:追記：開発側でも指摘されてましたが、今のところはlh_robustの例からjoint testを削除するという対応で、joint testは実装されてない模様。 lh_robust bug · Issue #320 · DeclareDesign/estimatr · GitHub