Building Cluster-Level Estimating Functions

From Marginal Means to Cluster Contributions

In a generalized estimating equation, observations from the same cluster are handled together. The cluster’s response vector supplies residuals, the marginal mean model describes the expected response vector, and a working covariance matrix weights those residuals. The contribution to the estimating function is built from these three ingredients, then added to the contributions from the other clusters.

The construction is useful because it makes the data structure visible in the algebra. It also provides a practical check: every matrix multiplication must have compatible dimensions, and the final contribution must have one entry for each regression coefficient.

Definition: For cluster $i$, let $m_i$ be its number of observations, $Y_i=(Y_{i1},\ldots,Y_{im_i})^T$ its response vector, and $y_i$ the observed value of that vector. Let $\mu_i(\beta)=E[Y_i\mid X_i]$ be the marginal mean vector under a model with parameter $\beta\in\mathbb{R}^p$. The mean derivative matrix is $$ D_i(\beta)=\frac{\partial \mu_i(\beta)}{\partial \beta^T}, $$ an $m_i\times p$ matrix. Let $V_i(\beta,\alpha)$ be an $m_i\times m_i$ symmetric positive-definite working covariance matrix, where $\alpha$ denotes any working covariance parameters. The cluster-level estimating function is $$ U_i(\beta)=D_i(\beta)^T V_i(\beta,\alpha)^{-1}\{y_i-\mu_i(\beta)\}, $$ and the overall estimating function is $U(\beta)=\sum_{i=1}^K U_i(\beta)$ for $K$ clusters.

The $j$th component of $\mu_i$ is the marginal mean for observation $j$ in cluster $i$. Row $j$ of $D_i$ is the derivative of that mean with respect to all $p$ components of $\beta$. The residual vector $y_i-\mu_i$ has length $m_i$. Multiplication by $V_i^{-1}$ weights and combines its entries; multiplication by $D_i^T$ maps the result to a vector of length $p$.

The covariance matrix is called a working covariance because it is part of the estimating-function construction. It must be invertible for this expression to be defined. It can depend on $\beta$ and on working parameters such as $\alpha$; in a calculation, use its value at the parameter values being considered. The derivative in $D_i$ is the derivative of the mean vector with respect to $\beta$, not a derivative of the working covariance.

Reading the Formula One Coefficient at a Time

The matrix expression can be expanded to show exactly how each observation contributes to each coefficient’s equation. This form is often the quickest way to check a calculation or identify which mean derivatives are zero.

Result: Suppose $\mu_i(\beta)$ is differentiable and $V_i$ is symmetric positive definite. For coefficient $r\in\{1,\ldots,p\}$, $$ U_{ir}(\beta) = \sum_{j=1}^{m_i} \frac{\partial \mu_{ij}(\beta)}{\partial \beta_r} \left[V_i(\beta,\alpha)^{-1}\{y_i-\mu_i(\beta)\}\right]_j. $$

Proof. Set $q_i=V_i^{-1}(y_i-\mu_i)$, an $m_i$-vector. By the definition of the derivative matrix, its entry in row $j$, column $r$ is $D_{i,jr}=\partial\mu_{ij}/\partial\beta_r$. The $r$th entry of $D_i^Tq_i$ is therefore the row-$r$, column-$1$ entry of that product: $$ (D_i^Tq_i)_r=\sum_{j=1}^{m_i}(D_i^T)_{rj}(q_i)_j =\sum_{j=1}^{m_i}\frac{\partial\mu_{ij}}{\partial\beta_r}(q_i)_j. $$ Substituting the definition of $q_i$ gives the stated expression. $\square$

This expansion distinguishes two roles. The $j$th component of the covariance-weighted residual is multiplied by how sensitive the $j$th mean is to coefficient $r$. An observation with zero derivative for that coefficient makes no direct contribution to its component of $U_i$, although it can still affect the weighted residual through the inverse covariance matrix.

Worked Example: A Binary Cluster With an Intercept

