5
$\begingroup$

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.

Improve this question
$\endgroup$
3
  • $\begingroup$ I've ran into this problem since I don't think predict methods are defined for gee objects. A quick cheat for getting around this is just to fit the analogous lm or glm model. The coefficients are exactly the same, the standard errors are what differ between the two. $\endgroup$
    – AdamO
    Jul 25, 2014 at 20:52
  • $\begingroup$ Also, one of the key relaxed assumptions of GEE is that of possible heteroscedasticity. Plugin estimates of residual variance are used to compute standard errors for linear models. So post-hoc tests of specific contrasts cannot necessarily be computed because you have to re-introduce assumptions about the distribution of errors around predicted values. $\endgroup$
    – AdamO
    Jul 25, 2014 at 21:00
  • $\begingroup$ FWIW, I am adding support for gee, geeglm, and geese objects in the next update of the lsmeans package. It'll be on CRAN in a couple of weeks. $\endgroup$
    – Russ Lenth
    Aug 30, 2014 at 15:02

3 Answers 3

Reset to default
3
$\begingroup$

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")
Improve this answer
$\endgroup$
5
$\begingroup$

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")
Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Thanks, it works great. Did not know doBy had those capabilities. $\endgroup$
    – user16263
    Jul 28, 2014 at 17:00
0
$\begingroup$

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,])
```
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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