Bayesian Posterior Construction for a GLMM

From a GLMM likelihood to a posterior

A generalized linear mixed model describes responses through a conditional distribution and a linear predictor containing both fixed and random effects. Bayesian posterior construction adds probability models for those unknown effects and for the variance components. The result is a joint distribution for the unknown quantities, conditional on the observed data.

The distinction between conditional and marginal likelihoods matters. The conditional likelihood treats the random effects as unknown quantities in the model. A marginal likelihood instead integrates them out. A Bayesian posterior that includes the random effects as unknowns uses the conditional likelihood together with their prior distribution; the integration appears only if the goal is a posterior that no longer includes those random effects.

Definition: For observed responses $y$, unknown quantities $\theta$, likelihood $p(y\mid\theta)$, and prior density $p(\theta)$, the posterior density is $p(\theta\mid y)=p(y\mid\theta)p(\theta)/p(y)$. The posterior kernel is the product $p(y\mid\theta)p(\theta)$, viewed as a function of $\theta$, with factors that are constant in $\theta$ omitted. It is proportional to the posterior density, not generally a probability by itself.

For a common clustered GLMM, let $j=1,\ldots,J$ index clusters and $i=1,\ldots,n_j$ index observations within cluster $j$. Write the linear predictor as $\eta_{ij}=x_{ij}^{\mathsf T}\beta+z_{ij}^{\mathsf T}b_j$, where $\beta$ is a fixed-effect vector and $b_j$ is the random-effect vector for cluster $j$. The response distribution and link determine the conditional mean $\mu_{ij}$, with $g(\mu_{ij})=\eta_{ij}$.

Assume that, conditional on all random effects and model parameters, the responses are independent and follow the specified GLM response family. Assume also that the cluster random effects are independent draws from a common distribution conditional on its variance-component parameters. These are modeling assumptions, not consequences of Bayes' rule. A scalar random intercept might satisfy $b_j\mid v\sim N(0,v)$; a vector of random effects might satisfy $b_j\mid G\sim N_q(0,G)$, where $G$ is a covariance matrix.

Constructing the joint posterior kernel

Under the conditional independence assumptions, the likelihood is the product of the observation-level response densities. In the vector random-effect case, combining this likelihood with the random-effect density and independent priors for $\beta$ and $G$ gives the following joint posterior kernel:

p(\beta,b_1,\ldots,b_J,G\mid y) \ \propto\ \left\{\prod_{j=1}^J\prod_{i=1}^{n_j} p(y_{ij}\mid \beta,b_j)\right\} \left\{\prod_{j=1}^J p(b_j\mid G)\right\} p(\beta)p(G).

This expression assumes prior independence between $\beta$ and $G$. If the prior specification instead makes them dependent, replace $p(\beta)p(G)$ by their joint prior density $p(\beta,G)$. Likewise, if a dispersion parameter or another unknown model parameter is included, its likelihood contribution and prior must be included. Any factor may be dropped from the kernel only when it is constant with respect to every unknown quantity being inferred.

For $b_j\mid G\sim N_q(0,G)$, the random-effect contribution is

\prod_{j=1}^J p(b_j\mid G) = (2\pi)^{-Jq/2}|G|^{-J/2} \exp\left\{-\frac{1}{2}\sum_{j=1}^J b_j^{\mathsf T}G^{-1}b_j\right\}.

The determinant term $|G|^{-J/2}$ depends on $G$, so it cannot be discarded when constructing a posterior for the variance components. The quadratic terms also depend on $G$ and the random effects. In a scalar random-intercept model, these contributions reduce to a product of normal densities involving $v^{-J/2}$ and $\sum_j b_j^2/v$, where $v$ is the random-intercept variance.

Key distinction: The joint posterior retains the random effects as unknowns. To obtain a posterior for $(\beta,G)$ alone, integrate the joint kernel over every random-effect vector: $$ p(\beta,G\mid y)\ \propto\ p(\beta)p(G) \int \left\{\prod_{j,i}p(y_{ij}\mid\beta,b_j)\right\} \left\{\prod_jp(b_j\mid G)\right\}\,db_1\cdots db_J. $$ The integral is generally not equal to the likelihood evaluated at $b_j=0$, nor to the likelihood evaluated at an estimated value of $b_j$.

Worked examples of posterior construction

Worked Example: A Binary Random-Intercept Model