Suppose a cluster has two binary responses and an intercept-only logistic mean model. At $\beta=0$, each marginal mean is $\mu_{ij}=1/2$, and the derivative of each mean with respect to $\beta$ is $\mu_{ij}(1-\mu_{ij})=1/4$. Let $y_i=(1,1)^T$, and use the working covariance $$ V_i=\frac14 \begin{pmatrix} 1&0.2\\ 0.2&1 \end{pmatrix}. $$ Calculate the cluster contribution.

The residual vector is $r_i=y_i-\mu_i=(1/2,1/2)^T$, and $D_i=(1/4,1/4)^T$. Since the determinant of the matrix with diagonal entries 1 and off-diagonal entries 0.2 is $1-0.2^2=0.96$, $$ V_i^{-1}=\frac{4}{0.96} \begin{pmatrix} 1&-0.2\\ -0.2&1 \end{pmatrix} =\frac{25}{6} \begin{pmatrix} 1&-0.2\\ -0.2&1 \end{pmatrix}. $$

Applying this inverse to the residual gives $$ V_i^{-1}r_i =\frac{25}{6} \begin{pmatrix} 0.5-0.1\\ -0.1+0.5 \end{pmatrix} = \begin{pmatrix}5/3\\5/3\end{pmatrix}. $$ Thus the scalar cluster contribution is $$ U_i(0)=D_i^TV_i^{-1}r_i =\frac14(5/3+5/3)=\frac56. $$ The calculation illustrates the prescribed order: first weight the residual vector with $V_i^{-1}$, then multiply by the transposed mean derivative.

Working Through Different Mean Models

The same assembly rule applies when the response distribution and link function change. What changes is the mean vector and its derivative matrix; the dimensions and order of the calculation remain the same.

Worked Example: A Continuous Response With a Linear Mean

For one two-observation cluster, let $X_i=\begin{pmatrix}1&0\\1&1\end{pmatrix}$, $\mu_i(\beta)=X_i\beta$, $\beta=(2,1)^T$, and $y_i=(3,1)^T$. Use $$ V_i=\begin{pmatrix}4&1\\1&4\end{pmatrix}. $$ Find $U_i(\beta)$.

The fitted mean vector and residual vector are $$ \mu_i=X_i\beta=\begin{pmatrix}2\\3\end{pmatrix}, \qquad r_i=y_i-\mu_i=\begin{pmatrix}1\\-2\end{pmatrix}. $$ Because the mean is linear in $\beta$, its derivative matrix is $D_i=X_i$. The determinant of $V_i$ is $16-1=15$, so $$ V_i^{-1}=\frac1{15}\begin{pmatrix}4&-1\\-1&4\end{pmatrix}, \qquad V_i^{-1}r_i=\frac1{15}\begin{pmatrix}6\\-9\end{pmatrix} =\begin{pmatrix}2/5\\-3/5\end{pmatrix}. $$ Therefore, $$ U_i(\beta)=X_i^TV_i^{-1}r_i =\begin{pmatrix}1&1\\0&1\end{pmatrix} \begin{pmatrix}2/5\\-3/5\end{pmatrix} =\begin{pmatrix}-1/5\\-3/5\end{pmatrix}. $$ The contribution has two entries, one for each coefficient. The first combines both weighted residuals because the intercept changes both fitted means; the second uses only the second row because the first observation’s predictor value for that coefficient is zero.

Worked Example: Poisson Counts With a Log Mean

Consider two count responses in one cluster with log means $\log\mu_{ij}=\beta+o_{ij}$, where the offsets are $o_{i1}=0$ and $o_{i2}=\log 4$. At $\beta=0$, the mean vector is $(1,4)^T$, and $D_i=(1,4)^T$, since $\partial\mu_{ij}/\partial\beta=\mu_{ij}$. Let $y_i=(2,1)^T$, and use the working covariance $$ V_i=\begin{pmatrix}1&1\\1&4\end{pmatrix}. $$

This covariance is positive definite because its leading diagonal entry is 1 and its determinant is $4-1=3>0$. The residual vector is $r_i=(1,-3)^T$, and $$ V_i^{-1}=\frac13\begin{pmatrix}4&-1\\-1&1\end{pmatrix}, \qquad V_i^{-1}r_i =\frac13\begin{pmatrix}7\\-4\end{pmatrix}. $$ The scalar contribution is consequently $$ U_i(0)=D_i^TV_i^{-1}r_i =(1,4)\begin{pmatrix}7/3\\-4/3\end{pmatrix} =\frac73-\frac{16}{3}=-3. $$ The negative value reflects the direction of the covariance-weighted residual relative to the mean’s change with $\beta$; it is not obtained by simply summing the raw residuals, which here sum to $-2$.

