# ancova analysis

#### 2024-05-27

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

    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("ancova_speed.csv", sep = ",", dec = ".")
    datafile <- datafile0
    colnames(datafile)

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

    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$speed=datafile[,6]
    datafile$speed <- scale(datafile$speed)
    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"))

# ancova analysis for skipping rate

    measure = "IA_SKIP"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)
    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    ancova_SKIP=glmer(depvar~readingmode*predictability+speed + (1|pp)+ (1|item),
                         datafile, family=binomial, control=glmerControl(optimizer="Nelder_Mead",optCtrl = list(maxfun = 200000)))
    summary(ancova_SKIP)

    ## Generalized linear mixed model fit by maximum likelihood (Laplace
    ##   Approximation) [glmerMod]
    ##  Family: binomial  ( logit )
    ## Formula: depvar ~ readingmode * predictability + speed + (1 | pp) + (1 |
    ##     item)
    ##    Data: datafile
    ## Control:
    ## glmerControl(optimizer = "Nelder_Mead", optCtrl = list(maxfun = 2e+05))
    ##
    ##      AIC      BIC   logLik deviance df.resid
    ##   5242.0   5288.2  -2614.0   5228.0     5465
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -2.2913 -0.5371 -0.3269  0.5079  5.8167
    ##
    ## Random effects:
    ##  Groups Name        Variance Std.Dev.
    ##  item   (Intercept) 0.0186   0.1364
    ##  pp     (Intercept) 0.3229   0.5683
    ## Number of obs: 5472, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                                  Estimate Std. Error z value Pr(>|z|)
    ## (Intercept)                      -1.32595    0.08239 -16.094  < 2e-16 ***
    ## readingmode2-1                   -1.30340    0.11104 -11.738  < 2e-16 ***
    ## predictability2-1                -0.16631    0.07485  -2.222   0.0263 *
    ## speed                             0.38856    0.05282   7.357 1.88e-13 ***
    ## readingmode2-1:predictability2-1 -0.10571    0.14964  -0.706   0.4799
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Correlation of Fixed Effects:
    ##             (Intr) rdn2-1 prd2-1 speed
    ## readngmd2-1  0.133
    ## prdctblt2-1  0.029  0.059
    ## speed        0.003  0.723  0.024
    ## rdng2-1:2-1  0.025  0.029  0.374 -0.014

    anova(ancova_SKIP)

    ## Analysis of Variance Table
    ##                            npar Sum Sq Mean Sq  F value
    ## readingmode                   1 608.54  608.54 608.5371
    ## predictability                1   5.37    5.37   5.3696
    ## speed                         1  57.31   57.31  57.3051
    ## readingmode:predictability    1   0.50    0.50   0.4960

    confint(ancova_SKIP)

    ## Computing profile confidence intervals ...

    ##                                       2.5 %     97.5 %
    ## .sig01                            0.0000000  0.2541007
    ## .sig02                            0.4510799  0.7186415
    ## (Intercept)                      -1.4927297 -1.1645627
    ## readingmode2-1                   -1.5234851 -1.0866246
    ## predictability2-1                -0.3138968 -0.0193142
    ## speed                             0.2851970  0.4924500
    ## readingmode2-1:predictability2-1 -0.4006721  0.1882777

