Why the Likelihood Is Marginal
A linear mixed model describes repeated biomarker measurements by combining a population-level mean trajectory with subject-level random effects and residual error. For likelihood-based inference, the observed data are the biomarker records, not the unobserved random effects. The likelihood for those records must therefore account for the random effects without treating them as known.
The key step is to integrate over the distribution of the random effects. Under Gaussian assumptions, that integration gives a multivariate normal distribution for each subject’s repeated measurements. Its mean comes from the fixed effects, while its covariance combines random-effects variation and residual variation.
Here \(D\) is the random-effects covariance matrix and \(R_i\) is the residual covariance matrix for subject \(i\). These matrices may depend on a vector of covariance parameters \(\theta\). The design matrices are treated as observed and fixed for the likelihood calculation. The assumptions include Gaussian random effects and Gaussian residual errors, as well as independence between the two sources of variation. The marginal covariance \(V_i\) must be positive definite for the usual nonsingular Gaussian density.
From the Model to the Marginal Likelihood
Conditional on \(b_i\), the observation vector has distribution \(y_i\mid b_i\sim N(X_i\beta+Z_i b_i,R_i)\). But \(b_i\) is not observed. The marginal density is obtained by averaging this conditional density over the random-effects distribution:
To see what this integration produces, note that \(Z_i b_i\) and \(\varepsilon_i\) are independent Gaussian vectors. Their sum is Gaussian. Its mean is zero, and its covariance is the sum of the two covariances: $$ \operatorname{Var}(Z_i b_i+\varepsilon_i) =Z_i\operatorname{Var}(b_i)Z_i^\mathsf{T}+\operatorname{Var}(\varepsilon_i) =Z_i D Z_i^\mathsf{T}+R_i. $$ Adding the fixed vector \(X_i\beta\) gives \(y_i\sim N(X_i\beta,V_i)\). Thus the integral is the multivariate normal density with that mean and covariance.
For \(N\) independent subjects, the observed-data likelihood is the product of their marginal densities. Equivalently, its log is a sum:
The determinant term \(\log|V_i|\) reflects the spread of the distribution: a more dispersed covariance spreads probability over a larger region. The quadratic term measures the subject’s residual vector \(y_i-X_i\beta\) relative to that covariance. A residual in a direction with high modeled variance contributes less to the quadratic term than an equally large residual in a direction with low variance.
Which Parameters Are Estimated?
The likelihood depends on two broad groups of unknowns. The fixed-effects parameters \(\beta\) determine the population-average mean at each subject’s observed times. The covariance parameters \(\theta\) determine the structure and size of subject-to-subject variation and residual variation through \(D\), \(R_i\), and therefore \(V_i\).
| Likelihood component | Typical parameters | Role in the model |
|---|---|---|
| Fixed-effects mean | Intercept, time coefficient, treatment coefficient, treatment-by-time coefficient | Determine \(X_i\beta\), the mean biomarker trajectory |
| Random-effects covariance | Random-intercept variance, random-slope variance, their covariance | Describe variation among subjects’ trajectories |
| Residual covariance | Residual variance and any parameters in \(R_i\) | Describe remaining within-subject variation |
The exact parameter list depends on the model specified. For example, if residual errors are independent with a common variance, then \(R_i=\sigma^2 I_{m_i}\), and \(\sigma^2\) is a covariance parameter. If \(D\) and \(R_i\) are specified as known, they are not estimated. Observed times and design-matrix entries are inputs to the likelihood, not parameters estimated from it.
Worked Example: A Random-Intercept Likelihood Contribution
Suppose one hypothetical subject has biomarker measurements \(y_i=(10,12)^\mathsf{T}\) at two visits. At these visits, the fixed-effects model gives means of 11 units, so \(X_i\beta=(11,11)^\mathsf{T}\). Consider a random-intercept model with random-intercept variance 4 and independent residual variance 1 at each visit. Calculate the subject’s marginal likelihood contribution at these parameter values.
For a random intercept, \(Z_i=(1,1)^\mathsf{T}\), \(D=(4)\), and \(R_i=I_2\). Hence $$ V_i=Z_i D Z_i^\mathsf{T}+R_i =\begin{pmatrix}4&4\\4&4\end{pmatrix} +\begin{pmatrix}1&0\\0&1\end{pmatrix} =\begin{pmatrix}5&4\\4&5\end{pmatrix}. $$ The determinant is \(5(5)-4(4)=9\), and $$ V_i^{-1}=\frac{1}{9}\begin{pmatrix}5&-4\\-4&5\end{pmatrix}. $$ The residual vector is \(r_i=y_i-X_i\beta=(-1,1)^\mathsf{T}\). Its quadratic form is $$ r_i^\mathsf{T}V_i^{-1}r_i =(-1,1)\frac{1}{9}\begin{pmatrix}5&-4\\-4&5\end{pmatrix} \begin{pmatrix}-1\\1\end{pmatrix} =2. $$ Since \(m_i=2\), the density contribution is $$ (2\pi)^{-1}|V_i|^{-1/2}\exp(-2/2) =\frac{1}{6\pi}e^{-1} \approx 0.01949. $$ This is a density value for a pair of measurements, not a probability that the exact pair occurs. The calculation uses the marginal covariance, so the shared random intercept induces covariance 4 between the two visits.
Worked Example: Combining Independent Subjects
Consider two hypothetical subjects, each measured twice. Let both subjects have fitted mean vector \((10,10)^\mathsf{T}\). Use a random-intercept variance of 1 and an independent residual variance of 1. Subject 1 has measurements \((8,10)^\mathsf{T}\), and subject 2 has measurements \((12,10)^\mathsf{T}\). Calculate the joint likelihood contribution for these subjects.
For each subject, \(V_i=\begin{pmatrix}2&1\\1&2\end{pmatrix}\), so \(|V_i|=3\) and \(V_i^{-1}=\frac{1}{3}\begin{pmatrix}2&-1\\-1&2\end{pmatrix}\). Subject 1 has residual \((-2,0)^\mathsf{T}\), giving quadratic form $$ (-2,0)V_i^{-1}(-2,0)^\mathsf{T}=\frac{8}{3}. $$ Subject 2 has residual \((2,0)^\mathsf{T}\), which gives the same quadratic form, \(8/3\). Therefore each density contribution is $$ (2\pi)^{-1}3^{-1/2}\exp\!\left(-\frac{1}{2}\cdot\frac{8}{3}\right) =\frac{1}{2\pi\sqrt{3}}e^{-4/3} \approx 0.02422. $$ Independence across subjects means their contributions multiply, giving a joint likelihood of approximately \(0.02422^2=0.0005866\). The multiplication is across subjects; the two measurements within each subject are handled together using that subject’s covariance matrix.
Random Slopes and the Covariance They Imply
The same likelihood formula applies when subjects have random slopes as well as random intercepts. In that case, \(Z_i\) has a column for each random effect, and the elements of \(D\) determine how the random intercepts and slopes vary and co-vary. The following calculation shows how those parameters enter a subject’s marginal likelihood.
Worked Example: A Random-Intercept-and-Slope Contribution
Suppose a hypothetical subject is measured at times 0 and 2 years, with observed values \((9,13)^\mathsf{T}\). The fixed-effects mean is 10 units at year 0 and 12 units at year 2. Use random intercept and slope effects with \(D=\begin{pmatrix}4&0.5\\0.5&1\end{pmatrix}\), and independent residual variance 1 at each visit. Calculate the marginal density contribution.
The random-effects design matrix is \(Z_i=\begin{pmatrix}1&0\\1&2\end{pmatrix}\). The covariance matrix \(D\) is valid because its first diagonal entry is positive and its determinant is \(4(1)-(0.5)^2=3.75>0\). The marginal covariance is $$ V_i=Z_i D Z_i^\mathsf{T}+I_2 =\begin{pmatrix}5&5\\5&11\end{pmatrix}, \qquad |V_i|=5(11)-5(5)=30. $$ Thus \(V_i^{-1}=\frac{1}{30}\begin{pmatrix}11&-5\\-5&5\end{pmatrix}\). The residual vector is \((-1,1)^\mathsf{T}\), and $$ (-1,1)V_i^{-1}(-1,1)^\mathsf{T} =\frac{26}{30} =\frac{13}{15} \approx 0.8667. $$ The likelihood contribution is therefore $$ \frac{1}{2\pi\sqrt{30}}\exp\!\left(-\frac{13}{30}\right) \approx 0.01884. $$ The off-diagonal covariance of 5 follows from the shared intercept plus the intercept-slope covariance multiplied by the later time: \(4+2(0.5)=5\). The later-visit variance also includes slope variation, which is why it is larger than the year-0 variance.
Likelihood Checks, Alternatives, and Limitations
The marginal likelihood is useful because it evaluates how well the model’s mean and covariance jointly account for the observed repeated records. It also makes explicit that the number of visits may differ across subjects: each subject contributes a density of dimension \(m_i\), and the contributions are multiplied under subject independence.
A conditional model and a marginal likelihood are related but are not interchangeable calculations. The conditional density \(f(y_i\mid b_i,\beta,\theta)\) describes measurements for a specified random effect. The observed-data likelihood instead integrates that density over \(f(b_i\mid\theta)\). Substituting a predicted or estimated \(b_i\) into the conditional density does not perform this integration and does not produce the marginal likelihood.
Before using this likelihood, check that the model’s covariance matrices are meaningful and that the assumptions are plausible for the data at hand. In particular, subject-level independence and the chosen within-subject covariance structure matter: an omitted dependence pattern can make \(V_i\) poorly represent the repeated measurements. A Gaussian likelihood can also be sensitive to unusual observations because the quadratic term grows with the squared, covariance-scaled residual.
The formula is conditional on the observed design and visit times. If the timing or occurrence of visits is informative about the biomarker process, this likelihood alone may not describe the full observed-data mechanism; the visit process may require separate consideration. Also, likelihood values are densities whose numerical scale depends on the measurement units, so comparing raw values from differently scaled outcomes is generally not meaningful.
A common mistake is to treat all rows in the long-format data as independent Gaussian observations. That uses only separate one-measurement contributions and discards the off-diagonal elements of \(V_i\). The correction is to form one covariance matrix for each subject’s observed measurement vector and use its determinant and quadratic form in that subject’s likelihood contribution. Another mistake is to call \(L(\beta,\theta)\) the probability that the parameters are correct. It is instead the density of the observed data as a function of candidate parameter values.
In practice, the likelihood is evaluated over candidate values of the mean and covariance parameters. Numerical difficulties can arise if a proposed covariance matrix is not positive definite, if the data contain little information about a variance component, or if the assumed covariance structure does not capture the repeated-measurement pattern. The likelihood formula identifies what is being evaluated; it does not by itself guarantee that all parameters are well determined or that the model is adequate.
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- For a random-intercept model with random-intercept variance \(\tau^2\) and independent residual variance \(\sigma^2\), what are the marginal variance and covariance of two measurements from the same subject?
- Suppose \(y_i\) has four observed measurements. What dimension is \(V_i\), and what is the coefficient of \(\log|V_i|\) in that subject’s log-likelihood contribution?
- Why is evaluating the conditional density at a predicted random intercept not equivalent to calculating the marginal likelihood?
- For one subject, suppose \(V_i=\begin{pmatrix}3&1\\1&3\end{pmatrix}\) and the residual vector is \((1,1)^\mathsf{T}\). Calculate the quadratic form.
Solutions. (1) Each marginal variance is \(\tau^2+\sigma^2\), and the covariance between the two measurements is \(\tau^2\), because they share the same random intercept. (2) \(V_i\) is \(4\times4\); the coefficient of \(\log|V_i|\) is \(-1/2\). (3) The marginal likelihood averages the conditional density over the full random-effects distribution by integration; a single predicted value does not do that. (4) The determinant is 8 and the inverse is \(\frac{1}{8}\begin{pmatrix}3&-1\\-1&3\end{pmatrix}\). Therefore the quadratic form is \((1,1)V_i^{-1}(1,1)^\mathsf{T}=4/8=0.5\).