# In R, when summary() and emmeans() give different results, which one should I trust?

I want to compare scores in the "control" condition to the "high" condition and to the "low" condition. The summary() and the emmeans() functions give different significance results for the "high" vs "control" contrast: the "high" vs "control" contrast is significant in the result given by summary() but not emmeans(). Why do they give different results? Which one should I trust?

Here is the dataset.

Here is code I ran to get the results:

# load library
library(tidyverse)
library(lmerTest)
library(emmeans)

df.mod <- df.csv

# set contrasts to compare "high" to "control", and compare "low" to "control" using summary()
highVScontrol <- c(-1/2, 1/2, 0)
lowVScontrol <- c(-1/2, 0, 1/2)
contrasts(df.mod$condition) <- cbind(highVScontrol, lowVScontrol) # fit model mod <- lmer(score ~ condition + (1|subject), data = df.mod, REML = FALSE) summary(mod) ######### # Output # Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [ # lmerModLmerTest] # Formula: score ~ condition + (1 | subject) # Data: df.mod # # AIC BIC logLik deviance df.resid # 283.0 295.5 -136.5 273.0 85 # # Scaled residuals: # Min 1Q Median 3Q Max # -3.09796 -0.51784 0.06943 0.64971 2.46609 # # Random effects: # Groups Name Variance Std.Dev. # subject (Intercept) 0.2618 0.5116 # Residual 1.0024 1.0012 # Number of obs: 90, groups: subject, 30 # # Fixed effects: # Estimate Std. Error df t value Pr(>|t|) # (Intercept) -1.0650 0.1409 30.0000 -7.557 2e-08 *** # conditionhighVScontrol 0.5342 0.2985 60.0000 1.790 0.0786 . # conditionlowVScontrol 0.6646 0.2985 60.0000 2.226 0.0298 * # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # # Correlation of Fixed Effects: # (Intr) cndtnhVS # cndtnhghVSc 0.000 # cndtnlwVScn 0.000 -0.500 ######### # comparison using emmeans() emmeans(mod, list(pairwise ~ condition), adjust = "holm") ######### # Output # $$emmeans of condition # condition emmean SE df lower.CL upper.CL # high -1.332 0.209 85.8 -1.747 -0.9170 # control -0.466 0.209 85.8 -0.881 -0.0506 # low -1.397 0.209 85.8 -1.812 -0.9822 # # Degrees-of-freedom method: kenward-roger # Confidence level used: 0.95 # #$$pairwise differences of condition # contrast estimate SE df t.ratio p.value # high - control -0.8664 0.263 62.1 -3.295 0.0033 # high - low 0.0652 0.263 62.1 0.248 0.8050 # control - low 0.9316 0.263 62.1 3.543 0.0023 # # P value adjustment: holm method for 3 tests #######  ## 2 Answers I think there might be something off about your contrast coding. This is not something I am terribly familiar with, so I ran the model without the contrast coding, and instead recoded condition into numerical values as such (note I called your dataset test): test <- test %>% mutate(condition2=ifelse(condition == "control", 0, ifelse(condition == "low", 1, 2)), condition2 = as.factor(condition2))  Then I ran the lmer model: mod <- lmer(score ~ condition2 + (1|subject), data = test, REML = FALSE)  with the following output: Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest'] Formula: score ~ condition2 + (1 | subject) Data: test AIC BIC logLik deviance df.resid 283.0 295.5 -136.5 273.0 85 Scaled residuals: Min 1Q Median 3Q Max -3.09796 -0.51784 0.06943 0.64971 2.46609 Random effects: Groups Name Variance Std.Dev. subject (Intercept) 0.2618 0.5116 Residual 1.0024 1.0012 Number of obs: 90, groups: subject, 30 Fixed effects: Estimate Std. Error df t value Pr(>|t|) (Intercept) -0.4657 0.2053 82.8910 -2.268 0.025901 * condition21 -0.9316 0.2585 60.0000 -3.604 0.000638 *** condition22 -0.8664 0.2585 60.0000 -3.352 0.001394 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Note that these coefficients are telling us about the mean difference between control and low (condition21) and control and high (condition22). Now emmeans: emmeans(mod, list(pairwise ~ condition2), adjust = "holm")$emmeans of condition2
condition2 emmean    SE   df lower.CL upper.CL
0          -0.466 0.209 85.8   -0.881  -0.0506
1          -1.397 0.209 85.8   -1.812  -0.9822
2          -1.332 0.209 85.8   -1.747  -0.9170

Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95

$pairwise differences of condition2 contrast estimate SE df t.ratio p.value 0 - 1 0.9316 0.263 62.1 3.543 0.0023 0 - 2 0.8664 0.263 62.1 3.295 0.0033 1 - 2 -0.0652 0.263 62.1 -0.248 0.8050  Notice how the 0-1 and 0-2 contrasts exactly match the output from lmer. This was not the case when comparing your lmer output to your emmeans output. What we get from emmeans is a direct test of the 1-2 contrast, which we did not get in lmer. The p-values are different, but not different enough that you would change your conclusion about statistical significance (i.e., if your alpha level is p<.05). I am not sure why they are different, but it has to do with the procedure difference in lmertest vs. emmeans. • Thank you, Erik! You reminded me of the problem of my contrast coding. Yes, now I realized the contrast coding I used is for orthogonal contrast (there is another error too but I corrected it), but I should have used the non-orthogonal contrast coding. I'll answer with the correct contrast coding. Jan 30, 2020 at 20:42 Now I realized the contrast coding in my question is for the orthogonal contrast. I should use the contrast coding for the non-orthogonal contrast. Here is the corrected contrast coding (assuming the condition levels are not changed): # specify contrast weights control_high <- c(-1/2, 1/2, 0) control_low <- c(-1/2, 0, 1/2) # create temporary matrix mat.temp <- rbind(constant=1/3, control_high, control_low) # get the inverse of that matrix mat <- solve(mat.temp) # drop the first column mat <- mat[, -1] # Now set the contrasts contrasts(df$condition) <- mat


Now the summary() gives the following output that matches the significance test given in emmeans() and @Erik Ruzek's answer:

# Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [lmerModLmerTest]
# Formula: score ~ condition + (1 | subject)
# Data: df
#
# AIC      BIC   logLik deviance df.resid
# 283.0    295.5   -136.5    273.0       85
#
# Scaled residuals:
# Min       1Q   Median       3Q      Max
# -3.09796 -0.51784  0.06943  0.64971  2.46609
#
# Random effects:
# Groups   Name        Variance Std.Dev.
# subject  (Intercept) 0.2618   0.5116
# Residual             1.0024   1.0012
# Number of obs: 90, groups:  subject, 30
#
# Fixed effects:
# Estimate Std. Error      df t value Pr(>|t|)
# (Intercept)            -1.0650     0.1409 30.0000  -7.557    2e-08 ***
# conditioncontrol_high  -0.4332     0.1293 60.0000  -3.352 0.001394 **
# conditioncontrol_low   -0.4658     0.1293 60.0000  -3.604 0.000638 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Correlation of Fixed Effects:
# (Intr) cndtncntrl_h
# cndtncntrl_h 0.000
# cndtncntrl_l 0.000  0.500