Common Working Correlations

Picture First: Three Correlation Shapes

Exchangeable Correlation

\[ \operatorname{Corr}_i(\alpha) = \begin{pmatrix} 1 & \alpha & \alpha & \alpha \\ \alpha & 1 & \alpha & \alpha \\ \alpha & \alpha & 1 & \alpha \\ \alpha & \alpha & \alpha & 1 \end{pmatrix} \]

Use when: No natural time ordering; cluster members similar.

AR(1) Correlation

\[ \operatorname{Corr}_i(\alpha) = \begin{pmatrix} 1 & \alpha & \alpha^2 & \alpha^3 \\ \alpha & 1 & \alpha & \alpha^2 \\ \alpha^2 & \alpha & 1 & \alpha \\ \alpha^3 & \alpha^2 & \alpha & 1 \end{pmatrix} \]

Use when: Correlation decays with time lag.

Unstructured Correlation

\[ \operatorname{Corr}_i(\alpha) = \begin{pmatrix} 1 & \alpha_{12} & \alpha_{13} & \alpha_{14} \\ \alpha_{12} & 1 & \alpha_{23} & \alpha_{24} \\ \alpha_{13} & \alpha_{23} & 1 & \alpha_{34} \\ \alpha_{14} & \alpha_{24} & \alpha_{34} & 1 \end{pmatrix} \]

Use when: Small \(n_i\), enough subjects, no assumed pattern.

GEE with Different Correlations

# Independence
gee_ind <- geeglm(
  seizures ~ trt + log_base + log_age + visit,
  id = id,
  data = epi,
  family = poisson(link = "log"),
  corstr = "independence"
)

# Exchangeable
gee_exch <- geeglm(
  seizures ~ trt + log_base + log_age + visit,
  id = id,
  data = epi,
  family = poisson(link = "log"),
  corstr = "exchangeable"
)

# AR(1)
gee_ar1 <- geeglm(
  seizures ~ trt + log_base + log_age + visit,
  id = id,
  data = epi,
  family = poisson(link = "log"),
  corstr = "ar1"
)

# Unstructured
gee_unstr <- geeglm(
  seizures ~ trt + log_base + log_age + visit,
  id = id,
  data = epi,
  family = poisson(link = "log"),
  corstr = "unstructured"
)

Compare Coefficient Estimates

extract_coefs <- function(fit, name) {
  coef(summary(fit)) |>
    as.data.frame() |>
    tibble::rownames_to_column("term") |>
    mutate(model = name) |>
    select(model, term, Estimate, `Std.err`)
}

bind_rows(
  extract_coefs(gee_ind, "Independence"),
  extract_coefs(gee_exch, "Exchangeable"),
  extract_coefs(gee_ar1, "AR(1)"),
  extract_coefs(gee_unstr, "Unstructured")
) |>
  filter(term == "trtProgabide") |>
  knitr::kable(digits = 4)

Key Observation

Treatment effect estimates are similar across working correlations.

Standard errors differ – efficiency depends on how close working correlation is to truth.

Robust SEs (default in geepack) protect against misspecification.

Estimated Working Correlations

cat("Exchangeable alpha:", gee_exch$geese$alpha, "\n")

Modeling the Association (binary data)

For binary outcomes, FLW Ch. 13 treats the within-subject association as something to MODEL, not just a nuisance:

• Correlations are awkward for binary data (bounded by the means)
• Better target: pairwise odds ratios \(\operatorname{OR}_{jk}\)
• Alternating Logistic Regression (ALR): a second estimating equation for the association parameters \(\alpha\), solved alongside the mean equation for \(\boldsymbol\beta\)

So \(\alpha\) can be a target of inference (how does dependence change over time?), not only a working-correlation nuisance.

Two Types of Standard Errors

Bread (sensitivity): \(\ \mathbf{B} = \sum_i \mathbf{D}_i^\top \mathbf{V}_i^{-1} \mathbf{D}_i\)

