2012/01/23 (Mon) R, Stata: Robust standard errors

Stataで回帰分析するにはregコマンドを使う。

reg y x1 x2 

デフォルトでは誤差項のhomoskedasticityを仮定して標準誤差が計算されるので、これをheteroskedasticity robustにするためには、robustオプションをつけて推計する。

reg y x1 x2, robust



これをRでどうやるかというと、まずlm関数がStataのregコマンドに対応する。

lm(y ~ x1 + x2, data = dat)

Rでは、Stataと違い複数のデータを同時に保持できるので（これがRのStataに対する大きなアドバンテージの1つ)、引数dataに使用するデータを指定する。

lm()関数はlmオブジェクトを返すので、変数に代入することができる。

obj <- lm(y ~ x1 + x2, data = dat)
summary(obj)

summary関数はlmオブジェクトに用いると、係数と標準誤差の表を表示する。仮説検定などもobjを通じて行う。

さて、heteroskedasticity robust 標準誤差の計算するには、以下のようにやる。

V <- sandwich(obj)
diag(V)^(1/2)

場合によっては、library(sandwich)としてsandwichパッケージを先に読み込む必要があるかもしれない。なにをやっているかというと、もとの回帰分析の結果を用いて、分散・共分散行列だけを計算し直しているわけだ（これはsandwich関数の仕事)。で、係数の標準誤差はその行列の対角要素の平方根ということ。

はじめは少し面食らうのだけど、要するに、Rでは推計を一からやり直すのではなく、その他の情報はそのままに、標準誤差だけを計算し直す。StataとRでのこの発想の違いに気づくまで、しばらく時間がかかった。

別に計算時間が節約できるわけでもないから、特にユーザーとしてのメリットはない。しかし、想像するに、これはオープンソースならではの仕組みではないだろうか。すでにlm関数があるところへ、誰かがその標準誤差だけを計算し直す関数を書き、公開する。別の修正方法が提案されればまた他の誰かが別の関数を書く・・・。この方法なら、もとのlm関数を書き直す必要がないから、多くの人が開発に参加しやすそうだ。


ちなみに、デメリットもないわけではない。Robust cluster standard errors というのがあるが、クラスターの定義がlmオブジェクトに含まれていないので、計算するには元データからクラスター定義を引っ張ってくる必要がある。
J Hardenさんがこれを計算する関数を書いているが（jjharden.web.unc.edu/archives/344)、この関数では欠損値がある場合にデータのサイズとクラスター定義のミスマッチが生じてエラーが生じる可能性がある。そもそもlmオブジェクトとは別に元データを引数に入れる時点で今ひとつスマートではない。どうもひと工夫必要そうである。