Can GEEs handle the Simpson's paradox? I'm a somewhat seasoned user of mixed effects/multilevel models but I don't know much about Generalized Estimating Equations (GEEs). Over time, I've been seen many people claim that GEEs can be used to handle the same kind of correlated/clustered data that mixed effects models are often used for.
However, one thing that always puzzled me is: can GEEs handle the correlated data if the trend across clusters is opposite to the trend within clusters? I.e. the situation that is often cited as one of the examples of Simpson's paradox?
I've simulated some data like this in R:
set.seed(123456)
beta_s <- rnorm(10, 0, 0.5)
alpha_s <- rnorm(10, 0, 15)
x <- round(runif(100, 1, 10), 2)

df <- data.frame(id = rep(1:10, each = 10),
             alpha = 100, beta = -1.25, 
             alpha_s = rep(alpha_s, each = 10),
             beta_s = rep(beta_s, each = 10))

df <- transform(df, x = rnorm(100, 0.2 * alpha_s, 1))
df <- transform(df, y = alpha + alpha_s + (beta + beta_s) * x + rnorm(100))

ggplot2::ggplot(df, aes(x, y, col = factor(id))) +
  geom_point()


When I model this data using a mixed effects model fitted with cmdstanr, I can recover the negative beta parameter very well. However, when I try fitting it with geepack::geeglm(y ~ x, id = id, data = df), I get a strong positive estimate for the slope of x (that's also highly statistically significant!).
Am I perhaps specifying the model incorrectly or is GEE the wrong sort of tool for this job? Can we use GEEs to estimate within-cluster trends at all?

The default in geepack::geelm() is to have an independence working correlation matrix. If you alter this, you obtain similar results, i.e., geeglm(y ~ x, data = df, id = id, corstr = "ar1").