Meat (observed variability): \(\ \mathbf{M} = \sum_i \mathbf{D}_i^\top \mathbf{V}_i^{-1} (\mathbf{Y}_i - \boldsymbol\mu_i)(\mathbf{Y}_i - \boldsymbol\mu_i)^\top \mathbf{V}_i^{-1} \mathbf{D}_i\)

When to Use Which?

Comparing SE Types in geepack

# Model-based (naive) SEs require fitting with std.err = "jack"
# then extracting the naive variance from the geese object

# Fit model - geepack stores both variance estimates internally
gee_exch_fit <- geeglm(

  seizures ~ trt + log_base + log_age + visit,
  id = id,
  data = epi,
  family = poisson(link = "log"),
  corstr = "exchangeable"
)

# Extract robust (sandwich) SEs - this is the default in summary()
se_robust <- sqrt(diag(gee_exch_fit$geese$vbeta))

# Extract model-based (naive) SEs from vbeta.naiv
# This is (D'V^{-1}D)^{-1}, assuming working covariance is correct
se_naive <- sqrt(diag(gee_exch_fit$geese$vbeta.naiv))

data.frame(
  term = names(coef(gee_exch_fit)),
  SE_robust = round(se_robust, 4),
  SE_model_based = round(se_naive, 4),
  ratio = round(se_robust / se_naive, 2)
)

Interpretation of SE Comparison

Read the ratio = robust / naive column you just saw:

Ratio > 1 (robust larger): the working correlation overstates the true correlation, so the naive SE is too small. Naive inference is anti-conservative (Type I error inflated). Trust the robust SE.

Ratio < 1 (robust smaller): the working correlation understates the true correlation, so the naive SE is too large. Naive inference is conservative; the robust SE is sharper.

Mechanism: the naive SE believes the working \(\mathbf{V}_i\); when \(\mathbf{V}_i\) overstates correlation it credits each subject with too little independent information, flipping which SE is bigger.

The Small-Sample Problem

Robust SEs rely on an asymptotic (large-\(N\)) approximation:

\[ \hat{\boldsymbol\beta} \ \dot\sim\ N\!\left(\boldsymbol\beta,\ \mathbf{B}^{-1}\mathbf{M}\mathbf{B}^{-1}\right) \]

With few clusters (\(N < 40\)), the sandwich estimator is biased downward (FLW pp. 360-361).

This leads to:

• Underestimated standard errors
• Inflated Type I error rates

Bias-Corrected Sandwich Estimators

\(\mathbf{H}_i\) is the cluster leverage (“hat”) matrix, the cluster analogue of the scalar OLS leverage \(h_{ii}\).

Small-Sample Correction in Practice

# geepack does not implement KC/FG/MD; those live in the geesmv package.
# We illustrate the small-sample idea with a model-agnostic cluster bootstrap.
set.seed(667)
n_boot <- 200
boot_coefs <- matrix(NA, n_boot, 5)

for (b in 1:n_boot) {
  # Resample whole CLUSTERS (subjects) with replacement
  ids <- unique(epi$id)
  boot_ids <- sample(ids, length(ids), replace = TRUE)
  boot_data <- do.call(rbind, lapply(seq_along(boot_ids), function(i) {
    d <- epi[epi$id == boot_ids[i], ]
    d$id <- i   # LOAD-BEARING: GEE needs UNIQUE cluster ids after resampling,
                # or a duplicated subject is merged and the SEs are wrong
    d
  }))

  fit_b <- tryCatch(
    geeglm(
      seizures ~ trt + log_base + log_age + visit,
      id = id,
      data = boot_data,
      family = poisson(link = "log"),
      corstr = "exchangeable"
    ),
    error = function(e) NULL
  )

  if (!is.null(fit_b)) {
    boot_coefs[b, ] <- coef(fit_b)
  }
}

boot_se <- apply(boot_coefs, 2, sd, na.rm = TRUE)
names(boot_se) <- names(coef(gee_exch))

data.frame(
  term = names(boot_se),
  SE_robust = round(se_robust, 4),
  SE_bootstrap = round(boot_se, 4)
)

Degrees of Freedom Adjustments

For small samples, also adjust test statistics:

