Why is converting emmeans contrasts to a data.frame not reporting correct p-values?

Question

The p-values for the contrasts I am running are not being converted correctly to a data.frame. Why is this and how do I fix it?

Console output for emmeans:

> pairs(emmeans(lmer.mod, ~ Status*Stim*Treatment), simple = "each")
$`simple contrasts for Status`
Stim = 1, Treatment = None:
 contrast               estimate     SE   df t.ratio p.value
 Control - Subclinical  -0.24213 0.0571 57.5 -4.241  0.0002 
 Control - Clinical     -0.16275 0.0571 57.5 -2.851  0.0164 
 Subclinical - Clinical  0.07938 0.0571 57.5  1.390  0.3526 

Console output for data.frame of emmeans:

> mod.EMM <- pairs(emmeans(lmer.mod, ~ Status*Stim*Treatment), simple = "each")
> as.data.frame(mod.EMM)
   Stim Treatment      Status               contrast     estimate         SE       df      t.ratio      p.value
1    1      None           .  Control - Subclinical -0.242125000 0.05709000 57.46544  -4.24111052 3.680551e-03
2    1      None           .     Control - Clinical -0.162750000 0.05709000 57.46544  -2.85076195 2.721389e-01
3    1      None           . Subclinical - Clinical  0.079375000 0.05709000 57.46544   1.39034857 1.000000e+00

Reproducible example:

model1 <- lm(uptake ~ Type + Treatment + conc + Type*Treatment, data=CO2)

library(emmeans)
pairs(emmeans(model1, ~ Type*Treatment), simple="each")
# $`simple contrasts for Type`
# Treatment = nonchilled:
#   contrast             estimate   SE df t.ratio p.value
# Quebec - Mississippi     9.38 1.85 79 5.068   <.0001 
# 
# Treatment = chilled:
#   contrast             estimate   SE df t.ratio p.value
# Quebec - Mississippi    15.94 1.85 79 8.610   <.0001 
# 
# 
# $`simple contrasts for Treatment`
# Type = Quebec:
#   contrast             estimate   SE df t.ratio p.value
# nonchilled - chilled     3.58 1.85 79 1.934   0.0566 
# 
# Type = Mississippi:
#   contrast             estimate   SE df t.ratio p.value
# nonchilled - chilled    10.14 1.85 79 5.477   <.0001 

as.data.frame(pairs(emmeans(model1, ~ Type*Treatment), simple="each"))
#    Treatment        Type             contrast  estimate       SE df  t.ratio      p.value
# 1 nonchilled           . Quebec - Mississippi  9.380952 1.851185 79 5.067538 1.036140e-05
# 2    chilled           . Quebec - Mississippi 15.938095 1.851185 79 8.609670 2.252161e-12
# 3          .      Quebec nonchilled - chilled  3.580952 1.851185 79 1.934410 2.265719e-01
# 4          . Mississippi nonchilled - chilled 10.138095 1.851185 79 5.476542 1.995066e-06

model1 <- lm(uptake ~ Type + Treatment + conc + Type*Treatment, data=CO2)
pairs(emmeans(model1, ~ Type*Treatment), simple="each")
# $`simple contrasts for Type`
# Treatment = nonchilled:
#   contrast             estimate   SE df t.ratio p.value
# Quebec - Mississippi     9.38 1.85 79 5.068   <.0001 
# 
# Treatment = chilled:
#   contrast             estimate   SE df t.ratio p.value
# Quebec - Mississippi    15.94 1.85 79 8.610   <.0001 
# 
# 
# $`simple contrasts for Treatment`
# Type = Quebec:
#   contrast             estimate   SE df t.ratio p.value
# nonchilled - chilled     3.58 1.85 79 1.934   0.0566 
# 
# Type = Mississippi:
#   contrast             estimate   SE df t.ratio p.value
# nonchilled - chilled    10.14 1.85 79 5.477   <.0001 

as.data.frame(pairs(emmeans(model1, ~ Type*Treatment), simple="each"))
#    Treatment        Type             contrast  estimate       SE df  t.ratio      p.value
# 1 nonchilled           . Quebec - Mississippi  9.380952 1.851185 79 5.067538 1.036140e-05
# 2    chilled           . Quebec - Mississippi 15.938095 1.851185 79 8.609670 2.252161e-12
# 3          .      Quebec nonchilled - chilled  3.580952 1.851185 79 1.934410 2.265719e-01
# 4          . Mississippi nonchilled - chilled 10.138095 1.851185 79 5.476542 1.995066e-06

Update from outside help:

"It seems like the result of pairs() is not itself an emmGrid object that can be turned into a data frame, but a list containing two emmGrid objects. If you extract either of those objects by position from the list, using [[]], like so,

pairs(emmeans(model1, ~ Type*Treatment), simple = "each")[[2]]

then you can data.frame() each result and it will be correct. You end up with two different dataframes to hold the contrasts involving the two different variables, but each of these dataframes has the correct p-values in it."

I'm hoping someone has a better work-around for this issue so I can combine all the contrasts into one data.frame.

Answer 1

The different p-values you are seeing reflect unadjusted p-values vs p-values that were adjusted for multiple comparisons.

The ?emmeans::pairs documentation tells us:

Ordinarily, when simple is a list or "each", the return value is an emm_list object with each entry in correspondence with the entries of simple. However, with combine = TRUE, the elements are all combined into one family of contrasts in a single emmGrid object using rbind.emmGrid.. In that case, the adjust argument sets the adjustment method for the combined set of contrasts.

So, with your reproducible example, you can combine all the simple main effects into one data frame with the combine argument set to TRUE. And you can choose between unadjusted vs adjusted p-values by setting the adjust argument.

model1 <- lm(uptake ~ Type + Treatment + conc + Type*Treatment, data=CO2)

> pairs(emmeans(model1, ~ Type*Treatment), simple = "each", combine = TRUE,
+               adjust = "none")
 Treatment  Type        contrast             estimate   SE df t.ratio p.value
 nonchilled .           Quebec - Mississippi     9.38 1.85 79 5.068   <.0001 
 chilled    .           Quebec - Mississippi    15.94 1.85 79 8.610   <.0001 
 .          Quebec      nonchilled - chilled     3.58 1.85 79 1.934   0.0566 
 .          Mississippi nonchilled - chilled    10.14 1.85 79 5.477   <.0001 

Here's one with Bonferroni adjustment:

> pairs(emmeans(model1, ~ Type*Treatment), simple = "each", combine = TRUE,
+               adjust = "bonferroni")
 Treatment  Type        contrast             estimate   SE df t.ratio p.value
 nonchilled .           Quebec - Mississippi     9.38 1.85 79 5.068   <.0001 
 chilled    .           Quebec - Mississippi    15.94 1.85 79 8.610   <.0001 
 .          Quebec      nonchilled - chilled     3.58 1.85 79 1.934   0.2266 
 .          Mississippi nonchilled - chilled    10.14 1.85 79 5.477   <.0001 

P value adjustment: bonferroni method for 4 tests