Rで学ぶ計量経済学と機械学習 8
計量経済学４：固定効果モデル
安藤道人（立教大学） 三田匡能 (株式会社 GA technologies)
2025-10-22
1 本講義の目的
- 固定効果モデルの推定の基礎を学ぶ。
2 パッケージの読み込み
今回の分析で使うパッケージを読み込む。
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データセットを用いる。
## 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 データの可視化
# 資本ストック(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:TWFE)を推定する。
#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 fixestによる推定
lm_robustだと、固定効果の数が3つ以上になる場合など、推定速度が著しく遅くなることがある。そこで、はやい速度で固定効果モデルを推定できるfixestの`feols’を紹介する。ここではTWFE(twoway
fixed-effect model)のみを実施する。
twoway_fe2_fixest <- fixest::feols(inv ~ value + capital | firm + year, # firmとyearが固定効果
data = Grunfeld,
cluster = ~ firm) #firmレベルのクラスタロバストSE
screenreg(list(twoway_fe2, twoway_fe2_fixest),
custom.model.names = c("TWFE, lm_robust", "TWFE, fixest"),
include.ci = FALSE,
digits = 3)##
## ====================================================
## TWFE, lm_robust TWFE, fixest
## ----------------------------------------------------
## value 0.118 *** 0.118 ***
## (0.011) (0.011)
## capital 0.358 *** 0.358 ***
## (0.049) (0.048)
## ----------------------------------------------------
## R^2 0.952
## Adj. R^2 0.943
## Num. obs. 200 200
## RMSE 51.725
## N Clusters 10
## Num. groups: firm 10
## Num. groups: year 20
## R^2 (full model) 0.952
## R^2 (proj model) 0.720
## Adj. R^2 (full model) 0.943
## Adj. R^2 (proj model) 0.717
## ====================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05Adj. R^2 (full model)は、固定効果も含めた全体の説明力、Adj. R^2 (proj model)は、固定効果を除去したwithinモデルにおける説明変数の説明力を示す。
8 参考文献
安藤道人「計量経済学２」あるいは「計量経済特論２」の「パネルデータ分析」の講義資料
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)