What is the Joint Density Function of a Three-Level Mixed-Effects Model?

Question

This is a follow-up question to a question I posted earlier. Obviously, maximum-likelihood estimation of mixed-effects models requires the joint density function.

Let us assume a two-level mixed effects model:

$y = X\beta+Zu+e$

where $y$ is a vector of n observable random variables, $\beta$ is a vector of $p$ fixed effects, $X$ and $Z$ are known matrices, and $u$ and $e$ re vectors of $q$ and $n$ random effects such that $E(u) = 0$ and $E(e) = 0$ and

$ Var \begin{bmatrix} u \\ e \\ \end{bmatrix} = \begin{bmatrix} G & 0 \\ 0 & R \\ \end{bmatrix}\sigma^2$

where $G$ and $R$ are known positive definite matrices and $\sigma^2$ is a positive constant. The joint density we require to estimate $u$ and $\beta$ is

$f(y,u) = (2\pi\sigma^2)^{-1/2n-1/2q}(det \begin{bmatrix} G & 0 \\ 0 & R \\ \end{bmatrix})^{-1/2}\times\\ exp(-\frac 1{2\sigma^2} \begin{pmatrix} u \\ y-x\beta-Zu \\ \end{pmatrix}' \begin{bmatrix} G & 0 \\ 0 & R \\ \end{bmatrix}^{-1} \begin{pmatrix} u \\ y-x\beta-Zu \\ \end{pmatrix}).$

Now, let us assume a three-level mixed effects model:

$y = X\beta+Zu_{2}+Wu_{3}+e$

where $y$ is a vector of n observable random variables, $\beta$ is a vector of $p$ fixed effects, $X$, $Z$, and $W$ are known matrices, and $u_{2}$, $u_{3}$ as well as $e$ are vectors of $q$, $k$, and $n$ random effects such that $E(u_{2}) = 0$, $E(u_{3}) = 0$ and $E(e) = 0$ and

$ Var \begin{bmatrix} u_{2} \\ u_{3} \\ e \\ \end{bmatrix} = \begin{bmatrix} G & 0 & 0 \\ 0 & K & 0\\ 0 & 0 & R \\ \end{bmatrix}\sigma^2$

where $G$, $K$, and $R$ are known positive definite matrices and $\sigma^2$ is a positive constant. The subscripts of $u$ indicate the level of the random effect.

Would the joint density then be

$f(y,u_{2},u_{3}) = (2\pi\sigma^2)^{-1/2n-1/2q-1/2k}(det \begin{bmatrix} G & 0 & 0 \\ 0 & K & 0\\ 0 & 0 & R \\ \end{bmatrix})^{-1/2}\times\\ exp(-\frac 1{2\sigma^2} \begin{pmatrix} u_{2}\\u_{3} \\ y-x\beta-Zu \\ \end{pmatrix}' \begin{bmatrix} G & 0 & 0 \\ 0 & K & 0\\ 0 & 0 & R \\ \end{bmatrix}^{-1} \begin{pmatrix} u_{2}\\ u_{3} \\ y-x\beta-Zu \\ \end{pmatrix})?$

Answer 1

Maximum likelihood estimation of mixed models typically works with the marginal likelihood of the observed response outcome data $y$. This marginal likelihood is obtained by integrating out the random effects from the joint density, i.e., for the $i$-th sample unit we have $$\left \{ \begin{array}{l} y_i = X_i\beta + Z_ib + \varepsilon_i,\\\\ b_i \sim \mathcal N(0, D), \quad \varepsilon_i \sim \mathcal N(0, \Sigma). \end{array} \right.$$

Based on this model, the log-likelihood function of the linear mixed model is $$\begin{eqnarray} \ell(\theta) & = & \sum_{i = 1}^n p(y_i; \theta)\\ & = & \sum_{i = 1}^n \int p(y_i \mid b_i; \theta) \, p(b_i; \theta) \; db_i, \end{eqnarray}$$ where $\theta$ denotes the parameters of the model, namely the fixed effects $\beta$ and the unique elements of the covariance matrices, $D$ and $\Sigma$. The first term in the second line is the multivariate normal density from the model $[y \mid b]$, and the second term the multivariate normal density for the random effects. Now, in the case of linear mixed models and because the random effects enter linearly in the mean of the model $[y \mid b]$, the two normal distributions 'work' together, and the integral has a closed-form solution. Namely, the marginal model for $[y_i]$ is $$y_i \sim \mathcal N(X_i \beta, \; Z_i D Z_i^\top + \Sigma).$$

A couple of notes:

• Even if you have nested random effects, the formulation remains the same. The $Z$ matrix only gets another form (often being sparse).
• Most software that fit mixed models under maximum likelihood actually work with the implied marginal model and give you estimates of $\theta$. That is, they do not actually fit the mixed model. You could construct situations in which the implied marginal model may come from two different mixed models. In this case, you cannot tell which of the two is the correct one. The random-effect estimates you obtain, for example using function ranef() in R, typically come from a second separate step using empirical Bayes methodology.
• The integrals work nicely together in the case of linear mixed models (and few other cases). If you have categorical outcome data and normal random effects, you need to approximate the integrals numerically.

Answer 2

Using the general mixed model formulation, the log-likelihood of the 3-level model will have the same form as that of the 2-level model. This is because the formulation of the linear mixed effects model as $\mathbf{y} = \mathbf{X\beta}+\mathbf{Z}\boldsymbol{b}+\mathbf{e}$ does not say anything about the number of grouping variables, and in the case of more than 1, whether they they are nested, crossed, or partially crossed. See here for a discussion of nested vs crossed random effects. The model matrix $\mathbf{Z}$ incorporates the structure of all the random effects, while the vector $\mathbf{b}$ contains the conditional modes of all the random effects.

Note that I am using $\mathbf{b}$ rather than $\mathbf{u}$ for the random effects vector above. The reason for this should be evident from the following example in R which may help to illustrate:

Let’s simulate some nested 3-level data. Here we envisage each row of data to be an observation of a students height. Students are nested within classes and classes within schools

set.seed(15)
df1 <- expand.grid(pupil = c(1,2,3), class = c(1,2,3), school = c(1,2))
df1$Height <- rnorm(nrow(df1), 1.5, 0.2) + df1$class + df1$school

Now we fit a linear mixed model, using lme4 with a single grouping variable (a 2-level model) of pupils nested in classes only:

library(lme4)
lmm0 <- lmer(Height ~ 1 + (1|class), data = df1)
getME(lmm0, "Zt")
getME(lmm0, "b")

The last 2 lines extract the transpose of $\mathbf{Z}$ and the conditional modes of the random effects (lme4 uses $\mathbf{b}$ rather than $\mathbf{u}$)

3 x 18 sparse Matrix of class "dgCMatrix"
   [[ suppressing 18 column names ‘1’, ‘2’, ‘3’ ... ]]

1 1 1 1 . . . . . . 1 1 1 . . . . . .
2 . . . 1 1 1 . . . . . . 1 1 1 . . .
3 . . . . . . 1 1 1 . . . . . . 1 1 1

3 x 1 Matrix of class "dgeMatrix"
           [,1]
[1,] -0.9421207
[2,] -0.0119456
[3,]  0.9540663

Now we fit the 3-level model:

lmm1 <- lmer(Height ~ 1 + (1|school/class), data = df1)
getME(lmm1, "Zt")
getME(lmm1, "b")

8 x 18 sparse Matrix of class "dgCMatrix"
   [[ suppressing 18 column names ‘1’, ‘2’, ‘3’ ... ]]

1:1 1 1 1 . . . . . . . . . . . . . . .
1:2 . . . . . . . . . 1 1 1 . . . . . .
2:1 . . . 1 1 1 . . . . . . . . . . . .
2:2 . . . . . . . . . . . . 1 1 1 . . .
3:1 . . . . . . 1 1 1 . . . . . . . . .
3:2 . . . . . . . . . . . . . . . 1 1 1
1   1 1 1 1 1 1 1 1 1 . . . . . . . . .
2   . . . . . . . . . 1 1 1 1 1 1 1 1 1

8 x 1 Matrix of class "dgeMatrix"
           [,1]
[1,] -1.2863150
[2,] -0.6696624
[3,] -0.4058639
[4,]  0.3810632
[5,]  0.6372665
[6,]  1.3435117
[7,] -0.1097978
[8,]  0.1097978

So we can see that the additional “level” is handled by the structure of $\mathbf{Z}$ and all the random effects are in $\mathbf{b}$

So this means that the general form of the log-likelihood is the same for the 3-level model as it is for the 2-level model.

For further information about the structure of $\mathbf{Z}$ I highly recommend this free book by the primary author of lme4 and this paper.

References:
Bates, D., Mächler, M., Bolker, B., & Walker, S. (2014). Fitting linear mixed-effects models using lme4. arXiv preprint arXiv:1406.5823.

Bates, Douglas M. "lme4: Mixed-effects modeling with R." (2010): 470-474.