# readmode_pred

#### 2024-05-11

version of R:4.3.2 version of Rstudio: 2023.12.0+369

We run the models for each dependent variables (on FFD,SFD,GD,TRT,RPD) in the following order:

Full_mdoel=lmer(depvar~readingmode * predictability+(1+readingmode * predictability|pp)+(1+readingmode*predictability|item), datafile)

model=lmer(depvar~readingmode * predictability+(1+readingmode*predictability|pp)+(1|item), datafile)

model=lmer(depvar~readingmode* predictability+(0+readingmode*predictability|pp)+(1|item), datafile)

model=lmer(depvar~readingmode*predictability+(1+readingmode+predictability|pp)+(1|item), datafile)

model=lmer(depvar~readingmode*predictability+(0+readingmode+predictability|pp)+(1|item), datafile)

model=lmer(depvar ~readingmode*predictability+ (1+readingmode|pp) + (1|item), datafile)

model=lmer(depvar ~readingmode*predictability+ (0+readingmode|pp) + (1|item), datafile)

model=lmer(depvar ~readingmode*predictability+ (1+predictability|pp) + (1|item), datafile)

model=lmer(depvar ~readingmode*predictability+ (0+predictability|pp) + (1|item), datafile)

model=lmer(depvar ~readingmode*predictability+ (1|pp) + (1|item), datafile)

Only the model that successfully converged for each dependent variables was retained.

    library(lmerTest) # version 3.1-3

    ## 载入需要的程序包：lme4

    ## 载入需要的程序包：Matrix

    ##
    ## 载入程序包：'lmerTest'

    ## The following object is masked from 'package:lme4':
    ##
    ##     lmer

    ## The following object is masked from 'package:stats':
    ##
    ##     step

    library(lattice)  # version 0.21-9
    library(MASS)     # version 7.3-60
    library(plyr)     # version 1.8.9
    library(plotrix)  # version 3.8-4

    rm(list=ls())
    datafile0 = read.csv("datanew.csv", sep = ",", dec = ".")
    datafile <- datafile0
    colnames(datafile)

    ##  [1] "pp"             "item"           "readingmode"    "predictability"
    ##  [5] "condition"      "TRT"            "IA_SKIP"        "FFD"
    ##  [9] "SFD"            "GD"             "RPD"            "sacinlength"
    ## [13] "accuracy"       "ifquestion"

    nrow(datafile)

    ## [1] 5472

    datafile$pp = datafile[,1]
    datafile$item = datafile[,2]
    datafile$readingmode = datafile[,3]
    datafile$predictability = datafile[,4]
    datafile$condition = datafile[,5]
    datafile$pp = as.factor(datafile$pp)
    datafile$item = as.factor(datafile$item)
    datafile$readingmode = factor(datafile$readingmode, levels = c("S","A"))
    datafile$predictability = factor(datafile$predictability, levels = c("H","L"))
    datafile$condition = factor(datafile$condition, levels = c("SH","SL","AH","AL"))

# FOR ACCURACY

    measure = "accuracy"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)

## descriptive statistic analyses for accuracy

    mean.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), mean, na.rm = T)
    sd.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), sd, na.rm = T)
    se.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), std.error, na.rm = T)
    grand.mean=apply(mean.matrix, 2, mean, na.rm = T)
    grand.sd=apply(sd.matrix,2,mean, na.rm = T)
    grand.se=apply(se.matrix,2, mean, na.rm = T)
    data.frame <- data.frame(grand.mean,grand.sd,grand.se)
    data.frame

    ##    grand.mean   grand.sd   grand.se
    ## SH  0.9653439 0.07014272 0.03465608
    ## SL  0.9476562 0.07259828 0.03730159
    ## AH  0.9562500 0.09045085 0.04375000
    ## AL  0.9356771 0.11707317 0.05771901

