2
$\begingroup$

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,
    adjust = "mvt",
  )

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

> 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

emmeans::emmip from linear and mixed effect models with significant contrasts indicated

And emmeans::pwpp as another visual

emmeans::pwpp from linear and mixed effect models

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

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

I will answer briefly in general terms.

Where do the contrast SEs come from?

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}}$.

Why are they smaller for the mixed model than for the linear model [in some cases]?

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for your response. Am I understanding this correctly? In my case the variation between participants (rather than blocks) is large, so the variance of the estimated means are large. But the covariance between means will be highly correlated, so effectively this cancels out when comparing the within-participant difference related to Condition? I guess that makes sense at part of the reason to model participants as a random effect in the first place. $\endgroup$
    – Jem Arnold
    Commented Aug 5, 2023 at 15:21
  • $\begingroup$ Yes but even if they were modeled as fixed effects, you still get about the same test statistics for the within-participant effects. Between-participant differences will be much more significant though, because with fixed participants, you narrow the scope of the inference to just those participants you have observed, rather than a population of participants. $\endgroup$
    – Russ Lenth
    Commented Sep 5, 2023 at 21:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.