2
$\begingroup$

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.

Improve this question
$\endgroup$
0

1 Answer 1

Reset to default
1
$\begingroup$

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(>|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
Improve this answer
$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – renethestudent
    Commented Jul 7, 2020 at 12:56
  • $\begingroup$ 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. $\endgroup$
    – Benjamin Christoffersen
    Commented Jul 7, 2020 at 14:34
  • $\begingroup$ 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? $\endgroup$
    – Benjamin Christoffersen
    Commented Jul 7, 2020 at 14:40
  • 1
    $\begingroup$ 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. $\endgroup$
    – renethestudent
    Commented Jul 7, 2020 at 16:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.