Deriving the Multivariate Longitudinal Likelihood

Why Integrate Over Random Effects?

With multiple longitudinal markers, measurements from the same person are related. Some of that dependence can be represented by subject-specific random effects: conditional on those effects, the marker measurements follow a specified multivariate model. But a subject’s random effects are not observed. The likelihood for the observed data must therefore average the conditional density over the random-effects distribution. This integration produces the subject-specific marginal longitudinal density.

The calculation matters because the marginal likelihood is the quantity used to estimate the longitudinal model parameters from observed subjects. Replacing a random effect with an estimate, or multiplying separate marker likelihoods when the markers are dependent, does not in general give this likelihood. The result here concerns the longitudinal measurements only; it contains no event-time or hazard contribution.

Model and Assumptions

For subject \(i\), stack all observed marker measurements into a vector \(\mathbf y_i\) of length \(m_i\). The stacking order can be chosen for convenience, but it must be used consistently in the design matrices and covariance matrices. Let \(\mathbf X_i\boldsymbol\beta\) be the fixed-effects mean, \(\mathbf Z_i\mathbf b_i\) the random-effects contribution, and \(\mathbf R_i\) the within-subject residual covariance. The multivariate linear mixed model is

\mathbf y_i\mid \mathbf b_i \sim N_{m_i}\!\left(\mathbf X_i\boldsymbol\beta+\mathbf Z_i\mathbf b_i,\ \mathbf R_i\right), \qquad \mathbf b_i\sim N_q(\mathbf 0,\mathbf D).

Here \(\mathbf b_i\) has dimension \(q\), \(\mathbf D\) is its covariance matrix, and \(N_m(\boldsymbol\mu,\boldsymbol\Sigma)\) denotes an \(m\)-variate normal distribution with mean \(\boldsymbol\mu\) and covariance \(\boldsymbol\Sigma\). The design matrices may encode different marker-specific means, time effects, and random-effect loadings. The residual covariance \(\mathbf R_i\) can represent dependence among measurements, including dependence across markers within a visit.

Assumptions for the marginal density: Conditional on \(\mathbf b_i\) and the observed design, \(\mathbf y_i\) has the stated multivariate normal distribution. The random effects have the stated normal distribution, with \(\mathbf D\) positive definite. The residual covariance \(\mathbf R_i\) is positive definite, and the residual vector is independent of \(\mathbf b_i\). Subjects are independent under the model. The likelihood is conditional on the observed visit times and marker measurement pattern; if that pattern is informative in a way not accounted for by the model, it requires additional modeling.

These assumptions specify a joint distribution for all of subject \(i\)’s recorded markers. In particular, the measurements are not treated as independent merely because they belong to different markers. Earlier material on multivariate linear mixed models and within-visit residual covariance across markers supplies the model components; the focus here is the integration that turns the conditional model into a density for the observed vector.

Integrating the Random Effects

Write \(\mathbf r_i=\mathbf y_i-\mathbf X_i\boldsymbol\beta\), the observed vector centered at its fixed-effects mean. The density for subject \(i\), conditional on the design, is the integral of the conditional density times the random-effects density:

f_i(\mathbf y_i) = \int_{\mathbb R^q} f(\mathbf y_i\mid\mathbf b_i)\,f(\mathbf b_i)\,d\mathbf b_i.

To evaluate it, define \(\mathbf A_i=\mathbf D^{-1}+\mathbf Z_i^\mathsf{T}\mathbf R_i^{-1}\mathbf Z_i\) and \(\mathbf h_i=\mathbf Z_i^\mathsf{T}\mathbf R_i^{-1}\mathbf r_i\). The part of the exponent involving \(\mathbf b_i\) is a quadratic form. Expanding and completing the square gives

\begin{aligned} &\mathbf r_i^\mathsf{T}\mathbf R_i^{-1}\mathbf r_i -2\mathbf b_i^\mathsf{T}\mathbf h_i +\mathbf b_i^\mathsf{T}\mathbf A_i\mathbf b_i\\ &\quad = (\mathbf b_i-\mathbf A_i^{-1}\mathbf h_i)^\mathsf{T} \mathbf A_i(\mathbf b_i-\mathbf A_i^{-1}\mathbf h_i) +\mathbf r_i^\mathsf{T}\mathbf R_i^{-1}\mathbf r_i -\mathbf h_i^\mathsf{T}\mathbf A_i^{-1}\mathbf h_i. \end{aligned}

