# Significance and confidence intervals from emmeans::contrasts on linear mixed model

I'm trying to understand the results from emmeans::contrast applied to a linear mixed model with continuous covariate (WR) and categorical fixed effect (Condition).

A simple version of my post-hoc analysis question is: is there a significant effect of Condition based on the model predictions at a particular set of work rates?

I understand contrast significance is not based on confidence intervals, but I can't quite wrap my head around why some of these contrasts are significant when the model CIs overlap by so much. I also (kinda) understand the limitations of relying on p-values, and I would say the contrasts are not meaningful in my context, regardless of significance. I also understand that I can compare emmeans::emtrends for the linear & quadratic coefficients at each WR value, but I feel like I also want to show a (non-)difference in estimated means at those WR values?

To help understand, I also ran a simple linear model and the results are more what I would expect. I have compared them below.

So I think I'm not understanding how the lmer model SE & CIs are generated for the estimated marginal means and contrasts? I might be misunderstanding something way more basic here. Would appreciate any insight. Sorry, this is a bit long as I'm trying to show my thought process.

library(lme4)
library(lmerTest)
library(emmeans)
library(tidyverse)

## data frame below
model_lm <-
lm(
yvar ~ poly(WR, 2, raw = TRUE) * Condition,
data = df
)

model_lmer <-
lmerTest::lmer(
yvar ~ poly(WR, 2, raw = TRUE) * Condition + (1 + WR | ID),
data = df
)

> summary(model_lm)
...
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)                         -0.631866   0.376649  -1.678   0.0948 .
poly(WR, 2, raw = TRUE)1             2.319495   0.559582   4.145 4.84e-05 ***
poly(WR, 2, raw = TRUE)2             0.147620   0.177503   0.832   0.4065
Condition1                          -0.061505   0.535203  -0.115   0.9086
poly(WR, 2, raw = TRUE)1:Condition1  0.445687   0.808018   0.552   0.5818
poly(WR, 2, raw = TRUE)2:Condition1 -0.001247   0.261320  -0.005   0.9962
---

> summary(model_lmer)
...
Fixed effects:
Estimate Std. Error        df t value Pr(>|t|)
(Intercept)                          -0.13939    0.22676  26.53710  -0.615    0.544
poly(WR, 2, raw = TRUE)1              0.34976    0.53158  15.83207   0.658    0.520
poly(WR, 2, raw = TRUE)2              1.20331    0.09402 211.98383  12.798   <2e-16 ***
Condition1                           -0.07443    0.25336 203.31092  -0.294    0.769
poly(WR, 2, raw = TRUE)1:Condition1   0.44427    0.38294 203.41994   1.160    0.247
poly(WR, 2, raw = TRUE)2:Condition1   0.01715    0.12401 203.53251   0.138    0.890
---

df <-
structure(list(
ID = structure(c(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, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L
), levels = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"
), class = "factor"),
Condition = structure(c(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, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 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, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 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, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L
), levels = c("0", "1"), class = "factor"),
WR = c(0, 0.25, 0.5,
0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 2.89, 0, 0.25,
0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.73, 0, 0.25, 0.5,
0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.35, 0, 0.25, 0.5, 0.75,
1, 1.25, 1.5, 1.75, 2, 2.25, 2.4, 0, 0.25, 0.5, 0.75, 1, 1.25,
1.5, 1.75, 1.98, 0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 1.9,
0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.65, 0,
0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.55, 0, 0.25,
0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5,
0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3,
3.24, 0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.11, 0, 0.25,
0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.05, 0, 0.25, 0.5, 0.75, 1,
1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.72, 0, 0.25, 0.5, 0.75, 1, 1.25,
1.5, 1.75, 2, 2.25, 2.5, 2.64, 0, 0.25, 0.5, 0.75, 1, 1.25, 1.5,
1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 0, 0.25, 0.5, 0.75, 1,
1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 0, 0.25, 0.5,
0.75, 1, 1.25, 1.5, 1.75, 1.99, 0, 0.25, 0.5, 0.75, 1, 1.25,
1.5, 1.75, 1.97, 0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 1.94,
0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2.03),
yvar = c(0, 0,
0, 0, 0, 0.5, 0.5, 2, 4, 5, 7, 8, 9, 0, 0, 0, 0.5, 0.5, 1, 2,
4, 5, 8, 9, 10, 0, 0, 0.5, 0.5, 2, 4, 4, 5, 7, 8, 9, 0, 0, 0,
0.5, 2, 3, 5, 7, 8, 8, 9, 0, 0.5, 0.5, 0.5, 2, 4, 7, 7, 9, 0.5,
0.5, 1, 3, 3, 5, 7, 8, 9, 0, 0.5, 0.5, 1, 2, 2, 3, 4, 4, 6, 8,
10, 0, 0, 0, 0.5, 1, 2, 3, 4, 6, 8, 9, 0, 0, 0, 0, 0.5, 0.5,
0.5, 1, 1, 3, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0.5, 1, 1, 2, 3, 4,
5, 7, 8, 8, 0, 0, 0.5, 0.5, 2, 3, 3, 4, 6, 7, 0, 0, 1, 1, 2,
3, 4, 4, 6, 6, 0, 0.5, 1, 1, 2, 2, 3, 4, 5, 6, 6, 8, 0, 0, 0.5,
0.5, 2, 3, 4, 4, 5, 6, 7, 7, 0, 0, 0, 0, 0, 0.5, 0.5, 2, 3, 4,
5, 6, 7, 9, 10, 0, 0, 0, 0, 0, 0, 0.5, 1, 2, 3, 4, 7, 9, 10,
10, 0, 0, 0.5, 0.5, 2, 3, 4, 7, 10, 0, 0, 0.5, 0.5, 1, 3, 3,
7, 10, 0, 0, 0.5, 0.5, 1, 1, 3, 5, 7, 0, 0, 0.5, 2, 4, 6, 8,
10, 10)), row.names = c(NA, -227L), class = c("tbl_df", "tbl",
"data.frame"))


