There are already a lot of good questions on this topic (e.g., here). But they address complexities that I am not interested in.
I have some simple data. I am using basic GLM and OLS, with robust variance estimators. In Stata, I'm inputting:
input deadspace asthma
43 0
44 0
45 0
56 0
56 0
57 0
58 0
64 0
31 1
78 1
79 1
88 1
92 1
101 1
112 1
end
And in R, I'm inputting:
a1 <- c(43,44,45,56,56,57,58,64)
a2 <- c(31,78,79,88,92,101,112)
deadspace <- c(a1,a2)
asthma <- c(rep(0,length(a1)),rep(1,length(a2)))
My objective is to understand how to equate results from Stata and R for simple GLMs. I have succeeded at this task when OLS is used (i.e., regress in Stata, and lm in R). But when I use GLMs, I get different SE estimates, even when I use the same robust variance estimator. For example, in Stata, running:
glm deadspace asthma, irls family(gaussian) link(identity) vce(robust)
Gives a robust standard error for the asthma coefficient of 9.74. Based on other findings, it seems that Stata's default robust variance estimator is HC1. So, using this in R with a corresponding GLM:
mod1 <- glm(deadspace ~ asthma,family=gaussian(link="identity"))
coeftest(mod1, vcov = vcovHC(mod1, type="HC1"))
I get a standard error estimate of 10.11. Trying different HC variations in R does not give me a result that matches the Stata estimate.
However, if I used the "unbiased" option in Stata:
glm deadspace asthma, irls family(gaussian) link(identity) vce(unbiased)
I get an SE estimate of 10.16, which matches the corresponding R estimate when the HC2 variance estimator is used:
mod1 <- glm(deadspace ~ asthma,family=gaussian(link="identity"))
coeftest(mod1, vcov = vcovHC(mod1, type="HC2"))[2,2]
But this is the only way I've been able to equate robust variance estimates with GLM from Stata and R. Why is that?
Based on an article linked here, it seems Stata implements a small sample adjustment by default. But it is not clear to me which adjustment is being used. How can I modify either:
- the Stata code to give the same HC1 estimate of 10.11 that R returns when HC1 is used, or
- the R code to get the SE estimate of 9.74 that Stata returns when
vce(robust)
is used?
Alternatively, I'd settle for the EXACT equation that is being used to compute the value of 9.74 in Stata with vce(robust)
.
coeftest(mod1, vcov = vcovCL)
will give you the same adjustment by default as in Stata, i.e., with $n - 1$ rather than $n$ or $n - k$. See also my comment to the accepted answer. $\endgroup$ – Achim Zeileis Apr 17 '20 at 1:41