Use \(t_{N-p}\) or \(F\) reference distributions instead of \(z\) or \(\chi^2\).

Rule of Thumb

With 59 patients here, standard robust SEs are reasonable.

Diagnostic Goals for GEE

Pearson Residuals vs Fitted

Residuals by Visit and Treatment

Cluster-Level Residuals

Identifying Influential Clusters

Subject 49 has extremely high seizure counts:

# Check extreme subjects
epi |>
  group_by(id, trt) |>
  summarize(
    total_seizures = sum(seizures),
    mean_seizures = mean(seizures),
    .groups = "drop"
  ) |>
  arrange(desc(total_seizures)) |>
  head(5)

Sensitivity Analysis: Exclude Outlier

epi_no49 <- epi |> filter(id != 49)

gee_no49 <- geeglm(
  seizures ~ trt + log_base + log_age + visit,
  id = id,
  data = epi_no49,
  family = poisson(link = "log"),
  corstr = "exchangeable"
)

# Compare treatment effects
data.frame(
  model = c("Full data", "Exclude ID 49"),
  trt_coef = c(coef(gee_exch)["trtProgabide"],
               coef(gee_no49)["trtProgabide"]),
  SE = c(summary(gee_exch)$coef["trtProgabide", "Std.err"],
         summary(gee_no49)$coef["trtProgabide", "Std.err"])
) |>
  mutate(
    RR = exp(trt_coef),
    across(c(trt_coef, SE, RR), ~round(.x, 4))
  )

Choosing a Working Correlation: FLW’s method

The fitted \(\hat\alpha\) values and the robust-vs-naive SE comparison (Part III) are the in-chapter diagnostics. QIC (next slide) is a useful add-on, but it is an instructor extension.

QIC: Model Selection for GEE

• QIC -> compare correlation structures (same mean model)
• QICu -> compare mean models (same correlation); QICu is NOT valid for selecting the correlation structure

One actionable rule: smaller is better; use QIC for correlation, QICu for covariates.

QIC: Code

# QIC() returns a named vector; index by NAME (positions are
# geepack-version-dependent). "QIC" ranks correlation structures;
# "QICu" ranks mean models (same correlation).
qic_row <- function(fit) {
  q <- geepack::QIC(fit)
  c(QIC = unname(q["QIC"]), QICu = unname(q["QICu"]))
}

data.frame(
  correlation = c("Independence", "Exchangeable", "AR(1)", "Unstructured"),
  rbind(qic_row(gee_ind), qic_row(gee_exch),
        qic_row(gee_ar1), qic_row(gee_unstr))
) |>
  arrange(QIC)

GEE vs Mixed Models

Important: Standard GEE requires MCAR for valid inference. IPW-GEE extends validity to MAR when the dropout model is correctly specified. (IPW-GEE itself is developed in FLW Ch. 17-18, not Ch. 13; we preview it in Part VII.)

Coefficient Interpretation: link matters

For an identity link (continuous data): GEE and mixed-model \(\boldsymbol\beta\) have the same interpretation.

For a nonlinear link:

• GEE \(\boldsymbol\beta\): population-average (marginal) effect
• Mixed \(\boldsymbol\beta\mid \mathbf{b}\): subject-specific (conditional) effect

Attenuation: the logistic picture

For a binary outcome with a subject random intercept \(b_i \sim N(0, \sigma_b^2)\), averaging the individual logistic curves over \(b_i\) flattens the population-average curve:

Mixed Model for Comparison

# GLMM with random intercept
glmm_fit <- glmer(
  seizures ~ trt + log_base + log_age + visit + (1 | id),
  data = epi,
  family = poisson(link = "log")
)

# Compare coefficients
data.frame(
  term = c("(Intercept)", "trtProgabide", "log_base", "log_age", "visit"),
  GEE = round(coef(gee_exch), 4),
  GLMM = round(fixef(glmm_fit), 4)
)

Marginal vs Conditional Effects

# Treatment effect comparison
trt_gee <- coef(gee_exch)["trtProgabide"]
trt_glmm <- fixef(glmm_fit)["trtProgabide"]

