Score-Type Tests Under GEE

Why Evaluate the Estimating Function Under the Null?

A Wald test begins with an unrestricted GEE fit and asks how far the estimated contrast lies from its null value, relative to its estimated covariance. A score-type test takes a different route: fit the model with the null restriction imposed, then assess whether the tested estimating equations still have a substantial residual at that constrained fit. The calculation can be useful when the null model is easier to fit than the unrestricted model.

The word “score” needs care here. GEE supplies estimating functions, not in general derivatives of a log-likelihood. A GEE score-type test therefore uses the estimating function and its robust covariance; it is not automatically the classical likelihood score test.

Definition and Nuisance Adjustment

Suppose independent, identically distributed (i.i.d.) clusters are indexed by \(i=1,\ldots,n\). Write the mean-model parameter as \(\theta=(\beta^T,\gamma^T)^T\), where \(\beta\) is the \(q\)-dimensional parameter of interest and \(\gamma\) is a nuisance parameter. We test \(H_0:\beta=\beta_0\). Let \(U_i(\theta)\) be the cluster’s GEE contribution and \(U(\theta)=\sum_{i=1}^n U_i(\theta)\). Partition these into \(U_{i\beta}\) and \(U_{i\gamma}\), and likewise for their sums.

Definition: A GEE score-type test evaluates \(U_\beta\) at the constrained fit \((\beta_0,\hat\gamma_0)\), where the nuisance estimate solves \(U_\gamma(\beta_0,\hat\gamma_0)=0\). Its covariance must account for estimating \(\gamma\) under the constraint. The resulting quadratic statistic is compared with a chi-square reference distribution under the regularity conditions stated below.

Let \(A(\theta)\) denote the sensitivity matrix, the negative derivative of the mean estimating function:

A(\theta)=-E\!\left[\frac{\partial U_i(\theta)}{\partial\theta^T}\right], \qquad A(\theta)= \begin{pmatrix} A_{\beta\beta}&A_{\beta\gamma}\\ A_{\gamma\beta}&A_{\gamma\gamma} \end{pmatrix}.

The expectation here is for one cluster; the same partition applies to the empirical sum of negative derivatives. Assume \(A_{\gamma\gamma}\) is nonsingular. Define \(C=A_{\beta\gamma}A_{\gamma\gamma}^{-1}\). The nuisance-adjusted cluster contribution is

\psi_i=U_{i\beta}-C U_{i\gamma}, \qquad \Omega=\operatorname{Var}(\psi_i).

In practice, estimate \(C\) from the fitted sensitivity matrix and estimate \(\Omega\) from the outer products of the fitted adjusted cluster contributions. At the constrained fit, the total nuisance equation is zero, so \(\sum_i\psi_i=U_\beta(\beta_0,\hat\gamma_0)\). The adjustment still matters: it determines the variability of that total after the nuisance parameter has been estimated.

Derivation and Asymptotic Result

Theorem: Suppose clusters are independent and identically distributed, the number of parameters is fixed, and the null model is correctly specified in the sense that \(E[U_i(\theta_0)]=0\), where \(\theta_0=(\beta_0^T,\gamma_0^T)^T\). Assume the estimating functions are differentiable near \(\theta_0\), the constrained nuisance estimator is consistent, \(A_{\gamma\gamma}\) is nonsingular, and the cluster estimating functions satisfy a central limit theorem with finite nonsingular covariance for \(\psi_i\). Also assume the first-order expansions in Steps 1 and 2 below hold, with \(n^{-1}A_{\gamma\gamma}^{(n)}\xrightarrow{p}A_{\gamma\gamma}\) and \(n^{-1}A_{\beta\gamma}^{(n)}\xrightarrow{p}A_{\beta\gamma}\). Also assume consistent estimators of \(C\) and \(\Omega\) are available, with \(\Omega\) positive definite. Then \(T=U_\beta(\beta_0,\hat\gamma_0)^T\widehat{\operatorname{Var}}\{U_\beta(\beta_0,\hat\gamma_0)\}^{-1}U_\beta(\beta_0,\hat\gamma_0)\) converges in distribution under \(H_0\) to \(\chi_q^2\).

