---
title: "threefactors"
output: html_document
date: "2024-05-27"
---

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~ group * readingmode* predictability+(1+readingmode* predictability|pp)+(1+group* readingmode*predictability|item), datafile)
model=lmer(depvar~group* readingmode* predictability+(1+readingmode*predictability|pp)+(1|item), datafile)

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

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

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

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

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

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

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

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


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


```{r warning=FALSE}
library(lmerTest) # version 3.1-3
library(lattice)  # version 0.21-9
library(MASS)     # version 7.3-60
library(plyr)     # version 1.8.9
library(plotrix)  # version 3.8-4
```



```{r}
rm(list=ls())
datafile0 = read.csv("threefactors.csv", sep = ",", dec = ".")
datafile <- datafile0
colnames(datafile)
nrow(datafile)
```



```{r}
datafile$pp = datafile[,1]
datafile$item = datafile[,2]
datafile$group = datafile[,3]
datafile$readingmode = datafile[,4]
datafile$predictability = datafile[,5]
datafile$rm_pred=datafile[,6]
datafile$group_rm=datafile[,7]
datafile$group_pred=datafile[,8]
datafile$condition = datafile[,9]
datafile$pp = as.factor(datafile$pp)
datafile$item = as.factor(datafile$item)
datafile$group = factor(datafile$group, levels = c("K","M"))
datafile$readingmode = factor(datafile$readingmode, levels = c("S","A"))
datafile$predictability = factor(datafile$predictability, levels = c("H","L"))
datafile$rm_pred = factor(datafile$rm_pred, levels = c("SH","SL","AH","AL"))
datafile$group_rm=factor(datafile$group_rm,levels=c("KS","KA","MS","MA"))
datafile$group_pred=factor(datafile$group_pred,levels=c("KH","KL","MH","ML"))
datafile$condition=factor(datafile$condition, levels=c("KSH","KSL","KAH","KAL","MSH","MSL","MAH","MAL"))
```




# FOR ACCURACY
```{r}
measure = "accuracy"
datafile$depvar = datafile[,measure]
datafile$depvar = as.numeric(datafile$depvar)
```


## regression analysis for accuracy
```{r warning=FALSE}
contrasts(datafile$group) <- contr.sdif(2) 
contrasts(datafile$predictability) <- contr.sdif(2) 
contrasts(datafile$readingmode) <- contr.sdif(2) 
model_accuracy = glm(depvar~group*readingmode*predictability,datafile,family=binomial)
summary(model_accuracy)
```





# FOR IA_SKIP RATE
```{r}
measure = "IA_SKIP"
datafile$depvar = datafile[,measure]
datafile$depvar = as.numeric(datafile$depvar)
```


## linear mixed effect model analysis for skipping rate
```{r warning=FALSE}
contrasts(datafile$group) <- contr.sdif(2)
contrasts(datafile$predictability) <- contr.sdif(2) 
contrasts(datafile$readingmode) <- contr.sdif(2) 
model_SKIP=glmer(depvar ~ group*readingmode*predictability + (1|pp)+ (1|item),
                     datafile, family=binomial, control=glmerControl(optimizer="Nelder_Mead",optCtrl = list(maxfun = 200000)))
summary(model_SKIP)
confint(model_SKIP)
```





# FOR FIRST FIXATION DURATION (FFD)
```{r}
measure = "FFD"
datafile$depvar = datafile[,measure]
datafile$depvar = as.numeric(datafile$depvar)
```


## linear mixed effect model analysis for FFD
```{r warning=FALSE}
contrasts(datafile$group) <- contr.sdif(2)
contrasts(datafile$predictability) <- contr.sdif(2) 
contrasts(datafile$readingmode) <- contr.sdif(2) 
model_FFD= lmer(depvar ~group*readingmode*predictability+ (1+readingmode|pp) + (1|item), datafile)
summary(model_FFD, corr = FALSE)
confint(model_FFD)
```





# FOR SINGLE FIXATION DURATION (SFD)
```{r}
measure = "SFD"
datafile$depvar = datafile[,measure]
datafile$depvar = as.numeric(datafile$depvar)
```


## linear mixed effect model analysis for SFD
```{r warning=FALSE}
contrasts(datafile$group) <- contr.sdif(2)
contrasts(datafile$predictability) <- contr.sdif(2) 
contrasts(datafile$readingmode) <- contr.sdif(2) 
model_SFD= lmer(depvar ~group*readingmode*predictability+ (1+readingmode|pp) + (1|item), datafile)
summary(model_SFD, corr = FALSE)
confint(model_SFD)
```





# FOR GAZE DURATION
```{r}
measure = "GD"
datafile$depvar = datafile[,measure]
datafile$depvar = as.numeric(datafile$depvar)
```


## linear mixed effect model analysis for GD
```{r warning=FALSE}
contrasts(datafile$group) <- contr.sdif(2)
contrasts(datafile$predictability) <- contr.sdif(2) 
contrasts(datafile$readingmode) <- contr.sdif(2) 
model_GD= lmer(depvar ~group*readingmode*predictability+ (1+readingmode+predictability|pp) + (1|item), datafile)
summary(model_GD, corr = FALSE)
confint(model_GD)
```





# FOR TOTAL READING TIME (TRT)
```{r}
measure = "TRT"
datafile$depvar = datafile[,measure]
datafile$depvar = as.numeric(datafile$depvar)
```


## linear mixed effect model analysis for TRT
```{r warning=FALSE}
contrasts(datafile$group) <- contr.sdif(2)
contrasts(datafile$predictability) <- contr.sdif(2) 
contrasts(datafile$readingmode) <- contr.sdif(2) 
model_TRT= lmer(depvar ~group*readingmode*predictability+ (1+readingmode+predictability|pp) + (1|item), datafile)
summary(model_TRT, corr = FALSE)
confint(model_TRT)
```





# FOR REGRESSION PATH DURATION (RPD)
```{r}
measure = "RPD"
datafile$depvar = datafile[,measure]
datafile$depvar = as.numeric(datafile$depvar)
```


## linear mixed effect model analysis for RPD
```{r warning=FALSE}
contrasts(datafile$group) <- contr.sdif(2)
contrasts(datafile$predictability) <- contr.sdif(2) 
contrasts(datafile$readingmode) <- contr.sdif(2) 
model_RPD= lmer(depvar ~group*readingmode*predictability+ (1+readingmode|pp) + (1|item), datafile)
summary(model_RPD, corr = FALSE)
confint(model_RPD)
```