Because \(\mathbf A_i\) is positive definite, the integral of the first exponential term over \(\mathbf b_i\) is the normalizing integral of a \(q\)-variate Gaussian kernel, namely \((2\pi)^{q/2}|\mathbf A_i|^{-1/2}\). Combining it with the normalizing constants in the two densities yields

f_i(\mathbf y_i)= (2\pi)^{-m_i/2} |\mathbf R_i|^{-1/2}|\mathbf D|^{-1/2}|\mathbf A_i|^{-1/2} \exp\!\left\{-\frac12\left( \mathbf r_i^\mathsf{T}\mathbf R_i^{-1}\mathbf r_i -\mathbf h_i^\mathsf{T}\mathbf A_i^{-1}\mathbf h_i \right)\right\}.

This expression can also be written as a multivariate normal density. Define \(\mathbf V_i=\mathbf Z_i\mathbf D\mathbf Z_i^\mathsf{T}+\mathbf R_i\). The determinant identity \(|\mathbf V_i|=|\mathbf R_i||\mathbf D||\mathbf A_i|\) gives the normalizing factor \(|\mathbf V_i|^{-1/2}\). The inverse identity \(\mathbf V_i^{-1}=\mathbf R_i^{-1}-\mathbf R_i^{-1}\mathbf Z_i\mathbf A_i^{-1}\mathbf Z_i^\mathsf{T}\mathbf R_i^{-1}\) shows that the quadratic form in the exponent is \(\mathbf r_i^\mathsf{T}\mathbf V_i^{-1}\mathbf r_i\). Thus

Marginal longitudinal density: Under the stated Gaussian mixed-model assumptions, \[ \mathbf y_i\sim N_{m_i}\!\left(\mathbf X_i\boldsymbol\beta,\ \mathbf V_i\right), \qquad \mathbf V_i=\mathbf Z_i\mathbf D\mathbf Z_i^\mathsf{T}+\mathbf R_i, \] and \[ f_i(\mathbf y_i)= \frac{\exp\!\left\{-\tfrac12\mathbf r_i^\mathsf{T}\mathbf V_i^{-1}\mathbf r_i\right\}} {(2\pi)^{m_i/2}|\mathbf V_i|^{1/2}}. \]

The covariance has a direct interpretation: \(\mathbf Z_i\mathbf D\mathbf Z_i^\mathsf{T}\) is the covariance induced by variation in random effects, while \(\mathbf R_i\) is the residual covariance remaining conditional on those effects. Their sum is the covariance of the observed vector after the random effects have been integrated out.

Worked Example: A Two-Marker Marginal Density

Consider a hypothetical subject with one measurement on each of two markers. Suppose the random-effect covariance and residual covariance are \(\mathbf D=\begin{pmatrix}0.8&0.4\\0.4&0.8\end{pmatrix}\) and \(\mathbf R_i=\begin{pmatrix}0.7&0.6\\0.6&0.7\end{pmatrix}\), with \(\mathbf Z_i=\mathbf I_2\). Both covariance matrices are positive definite: their eigenvalues are \(0.4,1.2\) and \(0.1,1.3\), respectively. Take the fixed-effects mean to be \(\mathbf 0\) and the observed centered vector to be \(\mathbf r_i=(1,2)^\mathsf{T}\).

The marginal covariance is \(\mathbf V_i=\mathbf D+\mathbf R_i=\begin{pmatrix}1.5&1\\1&1.5\end{pmatrix}\). Its determinant is \(1.5(1.5)-1(1)=1.25\), and \(\mathbf V_i^{-1}=\begin{pmatrix}1.2&-0.8\\-0.8&1.2\end{pmatrix}\). The quadratic form counts both off-diagonal terms: \[ \mathbf r_i^\mathsf{T}\mathbf V_i^{-1}\mathbf r_i =1.2(1^2)+2(-0.8)(1)(2)+1.2(2^2) =1.2-3.2+4.8=2.8. \]