The covariance in this theorem is the covariance of the total adjusted estimating function, estimated by summing cluster outer products. In particular, if \(\hat\psi_i=U_{i\beta}(\beta_0,\hat\gamma_0)-\hat C U_{i\gamma}(\beta_0,\hat\gamma_0)\), use \(\hat V=\sum_i\hat\psi_i\hat\psi_i^T\). With fixed \(q\), no extra division by \(n\) belongs in the quadratic form when both the total estimating function and its total covariance are used.

Expand the nuisance equation.
At the true null value, a first-order Taylor expansion gives \(0=U_\gamma(\beta_0,\hat\gamma_0)=U_\gamma(\theta_0)-A_{\gamma\gamma}^{(n)}(\hat\gamma_0-\gamma_0)+o_p(\sqrt n)\), where \(A_{\gamma\gamma}^{(n)}\) is the summed sensitivity block. Thus \(\hat\gamma_0-\gamma_0=(A_{\gamma\gamma}^{(n)})^{-1}U_\gamma(\theta_0)+o_p(n^{-1/2})\).

Expand the tested equation.
The corresponding expansion is \(U_\beta(\beta_0,\hat\gamma_0)=U_\beta(\theta_0)-A_{\beta\gamma}^{(n)}(\hat\gamma_0-\gamma_0)+o_p(\sqrt n)\). Substituting Step 1 gives the sum of \(U_{i\beta}-C U_{i\gamma}\), up to \(o_p(\sqrt n)\), where \(C=A_{\beta\gamma}A_{\gamma\gamma}^{-1}\).

Apply the central limit theorem.
The adjusted cluster contributions have mean zero under the null and covariance \(\Omega\). Independence and the stated moment and nonsingularity conditions imply \(n^{-1/2}\sum_i\psi_i\xrightarrow{d}N(0,\Omega)\).

Studentize the result.
Consistency gives \(n^{-1}\hat V\xrightarrow{p}\Omega\). Slutsky’s theorem then implies that the quadratic form in the theorem converges to the sum of squares of \(q\) independent standard normal variables, which is \(\chi_q^2\). \(\square\)

The first-order expansion also explains why simply using the unadjusted variability of \(U_\beta\) can be wrong: the constrained nuisance fit responds to random variation in \(U_\gamma\), and that response changes the tested equation. For a scalar test, the statistic is the squared constrained estimating function divided by its estimated variance. For a joint test, the full covariance matrix is essential.

Worked Calculations

Worked Example: A Scalar Test With a Nuisance Intercept

Suppose a restricted fit gives three cluster contributions \(U_{i\beta}=(2,-1,0)\) and \(U_{i\gamma}=(1,-1,0)\). The nuisance equation is satisfied because \(1-1+0=0\). Suppose the fitted sensitivity blocks are \(A_{\beta\gamma}=2\) and \(A_{\gamma\gamma}=4\), so \(C=2/4=0.5\).

The adjusted contributions are \(2-0.5(1)=1.5\), \(-1-0.5(-1)=-0.5\), and \(0-0.5(0)=0\). Their total is \(1\), equal to the constrained tested estimating function \(2-1+0=1\). The cluster outer-product estimate of its variance is \(1.5^2+(-0.5)^2+0^2=2.5\). The score-type statistic is therefore \(T=1^2/2.5=0.4\), compared with a \(\chi_1^2\) reference distribution.

This is an illustrative calculation, not a claim that three clusters suffice for reliable asymptotic inference. It shows the mechanics of nuisance adjustment and the total-covariance scaling.

Worked Example: A Joint Test With Two Tested Parameters

At a restricted fit, suppose three cluster contributions are \(U_{i\beta}=(2.5,-0.25)^T,\;(-0.5,1.25)^T,\;(0,-1)^T\) and \(U_{i\gamma}=1,-1,0\). The nuisance contributions sum to zero. Let \(C=(0.5,-0.25)^T\), so the nuisance-adjusted cluster contributions are \(U_{i\beta}-C U_{i\gamma}\).

For clusters 1, 2, and 3 these are, respectively, \((2,0)^T\), \((0,1)^T\), and \((0,-1)^T\). Their total is \(R=(2,0)^T\), which is also the sum of the tested contributions. The covariance estimate is

\hat V= \begin{pmatrix}2\\0\end{pmatrix} \begin{pmatrix}2&0\end{pmatrix} + \begin{pmatrix}0\\1\end{pmatrix} \begin{pmatrix}0&1\end{pmatrix} + \begin{pmatrix}0\\-1\end{pmatrix} \begin{pmatrix}0&-1\end{pmatrix} = \begin{pmatrix}4&0\\0&2\end{pmatrix}.