data.frame(
  Model = c("GEE (marginal)", "GLMM (conditional)"),
  log_RR = c(trt_gee, trt_glmm),
  RR = exp(c(trt_gee, trt_glmm)),
  Interpretation = c("Population-average", "Subject-specific")
) |>
  mutate(across(c(log_RR, RR), ~round(.x, 4)))

When to Choose GEE vs GLMM

Note on missing data: Standard GEE is only valid under MCAR. If dropout depends on observed data (MAR), use IPW-GEE or likelihood-based methods (GLMM).

The Missing Data Problem in GEE

Standard GEE assumes MCAR (Missing Completely at Random):

• Probability of dropout does not depend on outcomes or covariates
• Example: Random equipment failure, administrative censoring

Reality: Dropout often depends on observed data (MAR):

• Patients with more seizures may drop out
• Treatment side effects may cause attrition
• This violates MCAR and biases standard GEE

IPW-GEE Example: Simulated Dropout

# Create data with MAR dropout
set.seed(667)
epi_mar <- epi |>
  group_by(id) |>
  mutate(
    # Dropout probability depends on previous seizures
    prev_seizures = lag(seizures, default = 0),
    # Higher seizures -> more likely to drop out
    dropout_prob = plogis(-2 + 0.1 * prev_seizures),
    dropout = rbinom(n(), 1, dropout_prob),
    # Cumulative dropout
    ever_dropped = cummax(dropout),
    observed = 1 - ever_dropped
  ) |>
  ungroup()

# Keep only observed data
epi_obs <- epi_mar |> filter(observed == 1)

cat("Original observations:", nrow(epi), "\n")

Fit Dropout Model

# Model probability of being observed given history
dropout_fit <- glm(
  observed ~ trt + log_base + log_age + visit + prev_seizures,
  data = epi_mar,
  family = binomial()
)

# Get predicted probabilities for observed subjects
epi_obs$pi_hat <- predict(dropout_fit, newdata = epi_obs, type = "response")
epi_obs$ipw <- 1 / epi_obs$pi_hat

# Check weight distribution
summary(epi_obs$ipw)

Compare Standard GEE vs IPW-GEE

# Standard GEE (ignoring dropout mechanism)
gee_standard <- geeglm(
  seizures ~ trt + log_base + log_age + visit,
  id = id,
  data = epi_obs,
  family = poisson(link = "log"),
  corstr = "exchangeable"
)

# IPW-GEE (accounting for MAR dropout)
gee_ipw <- geeglm(
  seizures ~ trt + log_base + log_age + visit,
  id = id,
  data = epi_obs,
  family = poisson(link = "log"),
  corstr = "exchangeable",
  weights = ipw
)

# Compare treatment effects (save for inline interpretation next slide)
ipw_compare <- data.frame(
  Method = c("Full data (no dropout)", "Standard GEE (MAR ignored)", "IPW-GEE (MAR corrected)"),
  trt_coef = round(c(coef(gee_exch)["trtProgabide"],
                     coef(gee_standard)["trtProgabide"],
                     coef(gee_ipw)["trtProgabide"]), 4),
  SE = round(c(summary(gee_exch)$coef["trtProgabide", "Std.err"],
               summary(gee_standard)$coef["trtProgabide", "Std.err"],
               summary(gee_ipw)$coef["trtProgabide", "Std.err"]), 4)
)
ipw_compare

Interpreting the IPW-GEE Comparison

With the full data, the treatment log-rate-ratio is -0.042. After MAR dropout:

• Standard GEE (ignoring dropout) gives -0.08, a distance of 0.038 from the full-data value: dropout has biased it.
• IPW-GEE gives -0.19, a distance of 0.148 from full-data, not closer to the truth.

Takeaway: ignoring informative dropout pulls the estimate away from the full-data target; weighting by \(1/\hat\pi_{ij}\) moves it back (here, not all the way to) the full-data value.

IPW-GEE: Practical Considerations

Rule of thumb: If > 10% dropout and mechanism may be MAR, consider IPW-GEE.

