GEE Data Structures and Marginal Inputs

Why GEE Inputs Are Organized by Cluster

In a generalized estimating equations analysis, each cluster contributes a vector of responses rather than a collection of unrelated observations. A cluster might be a patient measured at several visits, a household containing several people, or a clinic enrolling multiple patients. The response measurements within a cluster may be dependent, while distinct clusters are treated as independent sampling units. Organizing the data this way makes the mean model and the assumed working covariance structure explicit.

This tutorial defines those data objects. It does not yet combine them into an estimating function. The key distinction is between the marginal mean, which specifies the expected response at each observation, and the working covariance matrix, which represents a chosen covariance structure for the responses in one cluster.

Definition: Let $i=1,\ldots,M$ index independent clusters, and let $j=1,\ldots,n_i$ index observations within cluster $i$. Write $Y_{ij}$ for the response and $x_{ij}$ for its $p$-dimensional covariate vector. The response vector is $\mathbf{Y}_i=(Y_{i1},\ldots,Y_{in_i})^\mathsf{T}$, and the covariate matrix $\mathbf{X}_i$ is the $n_i$-by-$p$ matrix whose $j$th row is $x_{ij}^\mathsf{T}$. The total number of observed responses is $N=\sum_{i=1}^M n_i$.

The superscript $\mathsf{T}$ denotes transpose. A covariate vector may include an intercept, treatment indicators, time, or other predictors. The response and covariate dimensions are tied to the cluster size: $\mathbf{Y}_i$ has $n_i$ entries and $\mathbf{X}_i$ has $n_i$ rows. Cluster sizes need not be equal.

Marginal Means and Their Derivatives

For a parameter vector $\beta$ of length $p$, a marginal mean model specifies $\mu_{ij}=E[Y_{ij}\mid x_{ij}]$. In the usual generalized linear model form, the linear predictor is $\eta_{ij}=x_{ij}^\mathsf{T}\beta$, and a link function $g$ relates it to the mean: $g(\mu_{ij})=\eta_{ij}$. Equivalently, $\mu_{ij}=g^{-1}(x_{ij}^\mathsf{T}\beta)$. The word marginal emphasizes that each expectation describes one response, not a conditional expectation given the other responses in its cluster.

Definition: The cluster marginal mean vector is $\boldsymbol{\mu}_i=(\mu_{i1},\ldots,\mu_{in_i})^\mathsf{T}$. Its derivative matrix is $\mathbf{D}_i=\partial\boldsymbol{\mu}_i/\partial\beta^\mathsf{T}$, an $n_i$-by-$p$ matrix. Row $j$ is the derivative of $\mu_{ij}$ with respect to the $p$ components of $\beta$.

When the inverse link is differentiable, the chain rule gives a useful row-by-row expression. Define $\dot{\mu}_{ij}=d g^{-1}(\eta_{ij})/d\eta_{ij}$. Then row $j$ of $\mathbf{D}_i$ is $\dot{\mu}_{ij}x_{ij}^\mathsf{T}$. This derivative records how a change in each regression parameter changes that observation’s marginal mean.

Result: If each $g^{-1}$ is differentiable at $\eta_{ij}=x_{ij}^\mathsf{T}\beta$, then $$ \mathbf{D}_i= \begin{pmatrix} \dot{\mu}_{i1}x_{i1}^\mathsf{T}\\ \vdots\\ \dot{\mu}_{in_i}x_{in_i}^\mathsf{T} \end{pmatrix}. $$ Derivation: For each $j$, $\mu_{ij}=g^{-1}(x_{ij}^\mathsf{T}\beta)$. Differentiating with respect to $\beta^\mathsf{T}$, the chain rule gives $$ \frac{\partial\mu_{ij}}{\partial\beta^\mathsf{T}} =\frac{d g^{-1}(\eta_{ij})}{d\eta_{ij}} \frac{\partial(x_{ij}^\mathsf{T}\beta)}{\partial\beta^\mathsf{T}} =\dot{\mu}_{ij}x_{ij}^\mathsf{T}. $$ Stacking these $n_i$ rows gives the displayed matrix and its stated dimensions. $\square$

