論理の流刑地

地獄の底を、爆笑しながら闊歩する

固定効果モデルでロバスト標準誤差 in R

備忘用メモ。

Rでパネルデータを扱うときに、固定効果モデルを行うときは、
plmパッケージを使うのが一般的だ。

最近の社会科学(経済学、社会学etc)で固定効果モデルを使うとき、
table上で報告されるのは、頑健標準誤差(robust standard error)である場合が多い。

従属変数がy , 独立変数x1 , x2 とする

Stataで

Stataを使って分析するときは、以下のようなコードをかくことになる。

xtreg y x1 x2, fe vce(id  time)

Rで

Rで固定効果モデルを行い、
上のStataコードと同じ標準誤差を得たい場合には、以下のようなコードになる

library(plm)
p_data <-  pdata.frame( test_data , index= c("id", "time") ) #plm用にデータ変換
fe_res <- plm( y ~ x1 + x2 , data = p_data , model="within") #plmで推定

### clusterロバスト標準誤差を手に入れる ###
cluster_n <- length(unique(test_data$id) ) #クラスタ数をみる

library(lmtest)
coeftest( fe_res , vcovHC(fe_res, type="HC1" , cluster = "id" ) * (cluster_n / (cluster_n - 1)) )

・plm()の戻り値に、クラスター変数を指定しつつvcovHC()を噛ませて、係数推定値のcluster-robustな分散共分散行列を得ている。
・coeftest() 関数の引数として指定する際には、clusterの数を取得しておき[クラスター数 / (クラスター数 - 1) ] を、vcov()の戻り値にかけて補正する。

なんでこうする必要があるのか

なぜこのような回りくどい補正が必要であるのかは、以下の記事(の最初の回答)に詳しい。

stata - Clustered standard errors in R using plm (with fixed effects) - Stack Overflow

・Stataでクラスタロバスト標準誤差を求めるとき、標準誤差を補正するためのスケールパラメータとして、
G_Nクラスター数、Kはパラメータ数、Nはケース数とすると、

  C =    \frac{G_N}{G_N - 1} \frac{N-1}{N - K}

を分散共分散行列にかけている。

・しかしRの場合は、スケールパラメータの第2項  \frac{N-1}{N - K}しか使っていないために、
Stataでの出力に合わせるためには、vcovHC()の出力に第1項 \frac{G_N}{G_N - 1}を外からかけてやる必要があるのである。