Summary: GEE and Missing Data

Key point: Standard GEE is only valid under MCAR. For MAR, use IPW-GEE with a correctly specified dropout model.

Full GEE Analysis: Epilepsy Data

# Final model with exchangeable correlation
summary(gee_exch)

Exponentiated Coefficients (Rate Ratios)

coefs <- coef(summary(gee_exch))

# Compute Wald-based 95% CIs manually
estimates <- coefs[, "Estimate"]
se <- coefs[, "Std.err"]
ci_lower <- estimates - 1.96 * se
ci_upper <- estimates + 1.96 * se

results <- data.frame(
  term = rownames(coefs),
  Estimate = estimates,
  SE = se,
  RR = exp(estimates),
  RR_lower = exp(ci_lower),
  RR_upper = exp(ci_upper),
  p_value = coefs[, "Pr(>|W|)"]
)

results |>
  mutate(across(c(Estimate, SE, RR, RR_lower, RR_upper), ~round(.x, 3))) |>
  mutate(p_value = format.pval(p_value, digits = 3))

Interpreting Treatment Effect

Treatment coefficient: \(\hat\beta =\) -0.042

Rate Ratio: \(\exp(\hat\beta) =\) 0.959

95% CI: (0.66, 1.38)

Interpretation: On average, progabide patients have about 4.1% lower seizure rates than placebo, but this effect is not statistically significant (p = 0.82).

Plain-Language Interpretation: GEE Output

• trtProgabide: rate ratio 0.96, about 4.1% lower (not significant)
• log_base: strongest predictor; higher baseline rate -> higher rate throughout
• visit: small decline per visit

Effect of Baseline Seizure Rate

log_base coefficient: \(\hat\beta =\) 1.223

For each doubling of baseline rate: \[\exp(\hat\beta \times \log 2) = \exp(0.848) = 2.33\]

Patients with higher baseline rates continue to have higher rates on treatment.

Predicted Rates

# Predictions for average patient
newdata <- expand.grid(
  trt = levels(epi$trt),
  log_base = log(median(epi_raw$base) / 4),
  log_age = log(median(epi_raw$age)),
  visit = 1:4
)

# Use the model matrix to show the linear-predictor mechanics
# (predict(gee_exch, newdata, type = "response") would also work)
X <- model.matrix(~ trt + log_base + log_age + visit, data = newdata)
newdata$pred_rate <- exp(X %*% coef(gee_exch))

ggplot(newdata, aes(visit, pred_rate, color = trt)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 3) +
  labs(
    title = "Predicted Seizure Rates (Median Patient)",
    x = "Visit",
    y = "Predicted Seizures per 2 Weeks",
    color = "Treatment"
  ) +
  theme_minimal()

GEE Extensions: Overdispersion

For count data, we may have overdispersion (\(\phi > 1\)):

# Estimate dispersion
pearson_stat <- sum(residuals(gee_exch, type = "pearson")^2)
df_resid <- nrow(epi) - length(coef(gee_exch))
phi_hat <- pearson_stat / df_resid

cat("Estimated dispersion:", round(phi_hat, 2), "\n")

Summary

Key Takeaways

Common Mistakes in GEE

GEE Recipe Card: A Step-by-Step Workflow

fit <- geeglm(y ~ x1 + x2 + time, id = subject_id,
              data = mydata, family = ..., corstr = "exchangeable")

QIC(fit_exch); QIC(fit_ar1); QIC(fit_unstr)

The Same GEE Fit in SAS

For SAS users, the exchangeable Poisson GEE is PROC GENMOD with a REPEATED statement (shown for reference, not run):

proc genmod data=epi;
   class id trt;
   model seizures = trt log_base log_age visit / dist=poisson link=log;
   repeated subject=id / type=exch corrw modelse;
run;

• class id trt; declares the cluster id and the categorical treatment
• dist=poisson link=log; is the marginal Poisson-log mean model
• repeated subject=id / type=exch requests the exchangeable working correlation (use type=ar(1) or type=un for the others)
• corrw prints the estimated working correlation; modelse adds the model-based SEs alongside the default empirical (sandwich) SEs