Worked Example: Binary Responses at Two Visits

Suppose one patient has binary responses at two visits, with a logistic marginal mean model using an intercept and a visit indicator $t$. Let $\beta=(0,-\log 2)^\mathsf{T}$, with $t=0$ at visit 1 and $t=1$ at visit 2. Thus $x_{i1}=(1,0)^\mathsf{T}$, $x_{i2}=(1,1)^\mathsf{T}$, and $$ \mathbf{X}_i= \begin{pmatrix} 1&0\\ 1&1 \end{pmatrix}. $$

The logistic inverse link is $g^{-1}(\eta)=e^\eta/(1+e^\eta)$. The linear predictors are $0$ and $-\log 2$, so the means are $\mu_{i1}=1/2$ and $\mu_{i2}=1/3$. Therefore $\boldsymbol{\mu}_i=(1/2,1/3)^\mathsf{T}$. For the logistic inverse link, the derivative with respect to $\eta$ is $\mu(1-\mu)$. The two derivative rows are consequently $(1/4)(1,0)$ and $(2/9)(1,1)$, giving $$ \mathbf{D}_i= \begin{pmatrix} 1/4&0\\ 2/9&2/9 \end{pmatrix}. $$

For a Bernoulli response the marginal variance is $\mu_{ij}(1-\mu_{ij})$, so the variance entries are $1/4$ and $2/9$. If the chosen working correlation between visits is $0.3$, the corresponding covariance is $0.3\sqrt{(1/4)(2/9)}=0.3/\sqrt{18}=1/(10\sqrt{2})$, approximately $0.07071$. Thus the working covariance is $$ \mathbf{V}_i= \begin{pmatrix} 1/4&1/(10\sqrt{2})\\ 1/(10\sqrt{2})&2/9 \end{pmatrix}. $$ The correlation matrix has determinant $1-0.3^2=0.91>0$, so it is positive definite; scaling its rows and columns by the positive standard deviations preserves that property.

Variance Functions and Working Covariance

A marginal mean model alone does not describe the variability of the responses. For a response family with variance function $v$, a common specification is $\operatorname{Var}(Y_{ij})=\phi v(\mu_{ij})$, where $\phi>0$ is a dispersion parameter. For Bernoulli and Poisson responses the usual variance functions are $v(\mu)=\mu(1-\mu)$ and $v(\mu)=\mu$, respectively, with $\phi=1$ in their standard forms. For continuous responses, the variance function and dispersion convention depend on the chosen model.

Let $\mathbf{A}_i$ be the diagonal matrix of marginal variances, with $j$th diagonal entry $\phi v(\mu_{ij})$. A working correlation matrix $\mathbf{R}_i(\alpha)$ specifies proposed pairwise correlations within cluster $i$, possibly through a low-dimensional parameter $\alpha$. The resulting working covariance matrix is $$ \mathbf{V}_i=\mathbf{A}_i^{1/2}\mathbf{R}_i(\alpha)\mathbf{A}_i^{1/2}. $$ Its diagonal entries are the specified marginal variances; its off-diagonal entries are those variances’ standard deviations multiplied by the corresponding working correlations. The word “working” matters: $\mathbf{V}_i$ is the covariance structure used for the analysis, and need not equal the true covariance matrix.