## regression analysis for accuracy

    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    model_accuracy = glm(depvar~readingmode*predictability,datafile,family=binomial)
    summary(model_accuracy)

    ##
    ## Call:
    ## glm(formula = depvar ~ readingmode * predictability, family = binomial,
    ##     data = datafile)
    ##
    ## Coefficients:
    ##                                  Estimate Std. Error z value Pr(>|z|)
    ## (Intercept)                        3.0124     0.1521  19.805   <2e-16 ***
    ## readingmode2-1                    -0.3903     0.3042  -1.283    0.200
    ## predictability2-1                 -0.2500     0.3042  -0.822    0.411
    ## readingmode2-1:predictability2-1  -0.2636     0.6084  -0.433    0.665
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## (Dispersion parameter for binomial family taken to be 1)
    ##
    ##     Null deviance: 384.67  on 994  degrees of freedom
    ## Residual deviance: 381.81  on 991  degrees of freedom
    ##   (因为不存在，4477个观察量被删除了)
    ## AIC: 389.81
    ##
    ## Number of Fisher Scoring iterations: 6

# FOR IA_SKIP RATE

    measure = "IA_SKIP"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)

## descriptive statistic analyses for skipping rate

    mean.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), mean, na.rm = T)
    sd.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), sd, na.rm = T)
    se.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), std.error, na.rm = T)
    grand.mean=apply(mean.matrix, 2, mean, na.rm = T)
    grand.sd=apply(sd.matrix,2,mean, na.rm = T)
    grand.se=apply(se.matrix,2, mean, na.rm = T)
    data.frame <- data.frame(grand.mean,grand.sd,grand.se)
    data.frame

    ##    grand.mean  grand.sd   grand.se
    ## SH  0.4356763 0.4575960 0.10112364
    ## SL  0.4063026 0.4505315 0.09901976
    ## AH  0.1193013 0.2702324 0.05776748
    ## AL  0.0988704 0.2368486 0.05088666

## linear mixed effect model analysis for skipping rate

    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    model_SKIP = glmer(depvar~readingmode*predictability+(1|pp)+(1|item),datafile,family=binomial)
    summary(model_SKIP)

    ## Generalized linear mixed model fit by maximum likelihood (Laplace
    ##   Approximation) [glmerMod]
    ##  Family: binomial  ( logit )
    ## Formula: depvar ~ readingmode * predictability + (1 | pp) + (1 | item)
    ##    Data: datafile
    ##
    ##      AIC      BIC   logLik deviance df.resid
    ##   5294.4   5334.1  -2641.2   5282.4     5466
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -2.1765 -0.5760 -0.3271  0.5215  7.2943
    ##
    ## Random effects:
    ##  Groups Name        Variance Std.Dev.
    ##  item   (Intercept) 0.0199   0.1411
    ##  pp     (Intercept) 0.5072   0.7122
    ## Number of obs: 5472, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                                  Estimate Std. Error z value Pr(>|z|)
    ## (Intercept)                      -1.34010    0.09855 -13.598   <2e-16 ***
    ## readingmode2-1                   -1.91795    0.07736 -24.792   <2e-16 ***
    ## predictability2-1                -0.18156    0.07485  -2.426   0.0153 *
    ## readingmode2-1:predictability2-1 -0.08871    0.14966  -0.593   0.5533
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Correlation of Fixed Effects:
    ##             (Intr) rdn2-1 prd2-1
    ## readngmd2-1 0.170
    ## prdctblt2-1 0.025  0.060
    ## rdng2-1:2-1 0.021  0.057  0.390

    confint(model_SKIP)

    ## Computing profile confidence intervals ...

    ##                                       2.5 %      97.5 %
    ## .sig01                            0.0000000  0.25702623
    ## .sig02                            0.5801780  0.88409561
    ## (Intercept)                      -1.5386555 -1.14650890
    ## readingmode2-1                   -2.0719422 -1.76747722
    ## predictability2-1                -0.3291933 -0.03451482
    ## readingmode2-1:predictability2-1 -0.3837663  0.20535466

# FOR FIRST FIXATION DURATION (FFD)

    measure = "FFD"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)

## descriptive statistic analyses for FFD

    mean.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), mean, na.rm = T)
    sd.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), sd, na.rm = T)
    se.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), std.error, na.rm = T)
    grand.mean=apply(mean.matrix, 2, mean, na.rm = T)
    grand.sd=apply(sd.matrix,2,mean, na.rm = T)
    grand.se=apply(se.matrix,2, mean, na.rm = T)
    data.frame <- data.frame(grand.mean,grand.sd,grand.se)
    data.frame

    ##    grand.mean  grand.sd grand.se
    ## SH   219.3804  62.99576 18.94725
    ## SL   219.0257  62.84565 17.67767
    ## AH   277.1583  99.07485 22.94028
    ## AL   292.1471 109.30093 24.94334

