A simple guide
In simple words, the rationale behind mixed models is the simultaneous presence of components that model the EXPECTED VALUE (fixed) as well as the VARIANCE (random).
Specifically, today we are going to cover LINEAR MIXED MODELS, which make the following assumptions (Bolker et al. 2022):
The expected values of the responses are linear combinations of the fixed predictor variables and the random effects.
The conditional distribution of the variable response is Gaussian (equivalently, the errors are “Normal”).
The random effects are normally distributed. Normally, we assume that the expected value of a random effect is equal to zero, but with a positive variance.
The most used packages and/or functions for frequentist LMMs:
For more information about mixed models in R, please, check out the new CRAN Task View: Mixed, Multilevel, and Hierarchical Models in R.
Well, from the hand of Tidyverse, the “tidy data” framework changed the way we code and work in R for data science. Tidy datasets are easy to manipulate, model and visualize, and have a specific structure (Wickham 2014):
Each variable is a column,
Each observation is a row, and
Each value have its own cell.
Tidy-data structure. Following three rules makes a dataset tidy: variables are in columns, observations are in rows, and values are in cells. Source: (Wickham and Grolemund 2017).
After planning the experimental design, identifying dependent variable, independent variable(s), conducting the experiment, and collecting the data…the expected path would be as follows:
flowchart LR
A[Dig the data] --> B{Model Selection}
B --> C[Significant\nEffects?]
subgraph Inference
C -->|Yes| D[Comparisons]
endflowchart LR
A[Dig the data] --> B{Model Selection}
B --> C[Significant\nEffects?]
subgraph Inference
C -->|Yes| D[Comparisons]
end
This workflow has been created using mermaid
For a complete experience, I recommend you to download and install the:
We are going to explore the new features offered by Quarto documents (*.qmd).
Quarto is a refined version (and successor) of R-markdown. It is an open-source scientific and technical publishing system built on Pandoc. It allows to combine R, Python, Julia, and Java (Observable JS) for the design of reports, applications, websites, blogs, scientific publications, and more…
#| echo: true
#| collapse: true
# install.packages("easypackages") # To load and/or install the packs
library(easypackages)
packages("readxl") # To open/save excel files
packages('dplyr', "tidyr","purrr", "forcats", "stringr") # Data wrangling
packages("nlme", "lme4") # Mixed models libraries
packages("broom", "broom.mixed") # Managing models results
packages("performance") # Check assumptions and performance
packages("emmeans","multcomp","multcompView",
"car", "multilevelmod") # Aov and mult comp
packages("ggplot2") # Figures
packages("agricolae") # Miscellaneous functionsAmong many packages we are going to use today, there is one from the tidymodels family that was specially designed to convert statistical objects into tidy format: the broom package (Robinson, Hayes, and Couch 2022). Particularly for mixed models, we will also need its spinoff: the broom.mixed package (Bolker and Robinson 2022) .
The broom package offer three key functions to manipulate models’ outcomes:
glance() to report information about the entire model
tidy() to summarize information about model components, and
augment() to add information about observations to a dataset
Now, we are going to reproduce the analysis I’ve done for one of my papers (Correndo et al. 2021). Particularly, we are going to reproduce Figure 2
For this paper, we have data from 4 different locations. We tested the levels of soil potassium fertility, hereinafter as soil test K (STK), in long-term experiments (2000-2009) where the treatments of interest were: (i) Control (unfertilized), (ii) NPS (fertilized with NPS), and (iii) Pristine conditions (No Ag-history).
At each plot/sample, the STK was measured at five-consecutive soil depths (0-20, 20-40, 40-60, 60-80, and 80-100 cm). Thus, they we took “repeated measurements” over the space.
We were NOT interested in comparing locations since they had very different previous history, and crop rotation, so confounding effects may have obscured the inference. Therefore, site was not a factor under consideration, and all the analysis were fitted independently by site.
# Read file
rm_data <-
read.csv("../data/repeated_measures_data.csv", header = TRUE) %>%
# Create PLOT column to identify subject (Exp. Unit for Rep. Measures)
unite(PLOT, BLOCK,TREAT, sep = "_", remove=FALSE) %>%
# OR
# Identify Subplot
ungroup() %>%
group_by(BLOCK, TREAT) %>%
# Create plot ID # Needed for Repeated Measures
mutate(plot = cur_group_id(), .after = PLOT) %>%
ungroup() %>%
mutate(DEPTH = as.factor(DEPTH),
depth = as.integer(DEPTH), # Needed for CorAR1
BLOCK = factor(BLOCK),
SITE = factor(SITE),
TREAT = factor(TREAT),
# Create a grouping variable (WHY?) # Needed for HetVar
GROUP = case_when(TREAT == "Pristine" ~ "Pristine",
TRUE ~ "Agriculture")
)
# File online? Try this...(remove "#")
# url_rm <- "https://raw.githubusercontent.com/adriancorrendo/univ6020/master/coding/data/repeated_measures_data.csv"
#rm_data <- read_excel(url_rm, col_names = TRUE)Now, let’s use several functions to explore the data.
First, the glimpse() function from dplyr
# Glimpse from dplyr
dplyr::glimpse(rm_data)Rows: 180
Columns: 9
$ SITE <fct> site_1, site_1, site_1, site_1, site_1, site_1, site_1, site_1, …
$ PLOT <chr> "I_Control", "I_Control", "I_Control", "I_Control", "I_Control",…
$ plot <int> 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 7, 7, 7, 7, 7, 2, 2, 2, 2, 2, 5, 5…
$ TREAT <fct> Control, Control, Control, Control, Control, Control, Control, C…
$ BLOCK <fct> I, I, I, I, I, II, II, II, II, II, III, III, III, III, III, I, I…
$ DEPTH <fct> 10, 30, 50, 70, 90, 10, 30, 50, 70, 90, 10, 30, 50, 70, 90, 10, …
$ STK <int> 105, 83, 103, 110, 127, 119, 98, 106, 107, 109, 132, 100, 96, 10…
$ depth <int> 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2…
$ GROUP <chr> "Agriculture", "Agriculture", "Agriculture", "Agriculture", "Agr…Then, the skim() function from skmir
# Skim from skimr
skimr::skim(rm_data)| Name | rm_data |
| Number of rows | 180 |
| Number of columns | 9 |
| _______________________ | |
| Column type frequency: | |
| character | 2 |
| factor | 4 |
| numeric | 3 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| PLOT | 0 | 1 | 5 | 12 | 0 | 9 | 0 |
| GROUP | 0 | 1 | 8 | 11 | 0 | 2 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| SITE | 0 | 1 | FALSE | 4 | sit: 45, sit: 45, sit: 45, sit: 45 |
| TREAT | 0 | 1 | FALSE | 3 | Con: 60, NPS: 60, Pri: 60 |
| BLOCK | 0 | 1 | FALSE | 3 | I: 60, II: 60, III: 60 |
| DEPTH | 0 | 1 | FALSE | 5 | 10: 36, 30: 36, 50: 36, 70: 36 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| plot | 0 | 1 | 5.00 | 2.59 | 1 | 3.00 | 5 | 7.00 | 9 | ▇▇▃▇▇ |
| STK | 0 | 1 | 181.84 | 63.09 | 74 | 138.75 | 174 | 221.25 | 406 | ▅▇▅▂▁ |
| depth | 0 | 1 | 3.00 | 1.42 | 1 | 2.00 | 3 | 4.00 | 5 | ▇▇▇▇▇ |
And let’s use ggplot2 for a better look
# Boxplot
rm_data %>%
dplyr::select(-depth) %>%
# Plot
ggplot() +
# Boxplots
geom_boxplot(aes(x = reorder(DEPTH, desc(DEPTH)), y = STK, fill = TREAT))+
# Axis labels
labs(x = "Soil depth (cm)", y = "STK (g/m2)")+
# Plot by site
facet_wrap(~SITE)+
# Flip axes
coord_flip()+
# Set scale type
scale_x_discrete()+
# Change theme
tidybayes::theme_tidybayes()
I’m sorry for this, but the most important step is ALWAYS to write down the model.
In a traditional approach blocks are defined as fixed, affecting the mean of the expected value. Yet there is no consensus about treating blocks as fixed or as random. For more information, read Dixon (2016).
In my opinion, we should treat blocks as random because the number of blocks we choose is a random sample of all the potential blocks we could have used. However, for a decent estimation of a random effect, the usual amount of blocks we use in agriculture research (3-5) is not sufficient.
Let’s define the model. For simplification (and avoid writing interaction terms), here we are going to consider that \(\tau_i\) is the “treatment”.
\[ y_{ij} = \mu + \tau_i + \beta_j + \epsilon_{ij} \]
\[ \epsilon_{ij} \sim N(0, \sigma^2_{e} )\] where \(\mu\) represents the overall mean (if intercept is used), \(\tau_i\) is the effect of treatment-j over \(\mu\), \(\beta_j\) is the effect of block-j over \(\mu\), and \(\epsilon_{ij}\) is the random effect of each experimental unit.
# SIMPLEST MODEL
fit_block_fixed <- function(x){
lm(# Response variable
STK ~
# Fixed
TREAT + DEPTH + TREAT:DEPTH + BLOCK,
# Data
data = x)
}An alternative approach is considering a MIXED MODEL, where blocks are considered “random”. Basically, we add a term to the model that it is expected to show a “null” overall effect over the mean of the variable of interest but introduces “noise”. By convention, a random effect is expected to have an expected value equal to zero but a positive variance as follows: \[ y_{ij} = \mu + \tau_i + \beta_j + \epsilon_{ij} \] \[ \beta_j \sim N(0, \sigma^2_{b} )\] \[ \epsilon_{ij} \sim N(0, \sigma^2_{e} )\] Similar than before, \(\mu\) represents the overall mean (if intercept is used), \(\tau_i\) is the effect of treatment-j over \(\mu\), \(\beta_j\) is the “random” effect of block-j over \(\mu\), and \(\epsilon_{ij}\) is the random effect of each experimental unit.
So what’s the difference? Simply specifying this component: \[ \beta_j \sim N(0, \sigma^2_b) \], which serves to model the variance.
How do we write that?
# RANDOM BLOCK
fit_block_random <- function(x){
nlme::lme(
# Fixed
STK ~ TREAT + DEPTH + TREAT:DEPTH,
# Random
random = ~1|BLOCK,
# Data
data = x)
}Run the candidate models
STK_models <-
rm_data %>%
# Let's group data to run multiple locations|datasets at once
group_by(SITE) %>%
# Store the data per location using nested arrangement
nest() %>%
# BLOCK as FIXED
mutate(model_0 = map(data, fit_block_fixed)) %>%
# BLOCK as RANDOM
mutate(model_1 = map(data, fit_block_random)) %>%
# COMPOUND SYMMETRY
mutate(model_2 = map(data, fit_corsymm)) %>%
# AUTO-REGRESSIVE ORDER 1
mutate(model_3 = map(data, fit_ar1)) %>%
# COMPOUND SYMMETRY + HETEROSKEDASTIC
mutate(model_4 = map(data, fit_corsymm_hetvar) ) %>%
# Data wrangling
pivot_longer(cols = c(model_0:model_4), # show alternative 'contains' model
names_to = "model_id",
values_to = "model") %>%
# Map over model column
mutate(results = map(model, broom.mixed::augment )) %>%
# Performance
mutate(performance = map(model, broom.mixed::glance )) %>%
# Extract AIC
mutate(AIC = map(performance, ~.x$AIC)) %>%
# Extract coefficients
mutate(coef = map(model, ~coef(.x))) %>%
# Visual-check plots
mutate(checks = map(model, ~performance::check_model(.))) %>%
ungroup()Checking assumptions is always important. To learn more about data exploration, tools to detect outliers, heterogeneity of variance, collinearity, dependence of observations, problems with interactions, among others, I highly recommend reading (Zuur, Ieno, and Elphick 2010).
(site_1_models %>% dplyr::filter(model_id == "model_0"))$checks[[1]]
(site_2_models %>% dplyr::filter(model_id == "model_0"))$checks[[1]]
(site_3_models %>% dplyr::filter(model_id == "model_0"))$checks[[1]]
(site_4_models %>% dplyr::filter(model_id == "model_0"))$checks[[1]]
Compare models performance
# Visual model selection
best_STK_models <-
STK_models %>%
group_by(SITE) %>%
# Use case_when to identify the best model
mutate(best_model =
case_when(AIC == min(as.numeric(AIC)) ~ "Yes",
TRUE ~ "No")) %>%
ungroup()
# Plot
best_STK_models %>%
ggplot()+
geom_point(aes(x = model_id, y = as.numeric(AIC),
color = best_model, shape = best_model),
size = 3)+
facet_wrap(~SITE)
# Final models
selected_models <- best_STK_models %>% dplyr::filter(best_model == "Yes")Estimate the effects of factors under study (and their interaction)
models_effects <-
selected_models %>%
# Type 3 Sum of Squares (Partial SS, when interactions are present)
mutate(ANOVA = map(model, ~Anova(., type = 3)) )
# Extract ANOVAS
models_effects$ANOVA[[1]]Analysis of Deviance Table (Type III tests)
Response: STK
Chisq Df Pr(>Chisq)
(Intercept) 326.261 1 < 2.2e-16 ***
TREAT 39.994 2 2.067e-09 ***
DEPTH 12.637 4 0.013191 *
TREAT:DEPTH 21.827 8 0.005247 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1models_effects$ANOVA[[2]]Analysis of Deviance Table (Type III tests)
Response: STK
Chisq Df Pr(>Chisq)
(Intercept) 794.721 1 < 2.2e-16 ***
TREAT 21.394 2 2.261e-05 ***
DEPTH 67.224 4 8.743e-14 ***
TREAT:DEPTH 52.957 8 1.099e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# MULTCOMPARISON
# emmeans and cld multcomp
# We need to specify ourselves the most important interaction to perform the comparisons
mult_comp <-
models_effects %>%
# Comparisons estimates (emmeans)
mutate(mc_estimates = map(model, ~emmeans(., ~ TREAT*DEPTH))) %>%
# Assign letters and p-value adjustment (multcomp)
mutate(mc_letters =
map(mc_estimates,
~as.data.frame(
# By specifies a strata or level to assign the letters
cld(., by = "DEPTH", decreasing = TRUE, details=FALSE,
reversed=TRUE, alpha=0.05, adjust = "tukey", Letters=LETTERS))))Now, we are going to reproduce Figure 2
# Create data frame for plot
plot_df <- mult_comp %>%
dplyr::select(SITE, mc_letters) %>%
unnest(mc_letters)
# Define your own colors
my_colors <- c("#ffaa00", "#7E5AA0", "#5c9c8c")
# Create the plot
STK_plot <-
plot_df %>%
# We need to re-express DEPTH from factor to character, and then to numeric
mutate(DEPTH = as.numeric(as.character(DEPTH))) %>%
# Re-order levels of the factors
mutate(TREAT = fct_relevel(TREAT,"Control", "NPS", "Pristine")) %>%
mutate(SITE = fct_relevel(SITE,"site_1", "site_2", "site_3", "site_4")) %>%
# Create plot
ggplot()+
# 01. LAYOUT
## Subplots
facet_wrap(~SITE, nrow = 2)+
## Axis titles
labs(x = "Soil depth (cm)", y = bquote(~NH[4]*'OAc-K (g' ~m^-2*')'))+
# 02. GEOMETRIES
## i. Points
geom_point(aes(x = DEPTH, y = emmean,
fill= TREAT,
shape = TREAT),
size = 3, col = "black")+
## Adjust shape aesthetics
scale_shape_manual(name="Fertilizer Scenario", values=c(24,23,21),
guide="legend")+
scale_colour_manual(name="Fertilizer Scenario",
values = my_colors,
guide='legend')+
scale_fill_manual(name="Fertilizer Scenario",
values = my_colors,
guide='legend')+
## ii. Add error bar
geom_errorbar(width = 0.25, aes(x = DEPTH, color = TREAT,
ymin = emmean-2*SE, ymax = emmean+2*SE))+
## iii. Add line
geom_line(size = 0.5,aes(x = DEPTH, y = emmean, color = TREAT))+
# 03. ADJUST XY AXIS
## Reverse the scale
scale_x_reverse(breaks=seq(0, 100, 20), limits = c(100,0))+
coord_flip()+
# 04. THEME
theme_bw()+
theme(strip.text = element_text(size = rel(1.25)),
strip.background = element_blank(),
# Grid
panel.grid = element_blank(),
# Axis
axis.title = element_text(size = rel(1.5)),
axis.text = element_text(size = rel(1.25), color = "black"),
# Legend
legend.position = "top", legend.title = element_blank(),
legend.text = element_text(size = rel(1.25)) )STK_plot
Figure 2. Soil profiles of STK (\(g~m^{-2}\)) under three different conditions: pristine soils (green circles), under grain cropping from 2000 to 2009 with no fertilizers added (Control, orange triangles), and under grain cropping from 2000 to 2009 with N, P, plus S fertilization (NPS, purple diamonds). Overlapping error bars indicate absence of significant differences between scenarios by soil depths combinations (Tukey’s HSD, p < 0.05).
☺
Adding references with Quarto has become quite easy. Let’s see…
This is a citation of the Tidyverse (Wickham et al. 2019) package.
Using a “*.bib” file: Figures were produced using the ggplot2 package (Wickham 2016). Models check were tested with the performance package (Lüdecke et al. 2021).
Using Footnotes: For example, this is some text with a Footnote1, this is a second Footnote2
Using visual editor: this option introduced by Quarto is simply awesome! 🤩. Let’s see. Insert -> Citation or “Crtl + Shift + F8”. With this option we can look for citations online via DOI, Crossref, etc… and insert them into our document (and to our *.bib file).
