Marginal Likelihood Integrals in Nonlinear Mixed-Effects Models

Why Integrate Over Random Effects?

In a nonlinear mixed-effects model, each subject has latent random effects that influence that subject’s observations. The conditional likelihood describes the data for a specified value of those effects. But when the effects have not been observed, the likelihood for the observed data must account for every possible value, weighted by how plausible that value is under the population distribution. This averaging is the individual-level marginal likelihood integral.

The integral is the target likelihood contribution before a computational approximation is chosen. Its definition also clarifies two different operations: integrate over the random effects for each subject, then combine the resulting subject-level contributions across subjects. Confusing those operations can change both the model and the likelihood being optimized.

Definition: For subject $i$, let $y_i$ denote the subject’s observed data, let $\eta_i$ be a $q$-dimensional random-effects vector, and let $\theta$ denote all model parameters, including fixed effects, residual-error parameters, and random-effects distribution parameters. Write $f(y_i\mid\eta_i;\theta)$ for the conditional data density or probability mass and $g(\eta_i;\theta)$ for the random-effects density. The individual marginal likelihood is $$ L_i(\theta)=p(y_i;\theta)=\int_{\mathbb{R}^q} f(y_i\mid\eta_i;\theta)\,g(\eta_i;\theta)\,d\eta_i. $$ The integrand is the product $f(y_i\mid\eta_i;\theta)g(\eta_i;\theta)$, and the integration dimension is $q$, the number of random effects for that subject.

The formula applies to continuous observations, discrete observations, or a mixture, provided the conditional likelihood is defined with respect to the appropriate measure and the random-effects density is normalized. A common model assumes $\eta_i\sim N(0,\Omega)$, where $\Omega$ is a positive-definite covariance matrix. In that case, $g$ is a multivariate normal density. Gaussian random effects are common but are not required by the definition: the integral uses whichever random-effects distribution the model specifies.

The conditional likelihood itself may factor over observations within a subject. For example, if the observations $Y_{ij}$, $j=1,\ldots,n_i$, are conditionally independent given $\eta_i$, then $f(y_i\mid\eta_i;\theta)=\prod_{j=1}^{n_i} f(y_{ij}\mid\eta_i;\theta)$. The integral is still over the subject’s random-effects vector, not over each observation separately. The conditional independence assumption must be part of the model; if the residual process has additional within-subject dependence, that dependence belongs in the conditional likelihood rather than being discarded.

From Subject Contributions to the Population Likelihood

Suppose there are $N$ subjects. To obtain the usual product likelihood, assume subjects are independent: their random effects are independently drawn from the specified population distribution, and their data are conditionally generated according to the model for their own effects. Then the joint density of the data and random effects factors across subjects. Integrating the effects subject by subject gives the population likelihood.

\begin{aligned} p(y_1,\ldots,y_N,\eta_1,\ldots,\eta_N;\theta) &=\prod_{i=1}^{N} f(y_i\mid\eta_i;\theta)g(\eta_i;\theta),\\ L(\theta) &=\int\cdots\int \prod_{i=1}^{N} f(y_i\mid\eta_i;\theta)g(\eta_i;\theta) \,d\eta_1\cdots d\eta_N\\ &=\prod_{i=1}^{N}\int f(y_i\mid\eta_i;\theta)g(\eta_i;\theta)\,d\eta_i =\prod_{i=1}^{N}L_i(\theta). \end{aligned}

The factorization in the last line follows because each factor in the integrand depends on only one subject’s random effects, and the integration region is the product of the subject-specific spaces. Equivalently, repeated integration of this nonnegative product separates into the product of the individual integrals. If subjects share latent effects or have other cross-subject dependence, this factorization does not hold as written; the joint model and its integral must represent that dependence.

The log likelihood is therefore a sum, $\ell(\theta)=\sum_{i=1}^{N}\log L_i(\theta)$. This is a computationally useful representation, but it does not mean that observations within a subject are independent after the random effects have been integrated out. A shared $\eta_i$ can induce marginal dependence among that subject’s repeated observations, even when those observations were conditionally independent given $\eta_i$.

Worked Example: A One-Dimensional Gaussian Integral

Consider one hypothetical subject with one observation $y=1$. Let $\eta\sim N(0,1)$ and, conditional on $\eta$, let $Y\mid\eta\sim N(\eta,1)$. The fixed parameters are treated as specified, so the subject-level integral has one integration variable.

