I want to test the effects of different variables on Reaction Times data. Following the recommendations of Lo & Andrews (2015) I compared the AIC/BIC of three GLMM with a Gaussian, a Gamma and an Inverse Gaussian distribution with an identity link (with the lme4 package)
Model_Mn_VAAST<- lmer(VAAST.RT~Conditionc* Type_Bloc_Encounteredc + (1|pp) + (1|Perso), Pilot5_FINAL_Mn2)
Model_Mn_VAAST_gamma<- glmer(VAAST.RT~Conditionc* Type_Bloc_Encounteredc + (1|pp) + (1|Perso), Pilot5_FINAL_Mn2,
family = Gamma(link = "identity"),
glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1000000)))
Model_Mn_VAAST_invGauss<- glmer(VAAST.RT~Conditionc* Type_Bloc_Encounteredc + (1|pp) + (1|Perso), Pilot5_FINAL_Mn2, family = inverse.gaussian(link = "identity"))
The model that best fits my data (with the lower AIC) is the Inverse Gaussian GLMM.
Model_Mn_VAAST: VAAST.RT ~ Conditionc * Type_Bloc_Encounteredc + (1 | pp) + (1 | Perso)
Model_Mn_VAAST_gamma: VAAST.RT ~ Conditionc * Type_Bloc_Encounteredc + (1 | pp) + (1 | Perso)
Model_Mn_VAAST_invGauss: VAAST.RT ~ Conditionc * Type_Bloc_Encounteredc + (1 | pp) + (1 | Perso)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
Model_Mn_VAAST 7 335238 335294 -167612 335224
Model_Mn_VAAST_gamma 7 328211 328268 -164099 328197 7026.4 0 < 2.2e-16 ***
Model_Mn_VAAST_invGauss 7 326411 326468 -163199 326397 1799.9 0 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Ok. I run the Inverse Gaussian GLMM and obtained the following output.
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [ glmerMod]
Family: inverse.gaussian ( identity )
Formula: VAAST.RT ~ Conditionc * Type_Bloc_Encounteredc + (1 | pp) + (1 | Perso)
Data: Pilot5_FINAL_Mn2
AIC BIC logLik deviance df.resid
326411.3 326468.2 -163198.7 326397.3 24842
Scaled residuals:
Min 1Q Median 3Q Max
-2.1492 -0.6180 -0.2170 0.3246 9.4115
Random effects:
Groups Name Variance Std.Dev.
pp (Intercept) 3.324e+03 57.650995
Perso (Intercept) 7.707e+01 8.779221
Residual 7.768e-05 0.008814
Number of obs: 24849, groups: pp, 201; Perso, 64
Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept) 856.078 6.629 129.145 <2e-16 ***
Conditionc 14.450 6.144 2.352 0.0187 *
Type_Bloc_Encounteredc -5.169 2.062 -2.507 0.0122 *
Conditionc:Type_Bloc_Encounteredc 6.647 3.378 1.968 0.0491 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Cndtnc Ty_B_E
Conditionc -0.109
Typ_Blc_Enc -0.023 0.055
Cndtn:T_B_E -0.110 0.080 0.038
I am a little bit confused concerning the obtained p-values.
Indeed, if I compare models with and without any of the fixed parameters (with the anova() finction) I did not obtain the same results.
For instance, concerning the interaction term I found the following results:
Models:
Model_Mn_VAAST_invGaussInt: VAAST.RT ~ Conditionc + Type_Bloc_Encounteredc + (1 | pp) + (1 | Perso)
Model_Mn_VAAST_invGauss: VAAST.RT ~ Conditionc * Type_Bloc_Encounteredc + (1 | pp) + (1 | Perso)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
Model_Mn_VAAST_invGaussInt 6 326412 326460 -163200 326400
Model_Mn_VAAST_invGauss 7 326411 326468 -163199 326397 2.3283 1 0.127
If I am right, from the summary() function the p-values are calculated based on Wald tests while there are calculated based on likelihood ratio tests with anova().
Which one is the best method ?
And how can I get the confidence intervals for the parameters estimates ? The only way I know is
confint(Model_Mn_VAAST_invGauss, method = "Wald")
But this doesn't apply if I choose the anova function for the p-values, no ?
Thank you for any help or information !
