---
title: "reading_analysis"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Overview

This notebook tests the effect of visual noise on reading performance of four groups with differen development. The readin material consists of a line semantically disconnected words or pseudowords with either 3 or 7 letters. The reading performance is evaluated by the time for reading a single item and the proportion of wrongly read items.

# Setup

```{r, warning=FALSE,message=FALSE,echo=FALSE}
library(plyr)
library(dplyr)
library(ggplot2)
library(lme4)
library(car)
library(glmmTMB)
library(sjPlot)
library(DHARMa)
library(ggpubr)

library(bbmle)
library(effects)
library(emmeans)
library(future)
library(readxl)
library(vioplot)

plan(multiprocess)

knitr::opts_chunk$set(fig.width=16, fig.height=16)

```
#Data preparation
```{r}

data_reading <- read_excel("data_reading4.xlsx")

# coerce variables that are categorical to factors:

data_reading$group<-as.factor(data_reading$group)
data_reading$word_type<-as.factor(data_reading$word_type)
data_reading$subject<-as.factor(data_reading$subject)
data_reading$letter_number<-as.factor(data_reading$letter_number)

#recode group names and word type levels
data_reading$group <- revalue(data_reading$group, c("1" = "ASD", "2" = "TD","3" = "ADHD", "4" = "SLD"))
data_reading$word_type <- revalue(data_reading$word_type, c("1" = "words", "2" = "pseudowords"))

#recode values with reading time = 0 (no line read) as missing

data_reading$time_per_word[which(data_reading$time_per_word==0)]<-NA
data_reading$time[which(data_reading$time==0)]<-NA

#calculate the proportion oferrorsfrom the read wordss
data_reading$prop_errors <- data_reading$errors/(data_reading$word_number-data_reading$misses)
data_reading$noise0 <- as.factor(data_reading$noise)
#exclusion of cases where no string was read
dat<-data_reading[!is.na(data_reading$time),]

```
#Test reading speed

#time_per_word

# using lme4 to select distribution

```{r,warning=FALSE,results='hide'}
#random intercept, inverse gaussian family, and identity link
fitt1<-glmer(time_per_word~group*letter_number*word_type*noise0+(1|subject),data=dat,family=inverse.gaussian(link="identity"),control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#log link
fitt2<-glmer(time_per_word~group*letter_number*word_type*noise0+(1|subject),dat=dat,family=inverse.gaussian(link="log"),na.action = na.omit,control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
car::Anova(fitt2)
anova(fitt1,fitt2)

#random intercept, Gamma family, and identity link

fitt3<-glmer(time_per_word~group*letter_number*word_type*noise0+(1|subject),data=dat,family=Gamma(link="identity"),na.action = na.omit,control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
anova(fitt2,fitt3)
fitt4<-glmer(time_per_word~group*letter_number*word_type*noise0+(1|subject),data=dat,family=Gamma(link="log"),na.action = na.omit,control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
```


