Simulate a 2-level Multilevel/HLM/linear mixed model

I think multilevel models are still “a thing”, right? Like they used to be all hot and sexy a few years ago but maybe they’ve started to become a little more common. The HLM software has made them very popular, although eventually I will rant about how much I dislike the whole “Level 1” and “Level 2” parlance that Raudenbush & Bryk helped disseminate. But for better or worse it’s here to stay and I need it if I want to communicate with my colleagues.

In any case, the model from which we will be simulating data in this post looks like this:

$y = X\mathbf{\beta} + Z\mathbf{u} + \mathbf{\epsilon}$

Where:
$y$ is an nX1 continuous vector of data (the dependent variable).
$X$ is an nXp model matrix with a pX1 vector of coefficients $\mathbf{\beta}$ (so this is the ‘fixed effects’ part of the model).
$Z$ is another nXq model matrix with a qX1 vector $\mathbf{u}$ (the ‘random effects’).
$\mathbf{\epsilon}$ are the residuals.

I know people reading this are probably more familiar withe the Level 1 and Level 2 approach of looking at things and I’m pretty sure it’s easy to make explicit its relationship with the way I stated the model, which is the Laird-Ware (1982) specification (this is, I think, a lot more common than the Raudenbush & Bryk one). Since I will be relying on matrices and matrix operations to simulate data, I thought it was a better idea to keep everything Laird-Ware.

In any case, by assuming that the random effects $\mathbf{u}$ are distributed $\mathbf{u} \sim N(0, \sigma^{2}D)$ (which is what we typically do when we fit a multilevel model) then $Var(y) = Var(Z\mathbf{u}) + Var(\mathbf{\epsilon})$. From regular OLS regression the assumption is that $Var(\mathbf{\epsilon})= \sigma^{2}I$ to which we simply need to add $Var(Z\mathbf{u})$. After doing some basic matrix algebra we now have $\epsilon = \sigma^{2}I+\sigma^{2}ZDZ^{t})$ so that $y \sim N(X\beta, \sigma^{2}(I+ZDZ^{t}))$

Because we’ve worked through the math a little bit here, now we know what we need to do in order to simulate data the way we want it: specify $X\beta$ for the fixed effects, the matrix $D$ for the covariance of the random effects and $\sigma^{2}$ for the residual.

::UPDATE::

Following the advice of quite a few people, I would like to direct you to the improved version of this function that is now in MY PUBLISHED ARTICLE. This way you get a nice function that has been vetted by peer-review AND you can cite something if you choose to use it.

The code is all the way at the end of the article. Enjoy!