1

How do I get the robust standard errors/sandwich variance estimators for GLM using a Gamma family with a log-link to match the robust standard errors from the GEE output?

library(sandwich)
library(lmtest)
data('CO2')

up_glm <- glm(uptake ~ Type, data = CO2, family = Gamma(link = 'log'))
up_gee <- gee::gee(uptake ~ Type , data = CO2, id = Plant, family = Gamma(link = 'log'))
up_gee_glm <- geepack::geeglm(uptake ~ Type , data = CO2, id = Plant, family = Gamma(link = 'log'))


summary(up_glm)$coefficients[,1:2]
coeftest(up_glm, vcov = vcovCL, cluster = CO2$Plant)
summary(up_gee)$coefficients[,c(1,2,4)]
summary(up_gee_glm)$coefficients[,1:2]

I've tried to follow this post and these directions but neither give me the robust SE that match the GEE output.

renethestudent
  • 13

1 Answers1

0

I've tried to follow this post and these directions but neither give me the correct robust SE.

I figure you are missing the cadjust argument. You can either create the robust covariance matrix yourself or pass the cluster and cadjust = FALSE argument to lmtest::coeftest (see help(coeftest)) like below:

library(sandwich)
library(lmtest)
options(digits = 3)
data('CO2')

up_glm <- glm(uptake ~ Type, data = CO2, family = Gamma(link = 'log'))
V <- vcovCL(up_glm, cluster = ~ Plant, cadjust = FALSE)


manually


co <- coef(up_glm)
cbind(Estimate = co, 
      Std. Error = sqrt(diag(V)), 
      t value = co / sqrt(diag(V)), 
      Pr(&gt;|z|) = 
          2 * pnorm(-abs(co / sqrt(diag(V)))))
#R>                 Estimate Std. Error t value Pr(>|z|)
#R> (Intercept)        3.513     0.0288  121.85 0.00e+00
#R> TypeMississippi   -0.474     0.1109   -4.27 1.94e-05


w/ coeftest


coeftest(up_glm, vcovCL, cluster = CO2$Plant, cadjust = FALSE)
#R> 
#R> z test of coefficients:
#R> 
#R>                 Estimate Std. Error z value Pr(>|z|)

#R> (Intercept)       3.5128     0.0288  121.85  < 2e-16 ***
#R> TypeMississippi  -0.4739     0.1109   -4.27  1.9e-05 ***
#R> ---
#R> Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Then you get the same output as from gee::gee

library(gee)
summary(
  gee(uptake ~ Type , data = CO2, id = Plant, 
      family = Gamma(link = 'log')))$coefficients
#R>     (Intercept) TypeMississippi 
#R>           3.513          -0.474 
#R>                 Estimate Naive S.E. Naive z Robust S.E. Robust z
#R> (Intercept)        3.513     0.0516    68.1      0.0288   121.85
#R> TypeMississippi   -0.474     0.0729    -6.5      0.1109    -4.27
Benjamin Christoffersen
  • 2,466
  • Thanks @benjamin, I have the cluster statement in there and even tried to compute it by hand. I get the same numbers as you. I guess my problem is that I want to match the robust SE from the gee outputs; the robust SEs from these outputs match to SAS and STATA. Do you have any suggestions for that? – renethestudent Jul 07 '20 at 12:56
  • It was not clear from your question that you wanted the same standard errors that you get from gee. You need to use the cadjust argument. – Benjamin Christoffersen Jul 07 '20 at 14:34
  • To clarify, do you want the same robust covariance matrix as gee::gee gives you or those that you get from SAS and STATA? This is not clear from your post or from the tile of your post. Could you edit your post to make this clear? – Benjamin Christoffersen Jul 07 '20 at 14:40
  • 1
    My apologies, I updated it to reflect that I would like the SE of the GLM to match the robust SE of the GEE outputs. And for clarification, the robust SE of the GEE outputs already match the robust SE outputs from Stata and SAS, so I'd like the GLM robust SE to match it. – renethestudent Jul 07 '20 at 16:51