Chapter 13 Mixed Models

13.1 To-Do

Look at the presentation (Google Slides) I made for Keith’s lab (the one Ashley asked me to make). It has a ton of information I can use here.

13.2 Preamble

I made mixed models its own section because it is not your typical analyses. I would highly recommend that any student starts with more traditional analyses such as ANOVA or linear regression before taking on mixed model analyses.

In order to do pipe friendly analyses broomExtra is a nice package to run at first. You can check out the GitHub repo to learn more about the functions (which also work with other statistical models beyond LMMs).

For a theoretical background I would recommend this chapter from the book “Just enough R”.

13.3 lmerControl

Within lme4 you can specify optimizers. The default is bobyqa the other common option is nloptwrap.

To learn more on optimizer checks see here

13.4 Code

# to speed up computation, let's use only 50% of the data
set.seed(123)

# linear model (model summaries across grouping combinations)
broomExtra::grouped_glance(
  data = dplyr::sample_frac(tbl = ggplot2::diamonds, size = 0.5),
  grouping.vars = c(cut, color),
  formula = price ~ carat - 1,
  ..f = stats::lm,
  na.action = na.omit
)

# linear mixed effects model (model summaries across grouping combinations)
broomExtra::grouped_glance(
  data = dplyr::sample_frac(tbl = ggplot2::diamonds, size = 0.5),
  grouping.vars = cut,
  ..f = lme4::lmer,4
  formula = price ~ carat + (carat | color) - 1,
  control = lme4::lmerControl(optimizer = "bobyqa")
  
  
  fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

# df.residual in `glance`
broom.mixed::tidy(fm1)

broom.mixed::glance(fm1)

# fetch the p-values
parameters::model_parameters(fm1) %>%
  broomExtra::easystats_to_tidy_names()
)


# p-values across groups in mixed-model
iris %>%
  group_by(Species) %>%
  group_modify(., ~ broomExtra::easystats_to_tidy_names(as.data.frame(parameters::model_parameters(
    lme4::lmer(Sepal.Length ~ Petal.Length + Petal.Width + Petal.Length * Petal.Width + (1 | Sepal.Width),
               data = .x,
               control = lme4::lmerControl(optimizer = "bobyqa")
    )
  ))))

easystats_to_tidy_names has been retired from everything I see online. Below is the only code I could find to replace it. It won’t extract p-values. Might need to check if its now part of ipmisc instead of broomExtra. I don’t think its a required argument anymore. You can see the original issue on GitHub here.

iris %>%
  group_by(Species) %>%
  group_modify(., ~ as.data.frame(parameters::model_parameters(
    lme4::lmer(Sepal.Length ~ Petal.Length + Petal.Width + Petal.Length * Petal.Width + (1 | Sepal.Width),
               data = .x,
               control = lme4::lmerControl(optimizer = "bobyqa")
    )
  )))

  # linear mixed effects model
  broomExtra::grouped_tidy(
    data = dplyr::mutate(MASS::Aids2, interval = death - diag),
    grouping.vars = sex,
    ..f = lme4::lmer,
    formula = interval ~ age + (1 | status),
    control = lme4::lmerControl(optimizer = "bobyqa"),
    tidy.args = list(conf.int = TRUE, conf.level = 0.99)
  )

Here is code I used from a recent project

tmp <- clear.labels(df.nirs) %>%
  group_by(metric) %>%
  group_modify(., ~ broomExtra::easystats_to_tidy_names(as.data.frame(parameters::model_parameters(
    lme4::lmer(value ~ Connection*Task + group*Connection + group*Task + group*Connection*Task + (1|Subject),
               data = .x,
               control = lme4::lmerControl(optimizer = "bobyqa")
    )
  )))) %>%
  filter(term != "(Intercept)") %>% # remove intercept from report
  mutate(term =  str_replace_all(term, ":", "×") %>% 
           str_replace_all("Connection", "") %>%
           str_replace_all("Task", "") %>%
           str_replace_all("group.L", "Group") %>%
           str_replace_all("->", "→")
  ) %>%
  janitor::clean_names("upper_camel") %>%
  dplyr::rename(
    "p" = "PValue",
    "dferror" = "DfError"
  ) %>%
  filter((str_detect(Term, "Group")) & p < .1) %>% # Give me only the results that show group differences, set p @.1
  mutate_if(is.numeric, round, 5) # round to 5 digits (makes it easier to read)

View(tmp)

To visualize your mixed model you can use sjPlot as shown below or here

pacman::p_load(sjPlot, sjmisc, ggplot2)
data(efc)
theme_set(theme_sjplot())

