GEE Reference

library(geepack)

# Binary — exchangeable correlation
fit <- geeglm(y ~ time + x, id = id, data = df,
              family = binomial, corstr = "exchangeable")

# Binary — AR(1)
fit <- geeglm(y ~ time + x, id = id, data = df,
              family = binomial, corstr = "ar1")

# Binary — independence
fit <- geeglm(y ~ time + x, id = id, data = df,
              family = binomial, corstr = "independence")

# Binary — unstructured
fit <- geeglm(y ~ time + x, id = id, data = df,
              family = binomial, corstr = "unstructured")

# Count — Poisson with exchangeable
fit <- geeglm(y ~ time + x, id = id, data = df,
              family = poisson, corstr = "exchangeable")

# Time × treatment
fit <- geeglm(y ~ time * treatment, id = id, data = df,
              family = binomial, corstr = "exchangeable")

# Polynomial time
fit <- geeglm(y ~ time + I(time^2) + treatment, id = id, data = df,
              family = binomial, corstr = "exchangeable")

fit <- geeglm(events ~ time + treatment, id = id, data = df,
              family = poisson, corstr = "exchangeable",
              offset = log(exposure_time))

Independence:       Exchangeable:       AR(1):              Unstructured:
[1 0 0 0]          [1 α α α]          [1  α  α² α³]       [1   r₁₂ r₁₃ r₁₄]
[0 1 0 0]          [α 1 α α]          [α  1  α  α²]       [r₁₂ 1   r₂₃ r₂₄]
[0 0 1 0]          [α α 1 α]          [α² α  1  α ]       [r₁₃ r₂₃ 1   r₃₄]
[0 0 0 1]          [α α α 1]          [α³ α² α  1 ]       [r₁₄ r₂₄ r₃₄ 1  ]

                    Start here
                        │
            ┌───────────┴───────────┐
            │  Few time points      │
            │  (≤ 5) and ≥ 100      │
            │  clusters?            │
            └───────┬───────┬───────┘
                Yes │       │ No
                    ▼       ▼
            ┌───────────┐ ┌─────────────────┐
            │ Consider  │ │ Exchangeable     │
            │ Unstruc-  │ │ (safe default)   │
            │ tured     │ └────────┬─────────┘
            └───────────┘          │
                          ┌────────┴────────┐
                          │ Expect temporal  │
                          │ decay?           │
                          └───┬─────────┬───┘
                          Yes │         │ No
                              ▼         ▼
                     ┌──────────┐  ┌──────────────┐
                     │ Try AR(1)│  │ Stay with    │
                     │ + compare│  │ exchangeable │
                     │ QIC      │  └──────────────┘
                     └──────────┘

# Correlation parameter(s)
fit$geese$alpha

# Full estimated working correlation matrix
# (not directly available; reconstruct from alpha and corstr)

coef(fit)
summary(fit)$coefficients

# ORs (binary) or IRRs (count)
exp(coef(fit))

# With confidence intervals
exp(confint(fit))

# Or manually from robust SEs
beta <- coef(fit)
se   <- summary(fit)$coefficients[, "Std.err"]
cbind(OR = exp(beta),
      lower = exp(beta - 1.96 * se),
      upper = exp(beta + 1.96 * se))

# Robust (sandwich) — reported by default in geepack
summary(fit)$coefficients[, "Std.err"]

# Model-based (naive) — only valid if corstr is correct
sqrt(diag(fit$geese$vbeta.naiv))

# Marginal predictions
predict(fit, type = "response")

# For new data
predict(fit, newdata = nd, type = "response")

# Link-scale predictions
predict(fit, type = "link")

fit$geese$alpha

fit$geese$gamma

# geepack provides QIC
QIC(fit_exch)
QIC(fit_ar1)
QIC(fit_indep)

# Compare in a table
data.frame(
  corstr = c("Exchangeable", "AR(1)", "Independence"),
  QIC = c(QIC(fit_exch)[1], QIC(fit_ar1)[1], QIC(fit_indep)[1])
)

robust <- summary(fit)$coefficients[, "Std.err"]
naive  <- sqrt(diag(fit$geese$vbeta.naiv))

data.frame(
  term  = names(coef(fit)),
  robust = round(robust, 4),
  naive  = round(naive, 4),
  ratio  = round(robust / naive, 3)
)

# Compare robust vs naive SEs (large discrepancy = misspecification)
robust <- summary(fit)$coefficients[, "Std.err"]
naive  <- sqrt(diag(fit$geese$vbeta.naiv))
max(abs(robust / naive - 1))  # Should be < 0.2

# Pearson residuals
res <- residuals(fit, type = "pearson")

# Plot residuals vs fitted
plot(fitted(fit), res, xlab = "Fitted", ylab = "Pearson Residuals")
abline(h = 0, lty = 2)

# Dispersion parameter (should be near 1 for Poisson)
fit$geese$gamma

# Cook's distance analogue — identify influential clusters
# (Not built into geepack; use leave-one-cluster-out)
ids <- unique(df$id)
betas <- matrix(NA, length(ids), length(coef(fit)))
for (i in seq_along(ids)) {
  fit_i <- update(fit, data = df[df$id != ids[i], ])
  betas[i, ] <- coef(fit_i)
}
# Check for outlier clusters
apply(betas, 2, sd)

df <- df %>% arrange(id, time)

# Switch from unstructured to exchangeable or AR(1)
fit <- geeglm(y ~ time + x, id = id, data = df,
              family = binomial, corstr = "exchangeable")

# Option 1: Use GLMM instead (handles MAR via FIML)
fit_glmm <- glmer(y ~ time + x + (1 | id), data = df, family = binomial)

# Option 2: Multiple imputation + GEE (requires additional packages)

OR <- exp(coef(fit)["treatment"])

IRR <- exp(coef(fit)["time"])

beta <- coef(fit)["treatment"]
se   <- summary(fit)$coefficients["treatment", "Std.err"]
exp(beta + c(-1.96, 1.96) * se)

# Given GLMM random intercept variance sigma2
c_star <- 1 / sqrt(1 + (16 * sqrt(3) / (15 * pi))^2 * sigma2)
marginal_beta <- conditional_beta * c_star

fit <- geeglm(y ~ time + current_stress + treatment, id = id,
              data = df, family = binomial, corstr = "exchangeable")

# Time × group interaction (different trends by group)
fit <- geeglm(y ~ time * treatment, id = id, data = df,
              family = binomial, corstr = "exchangeable")

fit <- geeglm(count ~ time + treatment, id = id, data = df,
              family = poisson, corstr = "exchangeable")
exp(coef(fit))  # IRRs

# Ordinal GEE
fit <- ordgee(ordered(y) ~ time + treatment, id = id, data = df,
              corstr = "exchangeable")