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

