RでStataと同じクラスターロバスト標準誤差を取り扱う
RでStataと同じ様にクラスターロバスト標準誤差を取り扱うには
以前、計量経済学の実証コースで、「実証論文を1本選んで自分でデータを取ってきてReplicateせよ（論文と同じ分析結果を得よ）」という課題があり、この課題で苦労したクラスターロバスト標準誤差の取り扱いについてこの記事を書こうと思う。
計量経済学の分析手法を使う研究を行う多くはStataユーザーであると思う（近年、RとPythonユーザーがかなり増えているようだが）。私もB3から計量分析を学び始め、Stata畑で育った人間だ。授業はRで実証分析を学ぶプログラムで、Replicationも当然Rで行い、自分の書いたコードとstargazerなどでまとめた結果を提出せよといったものであった。そこで母語がStataの私が長時間解かす羽目になったの原因が、Rでクラスターロバスト標準誤差を取り扱う方法が分からなかったからである。
ロバスト標準誤差
Stataのデフォルトでは、誤差項の均一分散を仮定して標準誤差が計算されるので、reg y x1 x2, robustとして標準誤差をheteroskedasticity robustにする。これで不均一分散に対して頑健な標準誤差（heteroskedasticity robust standard error）を考慮した回帰分析を行うことができる。Rで同じ分析を行う方法も少しググれば出てくる。
クラスターロバスト標準誤差
Cluster-robust standard error（Cluster-robust ES）をRで推定するにはどうすればいいのだろう？
Cluster-robust standard errorとは、村や地域のクラスター間の誤差項の相関を考慮したSEのことである。例えば、推定された通常の標準誤差が過小推定されていた場合、本来（Cluster-robust ESが正しいとき）は有意でない係数が有意であると判断されることを避けることができる。
Stataでは、reg y x1 x2, robustだけではクラスター間の誤差項の相関を無視してしまうため、reg y x1 x2, vce(cluster id)を指定する必要がある。確か、, cluster(cluster id)でも出来た気が...（記憶力ゴミなのでここはご自分でご確認下さい）。
これと同じ作業を行うのに少し苦労した。Rのestimatrパッケージを使えばStataで計算されたCluster-robust標準誤差と同じ物を得られる。スペルがestimatorではないことに注意したい。面倒な方法としては、estimatr::coeftest::vcovHCでクラスターを指定する方法があるが、出力も中途半端な気がするし、一度推定したものをcoeftest::vcovHCするのが究極的二度手間である（Stataでは一行で書ける）。ここが唯一Rが完全に劣っている点に思えるが、estimatr::lm_robustを使えば一行でクラスタリングするクラスターを指定してCluster-robust ESを計算することが可能になる。
Stataと同じ結果を得る場合には、以下のようにすればいい。
first2_1 <- lm_robust(y ~ x1 + x2, clusters = cluster id, se_type = "stata", data = data)
非常に単純なようだが、estimatr::iv_robustで操作変数法を用いる場合には、se_type = "stata"では微妙にStataのそれと異なるSEが計算される（私はそうだった）。
操作変数法の場合には、
first2_1 <- lm_robust(y ~ x1 + x2 | z + x2, clusters = cluster id, se_type = "CR0", data = data)
と、se_type = "CR0"のように変更すればStataで計算されたCluster-robust ESと同じ結果が得られるはずだ。推定式のzは操作変数、x1が内生変数の場合である。
何故Cluster-robust ESでも通常のlm_robustの場合とiv_robustの場合でse_typeの指定を変更しなければいけないのかは分からない。実際にStataとRの両方を開きながら推定を行ったが、Stataでは論文の分析（Cluster-robust ES）を完全にReplicate出来ているのに、Rではestimatr::iv_robustでse_type = "stata"と指定すると、微妙に数値が異なるSEが計算されたのである。
おそらく日本語と英語でこの内容（StataとRでの推定結果の違い）が書かれている記事はゼロなので、息の長い記事になることを祈る。次はRmarkdownやLaTexでの資料作成の際に使う参考文献用のBibTexの使い方の記事でも書こうと思う。
とりあえずReplicateには成功したので安心しているが、何故計算方法（結果）が異なるのか気になる。ご存知の方いらっしゃればコメント欄にて教えて下さい。
