I am working with a pre/post test structure, measuring dichotomous outcomes. I am using GEE to estimate the coefficient for time (also a binary variable, 0 representing pre and 1 representing post), given that most people participate in only one of the time points, but there are a number of people who are represented in both.

I am interested in performing a Monte Carlo simulation as a "sensitivity test" to assess how many dependent data points there need to be before ignoring dependencies in data causes a huge problem in inference (e.g., 3% of the data points are from people in both the pre and post data? 10%? 50%?).

Thus, I am interested in simulating clustered data that are meant to be estimated using GEE, but I will compare the GEE coefficients and standard errors to the regular maximum likelihood generalized linear model (GLM) estimates.

However, I am unsure how to simulate data for a GEE. I know the data take the same functional form as in a traditional GLM:

$$g(\mu_{ij}) = X'_{ij}\beta$$

where $g()$, in my case, is the logistic link function. I also know that GEE is basically a different estimation method (often contrasted against maximum likelihood) combined with a method of estimating robust variances (described well in this answer: https://stats.stackexchange.com/a/62924/130869).

We also specify a error correlation structure, which I am assuming as exchangeable. I do not know how to simulate this, as binary data using the logistic link function do not have an error term in the equation.

Most examples of how to generate correlated data assume a parametric form of a traditional multilevel model, where I would simulate something like this (note the lack of error term):

$$\text{logit}(\pi_{ij}) = \beta_{0j}$$ and $$\beta_{0j} = \gamma_{00} + \gamma_{01}X_j + u_{0j}$$

This would involve simulating an average intercept, an average treatment effect, and a subject-specific deviation from that intercept (from a normal distribution, which a GEE does not assume).

However, is there a way to simulate my data that are more explicitly linked to GEE? I would like to specify a range of correlations that represent the true correlation in the exchangeable error structure, the proportion of data points that are dependent, and the average effect of time, in line with the assumptions of a generalized estimating equation.

Update, taking into account @IsabellaGhement's suggestion in the comments below. I tried using the simstudy package, but the observations aren't acting as if they are dependent when throwing them in a regression. Below, I simulate paired data and then run a GEE and GLM on them. The distributions of the estimates are the same, whereas I've seen previous simulations showing that a consequence of non-independence is a wider distribution of parameters. Relatedly, I get better 95% CI coverage with a GLM rather than a GEE; again, I've seen simulations where dependence of observations harms coverage for GLM models, which assume iid. What am I missing in trying to simulate dependent observations? The same type of thing occurred when I used mvtnorm::rmvnorm() and dichotomized the variables.

logit <- function(p) log(p / (1 - p))
results <- lapply(1:2000, function(zzz) {
  dat <- genCorGen(n = 1000, nvars = 2, params1 = c(.50, .53), 
                   dist = "binary", rho = .80, corstr = "cs", wide = FALSE)
  b1_pop <- logit(.53) # leveraging the fact that logit(.50) is 0
  gee_mod <- geeglm(X ~ period, binomial, dat, id = id, corstr = "exchangeable")
  gee_ub <- summary(gee_mod)$coef[2, 1] * 1.96 + summary(gee_mod)$coef[2, 2]
  gee_lb <- summary(gee_mod)$coef[2, 1] * 1.96 - summary(gee_mod)$coef[2, 2]
  gee_cover <- b1_pop < gee_ub & b1_pop > gee_lb
  glm_mod <- glm(X ~ period, binomial, dat)
  glm_ub <- summary(glm_mod)$coef[2, 1] * 1.96 + summary(glm_mod)$coef[2, 2]
  glm_lb <- summary(glm_mod)$coef[2, 1] * 1.96 - summary(glm_mod)$coef[2, 2]
  glm_cover <- b1_pop < glm_ub & b1_pop > glm_lb

  c(b1_gee = coef(gee_mod)[[2]], b1_glm = coef(glm_mod)[[2]],
    gee_cover = gee_cover, glm_cover = glm_cover)
results <- as.data.frame(do.call(rbind, results))
results %>% 
  gather() %>% 
  ggplot(aes(x = value, fill = key)) +
  geom_density() +
  facet_wrap(~ key)

The call to colMeans returns:

   b1_gee    b1_glm gee_cover glm_cover 
    0.120     0.120     0.243     0.366 

Indicating that the mean parameter estimate for both GEE and GLM were the population parameter, and the GLM had 37% coverage, while the GEE had 24%.

The call to ggplot shows that the distribution of parameter estimates were essentially the same:

enter image description here

  • 1
    $\begingroup$ Mark, see if this link will help you get started: rdatagen.net/post/intro-to-simstudy. Then check out this link: rdatagen.net/post/simstudy-update-two-functions-for-correlation for an example of generating correlated (repeated) binary data in R. $\endgroup$ Commented Aug 19, 2018 at 1:36
  • $\begingroup$ As far as I know, simulating data from non-Normal GEEs is extremely difficult or even impossible because these models are based on quasi-likelihoods, not actual likelihoods. $\endgroup$ Commented Aug 22, 2018 at 3:40
  • $\begingroup$ @JakeWestfall makes sense. Even still, GEE should be outperforming GLM in the simulation I described in the post, no? It makes it look like ignoring the dependencies between observations is the better analysis option. $\endgroup$
    – Mark White
    Commented Aug 22, 2018 at 13:12
  • $\begingroup$ Why is the coverage so low for both methods? Can you show us the R code for mvtnorm::rmvnorm() and dichotomization of the variables? I was surprised to read that the discritization method did not work, considering that it allows you to specify the dependence structure you want, which I thought was the objective. Also, I didn't see X and period defined anywhere. $\endgroup$
    – papgeo
    Commented Aug 24, 2018 at 21:39

1 Answer 1


I just came across this post and this may be too late, but what you did was complete correct except that your formulae for calculating CIs was wrong. The correct formulae should be:

gee_ub <- summary(gee_mod)$coef[2, 1] + 1.96 * summary(gee_mod)$coef[2, 2]

gee_lb <- summary(gee_mod)$coef[2, 1] - 1.96 * summary(gee_mod)$coef[2, 2]

and same for GLM. I tried your codes with above correction and got below:

   b1_gee    b1_glm gee_cover glm_cover 
    0.120     0.120     0.948     0.998 

It shows GEE is correct while GLM results in too small p values thus coverage rate too high, which should be the case as your simulated data is correlated.

  • $\begingroup$ Aha! Simple coding error it was after all. Thank you! $\endgroup$
    – Mark White
    Commented Mar 10, 2021 at 19:04

Your Answer

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