## linear mixed effect model analysis for FFD

    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    model_FFD= lmer(depvar ~readingmode*predictability+ (1+readingmode|pp) + (1|item), datafile)
    summary(model_FFD, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ readingmode * predictability + (1 + readingmode | pp) +
    ##     (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 48782.4
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -3.5323 -0.6004 -0.1516  0.4210  7.7360
    ##
    ## Random effects:
    ##  Groups   Name           Variance Std.Dev. Corr
    ##  item     (Intercept)      50.51   7.107
    ##  pp       (Intercept)    1019.17  31.924
    ##           readingmode2-1  704.61  26.545   0.25
    ##  Residual                8991.17  94.822
    ## Number of obs: 4070, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                                  Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)                       253.486      4.359   65.700  58.153   <2e-16
    ## readingmode2-1                     62.402      4.577   57.509  13.633   <2e-16
    ## predictability2-1                   7.592      3.061 3920.526   2.480   0.0132
    ## readingmode2-1:predictability2-1   16.796      6.122 3921.832   2.744   0.0061
    ##
    ## (Intercept)                      ***
    ## readingmode2-1                   ***
    ## predictability2-1                *
    ## readingmode2-1:predictability2-1 **
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    confint(model_FFD)

    ## Computing profile confidence intervals ...

    ##                                       2.5 %      97.5 %
    ## .sig01                             0.000000  11.8569401
    ## .sig02                            26.037508  38.9618228
    ## .sig03                            -0.109267   0.5766166
    ## .sig04                            17.590495  35.9046646
    ## .sigma                            92.722595  96.9538036
    ## (Intercept)                      244.888169 262.1126165
    ## readingmode2-1                    53.415809  71.5690108
    ## predictability2-1                  1.592379  13.5923197
    ## readingmode2-1:predictability2-1   4.802486  28.8023877

## simple effect analysis for FFD

    contrast.matrix <- matrix(c(
      -1, 1, 0, 0,
      0, 0, -1, 1), 4, 2,
      dimnames=list(c("SH","SL","AH","AL"),
                    c("SL-SH", "AL-AH")))
    inv.contrast.matrix <- matrix(t(ginv(contrast.matrix)), 4, 2,
                                  dimnames=list(c("SH","SL","AH","AL"),
                                                c("SL-SH", "AL-AH")))
    (contrasts(datafile$condition) <- fractions(inv.contrast.matrix))

    ##    SL-SH AL-AH
    ## SH -1/2     0
    ## SL  1/2     0
    ## AH    0  -1/2
    ## AL    0   1/2

    depvar.lmer1 = lmer(depvar ~ condition + (1|pp)+(1|item), datafile)
    summary(depvar.lmer1, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ condition + (1 | pp) + (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 48807.5
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -3.3827 -0.6168 -0.1696  0.4374  7.7367
    ##
    ## Random effects:
    ##  Groups   Name        Variance Std.Dev.
    ##  item     (Intercept)   48.71   6.979
    ##  pp       (Intercept) 1103.05  33.212
    ##  Residual             9147.39  95.642
    ## Number of obs: 4070, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)     254.157      4.501   65.814  56.466  < 2e-16 ***
    ## conditionSL-SH   -1.304      4.824 3969.537  -0.270    0.787
    ## conditionAL-AH   16.184      3.841 3932.263   4.214 2.57e-05 ***
    ## condition        60.767      3.122 3975.770  19.467  < 2e-16 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# FOR SINGLE FIXATION DURATION (SFD)

    measure = "SFD"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)

## descriptive statistic analyses for SFD

    mean.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), mean, na.rm = T)
    sd.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), sd, na.rm = T)
    se.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), std.error, na.rm = T)
    grand.mean=apply(mean.matrix, 2, mean, na.rm = T)
    grand.sd=apply(sd.matrix,2,mean, na.rm = T)
    grand.se=apply(se.matrix,2, mean, na.rm = T)
    data.frame <- data.frame(grand.mean,grand.sd,grand.se)
    data.frame

    ##    grand.mean  grand.sd grand.se
    ## SH   219.0746  61.14582 18.99782
    ## SL   217.9791  62.62147 18.32226
    ## AH   281.5049 100.33102 26.51894
    ## AL   301.2408 109.33359 29.41948

