【R】固定効果モデルの推定がめっちゃ速いパッケージ { fixest }
以前から話題になっていたが、最近のアップデートで操作変数法にも対応したということで使ってみたら、めっちゃ早かったのでシェア。
Rで固定効果モデル(経済学でいう固定効果モデル)を推定する場合、方法はいろいろあるが、基本的なものはlm()にダミー変数をfactorで入れたり、
パネルデータ分析のパッケージであるplmパッケージのplm()関数を使うものがある。
estimatrパッケージでも固定効果を入れて推定することが可能である。
keita43a.hatenablog.com
ただ、問題はサンプルサイズが大きくなった場合と固定効果の変数の数が増えた場合に、推定に著しく時間がかかることである。
estimatrでも固定効果を3ついれたりすると、収束しなかったりすることがある。
そんな大きなデータセットにおける問題を解決すべく登場したのがlfeパッケージのfelm()関数なのであった。
しかし、それを速さで超えてくるパッケージがfixestなのである。
まずリンク先の最初の図を見てほしい。ベンチマークで速さを比較したものであるが、lfeだけでなくjuliaのFixedEffectModelやStataのreghdfeにも勝っている。*1
fixestの特徴
fixestパッケージの特徴としては大きく2点あると思う。
- めっちゃ速い。
- OLS以外にGLMの関数も用意されている。
めっちゃ速い
本家のグラフで一目瞭然だが、自分でも比較してみた。
以下のサンプルデータ(n = 1000000, 固定効果カテゴリ4つ)で推定した結果を見ていただこう。
# パッケージ library(tidyverse) library(lfe) library(fixest) set.seed(3) # データセットアップ n <- 10^6 data1 <- tibble::tibble(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n), x4 = rnorm(n), id = factor(sample(20,n,replace=TRUE)), firm = factor(sample(13,n,replace=TRUE)), year = factor(sample(2010:2020,n,replace=TRUE)), group = factor(sample(15,n,replace=TRUE)) ) ## 固定効果 id.eff <- rnorm(nlevels(data1$id)) firm.eff <- rnorm(nlevels(data1$firm)) year.eff <- rnorm(nlevels(data1$year)) group.eff <- rnorm(nlevels(data1$group)) ## 従属変数 data1 <- data1 %>% mutate(y = 1*x1 + 0.5*x2 + 2*x3 + 3*x4 + id.eff[id] + firm.eff[firm] + year.eff[year] + group.eff[group] + rnorm(n))
## 推定と比較 mbm <- microbenchmark::microbenchmark( felm = felm(y ~ x1 + x2 + x3 + x4 | id + firm + year + group, data1), fixest = feols(y ~ x1 + x2 + x3 + x4 | id + firm + year + group, data1) ) autoplot(mbm)
結果の比較がこちらである。
データが簡素なのでどちらも速いが、それでもfixestが速い。
GLMの関数もある
OLSだけではなく一般化線形モデル(GLM)のポワソン回帰やロジスティック回帰も使える。
| 線形回帰(OLS) | feols() |
| ポワソン回帰 | feglm(, family = "poisson") |
| 負の二項回帰 | fenegbin() |
| ロジスティック回帰 | feglm(, family = logit") |
使い方
基本的な使い方は、丁寧なIntroductionが英語で用意されているのでここでは簡単な紹介に留める。
1. 基本の関数
基本的な式の書き方は | で変数と固定効果を分ける書き方である。
model1 = feols(y ~ x1 + x2 + x3 + x4 | id + firm + year + group, data1)
GLMはデフォルトのファミリーがポワソンになっているので、ポワソン回帰を行う場合は指定は不要である。
data(trade) # sample data model2 <- feglm(Euros ~ log(dist_km)|Origin+Destination+Product+Year, trade)
2. クラスター標準誤差
クラスター標準誤差の計算もしっかりオプションが用意されている。
例えば、固定効果の最初の2つのカテゴリについてクラスター標準誤差を計算したい場合はsummary関数の中で se = "twoway"とする。
summary(model2, se = "twoway")
このseの指定は"cluster"で1つのカテゴリ、"threeway"で３つ、"fourway"で４つの指定が可能である。
また、特定のカテゴリでクラスター標準誤差を計算したい場合はclusterオプションを使う
summary(model2, cluster = "Product") summary(model2, cluster = ~Product) #同じ結果
さらにロバスト標準誤差（HC１）はse = "hetero"で計算できるので、固定効果を使わないモデルのOLSとしても使うことができる。
3. 固定効果の推定値
summary関数では固定効果を除いた変数の係数推定値が表示されるが、固定効果自体に興味がある場合はfixef()関数で取り出せる。
fixed_effect <- fixef(model1) summary(fixed_effect) fixed_effect$year # 特定のカテゴリの取り出し
またplot関数で、もっとも差の大きな値の比較が簡単に可能になっている。(すべての固定効果推定値が表示されてないことに留意）
plot(fixed_effect)
この他にも、カテゴリ同士の結合(e.g. 月-年効果）、グループごとに変わる係数の推定（e.g. 男女で異なるβ) 、パネルデータ作成による簡単なラグ変数の作成などの機能も盛り込まれている。
テーブル(表)への出力
論文やレポート作成において、推定した結果をエクスポートして活用するのは、手間が省ける以上に再現可能性やミス回避のために重要だ。
以前紹介の記事を書いたmodelsummaryだが、公式では関数を書くことで、固定効果を表示するテーブルを出力することが可能となっている。
https://vincentarelbundock.github.io/modelsummary/articles/newmodels.html#fixed-effects-regression-1vincentarelbundock.github.io
しかし、fixestはパッケージの中にetableというテーブルへの出力関数が含まれている！
いたれりつくせりだな！！！
model2 <- feglm(Euros ~ log(dist_km)|Origin+Destination+Product+Year, trade) summary(model2)
> etable(model1_noFE, model1, cluster = "firm") model1_noFE model1 (Intercept) -0.4849 (0.3075) x1 1.006*** (0.001831) 1.002*** (0.000965) x2 0.5029*** (0.002008) 0.5017*** (0.001129) x3 1.998*** (0.001486) 2*** (0.00103) x4 2.997*** (0.002595) 3*** (0.000737) Fixed-Effects: -------------------- -------------------- id No Yes firm No Yes year No Yes group No Yes ___________________ ____________________ ____________________ Observations 1,000,000 1,000,000 S.E. type: Clustered by: firm by: firm R2 0.74832 0.94743 Within R2 -- 0.93442
etable(model1_noFE, model1, cluster = "firm", tex = TRUE)
もちろんカスタマイズが可能になっている。
くわしい設定方法としても、Vignetteが準備されている。
Exporting estimation tables
*1:2020年2月時点のデータであることに留意。他のパッケージにアップデートが入った可能性はある。
