Solved – How to estimate model predicted means (a.k.a. marginal means, lsmeans, or EM means) from a GEE model fitted in R

generalized-estimating-equationsr

I am trying to obtain model-predicted means and CI's for a categorical predictor in a GEE model fitted with the geeglm function (geepack package). The model is fitted with no problem, but where I am stuck is when trying to estimate the model-predicted group means. This is something fairly easy to do in other software packages (i.e., SAS and SPSS), but my point is to try to do this in R (also I was a bit disappointed to find no direct way to obtain overall tests for a categorical predictor other than fitting reduced and full models separately and then comparing them; but on the other hand I was able to find an easy way to compute the QIC with the MESS package).
Anyway, I searched around to see if perhaps this was done with another package (there's a package called lsmeans, but it seems not compatible with geepack or gee packages). I am surprised I have not been able to find something as common as estimating model predicted means for a GEE model in R, so I was wondering if someone knows a solution for this. I tried contacting the geepack maintainer, but no answer.

Why GEE instead of a mixed model (marginal means for a mixed model can be computed using the lmerTest package)? well, it has fewer assumptions, and is more robust with small samples.

Best Answer

The emmeans package now provides estimated marginal means for GEE models:

library(geepack)
library(emmeans)
warp.gee <- geeglm(breaks ~ tension, id=wool, family=gaussian, data=warpbreaks)
emmeans(warp.gee, ~tension)
 tension   lsmean       SE  df asymp.LCL asymp.UCL
 L       36.38889 5.774705 Inf  25.07067  47.70710
 M       26.38889 1.689200 Inf  23.07812  29.69966
 H       21.66667 2.042753 Inf  17.66294  25.67039

plot(emmeans(warp.gee, ~tension), horizontal=FALSE, ylab="Estimated mean")
Related Question