## linear mixed effect model analysis for SFD

    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    model_SFD= lmer(depvar ~readingmode*predictability+ (1+readingmode|pp) + (1|item), datafile)
    summary(model_SFD, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ readingmode * predictability + (1 + readingmode | pp) +
    ##     (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 39818.7
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -3.5783 -0.5836 -0.1664  0.4176  7.7634
    ##
    ## Random effects:
    ##  Groups   Name           Variance Std.Dev. Corr
    ##  item     (Intercept)      24.56   4.956
    ##  pp       (Intercept)    1088.67  32.995
    ##           readingmode2-1  727.10  26.965   0.29
    ##  Residual                8907.82  94.381
    ## Number of obs: 3324, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                                  Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)                       256.380      4.491   64.239  57.085  < 2e-16
    ## readingmode2-1                     70.765      4.786   61.050  14.785  < 2e-16
    ## predictability2-1                   9.471      3.318 3195.421   2.854  0.00434
    ## readingmode2-1:predictability2-1   22.191      6.633 3192.178   3.346  0.00083
    ##
    ## (Intercept)                      ***
    ## readingmode2-1                   ***
    ## predictability2-1                **
    ## readingmode2-1:predictability2-1 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    confint(model_SFD)

    ## Computing profile confidence intervals ...

    ##                                         2.5 %      97.5 %
    ## .sig01                             0.00000000  11.0799605
    ## .sig02                            26.87161806  40.2893756
    ## .sig03                            -0.08006775   0.6122112
    ## .sig04                            17.44016383  36.7390669
    ## .sigma                            92.05616474  96.7470263
    ## (Intercept)                      247.51357903 265.2526130
    ## readingmode2-1                    61.35003470  80.2945173
    ## predictability2-1                  2.96758614  15.9771012
    ## readingmode2-1:predictability2-1   9.19259435  35.2020165

## simple effect analysis for SFD

    contrast.matrix <- matrix(c(
      -1, 1, 0, 0,
      0, 0, -1, 1), 4, 2,
      dimnames=list(c("SH","SL","AH","AL"),
                    c("SL-SH", "AL-AH")))
    inv.contrast.matrix <- matrix(t(ginv(contrast.matrix)), 4, 2,
                                  dimnames=list(c("SH","SL","AH","AL"),
                                                c("SL-SH", "AL-AH")))
    (contrasts(datafile$condition) <- fractions(inv.contrast.matrix))

    ##    SL-SH AL-AH
    ## SH -1/2     0
    ## SL  1/2     0
    ## AH    0  -1/2
    ## AL    0   1/2

    depvar.lmer1 = lmer(depvar ~ condition + (1|pp)+(1|item), datafile)
    summary(depvar.lmer1, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ condition + (1 | pp) + (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 39840.8
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -3.4429 -0.5953 -0.1708  0.4168  7.7733
    ##
    ## Random effects:
    ##  Groups   Name        Variance Std.Dev.
    ##  item     (Intercept)   19.23   4.385
    ##  pp       (Intercept) 1136.58  33.713
    ##  Residual             9083.87  95.309
    ## Number of obs: 3324, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)     256.821      4.565   63.630  56.257  < 2e-16 ***
    ## conditionSL-SH   -2.293      5.001 3234.051  -0.458    0.647
    ## conditionAL-AH   20.823      4.438 3217.428   4.692 2.81e-06 ***
    ## condition        70.003      3.391 3253.868  20.643  < 2e-16 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# FOR GAZE DURATION

    measure = "GD"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)

## descriptive statistic analyses for GD

    mean.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), mean, na.rm = T)
    sd.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), sd, na.rm = T)
    se.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), std.error, na.rm = T)
    grand.mean=apply(mean.matrix, 2, mean, na.rm = T)
    grand.sd=apply(sd.matrix,2,mean, na.rm = T)
    grand.se=apply(se.matrix,2, mean, na.rm = T)
    data.frame <- data.frame(grand.mean,grand.sd,grand.se)
    data.frame

    ##    grand.mean  grand.sd grand.se
    ## SH   230.5432  76.46853 22.43798
    ## SL   233.0348  81.19390 22.37793
    ## AH   336.5712 149.33960 34.37506
    ## AL   372.4391 168.44972 38.36808