# ancova analysis for FFD

    measure = "FFD"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)
    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    ancova_FFD= lmer(depvar ~readingmode*predictability+ speed+(1+readingmode|pp) + (1|item), datafile)
    summary(ancova_FFD, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ readingmode * predictability + speed + (1 + readingmode |
    ##     pp) + (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 48740.8
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -3.5346 -0.5939 -0.1487  0.4244  7.6942
    ##
    ## Random effects:
    ##  Groups   Name           Variance Std.Dev. Corr
    ##  item     (Intercept)      54.5    7.383
    ##  pp       (Intercept)     920.3   30.337
    ##           readingmode2-1  572.3   23.923   0.40
    ##  Residual                8936.2   94.532
    ## Number of obs: 4070, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                                  Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)                       252.888      4.180   65.800  60.498  < 2e-16
    ## readingmode2-1                     37.666      5.823  118.146   6.468 2.34e-09
    ## predictability2-1                   7.399      3.051 3926.310   2.425  0.01534
    ## speed                             -18.949      3.034  543.869  -6.245 8.57e-10
    ## readingmode2-1:predictability2-1   16.871      6.101 3926.964   2.765  0.00571
    ##
    ## (Intercept)                      ***
    ## readingmode2-1                   ***
    ## predictability2-1                *
    ## speed                            ***
    ## readingmode2-1:predictability2-1 **
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    anova(ancova_FFD)

    ## Type III Analysis of Variance Table with Satterthwaite's method
    ##                            Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)
    ## readingmode                373901  373901     1  118.1 41.8412 2.336e-09 ***
    ## predictability              52559   52559     1 3926.3  5.8816  0.015345 *
    ## speed                      348465  348465     1  543.9 38.9947 8.570e-10 ***
    ## readingmode:predictability  68340   68340     1 3927.0  7.6476  0.005711 **
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    confint(ancova_FFD)

    ## Computing profile confidence intervals ...

    ##                                         2.5 %      97.5 %
    ## .sig01                             0.00000000  12.0990825
    ## .sig02                            24.63215307  37.1010864
    ## .sig03                             0.03262684   0.7382762
    ## .sig04                            14.83796733  32.9152707
    ## .sigma                            92.43198470  96.6469833
    ## (Intercept)                      244.65299619 261.1493622
    ## readingmode2-1                    26.24367256  49.0974694
    ## predictability2-1                  1.42275559  13.3818155
    ## speed                            -24.90496025 -13.0043944
    ## readingmode2-1:predictability2-1   4.91368660  28.8295499

# ancova analysis for SFD

    measure = "SFD"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)
    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    ancova_SFD= lmer(depvar ~readingmode*predictability+ speed+(1+readingmode|pp) + (1|item), datafile)
    summary(ancova_SFD, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ readingmode * predictability + speed + (1 + readingmode |
    ##     pp) + (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 39780.9
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -3.5750 -0.5881 -0.1600  0.4189  7.7233
    ##
    ## Random effects:
    ##  Groups   Name           Variance Std.Dev. Corr
    ##  item     (Intercept)      31.29   5.594
    ##  pp       (Intercept)     976.41  31.248
    ##           readingmode2-1  621.86  24.937   0.42
    ##  Residual                8845.44  94.050
    ## Number of obs: 3324, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                                  Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)                       256.001      4.296   64.276  59.589  < 2e-16
    ## readingmode2-1                     46.448      6.138  128.595   7.567 6.47e-12
    ## predictability2-1                   9.346      3.306 3199.444   2.827 0.004728
    ## speed                             -18.338      3.106  569.544  -5.904 6.11e-09
    ## readingmode2-1:predictability2-1   22.180      6.607 3195.381   3.357 0.000797
    ##
    ## (Intercept)                      ***
    ## readingmode2-1                   ***
    ## predictability2-1                **
    ## speed                            ***
    ## readingmode2-1:predictability2-1 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    anova(ancova_SFD)

    ## Type III Analysis of Variance Table with Satterthwaite's method
    ##                            Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)
    ## readingmode                506484  506484     1  128.6  57.259 6.467e-12 ***
    ## predictability              70693   70693     1 3199.4   7.992 0.0047276 **
    ## speed                      308310  308310     1  569.5  34.855 6.106e-09 ***
    ## readingmode:predictability  99679   99679     1 3195.4  11.269 0.0007973 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    confint(ancova_SFD)

    ## Computing profile confidence intervals ...

    ##                                         2.5 %     97.5 %
    ## .sig01                             0.00000000  11.494140
    ## .sig02                            25.31760994  38.246681
    ## .sig03                             0.04277045   0.750282
    ## .sig04                            15.20186878  34.444538
    ## .sigma                            91.72513546  96.396167
    ## (Intercept)                      247.52528150 264.477656
    ## readingmode2-1                    34.38734368  58.490440
    ## predictability2-1                  2.86975722  15.830425
    ## speed                            -24.44673416 -12.247689
    ## readingmode2-1:predictability2-1   9.22358414  35.129417

# ancova analysis for GD

    measure = "GD"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)
    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    ancova_GD= lmer(depvar ~readingmode*predictability+speed+ (1+readingmode|pp) + (1|item), datafile)
    summary(ancova_GD, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ readingmode * predictability + speed + (1 + readingmode |
    ##     pp) + (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 51894.5
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -2.5939 -0.6239 -0.1896  0.4308  5.2785
    ##
    ## Random effects:
    ##  Groups   Name           Variance Std.Dev. Corr
    ##  item     (Intercept)      233.2   15.27
    ##  pp       (Intercept)      834.1   28.88
    ##           readingmode2-1  1614.3   40.18   0.56
    ##  Residual                19512.4  139.69
    ## Number of obs: 4070, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                                  Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)                       293.900      4.584   72.041  64.110  < 2e-16
    ## readingmode2-1                     77.231      8.685  122.791   8.892 6.31e-15
    ## predictability2-1                  20.209      4.507 3924.866   4.484 7.55e-06
    ## speed                             -32.280      4.290  420.554  -7.525 3.24e-13
    ## readingmode2-1:predictability2-1   34.000      9.014 3924.705   3.772 0.000164
    ##
    ## (Intercept)                      ***
    ## readingmode2-1                   ***
    ## predictability2-1                ***
    ## speed                            ***
    ## readingmode2-1:predictability2-1 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    anova(ancova_GD)

    ## Type III Analysis of Variance Table with Satterthwaite's method
    ##                             Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)
    ## readingmode                1542911 1542911     1  122.8  79.073 6.307e-15 ***
    ## predictability              392252  392252     1 3924.9  20.103 7.549e-06 ***
    ## speed                      1104802 1104802     1  420.6  56.621 3.239e-13 ***
    ## readingmode:predictability  277615  277615     1 3924.7  14.228 0.0001644 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    confint(ancova_GD)

    ## Computing profile confidence intervals ...

    ##                                        2.5 %      97.5 %
    ## .sig01                             8.2522629  21.6806769
    ## .sig02                            22.1767500  36.4988911
    ## .sig03                             0.1910964   0.8677244
    ## .sig04                            27.9686174  53.0737870
    ## .sigma                           136.5865290 142.8170090
    ## (Intercept)                      284.8485566 302.9148449
    ## readingmode2-1                    60.0120443  94.2808268
    ## predictability2-1                 11.3802288  29.0540498
    ## speed                            -40.7899471 -23.7911864
    ## readingmode2-1:predictability2-1  16.3195910  51.6681150

# ancova analysis for TRT

    measure = "TRT"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)
    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    ancova_TRT= lmer(depvar ~readingmode*predictability+speed+ (1+readingmode+predictability|pp) + (1|item), datafile)
    summary(ancova_TRT, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ readingmode * predictability + speed + (1 + readingmode +
    ##     predictability | pp) + (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 60512.1
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -2.9195 -0.6732 -0.1757  0.4885  8.7410
    ##
    ## Random effects:
    ##  Groups   Name              Variance Std.Dev. Corr
    ##  item     (Intercept)         733.4   27.08
    ##  pp       (Intercept)        1582.3   39.78
    ##           readingmode2-1     2627.3   51.26   0.65
    ##           predictability2-1   191.1   13.82   0.22 0.44
    ##  Residual                   32333.4  179.81
    ## Number of obs: 4563, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                                  Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)                       349.896      6.362   85.826  54.994  < 2e-16
    ## readingmode2-1                     21.303     10.517  112.170   2.026   0.0452
    ## predictability2-1                  27.581      5.672   61.697   4.863 8.32e-06
    ## speed                             -75.346      4.878  403.775 -15.447  < 2e-16
    ## readingmode2-1:predictability2-1   22.088     10.800 4395.336   2.045   0.0409
    ##
    ## (Intercept)                      ***
    ## readingmode2-1                   *
    ## predictability2-1                ***
    ## speed                            ***
    ## readingmode2-1:predictability2-1 *
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    anova(ancova_TRT)

    ## Type III Analysis of Variance Table with Satterthwaite's method
    ##                             Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)
    ## readingmode                 132654  132654     1  112.2   4.1027   0.04519 *
    ## predictability              764655  764655     1   61.7  23.6491 8.323e-06 ***
    ## speed                      7714585 7714585     1  403.8 238.5949 < 2.2e-16 ***
    ## readingmode:predictability  135248  135248     1 4395.3   4.1829   0.04089 *
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    confint(ancova_TRT)

    ## Computing profile confidence intervals ...

    ##                                         2.5 %      97.5 %
    ## .sig01                            19.82175855  34.9201538
    ## .sig02                            31.44802877  49.5235749
    ## .sig03                             0.31484329   0.9153431
    ## .sig04                            -0.96170870   1.0000000
    ## .sig05                            35.83540044  67.4592167
    ## .sig06                            -0.90532467   0.9583859
    ## .sig07                             0.60316519  32.5783145
    ## .sigma                           176.04670490 183.6688545
    ## (Intercept)                      337.33536367 362.4031078
    ## readingmode2-1                    -0.09922848  42.3231837
    ## predictability2-1                 16.31312315  38.7645989
    ## speed                            -85.37348783 -65.3640263
    ## readingmode2-1:predictability2-1   0.81542622  43.2453914

# ancova analysis for RPD

    measure = "RPD"
    datafile$depvar = datafile[,measure]
    datafile$depvar = as.numeric(datafile$depvar)
    contrasts(datafile$predictability) <- contr.sdif(2)
    contrasts(datafile$readingmode) <- contr.sdif(2)
    ancova_RPD= lmer(depvar ~readingmode*predictability+speed+ (1+readingmode|pp) + (1|item), datafile)
    summary(ancova_RPD, corr = FALSE)

    ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
    ## lmerModLmerTest]
    ## Formula: depvar ~ readingmode * predictability + speed + (1 + readingmode |
    ##     pp) + (1 | item)
    ##    Data: datafile
    ##
    ## REML criterion at convergence: 55555.8
    ##
    ## Scaled residuals:
    ##     Min      1Q  Median      3Q     Max
    ## -1.9420 -0.5613 -0.2250  0.2880 10.3258
    ##
    ## Random effects:
    ##  Groups   Name           Variance Std.Dev. Corr
    ##  item     (Intercept)      852.1   29.19
    ##  pp       (Intercept)     1119.3   33.46
    ##           readingmode2-1  3504.7   59.20   0.55
    ##  Residual                48181.5  219.50
    ## Number of obs: 4070, groups:  item, 88; pp, 64
    ##
    ## Fixed effects:
    ##                                  Estimate Std. Error       df t value Pr(>|t|)
    ## (Intercept)                       338.069      6.345   78.903  53.279  < 2e-16
    ## readingmode2-1                     47.697     13.160  103.992   3.624  0.00045
    ## predictability2-1                  17.826      7.082 3921.034   2.517  0.01187
    ## speed                             -60.343      6.570  316.032  -9.185  < 2e-16
    ## readingmode2-1:predictability2-1   33.565     14.164 3920.725   2.370  0.01785
    ##
    ## (Intercept)                      ***
    ## readingmode2-1                   ***
    ## predictability2-1                *
    ## speed                            ***
    ## readingmode2-1:predictability2-1 *
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    anova(ancova_RPD)

    ## Type III Analysis of Variance Table with Satterthwaite's method
    ##                             Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)
    ## readingmode                 632911  632911     1  104.0 13.1360 0.0004502 ***
    ## predictability              305253  305253     1 3921.0  6.3355 0.0118742 *
    ## speed                      4064926 4064926     1  316.0 84.3669 < 2.2e-16 ***
    ## readingmode:predictability  270586  270586     1 3920.7  5.6160 0.0178456 *
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    confint(ancova_RPD)

    ## Computing profile confidence intervals ...

    ##                                         2.5 %      97.5 %
    ## .sig01                            19.38213057  39.0379890
    ## .sig02                            23.46172832  44.2362031
    ## .sig03                             0.08733848   0.9664572
    ## .sig04                            38.56389725  80.0998571
    ## .sigma                           214.63421901 224.4382726
    ## (Intercept)                      325.53402031 350.5754779
    ## readingmode2-1                    21.48078780  73.5691896
    ## predictability2-1                  3.94803815  31.7140281
    ## speed                            -73.34145339 -47.4018591
    ## readingmode2-1:predictability2-1   5.80385427  61.3509279