Hand Calculation for Binary Repeated Measures

From a Binary Mean Model to a Cluster Contribution

A generalized estimating equation (GEE) calculation for repeated binary responses connects three pieces of arithmetic: the fitted marginal probabilities, a working covariance matrix for responses within a cluster, and the cluster’s contribution to the estimating equation. Keeping these pieces distinct makes a hand calculation easier to check. The covariance is a working weight; the observed responses enter the contribution through their residuals from the fitted marginal means.

For cluster $i$, let $Y_i=(Y_{i1},\ldots,Y_{im_i})^\mathsf{T}$ be the vector of binary responses, let $\mu_i=E[Y_i\mid X_i]$ be the vector of marginal means, and let $X_i$ be the cluster design matrix. With a logistic mean model, each component is $\mu_{ij}=\operatorname{expit}(x_{ij}^\mathsf{T}\beta)$, where $\operatorname{expit}(t)=1/(1+\exp(-t))$. The coefficient vector $\beta$ describes the marginal mean; the working covariance supplies weights for the residual vector $Y_i-\mu_i$.

Definition: At a specified coefficient vector $\beta$, the cluster estimating-function contribution is $$ U_i(\beta)=D_i^\mathsf{T}V_i^{-1}(Y_i-\mu_i), $$ where $D_i=\partial\mu_i/\partial\beta^\mathsf{T}$ and $V_i$ is the chosen working covariance matrix for cluster $i$. The estimating equation is $\sum_i U_i(\beta)=0$. A contribution is a vector, with one component for each coefficient.

For Bernoulli responses, the variance of response $j$ under the fitted mean is $v_{ij}=\mu_{ij}(1-\mu_{ij})$. Write $A_i=\operatorname{diag}(v_{i1},\ldots,v_{im_i})$. With dispersion fixed at one and an exchangeable working correlation having diagonal entries one and common off-diagonal entry $\alpha$, the working covariance is

V_i=A_i^{1/2}R_i(\alpha)A_i^{1/2}.

Thus its diagonal entries are the Bernoulli variances, and its off-diagonal entry for responses $j$ and $k$ is $\alpha\sqrt{v_{ij}v_{ik}}$. This covariance need not equal the true covariance. It is used to weight the residuals in the estimating function.

Calculate the Fitted Marginal Probabilities

The first calculation is always to form the linear predictor $x_{ij}^\mathsf{T}\beta$ and apply the inverse-logit function. A probability is then converted to a Bernoulli variance by multiplying it by its complement. The design row specifies which coefficient values enter each predictor; the observed response does not enter the fitted probability.

Worked Example: Two Fitted Probabilities

Suppose two visits have covariates $x_1=0$ and $x_2=1$, and the logistic marginal mean model is $\operatorname{logit}(\mu_j)=\beta_0+\beta_1x_j$. At $\beta_0=-1$ and $\beta_1=0.5$, calculate both fitted probabilities and their Bernoulli variances.

At visit 1, the linear predictor is $-1+(0.5)(0)=-1$. Therefore, $\mu_1=1/(1+\exp(1))=0.2689$, and $v_1=\mu_1(1-\mu_1)=0.2689(0.7311)=0.1966$.

At visit 2, the linear predictor is $-1+(0.5)(1)=-0.5$. Therefore, $\mu_2=1/(1+\exp(0.5))=0.3775$, and $v_2=\mu_2(1-\mu_2)=0.3775(0.6225)=0.2350$.

The fitted mean vector is approximately $(0.2689,0.3775)^\mathsf{T}$, and the diagonal of the binary variance matrix is approximately $(0.1966,0.2350)$. These are marginal means and variances at the specified coefficients, not estimates calculated from the two observed outcomes.

Build and Check an Exchangeable Working Covariance

For two responses in a cluster, an exchangeable correlation matrix has ones on the diagonal and the same working correlation $\alpha$ in both off-diagonal positions. The covariance matrix rescales that correlation by the standard deviations implied by the two fitted Bernoulli means. Consequently, its off-diagonal entry is not generally just $\alpha$, nor is it just the square root of the product of the variances.

Worked Example: A Two-Visit Working Covariance

Suppose the fitted probabilities are $\mu_1=0.20$ and $\mu_2=0.40$, and use exchangeable working correlation $\alpha=0.30$. Construct $V$, and calculate its inverse.

