library(tidyverse)
library(numDeriv)
library(lme4)
library(MuMIn)
collotask <- read_csv("data/collotask2.csv")
collotask$Accuracy <- factor(collotask$Accuracy)
# Read_DE: doesn't make sense as continuous variable (three classes?)
table(collotask$Read_DE) / 60
##
## 0 1 2 3 4 5 6 7 8 10
## 1 4 2 13 5 4 2 7 3 2
collotask$Read_DE <- cut(collotask$Read_DE,
breaks = c(0, 3.5, 6.5, 10),
include.lowest = TRUE)
table(collotask$Read_DE) / 60
##
## [0,3.5] (3.5,6.5] (6.5,10]
## 20 11 12
# Self_rated_reading: alike
table(collotask$Self_rated_reading) / 60
##
## 5 6 7 8 9 10
## 1 3 3 13 12 11
collotask$Self_rated_reading <- cut(collotask$Self_rated_reading,
breaks=c(0, 8.5, 9.5, 10),
include.lowest = TRUE)
table(collotask$Self_rated_reading) / 60
##
## [0,8.5] (8.5,9.5] (9.5,10]
## 20 12 11
Full model without scaling:
model1l_noscale_nointeraction <- glmer(
Accuracy ~ (1|Item) + (1|Subject) +
Word_class + Majors + Script + Read_DE + Self_rated_reading + # factors
SPR + Onset + Age + Residence, # continuous
family = binomial, data = collotask)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0078977
## (tol = 0.001, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
Fails to converge! But: “[g]lmer fits may produce convergence warnings; these do not necessarily mean the fit is incorrect” (see ?convergence)
If you want to know how formulas work in lme4
:
https://cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf (p. 7) So we haven’t modelled this as a hierarchical model; the intercept can vary among Item
and Subject
.
summary(model1l_noscale_nointeraction)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Accuracy ~ (1 | Item) + (1 | Subject) + Word_class + Majors +
## Script + Read_DE + Self_rated_reading + SPR + Onset + Age +
## Residence
## Data: collotask
##
## AIC BIC logLik deviance df.resid
## 2849.5 2949.0 -1407.7 2815.5 2563
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2646 -0.6684 0.2866 0.6217 4.0903
##
## Random effects:
## Groups Name Variance Std.Dev.
## Item (Intercept) 1.3724 1.1715
## Subject (Intercept) 0.2274 0.4769
## Number of obs: 2580, groups: Item, 60; Subject, 43
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7499199 0.6965972 1.077 0.281683
## Word_classADV + ADJ 0.7295050 0.4511677 1.617 0.105894
## Word_classV + ADV -0.3942169 0.4497639 -0.876 0.380760
## Word_classV + N 0.1939114 0.4508383 0.430 0.667114
## MajorsOther -0.3114654 0.1981437 -1.572 0.115970
## ScriptOverlapping script 0.5308012 0.2455482 2.162 0.030641 *
## ScriptSame script 0.2180700 0.2378999 0.917 0.359328
## Read_DE(3.5,6.5] 0.3589830 0.2406344 1.492 0.135747
## Read_DE(6.5,10] 0.3537253 0.2305179 1.534 0.124911
## Self_rated_reading(8.5,9.5] 0.4058718 0.2515524 1.613 0.106643
## Self_rated_reading(9.5,10] 0.0985052 0.3000132 0.328 0.742657
## SPR -0.0005537 0.0021830 -0.254 0.799757
## Onset -0.0356359 0.0184366 -1.933 0.053250 .
## Age -0.0193631 0.0213453 -0.907 0.364335
## Residence 0.0069099 0.0019195 3.600 0.000318 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## convergence code: 0
## Model failed to converge with max|grad| = 0.0078977 (tol = 0.001, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
With scaling of continuous predictors (seems to be best-practice)
model1l_nointeraction <- glmer(
Accuracy ~ (1|Item) + (1|Subject) +
Word_class + Majors + Script + Read_DE + Self_rated_reading + # factors
scale(SPR) + scale(Onset) + scale(Age) + scale(Residence), # continuous
family = binomial, data = collotask)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0027427
## (tol = 0.001, component 1)
# fails to converge!
summary(model1l_nointeraction)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Accuracy ~ (1 | Item) + (1 | Subject) + Word_class + Majors +
## Script + Read_DE + Self_rated_reading + scale(SPR) + scale(Onset) +
## scale(Age) + scale(Residence)
## Data: collotask
##
## AIC BIC logLik deviance df.resid
## 2849.5 2949.0 -1407.7 2815.5 2563
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2640 -0.6684 0.2866 0.6217 4.0895
##
## Random effects:
## Groups Name Variance Std.Dev.
## Item (Intercept) 1.3719 1.1713
## Subject (Intercept) 0.2273 0.4768
## Number of obs: 2580, groups: Item, 60; Subject, 43
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.001149 0.388588 0.003 0.997641
## Word_classADV + ADJ 0.729241 0.451212 1.616 0.106055
## Word_classV + ADV -0.394883 0.449777 -0.878 0.379970
## Word_classV + N 0.192956 0.450865 0.428 0.668674
## MajorsOther -0.311668 0.198122 -1.573 0.115694
## ScriptOverlapping script 0.530262 0.245528 2.160 0.030797 *
## ScriptSame script 0.217644 0.237876 0.915 0.360218
## Read_DE(3.5,6.5] 0.358907 0.240613 1.492 0.135794
## Read_DE(6.5,10] 0.353898 0.230498 1.535 0.124695
## Self_rated_reading(8.5,9.5] 0.405970 0.251528 1.614 0.106524
## Self_rated_reading(9.5,10] 0.098499 0.299987 0.328 0.742653
## scale(SPR) -0.025880 0.101637 -0.255 0.799005
## scale(Onset) -0.236736 0.122518 -1.932 0.053328 .
## scale(Age) -0.145687 0.160438 -0.908 0.363847
## scale(Residence) 0.571229 0.158648 3.601 0.000317 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## convergence code: 0
## Model failed to converge with max|grad| = 0.0027427 (tol = 0.001, component 1)
Pseudo-R² (conditional and marginal, see Nagakawa/Schielzeth (2013) or at least ?r.squaredGLMM
):
“theoretical”: value computed for binomial family (=> logistic regression) “delta”: value computed for all distributions and link functions
r.squaredGLMM(model1l_nointeraction)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help
## page.
## Warning: The null model is correct only if all variables used by the
## original model remain unchanged.
## R2m R2c
## theoretical 0.1322254 0.4160710
## delta 0.1150818 0.3621255
car::vif(model1l_nointeraction)
## GVIF Df GVIF^(1/(2*Df))
## Word_class 1.000364 3 1.000061
## Majors 1.281381 1 1.131981
## Script 1.910760 2 1.175714
## Read_DE 1.540247 2 1.114032
## Self_rated_reading 2.334398 2 1.236072
## scale(SPR) 1.362885 1 1.167427
## scale(Onset) 1.966491 1 1.402316
## scale(Age) 3.125145 1 1.767808
## scale(Residence) 2.929310 1 1.711523
“use allFit to try the fit with all available optimizers (e.g. several different implementations of BOBYQA and Nelder-Mead, L-BFGS-B from optim, nlminb, …). While this will of course be slow for large fits, we consider it the gold standard; if all optimizers converge to values that are practically equivalent, then we would consider the convergence warnings to be false positives.” (?convergence
)
model1l_nointeraction_all <- allFit(model1l_nointeraction, parallel = "multicore")
## Loading required namespace: dfoptim
## bobyqa : [OK]
## Nelder_Mead :
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00349833
## (tol = 0.001, component 1)
## [OK]
## nlminbwrap : [OK]
## nmkbw :
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00197771
## (tol = 0.001, component 1)
## [OK]
## optimx.L-BFGS-B : [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD :
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0186393
## (tol = 0.001, component 1)
## [OK]
## nloptwrap.NLOPT_LN_BOBYQA :
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0186393
## (tol = 0.001, component 1)
## [OK]
summary(model1l_nointeraction_all)
## $which.OK
## bobyqa Nelder_Mead
## TRUE TRUE
## nlminbwrap nmkbw
## TRUE TRUE
## optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD
## FALSE TRUE
## nloptwrap.NLOPT_LN_BOBYQA
## TRUE
##
## $msgs
## $msgs$bobyqa
## NULL
##
## $msgs$Nelder_Mead
## [1] "Model failed to converge with max|grad| = 0.00349833 (tol = 0.001, component 1)"
##
## $msgs$nlminbwrap
## NULL
##
## $msgs$nmkbw
## [1] "Model failed to converge with max|grad| = 0.00197771 (tol = 0.001, component 1)"
##
## $msgs$nloptwrap.NLOPT_LN_NELDERMEAD
## [1] "Model failed to converge with max|grad| = 0.0186393 (tol = 0.001, component 1)"
##
## $msgs$nloptwrap.NLOPT_LN_BOBYQA
## [1] "Model failed to converge with max|grad| = 0.0186393 (tol = 0.001, component 1)"
##
##
## $fixef
## (Intercept) Word_classADV + ADJ
## bobyqa 0.001002746 0.7290180
## Nelder_Mead 0.001162473 0.7294117
## nlminbwrap 0.001092633 0.7289064
## nmkbw 0.001120910 0.7287662
## nloptwrap.NLOPT_LN_NELDERMEAD 0.002048260 0.7294751
## nloptwrap.NLOPT_LN_BOBYQA 0.002048260 0.7294751
## Word_classV + ADV Word_classV + N
## bobyqa -0.3949318 0.1932674
## Nelder_Mead -0.3944329 0.1936554
## nlminbwrap -0.3950053 0.1931752
## nmkbw -0.3951023 0.1929884
## nloptwrap.NLOPT_LN_NELDERMEAD -0.3921874 0.1963306
## nloptwrap.NLOPT_LN_BOBYQA -0.3921874 0.1963306
## MajorsOther ScriptOverlapping script
## bobyqa -0.3118840 0.5302395
## Nelder_Mead -0.3119496 0.5297724
## nlminbwrap -0.3118899 0.5302308
## nmkbw -0.3118123 0.5305088
## nloptwrap.NLOPT_LN_NELDERMEAD -0.3107849 0.5302793
## nloptwrap.NLOPT_LN_BOBYQA -0.3107849 0.5302793
## ScriptSame script Read_DE(3.5,6.5]
## bobyqa 0.2177443 0.3587764
## Nelder_Mead 0.2174216 0.3583727
## nlminbwrap 0.2177353 0.3587573
## nmkbw 0.2177846 0.3588464
## nloptwrap.NLOPT_LN_NELDERMEAD 0.2191564 0.3577619
## nloptwrap.NLOPT_LN_BOBYQA 0.2191564 0.3577619
## Read_DE(6.5,10] Self_rated_reading(8.5,9.5]
## bobyqa 0.3539272 0.4063498
## Nelder_Mead 0.3536658 0.4059844
## nlminbwrap 0.3539121 0.4063444
## nmkbw 0.3539536 0.4062045
## nloptwrap.NLOPT_LN_NELDERMEAD 0.3523673 0.4031273
## nloptwrap.NLOPT_LN_BOBYQA 0.3523673 0.4031273
## Self_rated_reading(9.5,10] scale(SPR)
## bobyqa 0.09903520 -0.02591064
## Nelder_Mead 0.09891321 -0.02591543
## nlminbwrap 0.09904073 -0.02591149
## nmkbw 0.09875128 -0.02592940
## nloptwrap.NLOPT_LN_NELDERMEAD 0.09448333 -0.02428635
## nloptwrap.NLOPT_LN_BOBYQA 0.09448333 -0.02428635
## scale(Onset) scale(Age) scale(Residence)
## bobyqa -0.2365218 -0.1458172 0.5712211
## Nelder_Mead -0.2367057 -0.1459578 0.5714390
## nlminbwrap -0.2365241 -0.1458140 0.5712165
## nmkbw -0.2365430 -0.1457750 0.5712934
## nloptwrap.NLOPT_LN_NELDERMEAD -0.2379770 -0.1463040 0.5724397
## nloptwrap.NLOPT_LN_BOBYQA -0.2379770 -0.1463040 0.5724397
##
## $llik
## bobyqa Nelder_Mead
## -1407.744 -1407.744
## nlminbwrap nmkbw
## -1407.744 -1407.744
## nloptwrap.NLOPT_LN_NELDERMEAD nloptwrap.NLOPT_LN_BOBYQA
## -1407.744 -1407.744
##
## $sdcor
## Item.(Intercept) Subject.(Intercept)
## bobyqa 1.171202 0.4768522
## Nelder_Mead 1.171058 0.4767834
## nlminbwrap 1.171203 0.4768511
## nmkbw 1.171198 0.4768775
## nloptwrap.NLOPT_LN_NELDERMEAD 1.170828 0.4770673
## nloptwrap.NLOPT_LN_BOBYQA 1.170828 0.4770673
##
## $theta
## Item.(Intercept) Subject.(Intercept)
## bobyqa 1.171202 0.4768522
## Nelder_Mead 1.171058 0.4767834
## nlminbwrap 1.171203 0.4768511
## nmkbw 1.171198 0.4768775
## nloptwrap.NLOPT_LN_NELDERMEAD 1.170828 0.4770673
## nloptwrap.NLOPT_LN_BOBYQA 1.170828 0.4770673
##
## $times
## user.self sys.self elapsed user.child
## bobyqa 7.25 0.03 7.28 NA
## Nelder_Mead 11.23 0.00 11.24 NA
## nlminbwrap 9.87 0.00 9.88 NA
## nmkbw 9.99 0.00 9.98 NA
## nloptwrap.NLOPT_LN_NELDERMEAD 4.17 0.00 4.17 NA
## nloptwrap.NLOPT_LN_BOBYQA 4.12 0.00 4.15 NA
## sys.child
## bobyqa NA
## Nelder_Mead NA
## nlminbwrap NA
## nmkbw NA
## nloptwrap.NLOPT_LN_NELDERMEAD NA
## nloptwrap.NLOPT_LN_BOBYQA NA
##
## $feval
## bobyqa Nelder_Mead
## 859 1617
## nlminbwrap nmkbw
## NA 1304
## nloptwrap.NLOPT_LN_NELDERMEAD nloptwrap.NLOPT_LN_BOBYQA
## 230 230
##
## attr(,"class")
## [1] "summary.allFit"
=> practically indistinguishable => convergence warnings shouldn’t bother us too much
The influence of Age
seems to be insignificant, so let’s drop it:
model2l_nointeraction <- glmer(
Accuracy ~ (1|Item) + (1|Subject) +
Word_class + Majors + Script + Read_DE + Self_rated_reading + # factors
scale(SPR) + scale(Onset) + scale(Residence), # continuous
family = binomial, data = collotask)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00421895
## (tol = 0.001, component 1)
# fails to converge!
summary(model2l_nointeraction)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Accuracy ~ (1 | Item) + (1 | Subject) + Word_class + Majors +
## Script + Read_DE + Self_rated_reading + scale(SPR) + scale(Onset) +
## scale(Residence)
## Data: collotask
##
## AIC BIC logLik deviance df.resid
## 2848.3 2942.0 -1408.2 2816.3 2564
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2279 -0.6687 0.2866 0.6247 4.1154
##
## Random effects:
## Groups Name Variance Std.Dev.
## Item (Intercept) 1.3716 1.1712
## Subject (Intercept) 0.2347 0.4845
## Number of obs: 2580, groups: Item, 60; Subject, 43
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.009524 0.390063 0.024 0.9805
## Word_classADV + ADJ 0.729052 0.451229 1.616 0.1062
## Word_classV + ADV -0.395086 0.449824 -0.878 0.3798
## Word_classV + N 0.193245 0.450877 0.429 0.6682
## MajorsOther -0.298054 0.199818 -1.492 0.1358
## ScriptOverlapping script 0.564885 0.245434 2.302 0.0214 *
## ScriptSame script 0.189598 0.238509 0.795 0.4267
## Read_DE(3.5,6.5] 0.359637 0.243362 1.478 0.1395
## Read_DE(6.5,10] 0.354549 0.233148 1.521 0.1283
## Self_rated_reading(8.5,9.5] 0.396800 0.254175 1.561 0.1185
## Self_rated_reading(9.5,10] 0.048382 0.298518 0.162 0.8712
## scale(SPR) -0.026347 0.102814 -0.256 0.7978
## scale(Onset) -0.269547 0.118578 -2.273 0.0230 *
## scale(Residence) 0.470489 0.114118 4.123 3.74e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## convergence code: 0
## Model failed to converge with max|grad| = 0.00421895 (tol = 0.001, component 1)
r.squaredGLMM(model2l_nointeraction)
## Warning: The null model is correct only if all variables used by the
## original model remain unchanged.
## R2m R2c
## theoretical 0.1314365 0.4163983
## delta 0.1144035 0.3624367
Let’s also try to exclude Majors
instead:
model3l_nointeraction <- glmer(
Accuracy ~ Word_class + Script + Read_DE + Self_rated_reading +
scale(Onset) + scale(SPR) + scale(Residence) + scale(Age) +
(1|Item) + (1|Subject),
family = binomial,
data = collotask
)
# => Majors leads to non-convergence
summary(model3l_nointeraction)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Accuracy ~ Word_class + Script + Read_DE + Self_rated_reading +
## scale(Onset) + scale(SPR) + scale(Residence) + scale(Age) +
## (1 | Item) + (1 | Subject)
## Data: collotask
##
## AIC BIC logLik deviance df.resid
## 2849.9 2943.6 -1408.9 2817.9 2564
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2334 -0.6659 0.2840 0.6200 4.1989
##
## Random effects:
## Groups Name Variance Std.Dev.
## Item (Intercept) 1.3718 1.1712
## Subject (Intercept) 0.2496 0.4996
## Number of obs: 2580, groups: Item, 60; Subject, 43
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.15448 0.37942 -0.407 0.683900
## Word_classADV + ADJ 0.72953 0.45110 1.617 0.105831
## Word_classV + ADV -0.39426 0.44967 -0.877 0.380608
## Word_classV + N 0.19363 0.45075 0.430 0.667503
## ScriptOverlapping script 0.63084 0.24532 2.572 0.010124 *
## ScriptSame script 0.25125 0.24500 1.026 0.305123
## Read_DE(3.5,6.5] 0.41223 0.24628 1.674 0.094157 .
## Read_DE(6.5,10] 0.27303 0.23206 1.177 0.239375
## Self_rated_reading(8.5,9.5] 0.33950 0.25614 1.325 0.185014
## Self_rated_reading(9.5,10] 0.01024 0.30494 0.034 0.973215
## scale(Onset) -0.28567 0.12275 -2.327 0.019948 *
## scale(SPR) -0.03320 0.10497 -0.316 0.751779
## scale(Residence) 0.57953 0.16394 3.535 0.000408 ***
## scale(Age) -0.12645 0.16542 -0.764 0.444614
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
print(model3l_nointeraction, correlation = TRUE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Accuracy ~ Word_class + Script + Read_DE + Self_rated_reading +
## scale(Onset) + scale(SPR) + scale(Residence) + scale(Age) +
## (1 | Item) + (1 | Subject)
## Data: collotask
## AIC BIC logLik deviance df.resid
## 2849.868 2943.557 -1408.934 2817.868 2564
## Random effects:
## Groups Name Std.Dev.
## Item (Intercept) 1.1712
## Subject (Intercept) 0.4996
## Number of obs: 2580, groups: Item, 60; Subject, 43
## Fixed Effects:
## (Intercept) Word_classADV + ADJ
## -0.15448 0.72953
## Word_classV + ADV Word_classV + N
## -0.39426 0.19363
## ScriptOverlapping script ScriptSame script
## 0.63084 0.25125
## Read_DE(3.5,6.5] Read_DE(6.5,10]
## 0.41223 0.27303
## Self_rated_reading(8.5,9.5] Self_rated_reading(9.5,10]
## 0.33950 0.01024
## scale(Onset) scale(SPR)
## -0.28567 -0.03320
## scale(Residence) scale(Age)
## 0.57953 -0.12645
r.squaredGLMM(model3l_nointeraction)
## Warning: The null model is correct only if all variables used by the
## original model remain unchanged.
## R2m R2c
## theoretical 0.1295903 0.4169367
## delta 0.1128100 0.3629487
car::vif(model3l_nointeraction) # no apparent problems with multicollinearity
## GVIF Df GVIF^(1/(2*Df))
## Word_class 1.000328 3 1.000055
## Script 1.785485 2 1.155950
## Read_DE 1.375312 2 1.082930
## Self_rated_reading 2.239665 2 1.223336
## scale(Onset) 1.848544 1 1.359612
## scale(SPR) 1.360442 1 1.166380
## scale(Residence) 2.932942 1 1.712583
## scale(Age) 3.112630 1 1.764265
Let’s see if we can still drop Age
:
model4l_nointeraction <- glmer(
Accuracy ~ Word_class + Script + Read_DE + Self_rated_reading +
scale(Onset) + scale(SPR) + scale(Residence) +
(1|Subject) + (1|Item),
family = binomial,
data = collotask
)
summary(model4l_nointeraction)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Accuracy ~ Word_class + Script + Read_DE + Self_rated_reading +
## scale(Onset) + scale(SPR) + scale(Residence) + (1 | Subject) +
## (1 | Item)
## Data: collotask
##
## AIC BIC logLik deviance df.resid
## 2848.4 2936.3 -1409.2 2818.4 2565
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2056 -0.6671 0.2860 0.6206 4.2149
##
## Random effects:
## Groups Name Variance Std.Dev.
## Item (Intercept) 1.372 1.1712
## Subject (Intercept) 0.255 0.5049
## Number of obs: 2580, groups: Item, 60; Subject, 43
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.14114 0.38010 -0.371 0.71039
## Word_classADV + ADJ 0.72938 0.45125 1.616 0.10602
## Word_classV + ADV -0.39436 0.44980 -0.877 0.38063
## Word_classV + N 0.19375 0.45087 0.430 0.66739
## ScriptOverlapping script 0.65715 0.24485 2.684 0.00728 **
## ScriptSame script 0.22560 0.24452 0.923 0.35620
## Read_DE(3.5,6.5] 0.41096 0.24820 1.656 0.09777 .
## Read_DE(6.5,10] 0.27670 0.23387 1.183 0.23675
## Self_rated_reading(8.5,9.5] 0.33376 0.25802 1.294 0.19582
## Self_rated_reading(9.5,10] -0.03007 0.30288 -0.099 0.92091
## scale(Onset) -0.31244 0.11870 -2.632 0.00848 **
## scale(SPR) -0.03344 0.10581 -0.316 0.75194
## scale(Residence) 0.49128 0.11685 4.204 2.62e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
r.squaredGLMM(model4l_nointeraction)
## Warning: The null model is correct only if all variables used by the
## original model remain unchanged.
## R2m R2c
## theoretical 0.1290290 0.4172089
## delta 0.1123282 0.3632077
car::vif(model4l_nointeraction)
## GVIF Df GVIF^(1/(2*Df))
## Word_class 1.000319 3 1.000053
## Script 1.685995 2 1.139499
## Read_DE 1.375057 2 1.082880
## Self_rated_reading 2.165983 2 1.213148
## scale(Onset) 1.702813 1 1.304919
## scale(SPR) 1.360405 1 1.166364
## scale(Residence) 1.466393 1 1.210947
=> Slightly lower AIC, so it’s definitely not a worse model.
Drop SPR
:
model5l_nointeraction <- glmer(
Accuracy ~ Word_class + Script + Read_DE + Self_rated_reading +
scale(Onset) + scale(Residence) +
(1|Item) + (1|Subject),
family = binomial,
data = collotask
)
summary(model5l_nointeraction)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Accuracy ~ Word_class + Script + Read_DE + Self_rated_reading +
## scale(Onset) + scale(Residence) + (1 | Item) + (1 | Subject)
## Data: collotask
##
## AIC BIC logLik deviance df.resid
## 2846.5 2928.5 -1409.3 2818.5 2566
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2032 -0.6706 0.2864 0.6201 4.2097
##
## Random effects:
## Groups Name Variance Std.Dev.
## Item (Intercept) 1.3719 1.1713
## Subject (Intercept) 0.2561 0.5061
## Number of obs: 2580, groups: Item, 60; Subject, 43
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.14867 0.37955 -0.392 0.69528
## Word_classADV + ADJ 0.72961 0.45123 1.617 0.10589
## Word_classV + ADV -0.39399 0.44980 -0.876 0.38107
## Word_classV + N 0.19409 0.45087 0.430 0.66684
## ScriptOverlapping script 0.65503 0.24516 2.672 0.00754 **
## ScriptSame script 0.24527 0.23680 1.036 0.30031
## Read_DE(3.5,6.5] 0.41998 0.24694 1.701 0.08900 .
## Read_DE(6.5,10] 0.28113 0.23387 1.202 0.22933
## Self_rated_reading(8.5,9.5] 0.33267 0.25840 1.287 0.19795
## Self_rated_reading(9.5,10] -0.03843 0.30228 -0.127 0.89884
## scale(Onset) -0.32764 0.10880 -3.011 0.00260 **
## scale(Residence) 0.48956 0.11686 4.189 2.8e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) W_AD+A W_V+AD Wr_V+N ScrpOs ScrpSs R_DE(3 R_DE(6 S__(8.
## Wrd_ADV+ADJ -0.599
## Wrd_clV+ADV -0.598 0.503
## Wrd_clssV+N -0.598 0.503 0.504
## ScrptOvrlps -0.218 0.004 -0.002 0.001
## ScrptSmscrp -0.270 0.002 0.000 0.001 0.354
## R_DE(3.5,6. -0.182 0.002 -0.002 0.000 0.216 -0.132
## R_DE(6.5,10 -0.202 0.002 -0.001 0.001 0.112 -0.045 0.401
## S__(8.5,9.5 -0.231 0.002 -0.001 0.000 -0.226 -0.090 -0.046 0.035
## S__(9.5,10] -0.208 0.000 -0.001 0.000 -0.245 0.180 -0.270 -0.259 0.509
## scale(Onst) -0.195 -0.004 0.002 -0.001 0.054 0.052 -0.040 0.068 0.407
## scal(Rsdnc) 0.187 0.007 -0.001 0.003 0.085 -0.307 -0.033 0.128 -0.293
## S__(9. scl(O)
## Wrd_ADV+ADJ
## Wrd_clV+ADV
## Wrd_clssV+N
## ScrptOvrlps
## ScrptSmscrp
## R_DE(3.5,6.
## R_DE(6.5,10
## S__(8.5,9.5
## S__(9.5,10]
## scale(Onst) 0.436
## scal(Rsdnc) -0.435 -0.215
r.squaredGLMM(model5l_nointeraction)
## Warning: The null model is correct only if all variables used by the
## original model remain unchanged.
## R2m R2c
## theoretical 0.1288586 0.4172371
## delta 0.1121806 0.3632345
car::vif(model5l_nointeraction)
## GVIF Df GVIF^(1/(2*Df))
## Word_class 1.000317 3 1.000053
## Script 1.552560 2 1.116252
## Read_DE 1.356791 2 1.079265
## Self_rated_reading 2.148059 2 1.210630
## scale(Onset) 1.425920 1 1.194119
## scale(Residence) 1.463135 1 1.209601
=> AIC almost unchanged.
Drop Self_rated_reading
:
model6l_nointeraction <- glmer(
Accuracy ~ Word_class + Script + Read_DE +
scale(Onset) + scale(Residence) +
(1|Item) + (1|Subject),
family = binomial,
data = collotask
)
summary(model6l_nointeraction)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Accuracy ~ Word_class + Script + Read_DE + scale(Onset) + scale(Residence) +
## (1 | Item) + (1 | Subject)
## Data: collotask
##
## AIC BIC logLik deviance df.resid
## 2845.0 2915.2 -1410.5 2821.0 2568
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2817 -0.6693 0.2859 0.6194 4.1770
##
## Random effects:
## Groups Name Variance Std.Dev.
## Item (Intercept) 1.3720 1.1713
## Subject (Intercept) 0.2763 0.5256
## Number of obs: 2580, groups: Item, 60; Subject, 43
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.07132 0.36973 -0.193 0.847044
## Word_classADV + ADJ 0.72932 0.45108 1.617 0.105917
## Word_classV + ADV -0.39434 0.44971 -0.877 0.380558
## Word_classV + N 0.19374 0.45075 0.430 0.667327
## ScriptOverlapping script 0.69348 0.24254 2.859 0.004246 **
## ScriptSame script 0.33134 0.23392 1.416 0.156637
## Read_DE(3.5,6.5] 0.36912 0.24284 1.520 0.128513
## Read_DE(6.5,10] 0.20249 0.22738 0.891 0.373168
## scale(Onset) -0.35945 0.09768 -3.680 0.000234 ***
## scale(Residence) 0.49951 0.10744 4.649 3.33e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) W_AD+A W_V+AD Wr_V+N ScrpOs ScrpSs R_DE(3 R_DE(6 scl(O)
## Wrd_ADV+ADJ -0.614
## Wrd_clV+ADV -0.614 0.503
## Wrd_clssV+N -0.613 0.502 0.504
## ScrptOvrlps -0.313 0.004 -0.002 0.001
## ScrptSmscrp -0.290 0.003 -0.001 0.001 0.404
## R_DE(3.5,6. -0.245 0.002 -0.002 0.000 0.176 -0.066
## R_DE(6.5,10 -0.253 0.001 -0.001 0.000 0.077 0.047 0.342
## scale(Onst) -0.088 -0.005 0.003 -0.001 0.221 0.022 0.064 0.168
## scal(Rsdnc) 0.099 0.008 -0.001 0.003 -0.037 -0.287 -0.164 0.037 -0.008
r.squaredGLMM(model6l_nointeraction)
## Warning: The null model is correct only if all variables used by the
## original model remain unchanged.
## R2m R2c
## theoretical 0.1252057 0.4171957
## delta 0.1089994 0.3631951
car::vif(model6l_nointeraction)
## GVIF Df GVIF^(1/(2*Df))
## Word_class 1.000287 3 1.000048
## Script 1.270650 2 1.061711
## Read_DE 1.168096 2 1.039608
## scale(Onset) 1.087833 1 1.042992
## scale(Residence) 1.176450 1 1.084643
Drop Read_DE
:
model7l_nointeraction <- glmer(
Accuracy ~ Word_class + Script +
scale(Onset) + scale(Residence) +
(1|Item) + (1|Subject),
family = binomial,
data = collotask
)
summary(model7l_nointeraction)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Accuracy ~ Word_class + Script + scale(Onset) + scale(Residence) +
## (1 | Item) + (1 | Subject)
## Data: collotask
##
## AIC BIC logLik deviance df.resid
## 2843.3 2901.9 -1411.7 2823.3 2570
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2038 -0.6679 0.2874 0.6193 4.1188
##
## Random effects:
## Groups Name Variance Std.Dev.
## Item (Intercept) 1.3720 1.171
## Subject (Intercept) 0.3025 0.550
## Number of obs: 2580, groups: Item, 60; Subject, 43
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.09296 0.35465 0.262 0.79322
## Word_classADV + ADJ 0.72970 0.45112 1.618 0.10576
## Word_classV + ADV -0.39368 0.44974 -0.875 0.38138
## Word_classV + N 0.19426 0.45079 0.431 0.66651
## ScriptOverlapping script 0.62748 0.24698 2.541 0.01106 *
## ScriptSame script 0.34969 0.24082 1.452 0.14647
## scale(Onset) -0.37613 0.09962 -3.776 0.00016 ***
## scale(Residence) 0.52348 0.10899 4.803 1.56e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) W_AD+A W_V+AD Wr_V+N ScrpOs ScrpSs scl(O)
## Wrd_ADV+ADJ -0.639
## Wrd_clV+ADV -0.640 0.503
## Wrd_clssV+N -0.639 0.503 0.504
## ScrptOvrlps -0.293 0.004 -0.002 0.001
## ScrptSmscrp -0.317 0.003 -0.001 0.001 0.423
## scale(Onst) -0.049 -0.005 0.003 -0.001 0.214 0.015
## scal(Rsdnc) 0.085 0.008 -0.002 0.003 -0.011 -0.313 -0.015
r.squaredGLMM(model7l_nointeraction)
## Warning: The null model is correct only if all variables used by the
## original model remain unchanged.
## R2m R2c
## theoretical 0.1223061 0.4183543
## delta 0.1065025 0.3642974
car::vif(model7l_nointeraction)
## GVIF Df GVIF^(1/(2*Df))
## Word_class 1.000253 3 1.000042
## Script 1.195892 2 1.045738
## scale(Onset) 1.057494 1 1.028345
## scale(Residence) 1.133449 1 1.064636
We could still drop Word_class
…
So far, we have also included random intercepts. The model gets somewhat more complicated (but possibly better) if we also include random slopes.
“[R]esearchers in ecology (Schielzeth & Forstmeier, 2009), psycholinguistics (Barr, Levy, Scheepers, & Tilly, 2013) and other fields have shown via simulations that mixed models without random slopes are anti-conservative or, in other words, they have a relatively high Type I error rate (they tend to find a lot of significant results which are actually due to chance).” (Bodo Winter: http://www.bodowinter.com/tutorial/bw_LME_tutorial2.pdf)
Still to do here: Checking model assumptions.