# Interpretation of generalized linear mixed-effects models: How should I proceed with a significant interaction term?

My experiment has the following factors:

• "year" (fixed) with levels "1" and "2"
• "age" (fixed) with levels "old" and "new"
• "treatment" (fixed) with levels "control" and "treated"
• "block" (random) with levels ("1_1", "2_1", "3_1", "1_2", "2_2", "3_2")
• "variety" (random) with levels ("A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L")

"Year", "age", and "treatment" are crossed, "id" is nested within "age", and "block" is nested within "year":

$$year \times block(year) \times age \times treatment \times id(age)$$

A visualization of my data:

I am analyzing this experiment with generalized linear mixed effects models using the lmerTest library in R. In the following, I show this for an example of my data:

dat <- structure(list(value = c(579, 575, 740, 546, 692, 490, 535, 436,
962, 474, 867, 824, 549, 954, 511, 876, 556, 462, 443, 620, 573,
626, 664, 975, 545, 846, 868, 973, 744, 597, 685, 768, 535, 460,
585, 589, 1085, 755, 1004, 651, 971, 612, 682, 711, 652, 502,
603, 609, 677, 604, 604, 488, 441, 742, 753, 532, 551, 432, 614,
556, 615, 639, 718, 550, 642, 845, 426, 791, 499, 603, 606, 570,
518, 638, 594, 598, 401, 479, 503, 738, 929, 480, 891, 606, 705,
679, 941, 758, 623, 446, 619, 1077, 873, 826, 1259, 892, 890,
470, 793, 617, 611, 501, 729, 549, 851, 460, 1032, 802, 526,
928, 582, 790, 680, 759, 558, 367, 617, 589, 394, 529, 825, 536,
748, 757, 655, 593, 767, 639, 675, 872, 665, 787, 660, 924, 1031,
590, 813, 783, 736, 472, 629, 504, 898, 420, 607, 781, 687, 761,
1083, 770, 508, 598, 923, 542, 548, 591, 526, 537, 526, 570,
632, 919, 797, 637, 855, 811, 613), year = structure(c(2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L), levels = c("1", "2"), class = "factor"),
age = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L), levels = c("old", "new"), class = "factor"),
treatment = c("control", "treated", "control", "treated",
"control", "treated", "treated", "control", "control", "treated",
"treated", "control", "treated", "control", "treated", "control",
"control", "treated", "treated", "control", "treated", "control",
"treated", "control", "treated", "control", "control", "treated",
"control", "treated", "control", "treated", "control", "treated",
"treated", "control", "control", "treated", "treated", "control",
"control", "treated", "treated", "control", "control", "treated",
"treated", "control", "control", "treated", "control", "treated",
"treated", "control", "control", "treated", "control", "treated",
"treated", "control", "treated", "control", "control", "treated",
"treated", "control", "treated", "control", "control", "treated",
"treated", "control", "treated", "control", "treated", "control",
"treated", "control", "treated", "control", "control", "treated",
"treated", "treated", "control", "treated", "control", "control",
"treated", "treated", "control", "control", "treated", "treated",
"control", "control", "treated", "treated", "control", "treated",
"control", "treated", "control", "treated", "control", "treated",
"control", "control", "treated", "control", "treated", "control",
"treated", "control", "treated", "treated", "control", "control",
"treated", "treated", "control", "treated", "control", "control",
"treated", "treated", "control", "treated", "control", "treated",
"control", "treated", "control", "control", "treated", "treated",
"control", "treated", "control", "treated", "control", "treated",
"control", "treated", "control", "control", "treated", "treated",
"control", "control", "treated", "treated", "control", "treated",
"control", "control", "treated", "treated", "control", "treated",
"control", "control", "treated", "control", "treated", "treated",
"control"), block = c("1_2", "1_2", "2_2", "2_2", "3_2",
"3_2", "4_2", "4_2", "1_1_ru", "1_1_ru", "2_1_ru", "2_1_ru",
"3_1_ru", "3_1_ru", "1_2", "1_2", "2_2", "2_2", "3_2", "3_2",
"4_2", "4_2", "1_1_ru", "1_1_ru", "3_1_ru", "2_1_ru", "3_1_ru",
"3_1_ru", "1_2", "1_2", "2_2", "2_2", "3_2", "3_2", "4_2",
"4_2", "1_1_ru", "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru",
"3_1_ru", "1_2", "1_2", "2_2", "2_2", "3_2", "3_2", "4_2",
"4_2", "1_1_ru", "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru",
"3_1_ru", "1_2", "1_2", "2_2", "2_2", "3_2", "3_2", "4_2",
"4_2", "1_1_ru", "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru",
"3_1_ru", "1_2", "1_2", "2_2", "2_2", "3_2", "3_2", "4_2",
"4_2", "1_1_ru", "1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru",
"1_2", "1_2", "2_2", "2_2", "3_2", "3_2", "4_2", "4_2", "1_1_ru",
"1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", "3_1_ru", "1_2",
"1_2", "2_2", "2_2", "3_2", "3_2", "4_2", "4_2", "1_1_ru",
"1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", "3_1_ru", "1_2",
"1_2", "2_2", "2_2", "3_2", "3_2", "4_2", "4_2", "1_1_ru",
"1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", "3_1_ru", "1_2",
"1_2", "2_2", "2_2", "3_2", "3_2", "4_2", "4_2", "1_1_ru",
"1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", "3_1_ru", "1_2",
"1_2", "2_2", "2_2", "3_2", "3_2", "4_2", "4_2", "1_1_ru",
"1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", "3_1_ru", "1_2",
"1_2", "2_2", "2_2", "3_2", "3_2", "4_2", "4_2", "1_1_ru",
"1_1_ru", "2_1_ru", "2_1_ru", "3_1_ru", "3_1_ru"), id = c("A",
"A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A",
"A", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B",
"B", "B", "B", "C", "C", "C", "C", "C", "C", "C", "C", "C",
"C", "C", "C", "C", "C", "D", "D", "D", "D", "D", "D", "D",
"D", "D", "D", "D", "D", "D", "D", "E", "E", "E", "E", "E",
"E", "E", "E", "E", "E", "E", "E", "E", "E", "F", "F", "F",
"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "G", "G",
"G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G",
"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
"H", "H", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I",
"I", "I", "I", "I", "J", "J", "J", "J", "J", "J", "J", "J",
"J", "J", "J", "J", "J", "J", "K", "K", "K", "K", "K", "K",
"K", "K", "K", "K", "K", "K", "K", "K", "L", "L", "L", "L",
"L", "L", "L", "L", "L", "L", "L", "L", "L", "L")), row.names = c(NA,
167L), class = "data.frame")


