2
$\begingroup$

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.

$\endgroup$

1 Answer 1

1
$\begingroup$

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)
$\endgroup$
1
  • $\begingroup$ Thank you for this detailed answer, I have used glht on my GEE, but it`s prompted a wrong message, so I used to think that glht dose not fit GEE, because of that, I took a lot of detours. Now I can go on my work. $\endgroup$
    – dbcoffee
    May 2, 2020 at 4:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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