Deriving GEE from the Multivariate Quasi-Score

From Marginal Information to an Estimating Equation

For correlated responses, a cluster supplies a vector of observations, a model for its marginal means, and a working description of its covariance. Those ingredients suffice to construct a generalized estimating equation (GEE); a joint probability model for the response vector is not required. The key operation is to transform the vector of residuals using the derivative of the mean and the inverse working covariance.

Suppose there are $n$ independent clusters. Cluster $i$ has response vector $Y_i$ of length $m_i$, with marginal mean $\mu_i(\theta)=E[Y_i]$. The parameter $\theta$ has dimension $p$. The mean may be specified component by component through a link and covariates. For this derivation, assume the mean is differentiable in a neighborhood of the true value $\theta_0$.

Definition: The mean derivative matrix for cluster $i$ is $D_i(\theta)=\partial\mu_i(\theta)/\partial\theta'$, an $m_i\times p$ matrix. Let $V_i(\theta)$ be a specified $m_i\times m_i$ working covariance matrix, symmetric and positive definite. The cluster multivariate quasi-score is $$ U_i(\theta)=D_i(\theta)'V_i(\theta)^{-1}\{Y_i-\mu_i(\theta)\}. $$ The GEE is $\sum_{i=1}^n U_i(\theta)=0$.

The residual vector has length $m_i$. Multiplication by $V_i^{-1}$ transforms that vector using the working covariance, and multiplication by $D_i'$ maps it into $p$ parameter directions. Thus each cluster contributes a $p$-vector, just as many estimating equations as there are parameters. The covariance may include marginal variances and a working association structure; it need not be the true covariance.

Why This Form Appears

A direct route to the formula is local weighted least squares. Fix a current parameter value $\theta$, and consider a small change $\delta$ in the parameter. Differentiability gives the first-order approximation $\mu_i(\theta+\delta)\approx\mu_i(\theta)+D_i(\theta)\delta$. If the current residual is $r_i=Y_i-\mu_i(\theta)$, the residual after this change is approximately $r_i-D_i\delta$. Measure the discrepancy with the working covariance by the quadratic form $(r_i-D_i\delta)'V_i^{-1}(r_i-D_i\delta)$.

