Tutorials › Generalized Linear Models › Iterative Solution of the GEE

Home  /  Tutorials  /  Generalized Estimating Equations

Generalized Estimating Equations · Tutorial 853 of 981

Iterative Solution of the GEE

Learn to update GEE regression estimates systematically, refresh the quantities that define the working covariance, and assess whether an iteration has converged.

Advanced 9 min read

What You'll Learn

  • Write a Newton update for the summed GEE and distinguish it from a frozen-covariance scoring step
  • Recompute fitted means, derivatives, and working covariance matrices at each iterate
  • Carry out and check iterations for logistic and clustered Poisson mean models
  • Use the local Newton convergence result and identify when its conditions are not met
  • Choose practical stopping criteria based on both parameter changes and the estimating-equation residual

Why GEE Requires Iteration

A generalized estimating equation (GEE) is usually a system of nonlinear equations in the regression parameter \(\beta\). Even when the mean model is simple, its fitted means and working covariance depend on the current value of \(\beta\). Solving the equation therefore involves repeatedly updating the parameter and recalculating the quantities used to form the next update.

Let \(Y_i\) be the response vector for cluster \(i\), let \(\mu_i(\beta)=E[Y_i]\) be its marginal mean vector, and let \(D_i(\beta)=\partial\mu_i(\beta)/\partial\beta^T\). Write \(V_i(\beta)\) for the working covariance matrix, with any working-correlation or dispersion parameters held fixed during a regression-parameter step. The summed GEE is

$$ U(\beta)=\sum_{i=1}^m D_i(\beta)^T V_i(\beta)^{-1}\{Y_i-\mu_i(\beta)\}=0. $$

Here \(m\) is the number of independent clusters. The solution is generally not obtained by rearranging this equation: the mean, its derivative, and often the covariance all change with \(\beta\).

Definition: An iterative GEE solution is a sequence of regression estimates \(\beta^{(0)},\beta^{(1)},\ldots\) constructed from the current estimating function and its derivative or a specified approximation to that derivative. At each iterate, fitted means and the working covariance are evaluated at the current parameter before forming the next step.

Newton and Scoring Updates

For a vector-valued function \(U(\beta)\), define its Jacobian \(J(\beta)=\partial U(\beta)/\partial\beta^T\). If \(J(\beta)\) is nonsingular, a Newton step solves the linear approximation \(U(\beta+\Delta)\approx U(\beta)+J(\beta)\Delta=0\). Thus the exact Newton update is

$$ \beta^{(t+1)}=\beta^{(t)}-\{J(\beta^{(t)})\}^{-1}U(\beta^{(t)}). $$

The Jacobian can include derivatives of \(D_i\) and \(V_i\), as well as the derivative of the residual. A commonly used GEE scoring step instead freezes the working covariance while forming the step and uses the positive semidefinite matrix

$$ H(\beta)=\sum_{i=1}^m D_i(\beta)^T V_i(\beta)^{-1}D_i(\beta),\qquad \beta^{(t+1)}=\beta^{(t)}+H(\beta^{(t)})^{-1}U(\beta^{(t)}). $$

The displayed inverse-based step requires \(H\) to be nonsingular; for example, this holds if each \(V_i\) is positive definite and the vertically stacked \(D_i\) has full column rank. This is a Newton-style update, but it is not automatically the exact Newton step: the formula for \(H\) omits derivative terms that can occur when the covariance or mean derivative varies with \(\beta\). When \(H\) is used, it should be identified as the scoring or frozen-covariance step. In a scalar example below, the exact derivative and this matrix happen to differ, so the update is computed from the exact Jacobian.

Result: If \(U\) is continuously differentiable near a root \(\beta^\ast\), \(J(\beta^\ast)\) is nonsingular, and \(J\) is locally Lipschitz, exact Newton iteration converges quadratically when started sufficiently close to \(\beta^\ast\). More precisely, for some neighborhood and finite constant \(C\), its errors satisfy \(\|\beta^{(t+1)}-\beta^\ast\|\le C\|\beta^{(t)}-\beta^\ast\|^2\).

To verify the claim, set \(e=\beta-\beta^\ast\). Since \(U(\beta^\ast)=0\), the fundamental theorem of calculus gives \(U(\beta)=\int_0^1 J(\beta^\ast+s e)e\,ds\). The Newton error \(e_{\mathrm{new}}\) obeys

$$ e_{\mathrm{new}} =e-J(\beta)^{-1}U(\beta) =J(\beta)^{-1}\int_0^1\{J(\beta)-J(\beta^\ast+s e)\}e\,ds. $$

