Obtaining robust standard errors and odds ratios for logistic regression in R

I’ve always found it frustrating how it’s so easy to produce robust standard errors in Stata and in R it’s so complicated. First, we have to estimate the standard errors separately and then replace the previous standard errors with the new ones. Second, if you want to estimate odds ratios instead of logit coefficients, then the robust standard errors need to be scaled. All of that is as simple as adding robust or in the Stata logit command.

I decided to make it as simple in R.

First, I have to give credit to Achim Zeileis in this question because he provided part of code to generate the robust standard errors.

The function accepts a glm object and can return logit coefficients with robust standard errors, odd ratios with adjusted robust standard errors or probability scaled coefficients with adjusted robust standard errors.

Here’s the function:

# This function estimates robust standad error for glm objects and
# returns coefficients as either logit, odd ratios or probabilities.
# logits are default
# argument x must be glm model.


# Credit to Achim here:
# http://stackoverflow.com/questions/27367974/
# different-robust-standard-errors-of-logit-regression-in-stata-and-r
# for the code in line 14 and 15

robustse <- function(x, coef = c("logit", "odd.ratio", "probs")) {
suppressMessages(suppressWarnings(library(lmtest)))
suppressMessages(suppressWarnings(library(sandwich)))

    sandwich1 <- function(object, ...) sandwich(object) *
                                       nobs(object) / (nobs(object) - 1)
    # Function calculates SE's
    mod1 <- coeftest(x, vcov = sandwich1) 
    # apply the function over the variance-covariance matrix
    
    if (coef == "logit") {
    return(mod1) # return logit with robust SE's
    } else if (coef == "odd.ratio") {
    mod1[, 1] <- exp(mod1[, 1]) # return odd ratios with robust SE's
    mod1[, 2] <- mod1[, 1] * mod1[, 2]
    return(mod1)
    } else {
    mod1[, 1] <- (mod1[, 1]/4) # return probabilites with robust SE's
    mod1[, 2] <- mod1[, 2]/4
    return(mod1)
    }
}

Now let’s give it a try. Let’s estimate two models, one with logit coefficients and robust SE’s and the same for odds ratios. Just to make sure, let’s compare it with the Stata output.

# In R for logit coefficients and robust standard errors:
suppressMessages(suppressWarnings(library(haven)))

dat <- read_dta("http://www.stata-press.com/data/r9/quad1.dta")
mod1 <- glm(z ~ x1 + x2 + x3, dat, family = binomial)
robustse(mod1, coef = "logit")

## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## (Intercept) 0.091740   0.025869  3.5464 0.0003905 ***
## x1          0.024050   0.025869  0.9297 0.3525395    
## x2          0.157143   0.089715  1.7516 0.0798448 .  
## x3          0.190162   0.089418  2.1267 0.0334490 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# In Stata for logit coefficients and robust standard errors:
use http://www.stata-press.com/data/r9/quad1.dta, clear
logit z x1 x2 x3, robust

# ------------------------------------------------------------------------------
#              |               Robust
#            z |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
# -------------+----------------------------------------------------------------
#           x1 |     .02405   .0258693     0.93   0.353    -.0266528    .0747529
#           x2 |   .1571428   .0897145     1.75   0.080    -.0186944      .33298
#           x3 |    .190162   .0894185     2.13   0.033      .014905    .3654189
#        _cons |     .09174   .0258686     3.55   0.000     .0410386    .1424415
# ------------------------------------------------------------------------------

# In R for odd ratios with adjusted standard errors
robustse(mod1, coef = "odd.ratio")

## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## (Intercept) 1.096080   0.028354  3.5464 0.0003905 ***
## x1          1.024342   0.026499  0.9297 0.3525395    
## x2          1.170163   0.104981  1.7516 0.0798448 .  
## x3          1.209445   0.108147  2.1267 0.0334490 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# In Stata for odd ratios with adjusted standard errors
logit z x1 x2 x3, robust or

# ------------------------------------------------------------------------------
#              |               Robust
#            z | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
# -------------+----------------------------------------------------------------
#           x1 |   1.024342    .026499     0.93   0.353     .9736992    1.077618
#           x2 |   1.170163   .1049806     1.75   0.080     .9814792    1.395119
#           x3 |   1.209445   .1081467     2.13   0.033     1.015017    1.441118
#        _cons |    1.09608    .028354     3.55   0.000     1.041892    1.153086
# ------------------------------------------------------------------------------

You can also use the stargazer package to produce nicely formatted tables with these new estimates and it should be exactly the same.

---

UPDATE: I’ve included this function my personal package which you can install with devtools::install_github("cimentadaj/cimentadaj"). Feel free to make any pull requests in the github repo.