# make categorical
efc$c161sex <- to_factor(efc$c161sex)

# fit model with interaction
fit <- lm(neg_c_7 ~ c12hour + barthtot * c161sex, data = efc)

plot_model(fit, type = "pred", terms = c("barthtot", "c161sex"))



plot_model(fit, type = "int")

plot_model(fit, type = "pred", terms = c("c161sex", "barthtot [0, 100]"))

fit <- lm(neg_c_7 ~ c12hour + c161sex * barthtot, data = efc)
plot_model(fit, type = "int")

plot_model(fit, type = "int", mdrt.values = "meansd")


# fit model with 3-way-interaction
fit <- lm(neg_c_7 ~ c12hour * barthtot * c161sex, data = efc)

# select only levels 30, 50 and 70 from continuous variable Barthel-Index
plot_model(fit, type = "pred", terms = c("c12hour", "barthtot [30,50,70]", "c161sex"))

plot_model(fit, type = "int")

13.5 Post-hoc in mixed models

This can be a bit tricky because you want to include your random effect.

You could use

emmeans::emmeans()
lmerTest::difflsmeans()
multcomp::glht()

I would look at this StackOverflow post which describes some differences between each approach

This is a good breakdown.

13.5.1 Examples

# using glht package
summary(glht(YOUR MODEL, linfct=mcp(YOUR FIXED FACTOR="Tukey")))

13.6 Running your models

For starters I am a big fan of using something that allows me to run through several models to get an overall feel for the data. This usually involves either broom or rstatix package which allows for pipe friendly modelling. I will commonly do this in order to evaluate if random factors should be included in my final model. If you are running mixed models you can use broomExtra or broom.mixed. I would recommend broomExtra which depends on broom.mixed (so its hits 2 birds with one stone).

# Question: Is there a subject effect in our data?
# if there is then we introduce it as a random effect.

tbl.SubjectEffect <- df.nirs %>% 
  group_by(metric) %>% 
  do(tidy(lm(value ~ Subject, .))) # requires broom

tbl.SubjectEffect <- df.nirs %>% 
  group_by(metric) %>% 
  rstatix::anova_test(value~Subject)
# answer is yes. There is a subject effect that we must accomodate

For a very complete guide on running ANOVAs in R for your first time. This post is immaculate.

You can see some exercises here

13.7 Modeling your data

To test the normality of your data you can use a few different methods

Plot your data
Check skewness and kurtosis
Shapiro test.

In order to visualize your linear model I use a few packages

sjPlot and sjmisc
jtools
interactions
ggeffects (used less)
To visualize lmer mixed models see here
To follow up a mixed-model here

# Load typical packages
pacman::p_load(sjPlot, sjmisc, ggplot2, lme4)

model <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars)
r <- report::report(model)
table_short(r)

performance::check_model(Model) # will make some plots to check your model 

# Individually checking model. 
# Note: I will do this if I see something problematic after running `check_model`

parameters::model_parameters(model)

performance::check_collinearity(model)
performance::plot(check_collinearity(model))

I am still working on code that will do this when I run my mixed models in a loop. But the GitHub page should help

# Load typical packages
pacman::p_load(sjPlot, sjmisc, ggplot2, lme4)

model <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars)
r <- report::report(model)
table_short(r)

performance::check_model(Model) # will make some plots to check your model 

# Individually checking model. 
# Note: I will do this if I see something problematic after running `check_model`

parameters::model_parameters(model)

performance::check_collinearity(model)
performance::plot(check_collinearity(model))

13.7.1 Mixed Models: Misc Info

Contrast Analysis post-hoc mixed models
Visualize your mixed models
Good introduction: https://ourcodingclub.github.io/tutorials/mixed-models/
More intro: https://www.jaredknowles.com/journal/2013/11/25/getting-started-with-mixed-effect-models-in-r
https://dynamicecology.wordpress.com/2015/02/05/how-many-terms-in-your-model-before-statistical-machismo/
https://m-clark.github.io/mixed-models-with-R/random_slopes.html
Package that gives pseudo-R^2 for mixed models
- https://cran.r-project.org/web/packages/jtools/vignettes/summ.html#report_robust_standard_errors
- https://rinterested.github.io/statistics/mixed_effects_comparison.html
easystats blog. Has some great posts with code included.

13.8 Looking at interactions

src: https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html

13.8.1 Packages

mertool is a gui for modelling –> https://www.jaredknowles.com/journal/2015/8/12/announcing-mertools tidymodels: https://www.tidymodels.org/learn/statistics/tidy-analysis/

To assess model performance I use the performance package. You can visualize your model using this vignette