```{r,warning=FALSE}
fitt2a<-glmer(time_per_word~group*letter_number*word_type*noise0+(letter_number+word_type+noise0||subject),dat=dat,family=inverse.gaussian(link="log"),na.action = na.omit)

fitt2b<-update(fitt2a,formula = ~ . -group:letter_number:word_type )

fitt2c<-update(fitt2b,formula = ~ . +scale(AGE) )
AICtab(fitt2a,fitt2b,fitt2c,fitt2,fitt1,fitt3,fitt4)
#simulationOutput<-simulateResiduals(fitt2c)
#plot(simulationOutput)
#simulationOutput<-simulateResiduals(fitt4)
#plot(simulationOutput)

```
#Gamma distribution with log link is better than with identity link
#using glmmTNB with varying dispersion parameter
```{r,warning=FALSE}
fit0a<-glmmTMB(time_per_word~group*letter_number*word_type*noise0+scale(AGE)+scale(FSIQ)+(1|subject),data=dat,family=Gamma(link='log'),na.action = na.omit,dispformula=~letter_number+noise0+word_type)
fit0b<-glmmTMB(time_per_word~group*letter_number*word_type*noise0+scale(AGE)+(noise0 + letter_number + word_type|subject),data=dat,family=Gamma(link='log'),na.action = na.omit,dispformula=~letter_number+noise0+word_type)

summary(fit0b)
simulationOutput<-simulateResiduals(fit0b,plot = T)
testOutliers(simulationOutput)
testDispersion(fit0b)
Anova(fit0b)
joint_tests(fit0b)

AICtab(fit0a,fit0b)


fit0d<-update(fit0b,dispformula=~group+noise0+letter_number+word_type,control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))
fit0e<-update(fit0d,formula=~.-group:letter_number:word_type-letter_number:word_type:noise0)

AICtab(fit0a,fit0b,fit0d,fit0e)
Anova(fit0e)
joint_tests(fit0e)
tab_model(fit0e, show.reflvl = TRUE)
plot_model(fit0e,type = "pred",terms = "group")
```
#for Figure 6
```{r}
w<-as.data.frame(allEffects(fit0b))[[2]]
w$noise<-rep(c(0,50,55,60,62.5),each=16)
w$gr<-factor(w$noise>0)
ggplot(w,aes((noise0),fit,group=interaction(group,gr),color=group))+geom_point()+geom_line()+facet_grid(word_type~letter_number)+geom_errorbar(aes(ymin=fit-se,ymax=fit+se),width=.2)+
scale_x_discrete(name="Noise [deg]")+ scale_y_continuous(name="Time per string (fitted)[s]")+
theme(axis.title.x = element_text(size=15), axis.text.x  = element_text(size=14))+
theme(axis.title.y = element_text(size=16), axis.text.y  = element_text(size=14),strip.text.x = element_text(size = 12),strip.text.y = element_text(size = 12))+
theme(legend.text=element_text(size=10))+
theme(legend.title = element_text(size=14))+ theme(legend.title = element_blank())
pairs(emmeans(fit0b,~group))
pairs(emmeans(fit0b,~group|noise0))
pairs(emmeans(fit0b,~word_type|group))
```
#For figure 7
```{r}
t62<-ggplot(subset(w,noise==62.5),aes(word_type,fit,fill=group))+geom_bar(stat="identity", color="black", alpha=0.7,
position=position_dodge()) +theme(axis.title.x = element_text(size=15), axis.text.x  = element_text(size=14))+theme(axis.title.y = element_text(size=15), axis.text.y  = element_text(size=14))+scale_x_discrete(name="letter string type")+
geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=.2,
position=position_dodge(.9))+facet_grid(~letter_number)+scale_y_continuous(name="Time per string [s] (fitted)",limits=c(0, 5))+ggtitle("Noise=62.5")+ theme(legend.title = element_blank())
t62
t0<-ggplot(subset(w,noise==0),aes(word_type,fit,fill=group))+geom_bar(stat="identity", color="black", alpha=0.7,
position=position_dodge()) +theme(axis.title.x = element_text(size=15), axis.text.x  = element_text(size=14))+theme(axis.title.y = element_text(size=15), axis.text.y  = element_text(size=14))+scale_x_discrete(name="letter string type")+
geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=.2,
position=position_dodge(.9))+facet_grid(~letter_number)+scale_y_continuous(name="Time per string [s] (fitted)",limits=c(0, 5))+ggtitle("Noise=0")+ theme(legend.title = element_blank())
figure1 <- ggarrange(t0, t62,
labels = c("A", "B"),
common.legend = TRUE, legend = "bottom")
figure1


```
```
#Accuracy analysis
```{r}
f<-glmmTMB(errors~group*noise*word_type*letter_number+ (1|subject)+ offset(log(word_number-misses)),ziformula = ~1, data=dat, family =  poisson)
f1<-glmmTMB(errors~group*noise*word_type*letter_number+ (1|subject)+ offset(log(word_number-misses)),ziformula = ~noise0, data=dat, family =  poisson)
f2<-update(f,ziformula=~letter_number)
f3<-update(f,ziformula=~letter_number+word_type)
f4<-update(f,ziformula=~letter_number+noise0,control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS"),optCtrl=list(iter.max=1e3,eval.max=1e3)))
fa<-glmmTMB(errors~group*noise0*word_type*letter_number+ (1|subject)+ offset(log(word_number-misses)),ziformula = ~1, data=dat, family =  poisson)
f1a<-glmmTMB(errors~group*noise0*word_type*letter_number+ (1|subject)+ offset(log(word_number-misses)),ziformula = ~noise0, data=dat, family =  poisson)
f2a<-update(fa,ziformula=~letter_number)
f3a<-update(fa,ziformula=~letter_number+word_type)
f4a<-update(fa,ziformula=~letter_number+noise0,control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS"),optCtrl=list(iter.max=1e3,eval.max=1e3)))

w0<-as.data.frame(allEffects(f4a))[[1]]
w0$noise<-rep(c(0,50,55,60,62.5),each=4)
w0$gr<-factor(w0$noise>0)

ggplot(w0,aes((noise0),fit,group=interaction(gr,group),color=group))+geom_point()+geom_line()+facet_grid(word_type~letter_number)+geom_errorbar(aes(ymin=fit-se,ymax=fit+se),width=.2)+
  scale_x_discrete(name="Noise [deg]")+ scale_y_continuous(name="Error number")+
  theme(axis.title.x = element_text(size=15), axis.text.x  = element_text(size=14))+
  theme(axis.title.y = element_text(size=16), axis.text.y  = element_text(size=14),strip.text.x = element_text(size = 12),strip.text.y = element_text(size = 12))+ 
  theme(legend.text=element_text(size=10))+
  theme(legend.title = element_text(size=14))
  


```


#relation precision-speed
```{r}

b0 = glmer(prop_errors~scale(time_per_word)*group*noise*word_type*letter_number+ (1|subject), weight=word_number-misses, data=dat, family = binomial(link="logit"),control=glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun = 100000)))
b1 = glmer(prop_errors~scale(time_per_word)*group*noise*word_type*letter_number+ (scale(time_per_word)+1|subject), weight=word_number-misses, data=dat, family = binomial(link="logit"),control=glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun = 100000)))
anova(b0,b1)
simulationOutput<-simulateResiduals(b1,plot=T)
testOutliers(simulationOutput,type = "bootstrap")
testResiduals(simulationOutput)
b2 = glmer(prop_errors~log(time_per_word)*group*noise*word_type*letter_number+ (1|subject), weight=word_number-misses, data=dat, family = binomial(link="logit"),control=glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun = 100000)))
anova(b0,b1,b2)
simulationOutput<-simulateResiduals(b2,plot=T)
testResiduals(simulationOutput)
b3 = glmer(prop_errors~log(time_per_word)*group*noise*word_type*letter_number+ (log(time_per_word)+1|subject), weight=word_number-misses, data=dat, family = binomial(link="logit"),control=glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun = 100000)))

