Autocorrelation and Mixing in MCMC Chains

Why Autocorrelation Matters

An MCMC output is a sequence of draws, not simply a collection of independent values. In a Bayesian NLME analysis, successive draws of a population parameter such as typical clearance, or a variance parameter such as an interindividual variance, may be similar because the chain moves gradually through the posterior distribution. When this happens, a long-looking chain can contain less information about a posterior mean than an equally long independent sample. Autocorrelation measures this serial dependence; lag plots make some of its patterns visible.

Autocorrelation is not itself a verdict that a chain has converged or failed. It is a property of the ordered draws, and its interpretation depends on whether the chain is sampling a stable target distribution. The useful question is whether dependence is affecting exploration and the precision of the summaries, and whether other diagnostics tell a consistent story.

Definitions and Assumptions

Definition: Let $\theta_1,\ldots,\theta_M$ be $M$ post-warmup draws of one scalar parameter from a chain, with sample mean $\bar{\theta}$. For lag $k=0,1,\ldots,M-1$, define $$ \widehat{\gamma}_k = \frac{1}{M}\sum_{t=1}^{M-k} (\theta_t-\bar{\theta})(\theta_{t+k}-\bar{\theta}), \qquad \widehat{\rho}_k = \frac{\widehat{\gamma}_k}{\widehat{\gamma}_0}. $$ Here $\widehat{\rho}_k$ is the sample autocorrelation at lag $k$; lag $k$ compares draws $k$ iterations apart. The common denominator $M$ is used at every lag, so this is one conventional finite-chain estimator, not the only possible estimator. It is defined only for $0\leq k<M$, and requires $\widehat{\gamma}_0>0$.

At lag zero, $\widehat{\rho}_0=1$. A positive value at lag one indicates that adjacent draws tend to be on the same side of their mean together; a negative value indicates a tendency for adjacent deviations to have opposite signs. A lag-$k$ autocorrelation describes linear association at that separation, not every possible kind of dependence.

Definition: For a covariance-stationary chain with finite variance $\sigma_\theta^2$, the population autocorrelation at lag $k$ is $\rho_k=\operatorname{Cov}(\theta_t,\theta_{t+k})/\sigma_\theta^2$. Covariance stationarity means that the mean and variance do not change with $t$, and that the covariance depends on the separation $k$, not on the particular iteration.

Interpreting the sample autocorrelations as estimates of chain behavior presumes that post-warmup draws are approximately from a stable, stationary regime. The chain should also be exploring the intended posterior rather than remaining trapped in an unvisited region. For long-run variance and effective-sample-size calculations, additional regularity is needed, commonly including sufficiently decaying autocorrelations so their sum is meaningful. If the chain is still drifting, an autocorrelation computed around one overall sample mean can be misleading.

Reading Lag Plots

A lag-$k$ plot places $\theta_t$ on the horizontal axis and $\theta_{t+k}$ on the vertical axis for all available pairs. A narrow upward-sloping band suggests positive serial association; a downward-sloping band suggests negative association. A diffuse cloud with little linear pattern is consistent with weak autocorrelation at that lag, but does not establish independence. Curvature, clusters, or distinct bands can indicate structure that a single correlation coefficient does not summarize well.

For a variance parameter, the parameter is positive and its posterior distribution may be skewed. A lag plot on the original scale and one on a log scale can therefore emphasize different features. Diagnose the parameter scale relevant to the intended summary, and consider a transformation when it helps reveal structure; do not assume that autocorrelation is unchanged by nonlinear transformations.

Worked Example: Calculating Sample Autocorrelations

Suppose a hypothetical post-warmup chain for a population parameter has five consecutive draws $(1,2,3,4,5)$ — five draws, not six. To avoid ambiguity, the data are $\theta_1=1,\theta_2=2,\theta_3=3,\theta_4=4,\theta_5=5$, so $M=5$. Calculate the sample autocorrelations at lags one and two.

The mean is $\bar{\theta}=(1+2+3+4+5)/5=3$, and the deviations are $(-2,-1,0,1,2)$. Thus $\widehat{\gamma}_0=[4+1+0+1+4]/5=2$.

At lag one, the paired products are $2,0,0,2$, so $\widehat{\gamma}_1=(2+0+0+2)/5=0.8$ and $\widehat{\rho}_1=0.8/2=0.4$.

