# How to replicate Stata's robust binomial GLM for proportion data in R?

There is an example on how to run a GLM for proportion data in Stata here

The IV is the proportion of students receiving free or reduced priced meals at school. The stata model looks like this.:

glm meals yr_rnd parented api99, link(logit) family(binomial) robust nolog


I'm interested in learning how to replicate this results in R (ideally using the same robust approach). Lets imagine that I have data about the number of students receiving free meals (Successes) and the rest of the students (Failures). I'm guessing the model in R could look something like this:

fitglm <- glm(cbind(Successes,Failures) ~ yr_rnd + parented + api99, family=binomial)


Also, it was pointed out to me elsewhere (Penguin_Knight) that the error message "meals has non-integer values" could be bad. I'm clueless regarding this error...

• In Stata vce(robust) rather than robust is what is now documented, but robust should still work. Commented Mar 14, 2014 at 11:43

Using the R package sandwich and lmtest, you can replicate the results like that (I assume that you've already downloaded the dataset or access it over the internet):

#-----------------------------------------------------------------------------
#-----------------------------------------------------------------------------

require(foreign)
require(sandwich)

#-----------------------------------------------------------------------------
#-----------------------------------------------------------------------------

#-----------------------------------------------------------------------------
# Inspect dataset
#-----------------------------------------------------------------------------

str(dat)

#-----------------------------------------------------------------------------
# Fit the glm
#-----------------------------------------------------------------------------

fitglm <- glm(meals ~ yr_rnd + parented + api99, family = binomial(logit), data = dat)

#-----------------------------------------------------------------------------
# Output of the model
#-----------------------------------------------------------------------------

summary(fitglm)

#-----------------------------------------------------------------------------
# Calculate robust standard errors by hand
#-----------------------------------------------------------------------------

cov.m1 <- vcovHC(fitglm, type = "HC1")

std.err <- sqrt(diag(cov.m1))

q.val <- qnorm(0.975)

r.est <- cbind(
Estimate = coef(fitglm)
, "Robust SE" = std.err
, z = (coef(fitglm)/std.err)
, "Pr(>|z|) "= 2 * pnorm(abs(coef(fitglm)/std.err), lower.tail = FALSE)
, LL = coef(fitglm) - q.val  * std.err
, UL = coef(fitglm) + q.val  * std.err
)

r.est


The model output using robust standard errors is:

                Estimate    Robust SE          z     Pr(>|z|)            LL           UL
(Intercept)  6.801682703 0.0724029936  93.942009  0.000000e+00  6.659775443  6.943589963
yr_rndYes    0.048252657 0.0321827112   1.499335  1.337868e-01 -0.014824298  0.111329612
parented    -0.766259824 0.0390852844 -19.604816  1.406590e-85 -0.842865573 -0.689654074
api99       -0.007304603 0.0002156354 -33.874790 1.566480e-251 -0.007727241 -0.006881966


A much more convenient way is using the coeftest and coefci functions from the lmtest package (output not shown but is identical to output above):

coeftest(fitglm, vcov. = vcovHC(fitglm, type = "HC1"))
coefci(fitglm, vcov. = vcovHC(fitglm, type = "HC1"))


The estimates and standard errors are fairly similar to those calculated using Stata but not exactly. The reason is that Stata uses a finite-sample adjustment (see this post). The Stata-output is (caution: I enter the variable yr_rnd as categorical variable to replicate R's behaviour, unlike the UCLA page):

------------------------------------------------------------------------------
|               Robust
meals |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
yr_rnd |
Yes  |   .0482527   .0321714     1.50   0.134    -.0148021    .1113074
parented |  -.7662598   .0390715   -19.61   0.000    -.8428386   -.6896811
api99 |  -.0073046   .0002156   -33.89   0.000    -.0077271   -.0068821
_cons |   6.801683   .0723775    93.98   0.000     6.659825     6.94354
------------------------------------------------------------------------------


To exactly replicate Stata's standard errors, we have to use @AchimZeileis' function (found here):

sandwich1 <- function(object, ...) sandwich(object) * nobs(object) / (nobs(object) - 1)
coeftest(fitglm, vcov. = sandwich1)

Estimate  Std. Error  z value Pr(>|z|)
(Intercept)  6.80168270  0.07237747  93.9751   <2e-16 ***
yr_rndYes    0.04825266  0.03217137   1.4999   0.1336
parented    -0.76625982  0.03907151 -19.6117   <2e-16 ***
api99       -0.00730460  0.00021556 -33.8867   <2e-16 ***


There are several methods available for the function vcovHC. Consult the help file of vcovHC for the details.

Note that if you use the option family = quasibinomial(logit), there will be no error message (see here).

• +1. What would happen if you use glm() with family=quasibinomial? Isn't it supposed to estimate robust standard errors by itself, or at least do something conceptually similar by computing standard errors accounting for over-dispersion? Commented Sep 5, 2016 at 19:35
• @amoeba: while family = quasibinomial leads to a robust estimate of the variance ($\alpha \hat p (1- \hat p)$), it is not the same robust estimator as the sandwich estimator (in which for each observation, the variance is estimated as $(y_i - \hat y_i)^2$ rather than by some assumed functional form). Commented Sep 5, 2016 at 22:39
• @CliffAB Thanks a lot! This agrees with what I've been reading during the last hour, e.g. in this book or in this SO post and linked SO posts. (Do you have better references?) It is a pity we do not seem to have a good CV thread that would accurately explain different approaches. Commented Sep 5, 2016 at 22:43
• @amoeba: I'm afraid I don't have any good references off hand; I just learned about it in a class which did not use a book. I recall my professor mentioning that it was an idea that actually began in the field of economics before being readily accepted by statisticians. But that's all I've got. Commented Sep 5, 2016 at 23:34

You can replicate the UCLA FAQ on proportions (with a percentage as a dependent variable) as follows:

require(foreign);require(lmtest);require(sandwich)
fitperc <- glm(meals ~ yr_rnd + parented + api99, family = binomial, data=meals)
## Warning message:
## In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!


I don't know if the warning above is an issue here or not. For some reason the intercept don't match in R and Stata, but since we don't interpret it usually in logit/probit anyway it shouldn't matter much.

summary(fitperc)
##
## Call:
## glm(formula = meals ~ yr_rnd + parented + api99, family = binomial,
##     data = meals, na.action = na.exclude)
##
## Deviance Residuals:
##      Min        1Q    Median        3Q       Max
## -1.77722  -0.18995  -0.01649   0.18692   1.60959
##
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)  6.801683   0.231914  29.329   <2e-16 ***
## yr_rndYes    0.048253   0.104210   0.463    0.643
## parented    -0.766260   0.090733  -8.445   <2e-16 ***
## api99       -0.007305   0.000506 -14.435   <2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
##     Null deviance: 1953.94  on 4256  degrees of freedom
## Residual deviance:  395.81  on 4253  degrees of freedom
##   (164 observations deleted due to missingness)
## AIC: 2936.7
##
## Number of Fisher Scoring iterations: 5


In R the small-sample corrections used are different than those in Stata, but the robust SEs are fairly similar:

coeftest(fitperc, function(x) vcovHC(x, type = "HC1"))
##
## z test of coefficients:
##
##                Estimate  Std. Error  z value Pr(>|z|)
## (Intercept)  6.80168270  0.07240299  93.9420   <2e-16 ***
## yr_rndYes    0.04825266  0.03218271   1.4993   0.1338
## parented    -0.76625982  0.03908528 -19.6048   <2e-16 ***
## api99       -0.00730460  0.00021564 -33.8748   <2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


To use the exact same small-sample correction you need to follow this post:

sandwich1 <- function(object, ...) sandwich(object) * nobs(object) / (nobs(object) - 1)
coeftest(fitperc, vcov = sandwich1)
##
## z test of coefficients:
##
##                Estimate  Std. Error  z value Pr(>|z|)
## (Intercept)  6.80168270  0.07237747  93.9751   <2e-16 ***
## yr_rndYes    0.04825266  0.03217137   1.4999   0.1336
## parented    -0.76625982  0.03907151 -19.6117   <2e-16 ***
## api99       -0.00730460  0.00021556 -33.8867   <2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


The log likelihood and the confidence intervals (slightly different as the estimation procedure seems to be different):

logLik(fitperc)
## 'log Lik.' -1464.363 (df=4)
confint(fitperc)
## Waiting for profiling to be done...
##                    2.5 %       97.5 %
## (Intercept)  6.352788748  7.262067304
## yr_rndYes   -0.155529338  0.253123151
## parented    -0.944775733 -0.588903012
## api99       -0.008303668 -0.006319185


To obtain the predictions:

meals_pred <- data.frame(api99=rep(c(500,600,700), 2),
yr_rnd=rep(c("No", "Yes"), times=1, each=3),
parented=rep(2.5, 6))
cbind(meals_pred, pred=predict(fitperc, meals_pred, "response"))
##   api99 yr_rnd parented      pred
## 1   500     No      2.5 0.7744710
## 2   600     No      2.5 0.6232278
## 3   700     No      2.5 0.4434458
## 4   500    Yes      2.5 0.7827873
## 5   600    Yes      2.5 0.6344891
## 6   700    Yes      2.5 0.4553849


See this question for a related discussion: