GEE Models For Recurrent Event Counts

Why Use a GEE for Repeated Event Counts?

A recurrent-event data set can be organized as a count for each subject in each follow-up interval. For example, a study may record how many symptom episodes each participant experiences in each month. Counts from the same participant are generally dependent: a participant with many episodes in one month may also tend to have many in another. Treating all participant-months as independent can therefore give unreliable standard errors.

A generalized estimating equation (GEE) represents the mean count in each interval and uses a working covariance matrix to account for within-subject dependence. Its central target is a marginal mean: the expected count in an interval for subjects with specified covariates, averaged over the population variation not included in those covariates. The working correlation helps estimate that mean model efficiently and calculate appropriate uncertainty; it is not itself a model for the sequence of event times.

Mean Model and Estimating Equation

A common specification for counts is a log-link model with an exposure offset:

Here \(\beta\) is the vector of regression coefficients. The offset \(\log(E_{ij})\) has coefficient fixed at one; it adjusts the expected count for the amount of time or other exposure observed. If every interval has the same exposure, the offset may be omitted or absorbed into the intercept. With an offset, \(\exp(X_{ij}^{\top}\beta)\) is the modeled mean count per unit exposure. Without one, the model describes the expected count over the interval represented by the data.

For each subject, write \(\mu_i=(\mu_{i1},\ldots,\mu_{im_i})^\top\) and define \(D_i=\partial\mu_i/\partial\beta^\top\), the matrix of derivatives of the mean vector with respect to the coefficients. A working covariance matrix has the form \(V_i=A_i^{1/2}R_i(\alpha)A_i^{1/2}\), where \(A_i\) is diagonal and contains the marginal variance-function values, and \(R_i(\alpha)\) is a chosen working correlation matrix. Under a Poisson variance function, \(A_i=\operatorname{diag}(\mu_{i1},\ldots,\mu_{im_i})\). Other variance functions can be used when they better represent marginal dispersion.

This equation can be understood as a weighted residual condition. The residual vector \(Y_i-\mu_i\) is scaled and combined across intervals using \(V_i^{-1}\); \(D_i^\top\) then projects those weighted residuals onto the directions in which changing \(\beta\) changes the mean. The estimate is reached when the summed projected residuals are zero. If \(V_i\) is the true covariance, this weighting can improve precision. If it is only an approximation, the mean equation can still estimate \(\beta\) consistently under the assumptions below, while robust variance estimation accounts for the actual within-subject variability.

Assumptions and Inference

The usual large-sample justification requires independent subjects (clusters), a correctly specified marginal mean model, and enough independent subjects for the cluster-level estimating equation to be well approximated by its large-sample distribution. Within a subject, counts may be dependent. The chosen \(V_i\) must be invertible for the observed data and should be a reasonable working covariance; it need not equal the true covariance for consistency of the mean estimator. Covariates and exposure must be measured appropriately, and the observed intervals must represent the target follow-up. Informative loss of intervals or follow-up requires additional assumptions or methods; a working correlation alone does not correct it.

The model may use a Poisson variance function even when the counts are not Poisson distributed, provided the mean is correctly specified and robust inference is used. In that case the Poisson variance is a working choice, not a claim that the full count distribution is Poisson. A common large-sample robust covariance estimator has the sandwich form

where \(\hat U_i=D_i^\top\hat V_i^{-1}(Y_i-\hat\mu_i)\) is the fitted contribution from subject \(i\). This variance estimate allows the actual covariance within subjects to differ from the working covariance, provided the mean model and independent-cluster large-sample conditions hold. It does not protect against a misspecified mean model, dependent subjects, or a very small number of clusters. Robust variance in marginal rate models is covered in the earlier material on robust variance in marginal rate models.

Working Correlation: What It Does and Does Not Say

The working correlation describes the dependence pattern used in \(R_i(\alpha)\), conditional on the fitted marginal means and variance function. Common choices include:

• Independence: \(R_i\) is the identity matrix. It assumes no within-subject correlation for weighting, but the robust variance can still account for actual clustering.
• Exchangeable: all distinct intervals within a subject have the same working correlation \(\alpha\).
• AR(1): the working correlation between intervals decreases with their separation, typically as \(\alpha^{|j-k|}\).
• Unstructured: each pair of intervals has its own working correlation, requiring substantially more information and a positive-definite estimate.

These are alternative covariance approximations, not alternative marginal mean estimands. A working correlation is not a claim that one event causes another, nor does it encode a model for event order. A complicated structure is not automatically better: it can be unstable or non-positive-definite, especially with few subjects or many intervals. Independence is often a useful starting point, while a plausible simpler structure may improve efficiency.

Worked Rate Model With a Treatment Indicator

Diagnostics, Comparisons, and Limitations

Check whether the fitted means are plausible across covariate patterns and intervals, whether residuals show systematic trends, and whether estimated working correlations are stable and admissible. Large residuals or patterns over time may indicate an inadequate mean model, an unsuitable variance function, or influential subjects. A correlation estimate at a boundary, failed matrix inversion, or substantial changes across plausible working structures are warnings to investigate rather than reasons to report whichever fit is most favorable. With few independent subjects, sandwich standard errors can be downward biased; robust inference is not automatically reliable merely because it is labelled robust.

GEE is useful when the scientific target is a population-average mean count and one wants inference robust to a misspecified within-subject covariance. A fully specified count likelihood can be more efficient if its distribution and dependence assumptions are credible, but incorrect distributional assumptions can affect both estimation and inference. A subject-specific random-effects model targets a different conditional interpretation and is not interchangeable with a marginal GEE. For recurrent-event questions framed around event order or event-time risk sets, PWP and WLW models provide different structures and estimands; the choice should follow the scientific question, not just data format.

A frequent mistake is to read an exchangeable working correlation as evidence that each subject has a fixed common tendency causing their recurrences, or to believe that robust standard errors repair every modeling problem. The correction is to distinguish the components: the mean model defines the target, the working covariance supplies weights, and the sandwich variance addresses covariance misspecification under independent clusters and a correct mean model. It cannot repair a wrong mean, informative missingness, or dependence between subjects.