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)
::grouped_glance(
broomExtradata = 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)
::grouped_glance(
broomExtradata = 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")
<- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
fm1
# df.residual in `glance`
::tidy(fm1)
broom.mixed
::glance(fm1)
broom.mixed
# fetch the p-values
::model_parameters(fm1) %>%
parameters::easystats_to_tidy_names()
broomExtra
)
# p-values across groups in mixed-model
%>%
iris group_by(Species) %>%
group_modify(., ~ broomExtra::easystats_to_tidy_names(as.data.frame(parameters::model_parameters(
::lmer(Sepal.Length ~ Petal.Length + Petal.Width + Petal.Length * Petal.Width + (1 | Sepal.Width),
lme4data = .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(
::lmer(Sepal.Length ~ Petal.Length + Petal.Width + Petal.Length * Petal.Width + (1 | Sepal.Width),
lme4data = .x,
control = lme4::lmerControl(optimizer = "bobyqa")
) )))
# linear mixed effects model
::grouped_tidy(
broomExtradata = 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
<- clear.labels(df.nirs) %>%
tmp group_by(metric) %>%
group_modify(., ~ broomExtra::easystats_to_tidy_names(as.data.frame(parameters::model_parameters(
::lmer(value ~ Connection*Task + group*Connection + group*Task + group*Connection*Task + (1|Subject),
lme4data = .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("->", "→")
%>%
) ::clean_names("upper_camel") %>%
janitor::rename(
dplyr"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.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.
<- df.nirs %>%
tbl.SubjectEffect group_by(metric) %>%
do(tidy(lm(value ~ Subject, .))) # requires broom
<- df.nirs %>%
tbl.SubjectEffect group_by(metric) %>%
::anova_test(value~Subject)
rstatix# 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
::p_load(sjPlot, sjmisc, ggplot2, lme4)
pacman
<- lm(mpg ~ wt + cyl + gear + disp, data = mtcars)
model <- report::report(model)
r table_short(r)
::check_model(Model) # will make some plots to check your model
performance
# Individually checking model.
# Note: I will do this if I see something problematic after running `check_model`
::model_parameters(model)
parameters
::check_collinearity(model)
performance::plot(check_collinearity(model)) performance
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
::p_load(sjPlot, sjmisc, ggplot2, lme4)
pacman
<- lm(mpg ~ wt + cyl + gear + disp, data = mtcars)
model <- report::report(model)
r table_short(r)
::check_model(Model) # will make some plots to check your model
performance
# Individually checking model.
# Note: I will do this if I see something problematic after running `check_model`
::model_parameters(model)
parameters
::check_collinearity(model)
performance::plot(check_collinearity(model)) performance
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://m-clark.github.io/mixed-models-with-R/random_slopes.html
Package that gives pseudo-R^2 for mixed models
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