Let, \mathbf y_i = \mathbf X_i\mathbf\beta + \mathbf Z_i\mathbf b_i+ \mathbf\epsilon_i,
where
\mathbf y_i\sim N(\mathbf X_i\mathbf\beta, \Sigma_i=\sigma^2\mathbf I_{n_i}+\mathbf Z_i \mathbf G\mathbf Z_i'),
\mathbf b_i\sim N(\mathbf 0, \mathbf G),
\mathbf\epsilon_i\sim N(\mathbf 0...