盆暗の学習記録

データサイエンスを中心として,日々学んだことの備忘録としていく予定です。初心者であり独学なので内容には誤りが含まれる可能性が大いにあります。

[R]{stargazer}でロバスト標準誤差を使用した分析結果を載せる

問題

Rでロバスト標準誤差を使うときはsandwich::vcovHC()ロバスト標準誤差の分散共分散行列を算出し,lmtest::coeftest()を使ってRの標準的な結果表示のスタイルに似せた結果表示をする,というのが基本的・標準的なロバスト標準誤差の呼び出し方かと思っているのですが,

lm1 <- lm(Volume ~ Girth, data = trees)
coef1 <- lmtest::coeftest(lm1, .vcov = sandwich::vcovHC(lm1, type = "HC1"))
> coef1

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) -36.94346    3.36514 -10.978 7.621e-12 ***
Girth         5.06586    0.24738  20.478 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

この場合,stargazer()で綺麗に結果表示できません

> library(stargazer)
> stargazer(coef1, type = "text")

====================================
             Dependent variable:    
         ---------------------------
                                    
------------------------------------
Girth             5.066***          
                   (0.247)          
                                    
Constant         -36.943***         
                   (3.365)          
                                    
====================================
====================================
Note:    *p<0.1; **p<0.05; ***p<0.01

解決策

vcovHC()で分散共分散行列を作って,そこからロバスト標準誤差のベクトルを取得するところまでを行い,それをstargazer()seの引数に指定すればよいみたいです(stargazer()のモデルを指定する部分には通常のlm()の結果を入れる)

vcov <- sandwich::vcovHC(lm1, type = "HC1")
robust_se1 <- sqrt(diag(vcov))
stargazer(lm1, se = robust_se1, type = "text")
> stargazer(lm1, se = robust_se1, type = "text")

===============================================
                        Dependent variable:    
                    ---------------------------
                              Volume           
-----------------------------------------------
Girth                          5.066           
                                               
                                               
Constant                    -36.943***         
                              (3.951)          
                                               
-----------------------------------------------
Observations                    31             
R2                             0.935           
Adjusted R2                    0.933           
Residual Std. Error       4.252 (df = 29)      
F Statistic           419.360*** (df = 1; 29)  
===============================================
Note:               *p<0.1; **p<0.05; ***p<0.01

例(2つの回帰モデルの場合)

モデルが2つ以上のときは,lmの結果とロバスト標準誤差それぞれをリストに入れてstargazer()に入れる

library(stargazer)
# estimate
lm1 <- lm(Volume ~ Girth, data = trees)
lm2 <- lm(Volume ~ Height, data = trees)

# robust SE
vcov <- sandwich::vcovHC(lm1, type = "HC1")
robust_se1 <- sqrt(diag(vcov))
vcov <- sandwich::vcovHC(lm2, type = "HC1")
robust_se2 <- sqrt(diag(vcov))

# print
stargazer(list(lm1, lm2), se=list(robust_se1, robust_se2), type = "text")
> stargazer(list(lm1, lm2), se=list(robust_se1, robust_se2), type = "text")

==========================================================
                                  Dependent variable:     
                              ----------------------------
                                         Volume           
                                   (1)            (2)     
----------------------------------------------------------
Girth                            5.066***                 
                                 (0.317)                  
                                                          
Height                                         1.543***   
                                                (0.427)   
                                                          
Constant                        -36.943***    -87.124***  
                                 (3.951)       (31.008)   
                                                          
----------------------------------------------------------
Observations                        31            31      
R2                                0.935          0.358    
Adjusted R2                       0.933          0.336    
Residual Std. Error (df = 29)     4.252         13.397    
F Statistic (df = 1; 29)        419.360***     16.164***  
==========================================================
Note:                          *p<0.1; **p<0.05; ***p<0.01

参考

なお,stargazerのような表記じゃなくデータフレームの形でロバスト標準誤差を含む分析結果がほしい場合は{broom}パッケージが便利です。