Say we have a mixed-effects model with a single grouping factor, indexed with $i$:
$$ y_i = X_i\beta + Z_ib_i + \epsilon_i \\ \epsilon_i \sim \mathcal N(0, \sigma^2) \\ b_i \sim \mathcal N(0, \Sigma) $$
$X_i$ is each group's model-matrix of predictors for fixed effects, and $Z_i$ is each group's model-matrix of predictors for random-effects.
In this setting, we get a simple closed-form likelihood for each group, because each $y_i$ is multivariate normal:
$$ y_i \sim \mathcal N(X_i\beta,\ Z_i \Sigma Z^T_i + \sigma^2I) $$
In the above the $y_i$s are one-dimensional vectors -- i.e., our model has a univariate response.
Question: What is the likelihood when we instead have a multivariate response?
For example, if we have a bivariate response:
$$ y_i = (Y^T_{i1},Y^T_{i2})^T \\ y_{i1} = X_{i1}\beta_1 + Z_{i1}b_{i1} + \epsilon_{i1} \\ y_{i2} = X_{i2}\beta_2 + Z_{i2}b_{i2} + \epsilon_{i2} \\[1em] \epsilon_i \sim \mathcal N\left(0, \begin{bmatrix} \sigma_1^2 & \rho\sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 & \sigma_2^2 \end{bmatrix}\right) \\[2em] b_i \sim \mathcal N\left(0, \begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{12} & \Sigma_{22} \end{bmatrix}\right) \\ $$
Then I believe $y_i$ should still be multivariate normal, but I don't quite understand how to calculate the covariance:
$$ y_i \sim \mathcal N(\begin{bmatrix}X_{i1}\beta_1\\X_{i2}\beta_2\end{bmatrix},?) $$