For Next Time

• Chapter 14: Transition models
• Modeling dependence on previous outcomes
• Comparison: marginal vs transition interpretations

Reading: FLW Chapter 14

Further Reading

• FLW Ch. 13 (this chapter): marginal models, GEE, working correlation, sandwich SEs, the Section 13.4 case studies (Muscatine, leprosy, arthritis)
• FLW Ch. 17-18: inverse probability weighting and IPW-GEE for MAR dropout (the proper development of Part VII)
• Liang & Zeger (1986): the original GEE paper
• Pan (2001): QIC for GEE model selection (the Part V extension)

References

• Fitzmaurice, Laird, Ware (2011). Applied Longitudinal Analysis. Ch. 13 (GEE, pp. 354-393); IPW-GEE developed in Ch. 17-18.
• Liang, Zeger (1986). Longitudinal data analysis using generalized linear models. Biometrika.
• Pan (2001). Akaike’s information criterion in generalized estimating equations. Biometrics. (QIC; beyond Ch. 13)
• Kauermann, Carroll (2001); Fay, Graubard (2001); Mancl, DeRouen (2001). Small-sample sandwich corrections. JASA / Biometrics. (beyond Ch. 13)
• Robins, Rotnitzky, Zhao (1995). Semiparametric regression with missing data. JASA. (IPW)

Appendix A: GEE Derivation and Fisher Scoring

The generalized estimating equation:

\[ U(\boldsymbol\beta) = \sum_{i=1}^{N} \mathbf{D}_i^\top \mathbf{V}_i^{-1}(\mathbf{Y}_i - \boldsymbol\mu_i) = 0 \]

is solved by an iterative update (Fisher scoring, using the expected information):

\[ \hat{\boldsymbol\beta}^{(t+1)} = \hat{\boldsymbol\beta}^{(t)} + \left(\sum_i \mathbf{D}_i^\top \mathbf{V}_i^{-1} \mathbf{D}_i\right)^{-1} \sum_i \mathbf{D}_i^\top \mathbf{V}_i^{-1}(\mathbf{Y}_i - \boldsymbol\mu_i^{(t)}) \]

Fisher scoring vs Newton-Raphson: Fisher scoring uses the expected information (model-based, more stable); Newton-Raphson uses the observed information (data-based). They coincide for GLMs with a canonical link in the IRLS sense, but that exact equivalence does not carry over to GEE with a non-identity working correlation. The update above is iteratively reweighted least squares with the current weights.

Appendix B: Sandwich Variance

The robust variance estimator:

\[ \widehat{\operatorname{Var}}(\hat{\boldsymbol\beta}) = \mathbf{B}^{-1} \mathbf{M} \mathbf{B}^{-1} \]

where (bread):

\[ \mathbf{B} = \sum_i \mathbf{D}_i^\top \mathbf{V}_i^{-1} \mathbf{D}_i \]

and (meat):

\[ \mathbf{M} = \sum_i \mathbf{D}_i^\top \mathbf{V}_i^{-1} (\mathbf{Y}_i - \hat{\boldsymbol\mu}_i)(\mathbf{Y}_i - \hat{\boldsymbol\mu}_i)^\top \mathbf{V}_i^{-1} \mathbf{D}_i \]

Appendix C: QIC Formula (Pan 2001; beyond Ch. 13)

\[ \operatorname{QIC}(R) = -2Q(\hat{\boldsymbol\beta}_R; I) + 2\operatorname{trace}(\hat{\boldsymbol\Omega}_I^{-1} \hat{\mathbf{V}}_R) \]

For comparing mean models with the same correlation:

\[ \operatorname{QICu} = -2Q(\hat{\boldsymbol\beta}; I) + 2p \]

• \(Q(\hat{\boldsymbol\beta}; I)\) is the quasi-likelihood evaluated under the independence working model
• \(\hat{\boldsymbol\Omega}_I\) is the independence-model (naive) variance
• Smaller QIC indicates a better model