Marginalizing Random Effects in Bayesian Computation

Why Integrate Out Random Effects?

A joint model uses subject-specific random effects to connect repeated biomarker measurements with an event process. In Bayesian computation, those effects can be sampled along with the global parameters, or eliminated from the sampling target by integrating them out. The second strategy is called marginalization. It can reduce the dimension of the Markov chain and improve estimates of global parameters, but it replaces direct sampling with an integral that may be expensive or difficult to evaluate.

The key computational question is not simply whether random effects can be integrated out. It is whether the resulting marginal likelihood can be evaluated accurately and cheaply enough to justify doing so. A closed-form integral may make marginalization attractive; a high-dimensional integral for every subject may make sampling more practical.

Definition: Let $\theta$ denote all global parameters, including any random-effects distribution parameters, and let $b_i$ be the random effects for subject $i$. If $D_i$ denotes that subject’s observed data, $L_i(\theta,b_i)$ is its likelihood contribution conditional on $b_i$, and $p(b_i\mid\theta)$ is the random-effects distribution, the marginal likelihood contribution is $$ m_i(\theta)=\int L_i(\theta,b_i)\,p(b_i\mid\theta)\,db_i. $$ It averages the conditional likelihood over the random-effects distribution; it does not set the random effect to its mean.

Two Posterior Targets

For a cohort of $n$ subjects, a common model assumes that subjects’ data are conditionally independent given the global parameters and their own random effects, and that the random effects are conditionally independent across subjects given the global parameters. Under these assumptions, the posterior that retains all random effects has density proportional to

p(\theta,b_1,\ldots,b_n\mid D) \propto p(\theta)\prod_{i=1}^{n}p(b_i\mid\theta)L_i(\theta,b_i).

A sampler targeting this posterior updates $\theta$ and the $b_i$'s. Alternatively, integrate each $b_i$ out. The marginal posterior for the global parameters is proportional to

p(\theta\mid D) \propto p(\theta)\prod_{i=1}^{n} \int p(b_i\mid\theta)L_i(\theta,b_i)\,db_i.

This expression follows by integrating the joint posterior over all $b_i$'s. The product of integrals is valid because, conditional on $\theta$, each subject’s random-effects prior and data contribution involve that subject’s $b_i$. If subjects share a latent effect or have dependent random effects, the appropriate joint integral may not factor into subject-specific terms.

In a joint model, the conditional subject likelihood may be written $L_i(\theta,b_i)=L_{Yi}(\theta,b_i)L_{Ti}(\theta,b_i)$, where the factors represent longitudinal and event data. Both factors depend on the same random effect, so the marginal contribution integrates their product. Integrating each factor separately and multiplying the results generally gives a different likelihood.

Worked Example: An Analytic Marginal Likelihood

Suppose a hypothetical subject has one biomarker measurement $y=2.4$, in arbitrary biomarker units, and the conditional model is $Y\mid b\sim N(\theta+b,\sigma^2)$. Let $\theta=1.5$ units, $b\sim N(0,\tau^2)$, $\sigma=0.6$ units, and $\tau=0.8$ units. Assume the measurement error is independent of $b$. The normal random effect and normal measurement error make the marginal distribution available exactly.

Because $Y=\theta+b+\epsilon$, where $b\sim N(0,0.8^2)$ and $\epsilon\sim N(0,0.6^2)$ are independent, the marginal distribution is $Y\sim N(\theta,\tau^2+\sigma^2)$. Thus the marginal variance is $0.8^2+0.6^2=0.64+0.36=1.00$ square units. The likelihood density at the observed value is

m(\theta)=\frac{1}{\sqrt{2\pi(1.00)}}\exp\left\{-\frac{(2.4-1.5)^2}{2(1.00)}\right\} =\frac{1}{\sqrt{2\pi}}\exp(-0.405) \approx 0.2661

The value $0.2661$ is a density in inverse biomarker units, not a probability of observing exactly $2.4$. Here marginalization replaces a one-dimensional integral with a closed-form density evaluation. If the subject-specific effect is needed after fitting, its conditional posterior is still available: $b\mid y,\theta$ is normal with mean $\tau^2(\tau^2+\sigma^2)^{-1}(y-\theta)=0.64(0.9)=0.576$ units and variance $\tau^2\sigma^2/(\tau^2+\sigma^2)=0.64(0.36)=0.2304$ square units. The effect has been removed from the sampling target, not made unrecoverable.