Here the matrices are evaluated at the current parameter value and held fixed while choosing the local change. Differentiating this quadratic form with respect to $\delta$, setting the derivative to zero, and dropping the common factor $-2$ gives the normal equations \[ D_i'V_i^{-1}(r_i-D_i\delta)=0. \] The term multiplying the current residual is precisely the cluster quasi-score. Summing over clusters gives the GEE. This local argument motivates the equation; it does not assert that the quadratic form is a negative log-likelihood.

Key distinction: The GEE is defined from a mean, a mean derivative, and a working covariance. It is not obtained by differentiating a joint response log-likelihood. In particular, the construction does not require a joint density or mass function for $Y_i$.

Worked Examples

Worked Example: Two Marginal Binary Means

Consider a cluster with two binary responses and a saturated marginal mean model. Let $\theta=(\theta_1,\theta_2)'$, with $\operatorname{logit}\{\mu_{ij}(\theta)\}=\theta_j$ for occasion $j$. Evaluate the cluster contribution at $\theta=(0,0)'$, for observed responses $Y_i=(1,0)'$. Use a working correlation of $0.4$.

At $\theta=(0,0)'$, each marginal mean is $1/2$, and the derivative of a logistic mean with respect to its logit is $\mu(1-\mu)=1/4$. Therefore $\mu_i=(1/2,1/2)'$ and $D_i=(1/4)I_2$. The working variances are both $1/4$; the working covariance is $0.4\sqrt{(1/4)(1/4)}=1/10$. Hence

V_i= \begin{pmatrix} 1/4&1/10\\ 1/10&1/4 \end{pmatrix}, \qquad V_i^{-1}= \begin{pmatrix} 100/21&-40/21\\ -40/21&100/21 \end{pmatrix}.

The residual is $r_i=(1/2,-1/2)'$. Multiplying in order, $V_i^{-1}r_i=(10/3,-10/3)'$, and therefore $U_i=D_i'V_i^{-1}r_i=(5/6,-5/6)'$. The positive association changes how the two residuals combine, while $D_i'$ expresses that adjusted residual in the two marginal-logit parameter directions.

Worked Example: A Two-Occasion Count Model

Let $Y_i$ contain two counts, with log means $\log\mu_{i1}=\theta_1$ and $\log\mu_{i2}=\theta_1+\theta_2$. At $\theta=(\log 2,\log 2)'$, the mean vector is $(2,4)'$, and $D_i=\begin{pmatrix}2&0\\4&4\end{pmatrix}$. Suppose the observed counts are $(3,2)'$, and use the positive definite working covariance $V_i=\begin{pmatrix}2&1\\1&4\end{pmatrix}$.

The residual is $(1,-2)'$. Since the determinant of $V_i$ is $2(4)-1=7$,

V_i^{-1}=\frac{1}{7} \begin{pmatrix} 4&-1\\ -1&2 \end{pmatrix}, \qquad V_i^{-1}(Y_i-\mu_i)= \frac{1}{7} \begin{pmatrix} 6\\-5 \end{pmatrix}.

Consequently, $U_i=D_i'V_i^{-1}(Y_i-\mu_i)$, so $U_{i1}=2(6/7)+4(-5/7)=-8/7$ and $U_{i2}=4(-5/7)=-20/7$. These are the two parameter-direction contributions from this cluster. The calculation uses the marginal means and the chosen covariance matrix; it does not require specifying a joint count distribution.

Worked Example: Continuous Responses and a Linear Mean

Suppose a two-response cluster has mean $\mu_i=X_i\theta$, where $X_i=\begin{pmatrix}1&0\\1&1\end{pmatrix}$. At $\theta=(1,2)'$, the mean is $(1,3)'$, and $D_i=X_i$. For observed responses $Y_i=(2,5)'$, the residual is $(1,2)'$. Take $V_i=\begin{pmatrix}1&1/2\\1/2&4\end{pmatrix}$.

The determinant is $4-1/4=15/4$, so $V_i^{-1}=\begin{pmatrix}16/15&-2/15\\-2/15&4/15\end{pmatrix}$. It follows that $V_i^{-1}(Y_i-\mu_i)=(4/5,2/5)'$. Applying the transposed mean derivative gives

U_i= \begin{pmatrix}1&1\\0&1\end{pmatrix} \begin{pmatrix}4/5\\2/5\end{pmatrix} = \begin{pmatrix}6/5\\2/5\end{pmatrix}.

The first parameter shifts both means, so its contribution adds the two adjusted residuals. The second shifts only the second mean, so its contribution is the second adjusted residual alone.

Two Properties That Do Not Require a Joint Likelihood

The first property is unbiasedness at the true parameter. It requires the marginal mean to be correct and the working weight to be well-defined; it does not require the working covariance to equal the true covariance.

Result: If $E[Y_i]=\mu_i(\theta_0)$, and $D_i(\theta_0)$ and $V_i(\theta_0)^{-1}$ are finite nonrandom matrices, then $E[U_i(\theta_0)]=0$. If clusters are independent, the sum of the cluster contributions also has expectation zero.

Proof. At $\theta_0$, the matrices in the quasi-score are fixed, and the residual has mean zero. By linearity of expectation,

E[U_i(\theta_0)] =D_i(\theta_0)'V_i(\theta_0)^{-1} E[Y_i-\mu_i(\theta_0)] =D_i(\theta_0)'V_i(\theta_0)^{-1}0=0.

Adding these zero expectations over the finite number of clusters gives $E[\sum_i U_i(\theta_0)]=0$. No step uses the distribution of the residual beyond its mean. $\square$

A second property explains the covariance weighting when it is correct. It is a local optimality result within the class of linear estimating functions of the cluster residual. It does not say that every working GEE is globally optimal, nor that a misspecified working covariance has this optimality property.

Result (optimal linear estimating function): Fix a parameter value, and suppose $Y_i$ has mean $\mu_i$, positive definite covariance $V_i$, and mean derivative $D_i$. Among estimating functions $A_i(Y_i-\mu_i)$, where $A_i$ is a fixed $p\times m_i$ matrix and $A_iD_i$ is nonsingular, the choice $A_i=D_i'V_i^{-1}$ has Godambe information at least as large as any other such choice, in the positive-semidefinite matrix order.

Proof. Write $r_i=Y_i-\mu_i$. For $g_i=A_i r_i$, its sensitivity matrix is $H_i=A_iD_i$, since differentiating the mean in the residual gives $-D_i$. Its covariance is $J_i=A_iV_iA_i'$. The Godambe information for this estimating function is $H_i'J_i^{-1}H_i$, which is defined under the stated nonsingularity conditions.

Let $C=A_iV_i^{1/2}$ and $K=V_i^{-1/2}D_i$, using the symmetric positive definite square root. Then $J_i=CC'$, $H_i=CK$, and $H_i'J_i^{-1}H_i=K'C'(CC')^{-1}CK$. Because $C$ has full row rank, $P=C'(CC')^{-1}C$ is the orthogonal projection onto the row space of $C$. Thus $P$ and $I-P$ are positive semidefinite, and

D_i'V_i^{-1}D_i-H_i'J_i^{-1}H_i =K'(I-P)K \succeq 0.

For $A_i=D_i'V_i^{-1}$, we have $C=K'$, so the columns of $K$ lie in the row space of $C$, giving equality. This proves the stated comparison. $\square$

What the Derivative and Covariance Are Doing

The mean derivative tells the equation how a small parameter change moves each component of the marginal mean. The working covariance determines how the residual vector is transformed before those componentwise changes are mapped back to parameter space. Both are essential: without $D_i$, the equation has no direct connection to the parameters; without $V_i^{-1}$, it omits the chosen covariance weighting.

For a working covariance that varies with $\theta$, the GEE is still formed by evaluating $D_i(\theta)$, $V_i(\theta)^{-1}$, and $Y_i-\mu_i(\theta)$ at each candidate value and setting the sum of the resulting contributions to zero. This is a definition of an estimating equation, not a claim that one has differentiated every parameter-dependent term. A likelihood score can contain additional terms arising from the distribution and its covariance; the GEE does not acquire such terms merely because its working covariance depends on $\theta$.

Common Pitfall

A frequent error is to treat $D_i'V_i^{-1}(Y_i-\mu_i)$ as though it had to be the derivative of a joint log-likelihood. That reasoning may lead someone to demand a joint response distribution or to add covariance-derivative terms automatically. The correction is to identify the object being constructed: this is a moment-based estimating function. Its expectation is zero when the marginal mean is correct, and its form follows from the local linearized, covariance-weighted residual criterion. A joint likelihood is one possible route to other estimating functions, but it is not a prerequisite for the GEE.

Key takeaway: For each cluster, the quasi-score is the residual vector, adjusted by the inverse working covariance, then mapped through the transposed derivative of the marginal mean. Summing these contributions yields the GEE without specifying a joint response likelihood.

Check Your Understanding

Use the definitions and derivations above to answer the following questions.

For a cluster with $m_i$ responses and $p$ parameters, what are the dimensions of $D_i$, $V_i^{-1}$, and $U_i$?
Which assumption makes the expected cluster quasi-score zero at the true parameter?
In the local linearization, what role does the matrix $D_i$ play in the residual after a parameter change?
Does unbiasedness require the working covariance to equal the true covariance? Explain.
Why does the GEE formula not, by itself, require a joint probability distribution for the response vector?