Gibbs Sampling for Conjugate Variance Components

Why Gibbs Updates Can Be Direct

In a Gaussian hierarchical model, random effects and variance components are linked: the observed responses inform the random effects, and the random effects inform their variance. Gibbs sampling uses this structure by drawing one parameter block at a time from its full conditional distribution—the distribution of that block given the data and the current values of all other blocks. When the model and priors are conjugate, these conditional distributions have standard forms, so each block can be sampled directly rather than accepted or rejected.

The calculations below use a random-intercept model. The same logic applies to other Gaussian random-effect structures when the covariance form and priors support the required conditional distributions. Conjugacy is a property of a specified likelihood-prior pairing; it is not a general guarantee for every GLMM or every prior.

Model and Conjugate Full Conditionals

Suppose there are \(J\) groups and \(n_j\) Gaussian responses in group \(j\). Use the model

Conditional on the parameters, the errors are independent across observations, and the group effects are independent across groups. The mean \(\mu\) is assigned a normal prior, \(\mu\sim N(m_0,V_0)\). Independently, assign inverse-gamma priors to the two variances. We use the shape-rate convention

where \(a>0\) and \(b>0\). Thus, larger \(b\) changes the scale of the prior on the variance; \(b\) is not a rate multiplying \(v\) in an exponential density. Assume proper priors and a posterior that is proper.

To derive the update for \(b_j\), retain the likelihood terms from group \(j\) and multiply by its normal prior. Completing the square in \(b_j\) gives a normal conditional distribution:

The conditional precision is the sum of the data precision \(n_j/\sigma^2\) and prior precision \(1/\tau^2\). The conditional mean is a precision-weighted compromise between the group's observed deviations from \(\mu\) and the prior mean zero.

Given the random effects and variances, the conditional distribution of \(\mu\) is also normal. Its precision adds the prior precision to the observation precision:

Write \(N=\sum_j n_j\) for the total number of observations and define the residual sum of squares \(SSE=\sum_j\sum_i(y_{ij}-\mu-b_j)^2\). Define \(SSB=\sum_j b_j^2\), the sum of squared group effects. Multiplication of each Gaussian likelihood contribution by the corresponding inverse-gamma prior yields

Each observation contributes one Gaussian residual term to the update for \(\sigma^2\); each random effect contributes one normal term to the update for \(\tau^2\). Consequently, the shape increments are \(N/2\) and \(J/2\), respectively—not the number of groups for the residual variance or the number of observations for the random-effect variance.

Why a Full-Conditional Update Preserves the Posterior

The key justification is conditional resampling. Let \(\theta_k\) be the block being updated and let \(\theta_{-k}\) denote all remaining blocks. If the current joint state has posterior density \(\pi(\theta_k,\theta_{-k}\mid\mathbf y)\), the new value \(\theta_k'\) is drawn from \(\pi(\theta_k'\mid\theta_{-k},\mathbf y)\), holding \(\theta_{-k}\) fixed. The density of the new joint state is therefore

The integral is over the old value of the updated block, while the other blocks remain fixed. The resulting density is the target posterior density at the new value. Thus each exact full-conditional update leaves the posterior invariant, and a sequence of such updates also leaves it invariant. Invariance alone does not say that a finite run has reached the posterior; that requires adequate exploration of the posterior, a separate practical concern.

A Complete Numerical Gibbs Sweep

A sweep requires a specified starting state, a specified update order, and draws from the resulting conditional distributions. The next example uses two hypothetical groups with responses \(2,4\) and \(5,7\). Let \(\mu\sim N(0,100)\), and let both variance priors be \(\operatorname{IG}(2,1)\). Start at \(\mu=4\), \(\sigma^2=1\), and \(\tau^2=4\). Update \(b_1,b_2,\mu,\sigma^2,\tau^2\) in that order.

Assumptions, Diagnostics, and Alternatives

These direct updates depend on the stated model. The response likelihood must be Gaussian with common residual variance \(\sigma^2\); conditional on that variance, the errors must be independent. The random effects must be Gaussian with the stated variance structure, and the normal prior on \(\mu\) and inverse-gamma priors must have the specified parameterizations. If there are known offsets or predictors, they can be included in the linear predictor and the corresponding residuals, but the conditional formulas must then use those terms consistently.

Conjugacy can fail when the likelihood is not Gaussian, when random effects have a nonnormal prior, or when covariance parameters are assigned priors that do not yield standard full conditionals. An inverse-gamma prior is also not automatically a sensible prior merely because it is convenient; its implications should be considered on the variance scale and in the context of the application.

Useful implementation checks include verifying that all sampled variances are positive, that the residual and effect sums of squares use the current values specified by the update order, and that the shape increments count observations or effects correctly. Compare conditional formulas against simple cases, such as one group or zero residual sum of squares. In a run, inspect sampled values and trace behavior for signs of sticking, extreme sensitivity to starting values, or poor movement. Correct full conditionals do not by themselves ensure that a chain explores the posterior well.

Unlike a Metropolis-Hastings update, a Gibbs update from an exact full conditional needs no proposal distribution and has no accept-reject step. That advantage is available only where the full conditional can be sampled directly. When it cannot, other sampling methods may be needed. Neither direct sampling nor posterior invariance guarantees fast mixing, especially when effects and their variance components are strongly dependent.

Common Mistakes

• Integrating the wrong value in the invariance argument. The old value of the updated block is integrated out while the other blocks are held fixed. Integrating the other blocks instead gives a marginal density, not the joint density at the newly updated state.
• Using the wrong inverse-gamma convention. Under the shape-rate convention here, the kernel is \(\exp(-b/v)\). A formula from a source using a scale or a different parameterization must be translated before use.
• Using the wrong sample-size count. The residual variance shape increment is half the number of observations; the random-effect variance shape increment is half the number of effects.
• Treating a conditional draw as a point estimate. A Gibbs update is a random draw from a conditional distribution. Its value need not equal that distribution's mean or mode.