UPDATE: 2024-12-07 00:53:04.270224
ここでは固定効果モデルについて簡単にまとめておく。下記の書籍を参考にしている。
固定効果モデル式、Rでモデリングする部分が中心となる。固定効果モデルとマルチレベルは似たようなモデルではあるものの、少し異なる部分もあるので、その点はここでは触れない。下記にまとめているので、そちらを参照ください。
固定効果モデルは経済系、マルチレベルは社会心理、生物系とかで使われているイメージ。
固定効果モデル(Fixed Effect Model)はパネルデータで使用される分析モデル。パネルデータは同一の個体を複数の時点で観察したデータのこと。ここでの個体は、人間や地域のことを指す。クロスセクションデータとは異なるので、データの性質を先にまとめておく。パネルデータの変数\(X_{it}\)は、個人\(i\)の時点\(t\)における値を表す。\(X_{it}\)はいくつかのパターンを持つ。
これらを考慮すると、パネルデータに対する回帰モデルは下記のように書ける。誤差項は1-3の各変数に対応する形で存在する。これらの誤差項が変数と相関せず、内生生が存在しないのであれば、回帰モデルでパラメタを正しく推定できるが、普通は相関する。
\[ Y_{it} = \beta_{0} + \overbrace{X_{it}^{(1)}\beta_{1} + \cdots + X_{it}^{(k)}\beta_{k}}^{1} + \overbrace{Z_{i}^{(1)}\alpha_{1} + \cdots + Z_{i}^{(l)}\alpha_{l}}^{2} + \overbrace{W_{t}^{(1)}\gamma_{1} + \cdots + W_{it}^{(m)}\gamma_{m}}^{3} + \overbrace{\epsilon_{it} + e_{i} + u_{t}}^{\varepsilon_{it}} \]
このままだと見通しが良くないので、変数の数を減らして簡略化したモデルを使って、より理解を深めていく。\(\beta_{1}\)を正しく推定するためには\(X_{it}\)と\(\epsilon_{it} + e_{i} + u_{t}\)の内生性の問題を解決する必要がある。
\[ \begin{eqnarray} Y_{it} &=& \beta_{0} + X_{it}\beta_{1} + Z_{i}\alpha + W_{t}\gamma + \epsilon_{it} + e_{i} + u_{t} \\ i &=& 1,\dots,n \\ t &=& 1,\dots,T \end{eqnarray} \]
そもそも\(e_{i},u_{t}\)は何か。
個体を「地域」と考えるのであれば、\(e_{i}\)は地域の固有の特徴、魅力、能力などが該当する。つまり、時点は関係なく、各地域が持っている固有の観測できない要因が、\(e_{i}\)には含まれる。これは普通に考えるのであれば\(X_{it}\)と相関するため、内生性の問題が生じる。そのため、このような変数を細分化して観測できないので、下記のように\(\phi_{i}, \psi_{t}\)でまとめて肩代わりしてもらうことで、\(X_{it}\)の影響\(\beta_{0}\)を推定しようと試みる。これが固定効果モデル。
\(\phi_{i}\)は時点\(t\)について不変の\(i\)固有の効果全てで、\(\psi_{t}\)は時点\(t\)にのみ依存して個人間で共通の効果をまとめたもの。
\[ \begin{eqnarray} \phi_{i} &=& \beta_{0} + Z_{i}\alpha + e_{i} \\ \psi_{t} &=& W_{t}\gamma + u_{t} \end{eqnarray} \]
\(\phi_{i}, \psi_{t}\)に関して、ダミー変数を利用すれば、個人\(i\)や地域\(i\)は数多くあることが一般的で、時点\(t\)も数年分あることが一般的なので、下記のようにシグマでまとめることが出来る。
\[ \begin{eqnarray} Y_{it} &=& \beta_{0} + X_{it}\beta_{1} + \phi_{i} + \psi_{t} + \epsilon_{it} \\ &=& X_{it}\beta_{1} + \sum_{j=1}^{n} I_{ij}\phi_{j} + \sum_{k=1}^{T} \delta_{tk}\psi_{k} + \epsilon_{it} \\ i&=&1,...,n \\ t &=& 1,...,T \\ \end{eqnarray} \]
これが固定効果モデル(2方向固定効果モデル)であり、推定にバイアスがかからないように個体固有の効果や時点固有の効果をモデルに組み込んでいる。このモデルをみるとわかるが、変数としてモデルに組み込めるのは、個人\(i\)と時点\(t\)の両方が変化する変数のみである。なぜなら、それ以外の特徴をもつ変数はすでに固定効果\(\phi_{i}, \psi_{t}\)としてモデルに含まれている。
推定に関しては、最小2乗法を利用すればよい。ただ、within推定と呼ばれる方法も利用される。固定効果を1つ含む下記の固定効果モデルを考える。
\[ Y_{it} = X_{it}\beta_{1} + \sum_{j=1}^{n} I_{ij}\phi_{j} + \epsilon_{it} \]
この式に対して、平均を計算して、
\[ \bar{ Y_{i} } = \frac{1}{T}\sum_{t=1}^{T}Y_{it}, \quad \bar{ X_{i} } = \frac{1}{T}\sum_{t=1}^{T}X_{it}, \quad \bar{ \epsilon_{i} } = \frac{1}{T}\sum_{t=1}^{T}X\epsilon_{it} \]
先程の式から平均を引くことで、
\[ \bar{ Y_{i} } = \bar{ X_{i} }\beta_{1} + \sum_{j=1}^{n} I_{ij}\phi_{j} + \bar{ \epsilon_{i} } \]
固定効果を消去することが出来る。このモデルで\(\beta\)を推定するのと、さきほどの固定効果モデルを推定しても、パラメタの値は同じである。
\[ \begin{eqnarray} Y_{it} - \bar{ Y_{i} } &=& X_{it}\beta_{1} - \bar{ X_{i} }\beta_{1} + \sum_{j=1}^{n} I_{ij}\phi_{j} - \sum_{j=1}^{n} I_{ij}\phi_{j} + \epsilon_{it} - \bar{ \epsilon_{i} } \\ &=& (X_{it} - \bar{ X_{i} })\beta_{1} + (\epsilon_{it} - \bar{ \epsilon_{i} }) \end{eqnarray} \]
Rによる実証分析
回帰分析から因果分析へ第2版のp269で使用されるデータをお借りして固定効果モデルの推定を行う。まずは必要なライブラリを読み込む。書籍では、fixestパッケージを使用しているが、estimatrパッケージでも固定効果モデルは推定できるはずなので、estimatrパッケージも利用する。
このデータには、47都道府県(pref)の1997年から2019年(year)までの失業率(unemp)と自殺死亡率(suicide)が含まれる。このデータを利用して、固定効果モデルを用いることで、失業率が自殺死亡率に影響するかどうかを分析する。
パネルデータなので、下記の通り地域\(i\)を北海道に限定すると、時点\(t\)が数多く記録されている。
prefdata <- readr::read_csv("~/Desktop/prefecture.csv")
prefdata %>%
filter(pref == '北海道') %>%
arrange(pref)## # A tibble: 23 × 4
## pref year suicide unemp
## <chr> <dbl> <dbl> <dbl>
## 1 北海道 1997 19.6 3.7
## 2 北海道 1998 26.7 4.8
## 3 北海道 1999 26.2 5
## 4 北海道 2000 26.6 5.5
## 5 北海道 2001 23.6 5.8
## 6 北海道 2002 24.6 6.1
## 7 北海道 2003 27.1 6.5
## 8 北海道 2004 26.5 5.8
## 9 北海道 2005 27.3 5.3
## 10 北海道 2006 26.4 5.4
## # ℹ 13 more rowsまずは地域の固定効果をモデルに取り込んで、パラメタを推定する。モデルの結果を見る限り、失業率が1%増加すると、人口10万人あたりの自殺者数がおよそ3人増加することになる。
## # A tibble: 1 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 unemp 3.07 0.0862 35.6 1.09e-181within推定でも同様の数値が計算できる。
prefdata %>%
group_by(pref) %>%
mutate(suicidebar = mean(suicide),
unempbar = mean(unemp),
suicide2 = suicide - suicidebar,
unemp2 = unemp - unempbar) %>%
lm(suicide2 ~ -1 + unemp2,
data = .) %>%
tidy() %>%
filter(term == 'unemp2')## # A tibble: 1 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 unemp2 3.07 0.0843 36.4 7.27e-190さらに時点効果を含めたモデルを作成する(2方向固定効果モデル)。パラメタの値が小さくなり、失業率が1%増加すると、人口10万人あたりの自殺者数がおよそ0.77人増加する、という結果となった。
prefdata %>%
lm(suicide ~ -1 + unemp + pref + as.factor(year),
data = .) %>%
tidy() %>%
filter(term == 'unemp')## # A tibble: 1 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 unemp 0.771 0.159 4.85 0.00000142fixestパッケージのfeols関数では|の後に固定効果を記述すれば固定効果モデルとして推定できる。先程の結果と同じくパラメタの値は0.77となったが、std.errorの値が異なることがわかる。
## OLS estimation, Dep. Var.: suicide
## Observations: 1,081
## Fixed-effects: pref: 47, year: 23
## Standard-errors: Clustered (pref)
## Estimate Std. Error t value Pr(>|t|)
## unemp 0.770868 0.298724 2.58054 0.013121 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1.76097 Adj. R2: 0.8622
## Within R2: 0.022751これは、固定効果モデルで推定する際に、同じ個体\(i\)が何度も記録されることにより標準誤差の計算を調整する必要がある。それをfeols関数では考慮してくれている。標準誤差に関しては下記にまとめている。
固定効果を確認したければ、fixef関数を利用する。
## $pref
## 三重 京都 佐賀 兵庫 北海道 千葉 和歌山 埼玉
## 15.45397 14.13019 17.11202 15.45899 17.49214 14.90050 19.09119 15.04125
## 大分 大阪 奈良 宮城 宮崎 富山 山口 山形
## 17.06491 15.84755 13.11925 16.34422 22.15730 20.04691 18.36547 20.79473
## 山梨 岐阜 岡山 岩手 島根 広島 徳島 愛媛
## 18.33060 17.04139 14.04390 24.71817 22.27627 15.54064 13.88266 17.44879
## 愛知 新潟 東京 栃木 沖縄 滋賀 熊本 石川
## 14.07081 23.07379 14.77983 18.14798 15.48147 14.92670 16.94297 15.64029
## 神奈川 福井 福岡 福島 秋田 群馬 茨城 長崎
## 13.67033 16.38198 16.29855 18.89778 27.86661 18.77796 16.90630 16.83410
## 長野 青森 静岡 香川 高知 鳥取 鹿児島
## 17.49636 22.45245 15.31457 14.55124 19.28600 17.42461 18.57975
##
## $year
## 1997 1998 1999 2000 2001 2002 2003
## 0.0000000 5.5833995 4.5654306 3.7648943 2.9584970 3.5622088 5.5267073
## 2004 2005 2006 2007 2008 2009 2010
## 4.3136910 4.7456166 4.7956255 5.2002432 4.6308207 4.2217904 2.8488754
## 2011 2012 2013 2014 2015 2016 2017
## 2.4818250 0.8309494 0.8858371 -0.2012504 -0.9517742 -2.1754777 -2.3757769
## 2018 2019
## -2.6837893 -2.8346302
##
## attr(,"class")
## [1] "fixest.fixef" "list"
## attr(,"references")
## pref year
## 0 1
## attr(,"exponential")
## [1] FALSEestimatrパッケージのlm_robust関数でも固定効果モデルを推定する。先程の結果と同じくパラメタの値は0.77となり、std.errorの値も調整されていることがわかる。
lm_robust(
suicide ~ unemp,
fixed_effects = ~ pref + year, #固定効果と年固定効果
clusters = pref,
data = prefdata,
se_type = "CR0"
)## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## unemp 0.7708683 0.2923652 2.636662 0.01137947 0.1823677 1.359369 46