The Bernoulli variances are $v_1=0.20(0.80)=0.16$ and $v_2=0.40(0.60)=0.24$. The covariance off the diagonal is $\alpha\sqrt{v_1v_2}=0.30\sqrt{0.16(0.24)}$. Here $0.16(0.24)=0.0384$, so $\sqrt{0.0384}=0.195959$, and the off-diagonal entry is $0.30(0.195959)=0.058788$, to six decimal places.

V= \begin{pmatrix} 0.16 & 0.058788\\ 0.058788 & 0.24 \end{pmatrix}.

For a symmetric matrix with diagonal entries $a,d$ and off-diagonal entry $c$, the determinant is $ad-c^2$, and its inverse, when the determinant is nonzero, is $(ad-c^2)^{-1}\begin{pmatrix}d&-c\\-c&a\end{pmatrix}$. Here the determinant is $0.16(0.24)-(0.058788)^2$, approximately $0.0384-0.003456=0.034944$. Hence

V^{-1}\approx \frac{1}{0.034944} \begin{pmatrix} 0.24 & -0.058788\\ -0.058788 & 0.16 \end{pmatrix} \approx \begin{pmatrix} 6.868 & -1.682\\ -1.682 & 4.579 \end{pmatrix}.

The determinant is positive, and both diagonal entries of $V$ are positive, so this particular two-by-two matrix is invertible and positive definite. This check is useful before using the inverse to weight residuals.

The inverse formula follows directly by multiplication: for $M=\begin{pmatrix}a&c\\c&d\end{pmatrix}$, $$\begin{pmatrix}a&c\\c&d\end{pmatrix}\begin{pmatrix}d&-c\\-c&a\end{pmatrix} =\begin{pmatrix}ad-c^2&0\\0&ad-c^2\end{pmatrix}=(ad-c^2)I.$$ Thus division by $ad-c^2$ gives the inverse whenever $ad-c^2\ne0$. In a hand calculation, retaining the unrounded covariance entries until the inverse is formed reduces rounding error.

Evaluate the Cluster Estimating Function

For a logistic mean model, differentiating the inverse-logit gives $\partial\mu_{ij}/\partial(x_{ij}^\mathsf{T}\beta)=\mu_{ij}(1-\mu_{ij})$. By the chain rule, if $X_i$ has rows $x_{ij}^\mathsf{T}$, then $D_i=A_iX_i$. This derivative matrix is used at the same coefficient vector as the fitted means and covariance. The residual vector is observed minus fitted, $Y_i-\mu_i$, not fitted minus observed.

Worked Example: Contributions From Two Binary Clusters

Consider two clusters, each with two visits and design rows $(1,0)$ and $(1,1)$. Use the coefficients $\beta_0=\log(0.20/0.80)=-1.386294$ and $\beta_1=\log(0.40/0.60)-\log(0.20/0.80)=0.980829$. These give fitted probabilities $\mu=(0.20,0.40)^\mathsf{T}$. Take $\alpha=0.30$, so the covariance and its inverse are those calculated above. Let the observed responses be $Y_1=(1,0)^\mathsf{T}$ and $Y_2=(0,1)^\mathsf{T}$. Calculate both cluster contributions and their sum.

The design matrix and derivative matrix for each cluster are $X=\begin{pmatrix}1&0\\1&1\end{pmatrix}$ and $D=AX=\begin{pmatrix}0.16&0\\0.24&0.24\end{pmatrix}$. For cluster 1 the residual vector is $r_1=(0.80,-0.40)^\mathsf{T}$. Using the unrounded covariance calculation to the precision shown,

V^{-1}r_1\approx \begin{pmatrix}6.868&-1.682\\-1.682&4.579\end{pmatrix} \begin{pmatrix}0.80\\-0.40\end{pmatrix} \approx \begin{pmatrix}6.167\\-3.177\end{pmatrix}.

Multiplying by $D^\mathsf{T}$ gives $U_1=D^\mathsf{T}V^{-1}r_1$. Its intercept component is $0.16(6.167)+0.24(-3.177)=0.224$, and its slope component is $0(6.167)+0.24(-3.177)=-0.763$. Thus $U_1\approx(0.224,-0.763)^\mathsf{T}$.

