Why the Random-Effects Distribution Matters
In a shared random-effects joint model, each subject has latent effects that connect their longitudinal trajectory and event-time process. The conditional subject contribution describes the observed data for a specified value of those effects. The marginal likelihood must also account for how plausible that value is under the random-effects distribution. Thus, the distribution is not an optional description attached to the model: its density weights the conditional contribution at every possible random-effects value.
This tutorial focuses on that weighting, especially when a subject has both a random intercept and a random slope. The intercept and slope may vary together across subjects. Their covariance determines the orientation and shape of the random-effects distribution, and therefore influences both the likelihood integral and the interpretation of subject-level heterogeneity.
The mean is conventionally set to zero because population-level means are represented by fixed effects; \(b_{0i}\) and \(b_{1i}\) are deviations from those means. The random-effects vectors are modeled as independent across subjects with the same distribution \(N_2(0,D)\). The matrix \(D\) must be positive definite: in this two-dimensional case, \(\sigma_0^2>0\), \(\sigma_1^2>0\), and \(\sigma_{01}^2<\sigma_0^2\sigma_1^2\). These conditions ensure a valid, nondegenerate density. The distribution is a model assumption, not a claim that random effects are directly observed or that their empirical distribution must look exactly normal in a finite sample.
The Density and Its Role in the Likelihood
For a positive-definite \(q\times q\) covariance matrix \(D\), the density of \(b_i\sim N_q(0,D)\), evaluated at a vector \(b\), is
The determinant \(|D|\) controls the density’s overall scale, while the quadratic form \(b^\mathsf{T}D^{-1}b\) measures how far \(b\) lies from the mean relative to the covariance structure. A value that is far from zero in ordinary Euclidean distance need not be equally unusual in every direction: the covariance determines which directions have greater or lesser spread.
Let \(Y_i\) denote the longitudinal data for subject \(i\), and let \(T_i\) and \(\delta_i\) denote that subject’s event or censoring time and event indicator. Write \(f_i(Y_i,T_i,\delta_i\mid b_i;\theta)\) for the conditional subject likelihood contribution under the longitudinal and event submodels, with model parameters \(\theta\). Integrating over the unobserved random effects gives the marginal contribution:
For independent subjects, the full likelihood is the product of these marginal contributions over \(i\). The integral is an average of the conditional likelihood over the assumed random-effects distribution. In particular, it is not generally the conditional likelihood evaluated at a single estimated random effect. The earlier material on marginal likelihood over shared random effects establishes this integration step; here the emphasis is on the normal density that supplies its weights.
Where the Multivariate Normal Density Comes From
A short change-of-variables argument gives the density formula. Since \(D\) is positive definite, write \(D=AA^\mathsf{T}\) for a nonsingular matrix \(A\), and let \(Z\sim N_q(0,I)\). Then \(b=AZ\) has covariance \(D\). The standard normal density of \(Z\) is \((2\pi)^{-q/2}\exp(-Z^\mathsf{T}Z/2)\). Under the transformation \(Z=A^{-1}b\), the absolute Jacobian determinant for the density transformation is \(|A|^{-1}\). Also, \(Z^\mathsf{T}Z=b^\mathsf{T}(AA^\mathsf{T})^{-1}b=b^\mathsf{T}D^{-1}b\), and \(|A|=|D|^{1/2}\). Substitution gives the stated density.
This derivation also explains why both the determinant and the inverse covariance matrix appear. The determinant accounts for the volume change from standard-normal coordinates; the inverse covariance scales deviations according to the model’s spread and dependence.
Worked Example: Evaluating a Bivariate Random-Effects Density
Suppose the random intercept and slope have covariance matrix \(D=\begin{pmatrix}0.20&0.10\\0.10&0.40\end{pmatrix}\). Evaluate the density at \(b=(0.30,-0.20)^\mathsf{T}\). These values and this model are hypothetical.
First, \(|D|=(0.20)(0.40)-(0.10)^2=0.0700\), so \(D^{-1}=\frac{1}{0.07}\begin{pmatrix}0.40&-0.10\\-0.10&0.20\end{pmatrix}\). The quadratic form is \(b^\mathsf{T}D^{-1}b=\frac{1}{0.07}\{(0.40)(0.30)^2-2(0.10)(0.30)(-0.20)+(0.20)(-0.20)^2\}\). The expression in braces is \(0.036+0.012+0.008=0.056\), giving a quadratic form of \(0.800\).
For \(q=2\), the density is therefore \(\frac{\exp(-0.800/2)}{2\pi\sqrt{0.0700}}=\frac{\exp(-0.400)}{1.66237}\approx 0.4032\). This is a density value at the stated random-effects vector, not a probability that the vector equals that exact point.
Correlated Intercepts and Slopes
The covariance \(\sigma_{01}\) describes how subjects’ intercept deviations and slope deviations vary together. Its scale depends on the units of each effect, so the correlation is often easier to compare:
For example, a positive covariance means that subjects with intercepts above the population-level intercept tend, on average, also to have slopes above the population-level slope. It does not mean that every such subject has a positive slope, nor does it establish a causal relationship. The interpretation of the intercept also depends on the time origin: with time centered at a meaningful value, the random intercept is the deviation at that time.
Allowing a nonzero covariance changes the shape of the joint density. If random intercept and slope are assumed independent, \(D\) is diagonal and the density factors into separate univariate normal densities. That restriction can be useful when justified or needed for stability, but it is not implied merely because the two effects are listed separately in the model.
Worked Example: How Correlation Changes the Density
Consider \(D=\begin{pmatrix}0.25&0.15\\0.15&0.25\end{pmatrix}\), so both random effects have standard deviation \(0.50\) and their correlation is \(0.15/(0.50\cdot0.50)=0.600\). At \(b=(0.40,0.20)^\mathsf{T}\), the determinant is \(0.25(0.25)-0.15^2=0.0400\), and \(D^{-1}=\frac{1}{0.04}\begin{pmatrix}0.25&-0.15\\-0.15&0.25\end{pmatrix}\).
The quadratic form is \(\frac{1}{0.04}\{0.25(0.40)^2-2(0.15)(0.40)(0.20)+0.25(0.20)^2\} =\frac{0.040-0.024+0.010}{0.04}=0.650\). Thus the density is \(\exp(-0.650/2)/(2\pi\sqrt{0.0400})=\exp(-0.325)/1.25664\approx0.5750\).
If the covariance were instead set to zero while keeping both variances at \(0.25\), the quadratic form would be \(0.40^2/0.25+0.20^2/0.25=0.800\), and the density would be \(\exp(-0.400)/(2\pi\sqrt{0.0625})\approx0.4269\). The positive-covariance model assigns a different density to this same-sign pair of deviations. This comparison concerns the density at one point; it does not by itself establish that the correlated model fits better overall.
Integration in Practice
For two random effects, the marginal likelihood integral is two-dimensional. It is commonly evaluated numerically, for example by Gaussian quadrature or Monte Carlo methods. A convenient perspective is to write \(b_i=AZ\), where \(AA^\mathsf{T}=D\) and \(Z\sim N_q(0,I)\). Then the integral is an expectation of the conditional likelihood over standard-normal \(Z\). Numerical methods approximate that expectation, with the covariance structure built into the transformation \(A\).
Worked Example: A Quadrature Approximation to One Subject’s Contribution
Use the covariance matrix from the previous example, whose Cholesky factor is \(A=\begin{pmatrix}0.50&0\\0.30&0.40\end{pmatrix}\), so \(AA^\mathsf{T}=D\). A two-node-per-coordinate standard-normal quadrature uses \(z_1,z_2\in\{-1,+1\}\), with weight \(0.25\) for each of the four pairs. Suppose, hypothetically, the conditional likelihood evaluations at \(b=Az\) are as follows.
| \(z\) | \(b=Az\) | Conditional likelihood \(f_i\) |
|---|---|---|
| \((+1,+1)^\mathsf{T}\) | \((0.50,0.70)^\mathsf{T}\) | 0.08 |
| \((+1,-1)^\mathsf{T}\) | \((0.50,-0.10)^\mathsf{T}\) | 0.12 |
| \((-1,+1)^\mathsf{T}\) | \((-0.50,0.10)^\mathsf{T}\) | 0.20 |
| \((-1,-1)^\mathsf{T}\) | \((-0.50,-0.70)^\mathsf{T}\) | 0.16 |
The quadrature approximation to the marginal subject contribution is \(0.25(0.08+0.12+0.20+0.16)=0.25(0.56)=0.140\). This is a numerical approximation to an integral, not an exact evaluation. The example shows how a method that averages conditional contributions under standard-normal quadrature nodes incorporates the random-effects covariance through the transformation from \(z\) to \(b\).
Assumptions, Comparisons, and Diagnostics
The normal random-effects specification makes several substantive and computational assumptions: the random-effects mean is zero; \(D\) is common across subjects in the modeled population; random-effects vectors are independent across subjects; and \(D\) is positive definite. The conditional likelihood must also be correctly specified for the observed data given the random effects. A well-behaved numerical integral cannot compensate for a seriously misspecified conditional model or random-effects distribution.
A diagonal \(D\) is the main simple alternative: it assumes independent random intercepts and slopes. An unstructured \(D\), as above, estimates their covariance as well as their variances. More restrictive covariance structures may reduce the number of parameters, while more flexible distributions may be considered when normality is inadequate. These choices involve a trade-off: added flexibility may better represent heterogeneity but can make estimation less stable, especially with limited information per subject.
Useful warning signs include an estimated variance close to zero, a correlation near \(+1\) or \(-1\), unstable estimates under small changes to the model, or numerical integration results that change appreciably when quadrature is refined. Such signs can reflect weak information, an overly complex covariance structure, or a model that does not adequately represent the data. They are prompts for investigation, not automatic proof that the model is wrong. Checking residual patterns and sensitivity to reasonable random-effects assumptions can help distinguish these possibilities.
Time units and origin matter for interpreting the covariance. Changing the unit of time rescales the random slope and its covariance with the intercept. Centering time changes which point on the trajectory the intercept represents and can change the intercept-slope covariance. Therefore, a covariance estimate should always be interpreted alongside the time scale and centering used in the longitudinal submodel.
Common Pitfall
A related error is to treat \(\phi_q(b;D)\) as the probability of observing the exact vector \(b\). For continuous random effects, an individual point has probability zero. The density is a weight per unit volume, and probabilities are obtained by integrating it over regions.
Check Your Understanding
Try each question, then compare your work with the solutions below.
- For a bivariate normal random-effects model, what does the determinant term in the density contribute, and what does the quadratic form measure?
- If \(\sigma_{01}=0\) in a bivariate normal model, what additional assumption does this make about the random intercept and slope?
- For \(D=\begin{pmatrix}0.16&0\\0&0.09\end{pmatrix}\), calculate the density at \(b=(0,0)^\mathsf{T}\).
- Why is a conditional likelihood evaluated at one estimated random-effects vector not generally the marginal subject likelihood?
- Name one change in the analysis setup that can alter the interpretation of the intercept-slope covariance without changing the observed subjects.
Solutions. (1) The determinant gives the density’s scale through the volume of the covariance ellipsoid; \(b^\mathsf{T}D^{-1}b\) measures squared distance from the mean relative to that covariance. (2) Under the joint normal model, zero covariance implies independence of the two random effects. (3) Here \(|D|=0.16(0.09)=0.0144\), \(|D|^{1/2}=0.12\), and the quadratic form at zero is zero, so the density is \(1/(2\pi\cdot0.12)\approx1.326\) per unit area in the random-effects coordinates. (4) The marginal contribution averages conditional likelihood values over the full random-effects distribution, rather than using only one value. (5) Changing the time origin changes the point at which the intercept is defined; changing time units rescales the slope and its covariance with the intercept.