## linear mixed effect model analysis for GD

    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    model_GD= lmer(depvar ~readingmode*predictability+ (1+readingmode+predictability|pp) + (1|item), datafile)
    summary(model_GD, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula:
    ## depvar ~ readingmode * predictability + (1 + readingmode + predictability |
    ##     pp) + (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 51947.7
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -2.6197 -0.6159 -0.1959  0.4198  5.2024
    ##
    ## Random effects:
    ##  Groups   Name              Variance Std.Dev. Corr
    ##  item     (Intercept)         227.3   15.08
    ##  pp       (Intercept)        1190.8   34.51
    ##           readingmode2-1     1878.6   43.34   0.30
    ##           predictability2-1   155.5   12.47   0.41 0.63
    ##  Residual                   19601.0  140.00
    ## Number of obs: 4070, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                                  Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)                       294.970      5.158   72.032  57.183  < 2e-16
    ## readingmode2-1                    119.057      7.137   64.756  16.681  < 2e-16
    ## predictability2-1                  20.018      4.790   72.080   4.179  8.1e-05
    ## readingmode2-1:predictability2-1   34.396      9.052 3890.573   3.800 0.000147
    ##
    ## (Intercept)                      ***
    ## readingmode2-1                   ***
    ## predictability2-1                ***
    ## readingmode2-1:predictability2-1 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    confint(model_GD)

    ## Computing profile confidence intervals ...

    ##                                         2.5 %      97.5 %
    ## .sig01                             8.00361648  21.5135234
    ## .sig02                            27.36645977  42.8796063
    ## .sig03                            -0.06037386   0.9960387
    ## .sig04                            -0.61121795   0.9885329
    ## .sig05                            31.02331813  56.7247578
    ## .sig06                            -0.52800636   0.9996223
    ## .sig07                             1.31838676  26.8326195
    ## .sigma                           136.89872190 143.1955697
    ## (Intercept)                      284.75357066 305.1331272
    ## readingmode2-1                   104.94504062 133.2103173
    ## predictability2-1                 10.47920411  29.4430100
    ## readingmode2-1:predictability2-1  16.60987387  52.1664180

## simple effect analysis for GD

    contrast.matrix <- matrix(c(
      -1, 1, 0, 0,
      0, 0, -1, 1), 4, 2,
      dimnames=list(c("SH","SL","AH","AL"),
                    c("SL-SH", "AL-AH")))
    inv.contrast.matrix <- matrix(t(ginv(contrast.matrix)), 4, 2,
                                  dimnames=list(c("SH","SL","AH","AL"),
                                                c("SL-SH", "AL-AH")))
    (contrasts(datafile$condition) <- fractions(inv.contrast.matrix))

    ##    SL-SH AL-AH
    ## SH -1/2     0
    ## SL  1/2     0
    ## AH    0  -1/2
    ## AL    0   1/2

    depvar.lmer1 = lmer(depvar ~ condition + (1|pp)+(1|item), datafile)
    summary(depvar.lmer1, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ condition + (1 | pp) + (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 51991
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -2.5724 -0.6382 -0.1985  0.4322  5.3334
    ##
    ## Random effects:
    ##  Groups   Name        Variance Std.Dev.
    ##  item     (Intercept)   207.4   14.40
    ##  pp       (Intercept)  1313.2   36.24
    ##  Residual             20101.8  141.78
    ## Number of obs: 4070, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)     295.427      5.318   71.554  55.548  < 2e-16 ***
    ## conditionSL-SH    1.950      7.154 3966.522   0.272    0.785
    ## conditionAL-AH   37.807      5.694 3931.531   6.640 3.57e-11 ***
    ## condition       118.924      4.622 3985.979  25.729  < 2e-16 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# FOR TOTAL READING TIME (TRT)

    measure = "TRT"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)

## descriptive statistic analyses for TRT

    mean.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), mean, na.rm = T)
    sd.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), sd, na.rm = T)
    se.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), std.error, na.rm = T)
    grand.mean=apply(mean.matrix, 2, mean, na.rm = T)
    grand.sd=apply(sd.matrix,2,mean, na.rm = T)
    grand.se=apply(se.matrix,2, mean, na.rm = T)
    data.frame <- data.frame(grand.mean,grand.sd,grand.se)
    data.frame

    ##    grand.mean grand.sd grand.se
    ## SH   281.1793 140.1422 36.94758
    ## SL   297.6751 140.9746 35.45558
    ## AH   393.4750 187.4217 41.81485
    ## AL   431.4322 202.9030 44.88687