Consider one cluster with two binary observations, $y_1=1$ and $y_2=0$. Let $Y_i\mid\beta,b$ be Bernoulli with logit $ \operatorname{logit}(p_i)=\beta+b$. Use the proper priors $\beta\sim N(0,4)$, $b\mid v\sim N(0,v)$, and $v\sim\operatorname{Exponential}(1)$, whose density is $p(v)=e^{-v}$ for $v>0$. The priors are taken independent except for the stated dependence of $b$ on $v$.

For these observations the conditional likelihood is $p(1\mid\beta,b)p(0\mid\beta,b)=p(1-p)$, where $p=\operatorname{logit}^{-1}(\beta+b)$. Thus the posterior kernel is

L(\beta,b)\, \frac{1}{\sqrt{8\pi}}\exp\left(-\frac{\beta^2}{8}\right) \frac{1}{\sqrt{2\pi v}}\exp\left(-\frac{b^2}{2v}\right)e^{-v}, \qquad v>0.

Evaluate it at $\beta=\log 2$, $b=-\log 2$, and $v=1$. Then $\beta+b=0$, so $p=0.5$ and $L=0.5(1-0.5)=0.25$. The fixed-effect prior density is $(8\pi)^{-1/2}\exp\{-(\log 2)^2/8\}\approx0.1878$. The random-effect density is $(2\pi)^{-1/2}\exp\{-(\log 2)^2/2\}\approx0.313748$, and $p(v=1)=e^{-1}\approx0.367879$. Multiplying the unrounded factors gives

0.25(0.187847)(0.313748)(0.367879)\ \approx\ 0.005420.

This is a value of the joint posterior kernel at one point, not the probability that the parameters equal those values. The posterior density must be normalized over its parameter space to become a density.

Worked Example: A Poisson Random-Intercept Model

Suppose one cluster has counts $y_1=0$ and $y_2=2$, with $Y_i\mid\beta,b\sim\operatorname{Poisson}(\lambda_i)$ and $\log(\lambda_i)=\beta+b$. Use the same priors as in the preceding example. Conditional independence gives

L(\beta,b)= \frac{e^{-\lambda_1}\lambda_1^0}{0!} \frac{e^{-\lambda_2}\lambda_2^2}{2!}.

At $\beta=\log 2$, $b=0$, and $v=1$, both rates equal $2$. Therefore $L=e^{-2}(e^{-2}2^2/2)=2e^{-4}\approx0.0366313$. The fixed-effect prior density is approximately $0.187847$, the random-effect density at $b=0,v=1$ is $(2\pi)^{-1/2}\approx0.398942$, and the variance prior density is $e^{-1}\approx0.367879$. Thus the kernel value is

(2e^{-4})(0.187847)(0.398942)(0.367879) \ \approx\ 0.001010.

The Poisson factorials are included in the likelihood shown. For fixed observed counts they do not depend on $\beta$, $b$, or $v$, so they could be omitted from a kernel for these parameters. They should not be omitted if the likelihood is being reported in full or if the quantities being compared make those factors parameter-dependent.

Worked Example: A Grouped Binomial Observation

A cluster contributes one grouped observation: $y=2$ successes among $m=3$ trials. Let $Y\mid\beta,b\sim\operatorname{Binomial}(3,p)$, with $\operatorname{logit}(p)=\beta+b$, and again use $\beta\sim N(0,4)$, $b\mid v\sim N(0,v)$, and $v\sim\operatorname{Exponential}(1)$.

At $\beta=0$, $b=0$, and $v=1$, the success probability is $p=0.5$. The full binomial likelihood is $L={3\choose2}(0.5)^2(0.5)^1=3(0.125)=0.375$. The prior densities at this point are $1/\sqrt{8\pi}\approx0.199471$, $1/\sqrt{2\pi}\approx0.398942$, and $e^{-1}\approx0.367879$. Hence

0.375(0.199471)(0.398942)(0.367879)\ \approx\ 0.010976.

The binomial coefficient is constant in the model parameters here, so a kernel omitting that coefficient would be proportional to, but numerically different from, the full kernel. Either convention gives the same normalized posterior when the omitted factor is genuinely constant over the parameter space.

Why the construction works

The factorization follows directly from the conditional model and Bayes' rule. Given $\beta$ and the cluster effects, the observation densities multiply because the responses are conditionally independent. Given $G$, the cluster effects have their specified random-effect densities. The prior supplies the remaining factors. Bayes' rule divides this product by the marginal density $p(y)$, which does not depend on the unknown parameters; it is therefore omitted when writing a posterior kernel.