At lag two, the paired products are $0,-1,0$, so $\widehat{\gamma}_2=-1/5=-0.2$ and $\widehat{\rho}_2=-0.2/2=-0.1$. This tiny, deliberately simple sequence illustrates the arithmetic only: five draws are far too few to diagnose a real chain reliably. The estimated negative value at lag two should not be treated as evidence of a stable negative-dependence pattern.

How Dependence Changes Information

The connection to information can be derived by calculating the variance of the chain average. Suppose the chain is covariance-stationary, has mean $\mu$, variance $\sigma_\theta^2$, and autocorrelation $\rho_k$. The sample average is $\bar{\theta}=M^{-1}\sum_{t=1}^M\theta_t$. Expanding its variance gives the exact finite-$M$ expression:

\begin{aligned} \operatorname{Var}(\bar{\theta}) &= \frac{1}{M^2}\sum_{t=1}^M\sum_{s=1}^M \operatorname{Cov}(\theta_t,\theta_s)\\ &= \frac{\sigma_\theta^2}{M} \left[1+2\sum_{k=1}^{M-1} \left(1-\frac{k}{M}\right)\rho_k\right]. \end{aligned}

For each separation $k$, there are $M-k$ pairs above the diagonal and $M-k$ matching pairs below it. Their covariances are each $\sigma_\theta^2\rho_k$; the $M$ diagonal terms each contribute $\sigma_\theta^2$. Dividing this sum by $M^2$ yields the displayed result. If the draws were independent, every $\rho_k$ for $k>0$ would be zero and the variance would be $\sigma_\theta^2/M$. Positive serial dependence generally increases the variance of the average, whereas negative dependence can decrease it. The net effect depends on the weighted sum, not on the sign of just one lag.

When the chain is sufficiently well behaved that the autocorrelation series is summable, the large-$M$ variance is commonly expressed using the integrated autocorrelation time, $\tau=1+2\sum_{k=1}^{\infty}\rho_k$. The corresponding effective sample size for estimating a mean is approximately $M/\tau$. This is a variance-equivalence measure: it compares the chain average's precision with that of an independent sample from the same target distribution. It is not a count of distinct posterior states, and it does not make a poorly explored chain reliable.

Worked Example: Positive Autocorrelation in a Clearance Chain

Imagine a hypothetical chain with $M=900$ post-warmup draws of typical clearance, measured in L/h. Suppose its autocorrelations are well described by $\rho_k=0.8^k$, as in a stationary first-order autoregressive dependence pattern. This is an illustrative assumption, not a claim about a real analysis.

The integrated autocorrelation time is $\tau=1+2\sum_{k=1}^{\infty}0.8^k$. The geometric series sums to $0.8/(1-0.8)=4$, so $\tau=1+2(4)=9$. The approximate effective sample size is $900/9=100$ independent-draw equivalents for the chain average.

Suppose the posterior standard deviation of clearance is $2.4$ L/h. The independent-draw approximation to the Monte Carlo standard error of its mean would be $2.4/\sqrt{900}=0.0800$ L/h. Accounting for the dependence gives $2.4/\sqrt{100}=0.240$ L/h. Thus, under the stated dependence model, the chain average has about three times the Monte Carlo standard error of the independence calculation. This concerns numerical precision of the estimated posterior mean, not the posterior uncertainty in clearance itself.

Worked Example: Negative Autocorrelation and Effective Information

Consider a second hypothetical chain with $M=1{,}000$ draws of a population parameter, with $\rho_k=(-0.4)^k$. This alternating pattern has a summable autocorrelation series. Its integrated autocorrelation time is $\tau=1+2\sum_{k=1}^{\infty}(-0.4)^k$.

The geometric sum is $-0.4/[1-(-0.4)]=-0.4/1.4=-0.2857$, so $\tau=1+2(-0.2857)=0.4286$. The approximate mean ESS is $1{,}000/0.4286\approx2{,}333$. It exceeds the number of recorded draws because the negative serial dependence makes the chain average less variable than an independent average of 1,000 draws with the same marginal variance.

This result is not a contradiction: autocorrelation changes the information in an average, and negative dependence can increase it. An ESS larger than $M$ does not mean the chain generated more than 1,000 distinct values, nor does it by itself establish that the target distribution was explored correctly. It is an approximation tied to the mean and the chain's dependence pattern.

Using the Diagnostics Together

