I used glmer to analyze reaction times in a repeated-measures design with 2 faixed factors (3 levels each) and 1 random factor. As I have a significant interaction, I would like to do post-hoc to get contrasts. I first used relevel to change the reference of the model as many times as needed to get the contrasts I needed. I have seen that glht can also be used (and avoid relevelling) but I get some p-values that do not make any sense when looking at the data. I tried relevelling one of my factor just before glht and the p-values make more sense. Is my contrast matrix OK? Does this mean I need to use relevel also with glht? Can anyone help me please?
I have an experiment with 2 repeated fixed factors Verb (vt,vnt,vnta) and Delay (170,350,500), and Subjects as the random factor. As I measured reaction times, I used glmer (lme4 package) following the recommendation by Lo & Andrews (2015). The ANOVA of the model gives a significant interaction between the 2 fixed factors and I would to run post-hoc tests. I first used relevel to change the reference of the model and test all needed contrasts. I get significant differences between some of the conditions, however I'm not sure any correction is applied to the data. This is my first question.
I have seen that glht (multcomp package) can also be used for post-hoc (and avoids relevel). I therefore tried this but I'm not sure I'm doing it well as I get significant differences that look very weird when looking at the data. When I used relevel just before glht, then I get some more coherent results (and similar estimates as in the first relevel analysis). Does this mean that glht needs to use relevel sometimes? Can anyone help me?
> m2b <- glmer(RT ~ Delai * Verbe + (1|Sujet), data = Data_Tacti2, family=Gamma(link="identity"), glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
> Anova(m2b)
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: RT
Chisq Df Pr(>Chisq)
Delai 677.379 2 <2e-16 ***
Verbe 3.562 2 0.1685
Delai:Verbe 10.161 4 0.0378 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
model with reference VNT 170
> summary(m2b)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Gamma ( identity )
Formula: RT ~ Delai * Verbe + (1 | Sujet)
Data: Data_Tacti2
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
AIC BIC logLik deviance df.resid
51151.3 51220.9 -25564.6 51129.3 4122
Scaled residuals:
Min 1Q Median 3Q Max
-3.3223 -0.6339 -0.1768 0.4510 6.8590
Random effects:
Groups Name Variance Std.Dev.
Sujet (Intercept) 836.62309 28.9244
Residual 0.08911 0.2985
Number of obs: 4133, groups: Sujet, 18
Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept) 494.110 9.414 52.485 < 2e-16 ***
Delai350 -76.068 4.328 -17.576 < 2e-16 ***
Delai500 -93.102 4.313 -21.585 < 2e-16 ***
Verbevnta -7.559 4.901 -1.542 0.12298
Verbevt -16.837 5.224 -3.223 0.00127 **
Delai350:Verbevnta 7.926 6.657 1.191 0.23379
Delai500:Verbevnta 8.381 6.271 1.337 0.18136
Delai350:Verbevt 19.670 7.014 2.804 0.00504 **
Delai500:Verbevt 10.198 5.624 1.813 0.06980 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Del350 Del500 Vrbvnt Verbvt Dl350:Vrbvn Dl500:Vrbvn Dl350:Vrbvt
Delai350 -0.262
Delai500 -0.212 0.381
Verbevnta -0.232 0.206 0.139
Verbevt -0.297 0.196 0.168 0.184
Dl350:Vrbvn 0.209 -0.360 -0.019 -0.555 0.024
Dl500:Vrbvn 0.280 -0.182 -0.439 -0.515 -0.051 0.351
Dl350:Vrbvt 0.297 -0.408 -0.066 -0.011 -0.576 0.101 0.038
Dl500:Vrbvt 0.195 -0.083 -0.350 0.052 -0.571 -0.106 0.131 0.406
relevelling to VT 170
> Data_Tacti2$Delai <- relevel (Data_Tacti2$Delai, ref = "170")
> Data_Tacti2$Verbe <- relevel (Data_Tacti2$Verbe, ref = "vt")
> m2b_vt_170 <- glmer(RT ~ Delai * Verbe + (1|Sujet), data = Data_Tacti2, family=Gamma(link="identity"), glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
> summary(m2b_vt_170)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Gamma ( identity )
Formula: RT ~ Delai * Verbe + (1 | Sujet)
Data: Data_Tacti2
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
AIC BIC logLik deviance df.resid
51151.3 51220.9 -25564.6 51129.3 4122
Scaled residuals:
Min 1Q Median 3Q Max
-3.3223 -0.6339 -0.1768 0.4510 6.8590
Random effects:
Groups Name Variance Std.Dev.
Sujet (Intercept) 836.63773 28.9247
Residual 0.08911 0.2985
Number of obs: 4133, groups: Sujet, 18
Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept) 477.272 7.337 65.050 < 2e-16 ***
Delai350 -56.397 4.444 -12.692 < 2e-16 ***
Delai500 -82.904 4.670 -17.753 < 2e-16 ***
Verbevnt 16.838 4.189 4.019 5.83e-05 ***
Verbevnta 9.279 5.293 1.753 0.07962 .
Delai350:Verbevnt -19.671 5.984 -3.287 0.00101 **
Delai500:Verbevnt -10.198 7.293 -1.398 0.16200
Delai350:Verbevnta -11.745 6.581 -1.785 0.07433 .
Delai500:Verbevnta -1.817 5.805 -0.313 0.75428
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Del350 Del500 Verbevnt Verbevnta Delai350:Verbevnt Delai500:Verbevnt Delai350:Verbevnta
Delai350 -0.114
Delai500 -0.015 0.276
Verbevnt 0.012 -0.010 0.108
Verbevnta -0.157 0.220 0.070 0.139
Delai350:Verbevnt 0.006 -0.226 0.095 -0.420 0.091
Delai500:Verbevnt 0.003 0.107 -0.456 -0.438 0.152 0.183
Delai350:Verbevnta 0.198 -0.428 0.006 0.101 -0.570 0.007 -0.177
Delai500:Verbevnta 0.062 -0.038 -0.332 0.043 -0.511 -0.173 0.100 0.342
relevelling to VNT 350
> Data_Tacti2$Delai <- relevel (Data_Tacti2$Delai, ref = "350")
> Data_Tacti2$Verbe <- relevel (Data_Tacti2$Verbe, ref = "vnt")
> m2b_vnt_350 <- glmer(RT ~ Delai * Verbe + (1|Sujet), data = Data_Tacti2, family=Gamma(link="identity"), glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
> summary(m2b_vnt_350)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Gamma ( identity )
Formula: RT ~ Delai * Verbe + (1 | Sujet)
Data: Data_Tacti2
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
AIC BIC logLik deviance df.resid
51151.3 51220.9 -25564.6 51129.3 4122
Scaled residuals:
Min 1Q Median 3Q Max
-3.3223 -0.6339 -0.1768 0.4510 6.8590
Random effects:
Groups Name Variance Std.Dev.
Sujet (Intercept) 836.62212 28.9244
Residual 0.08911 0.2985
Number of obs: 4133, groups: Sujet, 18
Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept) 420.875 7.076 59.478 < 2e-16 ***
Delai170 56.398 5.221 10.802 < 2e-16 ***
Delai500 -26.507 4.087 -6.486 8.8e-11 ***
Verbevnt -2.833 4.184 -0.677 0.498329
Verbevnta -2.466 4.030 -0.612 0.540484
Delai170:Verbevnt 19.671 5.219 3.769 0.000164 ***
Delai500:Verbevnt 9.473 5.672 1.670 0.094900 .
Delai170:Verbevnta 11.745 5.454 2.153 0.031280 *
Delai500:Verbevnta 9.928 5.097 1.948 0.051439 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Del170 Del500 Verbevnt Verbevnta Delai170:Verbevnt Delai500:Verbevnt Delai170:Verbevnta
Delai170 0.032
Delai500 0.107 0.225
Verbevnt -0.022 -0.140 0.117
Verbevnta -0.010 -0.035 0.073 0.235
Delai170:Verbevnt 0.041 -0.310 -0.045 -0.229 0.066
Delai500:Verbevnt -0.096 0.200 -0.356 -0.468 0.029 0.137
Delai170:Verbevnta -0.087 -0.438 0.037 0.181 -0.115 0.115 -0.189
Delai500:Verbevnta -0.105 0.094 -0.336 -0.045 -0.383 -0.039 0.184 -0.019
I used relevel as many times as needed to get all the contrasts between the 3 levels of factor Verb and 3 levels of factor Delay (not reported here)
now using glht
model m2b with reference VT 170
> contrast.matrix <- rbind(
+ `VNT170 vs VT170` = c(0, 0, 0, 0, 1, 0, 0, 0, 0),
+ `VNTA170 vs VT170` = c(0, 0, 0, -1, 1, 0, 0, 0, 0),
+ `VNT170 vs VNTA170` = c(0, 0, 0, 1, 0, 0, 0, 0, 0),
+ `VNT350 vs VT350` = c(0, 1, 0, 0, 0, 0, 0, -1, 0),
+ `VNTA350 vs VT350` = c(0, 0, 0, 0, 0, 1, 0, -1, 0),
+ `VNT350 vs VNTA350` = c(0, 1, 0, 0, 0, -1, 0, 0, 0),
+ `VNT500 vs VT500` = c(0, 0, 1, 0, 0, 0, 0, 0, -1),
+ `VNTA500 vs VT500` = c(0, 0, 0, 0, 0, 0, 1, 0, -1),
+ `VNT500 vs VNTA500` = c(0, 0, 1, 0, 0, 0, -1, 0, 0)
+ )
> comps_refVNT170 <- glht(m2b, contrast.matrix)
> summary(comps_refVNT170)
Simultaneous Tests for General Linear Hypotheses
Fit: glmer(formula = RT ~ Delai * Verbe + (1 | Sujet), data = Data_Tacti2,
family = Gamma(link = "identity"), control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 1e+05)))
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
VNT170 vs VT170 == 0 -16.837 5.224 -3.223 0.010 *
VNTA170 vs VT170 == 0 -9.278 6.472 -1.434 0.621
VNT170 vs VNTA170 == 0 -7.559 4.901 -1.542 0.544
VNT350 vs VT350 == 0 -95.738 9.629 -9.943 <0.001 ***
VNTA350 vs VT350 == 0 -11.744 9.171 -1.281 0.726
VNT350 vs VNTA350 == 0 -83.994 9.152 -9.178 <0.001 ***
VNT500 vs VT500 == 0 -103.300 8.198 -12.601 <0.001 ***
VNTA500 vs VT500 == 0 -1.816 7.856 -0.231 1.000
VNT500 vs VNTA500 == 0 -101.484 9.038 -11.228 <0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
I get p-value <.001 for the contrast VNT350-VT350 which is very weird (same for VNT500-VT500 and VNT500-VNTA500 when looking at the data), and this difference is not significant in the first relevel analysis
Is the matrix OK?
if I relevel tp reference to VNT 350 (although I'm not supposed to) the estimate is similar to that in the first relevel analysis and the p-value for VNT350-VT350 makes more sense
Does this mean I need to relevel also for glht?
> Data_Tacti2$Delai <- relevel (Data_Tacti2$Delai, ref = "350")
> Data_Tacti2$Verbe <- relevel (Data_Tacti2$Verbe, ref = "vnt")
> m2b_vnt_350 <- glmer(RT ~ Delai * Verbe + (1|Sujet), data = Data_Tacti2, family=Gamma(link="identity"), glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
> summary(m2b_vnt_350)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Gamma ( identity )
Formula: RT ~ Delai * Verbe + (1 | Sujet)
Data: Data_Tacti2
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
AIC BIC logLik deviance df.resid
51151.3 51220.9 -25564.6 51129.3 4122
Scaled residuals:
Min 1Q Median 3Q Max
-3.3223 -0.6339 -0.1768 0.4510 6.8590
Random effects:
Groups Name Variance Std.Dev.
Sujet (Intercept) 836.62212 28.9244
Residual 0.08911 0.2985
Number of obs: 4133, groups: Sujet, 18
Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept) 420.875 7.076 59.478 < 2e-16 ***
Delai170 56.398 5.221 10.802 < 2e-16 ***
Delai500 -26.507 4.087 -6.486 8.8e-11 ***
Verbevnt -2.833 4.184 -0.677 0.498329
Verbevnta -2.466 4.030 -0.612 0.540484
Delai170:Verbevnt 19.671 5.219 3.769 0.000164 ***
Delai500:Verbevnt 9.473 5.672 1.670 0.094900 .
Delai170:Verbevnta 11.745 5.454 2.153 0.031280 *
Delai500:Verbevnta 9.928 5.097 1.948 0.051439 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Del170 Del500 Verbevnt Verbevnta Delai170:Verbevnt Delai500:Verbevnt Delai170:Verbevnta
Delai170 0.032
Delai500 0.107 0.225
Verbevnt -0.022 -0.140 0.117
Verbevnta -0.010 -0.035 0.073 0.235
Delai170:Verbevnt 0.041 -0.310 -0.045 -0.229 0.066
Delai500:Verbevnt -0.096 0.200 -0.356 -0.468 0.029 0.137
Delai170:Verbevnta -0.087 -0.438 0.037 0.181 -0.115 0.115 -0.189
Delai500:Verbevnta -0.105 0.094 -0.336 -0.045 -0.383 -0.039 0.184 -0.019
> contrast.matrix2 <- rbind(
+ `VNT170 vs VT170` = c(0, 1, 0, 0, 0, 0, 0, -1, 0),
+ `VNTA170 vs VT170` = c(0, 0, 0, 0, 0, 1, 0, -1, 0),
+ `VNT170 vs VNTA170` = c(0, 1, 0, 0, 0, -1, 0, 0, 0),
+ `VNT350 vs VT350` = c(0, 0, 0, 0, 1, 0, 0, 0, 0),
+ `VNTA350 vs VT350` = c(0, 0, 0, 1, -1, 0, 0, 0, 0),
+ `VNT350 vs VNTA350` = c(0, 0, 0, 1, 0, 0, 0, 0, 0),
+ `VNT500 vs VT500` = c(0, 0, 1, 0, 0, 0, 0, 0, -1),
+ `VNTA500 vs VT500` = c(0, 0, 0, 0, 0, 0, 1, 0, -1),
+ `VNT500 vs VNTA500` = c(0, 0, 1, 0, 0, 0, -1, 0, 0)
+ )
> comps_refVNT350 <- glht(m2b_vnt_350, contrast.matrix2)
> summary(comps_refVNT350)
Simultaneous Tests for General Linear Hypotheses
Fit: glmer(formula = RT ~ Delai * Verbe + (1 | Sujet), data = Data_Tacti2,
family = Gamma(link = "identity"), control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 1e+05)))
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
VNT170 vs VT170 == 0 44.6531 9.0536 4.932 < 1e-04 ***
VNTA170 vs VT170 == 0 7.9262 7.0999 1.116 0.858353
VNT170 vs VNTA170 == 0 36.7269 8.4492 4.347 0.000132 ***
VNT350 vs VT350 == 0 -2.4665 4.0296 -0.612 0.991632
VNTA350 vs VT350 == 0 -0.3665 5.0804 -0.072 1.000000
VNT350 vs VNTA350 == 0 -2.8330 4.1839 -0.677 0.985839
VNT500 vs VT500 == 0 -36.4355 7.5282 -4.840 < 1e-04 ***
VNTA500 vs VT500 == 0 -0.4557 6.8929 -0.066 1.000000
VNT500 vs VNTA500 == 0 -35.9798 8.0848 -4.450 < 1e-04 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
