Estimating AR(1) Dependence from Residual Pairs

Why Estimate AR(1) Dependence from Residual Pairs?

In a generalized estimating equation, the working correlation describes how responses from the same subject are related after accounting for the fitted mean and marginal variance. For an AR(1) working correlation, responses one visit apart have correlation \(\alpha\), responses two visits apart have correlation \(\alpha^2\), and so on. Estimating \(\alpha\) therefore means using residual pairs to learn about a pattern across visit lags, rather than treating every pair as if it had the same correlation.

The distinction between visit lag and position among observed responses matters. If a subject misses visit 2 but is observed at visits 1 and 3, that pair has lag 2, not lag 1. Likewise, a lag with many observed pairs can contribute more to a pair-weighted fit than a lag represented by only a few pairs.

Residual Products and the AR(1) Moment Pattern

Definition: For subject \(i\) at scheduled visit \(t\), let \(Y_{it}\) be the response, \(\hat{\mu}_{it}\) its fitted marginal mean, and \(\hat{v}_{it}\) its fitted marginal variance. A standardized residual is \(e_{it}=(Y_{it}-\hat{\mu}_{it})/\sqrt{\hat{v}_{it}}\). For an observed pair at visits \(s\) and \(t\), its standardized residual product is \(e_{is}e_{it}\), and its visit lag is \(k=|t-s|\). Under a discrete-time AR(1) working correlation, the correlation at lag \(k\) is \(\alpha^k\).

The fitted variance in this definition should include any dispersion factor used in the marginal variance model. If the means and variances were known and the standardized errors had mean zero and variance one, their product would have expectation equal to their correlation. In practice, fitted residuals are plug-in quantities: the resulting moments are estimates, not exact identities.

Result: Suppose a fixed set of observed pairs at lag \(k\) is selected independently of the response values. If the true standardized errors have mean zero, unit variance, and correlation \(\alpha^k\) for every pair in that set, then the average of their residual products has expectation \(\alpha^k\).

Proof. Write the true standardized errors for a pair \(j\) at lag \(k\) as \(Z_{j1}\) and \(Z_{j2}\), and let there be \(n_k>0\) such pairs. By the definition of correlation and the unit marginal variances, \(E[Z_{j1}Z_{j2}]=\operatorname{Corr}(Z_{j1},Z_{j2})=\alpha^k\). Linearity of expectation then gives \(E[n_k^{-1}\sum_{j=1}^{n_k}Z_{j1}Z_{j2}] =n_k^{-1}\sum_{j=1}^{n_k}E[Z_{j1}Z_{j2}] =n_k^{-1}n_k\alpha^k=\alpha^k\). No independence between different pairs is needed for this expectation calculation. Replacing true errors with fitted residuals yields a plug-in moment and does not, by itself, preserve this exact unbiasedness. \(\square\)

Pooling by Lag and Fitting the AR(1) Curve

Let \(\mathcal{P}_k\) be the set of observed within-subject pairs whose scheduled visits differ by \(k\), and let \(n_k\) be its size. When \(n_k>0\), pool the products at that lag:

\bar{p}_k=\frac{1}{n_k}\sum_{(i,s,t)\in\mathcal{P}_k} e_{is}e_{it}, \qquad |t-s|=k.

The AR(1) model predicts \(\bar{p}_k\) to be near \(\alpha^k\). One way to estimate \(\alpha\) is to minimize a stated least-squares criterion over a valid parameter range. For the examples below, residuals suggest nonnegative serial association, so we use \(0\leq a<1\), where \(a\) is a candidate value for \(\alpha\). Two choices of criterion are useful to distinguish:

Equal-lag fitting: \(Q_{\mathrm{equal}}(a)=\sum_{k\in\mathcal{K}}(\bar{p}_k-a^k)^2\). Each observed lag summary receives equal weight.
Pair-count-weighted fitting: \(Q_{\mathrm{pair}}(a)=\sum_{k\in\mathcal{K}}n_k(\bar{p}_k-a^k)^2\). A lag with more observed pairs receives more weight.

