備忘用メモ。
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でクラスターロバスト標準誤差を求めるとき、標準誤差を補正するためのスケールパラメータとして、
はクラスター数、Kはパラメータ数、Nはケース数とすると、
を分散共分散行列にかけている。
・しかしRの場合は、スケールパラメータの第2項しか使っていないために、
Stataでの出力に合わせるためには、vcovHC()の出力に第1項を外からかけてやる必要があるのである。