Replikacja „solidnej” opcji Staty w R.

39

Próbowałem zreplikować wyniki opcji Stata robustw R. Użyłem rlmpolecenia z pakietu MASS, a także polecenia lmrobz pakietu „robustbase”. W obu przypadkach wyniki są zupełnie inne niż „solidna” opcja w Stacie. Czy ktoś może zasugerować coś w tym kontekście?

Oto wyniki, które uzyskałem, gdy uruchomiłem solidną opcję w Stata:

. reg yb7 buildsqb7 no_bed no_bath rain_harv swim_pl pr_terrace, robust

Linear regression                                      Number of obs =    4451
                                                       F(  6,  4444) =  101.12
                                                       Prob > F      =  0.0000
                                                       R-squared     =  0.3682
                                                       Root MSE      =   .5721

------------------------------------------------------------------------------
             |               Robust
         yb7 |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   buildsqb7 |   .0046285   .0026486     1.75   0.081    -.0005639     .009821
      no_bed |   .3633841   .0684804     5.31   0.000     .2291284    .4976398
     no_bath |   .0832654   .0706737     1.18   0.239    -.0552904    .2218211
   rain_harv |   .3337906   .0395113     8.45   0.000     .2563289    .4112524
     swim_pl |   .1627587   .0601765     2.70   0.007     .0447829    .2807346
  pr_terrace |   .0032754   .0178881     0.18   0.855    -.0317941    .0383449
       _cons |   13.68136   .0827174   165.40   0.000     13.51919    13.84353

I to właśnie uzyskałem w R z opcją lmrob:

> modelb7<-lmrob(yb7~Buildsqb7+No_Bed+Rain_Harv+Swim_Pl+Gym+Pr_Terrace, data<-bang7)
> summary(modelb7)

Call:
lmrob(formula = yb7 ~ Buildsqb7 + No_Bed + Rain_Harv + Swim_Pl + Gym + Pr_Terrace, 
    data = data <- bang7)
 \--> method = "MM"
Residuals:
      Min        1Q    Median        3Q       Max 
-51.03802  -0.12240   0.02088   0.18199   8.96699 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.648261   0.055078 229.641   <2e-16 ***
Buildsqb7    0.060857   0.002050  29.693   <2e-16 ***
No_Bed       0.005629   0.019797   0.284   0.7762    
Rain_Harv    0.230816   0.018290  12.620   <2e-16 ***
Swim_Pl      0.065199   0.028121   2.319   0.0205 *  
Gym          0.023024   0.014655   1.571   0.1162    
Pr_Terrace   0.015045   0.013951   1.078   0.2809    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Robust residual standard error: 0.1678 
Multiple R-squared:  0.8062,    Adjusted R-squared:  0.8059 
użytkownik56579
źródło
3
Witamy w Cross Validated! Uczyniłem twój tytuł nieco bardziej opisowym i dodałem trochę formatowania. Ogólnie rzecz biorąc, pytania dotyczące programowania nie są tutaj na temat, ale myślę, że twoje, ponieważ dotyczą pewnych problemów statystycznych. Mam nadzieję, że się zobaczymy ...
Matt Krause,
3
Byłoby to bardzo pomocne, gdybyś przynajmniej wkleił kod użyty do oszacowania modeli w Stata i R (nawet lepiej, jeśli podasz całkowicie powtarzalny przykład). Kiedy mówisz „wyniki różnią się” - jeśli szacujesz ten sam model, tylko standardowe błędy powinny się różnić, a nie szacunki współczynników.
Andy W
okej ... Oto wyniki, które uzyskałem dzięki solidnej opcji w STATA:
user56579
5
wygląda na lmrobto, że to nie to samo co reg y x, robust. Google „błędy standardowe zgodne z heteroskedastycznością R”. Otrzymasz strony pokazujące, jak korzystać z bibliotek lmtesti sandwich.
generic_user
3
Stata używa małego współczynnika korekcji próbki wynoszącego n / (nk). R zwykle robi coś innego, więc dostosuj się do tego.
Dimitriy V. Masterov

Odpowiedzi:

47

Charles jest już prawie w swojej odpowiedzi, ale robustopcja regresspolecenia (i innych poleceń szacowania regresji) w Stacie umożliwia użycie wielu typów heteroskedastyczności i autokorelacji solidnych estymatorów macierzy wariancji-kowariancji, podobnie jak coeftestfunkcja w lmtestpakiecie, która w turn zależy od odpowiednich macierzy wariancji-kowariancji wytwarzanych przez vcovHCfunkcję w sandwichpakiecie.

