不均一分散に対して頑健な標準誤差をRでやってみる


## 記事の背景
日本語でアクセスできる不均一分散の対処を実装した説明があんまりない。  
自分の復習を兼ねて、均一分散の仮定が満たされない場合の問題と、標準誤差の修正による対処法をRによる実装つきでまとめてみた。
メインはRによる実装なので、不均一分散に関する詳しい説明は別の方の記事などを読んでください。

## OLSの仮定
最小二乗法による推定には、一般的に次のような仮定がおかれる場合が多い。

1. 誤差項の期待値が0
2. 説明変数と誤差項が無相関
3. 均一分散

以上が満たされる場合、推定量は不偏性と一致性、有効性を持つ（Gauss-Markovの定理）。
今回は、このうち均一分散の仮定が満たされない場合について取り扱う。

## 何が問題か
均一分散の仮定が満たされなかったとしても、推定量は依然として、不偏性と一致性を有する。そのため、不均一分散が起きていようとも、計算される推定量の値は変化しない。

しかし、計算される標準誤差の値にはバイアスが生じてしまう。

これはそれなりに重要な問題で、t値は推定量/標準誤差で算出されるので、誤った標準誤差によって実施される仮説検定は信頼することができない。
つまり、本当は有意でないにも関わらず、有意であるという結果が得られる危険性（あるいはその逆）がある。こうなると、的外れな意思決定を行ってしまうリスクが高まってしまう。

## 対処法
標準誤差が正しい値にならないので、そのバイアスを修正する考え方をとる。具体的には、標準誤差として、「不均一分散に対して頑健な標準誤差」(heteroskedasticity robust standard error)を使用すればよい[^1]。単にロバスト標準誤差ともいうが、ここでは以下HRSEとする[^2]。

## どちらの標準誤差を使うべきか

不均一分散かどうかを識別する仮説検定はいくつか存在する。なので、これらの検定を実行することで、どちらを使うべきか判断してもよい。
しかし、マーケティングデータのような人間行動が介在するデータの多くは、均一分散の仮定は現実的ではない。
したがって、一般的に社会科学やビジネスで使うようなデータについてはHRSEを常に使えばよい[^3]。

## Rによる実装
多くの統計ソフトがそうであるように、Rもlm()関数を使うと、デフォルトで均一分散の仮定がおかれてしまう。別に、sandwichパッケージをインストールする必要がある。

まずは、lm()関数による実行結果を書く。
例で用いるデータ等はWooldridgeの教科書を用いた。

```r
library("wooldride")

# wooldridgeのgpa3データセットを利用
data("gpa3")

# GPAをSATスコア、勉強時間、性別、人種で回帰
lm<-lm(formula = cumgpa ~  sat + hsperc + tothrs + female + black + white, data = gpa3)
summary(lm)
```

実行結果

```
Call:
lm(formula = cumgpa ~ sat + hsperc + tothrs + female + black + 
    white, data = gpa3)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.04543 -0.41107  0.03558  0.46175  2.96070 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.8791381  0.2978285   2.952 0.003261 ** 
sat          0.0008721  0.0002330   3.743 0.000196 ***
hsperc      -0.0055988  0.0016165  -3.463 0.000565 ***
tothrs       0.0120630  0.0009322  12.941  < 2e-16 ***
female       0.1666743  0.0770153   2.164 0.030777 *  
black       -0.0260697  0.1922434  -0.136 0.892169    
white        0.0142993  0.1839229   0.078 0.938052    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8659 on 725 degrees of freedom
Multiple R-squared:  0.2407,	Adjusted R-squared:  0.2344 
F-statistic: 38.31 on 6 and 725 DF,  p-value: < 2.2e-16

```


Rを使ってHRSEを計算する際には、sandwich パッケージが用意されている。
以下が実行例。

```r
library("sandwich")
library("lmtest")

# 上のlm()の結果を引数にする。
coeftest(lm, vcov = vcovHC(lm, type = "HC0"))

```

vcvHC()関数によって、lm()関数の実行結果に対して、分散不均一の場合の分散共分散行列が計算される。その結果を用いて、coeftest()関数によって仮説検定を実行している。
実行結果

```
t test of coefficients:

               Estimate  Std. Error t value  Pr(>|t|)    
(Intercept)  0.87913808  0.30162219  2.9147 0.0036700 ** 
sat          0.00087206  0.00023824  3.6604 0.0002701 ***
hsperc      -0.00559876  0.00165096 -3.3912 0.0007335 ***
tothrs       0.01206297  0.00112596 10.7135 < 2.2e-16 ***
female       0.16667434  0.07850527  2.1231 0.0340839 *  
black       -0.02606973  0.18100090 -0.1440 0.8855160    
white        0.01429927  0.16743329  0.0854 0.9319647    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

```


上の実行結果と比較して、推定量は変わらず、標準誤差の値が変化していることが確認できる。この例では、有意性の評価に大きな違いは見られない。




## おわりに
Stataとかだと、HRSEは簡単に計算できた気がするけど、R、SPSSだと少し面倒くさい。不勉強ゆえかRパッケージの存在も最近まで知らなかった。不均一分散への対処って計量経済学の外だとあんまりやられないのだろうか。。。

他方、計量経済学ではわりと早いうちに習う重要な話題だが、実践上、大きな問題になることは少ないかもしれない。例で示したように、通常の標準誤差とHRSEの値は多くの場合、大標本のもとではほとんど変わらない。
むしろ、大きく変わるようなことがあれば、モデルの誤りを疑ったり、あまりに偏りのあるデータなら対数値をとるなどして、過剰な分散不均一を抑えたりするなどの対処を検討すべきである。

## 関連記事、リンク
英語での同趣旨の記事 こっちのほうが詳しい
https://www.brodrigues.co/blog/2018-07-08-rob_stderr/

回帰分析に関する各種用語の説明
http://knknkn.hatenablog.com/entry/2018/03/25/140418


## 参考文献
Angirist, Joshua D. and Jorn-Steffen Pischke(2009) "Mostly Harmless Econometrics: An Empiricist's Companion" Princeton University Press.  
Wooldridge, Jeffrey M. "Introductory Econometrics: A Modern Approach, 5th edition" South-Western Cengage Learning 

[^1]: ただし、（あたりまえだけど）HRSEは万能ではない。通常の標準誤差と、HRSEの値があまりにもかけ離れる場合は、モデルの誤りを疑ったほうが良い。  
[^2]: 標準誤差に修正をかけるのではなく、モデルから変えてバイアスを取り除こうという発想もある。詳しくは、Feasible Generalized Least Squareで調べよう。
[^3]: 一般的に、と書いたのは、そうでない場合もあるという意味。小標本（有限標本）において、均一分散の仮定が満たされている場合は、HRSEの方に下方バイアスがかかる可能性がある。Angiristたちの教科書では、この場合、通常の標準誤差と、HRSEのうち、大きいほうを使うべきと提案されている。