Rで学ぶ計量経済学と機械学習 8
計量経済学４：固定効果モデル
安藤道人（立教大学） 三田匡能 (株式会社 GA technologies)
2023-06-10
1 本講義の目的
- 固定効果モデルの推定の基礎を学ぶ。
2 パッケージの読み込み
今回の分析で使うパッケージを読み込む。
pacman::p_load(tidyverse,
estimatr, # lm_robustによるパネルデータ分析(メイン)
stargazer, # lmやplmのときの結果表
texreg, # lmやplmやlm_robustのときの結果表
plm, # plmによるパネルデータ分析(参考)
lmtest) # lmやplmでロバスト標準誤差に修正（参考)3 パネルデータ
パネルデータ(panel data)とは、複数の個体を複数の期間にわたって観察したデータである。 （e.g. 同じ人を追跡調査したデータ、国✕年度の公的統計データ）
例えば以下のような形式のもの。
data(Grunfeld) # plmパッケージに入っているデータの読み込み
Grunfeld %>% rename(`企業ID`=firm,
`年度`=year,
`投資額`=inv,
`企業価値`=value,
`資本ストック`=capital) %>%
filter(1951 <= `年度` & `年度` <= 1953) %>%
head() %>% knitr::kable()| 企業ID | 年度 | 投資額 | 企業価値 | 資本ストック | |
|---|---|---|---|---|---|
| 17 | 1 | 1951 | 755.9 | 4833.0 | 1207.7 |
| 18 | 1 | 1952 | 891.2 | 4924.9 | 1430.5 |
| 19 | 1 | 1953 | 1304.4 | 6241.7 | 1777.3 |
| 37 | 2 | 1951 | 588.2 | 2289.5 | 342.1 |
| 38 | 2 | 1952 | 645.5 | 2159.4 | 444.2 |
| 39 | 2 | 1953 | 641.0 | 2031.3 | 623.6 |
\(k\)個の説明変数\(X\)と被説明変数\(Y\)の個体\(i\)・時点\(t\)のパネルデータは、次のように表記される。
\[ (X_{1it}, X_{2it},…,X_{kit}, Y_{it}), \hspace{1em} i = 1,…,n \text{ and } t = 1,…,T \]
- 欠損値の有無による分類がある
- balanced panel：各個体と各時点の観測値がすべて揃っているパネルデータ
- unbalanced panel：欠損のあるパネルデータ
4 固定効果モデル
固定効果モデルはパネルデータの分析で使われる代表的な分析手法のひとつである。とくに、欠落変数バイアスを引き起こす個体固有効果（時間によって変化しない個体に固有の要素）を除去し、個体内の変動(within variation)を活用して回帰分析を行うものがよく使われる。
4.1 one-way固定効果モデル
4.1.1 個体固定効果モデル（グループ内推定：within）
以下のような個体固定効果モデルを考える。
\[ Y_{it} = \beta_0 + \beta_1 X_{it} +\theta_i + \varepsilon_{it} \]
パネルデータを用いることができる場合、以下の３つの方法によって個体固定効果（entity fixed effects）\(\theta_i\)を除去することができる。
(1) “一回の階差モデル（first difference model）”によるOLS推定
\[ (Y_{i,t+1} - Y_{it}) = (\beta_0 - \beta_0) + \beta_1 (X_{i,t+1}-X_{it}) + (\varepsilon_{i,t+1} - \varepsilon_{it}) \]
記号を置き換えて、
\[ \Delta Y_{it} = \beta_1 \Delta X_{it}+ \Delta \varepsilon_{it} \]
- 推定方法：
- 説明変数、被説明変数それぞれ\(t+1\)期から\(t\)期を引く
- 上の式をOLS推定する
(2) “\(n-1\)個のダミー説明変数”を用いたOLS推定
最小二乗ダミー変数推定（Least Squares Dummy Variables (LSDV) 推定）とも呼ばれる。
\[ Y_{it} = \beta_0 + \beta_1 X_{it} + \gamma_2 D2_i + \cdots + \gamma_n Dn_i + \varepsilon_{it} \\ \text{where } D2_i = \begin{cases} 1 & \text{for } i = 2\\ 0 & \text{otherwise} \end{cases} \text{, etc.} \]
- 推定方法：
- 個体ダミー変数（個体\(i\)に該当する場合に1、それ以外は0となるダミー変数）\(D2_i, \cdots, Dn_i\)を作成する
- 上の式をOLS推定する
(3) ”平均差分法（Entity-demeaned）”を用いたOLS推定
\[ \begin{align} \tilde{Y}_{it} &= \beta_1 \tilde{X}_{it} + \tilde{\varepsilon}_{it}, \\ \text{where } \tilde{Y}_{it} &= Y_{it} - \bar{Y}_i, \hspace{1em} \bar{Y}_i = \frac{1}{T} \sum^T_{t=1} Y_{it}\\ \tilde{X}_{it} &= X_{it} - \bar{X}_i, \hspace{1em} \bar{X}_i = \frac{1}{T} \sum^T_{t=1} X_{it}\\ \tilde{\varepsilon}_{it} &= \varepsilon_{it}- \bar{\varepsilon}_i, \hspace{1em} \bar{\varepsilon}_i = \frac{1}{T} \sum_{t=1}^T \varepsilon_{it} \end{align} \]
- 推定方法：
- 説明変数・被説明変数について、変数から期間平均を引く
- 上の式をOLS推定する
- \(n-1\)個の個体ダミー説明変数による推定と同じ推定値が得られる
- 統計ソフトでは通常は平均差分法による推定が行われる
4.1.2 時間固定効果モデル（グループ間推定：between）
なお、時間の固定効果（time fixed effects）\(\pi_t\)を除去したい場合も似たような手法で分析できる。
\[ Y_{it} = \beta_0 + \beta_1 X_{it} +\pi_t + \varepsilon_{it} \]
- \(T-1\)個の時間ダミー説明変数（時間\(t\)に該当する場合に1、それ以外は0となるダミー変数）を用いたOLS推定
- time-demeanedを用いたOLS推定
ただし、計量経済学ではこのような時間固定効果のみを想定したモデルを推定することはほとんどない。
4.2 two-way固定効果モデル
個体と時間の固定効果モデル（entity and time fixed effects model）は、two-way固定効果モデルと呼ばれ、以下のように表される。
\[ Y_{it} =\beta_1 X_{it} +\theta_i + \pi_t + \varepsilon_{it} \]
個体の固定効果\(\theta_i\)と時間の固定効果\(\pi_t\)の両方を除去したい場合は、それぞれの推定方法の組み合わせになる。
- \(n-1\)個の個体ダミー変数と\(T-1\)個の時間ダミー変数を用いたOLS推定
- entity demeaningと\(T-1\)個の時間ダミー変数を用いたOLS推定
- time demeaningと\(n-1\)個の個体ダミー変数を用いたOLS推定
- entity & time demeaningを用いたOLS推定
- 説明変数と被説明変数について、個体と時間両方の平均を引いてOLS推定
なお、パネルデータを活用した計量経済分析では、時間固定効果がないと仮定できるケースはまれであるため、通常はone-way固定効果モデルではなくtwo-way固定効果モデルを用いる。
5 使用するパッケージとデータ
5.1 パッケージとデータ
Rでは、{estimatr}パッケージのlm_robust関数や{plm}パッケージのplmによってパネルデータ分析を行うことができる。
ここでは、ロバスト標準誤差やクラスタロバスト標準誤差を簡単に利用できるestimatr::lm_robustを用いた分析方法を紹介する。plmを用いた推定も参考に残している。
なお、パネルデータ分析については、誤差項の系列相関も考慮する必要があるため、クラスタロバスト標準誤差を用いることが望ましい。
今回は、{plm}パッケージ含まれているGrunfeldデータセットを用いる。
# データ読み込み
data("Grunfeld")
head(Grunfeld)## firm year inv value capital
## 1 1 1935 317.6 3078.5 2.8
## 2 1 1936 391.8 4661.7 52.6
## 3 1 1937 410.6 5387.1 156.9
## 4 1 1938 257.7 2792.2 209.2
## 5 1 1939 330.8 4313.2 203.4
## 6 1 1940 461.2 4643.9 207.2Grunfeldは1935～1954年にかけてのアメリカの10の企業のbalanced
panelデータである。
firm: 企業（ID的なもの）inv: 投資総額value: 企業価値capital: 資本ストック
5.2 データの可視化
# 投資総額(Y)
Grunfeld %>% ggplot(aes(x = year,
y = inv,
colour = as.factor(firm))) +
geom_line()
# 企業価値(X1)
Grunfeld %>% ggplot(aes(x = year,
y = value,
colour = as.factor(firm))) +
geom_line()
# 資本ストック(X2)
Grunfeld %>% ggplot(aes(x = year,
y = capital,
colour = as.factor(firm))) +
geom_line()
# 企業価値(X1)と投資総額(Y)
Grunfeld %>% ggplot(aes(x = value,
y = inv,
colour = as.factor(firm))) +
geom_point()
# 資本ストック(X2)と投資総額(Y)
Grunfeld %>% ggplot(aes(x = capital,
y = inv,
colour = as.factor(firm))) +
geom_point()
# 企業価値(X1)と 資本ストック(X2)
Grunfeld %>% ggplot(aes(x = value,
y = capital,
colour = as.factor(firm))) +
geom_point()
6 lm_robustによる推定
投資総額を企業価値、資本ストック、個体固定効果、時間固定効果に回帰する。
\[
Y_{it} =\beta_1 X1_{it} + \beta_2 X2_{it} +\theta_i + \pi_t +
\varepsilon_{it}
\]
ここでは、estimatr::lm_robustによる推定方法を説明する。
6.1 Pooled OLS
まず、lmおよびlm_robustを用いてプールドOLSモデルを推定する。
# lm
pooled_lm <- lm(inv ~ value + capital,
data = Grunfeld)
# lm_robust
pooled_robust <- lm_robust(inv ~ value + capital,
data = Grunfeld,
se_type = "HC1") # HC1はロバスト標準誤差の種類
#結果表の出力
screenreg(list(pooled_lm, pooled_robust),
custom.model.names = c("Pooled", "Pooled, robust"),
include.ci = FALSE,
digits = 3)##
## ========================================
## Pooled Pooled, robust
## ----------------------------------------
## (Intercept) -42.714 *** -42.714 ***
## (9.512) (11.575)
## value 0.116 *** 0.116 ***
## (0.006) (0.007)
## capital 0.231 *** 0.231 ***
## (0.025) (0.049)
## ----------------------------------------
## R^2 0.812 0.812
## Adj. R^2 0.811 0.811
## Num. obs. 200 200
## RMSE 94.408
## ========================================
## *** p < 0.001; ** p < 0.01; * p < 0.056.2 One-way固定効果モデル
次に、固定効果として個体固定効果のみを含む固定効果モデル(one-way fixed-effect model) を推定する。これには、n-1個の個体ダミー変数を説明変数に導入する方法と、平均差分法(entity-demeaned)を用いて固定効果を除去する方法がある。
なおlmやlm_robustでは、変数Xをfactor(X)にして右辺に加えると、Xのそれぞれの値をカテゴリとみなしてそれぞれのダミー変数を説明変数に加えてくれる。
#n-1個の個体ダミーを直接導入
oneway_fe1 <- lm_robust(inv ~ value + capital + factor(firm),
data = Grunfeld,
clusters = firm, #クラスタのレベルを指定
se_type = "stata") #クラスタSE: Stataと同一のものを使う
#lm_robustの"fixed_effects"オプションを利用
oneway_fe2 <- lm_robust(inv ~ value + capital,
data = Grunfeld,
fixed_effects = ~ firm, #個体固定効果の指定
clusters = firm, #クラスタのレベルを指定
se_type = "stata") # クラスタロバストはStataと同様の方法で計算
#結果表の出力
screenreg(list(pooled_robust, oneway_fe1, oneway_fe2),
custom.model.names = c("Pooled", "Oneway 1", "Oneway 2"),
include.ci = FALSE,
digits = 3)##
## ======================================================
## Pooled Oneway 1 Oneway 2
## ------------------------------------------------------
## (Intercept) -42.714 *** -70.297
## (11.575) (92.355)
## value 0.116 *** 0.110 *** 0.110 ***
## (0.007) (0.016) (0.016)
## capital 0.231 *** 0.310 *** 0.310 ***
## (0.049) (0.054) (0.054)
## factor(firm)2 172.203 **
## (50.343)
## factor(firm)3 -165.275 **
## (46.345)
## factor(firm)4 42.487
## (76.822)
## factor(firm)5 -44.320
## (69.275)
## factor(firm)6 47.135
## (81.614)
## factor(firm)7 3.743
## (77.005)
## factor(firm)8 12.751
## (78.701)
## factor(firm)9 -16.926
## (74.876)
## factor(firm)10 63.729
## (91.047)
## ------------------------------------------------------
## R^2 0.812 0.944 0.944
## Adj. R^2 0.811 0.941 0.941
## Num. obs. 200 200 200
## RMSE 94.408 52.768 52.768
## N Clusters 10 10
## ======================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05二つの固定効果モデル推定においては、どちらの方法を用いても、係数推定値および標準誤差は全く同じとなる。ただし、個体ダミー（factor(firm)）を直接モデルに入れると、定数項や個々の個体ダミー変数の係数値の推定結果が表示される。
6.3 Two-way固定効果モデル
次に、個体固定効果に加えて、時間固定効果を加えた固定効果モデル（twoway fixed-effect model)を推定する。
#n-1個の個体ダミーとT-1個の年ダミーを直接導入
twoway_fe1 <- lm_robust(inv ~ value + capital + factor(firm) + factor(year) ,
data = Grunfeld,
clusters = firm, #クラスタのレベルを指定
se_type = "stata") # クラスタロバストはStataと同様の方法で計算
#lm_robustの"fixed_effects"オプションを利用
twoway_fe2 <- lm_robust(inv ~ value + capital,
data = Grunfeld,
fixed_effects = ~ firm + year, #個体固定効果、時間固定効果の指定
clusters = firm, #クラスタのレベルを指定
se_type = "stata") # クラスタロバストはStataと同様の方法で計算
#結果表の出力
screenreg(list(pooled_robust, oneway_fe1, oneway_fe2, twoway_fe1, twoway_fe2),
custom.model.names = c("Pooled", "Oneway 1", "Oneway 2", "Twoway 1", "Twoway 2"),
include.ci = FALSE,
digits = 3)##
## ===================================================================================
## Pooled Oneway 1 Oneway 2 Twoway 1 Twoway 2
## -----------------------------------------------------------------------------------
## (Intercept) -42.714 *** -70.297 -86.900
## (11.575) (92.355) (62.010)
## value 0.116 *** 0.110 *** 0.110 *** 0.118 *** 0.118 ***
## (0.007) (0.016) (0.016) (0.011) (0.011)
## capital 0.231 *** 0.310 *** 0.310 *** 0.358 *** 0.358 ***
## (0.049) (0.054) (0.054) (0.049) (0.049)
## factor(firm)2 172.203 ** 207.054 ***
## (50.343) (39.027)
## factor(firm)3 -165.275 ** -135.231 **
## (46.345) (35.132)
## factor(firm)4 42.487 95.354
## (76.822) (59.423)
## factor(firm)5 -44.320 -5.439
## (69.275) (50.629)
## factor(firm)6 47.135 102.889
## (81.614) (62.961)
## factor(firm)7 3.743 51.467
## (77.005) (57.614)
## factor(firm)8 12.751 67.491
## (78.701) (61.119)
## factor(firm)9 -16.926 30.218
## (74.876) (56.268)
## factor(firm)10 63.729 126.837
## (91.047) (70.615)
## factor(year)1936 -19.197
## (21.243)
## factor(year)1937 -40.690
## (34.158)
## factor(year)1938 -39.226 *
## (16.150)
## factor(year)1939 -69.470 *
## (27.708)
## factor(year)1940 -44.235 *
## (17.829)
## factor(year)1941 -18.804
## (18.317)
## factor(year)1942 -21.140
## (14.537)
## factor(year)1943 -42.978 **
## (12.874)
## factor(year)1944 -43.099 **
## (11.285)
## factor(year)1945 -55.683 **
## (15.601)
## factor(year)1946 -31.169
## (21.467)
## factor(year)1947 -39.392
## (27.132)
## factor(year)1948 -43.717
## (39.900)
## factor(year)1949 -73.495
## (39.260)
## factor(year)1950 -75.896
## (37.766)
## factor(year)1951 -62.481
## (50.717)
## factor(year)1952 -64.632
## (52.917)
## factor(year)1953 -67.718
## (44.894)
## factor(year)1954 -93.526 *
## (32.560)
## -----------------------------------------------------------------------------------
## R^2 0.812 0.944 0.944 0.952 0.952
## Adj. R^2 0.811 0.941 0.941 0.943 0.943
## Num. obs. 200 200 200 200 200
## RMSE 94.408 52.768 52.768 51.725 51.725
## N Clusters 10 10 10 10
## ===================================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05ここでも、どちらの方法を用いても係数推定値および標準誤差は全く同じとなる。しかし、個体ダミー（factor(firm))および時間ダミー(factor(year))を用いた推定の場合、定数項や個々の個体ダミーや年ダミーの推定結果が表示される。
7
（参考）plmによる推定
ここでは、plmを用いて固定効果モデルを推定する方法を紹介する。ただし、標準誤差はロバスト標準誤差やクラスタロバスト標準誤差ではない点に注意が必要である。そのため、plmを用いた場合の標準誤差の修正方法も紹介する。
7.1 パネルデータフレームへの転換
まず、{plm}パッケージのpdata.frame()を使って、個体と時間を有するパネル・データフレームとする。
ただし、今回使うGrunfeldデータ(plmパッケージに所収）は最初からパネルデータフレームとなっているため、下記のコードは説明のための便宜的な処置である。
original_df <- Grunfeld # いったんdfとして読み込む（下記のpdata.frameの説明のための便宜的な処置）
panel_df <- pdata.frame(original_df, index = c("firm", "year")) # 個体=firm、時間=yearと認識させたdfここで、firmとyearが因子（factor)に変換されている点に注意が必要である。
7.2 Pooled OLS
まず、lm関数あるいはplm関数を使って、プールドOLSモデルを推定する。
# lmによるpooled OLS
pool_lm <- lm(inv ~ value + capital, data = panel_df)
# plmによるpooled OLS
pool_plm <- plm(inv ~ value + capital, data = panel_df,
model = "pooling")
#stargazer
stargazer(pool_lm, pool_plm,
type = "text",
digits = 3,
df = FALSE,
column.labels = c("Pooled lm", "Pooled plm"),
model.names = F,
model.numbers = F)##
## ================================================
## Dependent variable:
## ----------------------------
## inv
## Pooled lm Pooled plm
## ------------------------------------------------
## value 0.116*** 0.116***
## (0.006) (0.006)
##
## capital 0.231*** 0.231***
## (0.025) (0.025)
##
## Constant -42.714*** -42.714***
## (9.512) (9.512)
##
## ------------------------------------------------
## Observations 200 200
## R2 0.812 0.812
## Adjusted R2 0.811 0.811
## Residual Std. Error 94.408
## F Statistic 426.576*** 426.576***
## ================================================
## Note: *p<0.1; **p<0.05; ***p<0.017.3 One-way固定効果モデル
次に、固定効果として個体固定効果のみを含むone-way固定効果モデルを推定する。
# lmによる固定効果モデル推定(個体固定効果のみ）
# panel_dfにおけるfirmはすでに因子に変換されているので、
# そのままでも個体ダミー変数となるが、誤解のないようにfactor(firm)と書く。
fe_oneway_lm <- lm(inv ~ value + capital + factor(firm),
data = panel_df)
# plmによる固定効果モデル推定（個体固定効果のみ）
# yearも因子なので、factor(year)と書くと、ダミー変数として認識される。
fe_oneway_plm <- plm(inv ~ value + capital,
data = panel_df,
model = "within") # within=固定効果モデル
#stargazer
stargazer(pool_lm, pool_plm, fe_oneway_lm, fe_oneway_plm,
type = "text",
digits = 3,
df = FALSE,
column.labels = c("Pooled lm", "Pooled plm", "Oneway lm", "Oneway plm"),
model.names = F,
model.numbers = F)##
## ================================================================
## Dependent variable:
## --------------------------------------------
## inv
## Pooled lm Pooled plm Oneway lm Oneway plm
## ----------------------------------------------------------------
## value 0.116*** 0.116*** 0.110*** 0.110***
## (0.006) (0.006) (0.012) (0.012)
##
## capital 0.231*** 0.231*** 0.310*** 0.310***
## (0.025) (0.025) (0.017) (0.017)
##
## factor(firm)2 172.203***
## (31.161)
##
## factor(firm)3 -165.275***
## (31.776)
##
## factor(firm)4 42.487
## (43.910)
##
## factor(firm)5 -44.320
## (50.492)
##
## factor(firm)6 47.135
## (46.811)
##
## factor(firm)7 3.743
## (50.565)
##
## factor(firm)8 12.751
## (44.053)
##
## factor(firm)9 -16.926
## (48.453)
##
## factor(firm)10 63.729
## (50.330)
##
## Constant -42.714*** -42.714*** -70.297
## (9.512) (9.512) (49.708)
##
## ----------------------------------------------------------------
## Observations 200 200 200 200
## R2 0.812 0.812 0.944 0.767
## Adjusted R2 0.811 0.811 0.941 0.753
## Residual Std. Error 94.408 52.768
## F Statistic 426.576*** 426.576*** 288.500*** 309.014***
## ================================================================
## Note: *p<0.1; **p<0.05; ***p<0.017.4 Two-way固定効果モデル
最後に、時間固定効果も含むtwo-way固定効果モデルを推定する。
# lm関数による固定効果モデル推定(個体固定効果と時間固定効果)
# panel_dfにおけるfirmとyearはすでに因子に変換されているので、
# そのままでも個体ダミー変数となるが、誤解のないようにfactor(firm), factor(year)と書く。
fe_twoway_lm <- lm(inv ~ value + capital + factor(firm) + factor(year),
data = panel_df)
# plmの固定効果モデル推定（個体固定効果と時間固定効果）
fe_twoway_plm <- plm(inv ~ value + capital,
data = panel_df,
model = "within",
effect = "twoways") # within=固定効果モデル
#stargazer
stargazer(fe_oneway_lm,
fe_oneway_plm,
fe_twoway_lm,
fe_twoway_plm,
type = "text",
digits = 3,
df = FALSE,
column.labels = c("Oneway lm",
"Oneway plm",
"Twoway lm",
"Twoway plm"),
model.names = F, model.numbers = F)##
## =================================================================
## Dependent variable:
## ---------------------------------------------
## inv
## Oneway lm Oneway plm Twoway lm Twoway plm
## -----------------------------------------------------------------
## value 0.110*** 0.110*** 0.118*** 0.118***
## (0.012) (0.012) (0.014) (0.014)
##
## capital 0.310*** 0.310*** 0.358*** 0.358***
## (0.017) (0.017) (0.023) (0.023)
##
## factor(firm)2 172.203*** 207.054***
## (31.161) (35.173)
##
## factor(firm)3 -165.275*** -135.231***
## (31.776) (35.709)
##
## factor(firm)4 42.487 95.354*
## (43.910) (50.722)
##
## factor(firm)5 -44.320 -5.439
## (50.492) (57.831)
##
## factor(firm)6 47.135 102.889*
## (46.811) (54.174)
##
## factor(firm)7 3.743 51.467
## (50.565) (58.179)
##
## factor(firm)8 12.751 67.491
## (44.053) (50.971)
##
## factor(firm)9 -16.926 30.218
## (48.453) (55.723)
##
## factor(firm)10 63.729 126.837**
## (50.330) (58.525)
##
## factor(year)1936 -19.197
## (23.676)
##
## factor(year)1937 -40.690
## (24.695)
##
## factor(year)1938 -39.226*
## (23.236)
##
## factor(year)1939 -69.470***
## (23.656)
##
## factor(year)1940 -44.235*
## (23.810)
##
## factor(year)1941 -18.804
## (23.694)
##
## factor(year)1942 -21.140
## (23.382)
##
## factor(year)1943 -42.978*
## (23.553)
##
## factor(year)1944 -43.099*
## (23.610)
##
## factor(year)1945 -55.683**
## (23.896)
##
## factor(year)1946 -31.169
## (24.116)
##
## factor(year)1947 -39.392*
## (23.784)
##
## factor(year)1948 -43.717*
## (23.970)
##
## factor(year)1949 -73.495***
## (24.183)
##
## factor(year)1950 -75.896***
## (24.346)
##
## factor(year)1951 -62.481**
## (24.864)
##
## factor(year)1952 -64.632**
## (25.350)
##
## factor(year)1953 -67.718**
## (26.611)
##
## factor(year)1954 -93.526***
## (27.108)
##
## Constant -70.297 -86.900
## (49.708) (56.047)
##
## -----------------------------------------------------------------
## Observations 200 200 200 200
## R2 0.944 0.767 0.952 0.720
## Adjusted R2 0.941 0.753 0.943 0.670
## Residual Std. Error 52.768 51.725
## F Statistic 288.500*** 309.014*** 110.983*** 217.442***
## =================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01lmを用いてもplmを用いてもvalueおよびcapitalの係数推定値および標準誤差は全く同じとなる。しかし、lmと個体ダミー・時間ダミーを用いた推定の場合、個体固定効果や時間固定効果の値が推定される。それぞれの個体固定効果はfirm
1との差を、それぞれの時間固定効果は1935年からの差を推定している。
なお、twoway固定効果は、plm関数のoneway推定と時間ダミーを用いて、以下のように推定することもできる。
# plmの固定効果モデル推定2（個体固定効果と時間固定効果）
fe_twoway_plm2 <- plm(inv ~ value + capital + factor(year),
data = panel_df,
model = "within") # within =固定効果モデル
#stargazerで3つのtwoway FEを比べる。
stargazer(fe_twoway_lm, fe_twoway_plm, fe_twoway_plm2,
type = "text",
digits = 3,
df = FALSE,
column.labels = c("Twoway lm", "Twoway plm", "Twoway plm2"),
model.names = F,
model.numbers = F)##
## ======================================================
## Dependent variable:
## ----------------------------------
## inv
## Twoway lm Twoway plm Twoway plm2
## ------------------------------------------------------
## value 0.118*** 0.118*** 0.118***
## (0.014) (0.014) (0.014)
##
## capital 0.358*** 0.358*** 0.358***
## (0.023) (0.023) (0.023)
##
## factor(firm)2 207.054***
## (35.173)
##
## factor(firm)3 -135.231***
## (35.709)
##
## factor(firm)4 95.354*
## (50.722)
##
## factor(firm)5 -5.439
## (57.831)
##
## factor(firm)6 102.889*
## (54.174)
##
## factor(firm)7 51.467
## (58.179)
##
## factor(firm)8 67.491
## (50.971)
##
## factor(firm)9 30.218
## (55.723)
##
## factor(firm)10 126.837**
## (58.525)
##
## factor(year)1936 -19.197 -19.197
## (23.676) (23.676)
##
## factor(year)1937 -40.690 -40.690
## (24.695) (24.695)
##
## factor(year)1938 -39.226* -39.226*
## (23.236) (23.236)
##
## factor(year)1939 -69.470*** -69.470***
## (23.656) (23.656)
##
## factor(year)1940 -44.235* -44.235*
## (23.810) (23.810)
##
## factor(year)1941 -18.804 -18.804
## (23.694) (23.694)
##
## factor(year)1942 -21.140 -21.140
## (23.382) (23.382)
##
## factor(year)1943 -42.978* -42.978*
## (23.553) (23.553)
##
## factor(year)1944 -43.099* -43.099*
## (23.610) (23.610)
##
## factor(year)1945 -55.683** -55.683**
## (23.896) (23.896)
##
## factor(year)1946 -31.169 -31.169
## (24.116) (24.116)
##
## factor(year)1947 -39.392* -39.392*
## (23.784) (23.784)
##
## factor(year)1948 -43.717* -43.717*
## (23.970) (23.970)
##
## factor(year)1949 -73.495*** -73.495***
## (24.183) (24.183)
##
## factor(year)1950 -75.896*** -75.896***
## (24.346) (24.346)
##
## factor(year)1951 -62.481** -62.481**
## (24.864) (24.864)
##
## factor(year)1952 -64.632** -64.632**
## (25.350) (25.350)
##
## factor(year)1953 -67.718** -67.718**
## (26.611) (26.611)
##
## factor(year)1954 -93.526*** -93.526***
## (27.108) (27.108)
##
## Constant -86.900
## (56.047)
##
## ------------------------------------------------------
## Observations 200 200 200
## R2 0.952 0.720 0.799
## Adjusted R2 0.943 0.670 0.763
## Residual Std. Error 51.725
## F Statistic 110.983*** 217.442*** 31.899***
## ======================================================
## Note: *p<0.1; **p<0.05; ***p<0.017.5
coeftestによる（クラスタ）ロバスト標準誤差
plmによるOLS推定における標準誤差も、デフォルトでは均一分散を仮定している。したがって、ここでもロバスト標準誤差を用いる必要がある。さらに、パネルデータ分析については、誤差項の系列相関も考慮する必要があるため、クラスタロバスト標準誤差を用いることが望ましい。
上記の重回帰モデルの推定結果の標準誤差を、（クラスタ）ロバスト標準誤差に置き換える。
ここでは、lmtestおよびplmパッケージに含まれるcoeftest()を用いる。
# coeftest()を用いてロバスト標準誤差を計算する。
# lmで推定したものをロバスト標準誤差にする。
fe_twoway_lm_robust <- coeftest(fe_twoway_lm,
vcovHC(fe_twoway_lm, method = "arellano", type = "HC1"))
#plmで推定したものをクラスタロバスト標準誤差にする。
fe_twoway_plm_clustered <- coeftest(fe_twoway_plm,
vcovHC(fe_twoway_plm, method = "arellano", type = "HC1"))
# stargazerで比べる。
stargazer(fe_twoway_lm, fe_twoway_lm_robust, fe_twoway_plm, fe_twoway_plm_clustered,
type = "text",
digits = 3,
df = FALSE,
column.labels = c("lm", "lm,robust", "plm", "plm,clustered"),
model.names = F,
model.numbers = F)##
## ====================================================================
## Dependent variable:
## ------------------------------------------------
## inv inv
## lm lm,robust plm plm,clustered
## --------------------------------------------------------------------
## value 0.118*** 0.118*** 0.118*** 0.118***
## (0.014) (0.019) (0.014) (0.010)
##
## capital 0.358*** 0.358*** 0.358*** 0.358***
## (0.023) (0.054) (0.023) (0.043)
##
## factor(firm)2 207.054*** 207.054***
## (35.173) (45.879)
##
## factor(firm)3 -135.231*** -135.231***
## (35.709) (43.642)
##
## factor(firm)4 95.354* 95.354
## (50.722) (67.374)
##
## factor(firm)5 -5.439 -5.439
## (57.831) (71.355)
##
## factor(firm)6 102.889* 102.889
## (54.174) (72.398)
##
## factor(firm)7 51.467 51.467
## (58.179) (74.159)
##
## factor(firm)8 67.491 67.491
## (50.971) (68.425)
##
## factor(firm)9 30.218 30.218
## (55.723) (70.870)
##
## factor(firm)10 126.837** 126.837
## (58.525) (80.741)
##
## factor(year)1936 -19.197 -19.197
## (23.676) (20.420)
##
## factor(year)1937 -40.690 -40.690
## (24.695) (25.228)
##
## factor(year)1938 -39.226* -39.226*
## (23.236) (22.904)
##
## factor(year)1939 -69.470*** -69.470**
## (23.656) (29.377)
##
## factor(year)1940 -44.235* -44.235**
## (23.810) (20.424)
##
## factor(year)1941 -18.804 -18.804
## (23.694) (20.063)
##
## factor(year)1942 -21.140 -21.140
## (23.382) (20.795)
##
## factor(year)1943 -42.978* -42.978**
## (23.553) (21.009)
##
## factor(year)1944 -43.099* -43.099*
## (23.610) (24.369)
##
## factor(year)1945 -55.683** -55.683**
## (23.896) (23.219)
##
## factor(year)1946 -31.169 -31.169
## (24.116) (23.811)
##
## factor(year)1947 -39.392* -39.392*
## (23.784) (22.847)
##
## factor(year)1948 -43.717* -43.717
## (23.970) (27.348)
##
## factor(year)1949 -73.495*** -73.495***
## (24.183) (26.067)
##
## factor(year)1950 -75.896*** -75.896***
## (24.346) (26.746)
##
## factor(year)1951 -62.481** -62.481*
## (24.864) (32.884)
##
## factor(year)1952 -64.632** -64.632*
## (25.350) (37.156)
##
## factor(year)1953 -67.718** -67.718*
## (26.611) (38.983)
##
## factor(year)1954 -93.526*** -93.526***
## (27.108) (29.316)
##
## Constant -86.900 -86.900
## (56.047) (75.436)
##
## --------------------------------------------------------------------
## Observations 200 200
## R2 0.952 0.720
## Adjusted R2 0.943 0.670
## Residual Std. Error 51.725
## F Statistic 110.983*** 217.442***
## ====================================================================
## Note: *p<0.1; **p<0.05; ***p<0.018 参考文献
安藤道人「計量経済学２」あるいは「計量経済特論２」の「パネルデータ分析」の講義資料
Chapter 16 Fixed Effects (Nick Huntington-Klein (2021)“The Effect: An Introuction to Research Design and Causality”)
Chapter 8 Panel Data (Scott Cunningham (2021) “Causal Inference: The Mixtape)