Post hoc test after GEE with repeated measures using R

Question

I'm working on repeated measures design GEE and would like to get the result of post-Hoc test, but I don`t know how to do it.

I just used those methods to get the result of post-Hoc test, but it seems wrong.

My data table looks like this:

id  times   y   dose    doTimes
1   1      200  250     3
1   2      300  250     3
1   3      280  210     3
1   4      280  125     3
2   1      248  254     3
2   2      345  148     3
2   3      2654.4 73    3
3   1      912  223     3
3   2      523.5 98     3
......

Then I got the fitted values:

gee <- geeglm(y ~ doTimes + dose, data=myData, id=id, family = gaussian,corstr = "exchangeable")
myData$fitted.Values <- gee$fitted.values
myData$fitted.Values


id  times   y   dose    doTimes fitted.Values
6   2   313.5   100      3      1036.6992
6   3   1101.8  100      3      1036.6992
6   4   2501.8  100      3      1036.6992
9   1   246     350      3      954.9299
9   2   676.5   350      3      954.9299
9   3   1125.6  350      3      954.9299
9   4   2123.8  350      3      954.9299
10  1   794     300      3      971.2838
10  2   598.5   300     3       971.2838
10  3   1625.4  300     3       971.2838
10  4   4379.2  300     3       971.2838
11  1   200     200     3       1003.9915

After that, I run TukeyHSD:

f <- aov(myData$fitted.Values ~ myData$doTimes)
HSD <- TukeyHSD(f, conf.level = 0.95)
HSD$`myData$doTimes`

         diff       lwr       upr              p adj
1-0  546.4474  535.9451  556.9497 0.0000000003426653
2-0  414.4079  405.7935  423.0224 0.0000000003426653
3-0  161.3648  150.0971  172.6325 0.0000000003426653
2-1 -132.0395 -140.8439 -123.2350 0.0000000003426653
3-1 -385.0826 -396.4962 -373.6690 0.0000000003426653
3-2 -253.0432 -262.7478 -243.3385 0.0000000003426653

The whole code looks like:

gee <- geeglm(y ~ doTimes + dose, data=myData, id=id, family = gaussian,corstr = "exchangeable")
myData$fitted.Values <- gee$fitted.values
f <- aov(myData$fitted.Values ~ myData$doTimes)
HSD <- TukeyHSD(f, conf.level = 0.95)
HSD$`myData$doTimes`

The value of p.adj seems wrong, I know maybe I made a very low-level mistake... but if anyone could help me, it would be appreciated! Thanks.

Answer 1

The p.values look odd because you are taking the fitted values and using this to perform an anova, which doesn't capture the variation in your data. If you are interested in doing a posthoc on the effect of doTimes, you can try using glht from multcomp:

library(multcomp)
library(geepack)
data(dietox)
#dietox$Weight = as.integer(dietox$Weight)
dietox$Cu <- as.factor(dietox$Cu)
mf <- formula(Weight ~ Cu + Time + I(Time^2) + I(Time^3))
gee1 <- geeglm(mf, data=dietox, id=Pig, corstr="ar1")

summary(glht(gee1,mcp("Cu" = "Tukey")))

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: geeglm(formula = mf, data = dietox, id = Pig, corstr = "ar1")

Linear Hypotheses:
                   Estimate Std. Error z value Pr(>|z|)
Cu035 - Cu000 == 0  -0.4467     1.4342  -0.311    0.947
Cu175 - Cu000 == 0   1.1821     1.7856   0.662    0.784
Cu175 - Cu035 == 0   1.6288     1.8353   0.888    0.646
(Adjusted p values reported -- single-step method)