I know that when analyzing trial level data (i.e., treating every observation of a given individual as a separate data point), something must be done to account for this non-independence. I've always just taken for granted that including random subject terms in the model does this. But now I'm wondering exactly how they do this, and to what extent? Does having a random subject slope in addition to the intercept do a "better" job at modeling that non-independence? Or is it an all or nothing sort of thing? Is a random subject intercept always sufficient?



It seems like there are two main facets to this question. First, where in the math do these dependencies come from, and second, how do we make that useful for statistics. I'll begin with the math.

Let's say we've got a one-way random effects model. Then we have $$ y_{ij} = \mu + \alpha_i + \e_{ij} $$ where $\alpha \sim \mathcal N(0, \sigma^2_\alpha I) \perp \e \sim \mathcal N(0, \sigma^2_\e I)$ and $i= 1,\dots,m$ and $j=1,\dots,n_i$. Let $n=\sum_i n_i$. This means $$ \Cov(y_{ij}, y_{kl}) = \begin{cases}\sigma^2_\e + \sigma^2_\alpha & i=k, j=l \\ \sigma^2_\alpha & i=k, j\neq l \\ 0 & i\neq k\end{cases}. $$ so by adding a random intercept we are modeling a useful, but very narrow, type of correlation. We are saying that each pair of distinct observations from a given subject have the same covariance, and that this covariance is the same no matter what subject we are considering.

It's not hard to imagine situations where this model is wrong. Many real situations have the variability proportional to the mean (e.g. think about the variation in measuring the weight of a mouse vs a blue whale) and that violates this model. Similarly maybe the measures are more tightly clustered for some subjects than others so it doesn't make sense to use the same correlation structure for every subject. And etc. This is a very useful correlation structure relative to how easy it is to specify and implement, but it also is very easily wrong (quote George Box here).

Now let's add random slopes. Let's say we have $$ y_{ij} = \mu + \alpha_i + \eta x_{ij} + \tau_ix_{ij} + \e_{ij} $$ where $\alpha \sim \mathcal N(0, \sigma^2_\alpha I)$, $\tau \sim \mathcal N(0, \sigma^2_\tau I)$, and $\e \sim \mathcal N(0 ,\sigma^2 I)$. We'll say $\e$ is independent of $\alpha$ and $\tau$ but we'll allow $\tau$ and $\alpha$ to be correlated within a group to give us translation invariance in $x$ (see e.g. page 7 of the lme4 manual by Bates et al.). Let $\rho=\Cov(\alpha_i, \tau_i)$.

We can work out the covariance structure to be $$ \Cov(y_{ij}, y_{kl}) = \begin{cases} \sigma^2_\alpha + \sigma^2_\tau x_{ij}^2 + 2 x_{ij} \rho + \sigma^2_\e & i=k, j=l \\ \sigma^2_\alpha + x_{ij}x_{kl}\sigma^2_\tau + \rho(x_{ij} + x_{kl}) & i=k, j\neq l \\ 0 & \text{o.w.} \end{cases} $$

That's a much more complicated beast, but we can immediately see some things. First, we haven't changed the fact that we are modeling different subjects as independent. All we've done is added a fancier covariance structure for the data coming from within each subject.

Secondly, the variance parameters $\sigma^2_\alpha$, $\sigma^2_\tau$, $\rho$, and $\sigma^2_\e$ are still the same between subjects so our previous comments about the possible violations of that still apply. Also note how this directly generalizes the covariance structure of the 1-way random effects model: if $\sigma^2_\tau = 0$ (forcing $\tau_1 = \dots = \tau_m = 0$) and then necessarily $\rho = 0$ too, this will reduce to a random intercepts model (although we still have our fixed effects slope in there).

The main change is the appearance of $x_{ij}$ terms. The variance goes up with $|x_{ij}|$ which I think makes sense because surely for a slope term the variance ought to be related to the magnitudes of the values in question since it depends on the scale of the data. We also have that the covariance of distinct within-subject measurements depends on the values of $x_{ij}$ and $x_{il}$. If both $x_{ij}$ and $x_{il}$ have the same sign and are large in magnitude then this will be large, otherwise it will be smaller. or go negative (and for more moderate values the sign of $\rho$ matters too). That's pretty much taking into account the correlation between $x_{ij}$ and $x_{il}$ (in spirit).

So there you have it. The way that mixed models achieve this modeling of dependence is by adding random variables with carefully constructed covariance structures. The variance parameters are also estimated so it is fit to the data that way, and they can be estimated to be essentially (or exactly, depending on how they're estimated) zero, in which case you've added unnecessary complexity, squandered degrees of freedom, and lost power.

But I don't think this is really the helpful level at which to make modeling decisions. Look at the covariance in random slopes: if I've got a model with a random intercept and a fixed slope, and I'm considering adding a random slope, I wouldn't look at that covariance and decide if it really seems like a better model for my data with $\sigma^2_\tau > 0$. No, I would think scientifically about what it means to allow the slopes to change by subject. So if I think it makes sense that the rate of change for each person is not a constant $\beta$ but instead a random $\beta + \tau_i$, then that correlation structure is the way to do it. But the science comes first, and the exact details of the correlation are in my opinion almost an afterthought to make the scientific intent of the model perform as advertised.

To emphasize this point, consider the canonical sleep study data (again see the lme4 manual):

sleep (credit: page 4 of the lme4 manual)

I think it makes a lot of sense to use a random slope here. Statistically it's significant, but more importantly scientifically I would bet that some people handle sleep deprivation better than others so we'd expect to see a range in rates per subject. And I make this decision without any reference to model parameters.

  • $\begingroup$ +6. Very nice answer, detailed and pedagogical. I was disappointed to see it had zero upvotes so I decided to give it a bounty :) $\endgroup$ – amoeba says Reinstate Monica May 5 '18 at 20:12
  • $\begingroup$ @amoeba thank you! I've been trying to get better at mixed models recently so it's nice to get feedback that I'm not just making this stuff up :) $\endgroup$ – jld May 5 '18 at 20:34
  • $\begingroup$ @Chaconne very nice answer! To add a little perspective, the so-called MMRM (bit.ly/2wBonxN) used to be and still is to some extent the state-of-the-art model for modelling longitudinal clinical trial data. This model estimates the full subject-level covariance matrix and thus takes all possible correlations between repeated measures into account. The random intercept (compound-symmetry) and random coefficient models that you describe are simplifications of this "super model". $\endgroup$ – Rune H Christensen May 9 '18 at 20:45
  • $\begingroup$ @RuneHChristensen thank you, and thanks for the link to that paper. I hadn't heard of that before, it's very interesting $\endgroup$ – jld May 10 '18 at 16:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.