The conditional density and random-effects density are $f(1\mid\eta)=(2\pi)^{-1/2}\exp\{-(1-\eta)^2/2\}$ and $g(\eta)=(2\pi)^{-1/2}\exp\{-\eta^2/2\}$. Thus the integrand is their product and

L(1)=\int_{-\infty}^{\infty} \frac{1}{2\pi}\exp\left\{-\frac{(1-\eta)^2+\eta^2}{2}\right\}\,d\eta.

Complete the square: $(1-\eta)^2+\eta^2=2(\eta-\tfrac12)^2+\tfrac12$. Therefore,

\begin{aligned} L(1) &=\frac{e^{-1/4}}{2\pi} \int_{-\infty}^{\infty}e^{-(\eta-1/2)^2}\,d\eta\\ &=\frac{e^{-1/4}}{2\pi}\sqrt{\pi} =\frac{e^{-1/4}}{\sqrt{4\pi}} \approx 0.2197. \end{aligned}

This is a density value for the observed continuous outcome, not the probability that $Y$ equals exactly 1. It agrees with the convolution result $Y\sim N(0,2)$: the observation has marginal density $(4\pi)^{-1/2}\exp(-1/4)$. The random-effects integration is one-dimensional because there is one random effect, even though the model has an observation and a residual error as well.

Dimension and the Shape of the Integrand

The integration dimension is set by the number of latent random effects in $\eta_i$, not by the number of records $n_i$. More observations can make the conditional likelihood more informative about $\eta_i$, and can change the shape of the integrand substantially, but do not by themselves add integration coordinates. Conversely, a model with several random effects per subject requires a multivariate integral even with very few observations.

In a standard Gaussian random-effects model with $\eta_i\sim N(0,\Omega)$, the integrand can be written explicitly as

f(y_i\mid\eta_i;\theta) (2\pi)^{-q/2}|\Omega|^{-1/2} \exp\left(-\frac{1}{2}\eta_i^\mathsf{T}\Omega^{-1}\eta_i\right).

The conditional likelihood contributes the data fit, while the random-effects density weights each point according to its population plausibility. Neither factor alone is the marginal likelihood. In particular, evaluating the conditional likelihood only at an estimated or most likely individual effect is not the same operation as integrating it over the random-effects distribution.

Worked Example: Two Random Effects and Two Observations

For a hypothetical subject, observe $y_{i1}=0$ and $y_{i2}=0$ at covariate values $x_1=-1$ and $x_2=1$. Suppose $Y_{ij}\mid\eta_i\sim N(\eta_{i1}+x_j\eta_{i2},1)$, independently for $j=1,2$ conditional on $\eta_i$, and $\eta_{i1},\eta_{i2}$ are independent $N(0,1)$. Here $q=2$. Writing $\eta_i=(a,b)^\mathsf{T}$, the integral is

L_i=\int_{\mathbb{R}^2} \frac{1}{2\pi} \exp\left[-\frac{a^2+b^2}{2}\right] \prod_{x\in\{-1,1\}}\frac{1}{\sqrt{2\pi}} \exp\left[-\frac{(0-a-xb)^2}{2}\right]\,da\,db.

For a direct evaluation, the random contribution to the two observations is $(a-b,a+b)^\mathsf{T}$. Its covariance matrix is $ZZ^\mathsf{T}=2I_2$, where the rows of $Z$ are $(1,-1)$ and $(1,1)$. Adding the conditional residual covariance $I_2$ gives marginal covariance $3I_2$. The marginal data vector is therefore $N(0,3I_2)$, so its density at $(0,0)^\mathsf{T}$ is

L_i=\frac{1}{(2\pi)^{2/2}|3I_2|^{1/2}} =\frac{1}{2\pi\cdot 3} =\frac{1}{6\pi} \approx 0.05305.

The calculation confirms that the two-dimensional integral uses $a$ and $b$. The two observations determine the dimension neither by being counted one-for-one nor by requiring a separate random-effects vector for each observation. Their conditional independence is used in the product inside the integrand; their marginal covariance is not $I_2$, because the shared random effects also contribute variation.

Combining Contributions and Checking the Likelihood

Once each $L_i(\theta)$ has been defined, independent subjects contribute multiplicatively to $L(\theta)$. This construction is central to estimation: changing $\theta$ changes both the conditional data model and, generally, the random-effects distribution. The likelihood evaluates the probability density or mass of the observed data under that full population model.