Shared Effects Must Be Integrated Jointly

The same random effect can induce dependence between a subject’s biomarker data and event data, even when those data are conditionally independent given the effect and global parameters. Marginalization preserves that dependence only if the joint conditional likelihood is integrated as a whole.

Worked Example: Integrating a Shared Effect Across Two Outcomes

Consider a hypothetical subject with two binary observations: $Y=1$ for a longitudinal outcome and $D=1$ for an event by a fixed assessment time. For illustration, let the shared effect take values $b=-1$ and $b=1$, each with prior probability $0.5$. Given $b$, assume the two outcomes are independent, with $P(Y=1\mid b)=\operatorname{logit}^{-1}(b)$ and $P(D=1\mid b)=\operatorname{logit}^{-1}(-1+b)$, where $\operatorname{logit}^{-1}(x)=1/(1+e^{-x})$.

For $b=-1$, the two success probabilities are $\operatorname{logit}^{-1}(-1)=0.2689$ and $\operatorname{logit}^{-1}(-2)=0.1192$. For $b=1$, they are $\operatorname{logit}^{-1}(1)=0.7311$ and $\operatorname{logit}^{-1}(0)=0.5000$. The marginal likelihood of the observed pair is therefore

P(Y=1,D=1) =0.5(0.2689)(0.1192)+0.5(0.7311)(0.5000) \approx 0.1988.

If the outcomes were integrated separately and the two marginal probabilities multiplied, the result would be $[0.5(0.2689+0.7311)][0.5(0.1192+0.5000)]=0.1548$, not $0.1988$. That calculation wrongly treats the outcomes as independent after removing their shared effect. Conditional independence given $b$ does not imply marginal independence after $b$ is integrated out.

Why Marginalization Can Help

Removing the $b_i$'s reduces the dimension of the state sampled by an MCMC algorithm. It can also avoid slow movement along strong posterior dependence between a random effect and global parameters. A further benefit is the Rao–Blackwell principle: replacing a random quantity by its conditional expectation can reduce Monte Carlo variance.

For a posterior quantity $g(\theta,b)$, the law of total variance gives

\operatorname{Var}\{g(\theta,b)\mid D\} = \operatorname{Var}\{E[g(\theta,b)\mid\theta,D]\mid D\} + E[\operatorname{Var}\{g(\theta,b)\mid\theta,D\}\mid D].

The second term is nonnegative. Thus, when estimating the posterior expectation of $g$, averaging over $b$ conditionally rather than using a single noisy draw of $b$ can reduce variance. This is a statement about Monte Carlo precision, not a guarantee that the marginal sampler runs faster overall. Evaluating each integrated likelihood may be costly, and the marginal posterior for $\theta$ may still be difficult to explore.

A practical sampler can also combine strategies: sample global parameters from a marginal target, then draw subject effects from $p(b_i\mid\theta,D_i)$ when subject-specific summaries are needed. Alternatively, some effects may be integrated analytically while others remain explicit. Such choices must preserve the intended posterior and account for any approximation used in the integration.

Worked Example: When Tensor Quadrature Becomes Expensive

Suppose a hypothetical cohort has $n=200$ subjects, each with a five-dimensional random effect. Assume no analytic integral is available and a tensor-product Gaussian quadrature rule with 9 nodes per dimension is being considered. For one subject, the rule requires $9^5=59{,}049$ node combinations. Across the cohort, this is $200(59{,}049)=11{,}809{,}800$ likelihood evaluations for one evaluation of the marginal log posterior, before accounting for the other computation.

For comparison, with a one-dimensional effect and the same 9 nodes, the corresponding total is $200(9)=1{,}800$ node evaluations. The five-dimensional calculation is $11{,}809{,}800/1{,}800=6{,}561$ times as many evaluations. These counts follow directly from a full tensor product; adaptive or sparse-grid methods can change the cost, as can model structure. In this hypothetical setting, full tensor quadrature may be practical for the scalar effects but unattractive for the five-dimensional effects. Sampling the effects avoids enumerating all node combinations, although it introduces 1,000 latent scalar coordinates ($200\times5$) and can mix poorly. The arithmetic signals a trade-off; it does not decide which sampler will be faster or more accurate without implementation and diagnostic checks.

Choosing and Checking a Strategy