For population and variance parameters, inspect the chain on a scale that is interpretable and check whether dependence persists across several lags. A slow, smooth traversal in a trace plot often corresponds to positive autocorrelation over many lags. A lag plot can expose a tight diagonal band, alternating pattern, clusters, or nonlinear structure. The estimated autocorrelation function (ACF), a plot of $\widehat{\rho}_k$ against lag, summarizes linear dependence across a range of separations.

There is no universal lag at which an ACF must become exactly zero. At large lags, fewer pairs are available and sample autocorrelations are noisy; isolated distant spikes should not be overinterpreted. Nor does a small ACF at a few lags prove adequate mixing: dependence can be nonlinear, and separate regions of posterior mass may not have been visited. Compare these findings with the trace plots, rank-normalized Rhat, effective sample size, and Monte Carlo standard errors. Each addresses a different aspect of chain behavior or estimation precision.

Thinning keeps only every $q$th draw. It can reduce the number of stored draws, but does not automatically improve the information per unit of computation: the retained sequence may be less correlated while containing fewer draws. Assess the retained chain's effective information and Monte Carlo error rather than treating a smaller lag-one autocorrelation as proof that thinning helped.

Worked Example: A Lag Plot for a Variance Parameter

Suppose a hypothetical NLME analysis produces consecutive post-warmup draws of an interindividual variance parameter (in squared parameter units): $0.20, 0.22, 0.21, 0.24, 0.23, 0.22$. The lag-one plot uses the five pairs $(0.20,0.22)$, $(0.22,0.21)$, $(0.21,0.24)$, $(0.24,0.23)$, and $(0.23,0.22)$.

These points lie in a compact region and most successive values are close, suggesting short-step movement on the original scale. That visual pattern is compatible with positive lag-one dependence, but six draws cannot establish its strength or persistence. The appropriate follow-up is to examine a much longer post-warmup sequence, its trace plot and ACF, and the diagnostics from other chains. Because variance parameters are positive and may be skewed, inspecting log-variance draws as well can reveal whether apparent clustering is mainly a scale effect. Neither the lag plot nor a transformation replaces checks of convergence and exploration.

Common Mistakes and Limitations

Equating autocorrelation with failure. MCMC draws are ordinarily dependent. The issue is how dependence affects exploration and the precision of the desired summary, not whether every autocorrelation is zero.
Claiming dependence always reduces information. Positive serial dependence generally inflates the variance of a chain average; negative dependence can reduce it. Consider the combined autocorrelation pattern.
Treating ACF as a convergence test. ACF describes within-chain serial structure. It cannot alone show that chains reached the same posterior distribution or visited all important regions.
Trusting a noisy tail of the ACF. High-lag estimates use fewer pairs, so apparent spikes may be sampling noise. Look for sustained patterns and use multiple diagnostics.
Reading a lag plot as a complete dependence summary. Correlation captures linear association, while clusters or curved patterns can be missed by a single coefficient.
Using stationarity formulas during drift. If the chain has not settled into a stable regime, a global mean and ACF can conceal movement. Investigate the trace plot and chain agreement before interpreting long-run quantities.

Autocorrelation diagnostics matter because a large number of stored draws can create a false impression of precision. They help connect visible chain behavior to the variability of estimates such as posterior means, while leaving the broader assessment of convergence and model adequacy to the full set of diagnostics.

Check Your Understanding

Try each question, then compare your reasoning with the concise solutions below.

For a chain of $M=40$ draws, what is the largest lag for which the sample autocorrelation formula in this tutorial uses at least one available pair?
Suppose the integrated autocorrelation time is $4$ for a chain of $800$ draws. What is its approximate mean ESS, and what does that quantity represent?
Why is it incorrect to say that all autocorrelation reduces information in a chain average?
A lag plot shows a narrow upward band, but the trace plot shows a long-term drift. Should the ACF be interpreted as a reliable summary of stationary mixing? Why?
Why can thinning lower the apparent lag-one autocorrelation without increasing effective information per unit of computation?

Solutions. (1) The largest lag is $M-1=39$, where there is one available pair. (2) The approximate mean ESS is $800/4=200$ independent-draw equivalents for estimating a mean; it is not a count of distinct states. (3) Positive dependence generally increases the variance of the average, but negative dependence can decrease it. (4) No. The drift calls the stable-regime assumption into question, so a global ACF may not represent stationary mixing; investigate chain stability and other diagnostics first. (5) Keeping fewer draws can reduce correlation among retained values but also discards draws. The retained chain's effective information depends on both its length and its dependence, not on lag-one autocorrelation alone.