Jednak domyślna
macierz wariancji-kowariancji używana przez te dwa elementy jest inna: 1. Domyślna zwrócona macierz wariancji-kowariancji vcocHCto tak zwane HC3przyczyny opisane na stronie podręcznika vcovHC.
2. sandwichOpcja zastosowana przez Charlesa coeftestwykorzystuje HC0solidną macierz wariancji-kowariancji.
3. Aby odtworzyć domyślne zachowanie Stata dotyczące korzystania z robustopcji w wywołaniu regress, należy poprosić vcovHCo użycie HC1solidnej macierzy wariancji-kowariancji.

Przeczytaj więcej na ten temat tutaj .

Poniższy przykład, który pokazuje wszystkie powyższe punkty, oparty jest na przykładzie tutaj .

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)

Ostatni wiersz powyższego kodu odtwarza wyniki ze Staty:

use http://www.ats.ucla.edu/stat/stata/webbooks/reg/elemapi2
regress api00 acs_k3 acs_46 full enroll, robust
tchakravarty
źródło
Link do danych nie działa. Czy możesz zaktualizować link? Czy ten sam plik: faculty.smu.edu/tfomby/eco5350/data/Examples/elemapi2.dta ?
vasili111
Jak odtworzyć również przedziały ufności?
vasili111
29

Znalazłem opis na następującej stronie internetowej, który replikuje opcję „solidnej” Staty w R.

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

Postępując zgodnie z instrukcjami, wystarczy załadować funkcję do sesji R, a następnie ustawić parametr „solidny” w funkcji podsumowania na PRAWDA.

summary(lm.object, robust=TRUE)
Alex Rato
źródło
10

Według stanu na kwiecień 2018 roku uważam, że chcesz estimatrpakiet , który zapewnia prawie spadek w zastępstwie lm. Kilka przykładów wyciągnięto prawie z dokumentacji:

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

carPakiet następnie ułatwia wykonywanie zbiorczych testy hipotez tych modeli:

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
alexpghayes
źródło
4

Zredagowałbym pytanie. Mylisz solidną regresję z solidnym poleceniem Staty. Wydaje się, że wprowadzenie tego zamieszania nie przynosi korzyści.

Myślę, że istnieje kilka podejść. Nie obejrzałem ich wszystkich i nie jestem pewien, co jest najlepsze:

Pakiet kanapkowy:

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

Ale z jakiegoś powodu nie daje mi to takich samych odpowiedzi, jakie otrzymuję od Staty. Nigdy nie próbowałem dowiedzieć się, dlaczego - ale powyżej w komentarzach jest sugerowana odpowiedź - po prostu nie używam tego pakietu.

Pakiet rms:

Uważam, że to trochę kłopotliwe w pracy, ale zwykle uzyskuję dobre odpowiedzi z pewnym wysiłkiem. I to jest dla mnie najbardziej przydatne.

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

Możesz go kodować od zera

Zobacz ten post na blogu ( http://thetarzan.wordpress.com/2011/05/28/heteroskedasticity-robust-and-clustered-standard-errors-in-r/ ). Wygląda na najbardziej bolesną opcję, ale niezwykle łatwą i ta opcja często działa najlepiej.

Charles
źródło
4
Charles ma rację w głównej kwestii, ale aby wyjaśnić, co sugeruje gdzie indziej, zauważ, że Stata nie ma robustpolecenia! (Istnieje polecenie programisty _robust, które nie ma bezpośredniego związku z tym.) Zamiast tego, aby uzyskać solidne błędy standardowe (Huber-Eicker-White-sandwich), nowoczesne podejście w Stata jest określane vce(robust)jako opcja. Starsze podejście do określenia robustopcji nadal działa. Mówiąc szerzej, zamieszanie spowodowane różnicą między solidną regresją (itp.) A „solidną” SE jest niefortunne.
Nick Cox
Hej. Wielkie dzięki. Kody działają i rzeczywiście zapewniają wyniki, które robi Stata. Tylko pytanie. Rozumiem, że solidna regresja różni się od solidnych błędów standardowych i że silna regresja jest stosowana, gdy dane zawierają wartości odstające. Ale rozwiązuje również problem heteroskedastyczności. Czy ktoś mógłby mi powiedzieć, czy oszacowanie rodzaju MM dostarczone przez polecenie „lmrob” z pakietu „robustbase” może być użyte jako rozwiązanie problemu wartości odstających i heteroskedastyczności jednocześnie?
user56579,
@ user56579 Domyślam się, że chcesz zadać osobne pytanie na ten temat.
tchakravarty