testOutliers(simulationOutput,type = "bootstrap")

b4 = glmmTMB(prop_errors~log(time_per_word)*group*noise*word_type*letter_number+ (log(time_per_word)+1|subject), weight=word_number-misses, data=dat, family = betabinomial,control=glmmTMBControl(optCtrl=list(iter.max=1e3,eval.max=1e3)))
AICtab(b0,b1,b2,b3,b4)
simulationOutput<-simulateResiduals(b4,plot=T)

testOutliers(simulationOutput,type = "bootstrap")
summary(b4)
tab_model(b4)
anova(b3,b4)
testZeroInflation(b4)

Anova(b4)
b5<-glmmTMB(prop_errors~log(time_per_word)*group*noise0+log(time_per_word)*group*word_type+log(time_per_word)*group*letter_number+ group*noise0*word_type+group*noise0*letter_number+(log(time_per_word)+1|subject), weight=word_number-misses, data=dat, family = betabinomial,control=glmmTMBControl(optCtrl=list(iter.max=1e3,eval.max=1e3)))

#for Figure 7 (restricting ploting values to maximal time per word + 1 sec)
wpt<-as.data.frame(allEffects(b4,xlevels=list(time_per_word=1:20)))[[1]]
a<-aggregate(time_per_word~group+noise0+word_type+letter_number,dat,max)
a$time_per_word<-a$time_per_word+1
a<-(a %>% slice(rep(row_number(), each=20)))
wpt0<-(filter(wpt,wpt$time_per_word<a$time_per_word))
wpt0$noise0<-as.factor(wpt0$noise)
pt1<-ggplot(subset(wpt0,word_type=="words"),aes((time_per_word),1-fit,color=noise0))+geom_line(size=1)+facet_grid(group~letter_number)+
   scale_x_continuous(name="Time per word [s]")+ scale_y_continuous(name="Proportion of correct words (fitted)")+
   theme(axis.title.x = element_text(size=15), axis.text.x  = element_text(size=14))+
   theme(axis.title.y = element_text(size=16), axis.text.y  = element_text(size=14),strip.text.x = element_text(size = 12),strip.text.y = element_text(size = 12))+
   theme(legend.text=element_text(size=10))+
   theme(legend.title = element_text(size=14))+ 
  scale_color_discrete(name = "Noise")
pt2<-ggplot(subset(wpt0,word_type=="pseudowords"),aes((time_per_word),1-fit,color=noise0))+geom_line(size=1)+facet_grid(group~letter_number)+
  scale_x_continuous(name="Time per word [s]")+ scale_y_continuous(name="Proportion of correct pseudowords (fitted)")+
  theme(axis.title.x = element_text(size=15), axis.text.x  = element_text(size=14))+
  theme(axis.title.y = element_text(size=16), axis.text.y  = element_text(size=14),strip.text.x = element_text(size = 12),strip.text.y = element_text(size = 12))+
  theme(legend.text=element_text(size=10))+
  theme(legend.title = element_text(size=14))+ 
  scale_color_discrete(name = "Noise")
figure <- ggarrange(pt1, pt2,
labels = c("A", "B"),
common.legend = TRUE, legend = "bottom")
figure



```

