rlinear-regressionlme4mixed-modelsglmm

How to report random effect in the mixed effects model

I use the iris dataset in R as the example below.

library(lme4)
mixed.fit.b <- lmer(Sepal.Width ~ Sepal.Length + (1+ Sepal.Length|Species), data = iris)
summary(mixed.fit.b)

ranef(mixed.fit.b)$Species
coef(mixed.fit.b)$Species
predict(mixed.fit.b)

The random intercept and slope are shown below.

           (Intercept) Sepal.Length
setosa     -0.49549054   0.78331501
versicolor  1.19197858   0.26689317
virginica   1.17260303   0.27282273

That means the intercept is -0.49549054 (fixed + random intercept) and slope is 0.78331501 (fixed + random slope) for setosa right? So, there are three couples of intercepts and slopes. In a general linear model, we can say the y = intercept + slope and the y changed a slope per x. But in mixed models, there are three three couples of intercepts and slopes. How to report them?

Solution

  • Recommendation

    I agree with Roland that using a mixed model with three clusters is not advisable (see Gelman & Hill, 2007, linked below). Speaking to your main question, we would like to know how much random variation is being teased apart by the model, and in fact this is really supposed to be used as justification to use the model in the first place (see Meteyard & Davies, 2020, where they propose fitting random effects-only models as a precursor to your full model).

    With this in mind, what I would recommend is reporting these parts of your model to account for the random effects (including what you mentioned). While they are not all necessary, they are certainly informative:

    Worked Example: Fitting the Model

    I will continue to use the data you have as an example, but I will change it to something that converges, as your first model was singular, which is not usable when reporting. I have loaded four libraries necessary for the functions used and fit the model.

    #### Load Libraries ####
library(lmerTest)
library(lattice)
library(broom.mixed)
library(performance)

#### Fit Model ####
fit <- lmer(Petal.Length
                    ~ Petal.Width
                    + (1 + Petal.Width |Species), 
                    data = iris)

    Summarizing the Random Effect Variance

    Next, we can summarize the model, but I prefer to use the tidy function from the broom.mixed package for quickly obtaining RE estimates.

    ##### Summarise Model ####
tidy(fit)

    Shown below:

    # A tibble: 6 × 8
  effect   group    term                      estim…¹ std.e…² stati…³    df p.value
  <chr>    <chr>    <chr>                       <dbl>   <dbl>   <dbl> <dbl>   <dbl>
1 fixed    NA       (Intercept)                 2.43    0.898    2.71  1.97   0.115
2 fixed    NA       Petal.Width                 1.10    0.414    2.65  2.00   0.118
3 ran_pars Species  sd__(Intercept)             1.52   NA       NA    NA     NA    
4 ran_pars Species  cor__(Intercept).Petal.W…  -0.408  NA       NA    NA     NA    
5 ran_pars Species  sd__Petal.Width             0.646  NA       NA    NA     NA    
6 ran_pars Residual sd__Observation             0.362  NA       NA    NA     NA    
# … with abbreviated variable names ¹​estimate, ²​std.error, ³​statistic

    Based off this information, we now know:

    The appendix for one of the cited articles above shows how you can report this information, which is the main question you had if I am not mistaken.

    Caterpillar Plot of Random Effects

    We can report this information in an article if necessary, but we can also plot these effects with a caterpillar plot:

    #### Plot Dotplot of Random Effects ####
dotplot(ranef(fit))

    enter image description here

    This is partly why Roland warned you about having a higher number of random effect clusters. You only see three intercepts and slopes here, which doesn't capture a lot of variation that you would want to get from a mixed model. However, we can make two key insights here:

    1. The virginica plant has much higher petal lengths on average (random intercepts) than the other species.
    2. The versicolor plant has a higher slope term for petal widths than the others.

    Other Random Effects Estimates

    We can also obtain the ICC and pseudo R2 estimates with this code:

    #### Check Performance ####
performance(fit)

    Shown below:

    # Indices of model performance

AIC     |    AICc |     BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma
-----------------------------------------------------------------------------
154.478 | 155.066 | 172.542 |      0.957 |      0.231 | 0.944 | 0.355 | 0.362

    What to note from this summary:

    Citations

    Below are the citations I mentioned earlier. Gelman & Hill is a canonical source for learning about mixed models. The article by Meteyard & Davies is a best-practice guide for running mixed models. Let me know if you found this answer helpful.