I am running a GLMM using the glmer function from the lme4 package in R. I would like to use a Helmert contrast to test a specific difference for the factor condition. If I run that in the glmer function, that works well:
contrasts(data_all$condition_f) = contr.helmert(3)
model_vector <- glmer(cbind(normCor, normIncor) ~
condition_f*sorting_f*time_f + rescaled_VM + rescaled_L +
(1+condition_f + time_f + sorting_f|id) + (1+time_f|word_id),
data = data_all, family = binomial(link = 'logit'),
glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
summary(model_vector)
results in:
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: cbind(normCor, normIncor) ~ condition_f * sorting_f * time_f
+ rescaled_VM + rescaled_L + (1 + condition_f + time_f +
sorting_f | id) + (1 + time_f | word_id)
Data: data_all
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
AIC BIC logLik deviance df.resid
35224.6 35450.3 -17580.3 35160.6 8509
Scaled residuals:
Min 1Q Median 3Q Max
-13.5117 -1.1089 -0.1801 1.0330 8.2894
Random effects:
Groups Name Variance Std.Dev. Corr
word_id (Intercept) 0.72186 0.8496
time_fretention 0.22858 0.4781 -0.75
id (Intercept) 1.68692 1.2988
condition_f1 0.02358 0.1536 -0.11
condition_f2 0.01460 0.1208 0.16 0.01
time_fretention 0.91416 0.9561 -0.68 -0.03 0.08
sorting_fsorted 0.83799 0.9154 -0.46 0.11 -0.37 0.13
Number of obs: 8541, groups: word_id, 73; id, 59
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.206778 0.254561 0.812 0.41663
condition_f1 -0.070279 0.127098 -0.553 0.58030
condition_f2 -0.006019 0.075334 -0.080 0.93632
sorting_fsorted 0.448982 0.326037 1.377 0.16848
time_fretention -1.174478 0.183483 -6.401 1.54e-10 ***
rescaled_VM -0.008824 0.131198 -0.067 0.94637
rescaled_L 0.255718 0.142523 1.794 0.07278 .
condition_f1:sorting_fsorted 0.168045 0.053817 3.123 0.00179 **
condition_f2:sorting_fsorted 0.099171 0.038233 2.594 0.00949 **
condition_f1:time_fretention 0.146651 0.077239 1.899 0.05761 .
condition_f2:time_fretention -0.035498 0.045058 -0.788 0.43080
sorting_fsorted:time_fretention -0.066772 0.253215 -0.264 0.79201
condition_f1:sorting_fsorted:time_fretention -0.110543 0.050520 -2.188 0.02866 *
condition_f2:sorting_fsorted:time_fretention -0.082406 0.030258 -2.723 0.00646 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Where the two Helmert levels can be seen in condition_f1 & condition_f2. Note the significant interactions with sorting (as well as the three-way interactions).
However, when I want to run the car Anova function to calculate more accurate p-values, it does not show the Helmert contrasts anymore:
Anova(model_vector, test = "Chi")
results in:
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: cbind(normCor, normIncor)
Chisq Df Pr(>Chisq)
condition_f 3.0183 2 0.221102
sorting_f 2.9130 1 0.087870 .
time_f 75.3887 1 < 2.2e-16 ***
rescaled_VM 0.0045 1 0.946375
rescaled_L 3.2192 1 0.072778 .
condition_f:sorting_f 8.2811 2 0.015914 *
condition_f:time_f 4.7222 2 0.094318 .
sorting_f:time_f 0.0663 1 0.796731
condition_f:sorting_f:time_f 11.8321 2 0.002696 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
As can be seen, now there is no separation in the Helmert contrast anymore, while the relevant interactions remain significant.
I also tried to calculate the contrasts using emmeans post-hoc, by not specifying the contrast in the model and, after running the model, run:
helmert.emmc <- function(levs, ...) {
M <- as.data.frame(contr.helmert(levs))
names(M) <- paste(levs[-1],"vs earlier")
attr(M, "desc") <- "Helmert contrasts"
M
}
model_vector_emmeans <- emmeans(model_vector, ~condition_f|sorting_f)
contrast(model_vector_emmeans, "helmert")
This gives as output
sorting_f = shuffled:
contrast estimate SE df z.ratio p.value
Somewhat similar vs earlier 0.0743 0.213 Inf 0.350 0.7267
Dissimilar vs earlier 0.0621 0.366 Inf 0.170 0.8651
sorting_f = sorted:
contrast estimate SE df z.ratio p.value
Somewhat similar vs earlier 0.0132 0.214 Inf 0.061 0.9510
Dissimilar vs earlier -0.4501 0.368 Inf -1.224 0.2210
This gives four p values, in line with the Helmert contrast per level of sorting, but does not give me the overall p-value for the interaction.
Instead running
model_vector_emmeans <- emmeans(model_vector, ~condition_f*sorting_f)
contrast(model_vector_emmeans, "helmert")
Gives me five levels, of which I am not sure where they come from, I think from creating the interaction (2 * 3 levels) then applying the contrast (resulting in n-1 = 5 levels), also not what I intended to do:
contrast estimate SE df z.ratio p.value
Somewhat similar shuffled vs earlier 0.0743 0.213 Inf 0.350 0.7267
Dissimilar shuffled vs earlier 0.0621 0.366 Inf 0.170 0.8651
Very similar sorted vs earlier 1.4521 0.863 Inf 1.682 0.0925
Somewhat similar sorted vs earlier 1.5047 0.968 Inf 1.554 0.1201
Dissimilar sorted vs earlier 0.3466 1.089 Inf 0.318 0.7502
Hence, my question is: how can I calculate accurate significance values for Helmert contrasts using the R glmer function, using Anova, emmeans or some other means?