Hence \(T=R^T\hat V^{-1}R=(2,0)\begin{pmatrix}1/4&0\\0&1/2\end{pmatrix}(2,0)^T=1\). Under the null and the theorem’s conditions, the reference distribution is \(\chi_2^2\). The two-dimensional covariance, rather than two separate univariate calculations, is what defines the joint test.

Worked Example: A Treatment Test in a Clustered Linear GEE

Consider four independent clusters of two observations each. Fit an identity-link mean model with mean \(\gamma+\beta x_i\), where \(x_i=0\) for the first two clusters and \(x_i=1\) for the last two. Use working independence and unit variance. Test \(H_0:\beta=0\), treating the intercept \(\gamma\) as nuisance. The observed responses in the four clusters are \((9,9),(10,10),(10.5,10.5),(10.5,10.5)\).

Under the null, the restricted intercept is the grand mean, \(\hat\gamma_0=80/8=10\). The cluster residual sums are \(-2,0,1,1\), so the nuisance equation sums to zero. The treatment estimating-function contributions are \(0,0,1,1\), giving constrained tested total \(R=2\). The sensitivity blocks are \(A_{\beta\gamma}=4\) (four treated observations) and \(A_{\gamma\gamma}=8\) (eight observations), hence \(C=4/8=0.5\).

The adjusted cluster contributions are \(0-0.5(-2)=1\), \(0-0.5(0)=0\), \(1-0.5(1)=0.5\), and \(1-0.5(1)=0.5\). Their squared sum is \(1+0+0.25+0.25=1.5\). Thus \(T=2^2/1.5=8/3\approx2.667\), with a \(\chi_1^2\) reference distribution. This example makes explicit that the nuisance adjustment can alter the robust variance even though the total adjusted contribution equals the tested total at the restricted fit.

How This Differs From a Likelihood Score Test

In a regular parametric likelihood model, the classical score test evaluates the derivative of the log-likelihood at the null-constrained maximum likelihood estimate of nuisance parameters. Its covariance is based on Fisher information, with likelihood identities linking score variance and expected curvature under the model. Those facts yield the familiar likelihood score statistic.

A GEE contribution is generally not a log-likelihood derivative. Its sensitivity and its covariance need not be equal, particularly when the working correlation is misspecified. The GEE test therefore uses the sandwich covariance of the nuisance-adjusted estimating functions. Calling it “score-type” describes its constrained-fit construction; it does not supply a likelihood or justify replacing robust covariance with model-based information.

The theorem treats the working-correlation structure as fixed or as a plug-in whose estimation has no first-order effect on the mean estimating function under the usual GEE conditions. If estimation of additional nuisance parameters does affect the first-order behavior, include their estimating equations and sensitivity blocks in the joint expansion before forming the adjusted covariance.

Common Pitfall

Wrong reasoning: “The nuisance equation is zero at the constrained fit, so I can ignore nuisance estimation and use only the variability of the tested contributions.”
Correction: The nuisance equation’s total is zero, but its cluster contributions are not generally zero. Estimation of nuisance parameters changes the tested equation through \(C U_{i\gamma}\). Use adjusted cluster contributions to estimate the covariance, while recognizing that their total equals the constrained tested equation.
Another error: Treating a GEE score-type statistic as a likelihood score statistic and using a likelihood-information covariance. GEE generally has no likelihood score or information identity; use the appropriate robust covariance.

Key takeaway: Fit the null-constrained GEE, evaluate the tested estimating function there, and studentize it with the cluster-robust covariance of nuisance-adjusted contributions. The chi-square calibration is asymptotic and depends on independent clusters, regularity, and consistent covariance estimation.

Check Your Understanding

Use the definitions and derivation above to answer these questions.

Why does the total adjusted estimating function equal the constrained tested estimating function when the nuisance equation is solved?
For scalar \(\beta\), write the score-type statistic using the constrained tested total and its estimated total variance.
What role does the sensitivity block \(A_{\beta\gamma}\) play in nuisance adjustment?
Why is a likelihood-information covariance generally not justified for a GEE score-type test?
In the treatment example, what would be lost by calculating covariance from the unadjusted treatment contributions alone?

What You'll Learn