I am trying to understand how to stratify two different mixed effect models. For this, I am reading this guide for mixed effect models.


In this article, they mentioned that in a Conditional three-level growth model of the form.

m1 <- lmer(y ~ time*tx + 
             (time | subjects:therapist) + 
             (1 | therapist) + 
             (0 + tx | therapist) +
             (0 + time | therapist) +
             (0 + time:tx | therapist), data=df)

enter image description here

enter image description here

I want to apply a similar reasoning for my two models :

fit 1 <- lmer(lVolume ~ Country * Day + (1 | Rat.ID))
fit 2 <- lmer(lVolume ~ Country * Day + (1 + Day| Rat.ID))

In both models, I am interested in the prediction of exp(fitted) (I fit an exponential model, that is why I took lVolume).

For model fit1, if I am not wrong, what I am estimating is :

$$ \begin{aligned} \operatorname{lVolume}_{i} &\sim N \left(\alpha_{j[i]} + \beta_{1}(\operatorname{Day}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\gamma_{0}^{\alpha} + \gamma_{1}^{\alpha}(\operatorname{Country}_{\operatorname{England}}) + \gamma_{2}^{\alpha}(\operatorname{Day} \times \operatorname{Country}_{\operatorname{England}}), \sigma^2_{\alpha_{j}} \right) \text{, for Rat.ID j = 1,} \dots \text{,J} \end{aligned} $$

For model fit2, if I am not wrong, what I am estimating is : $$ \begin{aligned} \operatorname{lVolume}_{i} &\sim N \left(\alpha_{j[i]} + \beta_{1j[i]}(\operatorname{Day}), \sigma^2 \right) \\ \left( \begin{array}{c} \begin{aligned} &\alpha_{j} \\ &\beta_{1j} \end{aligned} \end{array} \right) &\sim N \left( \left( \begin{array}{c} \begin{aligned} &\gamma_{0}^{\alpha} + \gamma_{1}^{\alpha}(\operatorname{Country}_{\operatorname{England}}) \\ &\gamma^{\beta_{1}}_{0} + \gamma^{\beta_{1}}_{1}(\operatorname{Country}_{\operatorname{England}}) \end{aligned} \end{array} \right) , \left( \begin{array}{cc} \sigma^2_{\alpha_{j}} & \rho_{\alpha_{j}\beta_{1j}} \\ \rho_{\beta_{1j}\alpha_{j}} & \sigma^2_{\beta_{1j}} \end{array} \right) \right) \text{, for Rat.ID j = 1,} \dots \text{,J} \end{aligned} $$

So, I doubt how I can interpret my two models in terms of the organization chart from the guide and get the level equations similar to the guide lecture. I really don't understand the differences in (1 + Day | Rat.ID) and (1 | Rat.ID) because, to me, both models estimate the amount of variation in the rate of change of volume between rats. Is this correct?

I really appreciate any kind of thoughts and help/hints for this. Thanks.


1 Answer 1


Your situation is much simpler than the multi-level model that you first cite as an example. I'd recommend that you start by looking at this site's lmer cheat sheet, and its links, about mixed models. Only then should you look at more complex models.

As the notation in your question is a bit confusing, I'll limit myself to trying to clarify some conceptual issues.

Both your fit1 and fit2 include fixed effects of Country and Day and their interaction. The models only differ in the random effects. Your interpretation of those models seems to be incorrect in a few ways.

First, the coefficients not only for Day but also for Country and for their interaction should appear in the fixed-effect part of the model. You only show Day there explicitly in both attempts at interpreting the lmer models. Perhaps that's implied in the way that you wrote the random-effects parts of the model, but it's hard to say.* The random effects are usually presented as having means of 0, which is not how you wrote them for either of your models.

Second, there is no separate modeling of $\sigma_{a_j}^2$ for each Rat, as one might interpret the way you wrote the random effects for both models. In fit1 a single Gaussian distribution is fit to all of the individual Rat intercepts. (Perhaps that's what you intended, but the extra $j$ subscript is confusing.**) Those random intercepts just allow for variation among Rats about the fixed-effect intercept. To be clear, there is no modeling of Country or the Country:Day interaction within individual Rats (although it's hard to know if that's what you intended to convey). That model does NOT "estimate the amount of variation in the rate of change of volume between rats," just their hypothetical log-volumes at Day = 0 and your baseline value of Country.

In fit2 the (1 + Day| Rat.ID) term does also allow for differences among rats in the Day coefficient, which you can consider the slope with respect to time. Again, the extra $j$ subscripts in your random-effect variance terms is a bit confusing;** there is just one Gaussian distribution of slopes encompassing all Rats. The fit2 also models correlations between the individual intercepts and the slopes with respect to time in a way that is similar to what you write; see the cheat sheet linked above.

In response to comment and related question:

The equations you show are evidently from the extract_eq() function of the equatiomatic package. That can include fixed-effect terms in the mean values of expressions for random-effect terms, as explained for your fit1 in this answer. Furthermore, its choice of symbols used to describe a model can be more complicated than is needed for simple models like yours. That led to some of the confusion in my original interpretation (above) of those equations.

With respect to your fit2, recognize that all terms in mean values for expressions of random effects that don't involve your random effects of rats, indexed by $j$, are part of the fixed-effect portion of the model. I find it simplest to first assume that the coefficients involved in random effects are at their mean values, to obtain the full fixed-effect part of the model. Then evaluate the random-effect (co)variance terms as having zero means, adding to the fixed effects the random-effect terms for each rat $j$.

In your fit2, at the mean values of $\alpha_j$ and $\beta_{1_j}$ for the first term in the extract_eq() output, you would have for the distribution of fixed effects plus error variance $\sigma^2$:

$$lVolume \sim N\left(\gamma_{0}^{\alpha} + \gamma_{1}^{\alpha}(\operatorname{Country}_{\operatorname{England}}) + \left(\gamma^{\beta_{1}}_{0} + \gamma^{\beta_{1}}_{1}(\operatorname{Country}_{\operatorname{England}})\right)(\operatorname{Day}), \sigma^2 \right). $$

Multiplying out the terms in the mean value for this normal distribution, the overall fixed-effect intercept is thus $\gamma_{0}^{\alpha}$; $\gamma_{1}^{\alpha}$ is the fixed-effect coefficient for Country; $\gamma^{\beta_{1}}_{0}$ is the fixed-effect coefficient for Day; $\gamma^{\beta_{1}}_{1}$ is the fixed-effect coefficient for the Day by Country interaction.

Then you can simply think of the random-effect coefficients as distributed bivariate normal with means of 0 and the covariance structure indicated in your answer. Those add extra terms to what's predicted based just on fixed effects. The random intercept term is just added directly. The random slope is first multiplied by the corresponding Day value, then added.

*I now understand that this display is based on the extract_eq() function of the equatiomatic package. You can extract the fixed-effects portion of the model by substituting the mean values of the random effects into the corresponding portions of the model. Then the random effect coefficients have 0 means, to incorporate into terms that include the random effects.

**This is also from extract_eq(). It indexes the rat random effects with $j$ and calls the intercept $\alpha$. So $\sigma^2_{\alpha_j}$ is meant to represent the variance of the random intercept $\alpha_j$; $\sigma^2_{\beta_{1_j}}$ represents the random-effect variance of the $\beta_1$ coefficient.

  • $\begingroup$ Now I understood a little bit more what you ment, it took me a while sorry, which it is still a bit confusing is the mathematical interpretation, let's just say fit1. We have only the Passage*Day interaction. But how can I represent the model fit1 : fit 1 <- lmer(lVolume ~ Country * Day + (1 | Rat.ID)) in terms of expected values and variances if I can ask. $\endgroup$
    – Rosa Maria
    Mar 18, 2022 at 21:32
  • 1
    $\begingroup$ @RosaMaria I think that my answer to another of your questions explains that for fit1. I hadn't appreciated that your equations come from the extract_eq() function in the equatiomatic package. The terminology used by that function can be hard to understand, even in simple cases, as it can include things that are fundamentally fixed effects as part of its formulas for random effects. You seem to have several related questions on this site; please see if you might be able to delete some of them now. $\endgroup$
    – EdM
    Mar 19, 2022 at 18:03
  • $\begingroup$ Hi sorry to bother you again, a technical question: the distribution of the predicted values by the lmer model follows a normal distribution approximately or almost a normal distribution or can it be non normal ? Thanks $\endgroup$
    – Rosa Maria
    Aug 25, 2022 at 12:54
  • 1
    $\begingroup$ @RosaMaria fitting an lmer model assumes that each random effect comes from a normal distribution and that the residual variance is also normally distributed. The actual estimates of the individual random effects need not strictly follow a normal distribution, as the model must fit both the fixed and the random-effect terms together so there will be tradeoffs from ideal normality. The distribution of "predicted values" for the entire model depend on the distributions assumed for the fixed-effect predictors and the choice of values for the random effects. $\endgroup$
    – EdM
    Aug 25, 2022 at 13:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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