40
$\begingroup$

I have been trying to replicate the results of the Stata option robust in R. I have used the rlm command form the MASS package and also the command lmrob from the package "robustbase". In both cases the results are quite different from the "robust" option in Stata. Can anybody please suggest something in this context?

Here are the results I obtained when I ran the robust option in 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

And this is what I obtained in R with the lmrob option:

> 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 
$\endgroup$
  • 3
    $\begingroup$ Welcome to Cross Validated! I made your title a little bit more descriptive and added some formatting. In general, programming questions aren't on topic here, but I think yours is because it involves some statistical issues. Hope to see you around.... $\endgroup$ – Matt Krause Sep 28 '14 at 13:14
  • 3
    $\begingroup$ It would help tremendously if you at least pasted the code used to estimate the models in Stata and R (even better if you provide an entirely reproducible example). When you say "results differ" - if you are estimating the same model only the standard errors should differ, not the coefficient estimates. $\endgroup$ – Andy W Sep 28 '14 at 13:23
  • $\begingroup$ okay... These are the results I obtained by the robust option in STATA: $\endgroup$ – user56579 Sep 28 '14 at 13:47
  • 5
    $\begingroup$ looks like lmrob is not the same as reg y x, robust. Google "heteroskedasticity-consistent standard errors R". You'll get pages showing you how to use the lmtest and sandwich libraries. $\endgroup$ – generic_user Sep 28 '14 at 14:12
  • 3
    $\begingroup$ Stata uses a small sample correction factor of n/(n-k). R usually does something else, so make sure you adjust for that. $\endgroup$ – Dimitriy V. Masterov Sep 28 '14 at 15:36
48
$\begingroup$

Charles is nearly there in his answer, but robust option of the regress command (and other regression estimation commands) in Stata makes it possible to use multiple types of heteroskedasticity and autocorrelation robust variance-covariance matrix estimators, as does the coeftest function in the lmtest package, which in turn depends on the respective variance-covariance matrices produced by the vcovHC function in the sandwich package.

However, the default variance-covariance matrices used by the two is different:
1. The default variance-covariance matrix returned by vcocHC is the so-called HC3 for reasons described in the man page for vcovHC.
2. The sandwich option used by Charles makes coeftest use the HC0 robust variance-covariance matrix.
3. To reproduce the Stata default behavior of using the robust option in a call to regress you need to request vcovHC to use the HC1 robust variance-covariance matrix.

Read more about it here.

The following example that demonstrates all the points made above is based on the example here.

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)

The last line of code above reproduces results from Stata:

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

I found a description on the following website that replicates Stata's ''robust'' option in R.

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

Following the instructions, all you need to do is load a function into your R session and then set the parameter ''robust'' in you summary function to TRUE.

summary(lm.object, robust=TRUE)
$\endgroup$
10
$\begingroup$

As of April 2018 I believe you want the estimatr package, which provides a near drop in replacement for lm. Several examples pulled nearly from the documentation:

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

The car package then makes it easy to perform omnibus hypothesis tests for these models:

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
$\endgroup$
4
$\begingroup$

I'd edit the question. You're confusing robust regression with Stata's robust command. There seems to be no benefit to introducing this confusion.

I think there are a few approaches. I haven't looked at them all and not sure which is the best:

The sandwich package:

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

But this doesn't give me the same answers I get from Stata for some reason. I've never tried to work out why - but above in comments there is a suggested answer - I just don't use this package.

The rms package:

I find this a bit of a pain to work with but usually get good answers with some effort. And it is the most useful for me.

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

You can code it from scratch

See this blog post (http://thetarzan.wordpress.com/2011/05/28/heteroskedasticity-robust-and-clustered-standard-errors-in-r/). It looks like the most painful option, but remarkably easy and this option often works the best.

$\endgroup$
  • 4
    $\begingroup$ Charles is correct on the main point, but to make explicit what is implied elsewhere note that Stata has no robust command! (There is a programmer's command _robust, not directly relevant here.) Rather, to get robust (Huber-Eicker-White-sandwich) standard errors, the modern approach in Stata is to specify vce(robust) as an option. The older approach to specify a robust option still works. More broadly, the confusion caused by the difference between robust regression (etc.) and "robust" SEs is unfortunate. $\endgroup$ – Nick Cox Sep 29 '14 at 12:32
  • $\begingroup$ Hey. Thanks a lot. The codes work and it does indeed provide with the results that Stata does. Just a question. I understand that robust regression is different from robust standard errors, and that robust regression is used when your data contains outliers. But it also solves the problem of heteroskedasticity. Could anyone please tell me whether the MM kind estimation provided by the "lmrob" command from package "robustbase" be used as a solution to the problem of outliers and heteroskedasticity simultaneuosly? $\endgroup$ – user56579 Sep 29 '14 at 14:49
  • $\begingroup$ @user56579 My guess is that you want to ask a separate question about this. $\endgroup$ – tchakravarty Oct 4 '14 at 2:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.