Autocorrelation and Slow-Mixing Random Effects

Why Autocorrelation Matters in Joint-Model Chains

An MCMC chain is useful when its draws explore the posterior distribution well enough to estimate quantities of interest. Consecutive draws are often dependent: if a parameter is currently high, the next draw may also be high. This persistence is autocorrelation. When it is strong, a long chain can contain much less information about a posterior mean than the same number of nearly independent draws.

This issue can be especially visible for subject-specific random effects. In a joint model, a subject’s random effect can be informed by that subject’s longitudinal measurements and event history, as well as by population-level parameters. Some subjects have sparse or weakly informative records, so their effects may be difficult to explore efficiently. The number of stored draws alone does not show whether those effects mix well.

Here, \(x^{(m)}\) denotes the value of a sampled scalar quantity at iteration \(m\) within one chain. It may represent a population parameter or a subject-specific effect such as \(b_i\), the random effect for subject \(i\). Autocorrelation is calculated within a chain; it is not the correlation between different subjects’ effects or the correlation between repeated biomarker measurements.

Definition: For a stationary chain with finite variance, the lag-\(k\) autocorrelation of \(x\) is \(\rho_k=\operatorname{Corr}(x^{(m)},x^{(m+k)})\). It describes the linear dependence between draws \(k\) iterations apart. In a finite chain, a common estimate is

\(\hat{\rho}_k=\dfrac{\sum_{m=1}^{n-k}(x^{(m)}-\bar{x})(x^{(m+k)}-\bar{x})}{\sum_{m=1}^{n}(x^{(m)}-\bar{x})^2}\),

where \(n\) is the number of retained draws from that chain and \(\bar{x}\) is their mean. Estimates at large lags are noisy, so the pattern across early lags is usually more informative than any single late-lag estimate.

From Persistent Draws to Less Monte Carlo Information

The connection between autocorrelation and Monte Carlo precision can be derived directly. Suppose a stationary chain has variance \(\sigma_x^2\), and let \(\bar{x}\) be the mean of \(n\) draws. Since the variance of a sum includes the covariances between every pair of draws,

\operatorname{Var}(\bar{x}) =\frac{\sigma_x^2}{n} \left[1+2\sum_{k=1}^{n-1}\left(1-\frac{k}{n}\right)\rho_k\right].

For a sufficiently long chain with a summable autocorrelation sequence, the factor in brackets is approximately \(1+2\sum_{k=1}^{\infty}\rho_k\). This factor is the integrated autocorrelation time, often written \(\tau\). The effective sample size is approximately \(n/\tau\): it is the size of an independent sample that would yield roughly the same Monte Carlo variance for the mean. It does not change the number of draws actually generated.

Key result: If autocorrelation decreases geometrically as in a stationary AR(1) process, with \(\rho_k=\rho^k\) and \(0\leq\rho<1\), then \(\tau=(1+\rho)/(1-\rho)\). Thus \(N_{\mathrm{eff}}\approx n(1-\rho)/(1+\rho)\). This is a useful illustration and rough diagnostic, not a replacement for the autocorrelation-based effective sample size calculation used by modern MCMC software.

To see the result, substitute \(\rho_k=\rho^k\) into the long-chain approximation. The geometric series gives \(\sum_{k=1}^{\infty}\rho^k=\rho/(1-\rho)\), so \(\tau=1+2\rho/(1-\rho)=(1+\rho)/(1-\rho)\). Equating the approximate variance \(\sigma_x^2\tau/n\) with the variance \(\sigma_x^2/N_{\mathrm{eff}}\) of an independent sample gives \(N_{\mathrm{eff}}\approx n/\tau\).

Autocorrelation is not the same as bias, and a low autocorrelation does not establish convergence. A chain can move rapidly within the wrong region, or different chains can explore different regions. Use autocorrelation alongside trace plots, between-chain agreement, and effective sample size; the earlier material on those diagnostics gives their complementary roles.

Worked Example: Interpreting a High Lag-One Autocorrelation