These are fitting choices, not two names for the same estimator. Equal-lag fitting gives each lag summary equal influence; pair-count weighting gives more influence to a lag represented by more pairs. Pair-count weighting also has a direct least-squares interpretation. For any lag \(k\), if its individual products are \(p_1,\ldots,p_{n_k}\) and their mean is \(\bar{p}_k\), then \(\sum_{j=1}^{n_k}(p_j-a^k)^2 =\sum_{j=1}^{n_k}(p_j-\bar{p}_k)^2+n_k(\bar{p}_k-a^k)^2\). To verify this, write \(p_j-a^k=(p_j-\bar{p}_k)+(\bar{p}_k-a^k)\), square and sum; the cross term is zero because \(\sum_j(p_j-\bar{p}_k)=0\). The first term does not depend on \(a\), so minimizing squared discrepancies across individual products gives the same minimizing \(a\) as \(Q_{\mathrm{pair}}\). This does not mean that pair-count weighting accounts for dependence among products or is automatically optimal.

Worked Example: Fitting Two Pooled Lag Moments

Suppose the four products at lag 1 are \(0.30, 0.50, 0.40,\) and \(0.20\), and the two products at lag 2 are \(0.10\) and \(0.20\). First calculate the pooled moments: \(\bar{p}_1=(0.30+0.50+0.40+0.20)/4=1.40/4=0.35\); \(\bar{p}_2=(0.10+0.20)/2=0.30/2=0.15\).

