MCMC Targets and Posterior Draws in Bayesian GLMMs

From a Bayesian GLMM to a Posterior Target

A Bayesian GLMM describes the data with a likelihood and the unknown quantities with a prior distribution. Computation by Markov chain Monte Carlo (MCMC) does not change that model: its purpose is to generate draws that can approximate summaries of the model’s posterior distribution. The central question is therefore not simply “What numbers did the algorithm produce?” but which posterior distribution are those numbers intended to represent?

Write $y_{ij}$ for the response from observation $i$ in group $j$, $x_{ij}$ for its fixed-effect design vector, and $z_{ij}$ for its random-effect design vector. A GLMM specifies a conditional response distribution given the random effects, with a linear predictor such as $\eta_{ij}=x_{ij}^{T}\beta+z_{ij}^{T}b_j$, linked to the conditional mean by a link function. The random effects $b_j$ have a distribution governed by variance or covariance parameters. In Bayesian analysis, prior distributions are also assigned to the unknown parameters.

Definition: The joint posterior target is the posterior distribution of all unknown quantities represented together, such as fixed effects $\beta$, random effects $b=(b_1,\ldots,b_J)$, and variance or covariance parameters $D$. Its density is proportional to the likelihood times the random-effects density and the parameter prior: $$ p(\beta,b,D\mid y)\ \propto\ p(y\mid \beta,b)\,p(b\mid D)\,p(\beta,D). $$ A marginal posterior, such as $p(\beta\mid y)$, is obtained by integrating the joint posterior over the other unknown quantities.

This expression assumes that the model’s likelihood, random-effects distribution, and priors are specified, and that their product can be normalized to a proper posterior density. If the model includes additional unknown quantities, they belong in the target and in the prior factorization. A sampler may represent the joint target, or a marginal target if other quantities have been integrated out. A posterior summary is meaningful only when its target is stated clearly.

Worked Example: Writing the Joint Target for a Logistic Random-Intercept GLMM

Suppose a hypothetical data set has three groups, each with five binary outcomes. The numbers of successes are $r_1=3$, $r_2=1$, and $r_3=4$. In a logistic random-intercept model, let $b_j$ be the intercept deviation for group $j$, and let $\operatorname{logit}^{-1}(u)=1/(1+e^{-u})$. Use the priors $\beta\sim N(0,2.5^2)$, $b_j\mid\sigma^2\sim N(0,\sigma^2)$, and a half-normal prior with scale 1 for $\sigma>0$.

Conditional on $\beta$ and $b_j$, the success probability in group $j$ is $q_j=\operatorname{logit}^{-1}(\beta+b_j)$. The binomial likelihood is $$ p(y\mid\beta,b)=\prod_{j=1}^{3}{5\choose r_j}q_j^{r_j}(1-q_j)^{5-r_j}. $$

The random-effects density contributes $\prod_{j=1}^{3}\phi(b_j;0,\sigma^2)$, where $\phi(\,\cdot\,;0,\sigma^2)$ is a normal density, and the priors contribute $p(\beta)p(\sigma)$. Thus the joint posterior target is $$ p(\beta,b_1,b_2,b_3,\sigma\mid y)\ \propto\ \prod_{j=1}^{3}{5\choose r_j}q_j^{r_j}(1-q_j)^{5-r_j} \prod_{j=1}^{3}\phi(b_j;0,\sigma^2)\,p(\beta)p(\sigma). $$

The target includes the group effects as well as $\beta$ and $\sigma$. For instance, a posterior summary of $\sigma$ based on this joint target averages over uncertainty about the $b_j$ values rather than treating their unknown values as fixed and known. A marginal target for $\sigma$ would integrate $\beta$ and all three group effects out of the joint posterior.

Why Markov Chain Draws Can Approximate Posterior Summaries

An MCMC algorithm produces a sequence of states $X_1,X_2,\ldots$, where a state may contain all unknown quantities in the joint target. The Markov property means that the distribution of the next state, given the chain’s history, depends on the current state. Consecutive states are generally dependent. That dependence is not itself an error: the chain is constructed so that its intended stationary distribution is the posterior target.

Key result: If a Markov chain has the posterior target $\pi$ as its invariant distribution and satisfies suitable ergodicity conditions (including irreducibility, aperiodicity, and positive recurrence in standard settings), then, for a function $h$ with a finite posterior expectation, $$ \frac{1}{M}\sum_{m=1}^{M}h(X_m)\ \longrightarrow\ E_{\pi}\{h(X)\} $$ as $M$ increases, under the usual Markov chain law of large numbers. This result justifies estimating posterior expectations by averages of chain draws; it does not make the draws independent.