Suppose a hypothetical joint-model chain retains \(n=1{,}800\) draws for a log-hazard association parameter. Its estimated lag-one autocorrelation is about \(0.80\), and the autocorrelation pattern appears roughly geometric. Use the AR(1) approximation to assess the information in the chain.

The approximate integrated autocorrelation time is \(\tau=(1+0.80)/(1-0.80)=1.80/0.20=9.00\). Therefore \(N_{\mathrm{eff}}\approx1{,}800/9.00=200\) effective draws. If the posterior standard deviation of the parameter is \(0.60\) per unit increase in the latent association scale, the approximate Monte Carlo standard error of its posterior mean is \(0.60/\sqrt{200}=0.0424\) on that scale.

The chain contains 1,800 retained draws but, under this approximation, the mean has Monte Carlo precision similar to that from about 200 independent draws. This does not mean that the posterior standard deviation should be divided by \(\sqrt{1{,}800}\) when assessing the simulation error of the mean; the dependence reduces the effective information.

Inspecting Subject Effects as Well as Population Parameters

A single summary for a model is not enough. Inspect autocorrelation for parameters central to the scientific question and for selected subject effects. In a shared-random-effects joint model, useful targets may include a population association parameter, a variance parameter, and \(b_i\) for subjects with notably sparse longitudinal or event information. The posterior scales and mixing behavior can differ substantially across these quantities.

For each target, examine the chain trace and autocorrelation by chain. A slowly decaying autocorrelation pattern indicates that nearby draws remain similar over many iterations. A trace that drifts, sticks for long stretches, or shows abrupt changes may indicate a more serious problem than ordinary persistence. Comparing chains separately matters: pooling their draws before inspecting dependence can conceal differences in where chains have spent their time.

Worked Example: One Subject Effect Mixes More Slowly

Consider a hypothetical fit with 2,000 retained draws for each of two subject-specific random intercepts. Subject 17 has sparse biomarker measurements and a censored event time; subject 42 has more longitudinal measurements. Suppose their lag-one autocorrelations are approximately \(0.95\) and \(0.35\), respectively, and both patterns are roughly geometric. Their posterior standard deviations are \(0.80\) and \(0.50\) log-biomarker units.

For subject 17, \(\tau\approx(1+0.95)/(1-0.95)=1.95/0.05=39.0\), so \(N_{\mathrm{eff}}\approx2{,}000/39.0=51.3\). Its posterior-mean Monte Carlo standard error is approximately \(0.80/\sqrt{51.3}=0.112\) log-biomarker units.

For subject 42, \(\tau\approx(1+0.35)/(1-0.35)=1.35/0.65=2.08\), so \(N_{\mathrm{eff}}\approx2{,}000/2.08=962\). Its corresponding Monte Carlo standard error is approximately \(0.50/\sqrt{962}=0.0161\) log-biomarker units.

These are illustrative AR(1)-based approximations, not exact software estimates. They show why equal draw counts do not imply equal precision: the sparse subject’s effect is much more persistent and has a substantially larger estimated simulation error. The different posterior standard deviations also reflect posterior uncertainty, which should not be confused with Monte Carlo error.

Assessing Practical Remedies

A remedy should be judged by whether it improves exploration and the precision of the quantities that matter, while preserving the intended posterior target. Use comparable runs and chain summaries: retain the same post-warmup draw count where possible, examine traces and autocorrelation by chain, and compare effective sample sizes and Monte Carlo standard errors. Posterior summaries should also be compatible across the runs. An apparently improved autocorrelation estimate is not sufficient if chains disagree or traces show nonstationarity.

For random-effects models, centered and noncentered parameterizations can have different mixing behavior; which works better depends on the model and data. Blocking parameters or marginalizing random effects may also help in appropriate settings. These are practical options, not guarantees. Change one feature at a time when feasible, then rerun the diagnostics. Do not assume that a particular parameterization is universally superior.

Worked Example: Comparing Two Parameterizations by Chain Summaries

In a hypothetical joint-model analysis, four chains each retain 1,000 draws for an association parameter. Under one parameterization, the lag-one autocorrelation is about \(0.92\); after switching parameterization and rerunning the same model, it is about \(0.45\). Suppose the autocorrelation patterns are roughly geometric, the chains in both runs agree on a posterior mean of about \(-0.72\), and the posterior standard deviation is \(0.30\).