How Contributions Combine

Each cluster produces a $p$-vector, regardless of its size. The overall estimating function is the ordinary vector sum of those contributions. It is not necessary for clusters to have equal sizes: $D_i$, $V_i$, and $y_i-\mu_i$ are formed at the size $m_i$ of cluster $i$, while the resulting $U_i$ always has length $p$.

Result: Suppose a cluster’s observations can be divided into blocks, and, in the same ordering, its mean derivative matrix and working covariance have forms $$ D_i=\begin{pmatrix}D_{i1}\\ \vdots\\D_{iB}\end{pmatrix}, \qquad V_i=\operatorname{blockdiag}(V_{i1},\ldots,V_{iB}), $$ where each $V_{ib}$ is invertible. If $r_{ib}$ is the residual vector for block $b$, then $$ D_i^TV_i^{-1}(y_i-\mu_i) =\sum_{b=1}^{B}D_{ib}^TV_{ib}^{-1}r_{ib}. $$

Proof. The inverse of a block-diagonal matrix with invertible diagonal blocks is block diagonal with the inverse of each block. This follows by multiplying the two block-diagonal matrices: each diagonal block of their product is $V_{ib}V_{ib}^{-1}=I$, and every off-diagonal block is zero. Hence $$ V_i^{-1}(y_i-\mu_i) = \begin{pmatrix} V_{i1}^{-1}r_{i1}\\ \vdots\\ V_{iB}^{-1}r_{iB} \end{pmatrix}. $$ Multiplying on the left by $D_i^T=(D_{i1}^T,\ldots,D_{iB}^T)$ gives $$ D_i^TV_i^{-1}(y_i-\mu_i) =\sum_{b=1}^{B}D_{ib}^TV_{ib}^{-1}r_{ib}, $$ as claimed. $\square$

This result is an algebraic statement about the specified block structure. When the working covariance has nonzero entries linking blocks, its inverse may also link them, and the contribution generally cannot be calculated as a sum of separate block contributions using only the corresponding block covariance matrices.

Common Mistakes and a Practical Check

Reversing the matrix order. The expression $D_i^TV_i^{-1}r_i$ ends in a $p$-vector. The order $V_i^{-1}D_i^Tr_i$ is generally not even dimensionally defined.
Using the covariance instead of its inverse. The residual weighting step is $V_i^{-1}r_i$, not $V_ir_i$. Check the dimensions and the specified formula before multiplying.
Differentiating the working covariance in this construction. $D_i$ contains derivatives of the mean vector with respect to $\beta$. The cluster contribution is defined using that matrix and the current inverse working covariance; it does not add a derivative of $V_i$.
Adding raw residuals across observations or clusters. Residuals are first weighted within their cluster and mapped through $D_i^T$. Only the resulting $p$-vectors are summed across clusters.

A reliable calculation proceeds in four checks: write down $\mu_i$ and $D_i$; form the residual $y_i-\mu_i$; compute $V_i^{-1}(y_i-\mu_i)$; and multiply by $D_i^T$. Finally, verify that every $U_i$ has $p$ entries before summing over $i$. This keeps the cluster structure explicit without requiring equal cluster sizes or a common response type.

Check Your Understanding

Use the definitions and calculations above to answer the following.

If cluster $i$ has $m_i$ observations and the model has $p$ coefficients, what are the dimensions of $D_i$, $V_i^{-1}$, and $U_i$?
In the coordinate expansion of $U_{ir}$, what does the derivative $\partial\mu_{ij}/\partial\beta_r$ measure?
Why can an observation affect a coefficient’s contribution even when its own mean derivative for that coefficient is zero?
For the Poisson example, verify that the stated working covariance is positive definite using its determinant and leading diagonal entry.
Under what covariance structure can a cluster contribution be written as a sum of separate block contributions?