## linear mixed effect model analysis for TRT

    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    model_TRT= lmer(depvar ~readingmode*predictability+ (1+readingmode+predictability|pp) + (1|item), datafile)
    summary(model_TRT, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula:
    ## depvar ~ readingmode * predictability + (1 + readingmode + predictability |
    ##     pp) + (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 60734.4
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -2.9131 -0.6511 -0.2148  0.4682  8.8914
    ##
    ## Random effects:
    ##  Groups   Name              Variance Std.Dev. Corr
    ##  item     (Intercept)         835.0   28.90
    ##  pp       (Intercept)        2693.0   51.89
    ##           readingmode2-1     3048.2   55.21    0.32
    ##           predictability2-1   178.6   13.36    0.58 -0.09
    ##  Residual                   33608.4  183.33
    ## Number of obs: 4563, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                                  Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)                       351.862      7.702   83.183  45.687  < 2e-16
    ## readingmode2-1                    121.049      8.864   60.937  13.657  < 2e-16
    ## predictability2-1                  27.918      5.757   61.522   4.850 8.78e-06
    ## readingmode2-1:predictability2-1   22.665     11.016 4384.806   2.058   0.0397
    ##
    ## (Intercept)                      ***
    ## readingmode2-1                   ***
    ## predictability2-1                ***
    ## readingmode2-1:predictability2-1 *
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    confint(model_TRT)

    ## Computing profile confidence intervals ...

    ##                                         2.5 %     97.5 %
    ## .sig01                            21.55719339  37.049278
    ## .sig02                            42.16713902  63.560564
    ## .sig03                            -0.01862441   0.609707
    ## .sig04                            -1.00000000   1.000000
    ## .sig05                            39.63964585  72.110993
    ## .sig06                            -1.00000000   1.000000
    ## .sig07                             0.00000000  32.884395
    ## .sigma                           179.49515281 187.257760
    ## (Intercept)                      336.66156202 367.030158
    ## readingmode2-1                   103.54642847 138.609245
    ## predictability2-1                 16.47933690  39.246365
    ## readingmode2-1:predictability2-1   1.03545777  44.277834

## simple effect analysis for TRT

    contrast.matrix <- matrix(c(
      -1, 1, 0, 0,
      0, 0, -1, 1), 4, 2,
      dimnames=list(c("SH","SL","AH","AL"),
                    c("SL-SH", "AL-AH")))
    inv.contrast.matrix <- matrix(t(ginv(contrast.matrix)), 4, 2,
                                  dimnames=list(c("SH","SL","AH","AL"),
                                                c("SL-SH", "AL-AH")))
    (contrasts(datafile$condition) <- fractions(inv.contrast.matrix))

    ##    SL-SH AL-AH
    ## SH -1/2     0
    ## SL  1/2     0
    ## AH    0  -1/2
    ## AL    0   1/2

    depvar.lmer1 = lmer(depvar ~ condition + (1|pp)+(1|item), datafile)
    summary(depvar.lmer1, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ condition + (1 | pp) + (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 60777.3
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -2.5919 -0.6669 -0.2207  0.4728  8.7966
    ##
    ## Random effects:
    ##  Groups   Name        Variance Std.Dev.
    ##  item     (Intercept)   819.2   28.62
    ##  pp       (Intercept)  2842.6   53.32
    ##  Residual             34383.4  185.43
    ## Number of obs: 4563, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)     352.403      7.846   82.238  44.913  < 2e-16 ***
    ## conditionSL-SH   16.892      8.428 4441.558   2.004   0.0451 *
    ## conditionAL-AH   39.827      7.266 4416.623   5.481 4.47e-08 ***
    ## condition       120.714      5.589 4443.426  21.598  < 2e-16 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# FOR REGRESSION PATH DURATION (RPD)

    measure = "RPD"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)

## descriptive statistic analyses for RPD

    mean.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), mean, na.rm = T)
    sd.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), sd, na.rm = T)
    se.matrix = tapply(datafile$depvar, list(datafile$pp, datafile$condition), std.error, na.rm = T)
    grand.mean=apply(mean.matrix, 2, mean, na.rm = T)
    grand.sd=apply(sd.matrix,2,mean, na.rm = T)
    grand.se=apply(se.matrix,2, mean, na.rm = T)
    data.frame <- data.frame(grand.mean,grand.sd,grand.se)
    data.frame

    ##    grand.mean grand.sd grand.se
    ## SH   283.5736 156.5563 50.84719
    ## SL   285.5320 155.1160 46.37251
    ## AH   386.2382 213.9354 49.13321
    ## AL   421.2266 234.2873 54.23768

