题
我想找到最类似于STATA输出的R实现,将最小二乘回归函数与HeteroSkedastic校正的标准误差拟合。具体而言,我希望校正后的标准误差在“摘要”中,而不必为我的初始假设测试进行其他计算。我正在寻找一种与Eviews和Stata提供的“干净”的解决方案。
到目前为止,使用“ lmtest”软件包,我能想到的最好的是:
model <- lm(...)
coeftest(model, vcov = hccm)
这为我提供了想要的输出,但是它似乎并没有用于其陈述的目的。我还必须将摘要带有错误的标准错误来读取R^2和f stat,等等。考虑到动态R的动态性,我觉得应该存在“一行”解决方案。
谢谢
解决方案
我认为你在正确的轨道上 coeftest
在lmtest包中。看一下 三明治包 其中包括此功能,旨在与您已经找到的LMTest软件包并驾齐驱。
> # generate linear regression relationship
> # with Homoskedastic variances
> x <- sin(1:100)
> y <- 1 + x + rnorm(100)
> ## model fit and HC3 covariance
> fm <- lm(y ~ x)
> vcovHC(fm)
(Intercept) x
(Intercept) 0.010809366 0.001209603
x 0.001209603 0.018353076
> coeftest(fm, vcov. = vcovHC)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.01973 0.10397 9.8081 3.159e-16 ***
x 0.93992 0.13547 6.9381 4.313e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
要进行F测试,请查看功能 waldtest()
:
> waldtest(fm, vcov = vcovHC)
Wald test
Model 1: y ~ x
Model 2: y ~ 1
Res.Df Df F Pr(>F)
1 98
2 99 -1 48.137 4.313e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
如果您想要单线...
有很多例子 使用HC和HAC协方差矩阵估计器的计量计算器 Vignette附带了链接LMTest和Sandwich的三明治包装,以完成您想要的工作。
编辑: 单线可能很简单:
mySummary <- function(model, VCOV) {
print(coeftest(model, vcov. = VCOV))
print(waldtest(model, vcov = VCOV))
}
我们可以这样使用(在上面的示例上):
> mySummary(fm, vcovHC)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.01973 0.10397 9.8081 3.159e-16 ***
x 0.93992 0.13547 6.9381 4.313e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Wald test
Model 1: y ~ x
Model 2: y ~ 1
Res.Df Df F Pr(>F)
1 98
2 99 -1 48.137 4.313e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
其他提示
我找到了一个r函数,可以正是您想要的。它为您提供了强大的标准错误,而无需进行其他计算。你跑 summary()
在lm.Object上,如果您设置参数 robust=T
它为您提供了类似Stata的异方差的一致标准错误。
summary(lm.object, robust=T)
您可以在 https://econonictheoryblog.com/2016/08/08/robust-standard-andard-errors-in-r/
现在有一个单行解决方案 lm_robust
来自 estimatr
包裹, ,您可以从Cran中安装 install.packages(estimatr)
.
> library(estimatr)
> lmro <- lm_robust(mpg ~ hp, data = mtcars, se_type = "stata")
> summary(lmro)
Call:
lm_robust(formula = mpg ~ hp, data = mtcars, se_type = "stata")
Standard error type: HC1
Coefficients:
Estimate Std. Error Pr(>|t|) CI Lower CI Upper DF
(Intercept) 30.09886 2.07661 4.348e-15 25.85785 34.33987 30
hp -0.06823 0.01356 2.132e-05 -0.09592 -0.04053 30
Multiple R-squared: 0.6024 , Adjusted R-squared: 0.5892
F-statistic: 45.46 on 1 and 30 DF, p-value: 1.788e-07
您还可以获取整理的输出:
> tidy(lmro)
term estimate std.error p.value ci.lower
1 (Intercept) 30.09886054 2.0766149 4.347723e-15 25.85784704
2 hp -0.06822828 0.0135604 2.131785e-05 -0.09592231
ci.upper df outcome
1 34.33987404 30 mpg
2 -0.04053425 30 mpg
这 "stata"
标准错误默认为“ HC1”标准错误,这是默认错误 rob
Stata中的标准错误。你也可以得到 "classical", "HC0", "HC1", "HC2", "HC3"
以及各种聚类的标准错误(包括与Stata匹配的标准错误)。