No significant terms for Condition1 or the interactions. (This model is probably overkill for the current MRE data, but it is more appropriate for my full dataset).

emm_lm <-
emmeans::emmeans(
model_lm,
~ Condition | WR,
at = list(WR = c(0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5)),
)

emm_lmer <-
emmeans::emmeans(
model_lmer,
~ Condition | WR,
at = list(WR = c(0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5)),
)

contrasts_lm <-
emmeans::contrast(
emm_lm,
method = "consec",
simple = list("Condition"),
combine = TRUE,
)

contrasts_lmer <-
emmeans::contrast(
emm_lmer,
method = "consec",
simple = list("Condition"),
combine = TRUE,
)

> contrasts_lm
WR  contrast                estimate    SE  df t.ratio p.value
0   Condition1 - Condition0  -0.0615 0.535 221  -0.115  0.9994
0.5 Condition1 - Condition0   0.1610 0.309 221   0.521  0.9523
1   Condition1 - Condition0   0.3829 0.285 221   1.342  0.5290
1.5 Condition1 - Condition0   0.6042 0.312 221   1.939  0.2098
2   Condition1 - Condition0   0.8249 0.312 221   2.641  0.0432
2.5 Condition1 - Condition0   1.0449 0.375 221   2.790  0.0289
3   Condition1 - Condition0   1.2643 0.621 221   2.037  0.1743
3.5 Condition1 - Condition0   1.4831 1.046 221   1.418  0.4821

P value adjustment: mvt method for 8 tests

> contrasts_lmer
WR  contrast                estimate    SE  df t.ratio p.value
0   Condition1 - Condition0  -0.0744 0.253 203  -0.294  0.9908
0.5 Condition1 - Condition0   0.1520 0.146 203   1.040  0.7192
1   Condition1 - Condition0   0.3870 0.135 203   2.864  0.0235
1.5 Condition1 - Condition0   0.6306 0.148 203   4.274  0.0001
2   Condition1 - Condition0   0.8827 0.148 203   5.968  <.0001
2.5 Condition1 - Condition0   1.1434 0.178 203   6.429  <.0001
3   Condition1 - Condition0   1.4127 0.295 203   4.784  <.0001
3.5 Condition1 - Condition0   1.6906 0.498 203   3.397  0.0044

Degrees-of-freedom method: kenward-roger
P value adjustment: mvt method for 8 tests


The significant contrasts from model_lm make sense to me (and might be spurious anyway), but those from model_lmer I don't fully understand, given the model Condition terms are non-significant and the emmeans SEs & CIs are large at the particular WR values.

Where do the contrast SEs come from? Why are they smaller for the mixed model than for the linear model, when the emmean SEs are larger for the lmem and smaller for the lm? (not shown above, plotted below)

To visualise why I'm confused, I plotted emmeans::emmip and added the contrast significance values. CIs are from the emm_lm & emm_lmer. My intuition would be that the lmer contrasts should be less significant given the larger confidence intervals around the model predictions, but that is obviously wrong?

p.value_lmer <-
tibble(summary(contrasts_lmer)) |>
mutate(
WR = as.numeric(WR),
model = "lmer",
p.format = case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "*",
TRUE ~ "")
) |>
left_join(
tibble(summary(emm_lmer)) |>
group_by(WR) |>
summarise(
y.position = max(emmean) + max(SE) * 1.96
),
by = "WR"
)

plot_lmer <-
emmeans::emmip(
model_lmer,
Condition ~ WR,
CIs = TRUE,
at = list(WR = c(0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5))
) +
geom_point(position = position_dodge(width = 0.1)) +
geom_text(
data = p.value_lmer,
inherit.aes = FALSE,
aes(x = WR, y = y.position, label = p.format),
size = 5) +
coord_cartesian(ylim = c(-1, 20)) +
labs(
title = "Model_LMER",
caption = "* indicates significant contrast Condition1 - Condition0")

# and similar for plot_lm, both figures attached below


And emmeans::pwpp as another visual

If any of that made any sense, I would appreciate any advice!

Each contrast is a linear combination of the regression coefficients; notationally, $$e = co_0 b_0 + c_1b_1 + c2_b2 + \cdots = \mathbf{c}'\mathbf{b}$$ in vector notation, where $$\mathbf{b}$$ is the regression coefficients and $$\mathbf{c}$$ are the weights for the linear combination. The variance-covariance matrix for $$\mathbf{b}$$, denoted $$\mathbf{V}$$, is obtained from the model using the vcov() function, and by standard linear-model theory, $$SE(e) = \sqrt{\mathbf{c}'\mathbf{Vc}}$$.
Consider two estimated means $$m_1$$ and $$m_2$$, then $$Var(m_1 - m_2) = Var(m_1) + Var(m_2) - 2Cov(m_1,m_2)$$. In a fixed-effects model, these estimates are often uncorrelated, making the third term zero. In a mixed model, they can be correlated, and if they are positively correlated then the hird term subtracts from the sum of the variances, making the variance of the difference possibly smaller than is obtained from the fixed-effects model. In fact this often is the case when you have a random block effect and there is substantial variation between the blocks. In such cases, accounting for blocks is a very powerful technique because the between-block variation cancels out when you look at within-block differences.