## linear mixed effect model analysis for RPD

    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    model_RPD= lmer(depvar ~readingmode*predictability+ (1+readingmode|pp) + (1|item), datafile)
    summary(model_RPD, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ readingmode * predictability + (1 + readingmode | pp) +
    ##     (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 55639.4
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -1.9542 -0.5366 -0.2376  0.2876 10.2545
    ##
    ## Random effects:
    ##  Groups   Name           Variance Std.Dev. Corr
    ##  item     (Intercept)      833.6   28.87
    ##  pp       (Intercept)     1829.4   42.77
    ##           readingmode2-1  4459.4   66.78   0.26
    ##  Residual                48776.1  220.85
    ## Number of obs: 4070, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                                  Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)                       341.600      7.173   79.700  47.626  < 2e-16
    ## readingmode2-1                    123.401     11.094   58.314  11.123 4.82e-16
    ## predictability2-1                  18.214      7.132 3914.594   2.554   0.0107
    ## readingmode2-1:predictability2-1   33.789     14.263 3915.269   2.369   0.0179
    ##
    ## (Intercept)                      ***
    ## readingmode2-1                   ***
    ## predictability2-1                *
    ## readingmode2-1:predictability2-1 *
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    confint(model_RPD)

    ## Computing profile confidence intervals ...

    ##                                        2.5 %      97.5 %
    ## .sig01                            18.9982454  38.8658562
    ## .sig02                            32.6169564  54.4924377
    ## .sig03                            -0.1302937   0.6278195
    ## .sig04                            46.1856411  88.7462786
    ## .sigma                           215.9667289 225.8259936
    ## (Intercept)                      327.4431633 355.7235455
    ## readingmode2-1                   101.4705135 145.2997816
    ## predictability2-1                  4.2271727  32.1890258
    ## readingmode2-1:predictability2-1   5.8541628  61.7918104

## simple effect analysis for RPD

    contrast.matrix <- matrix(c(
      -1, 1, 0, 0,
      0, 0, -1, 1), 4, 2,
      dimnames=list(c("SH","SL","AH","AL"),
                    c("SL-SH", "AL-AH")))
    inv.contrast.matrix <- matrix(t(ginv(contrast.matrix)), 4, 2,
                                  dimnames=list(c("SH","SL","AH","AL"),
                                                c("SL-SH", "AL-AH")))
    (contrasts(datafile$condition) <- fractions(inv.contrast.matrix))

    ##    SL-SH AL-AH
    ## SH -1/2     0
    ## SL  1/2     0
    ## AH    0  -1/2
    ## AL    0   1/2

    depvar.lmer1 = lmer(depvar ~ condition + (1|pp)+(1|item), datafile)
    summary(depvar.lmer1, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ condition + (1 | pp) + (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 55672.4
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -1.8395 -0.5578 -0.2374  0.2793 10.3235
    ##
    ## Random effects:
    ##  Groups   Name        Variance Std.Dev.
    ##  item     (Intercept)   770     27.75
    ##  pp       (Intercept)  2111     45.95
    ##  Residual             49807    223.18
    ## Number of obs: 4070, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                 Estimate Std. Error        df t value Pr(>|t|)
    ## (Intercept)     341.7877     7.4226   78.8630  46.047  < 2e-16 ***
    ## conditionSL-SH   -0.4522    11.2641 3966.2383  -0.040    0.968
    ## conditionAL-AH   35.8975     8.9631 3933.2888   4.005 6.32e-05 ***
    ## condition       124.0698     7.2674 3993.4780  17.072  < 2e-16 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1