Learn to write and simplify the joint log-likelihood of observations and latent Gaussian random effects while retaining every term that depends on population parameters.
Advanced
11 min read
What You'll Learn
Define the complete data and its joint density in an NLME model
Factor the likelihood across independent subjects
Derive the Gaussian random-effect log-density and identify its parameter-dependent terms
Construct complete-data likelihoods for Gaussian and Poisson observation models
Distinguish removable constants from terms needed for population-parameter estimation
Recognize common errors in forming or simplifying the complete-data log-likelihood
Why the Complete-Data Log-Likelihood Is a Joint Density
In a nonlinear mixed-effects model, each subject has observations and an unobserved subject-specific random effect. The complete-data log-likelihood treats both as if they were available: it is the log of the joint density of the observations and the random effects. This is different from the likelihood based on observations alone, which requires integrating out the random effects. Writing the joint form correctly is essential because it makes explicit which terms involve the population parameters.
Let \(i=1,\ldots,N\) index independent subjects. Write \(y_i\) for the vector of observations on subject \(i\), and \(\eta_i\) for that subject’s \(q\)-dimensional random-effect vector. Let \(\beta\) denote fixed effects, \(\Omega\) the random-effect covariance matrix, and \(\phi\) any additional observation-model parameters, such as residual variance parameters. Collect the population parameters as \(\theta=(\beta,\phi,\Omega)\). The conditional observation density is \(p(y_i\mid\eta_i;\beta,\phi)\), and the random effects have distribution \(\eta_i\sim N(0,\Omega)\).
Definition: The complete data are \((y,\eta)=\{(y_i,\eta_i):i=1,\ldots,N\}\). Under independent subjects, the complete-data likelihood and log-likelihood are
$$
L_c(\theta;y,\eta)=\prod_{i=1}^N p(y_i\mid\eta_i;\beta,\phi)\,p(\eta_i;\Omega),
\qquad
\ell_c(\theta;y,\eta)=\sum_{i=1}^N\left\{\log p(y_i\mid\eta_i;\beta,\phi)+\log p(\eta_i;\Omega)\right\}.
$$
Here, \(\eta\) is part of the complete-data argument, not an additional population parameter. The model assumes subjects are independent and that each subject’s random effect has the specified Gaussian distribution. The conditional density \(p(y_i\mid\eta_i;\beta,\phi)\) must describe the observation process correctly; it may factor across records within a subject when those records are conditionally independent, or it may be a joint density that includes within-subject dependence. Do not assume record-level independence unless the model specifies it.
The marginal likelihood is instead \(L(\theta;y)=\prod_i\int p(y_i\mid\eta_i;\beta,\phi)p(\eta_i;\Omega)\,d\eta_i\). That integration is precisely what makes the observed-data likelihood different. The complete-data expression keeps the random effects explicit; it does not replace them with fitted values or condition on their estimated modes.
The Gaussian Random-Effect Contribution
For a \(q\)-dimensional Gaussian random effect with positive-definite covariance matrix \(\Omega\), the density is
$$
p(\eta_i;\Omega)=(2\pi)^{-q/2}|\Omega|^{-1/2}\exp\left(-\frac{1}{2}\eta_i^\mathsf{T}\Omega^{-1}\eta_i\right).
$$
Taking the logarithm gives a contribution for subject \(i\) of
$$
\log p(\eta_i;\Omega)=-\frac{q}{2}\log(2\pi)-\frac{1}{2}\log|\Omega|-\frac{1}{2}\eta_i^\mathsf{T}\Omega^{-1}\eta_i.
$$
Summing over subjects, the complete-data log-likelihood can therefore be written
$$
\ell_c(\theta;y,\eta)=\sum_{i=1}^N\log p(y_i\mid\eta_i;\beta,\phi)-\frac{Nq}{2}\log(2\pi)-\frac{N}{2}\log|\Omega|-\frac{1}{2}\sum_{i=1}^N\eta_i^\mathsf{T}\Omega^{-1}\eta_i.
$$
This expression separates the observation-model contribution from the Gaussian random-effect contribution. The latter depends on \(\Omega\) through both \(\log|\Omega|\) and \(\Omega^{-1}\). Both are essential when estimating the covariance. The observation contribution depends on \(\beta\), \(\phi\), and the random effects wherever they enter the conditional model.
A term may be discarded from an objective only if it is constant with respect to every parameter being optimized. For example, \(-Nq\log(2\pi)/2\) is constant when \(N\) and \(q\) are fixed. By contrast, \(-N\log|\Omega|/2\) is not constant when \(\Omega\) is estimated. A data-only term such as \(-\log(y_{ij}!)\) may be constant with respect to the model parameters, but it remains part of the exact log-density.
Deriving the Form and Checking Parameter Dependence
The factorization follows from the joint-density rule and the model’s independence assumptions: for each subject, the joint density of observations and random effect is \(p(y_i,\eta_i)=p(y_i\mid\eta_i;\beta,\phi)p(\eta_i;\Omega)\); independence across subjects makes the full density the product of these terms. Taking a logarithm converts the product into a sum. The Gaussian contribution follows by taking the logarithm of its normalized density, so the determinant and quadratic terms cannot be omitted without changing the parameter-dependent objective.
1
Write the subject-level joint density. Use \(p(y_i,\eta_i)=p(y_i\mid\eta_i;\beta,\phi)p(\eta_i;\Omega)\), not just the conditional observation density.
2
Multiply across independent subjects. The model’s subject-level independence gives a product over \(i=1,\ldots,N\).
3
Take logarithms and expand the Gaussian term. Retain every term involving a parameter being estimated, including the covariance determinant and quadratic form.
4
Mark parameter-free constants only after expansion. Whether a term can be omitted depends on the optimization target, not on whether it looks numerically constant at the current parameter values.
Worked Example: Gaussian Random-Intercept Model
Consider a hypothetical model with two independent subjects and two observations per subject:
\(y_{ij}=\beta+\eta_i+\epsilon_{ij}\), where \(\eta_i\sim N(0,\omega^2)\), \(\epsilon_{ij}\sim N(0,\sigma^2)\), and the residuals are conditionally independent given \(\eta_i\). Suppose the observations are \(y_1=(10,12)\) and \(y_2=(8,9)\). To evaluate the complete-data log-likelihood at a particular point, use \(\beta=9\), \(\sigma^2=4\), \(\omega^2=1\), \(\eta_1=1\), and \(\eta_2=-1\). These are hypothetical values supplied for evaluating the expression, not estimates from a study.
For each observation, the conditional residual is \(y_{ij}-(\beta+\eta_i)\). The four residuals are \(0,2,0,1\), so their squared sum is \(0^2+2^2+0^2+1^2=5\). With \(M=4\) total observations, the conditional Gaussian contribution is
$$
-\frac{M}{2}\log(2\pi\sigma^2)-\frac{5}{2\sigma^2}
=-2\log(8\pi)-\frac{5}{8}
\approx -7.073.
$$
The random-effect contribution is
$$
-\frac{N}{2}\log(2\pi\omega^2)-\frac{\eta_1^2+\eta_2^2}{2\omega^2}
=-\log(2\pi)-\frac{1^2+(-1)^2}{2}
\approx -2.838.
$$
Thus \(\ell_c\approx -7.073-2.838=-9.911\). The residual term involves \(\beta\) and \(\sigma^2\) through the residuals and variance; the random-effect term involves \(\omega^2\) through its logarithm and the scaled squared effects. The values \(\eta_1,\eta_2\) are held fixed for this evaluation of the complete-data function.
Worked Example: A Correlated Two-Dimensional Random Effect
Suppose a hypothetical subject has a two-dimensional random effect \(\eta=(0.5,-1)^\mathsf{T}\) and a candidate covariance matrix
$$
\Omega=\begin{pmatrix}4&1\\1&1\end{pmatrix}.
$$
The determinant is \(4(1)-1(1)=3\), and the inverse is
$$
\Omega^{-1}=\frac{1}{3}\begin{pmatrix}1&-1\\-1&4\end{pmatrix}.
$$
The quadratic form is
$$
\eta^\mathsf{T}\Omega^{-1}\eta
=(0.5,-1)\begin{pmatrix}0.5\\-1.5\end{pmatrix}
=0.25+1.5=1.75.
$$
For one subject and \(q=2\), the random-effect log-density is therefore
$$
-\log(2\pi)-\frac{1}{2}\log(3)-\frac{1}{2}(1.75)
\approx -3.262.
$$
If this subject’s conditional observation log-density at the same parameter values is \(-2.400\), the complete-data log-likelihood contribution is \(-2.400-3.262=-5.662\). In an optimization over \(\Omega\), neither the log-determinant nor the quadratic form is a constant: both change as the covariance changes. Keeping only the quadratic penalty would therefore give the wrong covariance objective.
Worked Example: Poisson Observations and a Gaussian Random Effect
As a different hypothetical observation model, let one subject have counts \(y=(2,0)\), with conditionally independent observations \(Y_j\mid\eta\sim\operatorname{Poisson}(\lambda)\) and \(\lambda=\exp(\beta+\eta)\). Let \(\eta\sim N(0,\omega^2)\). Evaluate the complete-data log-likelihood at \(\beta=0\), \(\eta=\log 2\), and \(\omega^2=0.25\). Then \(\lambda=\exp(\log 2)=2\).
The conditional log-likelihood is \(\sum_j\{y_j\log\lambda-\lambda-\log(y_j!)\}\). Substitution gives
$$
[2\log(2)-2-\log(2!)]+[0\log(2)-2-\log(0!)]
=2\log(2)-4-\log(2)
\approx -3.307.
$$
Here \(\log(2!)=\log 2\) and \(\log(0!)=0\). The Gaussian prior contribution is
$$
-\frac{1}{2}\log(2\pi\omega^2)-\frac{\eta^2}{2\omega^2}
=-\frac{1}{2}\log(\pi/2)-\frac{(\log 2)^2}{0.5}
\approx -1.187.
$$
The resulting complete-data log-likelihood is approximately \(-3.307-1.187=-4.494\). The factorial terms are constant with respect to \(\beta\) and \(\omega^2\) for these fixed counts, but the Poisson rate terms and both Gaussian random-effect terms must be retained when their corresponding parameters are optimized.
Diagnostics, Alternatives, and Common Mistakes
Several checks help catch errors when constructing the expression. First, the Gaussian covariance contribution should contain a determinant term and a quadratic form. Second, the observation contribution should be the correct normalized conditional density; for Gaussian observations, check that the variance normalizer is included if that variance is estimated. Third, verify that all subject-specific terms are summed over the correct subjects and that the random-effect dimension matches the dimension used in the Gaussian normalizing constant. Finally, classify terms relative to the parameters being optimized. A term may be constant for one optimization but not another.
The complete-data likelihood and the marginal likelihood serve different roles. The complete-data expression is a joint density evaluated at \(y\) and \(\eta\); the marginal likelihood is based only on \(y\), with each \(\eta_i\) integrated out. A conditional mode or point estimate for \(\eta_i\) is not a substitute for this integration. In nonlinear models, the integral generally cannot be replaced by evaluating the integrand at one random-effect value.
A frequent mistake is to write only \(\sum_i\log p(y_i\mid\eta_i;\beta,\phi)\), reasoning that the random effects are already included in the complete data. That omits their model density. The correction is to include \(\sum_i\log p(\eta_i;\Omega)\) as well: the random effects are latent variables in the complete data, and their distribution is part of the joint model. Another mistake is to discard \(-N\log|\Omega|/2\) as a “normalizing constant.” It is not constant in an estimation of \(\Omega\).
The form above assumes a zero-mean Gaussian random-effect distribution with common covariance \(\Omega\). If the model instead specifies a nonzero random-effect mean or a different covariance structure, its density must be written using that specified distribution; the zero-mean formula cannot be reused unchanged. Likewise, dependence between subjects requires a joint random-effect density rather than a sum of independent subject contributions. The introductory material on the EM view of latent random effects motivates treating these effects as missing data; here the focus is the exact joint log-likelihood they produce.
Key takeaway: Write the complete-data log-likelihood as the sum of conditional observation log-densities and random-effect log-densities. Keep every term that depends on a population parameter, especially the Gaussian covariance log-determinant and quadratic form.
Check Your Understanding
Try each question, then compare your reasoning with the solutions below.
For independent subjects with \(p(y_i\mid\eta_i;\beta,\phi)\) and \(\eta_i\sim N(0,\Omega)\), write the complete-data log-likelihood.
Why can the Gaussian random-effect term \(-N\log|\Omega|/2\) not be dropped when estimating \(\Omega\)?
For a scalar random effect with variance \(\omega^2\), identify the two parts of its summed log-density that depend on \(\omega^2\).
When is a term such as \(-\log(y!)\) removable from an optimization objective, and what should be remembered about it?
Solutions. (1) It is \(\sum_i\{\log p(y_i\mid\eta_i;\beta,\phi)+\log p(\eta_i;\Omega)\}\), because subject independence gives a product of joint densities and the logarithm turns that product into a sum. (2) The determinant changes with \(\Omega\), so dropping it changes the objective being optimized. (3) The terms are \(-N\log(\omega^2)/2\) and \(-\sum_i\eta_i^2/(2\omega^2)\), in addition to constants independent of \(\omega^2\). (4) For fixed observed \(y\), \(-\log(y!)\) is removable if optimizing only parameters that do not enter it; it remains part of the exact log-density and should not be removed when an exact likelihood value is needed.