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

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}B_1 + Z_{i1}b_{i1} + \epsilon_{i1}
\\
y_{i2} = X_{i2}B_2 + Z_{i2}b_{i2} + \epsilon_{i2}
\\
\epsilon_i \sim \mathcal N(0, \begin{bmatrix}
    \sigma_1^2 & \rho\sigma_1\sigma_2  \\
    \rho\sigma_1\sigma_2 & \sigma_2^2 
  \end{bmatrix})
\\
b_i \sim \mathcal N(0, \begin{bmatrix}
    \Sigma_{11} & \Sigma_{12}  \\
    \Sigma_{12} & \Sigma_{22}
  \end{bmatrix})
\\
$$

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(X_i\beta,?)
$$