For example, setting $h(X)=\beta_k$ estimates a posterior mean coefficient, while setting $h(X)=I(\beta_k>0)$ estimates the posterior probability that the coefficient is positive. Empirical quantiles of the sampled values estimate posterior quantiles. These are numerical approximations to features of the posterior distribution, not replacements for the model or its likelihood.

A chain started at an arbitrary state need not initially behave as if it were sampling from its stationary distribution. In practice, analysts assess whether the chain has explored the relevant posterior region and whether its summaries are stable. A trace plot, which plots sampled values against iteration, can reveal prolonged sticking, drifting, or abrupt changes. These are warning signs, not a proof that a chain has or has not converged. Details of formal diagnostics and sampling algorithms require separate treatment.

Posterior Uncertainty and Monte Carlo Error

Two kinds of uncertainty must not be confused. Posterior uncertainty describes uncertainty about a parameter conditional on the observed data and model; the posterior standard deviation is one measure of it. Monte Carlo error is numerical error from estimating a posterior quantity with a finite, dependent chain. A posterior may be broad even when its mean is estimated precisely, or narrow while its estimated mean is numerically unstable.

If $h(X_m)$ has posterior variance $\gamma_0$ and lag-$k$ autocovariance $\gamma_k$, then under standard conditions the long-run variance of its sample average depends on $$ \gamma_0+2\sum_{k=1}^{\infty}\gamma_k. $$ Consequently, the independent-draw standard error $\sqrt{\gamma_0/M}$ generally does not apply when the draws are autocorrelated. Positive autocorrelation typically makes the actual numerical uncertainty larger than that naive calculation suggests. With a sufficiently long, well-behaved chain, Monte Carlo error can be estimated using methods that account for dependence; a small chain alone does not support a reliable estimate.

Worked Example: A Chain Average and the Effect of Dependence

Consider eight illustrative draws of a fixed-effect coefficient, in log-odds units: $0.5,0.5,0.5,0.5,1.5,1.5,1.5,1.5$. These numbers are hypothetical MCMC output, not data from a real study. The draw average is $$ \bar{\beta}=\frac{4(0.5)+4(1.5)}{8}=\frac{8}{8}=1.0. $$

The sample variance of these eight draws is $$ s^2=\frac{4(0.5-1)^2+4(1.5-1)^2}{8-1} =\frac{2}{7}=0.2857, $$ so their sample standard deviation is $s=\sqrt{2/7}\approx0.5345$ log-odds units. If they were independent, a naive estimated standard error for the draw average would be $$ \frac{s}{\sqrt{8}}=\sqrt{\frac{2}{56}}\approx0.1890 $$ log-odds units.

But these draws occur in long runs of the same value. Their strong serial pattern means that the independent-draw calculation is not a justified estimate of the average’s Monte Carlo error. In this short sequence, there is also too little information for a dependable dependence-adjusted error estimate. The value $0.1890$ is therefore the result of an explicitly hypothetical independence calculation—not a trustworthy precision statement for these chain draws. The draw average may still be a useful illustration of the estimator, but a real analysis needs enough well-exploring draws to support a numerical precision assessment.

Estimating Probabilities and Variance Summaries

A posterior probability is an expectation of an indicator function. For a parameter $\theta$ and a set $A$, $P(\theta\in A\mid y)=E\{I(\theta\in A)\mid y\}$. Its MCMC estimate is the fraction of draws in $A$. This fraction estimates a probability under the posterior; it is not a frequentist proportion in the original sample, and it does not measure the probability that the model or a null hypothesis is true.

Worked Example: Estimating a Posterior Probability

Suppose ten hypothetical, post-warmup draws of a random-intercept standard deviation $\sigma$, in log-odds units, are already sorted: $0.12,0.18,0.20,0.22,0.25,0.28,0.31,0.45,0.62,0.90$. Estimate $P(\sigma>0.30\mid y)$ using these draws.

Four of the ten values exceed $0.30$: $0.31,0.45,0.62,$ and $0.90$. Thus the empirical posterior probability estimate is $$ \widehat{P}(\sigma>0.30\mid y) =\frac{1}{10}\sum_{m=1}^{10}I(\sigma^{(m)}>0.30) =\frac{4}{10}=0.40. $$

The interpretation is that the sampled draws estimate a 0.40 posterior probability that the random-intercept standard deviation exceeds 0.30 log-odds units, conditional on this model and data. It is an approximation based on only ten draws; dependence or incomplete exploration could make it imprecise. It is not a statement that there is a 40% chance the fixed, unknown parameter changes value after the analysis.