Write the conditional response model.
Specify the response family and linear predictor, then form the product of conditional response densities for the observed data.

Write the random-effect distribution.
Include its full density, including determinant or variance-scale terms that depend on the unknown variance components.

Multiply by the prior.
Add the fixed-effect and variance-component prior densities, using a joint prior instead if the parameters are not prior-independent.

Remove only constants in all unknowns.
The remaining product is the posterior kernel. Normalize it only when a normalized posterior density or posterior probability is needed.

Checks, alternatives, and limitations

A correct algebraic kernel cannot repair a misspecified conditional response distribution, link, random-effect structure, or independence assumption. Before computation, verify that the response likelihood matches the data type, the predictor uses the intended covariates, and each random effect is assigned to the correct cluster. Check that each prior has the intended support and that the variance-component density is expressed in the same parameterization as the model. A prior on a variance $v$ is not interchangeable with a prior on a standard deviation without accounting for the change of variable.

A useful computational check is to evaluate the log kernel rather than multiply many small densities directly. The log kernel is the sum of the log likelihood and log prior contributions. Compare its individual terms at a few valid parameter values; an omitted determinant term, sign error in a quadratic form, or use of the wrong conditional mean can then be easier to detect. Also check that the posterior is proper: a finite-looking kernel at selected points does not establish that its integral over the full parameter space is finite.

There are several related approaches, but they answer different questions. A frequentist marginal likelihood integrates the random effects out of the conditional likelihood and is used for likelihood-based estimation. The joint Bayesian posterior retains random effects and combines their conditional distribution with the likelihood and priors; integrating the posterior over the random effects gives the corresponding marginal posterior for population-level parameters. Penalized quasi-likelihood uses an approximation rather than this exact conditional response likelihood in its construction, so a kernel built from the exact likelihood should not be described as a PQL posterior.

Posterior construction specifies the target distribution, not the numerical method for exploring it. Integration and simulation may be difficult in high dimensions, especially with weakly informed variance components or sparse outcomes. A computed chain or approximation should be assessed for convergence and sensitivity to modeling choices; a well-written kernel alone does not demonstrate that computation has succeeded. Prior choice also matters: the examples use simple proper priors to make the multiplication transparent, not to claim that those priors are universally suitable.

Common mistakes

Using only the likelihood. A likelihood is not a Bayesian posterior kernel. Include the prior density for every unknown parameter.
Dropping the random-effect density. If $b_j$ is included in the joint posterior, its distribution conditional on the variance components is part of the model and must appear.
Dropping variance-dependent normalizing terms. Terms such as $|G|^{-J/2}$ depend on the variance components and cannot be treated as constants.
Calling a kernel value a probability. A density value is not the probability of a single parameter point. Posterior probabilities require integration over a specified set.
Substituting an estimated random effect for integration. A likelihood evaluated at a fitted random effect is not generally the marginal likelihood or the marginal posterior.

Key takeaway: A GLMM joint posterior kernel is the product of the conditional response likelihood, the random-effect density, and the prior densities for the remaining unknowns. State the conditioning clearly, preserve every parameter-dependent factor, and distinguish the kernel from both a normalized posterior and a marginal likelihood.

Check Your Understanding

For each question, identify the relevant factorization or calculation before reading the solution.

In a joint posterior for $(\beta,b,G)$, why must the density $p(b\mid G)$ appear even when $b$ is called a random effect?
Suppose the likelihood is multiplied by $p(\beta)p(G)$ but no random-effect density is included. What model quantity is missing?
For a scalar random intercept with $J$ clusters and variance $v$, which part of the normal random-effect density shows that $v^{-J/2}$ cannot be dropped from a posterior kernel for $v$?
Explain why the value of a posterior kernel at one point is not the posterior probability of that parameter value.

Solutions. (1) In the joint posterior, $b$ is an unknown quantity whose conditional distribution is part of the hierarchical model. Bayes' rule therefore contributes $p(b\mid G)$ alongside the response likelihood and other priors. (2) The model's random-effect distribution $p(b\mid G)$ is missing; without it, the expression does not specify the intended joint hierarchical posterior. (3) The product of $J$ normal densities contains the factor $(2\pi v)^{-J/2}$, including $v^{-J/2}$, which varies with $v$. (4) A point has probability zero under a continuous posterior distribution. The kernel is an unnormalized density value; a probability is obtained by integrating a normalized density over a set of parameter values.