I fit the full design as:

model_full <- lmerTest::lmer(value ~ treatment*age*year*(1|block)*(1|id),
data = dat,
REML = TRUE)

summary(model_full)


Is this the correct way to write the design formula? Can I indicate that the random factor "id" is nested in the fixed factor "age"?

My model returns:

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: value ~ treatment * age * year * (1 | block) * (1 | id)
Data: dat

REML criterion at convergence: 2012.2

Scaled residuals:
Min       1Q   Median       3Q      Max
-2.42972 -0.63655 -0.08531  0.65485  3.06605

Random effects:
Groups   Name        Variance Std.Dev.
id       (Intercept)  1881.3   43.37
block    (Intercept)   126.7   11.26
Residual             14691.6  121.21
Number of obs: 167, groups:  id, 12; block, 7

Fixed effects:
Estimate Std. Error      df t value Pr(>|t|)
(Intercept)                     821.72      34.23   36.46  24.003  < 2e-16 ***
treatmenttreated               -240.38      40.41  144.50  -5.948 1.95e-08 ***
agenew                           56.86      48.06   48.99   1.183 0.242457
year2                          -136.51      38.76   49.24  -3.522 0.000934 ***
treatmenttreated:agenew         116.85      57.58  144.14   2.029 0.044248 *
treatmenttreated:year2           88.92      53.45  144.27   1.664 0.098378 .
agenew:year2                   -101.57      53.92  144.20  -1.884 0.061597 .
treatmenttreated:agenew:year2   -18.60      75.92  144.06  -0.245 0.806794
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) trtmnt agenew year2  trtmn: trtm:2 agnw:2
tretmnttrtd -0.590
agenew      -0.687  0.420
year2       -0.647  0.521  0.438
trtmnttrtd:  0.414 -0.702 -0.608 -0.366
trtmnttrt:2  0.446 -0.756 -0.318 -0.689  0.530
agenew:yer2  0.442 -0.374 -0.649 -0.684  0.542  0.495
trtmnttr::2 -0.314  0.532  0.461  0.485 -0.758 -0.704 -0.710


I interpret the significant (P = 0.044) interaction between "treatment" and "age" as that the change of "value" between "control" and "treated" conditions depends on whether a sample belongs to the "old" or "new" level. How would I compare this change in response to treatment between the "old" and "new" levels? I am interested in the direction of the change and whether it is smaller/larger. Could I generate subsets of my data for the "old" and "new" conditions and use generalized linear mixed-effects models on these to test the effect of "treatment"?

dat_old <- subset(dat, age == "old")
dat_new <- subset(dat, age == "new")

model_old <- lmerTest::lmer(value ~ treatment*year*(1|block)*(1|id),
data = dat_old,
REML = TRUE)

model_new <- lmerTest::lmer(value ~ treatment*year*(1|block)*(1|id),
data = dat_new,
REML = TRUE)


Also, is there a good way to add the results of these generalized linear mixed-effects models to a plot like the one above?

I have some comments before looking at the treatment effect comparison between "old" and "new".

This is a linear mixed model (a linear regression with random effects), not a generalized linear mixed model:

• linear model: $$\operatorname{E}(Y | X) = X\beta$$
• generalized linear model: $$g(\operatorname{E}(Y | X)) = X\beta$$ where $$g$$ is the (non-linear) link function