Using equal-lag fitting gives \(Q_{\mathrm{equal}}(a)=(a-0.35)^2+(a^2-0.15)^2\). Its derivative is \(Q'_{\mathrm{equal}}(a)=2(a-0.35)+4a(a^2-0.15)\). Setting the derivative to zero and dividing by 2 gives \(2a^3+0.7a-0.35=0\). The derivative of the left side is \(6a^2+0.7>0\), so it is strictly increasing; it is negative at 0 and positive at 1. Thus there is exactly one root on \([0,1]\), approximately \(a=0.36316\). At this value, \(2a^3+0.7a-0.35\) is approximately \(0.000003\), confirming the numerical solution. This estimate balances the lag-1 moment \(0.35\) against the lag-2 moment \(0.15\), which the model represents as \(a\) and \(a^2\).

How Unequal Lag Counts Change the Estimate

When pair counts differ, changing the weights can change the fitted parameter, even though the pooled moments themselves have not changed. The following example isolates that effect. For a differentiable fitting criterion, an interior minimum must have derivative zero; checking the derivative and the chosen parameter range is essential before reporting a numerical estimate.

Worked Example: Equal-Lag Versus Pair-Weighted Fitting

Suppose there are eight lag-1 products with mean \(0.40\), and two lag-2 products with mean \(0.10\). Equal-lag fitting uses \(Q_{\mathrm{equal}}(a)=(a-0.40)^2+(a^2-0.10)^2\), whose derivative is \(2(a-0.40)+4a(a^2-0.10)\). At \(a=0.3716\), this derivative is \(-0.0568+4(0.3716)(0.03808656)=-0.0568+0.056611862784=-0.000188137216\). At \(a=0.3717\), it is \(-0.0566+4(0.3717)(0.03816089)=-0.0566+0.056737611252=0.000137611252\). The derivative changes sign between these values, and the minimizing value is approximately \(a=0.37166\).

With pair-count weighting, the criterion is \(Q_{\mathrm{pair}}(a)=8(a-0.40)^2+2(a^2-0.10)^2\). Its derivative is \(16(a-0.40)+8a(a^2-0.10)\). Dividing the stationary equation by 8 gives \(a^3+1.9a-0.8=0\). The left side is strictly increasing for \(a\geq0\), since its derivative is \(3a^2+1.9>0\); its values at 0 and 1 have opposite signs. Hence it has exactly one root in \([0,1]\), approximately \(a=0.389865\). Pair weighting moves the estimate toward the more numerous lag-1 moment of \(0.40\). The lag-2 moment still matters, but it has fewer pairs and therefore less weight under this criterion.

A simpler moment estimate uses only lag 1: if there are observed adjacent-visit pairs, take \(\hat{\alpha}=\bar{p}_1\), subject to any parameter-range convention used for the working correlation. This follows because the AR(1) correlation at lag 1 is \(\alpha\). An all-lag fit uses more of the pattern \(\alpha^k\), but requires choosing how the lag summaries are weighted. In either case, the counts \(n_k\) describe which pairs actually contribute; they should not be replaced by the number of subjects or by a presumed complete schedule.

Missing Visits and the Meaning of Lag

The scheduled visit labels determine the lag, not the order of the recorded values after omissions. This is especially important when visit intervals are equally spaced and a response is missing at an intermediate visit. Keeping the original labels preserves the intended discrete-time AR(1) pattern: a gap of two scheduled intervals corresponds to \(\alpha^2\), not \(\alpha\).

Worked Example: A Missed Visit Changes Which Moment a Pair Enters

Subject A is observed at visits 1, 2, and 3, with standardized residuals \(0.2, 0.4,\) and \(0.5\), respectively. Subject B is observed at visits 1, 3, and 4; visit 2 is missing, and the residuals at visits 1, 3, and 4 are \(0.3, 0.6,\) and \(0.5\).

For subject A, the lag-1 products are \(0.2(0.4)=0.08\) and \(0.4(0.5)=0.20\); the lag-2 product is \(0.2(0.5)=0.10\). For subject B, the pair at visits 3 and 4 contributes a lag-1 product \(0.6(0.5)=0.30\); the pair at visits 1 and 3 contributes a lag-2 product \(0.3(0.6)=0.18\); and the pair at visits 1 and 4 contributes a lag-3 product \(0.3(0.5)=0.15\).

Pooling by actual scheduled lag gives \(n_1=3\) and \(\bar{p}_1=(0.08+0.20+0.30)/3=0.58/3\approx0.1933\); \(n_2=2\) and \(\bar{p}_2=(0.10+0.18)/2=0.14\); and \(n_3=1\) and \(\bar{p}_3=0.15\). The lag-1 moment estimate is therefore approximately \(0.1933\). If subject B’s three observed responses were incorrectly renumbered as visits 1, 2, and 3, the product \(0.3(0.6)=0.18\) would be misclassified as lag 1, and the product \(0.3(0.5)=0.15\) would be misclassified as lag 2. That would change both the pooled lag summaries and any AR(1) fit based on them.

Assumptions, Interpretation, and a Common Mistake

The moment interpretation is most straightforward when the scheduled visit gaps are meaningful, the marginal mean and variance model are adequate, and the observed pair set is not selected according to the residual outcomes in a way that distorts their products. If missingness depends on unobserved responses, the observed pooled products need not represent the corresponding population correlations. Even with noninformative missingness, estimated means and variances make the residual products approximate moment information. Also, products from pairs in the same subject can be dependent, so the number of pairs is not automatically the number of independent observations.

A common mistake is to treat consecutive recorded responses as consecutive scheduled visits. The mistaken reasoning is: “These are the next two observed measurements for this subject, so their product estimates the lag-1 correlation.” The correction is to use the scheduled visit labels: if the measurements occurred at visits 1 and 3, their lag is 2, so their AR(1) target is \(\alpha^2\). Another mistake is to report a fitted value without saying whether lag summaries were weighted equally or by their pair counts. Those criteria can produce different answers, as the calculations above show.

Key takeaway: Pool standardized residual products using actual within-subject visit lags. Then state how the lag-specific moments are combined: a lag-1 moment estimate, equal-lag fitting, and pair-count-weighted fitting are distinct procedures. Unequal counts and missing visits affect the calculation through the pairs that are observed at each true lag.

Check Your Understanding

Use the definitions and calculations above to answer these questions.

Under an AR(1) working correlation, what correlation does the model assign to a pair separated by three scheduled visits?
There are six observed products at lag 1 with mean \(0.25\), and two at lag 2 with mean \(0.10\). Write the equal-lag and pair-count-weighted fitting criteria.
A subject is observed at visits 2 and 5, with no observations at visits 3 or 4. What is the pair’s lag, and what AR(1) correlation does the model assign it?
Why does the expectation argument for pooled residual products not require independence between different pairs?
Give one reason the observed pooled product at a lag may fail to estimate the population correlation at that lag reliably.