Strategy	Potential advantage	Main computational risk
Sample each $b_i$	Avoids an explicit high-dimensional integral; subject effects are directly available in posterior draws.	Many latent coordinates and potentially strong dependence can slow exploration.
Integrate each $b_i$ analytically	Removes latent coordinates without numerical integration error.	Closed forms are available only for particular model structures.
Integrate numerically	Can reduce sampling dimension while retaining a likelihood-based target.	Repeated integration can be costly; inaccurate integration can distort the posterior.

Before relying on numerical marginalization, check that the integration method accurately represents the full contribution $p(b_i\mid\theta)L_i(\theta,b_i)$, including its tails and any sharp concentration. Increase quadrature accuracy and compare important posterior summaries; substantial changes indicate that the chosen approximation is not yet reliable. For Monte Carlo integration, assess the stability and precision of the estimated integral, not just the apparent stability of the outer chain. Calculations on the log scale should combine weighted terms with numerically stable log-sum-exp methods when needed.

Also verify the model assumptions behind the factorization. Separate subject integrals require the stated conditional independence structure. Longitudinal and event components sharing $b_i$ must remain together inside each subject’s integral. If the integral is approximated, the approximation error is distinct from ordinary MCMC error: a long chain cannot correct a systematically inaccurate likelihood.

Common Pitfalls and Limitations

Replacing integration with the mean effect. Using $L_i(\theta,E[b_i\mid\theta])$ is generally not equal to $m_i(\theta)$. In general, the likelihood is nonlinear in $b_i$, so averaging the likelihood is not the same as evaluating it at an average effect.
Integrating likelihood components separately. If both components depend on the same random effect, integrate their product. The discrete example shows that multiplying separately marginalized contributions changes the model’s implied dependence.
Assuming fewer sampled parameters means less total computation. Marginalization may lower chain dimension but make every target evaluation much more expensive. Compare the cost and accuracy of the integral with the cost of updating the effects.
Discarding random effects permanently. Marginalization removes them from the main sampling target; conditional posterior draws or summaries can still be obtained given draws of the global parameters.

There is no universally best strategy. Analytic integration is especially attractive when available and tractable. Numerical integration is most plausible when the effective integration dimension is modest and its accuracy can be checked. Sampling is often preferable when integration is high-dimensional or costly, provided the resulting chain can be explored reliably. Hybrid approaches can be useful, but they do not eliminate the need to verify the target distribution and the accuracy of any approximations.

Key takeaway: Marginalizing random effects replaces sampling them with a likelihood averaged over their distribution. This can reduce dimension and Monte Carlo noise, but only pays off when the resulting integrals are tractable and accurate. Preserve shared-effect dependence by integrating the complete subject likelihood.

Check Your Understanding

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

State the assumptions that allow the marginal posterior to factor into a product of subject-specific integrals.
Why is the product of separately integrated longitudinal and event likelihoods generally incorrect when both use the same random effect?
A one-dimensional random effect has a closed-form marginal likelihood. What is one computational advantage of integrating it out, and what is one reason to retain it in the sampler?
With 6 quadrature nodes per dimension and a four-dimensional effect, how many tensor-product nodes are required for one subject?
Does a stable MCMC trace establish that a numerically integrated likelihood is accurate? Explain.

Solutions. (1) Conditional on global parameters, subjects’ data and random effects must factor by subject, with each contribution depending on its own effect; shared latent structure can prevent this factorization. (2) The shared effect induces marginal dependence. The correct contribution integrates the product of the conditional likelihood factors, rather than multiplying their separate averages. (3) Integration removes that latent coordinate and may reduce posterior dependence or Monte Carlo noise; retaining it can be more practical if the integral is expensive, and it makes subject-specific effects directly available in draws. (4) $6^4=1{,}296$ nodes. (5) No. A chain can be stable for an inaccurately approximated target. Integration accuracy must be checked separately, for example by increasing quadrature accuracy and comparing posterior results.

Strategy	Potential advantage	Main computational risk
Sample each \(b_i\)	Avoids an explicit high-dimensional integral; subject effects are directly available in posterior draws.	Many latent coordinates and potentially strong dependence can slow exploration.
Integrate each \(b_i\) analytically	Removes latent coordinates without numerical integration error.	Closed forms are available only for particular model structures.
Integrate numerically	Can reduce sampling dimension while retaining a likelihood-based target.	Repeated integration can be costly; inaccurate integration can distort the posterior.