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.

#> 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)
#> 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).

  • $\begingroup$ Seems to me this question is posted in 2 places. You should post in only one. I answered it elsewhere. $\endgroup$ – rvl Feb 9 at 17:47
  • $\begingroup$ 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. $\endgroup$ – Frederick Feb 13 at 9:19

Your Answer

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

Browse other questions tagged or ask your own question.