How can I estimate model predicted means (a.k.a. marginal means, lsmeans, or EM means) from a GEE model fitted in R? 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.  
 A: The LSmeans function in the doBy package may be helpful.
Here is a simple modification of an example in the vignette.
library(doBy)
library(geepack)
warp.gee <- geeglm(breaks ~ tension, id=wool, family=gaussian, data=warpbreaks)
LSmeans(warp.gee,effect="tension")

A: 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")

A: The spind package offers a predict function similar to other uses of predict:
https://www.rdocumentation.org/packages/spind/versions/2.1.3/topics/predict.GEE
# load packages
library(geepack)
library(spind)

n <- nrow(warpbreaks) # number of cases
trainIndex <- sample(1:n, n*.60) # random subset of 60% of cases (training)
testIndex <- setdiff(1:n, trainIndex) # rest of cases not in training (testing)

# model fit
warp.gee.fit <- geeglm(breaks ~ tension, id=wool, family=gaussian, data=warpbreaks[trainIndex,])

# predict test cases
predict(warp.gee.fit, warpbreaks[testIndex,])
```