Furthermore, {lme4} ignores the product terms between the fixed effects and the random effects, and between the two random effects. This is clear in the summary output; see the "Random effects" and the "Fixed effects" parts. The actual LMM formula is:

value ~ treatment * age * year + (1 | block) + (1 | id)

• The mean structure, $$\operatorname{E}(Y | X)$$, is treatment * age * year.
• There are two crossed random effects, block and id.

Nesting describes a structuring of random effects (actually, of the variability attributed to these random effects). A random factor (id) cannot be nested within a fixed factor (age). You can let old and new id's have different variances by replacing (1 | id) with (0 + age || id). This doesn't improve the model. Similarly, block is not nested in year.

Most importantly, you don't need to split the data by age and fit two separate models to answer questions about the treatment effect for new vs old individuals. I use the {marginaleffects} package to illustrate how to do this analysis. You'll need to decide how to deal with the the year covariate as it interacts with treatment and age. Do you want to look at the treatment effect for years 1 and 2 separately, average over the two years, or compare the treatment effect between years?

I interpret the significant (P = 0.044) interaction between "treatment" and "age" as that the change of "value" between "control" and "treated" conditions depends on whether a sample belongs to the "old" or "new" level.

Not quite. The interaction term treatmenttreated:agenew corresponds to this comparison:

\begin{aligned} \left(\operatorname{E}Y_{\text{treated, new, year=1}} - \operatorname{E}Y_{\text{control, new, year=1}}\right) - \left(\operatorname{E}Y_{\text{treated, old, year=1}} - \operatorname{E}Y_{\text{control, old, year=1}}\right) \end{aligned}

The comparison is complex due to the three-way interaction. First, it concerns the expected response in year 1 only. The interaction itself is a difference in differences: it tells you whether there is a difference in the expected improvement due to treatment when we compare new and old individuals.

Let's do this calculation explicitly with the help of {marginaleffects}. See also the tutorial on Comparisons in the Marginal Effects Zoo (the package documentation).

# Calculate the treatment effect, treated - control,
# for each age in year 1.
comparisons(
fit,
variables = "treatment",
newdata = datagrid(year = "1", age = c("old", "new")),
re.form = NA
)
#>       Term          Contrast year age Estimate Std. Error     z Pr(>|z|)
#>  treatment treated - control    1 old     -240       40.4 -5.95   <0.001
#>  treatment treated - control    1 new     -124       41.0 -3.01   0.0026

# Compare the treatment effect between old (b1) and new (b2) in year 1.
comparisons(
fit,
hypothesis = "b2 - b1 = 0",
variables = "treatment",
newdata = datagrid(year = "1", age = c("new", "old")),
re.form = NA
)
#>     Term Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
#>  b2-b1=0     -117       57.6 -2.03   0.0424 4.6  -230     -4


Same as the treatmenttreated:agenew row in the summary output.

Let's also look at how the change in response to treatment between old and new changes itself from year 1 to year 2.

# Compare the treatment and the control (ie. estimate the difference treated - control
# for each combination of year and age.
comparisons(
fit,
variables = "treatment",
newdata = datagrid(year = c("1", "2"), age = c("old", "new")),
re.form = NA
)
#>       Term          Contrast year age Estimate Std. Error     z Pr(>|z|)
#>  treatment treated - control    1 old   -240.4       40.4 -5.95   <0.001
#>  treatment treated - control    1 new   -123.5       41.0 -3.01   0.0026
#>  treatment treated - control    2 old   -151.5       35.0 -4.33   <0.001
#>  treatment treated - control    2 new    -53.2       35.0 -1.52   0.1283

# Suppose the four rows/comparisons in the table above are labeled b1 to b4, and
# formulate a more complex contrast as a difference in differences.
# Of course, each b term itself is a difference between treated and controls.
comparisons(
fit,
hypothesis = "(b4 - b3) - (b2 - b1) = 0",
variables = "treatment",
newdata = datagrid(year = c("1", "2"), age = c("old", "new")),
re.form = NA
)
#>               Term Estimate Std. Error      z Pr(>|z|)   S 2.5 % 97.5 %
#>  (b4-b3)-(b2-b1)=0    -18.6       75.9 -0.245    0.806 0.3  -167    130


Hopefully, you are not surprised that this computation gives the same result as the treatmenttreated:agenew:year2 term in the summary table. If you are surprised and want to understand this better, you'll have to read up on design matrices and the different ways to parametrize a (generalized) linear model. We're using a treatment contrast matrix, which is the default in R.

However, there is no need to do this reading! Instead, learn to express the comparisons you are interested in as a combination of treatment, age and year. For example, you might be interested in the treatment effect for "old" individuals averaged over years 1 and 2. (This comparison doesn't correspond exactly to one of the coefficients in the summary table.)

Now, with this great power comes great responsibility. Unless the comparisons you make are specified in advance (ie. the analysis is pre-registered), if you look at multiple contrasts in a kind of exploratory analysis to discover "significant" ones, you will have to apply a correction for multiple hypothesis testing.