問題

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

例（２つの回帰モデルの場合）

モデルが２つ以上のときは，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