固定効果モデル(Fixed Effect)を考える。

固定効果モデルとは

固定効果モデルとは、ざっくり書くと、パネルデータに対するOLSの際にパネルデータ内の個人毎に異なる「個人差」のような部分を除去してOLS推定ができるようになるモデル。

固定効果モデルイメージ

これの何が嬉しいかというと、「個人の才能」みたいな定量的に測れない説明変数が関係ありそうなモデルを考えたいとする。その場合、数値化されたデータは手に入らないため、説明変数を入れることができない。しかし、これを入れないと欠落変数バイアスがかかりα, βなどで誤ったOLS推定量を計算してしまう。
このとき、「個人の才能」のように時間経過によって変わらないであろうものであれば固定効果モデルを考えることでこの問題を解決できる。

固定効果モデルのややちゃんとした説明

ちゃんと書くならば、個人i、時間tでのYを考えるにあたり個人i,時間tによって変わる説明変数X, 誤差項vと、個人iのみによって変わる時間不変のF(固定効果）からなる回帰モデル $Y_{i,t} = αX_{i,t} + F_i + v_{i,t}$ を考える。

このとき、固定効果Fは個人iの才能や能力といったような、個人i固有の定量化しづらいものとなる。定量化できないということは説明変数に入れることができない。説明変数に入れない場合は、Fは誤差項に吸収されるので、 $Y_i,t = αX_{i,t} + ( F_i + v_{i,t} ) = αX_{i,t} + ε_{i,t}$ となる。
しかし、このFは多くの場合でXと相関が起きるため、説明変数Xと誤差項 $ε_{i,t} = F_i + v_{i,t}$ で相関が生まれることでOLSをしても推定量に一致性を持たない結果となる（ちなみに、FとXで相関がない場合はFは変量効果と呼ばれ、変量効果モデルという別の問題を考える必要がある）。

このようなときには、一般的にはFを除去することが考えられる。回帰モデル $Y_{i,t} = αX_{i,t} + F_i + v_{i,t}$ の時間平均は $Y'_i = αX'_i + F_i + v'_i$ となる。式からわかるように、 $F_i$ は時間不変なので時間平均も $F_i$ となる。ここから、 $Y_{i,t} - Y'_i = α(X_{i,t} - X'_i) + (v_{i,t} - v'_i) = αX''_{i,t} + v''_{i,t}$ が導出される。このとき、 $X'', v''$ の間に相関はないのでαのOLS推定量は一致性を持つことができる。

また、式からわかるように固定効果Fは、Yに関係ある 特定の 時間不変説明変数ではなく、 すべての 時間不変要素となるので定量化しづらい要素が全て入ることになるので欠落変数バイアスの影響を減らす効果がある（逆に、全てがまとめて消去されるので時間不変の個々の要素の影響を考えることはできない）。

以上はややざっくりとした概要なので、もっとちゃんとした理論的な部分は書籍などを参照してください。ブログだと下記がわかりやすくまとまっていた。

［R］Rでパネルデータ分析：固定効果モデル - 盆暗の学習記録

lfeパッケージ

概要

この固定効果モデルをRで導入するには色々とパッケージがあるようだが、lfeパッケージのfelm関数を使ってみる（ちなみにplmパッケージが一般的らしい

CRAN - Package lfe

基本的な使い方はbase::lm関数と同じだがformulaは|でブロックを区切って以下のような書き方をする。
y ~ x1 + x2 | f1 + f2 | (Q|W ~ x3+x4) | clu1 + clu2

第1ブロックは、推定したい回帰式(lm関数と同じ)。
第2ブロックは、固定効果として除去したい時不変の変数。
第3ブロックは、操作変数の設定。例だとy ~ x1 + x2 + Q + Wという回帰式のうち、操作変数Q,WがQ ~ x3+x4, W ~ x3+x4に従うという意味。
第4ブロックは、クラスターロバスト標準誤差(robust standard error)のためのクラスタリング対象変数

これらを指定し、あとはlm関数とほぼ同様に重みなどを指定します。

サンプルコード

ドキュメントにあるサンプルコードは以下。idとfirmの固定効果を考えるモデル。

oldopts <- options(lfe.threads=1)

## create covariates
x <- rnorm(1000)
x2 <- rnorm(length(x))

## individual and firm
id <- factor(sample(20,length(x),replace=TRUE))
firm <- factor(sample(13,length(x),replace=TRUE))

## effects for them
id.eff <- rnorm(nlevels(id))
firm.eff <- rnorm(nlevels(firm))

## left hand side
u <- rnorm(length(x))
y <- x + 0.5*x2 + id.eff[id] + firm.eff[firm] + u

## estimate and print result
est <- felm(y ~ x+x2| id + firm)
summary(est)

# =>

Call:
   felm(formula = y ~ x + x2 | id + firm) 

Residuals:
    Min      1Q  Median      3Q     Max 
-4.4196 -0.9338 -0.0152  0.9765  4.6023 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)    
x  -0.48704    0.04616  -10.55   <2e-16 ***
x2  0.61937    0.04671   13.26   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.411 on 966 degrees of freedom
Multiple R-squared(full model): 0.4624   Adjusted R-squared: 0.4441 
Multiple R-squared(proj model): 0.2294   Adjusted R-squared: 0.203 
F-statistic(full model):25.18 on 33 and 966 DF, p-value: < 2.2e-16 
F-statistic(proj model): 143.8 on 2 and 966 DF, p-value: < 2.2e-16

以下では、固定効果モデルと操作変数法を組み合わせた固定効果操作変数法をおこなうことで、逆の因果関係(reverse causation)をもたらす時間不変の固定効果と時変の変数を除去した上での推定値を求めている。

# make an example with 'reverse causation'
# Q and W are instrumented by x3 and the factor x4. Report robust s.e.
x3 <- rnorm(length(x))
x4 <- sample(12,length(x),replace=TRUE)
Q <- 0.3*x3 + x + 0.2*x2 + id.eff[id] + 0.3*log(x4) - 0.3*y + rnorm(length(x),sd=0.3)
W <- 0.7*x3 - 2*x + 0.1*x2 - 0.7*id.eff[id] + 0.8*cos(x4) - 0.2*y+ rnorm(length(x),sd=0.6)

# add them to the outcome
y <- y + Q + W
ivest <- felm(y ~ x + x2 | id+firm | (Q|W ~x3+factor(x4)))
summary(ivest,robust=TRUE)

# =>

Call:
   felm(formula = y ~ x + x2 | id + firm | (Q | W ~ x3 + factor(x4))) 

Residuals:
     Min       1Q   Median       3Q      Max 
-2.81299 -0.70052  0.03803  0.65092  2.81178 

Coefficients:
         Estimate Robust s.e t value Pr(>|t|)    
x         1.14222    0.17483   6.533 1.04e-10 ***
x2        0.55946    0.03365  16.626  < 2e-16 ***
`Q(fit)`  0.94039    0.11730   8.017 3.13e-15 ***
`W(fit)`  1.03717    0.04922  21.071  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9975 on 964 degrees of freedom
Multiple R-squared(full model): 0.7319   Adjusted R-squared: 0.7222 
Multiple R-squared(proj model): 0.6157   Adjusted R-squared: 0.6018 
F-statistic(full model, *iid*):82.19 on 35 and 964 DF, p-value: < 2.2e-16 
F-statistic(proj model): 448.1 on 4 and 964 DF, p-value: < 2.2e-16 
F-statistic(endog. vars):587.9 on 2 and 964 DF, p-value: < 2.2e-16