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.