Регрессия с исправленными гетероскедастичностью стандартными ошибками - PullRequest
30 голосов
/ 08 декабря 2010

Я хотел бы найти реализацию R, которая больше всего похожа на вывод Stata для подгонки функции регрессии наименьших квадратов с исправленными гетероскедастическими стандартными ошибками.В частности, я хотел бы, чтобы исправленные стандартные ошибки были в «сводке» и не требовали дополнительных вычислений для моего начального этапа проверки гипотез.Я ищу решение, столь же «чистое», как то, что предоставляют Eviews и Stata.

Пока что, используя пакет «lmtest», лучшее, что я могу предложить:1005 * Это дает мне вывод, который я хочу, но он не использует "coeftest" для своей заявленной цели.Мне также пришлось бы использовать сводку с неверными стандартными ошибками для считывания R ^ 2 и F stat и т. Д. Я чувствую, что должно существовать решение этой проблемы в «одну строку», учитывая динамический R.

Спасибо

Ответы [ 3 ]

39 голосов
/ 08 декабря 2010

Я думаю, что вы на правильном пути с 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 виньетка, которая поставляется с пакетом сэндвичей, связывающих lmtest и бутерброд, для того, что вы хотите.

Редактировать: Однострочник может быть простым:

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
10 голосов
/ 23 августа 2016

Я нашел функцию R, которая делает именно то, что вы ищете. Это дает вам надежные стандартные ошибки без необходимости делать дополнительные вычисления. Вы запускаете summary() для объекта lm.object, и если вы устанавливаете параметр robust=T, он возвращает вам стандартные ошибки, подобные гетероскедастичности Stata.

summary(lm.object, robust=T)

Вы можете найти функцию на https://economictheoryblog.com/2016/08/08/robust-standard-errors-in-r/

2 голосов
/ 18 апреля 2018

В настоящее время существует однострочное решение с использованием 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).

...