str(AY_rt)
## 'data.frame': 3116 obs. of 6 variables:
## $ subID : Factor w/ 164 levels "sub01","sub04",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ groupID: Factor w/ 2 levels "AE","UAA": 2 2 2 2 2 2 2 2 2 2 ...
## $ timeID : Factor w/ 2 levels "post","pre": 2 2 2 2 2 2 2 2 2 2 ...
## $ trialID: int 4 30 36 43 57 66 68 89 93 105 ...
## $ condnID: Factor w/ 1 level "AY": 1 1 1 1 1 1 1 1 1 1 ...
## $ probeRT: int 1191 919 542 599 1115 779 767 614 986 683 ...
fit <- lmer(probeRT ~ groupID * timeID + (1|trialID), data = AY_rt, REML = TRUE)
# REML = TRUE following the link and advice here: https://stats.stackexchange.com/questions/48671/what-is-restricted-maximum-likelihood-and-when-should-it-be-used (Reference #8)
qqpoints<- qqnorm(resid(fit))
qqline(resid(fit))
shapiro.test(resid(fit))
##
## Shapiro-Wilk normality test
##
## data: resid(fit)
## W = 0.9313, p-value < 2.2e-16
#non-normal residuals
Fails Shapiro-Wilkes test. Following advice from reference 13, we have the following GLMM model:
summary(GHQ <- glmer(probeRT ~ groupID + timeID + trialID+groupID*timeID*trialID + (1+groupID*timeID|subID) + (1|trialID), data = AY_rt,family = Gamma(link="identity"), nAGQ = 0))
## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
## Family: Gamma ( identity )
## Formula: probeRT ~ groupID + timeID + trialID + groupID * timeID * trialID +
## (1 + groupID * timeID | subID) + (1 | trialID)
## Data: AY_rt
##
## AIC BIC logLik deviance df.resid
## 41388.8 41509.6 -20674.4 41348.8 3096
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8002 -0.5988 -0.1165 0.4611 4.9158
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subID (Intercept) 8.631e+03 92.9013
## groupIDUAA 7.122e+01 8.4390 -1.00
## timeIDpre 1.137e+04 106.6339 -0.43 0.43
## groupIDUAA:timeIDpre 1.725e+03 41.5293 0.07 -0.07 -0.76
## trialID (Intercept) 1.038e+03 32.2130
## Residual 4.915e-02 0.2217
## Number of obs: 3116, groups: subID, 164; trialID, 120
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 845.14670 18.59450 45.451 <2e-16 ***
## groupIDUAA -36.85759 24.14073 -1.527 0.127
## timeIDpre 19.01885 22.95115 0.829 0.407
## trialID 0.47808 0.21896 2.183 0.029 *
## groupIDUAA:timeIDpre -14.78429 30.97153 -0.477 0.633
## groupIDUAA:trialID 0.19580 0.28259 0.693 0.488
## timeIDpre:trialID -0.06861 0.28347 -0.242 0.809
## groupIDUAA:timeIDpre:trialID -0.16343 0.39465 -0.414 0.679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) grIDUAA tmIDpr trilID grpIDUAA:tmID grpIDUAA:trID tID:ID
## groupIDUAA -0.692
## timeIDpre -0.605 0.465
## trialID -0.703 0.458 0.481
## grpIDUAA:tmID 0.448 -0.648 -0.739 -0.356
## grpIDUAA:trID 0.461 -0.683 -0.371 -0.657 0.531
## tmIDpr:trID 0.457 -0.351 -0.717 -0.654 0.530 0.505
## gIDUAA:ID:I -0.328 0.487 0.513 0.469 -0.737 -0.714 -0.716
#getME(GHQ,'b')