GEE vs mixed model for time-varying covariate

Question

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.

Answer 1

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.

Answer 2

library(tidyverse)
library(geepack)
library(emmeans)
library(MESS)
#> Loading required package: geeM
#> Loading required package: Matrix
#> 
#> 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, 
                    family=binomial(link="logit"),
                    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)}