Result: Suppose all specified marginal variances are positive and $\mathbf{R}_i(\alpha)$ is a positive semidefinite correlation matrix. Then $\mathbf{V}_i$ is positive semidefinite and has diagonal entries $\phi v(\mu_{ij})$. If $\mathbf{R}_i(\alpha)$ is positive definite, then $\mathbf{V}_i$ is positive definite. Proof: The matrix $\mathbf{A}_i^{1/2}$ is diagonal with entries $\sqrt{\phi v(\mu_{ij})}$. For any $n_i$-vector $a$, $$ a^\mathsf{T}\mathbf{V}_i a =(\mathbf{A}_i^{1/2}a)^\mathsf{T} \mathbf{R}_i(\alpha)(\mathbf{A}_i^{1/2}a)\geq 0, $$ by positive semidefiniteness of $\mathbf{R}_i(\alpha)$. The $j$th diagonal entry of $\mathbf{V}_i$ is $\sqrt{\phi v(\mu_{ij})}\times 1\times\sqrt{\phi v(\mu_{ij})}=\phi v(\mu_{ij})$. If $\mathbf{R}_i(\alpha)$ is positive definite, then for nonzero $a$, $\mathbf{A}_i^{1/2}a\ne0$, so the same quadratic form is strictly positive. $\square$

Worked Example: Working Covariance for Clustered Counts

Consider a cluster with three count responses and a Poisson log-link marginal model. Let the covariate be time $t=0,1,2$, with $x_{ij}=(1,t)^\mathsf{T}$ and $\beta=(\log 2,\log 2)^\mathsf{T}$. The covariate matrix and linear predictors are $$ \mathbf{X}_i= \begin{pmatrix} 1&0\\ 1&1\\ 1&2 \end{pmatrix}, \qquad \eta_{i1}=\log 2,\quad \eta_{i2}=\log 4,\quad \eta_{i3}=\log 8. $$

The means are $2,4,8$, so $\boldsymbol{\mu}_i=(2,4,8)^\mathsf{T}$. The Poisson variance function is $v(\mu)=\mu$, hence $\mathbf{A}_i=\operatorname{diag}(2,4,8)$ when $\phi=1$. Because the derivative of the exponential inverse link is $\mu$, the derivative matrix has rows $\mu_{ij}x_{ij}^\mathsf{T}$: $$ \mathbf{D}_i= \begin{pmatrix} 2&0\\ 4&4\\ 8&16 \end{pmatrix}. $$

Choose an exchangeable working correlation with all distinct-time correlations equal to $0.2$. Its matrix is $$ \mathbf{R}_i= \begin{pmatrix} 1&0.2&0.2\\ 0.2&1&0.2\\ 0.2&0.2&1 \end{pmatrix}. $$ Its eigenvalues are $0.8,0.8,1.4$, so it is positive definite. The working covariance has diagonal $2,4,8$. Its off-diagonal entries are $0.2\sqrt{2\cdot4}=0.4\sqrt{2}\approx0.5657$, $0.2\sqrt{2\cdot8}=0.8$, and $0.2\sqrt{4\cdot8}=0.8\sqrt{2}\approx1.1314$. Thus $$ \mathbf{V}_i= \begin{pmatrix} 2&0.4\sqrt{2}&0.8\\ 0.4\sqrt{2}&4&0.8\sqrt{2}\\ 0.8&0.8\sqrt{2}&8 \end{pmatrix}. $$ The off-diagonal entries are covariances, not correlations: for example, the correlation for the first and third responses is $0.8/\sqrt{2\cdot8}=0.2$.

Independence Across Clusters, Dependence Within Clusters

The usual GEE sampling structure treats $\mathbf{Y}_1,\ldots,\mathbf{Y}_M$ as independent random vectors. This is a statement about the joint distributions of whole clusters. It does not assert that $Y_{i1},\ldots,Y_{in_i}$ are independent within a cluster. Indeed, the purpose of a non-diagonal working covariance is to represent their within-cluster dependence.