For cluster 2 the residual vector is $r_2=(-0.20,0.60)^\mathsf{T}$. Then

V^{-1}r_2\approx \begin{pmatrix}6.868&-1.682\\-1.682&4.579\end{pmatrix} \begin{pmatrix}-0.20\\0.60\end{pmatrix} \approx \begin{pmatrix}-2.383\\3.084\end{pmatrix}.

The intercept component is $0.16(-2.383)+0.24(3.084)=0.359$, and the slope component is $0.24(3.084)=0.740$. Therefore $U_2\approx(0.359,0.740)^\mathsf{T}$. The total estimating function at these coefficients is $U_1+U_2\approx(0.583,-0.022)^\mathsf{T}$. The displayed cluster contributions are rounded separately, so their displayed slope components sum to $-0.023$; the coefficients have not been updated, and this total is not expected to be exactly zero at an arbitrary specified coefficient vector.

Why the Cluster Contributions Are Centered

The cluster formula has a useful expectation property. For fixed covariates and fixed working covariance parameters, if the marginal mean model is correct, each cluster contribution has expectation zero. This is a statement about the estimating function, not a claim that the working covariance is the true covariance.

Result: Suppose $E[Y_i\mid X_i]=\mu_i(\beta_0)$, and $D_i$ and $V_i$ are finite matrices determined by the covariates and fixed parameter values, with $V_i$ invertible. Then $E[U_i(\beta_0)\mid X_i]=0$.

To prove the result, condition on $X_i$. At $\beta_0$, $D_i$ and $V_i^{-1}$ are fixed with respect to the response vector, so linearity of conditional expectation gives

E[U_i(\beta_0)\mid X_i] =D_i^\mathsf{T}V_i^{-1}E[Y_i-\mu_i(\beta_0)\mid X_i] =D_i^\mathsf{T}V_i^{-1} \{E[Y_i\mid X_i]-\mu_i(\beta_0)\}=0.

The final equality uses the assumed correct marginal mean. No assumption that $V_i$ equals the true conditional covariance was used. This calculation concerns fixed working covariance parameters; it does not treat an outcome-dependent plug-in covariance as fixed inside the conditional expectation. $\square$

Common Mistakes and Practical Checks

The arithmetic is short, but several small errors can change the direction or scale of a contribution. Use the following checks before interpreting a hand calculation.

Omitting the square root scaling. The exchangeable covariance is $\alpha\sqrt{v_{ij}v_{ik}}$. For variances $0.16$ and $0.24$, the square root is about $0.196$; with $\alpha=0.30$, the covariance is about $0.0588$. The square root itself is not the covariance.
Using the wrong residual order. The stated GEE uses $Y_i-\mu_i$. Reversing the subtraction changes the sign of every contribution.
Mixing up $D_i$ and $D_i^\mathsf{T}$. With $m_i$ responses and $p$ coefficients, $D_i$ is $m_i$-by-$p$, so $D_i^\mathsf{T}V_i^{-1}(Y_i-\mu_i)$ is $p$-by-one, as a coefficient estimating function must be.
Calling a contribution an estimate. A nonzero $U_i(\beta)$ is one cluster’s contribution at the supplied coefficients. Estimating coefficients requires solving the summed estimating equation; this calculation alone does not perform that solution.
Rounding too early. Keep extra digits in fitted probabilities, variances, covariance entries, and the inverse, then round the final contribution to a sensible precision.

Key takeaway: At the specified coefficients, calculate the marginal probabilities first, form the Bernoulli variances and working covariance, compute the observed-minus-fitted residual vector, and apply $D_i^\mathsf{T}V_i^{-1}$. Keeping this order explicit makes each cluster contribution auditable.

Check Your Understanding

Use the definitions and calculation sequence above to check your understanding.

For a fitted binary probability of $0.65$, what is the Bernoulli variance?
If two fitted variances are $0.09$ and $0.25$, and the exchangeable working correlation is $0.20$, what is the off-diagonal working covariance?
Why does the fitted marginal probability depend on the covariates and coefficients, but not on the observed binary response?
For a two-visit cluster with residual vector $(0.3,-0.1)^\mathsf{T}$, what does the inverse working covariance do before multiplication by $D_i^\mathsf{T}$?
Under what mean-model condition does the fixed-working-covariance cluster contribution have conditional expectation zero?