Therefore the density at this observed vector is \[ f_i(\mathbf y_i) =\frac{\exp(-2.8/2)}{(2\pi)^{2/2}\sqrt{1.25}} =\frac{\exp(-1.4)}{2\pi\sqrt{1.25}} \approx 0.03510. \] This is a density for the two-dimensional measurement vector, not the probability of observing exactly \((1,2)^\mathsf{T}\). Its units are the reciprocal of the product of the two measurement units.

Worked Example: Shared Random Intercept Across Two Markers

In another hypothetical model, one random intercept \(b_i\sim N(0,0.5)\) contributes equally to one measurement on each of two markers. Conditional on \(b_i\), the residuals are independent with variance \(0.5\), so \(\mathbf Z_i=(1,1)^\mathsf{T}\) and \(\mathbf R_i=0.5\mathbf I_2\). These assumptions give a positive-definite residual covariance and make the shared random intercept the source of cross-marker covariance.

The marginal covariance is \[ \mathbf V_i =0.5\begin{pmatrix}1&1\\1&1\end{pmatrix} +0.5\begin{pmatrix}1&0\\0&1\end{pmatrix} =\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}. \] For centered observations \(\mathbf r_i=(1,2)^\mathsf{T}\), \(|\mathbf V_i|=1-0.25=0.75\) and \(\mathbf V_i^{-1}=\frac{1}{0.75}\begin{pmatrix}1&-0.5\\-0.5&1\end{pmatrix}\). Multiplying gives \(\mathbf V_i^{-1}\mathbf r_i=(0,2)^\mathsf{T}\), so the quadratic form is \(1(0)+2(2)=4\). Hence \[ f_i(\mathbf y_i) =\frac{\exp(-4/2)}{2\pi\sqrt{0.75}} \approx 0.02487. \] The off-diagonal covariance \(0.5\) is present even though the residuals are independent: it arises from the shared random intercept.

From Subject Densities to the Longitudinal Likelihood

Given independent subjects, the likelihood for the longitudinal parameters is the product of their marginal densities. If \(\boldsymbol\theta\) collects the fixed effects and covariance parameters, then

L(\boldsymbol\theta)=\prod_{i=1}^{N} f_i(\mathbf y_i;\boldsymbol\theta), \qquad \ell(\boldsymbol\theta)=\log L(\boldsymbol\theta) =\sum_{i=1}^{N}\log f_i(\mathbf y_i;\boldsymbol\theta).

This product is over subjects, not over individual measurements treated as unrelated observations. The determinant and quadratic-form terms in each subject’s density account for the dependence among that subject’s measurements.

Worked Example: Combining Two Subject Contributions

Suppose two hypothetical independent subjects share the covariance from the first worked example, \(\mathbf V=\begin{pmatrix}1.5&1\\1&1.5\end{pmatrix}\), with fixed-effects means already subtracted. Their centered vectors are \(\mathbf r_1=(1,2)^\mathsf{T}\) and \(\mathbf r_2=(0,1)^\mathsf{T}\). The first quadratic form is \(2.8\), as calculated above. For the second, \[ \mathbf r_2^\mathsf{T}\mathbf V^{-1}\mathbf r_2 =1.2(0^2)+2(-0.8)(0)(1)+1.2(1^2)=1.2. \]

Each subject has determinant \(1.25\), so their joint likelihood contribution is \[ L=f_1f_2 =\frac{\exp[-(2.8+1.2)/2]}{(2\pi)^2(1.25)} =\frac{\exp(-2)}{49.348} \approx 0.002743. \] Equivalently, the log likelihood contribution is \[ \ell=\log L=-2-\log\!\left((2\pi)^2(1.25)\right)\approx -5.898. \] This calculation combines the two subject-level densities under the independence assumption; it does not multiply the markers within a subject as if their covariance were zero.

Interpretation, Diagnostics, and Alternatives

The marginal form is useful because it gives a density for the observed measurements alone, with random effects integrated out rather than observed or fixed at estimated values. It also makes clear how random-effect covariance and residual covariance jointly determine dependence. The same subject may have measurements at unequal times or a different number of recorded measurements: the dimensions of \(\mathbf y_i\), \(\mathbf Z_i\), and \(\mathbf R_i\) can vary across subjects, provided each subject’s model is constructed consistently.

