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

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})?$$

• What makes you think that an additional model matrix ? – Robert Long Jan 10 '19 at 10:03
• @RobertLong, I am not sure I get your question. Seems like you got cut off. Would you mind explaining what you mean? Cheers – DomB Jan 10 '19 at 10:19
• Well, you added a new "level" to a basic 2-level model. Why do you think that requires an additional model matrix ? This will be handled within $\mathbf{Z}$, and what you have written as $u_2$ and $u_3$ will be the components of $\mathbf{u}$. If this doesn't make sense I will expand with a full answer with an example a bit later today. I suspect you may be getting confused between the conventions used in "multilevel" analysis and those in mixed model methodology (which are entirely equivalent but often presented differently. – Robert Long Jan 10 '19 at 10:37

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. • Hence, compared to density of a simple fixed-effect model, the joint density of a mixed-model changes because we now incorporate the model matrix$Z$to account for the random effects. However,$Z\$ may contain different levels such as class, school etc. Super answer! – DomB Jan 10 '19 at 13:58

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.