# How to estimate model predicted means per group from a GEE model fitted in R?

This question is related to a similar post: How can I estimate model predicted means (a.k.a. marginal means, lsmeans, or EM means) from a GEE model fitted in R?

The difference is that I do not only fit categorical predictors but also continuous. This problem could also be a pure coding issue but I would like to use this question to cross-validate my approach as well. Hence, my problem is twofold.

Part 1
First, I fitted a GEE for longitudinal data with multiple observations per ID. In R I use the geepack::geeglm function and update my model with all possible corstr except for "userdefined". Following this post: QIC or QICu for variable selection in GEE-GLM variable selection I compute the QIC of every model and choose the model with the lowest QIC.

Is this a valid approach of selecting a GEE model?

Part 2
Second, and most important to me: I have fitted a GEE using the approach above. Explanatory variables are one categorical variable (i.e. treatment 1/0 [yes/no]) and one continues (between 0 and 1).

Fitting the model is no problem. But now I want to find the estimated mean per group (that is per treatment group).

To calculate the estimated means I am using emmeans::emmeans()because this also returns both LCL and UCL. What I don't understand is: why does emmeans()is using the overall observed mean instead of the group mean? How can I make sure to compute the estimated group mean (including LCL/UCL)?

I tried to make a minimal reproducible example below. Hope someone can help.

library("geepack")
library("emmeans")
library("dplyr")
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#>     filter, lag
#> The following objects are masked from 'package:base':
#>
#>     intersect, setdiff, setequal, union

# Extending to pigs data to create a fourth variable representing the
#treatment group
pigs.extended <- pigs %>% mutate(trt = c(rep("a", 20), rep("b", 9)))

# The observed means for percent by treatment group
pigs.extended %>% group_by(trt) %>% summarise(mean = mean(percent))
#> # A tibble: 2 x 2
#>   trt    mean
#>   <chr> <dbl>
#> 1 a      13.2
#> 2 b      12.3

# The overall mean (wihtout grouping)
pigs.extended %>% summarise(mean = mean(percent))
#>       mean
#> 1 12.93103

# A simple example model
model <- geepack::geeglm(conc ~ percent + trt, id = source, data = pigs.extended)
model
#>
#> Call:
#> geepack::geeglm(formula = conc ~ percent + trt, data = pigs.extended,
#>     id = source)
#>
#> Coefficients:
#> (Intercept)     percent        trtb
#>   20.498948    1.049322   10.703857
#>
#> Degrees of Freedom: 29 Total (i.e. Null);  26 Residual
#>
#> Scale Link:                   identity
#> Estimated Scale Parameters:  [1] 36.67949
#>
#> Correlation:  Structure = independence
#> Number of clusters:   3   Maximum cluster size: 10

# Estimating (overall) mean with emmeans
emmeans(model, "percent")
#>  percent emmean   SE  df asymp.LCL asymp.UCL
#>     12.9   39.4 1.65 Inf      36.2      42.6
#>
#> Results are averaged over the levels of: trt
#> Covariance estimate used: vbeta
#> Confidence level used: 0.95

# why is the percent mean for each group 12.9 and not 13.2 and 12.3?
emmeans(model, "percent", by = "trt")
#> trt = a:
#>  percent emmean    SE  df asymp.LCL asymp.UCL
#>     12.9   34.1 3.252 Inf      27.7      40.4
#>
#> trt = b:
#>  percent emmean    SE  df asymp.LCL asymp.UCL
#>     12.9   44.8 0.327 Inf      44.1      45.4
#>
#> Covariance estimate used: vbeta
#> Confidence level used: 0.95


Created on 2019-02-07 by the reprex package (v0.2.0).

• Seems to me this question is posted in 2 places. You should post in only one. I answered it elsewhere. – rvl Feb 9 at 17:47
• The second part is indeed posted on stackoverflow in a similar form since I think that this could be a pure coding issue. Part 1 belongs here in in my opinion. – Frederick Feb 13 at 9:19