By continuity and nonsingularity, \(J(\beta)^{-1}\) is bounded in a sufficiently small neighborhood of \(\beta^\ast\). If \(J\) is Lipschitz there with constant \(L\), then \(\|J(\beta)-J(\beta^\ast+s e)\|\le L(1-s)\|e\|\). Consequently, if \(\|J(\beta)^{-1}\|\le M\), \(\|e_{\mathrm{new}}\|\le ML\|e\|^2\int_0^1(1-s)\,ds=(ML/2)\|e\|^2\). Choosing the starting error small enough keeps the iterates in that neighborhood and makes the errors decrease quadratically. This is a local result; it does not guarantee convergence from an arbitrary starting value. \(\square\)

Refresh the Mean and Working Covariance

At iteration \(t\), evaluate the linear predictors and fitted means from \(\beta^{(t)}\). For a GLM mean model, this means calculating \(\eta_i^{(t)}=X_i\beta^{(t)}\) and applying the inverse link to obtain \(\mu_i^{(t)}\). Then calculate \(D_i^{(t)}\), the variance function, and the working covariance used for that step. In the common representation \(V_i=\phi A_i^{1/2}R_i(\alpha)A_i^{1/2}\), \(A_i\) is determined by the current fitted marginal variances, while \(R_i(\alpha)\) is the working correlation matrix.

The distinction between updating \(\beta\) and refreshing the working covariance matters. If \(\beta\) changes, the fitted means and usually \(A_i\) change too; using an old \(V_i\) in the next equation is not the same as evaluating the specified GEE at the new parameter. Any working-correlation or dispersion estimates used in the procedure should likewise be refreshed according to the chosen algorithm, with their values held fixed for the regression step if that is how the step is defined.

1
Initialize.
Choose \(\beta^{(0)}\) and the working-correlation and dispersion values to use in the first step.
2
Evaluate at the current estimate.
Calculate every \(\mu_i^{(t)}\), \(D_i^{(t)}\), \(V_i^{(t)}\), and the summed \(U(\beta^{(t)})\).
3
Form and apply a step.
Use the exact Jacobian for Newton, or state and use a chosen scoring approximation. If the system is ill-conditioned, consider a damped step rather than taking the full increment.
4
Refresh and check.
Recalculate the fitted quantities at the new estimate. Check both the parameter change and the estimating-equation residual before declaring convergence.

Worked Iterations

Worked Example: Newton Iteration for a Logistic Intercept

Suppose 14 of 20 independent binary observations are successes, and the model gives each observation probability \(p=\{1+\exp(-\beta)\}^{-1}\). The summed equation is \(U(\beta)=14-20p\). Its derivative is \(J(\beta)=-20p(1-p)\), so the exact Newton increment is \(\Delta=U/\{20p(1-p)\}\).

Start at \(\beta^{(0)}=0\), giving \(p^{(0)}=0.5\), \(U^{(0)}=14-10=4\), and \(H^{(0)}=20(0.5)(0.5)=5\). Therefore \(\Delta^{(0)}=4/5=0.8\), and \(\beta^{(1)}=0.8\). Refresh the probability: \(p^{(1)}=1/(1+\exp(-0.8))\approx0.689974\). Then \(U^{(1)}=14-20(0.689974)\approx0.20052\), and \(H^{(1)}=20(0.689974)(0.310026)\approx4.2782\). The next increment is about \(0.04687\), giving \(\beta^{(2)}\approx0.84687\).

At \(\beta^{(2)}\), \(p\) is about \(0.69991\), so \(U\approx14-20(0.69991)=0.0018\). One more increment is about \(0.0018/4.200\approx0.00043\), bringing \(\beta\) close to \(\log(0.7/0.3)\approx0.84730\). The iteration is checked against the equation residual, not merely against a small-looking coefficient change.

Worked Example: Correct Newton Step for Clustered Poisson Counts

There are six clusters of size two, with a common Poisson mean \(\mu=\exp(\beta)\), working dispersion one, and exchangeable working correlation \(R\) with \(\rho=0.2\). The total of the twelve observed counts is 12. Thus \(V_i=\mu R\), \(D_i=\mu(1,1)^T\), and \(R^{-1}(1,1)^T=(1,1)^T/1.2\). The summed equation is

$$ U(\beta)=\frac{12-12\mu}{1.2}=10-10\mu. $$

