3
$\begingroup$

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)

# read in data
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 
#######
$\endgroup$

2 Answers 2

1
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – chaoh
    Jan 30, 2020 at 20:42
1
$\begingroup$

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       
$\endgroup$

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.