There are \(n=4{,}000\) retained draws across the chains. For the first run, the AR(1) proxy gives \(N_{\mathrm{eff}}\approx4{,}000(1-0.92)/(1+0.92)=4{,}000(0.08/1.92)=167\). The approximate Monte Carlo standard error of the mean is \(0.30/\sqrt{167}=0.0232\).

For the second run, \(N_{\mathrm{eff}}\approx4{,}000(1-0.45)/(1+0.45)=4{,}000(0.55/1.45)=1{,}517\), and the approximate Monte Carlo standard error is \(0.30/\sqrt{1{,}517}=0.00770\). In this example, the second run has lower estimated persistence and a smaller estimated simulation error, while the posterior mean is similar. That is encouraging evidence of improved efficiency, but it does not alone prove that either run converged or that the model is appropriate.

Failure Modes, Alternatives, and Limits

High autocorrelation is a diagnostic clue, not a diagnosis of its own cause. Dependence can arise from the posterior structure, the parameterization, or the sampling algorithm. A very high lag-one estimate can also be unstable in a short chain. Interpret it with the trace and the full early-lag pattern rather than treating a threshold as a pass-or-fail rule.

Extending a chain can increase effective sample size if the chain is already exploring the posterior adequately. It does not repair chains that have not reached a stable target or that remain trapped in a region. Thinning by keeping only every \(q\)th draw may lower the observed correlation among saved draws, but it also discards draws and generally does not create more information per computation. It is mainly a storage convenience, not a default solution to slow mixing.

Autocorrelation-based summaries have limits. They concern Monte Carlo error for specific quantities, and good precision for a population mean does not guarantee adequate exploration of tails or rare regions. Subject effects may require separate attention because a few poorly mixing effects can be hidden by summaries over many subjects. Conversely, a high autocorrelation does not by itself show that posterior estimates are unusable: the relevant question is whether the chains adequately explore the target and provide sufficient effective information for the intended summaries.

Mistake: “The chain has 10,000 draws, so its posterior mean must be precise.”
Correction: Precision depends on dependence as well as draw count. Examine effective sample size and Monte Carlo error for the parameter or subject effect being reported.
Mistake: “The lag-one autocorrelation fell after thinning, so the chain is fixed.”
Correction: Thinning can reduce the correlation among stored draws while discarding information. Check exploration and effective sample size per computation instead.
Mistake: “A low autocorrelation proves convergence.”
Correction: Autocorrelation describes within-chain dependence. It does not establish that chains reached or adequately explored the posterior.

Key takeaway: Diagnose slow mixing at the level of individual sampled quantities, including subject effects. Read autocorrelation together with traces and between-chain agreement, and assess remedies through effective sample size and Monte Carlo error rather than raw draw count alone.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

A chain has 3,000 retained draws and an approximate AR(1) lag-one autocorrelation of \(0.50\). Find the approximate effective sample size.
Why can a subject-specific random effect have a much smaller effective sample size than a population parameter in the same fit?
Does an autocorrelation estimate near zero establish that a chain has converged? Explain.
A researcher thins a chain and observes lower lag-one autocorrelation among saved draws. What additional comparison is needed before claiming the change improved sampling?

Solutions. (1) \(N_{\mathrm{eff}}\approx3{,}000(1-0.50)/(1+0.50)=3{,}000(0.50/1.50)=1{,}000\). This is the AR(1) approximation, not necessarily the software’s final effective sample size. (2) The subject’s data may provide less information about its effect, or its draws may be more strongly dependent on other sampled quantities. Its own trace and autocorrelation should therefore be inspected rather than inferred from population-level summaries. (3) No. Low within-chain autocorrelation does not show that the chain reached the posterior or that separate chains agree; inspect traces and between-chain diagnostics as well. (4) Compare effective information and Monte Carlo error for the quantities of interest, accounting for the discarded draws and computational cost, and check that the chains explore consistently. A lower saved-draw autocorrelation alone is not evidence of increased sampling efficiency.