A direct diagnostic is to check that each covariance matrix is symmetric and positive definite and that its dimensions match the stacked data. Numerical calculations should also be checked for positive determinants and sensible quadratic forms: for positive-definite \(\mathbf V_i\), \(\mathbf r_i^\mathsf{T}\mathbf V_i^{-1}\mathbf r_i\) cannot be negative. When the marginal covariance is built from blocks, verify that the blocks follow the stated stacking order. A mistaken ordering can create apparently plausible but incorrect cross-marker covariances.

Two common alternatives are useful computationally, but they do not change the target likelihood. One can evaluate the integral over random effects directly, or use the equivalent marginal multivariate normal density involving \(\mathbf V_i\). The latter is convenient when the marginal covariance can be formed and factored efficiently; the completed-square form can be useful when the random-effects dimension is smaller than the number of measurements. In contrast, evaluating the conditional density at a fitted or predicted random effect is not integration and does not give \(f_i(\mathbf y_i)\). Likewise, multiplying marker-specific marginal densities is valid only when the model makes the relevant marker vectors independent; shared random effects or nonzero cross-marker residual covariance generally invalidate that factorization.

The result depends on the Gaussian assumptions. If the random-effects distribution or conditional measurement distribution is non-Gaussian, the integrated density need not have this multivariate normal form, even though the integral definition remains meaningful. And if the schedule or presence of measurements depends informatively on unmodeled processes, conditioning on the observed pattern alone may not describe the full data-generating mechanism. These are limits of the stated likelihood, not reasons to alter the Gaussian integration.

Common mistake: counting an off-diagonal term once. In a quadratic form \(\mathbf r^\mathsf{T}\mathbf V^{-1}\mathbf r\), each off-diagonal element contributes twice. For example, the cross-term for a two-dimensional vector is \(2(\mathbf V^{-1})_{12}r_1r_2\), not \((\mathbf V^{-1})_{12}r_1r_2\).
Correction: Expand the full matrix product, or multiply \(\mathbf V^{-1}\mathbf r\) first and then take its inner product with \(\mathbf r\). Check that the result is nonnegative.
Common mistake: treating a plug-in random effect as the marginal model. A conditional density evaluated at one estimated \(\mathbf b_i\) omits uncertainty and integration over the random-effects distribution.
Correction: Use the integral over \(\mathbf b_i\), or the equivalent marginal density with covariance \(\mathbf Z_i\mathbf D\mathbf Z_i^\mathsf{T}+\mathbf R_i\).

Key takeaway: Integrating a Gaussian multivariate mixed model over its Gaussian random effects gives a multivariate normal density for each subject’s stacked marker vector. Its mean is the fixed-effects mean, its covariance is the sum of random-effect and residual covariance contributions, and independent subject densities multiply to form the longitudinal likelihood.

Check Your Understanding

Answer each question, then compare with the brief solutions.

If \(\mathbf Z_i=\mathbf I\), \(\mathbf D=\begin{pmatrix}1&0.3\\0.3&1\end{pmatrix}\), and \(\mathbf R_i=\begin{pmatrix}0.5&0\\0&0.5\end{pmatrix}\), what is the marginal covariance?
For a two-dimensional residual vector \((1,2)^\mathsf{T}\) and inverse covariance \(\begin{pmatrix}1.2&-0.8\\-0.8&1.2\end{pmatrix}\), calculate the quadratic form.
Why is the likelihood for independent subjects a product of subject densities rather than a product of independent measurement densities?
Does evaluating \(f(\mathbf y_i\mid\mathbf b_i)\) at a predicted random effect produce the marginal density? Explain.
Name two assumptions that permit the closed-form multivariate normal marginal density.

Solutions. (1) Add the matrices: \(\mathbf V_i=\begin{pmatrix}1.5&0.3\\0.3&1.5\end{pmatrix}\). (2) \(1.2+2(-0.8)(1)(2)+1.2(2^2)=1.2-3.2+4.8=2.8\). (3) Measurements within a subject may be dependent through random effects or residual covariance; the subject density accounts for that dependence, while independence is assumed across subjects. (4) No. It is a conditional density at one value of the random effect; the marginal density averages the conditional density over the random-effects distribution. (5) For example, the conditional measurement distribution and random-effects distribution are Gaussian, and the residual vector is independent of the random effects. The stated covariance matrices must also be positive definite for the displayed density to apply.