Robust SE clustered GLM Gamma Log Link to match GEE Robust SE 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.
 A: 
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