Because \(d\mu/d\beta=\mu\), the exact negative derivative is \(H(\beta)=-dU/d\beta=10\mu\). Starting at \(\mu^{(0)}=1.5\), the equation gives \(U=10-15=-5\), \(H=15\), and \(\Delta=U/H=-1/3\). Therefore \(\beta^{(1)}=\beta^{(0)}-1/3\) and \(\mu^{(1)}=1.5\exp(-1/3)\approx1.07480\).

At the refreshed mean, \(U^{(1)}=10-10(1.07480)\approx-0.74797\), \(H^{(1)}=10(1.07480)\approx10.74797\), and \(\Delta^{(1)}\approx-0.06959\). This gives \(\mu^{(2)}\approx1.00254\). The next equation residual is \(U^{(2)}=10-10(1.00254)\approx-0.0254\), with \(H^{(2)}\approx10.0254\) and \(\Delta^{(2)}\approx-0.002534\); the next mean is approximately 1.00000. As a direct arithmetic check, at \(\mu=1.99969\) the equation gives \(U=10-10(1.99969)=-9.9969\).

The factor of \(\mu\) cancels between \(D_i\) and \(V_i^{-1}\) in the GEE contribution, but the residual still depends on \(\mu\). Differentiating that residual with respect to \(\beta\) gives the \(10\mu\) above. This calculation should not be replaced by inserting an extra \(1/\mu\) into \(U\).

Worked Example: Refreshing the Working Covariance

Continue with the clustered Poisson model above. For a size-two cluster and \(\rho=0.2\), \(R=\begin{pmatrix}1&0.2\\0.2&1\end{pmatrix}\). At the initial mean \(\mu=1.5\), the working covariance is \(V_i=1.5R=\begin{pmatrix}1.5&0.3\\0.3&1.5\end{pmatrix}\). After the first update, the mean is \(\mu\approx1.07480\). The covariance for evaluating the next step is therefore \(V_i\approx1.07480R=\begin{pmatrix}1.07480&0.21496\\0.21496&1.07480\end{pmatrix}\).

The correlation parameter remains \(0.2\) in this calculation, but the marginal variances and hence the covariance matrix change with the fitted mean. To form the next GEE, calculate the new residuals and use this refreshed matrix, rather than carrying the initial covariance forward. This makes the numerical iteration correspond to the stated Poisson mean and working-correlation model.

Convergence Checks and Common Pitfalls

A practical convergence check uses at least two quantities: the size of the parameter update and the size of the estimating-equation residual. For example, one may stop only when \(\|\beta^{(t+1)}-\beta^{(t)}\|\) is small relative to the scale of \(\beta\), and a suitably scaled \(\|U(\beta^{(t+1)})\|\) is also small. Tolerances depend on numerical precision and the problem scale; neither criterion alone is universally sufficient. The residual must be evaluated using the refreshed quantities at the candidate estimate.

A small increment can result from a nearly singular update matrix, not from a root. Conversely, a small residual in one parameterization may not indicate that every component of a vector estimating equation is adequately solved. If an update increases the residual sharply or produces invalid fitted means or covariance matrices, use a safer starting value, inspect the computations, or damp the step. The local Newton result does not promise that every full step improves the residual.

  • Confusing Newton with scoring. The exact Jacobian may contain derivatives of the working covariance; \(H=\sum D_i^TV_i^{-1}D_i\) is a specified approximation, not an identity in every GEE.
  • Using stale fitted quantities. After updating \(\beta\), recalculate the means, derivatives, and covariance matrices before evaluating the next equation.
  • Checking only coefficient changes. Verify that the summed equation is close to zero as well as checking the parameter increment.
Key takeaway: A GEE iteration is defined by both its parameter update and how its fitted quantities are evaluated. State whether the step uses the exact Jacobian or a scoring approximation, refresh means and working covariance at each iterate, and assess convergence using both the update and the equation residual.

Check Your Understanding

Use the equations and update logic in this tutorial to answer the following.

  1. What quantities in the summed GEE generally change when \(\beta\) changes?
  2. Why can \(H=\sum_iD_i^TV_i^{-1}D_i\) differ from the exact negative derivative of \(U\)?
  3. In the clustered Poisson example, calculate \(U\) at \(\mu=2.1\) and state whether it is positive or negative.
  4. Why is a small parameter increment alone not enough to establish convergence?
  5. Under what conditions does the local Newton result give quadratic convergence, and why does it not guarantee convergence from every starting value?

Previous

Next

Back to Tutorials
Clinical Biostats