Worked Example: A Two-Subject Population Likelihood

Use the one-random-effect model from the first example for two independent hypothetical subjects, with observed values $y_1=0$ and $y_2=1$. The marginal distribution is $Y_i\sim N(0,2)$, so

L_1=\frac{1}{\sqrt{4\pi}}\approx 0.2820948,\qquad L_2=\frac{e^{-1/4}}{\sqrt{4\pi}}\approx 0.2196956.

Independence across subjects gives

L(\theta)=L_1L_2 =\frac{e^{-1/4}}{4\pi} \approx 0.0619750, \qquad \ell(\theta)=\log L(\theta) =-\frac14-\log(4\pi) \approx -2.781.

The product is the joint density at the two observed outcomes under the specified model and parameter values. The corresponding log likelihood is a sum of the two individual log densities. Using unrounded values for the product gives approximately $0.0619750$; multiplying the displayed rounded values gives the same result to the shown precision.

Assumptions, Alternatives, and Failure Modes

The marginal integral is a definition of the likelihood under the model, not a numerical method. To use the product form, the subjects must be independent under the model. Within each subject, the conditional likelihood must correctly encode the observation process and any residual dependence. The random-effects density must be proper and normalized, and the integral must exist and be finite at the parameter values considered. For continuous observations, a likelihood contribution is a density and can exceed 1; it is not required to be a probability bounded by 1.

An exact analytic integral is available for some models, as in the Gaussian examples. In nonlinear models, the integral often has no closed form. One can then use an approximation to the integral, or a numerical integration approach. These approaches target the same marginal likelihood when correctly implemented, but their finite computational results can differ. The earlier material on choosing an approximation concerns that distinction; the defining integral remains the reference quantity. Replacing the integral with a conditional likelihood evaluated at a single estimated effect is a different objective unless a specific approximation justifies that replacement.

Several checks help identify a malformed setup before considering numerical accuracy:

Check the dimension. Count the components of $\eta_i$. If there are $q$ random effects, the individual integral is over $q$ coordinates, even when the subject has many observations.
Check both factors. The integrand must include the conditional likelihood and the random-effects density. Omitting the density fails to weight effects by their population distribution and generally does not define the marginal likelihood.
Check the level of independence. Conditional independence of repeated observations does not imply their marginal independence after a shared random effect is integrated out. Independence across subjects is a separate assumption used to obtain the population product.
Check probability versus density. For continuous outcomes, $L_i$ is a density value at the observed data. Calling it the probability of observing that exact value is incorrect.
Check the parameter dependence. The random-effects distribution may depend on estimated covariance parameters. Treating its density as a fixed weight while estimating those parameters would not evaluate the stated model likelihood.

A common mistake is to multiply all conditional observation densities over all subjects and then integrate over one shared random effect, or to multiply conditional densities evaluated at subject-specific estimates. The first construction imposes a common latent effect rather than independent subject effects; the second omits the population averaging. The correction is to form one integral per subject, using that subject’s data and random-effects distribution, and then multiply those marginal contributions only under the across-subject independence assumption.

Key takeaway: The individual marginal likelihood is an integral of conditional data fit times the random-effects density over the subject’s random-effects vector. Its dimension is the number of random effects. Under independent subjects, the population likelihood is the product of these subject-level integrals.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

A model gives each subject three random effects and records 12 observations per subject. What is the dimension of the subject-level marginal likelihood integral?
State the two factors in the integrand and explain what each contributes.
Under what assumption can the population likelihood be written as a product of individual marginal likelihoods?
Why can repeated observations be conditionally independent but marginally dependent?
In the Gaussian model $Y\mid\eta\sim N(\eta,1)$, $\eta\sim N(0,1)$, what is the marginal density at $y=0$?

Solutions. (1) The integral is three-dimensional, because there are three random effects; 12 observations affect the conditional likelihood but not the integration dimension. (2) The factors are $f(y_i\mid\eta_i;\theta)$, the data likelihood conditional on the effects, and $g(\eta_i;\theta)$, the population density that weights effect values. (3) The subjects must be independent under the model, including independent random effects and subject-specific data generation. (4) A shared random effect links the observations; averaging over its uncertainty induces marginal dependence even if the conditional likelihood factors. (5) The marginal distribution is $N(0,2)$, so the density at zero is $1/\sqrt{4\pi}\approx 0.2821$.