查尔斯几乎在他的答案中，但是Statarobust中的命令选项regress（和其他回归估计命令）使得可以使用多种类型的异方差和自相关鲁棒方差 - 协方差矩阵估计器，包中的coeftest函数也是如此lmtest，其中turn 取决于包中vcovHC函数产生的各自的方差-协方差矩阵sandwich。

但是，两者使用的默认方差-协方差矩阵是不同的：
1. 返回的默认方差-协方差矩阵vcocHC就是所谓HC3的man page for 中描述的原因vcovHC。
2. sandwichCharles 使用的选项coeftest使用了HC0稳健的方差-协方差矩阵。3. 要重现 Stata在调用
中使用选项的默认行为，您需要请求使用稳健的方差-协方差矩阵。robustregressvcovHCHC1

在此处阅读更多相关信息。

以下示例演示了上述所有要点，基于此处的示例。

library(foreign)
library(sandwich)
library(lmtest)

dfAPI = read.dta("http://www.ats.ucla.edu/stat/stata/webbooks/reg/elemapi2.dta")
lmAPI = lm(api00 ~ acs_k3 + acs_46 + full + enroll, data= dfAPI)
summary(lmAPI)                                  # non-robust

# check that "sandwich" returns HC0
coeftest(lmAPI, vcov = sandwich)                # robust; sandwich
coeftest(lmAPI, vcov = vcovHC(lmAPI, "HC0"))    # robust; HC0 

# check that the default robust var-cov matrix is HC3
coeftest(lmAPI, vcov = vcovHC(lmAPI))           # robust; HC3 
coeftest(lmAPI, vcov = vcovHC(lmAPI, "HC3"))    # robust; HC3 (default)

# reproduce the Stata default
coeftest(lmAPI, vcov = vcovHC(lmAPI, "HC1"))    # robust; HC1 (Stata default)

上面的最后一行代码复制了来自 Stata 的结果：

use http://www.ats.ucla.edu/stat/stata/webbooks/reg/elemapi2
regress api00 acs_k3 acs_46 full enroll, robust

我在以下网站上找到了一个描述，它在 R 中复制了 Stata 的“强大”选项。

https://economictheoryblog.com/2016/08/08/robust-standard-errors-in-r

按照说明操作，您只需将函数加载到 R 会话中，然后将摘要函数中的参数“robust”设置为 TRUE。

summary(lm.object, robust=TRUE)

截至 2018 年 4 月，我相信你想要这个estimatr包，它几乎可以替代lm. 几乎从文档中提取了几个示例：

library(estimatr)
library(car)

# HC1 robust standard errors
model <- lm_robust(GPA_year2 ~ gpa0 + ssp, data = alo_star_men,
                   se_type = "stata")
summary(model)
#> 
#> Call:
#> lm_robust(formula = GPA_year2 ~ gpa0 + ssp, data = alo_star_men, 
#>     se_type = "stata")
#> 
#> Standard error type:  HC1 
#> 
#> Coefficients:
#>             Estimate Std. Error  Pr(>|t|) CI Lower CI Upper  DF
#> (Intercept) -3.60625    1.60084 0.0258665 -6.77180  -0.4407 137
#> gpa0         0.06814    0.02024 0.0009868  0.02812   0.1082 137
#> ssp          0.31917    0.18202 0.0817589 -0.04077   0.6791 137
#> 
#> Multiple R-squared:  0.09262 ,   Adjusted R-squared:  0.07937 
#> F-statistic: 6.992 on 2 and 137 DF,  p-value: 0.001284

# HC1 cluster robust standard errors
model2 <- lm_robust(GPA_year2 ~ gpa0 + ssp, cluster = ssp,
                   data = alo_star_men, se_type = "stata")
summary(model2)
#> 
#> Call:
#> lm_robust(formula = GPA_year2 ~ gpa0 + ssp, data = alo_star_men, 
#>     clusters = ssp, se_type = "stata")
#> 
#> Standard error type:  stata 
#> 
#> Coefficients:
#>             Estimate Std. Error Pr(>|t|) CI Lower CI Upper DF
#> (Intercept) -3.60625   1.433195 0.240821 -21.8167  14.6042  1
#> gpa0         0.06814   0.018122 0.165482  -0.1621   0.2984  1
#> ssp          0.31917   0.004768 0.009509   0.2586   0.3798  1
#> 
#> Multiple R-squared:  0.09262 ,   Adjusted R-squared:  0.07937 
#> F-statistic: 6.992 on 2 and 137 DF,  p-value: 0.001284

然后，该car软件包可以轻松地为这些模型执行综合假设检验：

linearHypothesis(model, c("gpa0 = ssp"))
#> Linear hypothesis test
#> 
#> Hypothesis:
#> gpa0 - ssp = 0
#> 
#> Model 1: restricted model
#> Model 2: GPA_year2 ~ gpa0 + ssp
#> 
#>   Res.Df Df  Chisq Pr(>Chisq)
#> 1    138                     
#> 2    137  1 1.8859     0.1697

我会编辑这个问题。您将稳健回归与 Stata 的稳健命令混淆了。引入这种混淆似乎没有任何好处。

我认为有几种方法。我没有看过它们，不确定哪个是最好的：

三明治包装：

library(sandwich)    
coeftest(model, vcov=sandwich)

但是由于某种原因，这并没有给我从 Stata 得到的相同答案。我从来没有试图找出原因 - 但在评论上面有一个建议的答案 - 我只是不使用这个包。

rms 包：

我发现这有点痛苦，但通常会通过一些努力得到好的答案。它对我来说是最有用的。

model = ols(a~b, x=TRUE)    
robcov(model)

您可以从头开始编写代码

请参阅此博客文章（http://thetarzan.wordpress.com/2011/05/28/heteroskedasticity-robust-and-clustered-standard-errors-in-r/）。它看起来是最痛苦的选择，但非常容易，而且这种选择通常效果最好。