Poisson GEE for Clustered Counts

Why a Poisson GEE for Clustered Counts?

Counts such as clinic visits, infections, or emergency-department arrivals are often recorded repeatedly for the same person, or for people grouped within a clinic. A Poisson mean model gives a useful starting point when the response is a count and its mean should remain positive. But repeated observations within a cluster may be associated, so treating every count as independent can use the wrong residual weighting.

A Poisson generalized estimating equation (GEE) combines a marginal mean model with a working covariance matrix for each cluster. The mean model determines the population-average relationship of interest; the working covariance determines how residuals within a cluster contribute to estimating its coefficients. The working covariance need not be the true covariance for the estimating equation to have mean zero when the marginal mean is correctly specified, though a different working covariance can produce a different finite-sample estimate.

Mean Model and Cluster Estimating Equation

Let \(Y_{ij}\) be the count for observation \(j\) in cluster \(i\), where \(i=1,\ldots,K\) and \(j=1,\ldots,m_i\). Clusters are independent of one another; observations within a cluster may be dependent. Let \(x_{ij}\) be a \(p\)-component covariate vector, including an intercept if appropriate. The log-link marginal mean model is

For a Poisson marginal variance specification, \(\operatorname{Var}(Y_{ij}\mid x_{ij})=\mu_{ij}\). This specifies the one-observation variance used in the GEE; it does not assert that observations within a cluster are independent. Write \(Y_i=(Y_{i1},\ldots,Y_{im_i})'\), \(\mu_i=(\mu_{i1},\ldots,\mu_{im_i})'\), and let \(X_i\) be the \(m_i\) by \(p\) matrix whose rows are \(x_{ij}'\). Define \(A_i=\operatorname{diag}(\mu_{i1},\ldots,\mu_{im_i})\).

The matrix \(R_i(\alpha)\) is a working correlation matrix: it is used to weight residuals, and its parameters may be estimated from the data. For a two-observation exchangeable working correlation with correlation \(\rho\), \(R_i=\begin{pmatrix}1&\rho\\\rho&1\end{pmatrix}\), where \(|\rho|<1\) ensures positive definiteness. Setting \(R_i=I\) gives the independence working correlation.

Deriving the Poisson Cluster Score

The log link gives \(\partial\mu_{ij}/\partial\beta=\mu_{ij}x_{ij}\), so \(D_i=A_iX_i\). Substituting this derivative and the working covariance into the GEE yields a convenient expression for the contribution \(U_i(\beta)\) from cluster \(i\).

To verify the expression, use \(V_i^{-1}=\phi^{-1}A_i^{-1/2}R_i^{-1}A_i^{-1/2}\). Then

This derivation shows exactly where the working correlation enters: it acts on standardized residuals \(A_i^{-1/2}(Y_i-\mu_i)\), and the resulting vector is scaled by \(A_i^{1/2}\) and projected onto the covariates. Under independence, the expression simplifies to \(U_i=X_i'A_i^{-1}(Y_i-\mu_i)\) only if read incorrectly; the correct simplification is \(U_i=X_i'(Y_i-\mu_i)/\phi\), since \(A_i^{1/2}I A_i^{-1/2}=I\). Thus each observation contributes its covariate vector times its raw residual, up to the common dispersion factor.

What Exchangeable Correlation Does to an Intercept Score

A particularly transparent case is an intercept-only model, so every observation has the same mean \(\mu=\exp(\beta)\). Suppose a cluster has size \(m\) and exchangeable working correlation \(R=(1-\rho)I+\rho\mathbf{1}\mathbf{1}'\), with \(-1/(m-1)<\rho<1\) when \(m>1\). Here \(\mathbf{1}\) is the \(m\)-vector of ones. The cluster contribution can be reduced to a scalar multiple of its total residual.

To derive it, \(A_i=\mu I\), so the factors \(A_i^{1/2}\) and \(A_i^{-1/2}\) cancel, leaving \(U_i=\phi^{-1}\mathbf{1}'R_i^{-1}(Y_i-\mu\mathbf{1})\). Since \(R_i\mathbf{1}=\{1+(m_i-1)\rho\}\mathbf{1}\), it follows that \(R_i^{-1}\mathbf{1}=\{1+(m_i-1)\rho\}^{-1}\mathbf{1}\). Symmetry of \(R_i^{-1}\) then gives \(\mathbf{1}'R_i^{-1}=\{1+(m_i-1)\rho\}^{-1}\mathbf{1}'\), proving the formula. The denominator is positive under the stated exchangeable validity bounds. □

Thus a cluster total residual is divided by a size-dependent factor. For \(\rho>0\), larger clusters receive less weight relative to their number of observations than under independence. If all clusters have the same size, this factor is common to every cluster and cancels from the estimating equation; with unequal sizes, it can change the estimate.

A Direct Check of the Inverse-Correlation Weight

Even in a simple cluster, the inverse correlation matrix can change a score in ways that are easy to miscalculate. For two observations with equal fitted means, the standardized residuals are the raw residuals divided by the common square root of the mean. The following calculation checks the inverse and the resulting intercept contribution explicitly.

Interpretation, Assumptions, and a Common Pitfall

In a log-link marginal model, \(\beta_r\) describes a population-average log-mean change per unit increase in covariate \(r\), holding the other covariates fixed. Exponentiating a coefficient gives the corresponding multiplicative mean ratio. The GEE uses the marginal means and the cluster residual weights to estimate these coefficients. It does not require a fully specified joint distribution for all counts within a cluster.

For the usual unbiasedness argument, assume independent clusters, finite relevant expectations, a correctly specified marginal mean conditional on the covariates, and a working covariance/inverse that is well-defined and determined by covariates or by nuisance quantities treated as fixed in the estimating equation. Then \(E[Y_i-\mu_i\mid X_i]=0\), and consequently \(E[U_i(\beta_0)\mid X_i]=0\). The choice of working correlation can affect efficiency and the realized estimate; it does not repair a misspecified mean model. Robust sandwich inference addresses covariance misspecification under its regularity conditions, but does not make a wrong mean model correct.

• Pitfall: treating the Poisson working variance as proof of independence. The marginal Poisson variance specifies the diagonal entries of the working covariance. The off-diagonal entries are supplied by \(R_i(\alpha)\); a nonidentity working correlation explicitly allows within-cluster dependence.
• Pitfall: multiplying by \(R_i^{-1}\) but forgetting the variance scaling. For unequal fitted means, the score weights standardized residuals, not raw residuals alone. The factors \(A_i^{1/2}\) and \(A_i^{-1/2}\) must remain in the expression.
• Pitfall: calling the GEE a Poisson likelihood score. Marginal Poisson distributions do not specify a joint likelihood for dependent observations. The GEE is an estimating equation built from marginal means and a working covariance.
• Pitfall: claiming correlation always changes the coefficient estimate. In the intercept-only exchangeable case with equal cluster sizes, the common cluster weight cancels from the equation. Unequal sizes, covariate patterns, or means can make the weighting consequential, but a change is not automatic.