For a variance component such as $\sigma^2$, draw-by-draw transformation matters. If the sampler stores $\sigma$, then the draws $(\sigma^{(m)})^2$ represent draws of the variance, and averaging them estimates $E(\sigma^2\mid y)$. Squaring the posterior mean of $\sigma$ instead estimates $\{E(\sigma\mid y)\}^2$, which is generally different. Similarly, applying a nonlinear transformation to a posterior mean does not generally give the posterior mean of the transformed parameter.

Worked Example: Mean Standard Deviation Versus Mean Variance

Suppose four illustrative draws of a random-intercept standard deviation, in log-odds units, are $0.2,0.4,0.6,0.8$. Their mean standard deviation is $$ \overline{\sigma}=\frac{0.2+0.4+0.6+0.8}{4}=\frac{2.0}{4}=0.5. $$ The corresponding average of the variance draws is $$ \overline{\sigma^2} =\frac{0.2^2+0.4^2+0.6^2+0.8^2}{4} =\frac{0.04+0.16+0.36+0.64}{4} =\frac{1.20}{4}=0.30 $$ in squared log-odds units. Squaring the mean standard deviation instead gives $(0.5)^2=0.25$, not $0.30$. Thus, the posterior mean variance is estimated by averaging the squared draws, not by squaring the estimated posterior mean standard deviation.

Checking the Target and Recognizing Failure Modes

A posterior summary depends on both the target and the chain’s exploration of it. Before reporting one, check that the stored variables correspond to the intended parameterization: standard deviation versus variance, covariance versus correlation, and conditional versus marginal quantity. Also verify that transformations and inequalities used in summaries are applied to each draw on the appropriate scale.

Warning signs include chains that remain in a narrow range for long periods, chains that show different persistent regions, or summaries that change substantially when later draws are included. A chain can appear visually stable yet fail to explore an important region, so plots and stability checks are evidence to consider rather than guarantees. If the posterior has multiple separated regions, for example, a chain trapped in one region can produce a precise-looking summary for the wrong part of the target.

The main alternatives depend on the structure of the posterior. When a posterior summary can be derived analytically, direct calculation avoids Monte Carlo error. Numerical integration can be effective for low-dimensional targets but becomes difficult as the number of unknowns grows. MCMC is useful for complex, high-dimensional targets because it can approximate posterior quantities without requiring independent draws, but its accuracy depends on correct targeting and adequate exploration. None of these computational approaches removes the need to assess the model and prior assumptions.

Common Mistakes and Practical Limits

Treating every draw as independent. A long chain is not equivalent to the same number of independent posterior draws. Account for dependence when judging numerical precision.
Confusing posterior spread with Monte Carlo error. Posterior spread describes uncertainty about a parameter; Monte Carlo error describes uncertainty in a computed summary.
Reporting a draw fraction without naming its event. State the probability being estimated and its context, such as $P(\sigma>0.30\mid y)$ for a standard deviation measured in log-odds units.
Summarizing the wrong scale. For a posterior mean variance, average the squared standard-deviation draws; do not square their average.
Assuming a stable-looking sequence proves convergence. Visual stability is helpful but cannot alone establish that all relevant posterior regions were explored.

MCMC results are approximations, and their reliability can be limited by dependence, slow exploration, or a poorly specified model. The posterior itself also depends on the likelihood and priors. A very large number of draws cannot correct a target that does not represent the intended scientific model.

Key takeaway: First specify the posterior target. Then interpret each chain summary as an approximation to a quantity under that target, and distinguish posterior uncertainty from the additional numerical uncertainty caused by finite, dependent draws.

Check Your Understanding

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

In a GLMM with fixed effects, group-level random effects, and a covariance parameter, what quantities belong in the joint posterior target?
A chain gives draws $\theta^{(1)},\ldots,\theta^{(M)}$. Write an estimator of $P(\theta>c\mid y)$.
Why is the usual sample standard deviation divided by $\sqrt{M}$ not generally a reliable Monte Carlo standard error for dependent draws?
If the sampler stores standard-deviation draws $\sigma^{(m)}$, how would you estimate $E(\sigma^2\mid y)$?
Ten draws contain six values above a specified threshold. State the empirical posterior probability estimate and one limitation of that estimate.

Solutions. (1) Include the fixed effects, random effects, covariance parameter, and any other unknown model quantities represented by the sampler; the target density is proportional to the likelihood times the random-effects density and prior. (2) Use $\frac{1}{M}\sum_{m=1}^{M}I(\theta^{(m)}>c)$. (3) That calculation assumes independent draws, whereas serial dependence changes the variance of the average and can make its error larger. (4) Calculate $\frac{1}{M}\sum_{m=1}^{M}(\sigma^{(m)})^2$, averaging the squared draws. (5) The estimate is $6/10=0.60$; with only ten draws it may be imprecise, particularly if the draws are dependent or the chain did not explore the target adequately.