I am doing a randomised controlled trial and I need to use a GEE model for my statistical analysis.

I have two groups (intervention and control) and I collect injured subjects every month (during 7 months) across a whole season (either injured or not. i.e. injury prevalence), and in the end I want to compare both groups throughout the season. I have read two papers that use GEE with similar study designs compared to my study (Harøy et al. 2017: https://doi.org/10.1136/bjsports-2017-098937 and Åkerlund et al. 2022: https://doi.org/10.1007/s00167-021-06644-2), although they use SPSS for statistical analysis, also idk how to perform GEE there.

I started reading the documents from the "gee" package and end up running this code:

migee <- gee(Injury ~ GRUPO * Tiempo, family = poisson(link = "log"), 
           data = GEE_trial, corstr = "exchangeable", id = id)

"Injury" is a dichotomous variable that means either injured player or not "GRUPO" is the group, either control or intervention "Tiempo" is each month GEE_trial is the dataframe "id" is a number assigned to each subject

I used log poisson and an exchangeable correlation structure since I read in the aforementioned papers that they used this format. Although I am not sure whether this is properly performed or not.

After running the analysis and the "summary" function, various coefficients appeared (note that i only run the analysis with data from three months instead of the entire season) with standard errors and z-scores. Coefficients:

                 Estimate Naive S.E.    Naive z Robust S.E.   Robust z
(Intercept)    -2.1400662  0.6492832 -3.2960444   0.6642112 -3.2219666
GRUPO1          1.1415373  0.7362179  1.5505427   0.7289733  1.5659522
Tiempo2         0.6931472  0.7952062  0.8716571   0.7952062  0.8716571
Tiempo3         0.6931472  0.7952062  0.8716571   0.7952062  0.8716571
GRUPO1:Tiempo2 -1.9459101  1.0836834 -1.7956445   1.0816408 -1.7990354
GRUPO1:Tiempo3 -1.5404450  1.0167826 -1.5150191   1.0016388 -1.5379247          

I read that by exponentiating the "GRUPO1" coefficient, I would get the Odds ratio, since this is the effect size that I need to present in the paper. However, in the mentioned papers they present the Odds ratio (or Rate ratio (Åkerlund 2022) with the confidence interval.

First, I have no idea whether I have performed the analysis properly. Second, if correct, how can i calculate the confidence interval for this effect size (Odds ratio or Rate ratio)?

Third, if incorrect, how can I perform the analysis correctly?


1 Answer 1


We can't help with the "is this correct" because that is predicated on the scientific question at hand and something only you and subject-matter experts can determine. If you are naive to statistics and their role in your analysis, I suggest you consult a statistician within the institution from where you are working on your PhD to guide you. These types of discussion are also off-topic for stack overflow (you may consider Cross Validated)

As far as the objective, coding-related question goes, I dont believe gee objects can use typical functions to gather the confidence intervals (ie, confint()). However, you can calculate the Wald confidence interval directly using the results of the model.

You didnt provide data, so using this reproducible data:

n <- 1e3
GEE_trial <- data.frame(Injury = sample(0:1, n, replace = TRUE),
                 id = factor(LETTERS[1:5]),
                 GRUPO =  factor(sample(0:1, n, replace = TRUE)),
                 Tiempo = factor(sample(1:3,  n, replace = TRUE)))

You can run the model (migee), extract the results (mdlresults), then calculate the confidence intervals and organize the results in a data frame (finalres):

# model
migee <- gee::gee(Injury ~ GRUPO * Tiempo,
                  family = poisson(link = "log"),
                  data = GEE_trial,
                  corstr = "exchangeable",
                  id = id)

# 95% Wald confidence interval
ci <- 0.95
xx <- qnorm((1 + ci) / 2)

# final results
finalres <- data.frame(param = rownames(mdlresults),
                       est = mdlresults$Estimate,
                       lwr =  mdlresults$Estimate - xx * mdlresults$Robust.S.E.,
                       upr = mdlresults$Estimate + xx * mdlresults$Robust.S.E.)


# > finalres
#            param         est         lwr         upr
# 1    (Intercept) -0.79094492 -0.96054071 -0.62134914
# 2         GRUPO1  0.14851284 -0.07233009  0.36935577
# 3        Tiempo2  0.07898841 -0.15237451  0.31035133
# 4        Tiempo3  0.12337459 -0.09494036  0.34168953
# 5 GRUPO1:Tiempo2 -0.24411386 -0.56936113  0.08113342
# 6 GRUPO1:Tiempo3 -0.12646163 -0.42464761  0.17172435

Your Answer

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

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