A consequence of independence across clusters is that responses from different clusters have zero cross-covariance, provided their second moments are finite. To see this, take $i\ne k$. Independence gives $E[Y_{ij}Y_{k\ell}]=E[Y_{ij}]E[Y_{k\ell}]$, and therefore $$ \operatorname{Cov}(Y_{ij},Y_{k\ell}) =E[Y_{ij}Y_{k\ell}]-E[Y_{ij}]E[Y_{k\ell}]=0. $$ Accordingly, the covariance matrix of the stacked responses has cluster covariance blocks and zero off-block entries when the true second moments exist. Zero cross-covariance alone does not establish independence in general; the independent-cluster assumption is stronger.

Worked Example: Unequal Cluster Sizes and Block Structure

Suppose there are two independent clusters: a patient observed twice and a patient observed once. Their response vectors have dimensions two and one, respectively. For illustration, assume the first cluster’s true covariance is $$ \boldsymbol{\Sigma}_1= \begin{pmatrix} 1&0.4\\ 0.4&1 \end{pmatrix}, $$ and the second cluster’s variance is $\boldsymbol{\Sigma}_2=(2)$. If the responses are stacked cluster by cluster, independence makes the full covariance matrix $$ \operatorname{Cov} \begin{pmatrix} \mathbf{Y}_1\\ \mathbf{Y}_2 \end{pmatrix} = \begin{pmatrix} 1&0.4&0\\ 0.4&1&0\\ 0&0&2 \end{pmatrix}. $$

The $0.4$ entries describe within-cluster covariance for the first patient; the zeros between blocks follow from independence of the patients’ response vectors. The second patient’s cluster size is one, so there is no within-cluster off-diagonal entry to specify. This example also shows why a single common response-vector dimension is unnecessary: each block is sized according to its own $n_i$.

How to Check the Inputs

Before proceeding to any cluster-level calculation, check the objects in a fixed order. The response vector and mean vector should each have $n_i$ entries; the covariate matrix should have $n_i$ rows and $p$ columns; the derivative matrix should be $n_i$-by-$p$; and the working covariance should be $n_i$-by-$n_i$. The covariance diagonal should match the stated marginal variance function, and the working correlation should have unit diagonal and be a valid positive semidefinite correlation matrix.

Index the observations.
Record the cluster index $i$ and the within-cluster index $j$, allowing $n_i$ to vary across clusters.

Specify mean inputs.
Write the covariate rows, apply the link-scale model, and collect the fitted marginal means and their derivatives.

Specify variability.
Evaluate the marginal variance function at each mean, then combine those variances with a valid within-cluster working correlation.

Check cluster assumptions.
Confirm that the independent units are clusters and that any proposed dependence is represented within, not between, clusters.

Common Mistakes

Treating every response as an independent sampling unit. Repeated responses from the same patient can be dependent. The independent units in the standard setup are the clusters.
Confusing a correlation with a covariance. A working correlation has unit diagonal; the working covariance rescales it by the marginal standard deviations.
Using a single common cluster size without checking. The matrices for cluster $i$ must use its actual $n_i$, which may differ from the size of another cluster.
Assuming a working covariance is automatically the true covariance. It is a chosen structure for the analysis; its diagonal and correlation interpretation should be explicit, but it need not describe the true within-cluster dependence exactly.

Key takeaway: For each cluster, organize the responses, covariates, marginal means, derivatives, and working covariance with dimensions determined by that cluster’s size. Independence applies across the response vectors; dependence may remain within each vector.

Check Your Understanding

Use the definitions and dimension checks in this tutorial to answer the following.

If cluster $i$ has four responses and the model has three regression parameters, what are the dimensions of $\mathbf{Y}_i$, $\mathbf{X}_i$, $\mathbf{D}_i$, and $\mathbf{V}_i$?
For a logistic mean model, what is the derivative of $\mu_{ij}$ with respect to the linear predictor, expressed using $\mu_{ij}$?
Explain why an off-diagonal working covariance entry is not generally equal to its working correlation.
What does independence across clusters imply about cross-covariances, and why does zero cross-covariance alone not prove independence?
For a Poisson marginal model with mean $\mu_{ij}$ and $\phi=1$, what is the specified marginal variance?