GEE vs mixed model for time-varying covariate

Assuming the attached dataset, I am looking to examine whether participants who are treated are more/less likely to have high addiction severity compared to non-treated.

Both treatment and high addiction severity are time-varying (ie, re-assessed at each visit).

Given the correlated nature of the data, I think GEE would be a good fit. If my understanding in right, this approach would tell me if, on average, treated people are more/less likely to have high addiction severity compared to non-treated people.

However, someone has suggested that I use a mixed model to explore differences among treated and non-treated individuals with respect to this covariate (and other ones).

I have read a bit about the distinction between the 2 methods and understand that GEE gives us population-average effects whereas a mixed model would give us subject-specific effects (between-subject differences are treated as random). This distinction makes sense to me in the context of a time-fixed covariate but cannot really understand how it would work in the context of a time-varying covariate …Would this analysis make sense, and if so, how would the results be different compared to a GEE?

There is an issue of temporality/reverse causality here but please ignore for now as I am trying to understand the distinction between the 2 methods.

Thanks so much in advance.

• See my answer here and the referenced work stats.stackexchange.com/questions/163216/… GEE will work with time-varying covariates. This article may shed some light: academic.oup.com/biostatistics/article/6/4/633/228809 – StatsStudent Jan 27 at 9:08
• The four subjects for whom data is shown each have 3, 4, 2 and 3 visits, respectively. If you are going to have a "time" factor then Dimitris Rizopoulos correctly points out that GEE is problematic with missing values (presumably the 3rd, 4th and/or 5th visits are missing in the example data). If OTOH your design simply allows for a variable number of visits then I'm pretty sure you can use GEE which does allow variable cluster sizes. Each subject is a cluster with anywhere from 2 to 5 (or whatever the range is in the full data) observations. – Brent Hutto Jan 30 at 11:20

Indeed, for dichotomous outcomes, as you seem to have here, the corresponding mixed effects model, namely a mixed effects logistic regression gives you fixed effects coefficients that have an intepretation conditional on the random effects. A detailed explanation can be found here. Most often, this is not the interpretation you want. The GEE approach does give you coefficients with a marginal / population-averaged interpretation.

However, an additional practical point that you also need to consider is missing data. You have not given us enough details with regard to this point for your application, but almost always we have to deal with incomplete data. With regard to this point, mixed models give you valid results under the less stringent missing at random assmuption compared to the (standard not weighted) GEE that give you valid results under the less realistic missing completely at random assumption.

Taking both points (i.e., interpretation and missing data) into account, you would most often like to fit a mixed model to be more protected for the missing data but want to obtain parameters that have a population averaged intepretation. An early solution towards this direction was the marginalized mixed models propoposed by Heagerty, but, in general, these are computionally intensive to fit. A more recent approach that seems to solve the problem has been proposed by Hedeker et al. This is implemented in the function marginal_coefs() of the R package GLMMadaptive. You can find an example on how to use this function here.

• Thank you for your answer! I think what I am wondering, first of all, is whether a mixed model would answer the same question, on a conceptual level? The one mentioned in my post above. I am a bit familiar with GEE and thought that it would be the best fit here, so I wonder why/if a mixed model would be the better approach? – blackandwhite Jan 31 at 15:33
• Marginalized mixed models are suitable to answer the same question as GEEs, protecting you at the same time for missing at random missing data. – Dimitris Rizopoulos Jan 31 at 19:20
library(tidyverse)
library(geepack)
library(emmeans)
library(MESS)
#>
#> Attaching package: 'Matrix'
#> The following object is masked from 'package:tidyr':
#>
#>     expand

unlogit <- function(x) {exp(x)/(1+exp(x))}

# Simulate a dataset with 40 subjects, each with 2-5 visits
# where p(outcome)=0.6 in untreated, p(outcome)=0.2 in treated

set.seed(279)

subj <- c(rep(1:10, each=2), rep(11:20, each=3), rep(21:30, each=4), rep(31:40, each=5))
visit <- c(rep(1:2, 10), rep(1:3, 10), rep(1:4, 10), rep(1:5, 10))
unique <- runif(40,0,1)
dat <- data.frame(subj,visit) %>% mutate(trt=rbinom(n(), 1, unique[subj]),
outcome=rbinom(n(), 1, 0.6-(0.4*trt)))

# Fit the simplest possible GEE model, treating "visit" as a nuisance clustering effect

fit.exch <- geeglm(outcome ~ trt, data=dat,
id=subj, corstr = "exchangeable", std.err="san.se")
summary(fit.exch)
#>
#> Call:
#> geeglm(formula = outcome ~ trt, family = binomial(link = "logit"),
#>     data = dat, id = subj, corstr = "exchangeable", std.err = "san.se")
#>
#>  Coefficients:
#>              Estimate   Std.err  Wald Pr(>|W|)
#> (Intercept)  0.001003  0.207653 0.000   0.9961
#> trt         -0.888043  0.358173 6.147   0.0132 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Estimated Scale Parameters:
#>             Estimate Std.err
#> (Intercept)   0.9999 0.07357
#>
#> Correlation: Structure = exchangeable  Link = identity
#>
#> Estimated Correlation Parameters:
#>       Estimate Std.err
#> alpha -0.01807 0.07048
#> Number of clusters:   40   Maximum cluster size: 5
QIC(fit.exch)
#>       QIC      QICu Quasi Lik       CIC    params      QICC
#>   185.220   185.192   -90.596     2.014     2.000   185.396

#Estimate p(outcome) in untreated
unlogit(summary(fit.exch)\$coef["(Intercept)", "Estimate"])
#> [1] 0.5003

#Estimate p(outcome) in treated
unlogit(summary(fit.exch)$$coef["(Intercept)", "Estimate"] + summary(fit.exch)$$coef["trt", "Estimate"])
#> [1] 0.2917


Created on 2019-01-30 by the reprex package (v0.2.1)