Robust Variance in Stata and R?

Question

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:

1. the Stata code to give the same HC1 estimate of 10.11 that R returns when HC1 is used, or
2. 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).

Answer 1

Stata uses a small sample adjustment for their GLM function, in which the standard sandwich variance (obtained via the vce(robust) command) is scaled according to the sample size and the number of parameters in the model. R uses no such adjustment. Thus, the specific relation between the sandwich SE from Stata and R (specifying HC1 when using the latter) is:

$$ SE_{sandwich, Stata} = \left (SE_{sandwich, R} \right ) \times \sqrt{ (n-k)/(n-1) } $$

where $n$ is the sample size, and $k$ is the number of parameters in the model. Thus, for the specific problem, the sandwich standard error estimate obtained in R was 10.11. Using the above equation with two parameters (intercept and asthma coefficient) and a sample size of 15, we can convert the R estimate (10.11) to it's corresponding Stata value (9.74) as:

$$ 10.11 \times \sqrt{(15-2)/(15-1)} = 9.74$$