---
title: "007_endAnalysis_final"
author: "Katharine Aveni"
date: "October 1, 2021"
output: html_document
---

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

##Load Libraries & Data
```{r inital data loading}

#sentence_window_clean <- readRDS("007_AVT_data_clean.rds", refhook = NULL)
sentence_window_clean <- readRDS("007_AVT_data_clean_with_anticipated.rds", refhook = NULL)

#load required packages
library("Matrix")
library("plyr")
library("tidyr")
library("lme4")
library("ggplot2")
library("eyetrackingR")
library("dplyr")
```

##Verb + Object (end) window
### zoom in to end window (already offset by 200ms)
```{r}
end_window <- subset_by_window(sentence_window_clean,
                                window_start_time = 658, 
                                window_end_time = 1908,
                                rezero = TRUE, remove = TRUE)
```
Reduced end window by 42ms in order to avoid final small bin.
Alternatively, could start the window a little bit later to capture end of target at expense of onset, or could re-export data to be able to extend the window another 46ms.

###create object window time series
```{r Create verb window time series}

end_time_seq <-make_time_sequence_data(end_window, time_bin_size = 50,
                                   predictor_columns = c("Group", "trial_name"),
                                   aois = c("Target", "ActionDis", "AgentDis", "Unrelated")
)


```

###Interpolate values for set 9 trials for subjects 2-6.

```{r Create vectors for interpolated trials (problem trials)}

#SUBJECT 007_05
#interpolate trial 9_3
`ei_05_3_2` <- subset(end_time_seq, Subject == "007_05" & trial_name == "ExpA_Set3_Item2", 
                    select = c("Prop"))
`ei_05_7_3` <- subset(end_time_seq, Subject == "007_05" & trial_name == "ExpA_Set7_Item3", 
                    select = c("Prop"))
end_time_seq[which(end_time_seq$Subject == "007_05" & end_time_seq$trial_name == "ExpA_Set9_Item3"), "Prop"] <- ((`ei_05_3_2` + `ei_05_7_3`)/2)


#interpolate trial 9_2
"ei_05_7_2" <- subset(end_time_seq, Subject == "007_05" & trial_name == "ExpA_Set7_Item2", 
                    select = c("Prop"))
"ei_05_12_1" <- subset(end_time_seq, Subject == "007_05" & trial_name == "ExpA_Set12_Item1", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_05" & end_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((ei_05_7_2 + ei_05_12_1) / 2)

#interpolate trial 9_4
"ei_05_11_4" <- subset(end_time_seq, Subject == "007_05" & trial_name == "ExpA_Set11_Item4", 
                    select = c("Prop"))
"ei_05_12_2" <- subset(end_time_seq, Subject == "007_05" & trial_name == "ExpA_Set12_Item2", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_05" & end_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((ei_05_11_4 + ei_05_12_2) / 2)

#SUBJECT 007_02
#interpolate trial 9_4
"ei_02_10_2" <- subset(end_time_seq, Subject == "007_02" & trial_name == "ExpA_Set10_Item2", 
                    select = c("Prop"))
"ei_02_3_2" <- subset(end_time_seq, Subject == "007_02" & trial_name == "ExpA_Set3_Item2", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_02" & end_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((ei_02_10_2 + ei_02_3_2) / 2)

#interpolate trial 9_1
"ei_02_3_2" <- subset(end_time_seq, Subject == "007_02" & trial_name == "ExpA_Set3_Item2", 
                    select = c("Prop"))
"ei_02_5_3" <- subset(end_time_seq, Subject == "007_02" & trial_name == "ExpA_Set5_Item3", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_02" & end_time_seq$trial_name == "ExpA_Set9_Item1"), "Prop"] <- ((ei_02_3_2 + ei_02_5_3) / 2)

#interpolate trial 9_2
"ei_02_3_4" <- subset(end_time_seq, Subject == "007_02" & trial_name == "ExpA_Set3_Item4", 
                    select = c("Prop"))
"ei_02_5_1" <- subset(end_time_seq, Subject == "007_02" & trial_name == "ExpA_Set5_Item1", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_02" & end_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((ei_02_3_4 + ei_02_5_1) / 2)

#SUBJECT 007_06
#interpolate trial 9_4
"ei_06_2_1" <- subset(end_time_seq, Subject == "007_06" & trial_name == "ExpA_Set2_Item1", 
                    select = c("Prop"))
"ei_06_7_2" <- subset(end_time_seq, Subject == "007_06" & trial_name == "ExpA_Set7_Item2", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_06" & end_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((ei_06_2_1 + ei_06_7_2) / 2)

#interpolate trial 9_2
"ei_06_7_4" <- subset(end_time_seq, Subject == "007_06" & trial_name == "ExpA_Set7_Item4", 
                    select = c("Prop"))
"ei_06_4_2" <- subset(end_time_seq, Subject == "007_06" & trial_name == "ExpA_Set4_Item2", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_06" & end_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((ei_06_7_4 + ei_06_4_2) / 2)

#SUBJECT 007_03
#interpolate trial 9_4
"ei_03_11_2" <- subset(end_time_seq, Subject == "007_03" & trial_name == "ExpA_Set11_Item2", 
                    select = c("Prop"))
"ei_03_12_1" <- subset(end_time_seq, Subject == "007_03" & trial_name == "ExpA_Set12_Item1", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_03" & end_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((ei_03_11_2 + ei_03_12_1) / 2)

#interpolate trial 9_2
"ei_03_6_2" <- subset(end_time_seq, Subject == "007_03" & trial_name == "ExpA_Set6_Item2", 
                    select = c("Prop"))
"ei_03_5_1" <- subset(end_time_seq, Subject == "007_03" & trial_name == "ExpA_Set5_Item1", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_03" & end_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((ei_03_6_2 + ei_03_5_1) / 2)

#interpolate trial 9_3
"ei_03_7_4" <- subset(end_time_seq, Subject == "007_03" & trial_name == "ExpA_Set7_Item4", 
                    select = c("Prop"))
"ei_03_4_2" <- subset(end_time_seq, Subject == "007_03" & trial_name == "ExpA_Set4_Item2", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_03" & end_time_seq$trial_name == "ExpA_Set9_Item3"), "Prop"] <- ((ei_03_7_4 + ei_03_4_2) / 2)

#SUBJECT 007_04
#interpolate trial 9_4
"ei_04_11_2" <- subset(end_time_seq, Subject == "007_04" & trial_name == "ExpA_Set11_Item2", 
                    select = c("Prop"))
"ei_04_7_3" <- subset(end_time_seq, Subject == "007_04" & trial_name == "ExpA_Set7_Item3", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_04" & end_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((ei_04_11_2 + ei_04_7_3) / 2)

#interpolate trial 9_1
"ei_04_10_4" <- subset(end_time_seq, Subject == "007_04" & trial_name == "ExpA_Set10_Item4", 
                    select = c("Prop"))
"ei_04_1_2" <- subset(end_time_seq, Subject == "007_04" & trial_name == "ExpA_Set1_Item2", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_04" & end_time_seq$trial_name == "ExpA_Set9_Item1"), "Prop"] <- ((ei_04_10_4 + ei_04_1_2) / 2)

#interpolate trial 9_2
"ei_04_6_2" <- subset(end_time_seq, Subject == "007_04" & trial_name == "ExpA_Set6_Item2", 
                    select = c("Prop"))
"ei_04_5_1" <- subset(end_time_seq, Subject == "007_04" & trial_name == "ExpA_Set5_Item1", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_04" & end_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((ei_04_6_2 + ei_04_5_1) / 2)

#interpolate trial 9_3
"ei_04_1_3" <- subset(end_time_seq, Subject == "007_04" & trial_name == "ExpA_Set1_Item3", 
                    select = c("Prop"))
"ei_04_4_2" <- subset(end_time_seq, Subject == "007_04" & trial_name == "ExpA_Set4_Item2", 
                    select = c("Prop"))
end_time_seq[ which(end_time_seq$Subject == "007_04" & end_time_seq$trial_name == "ExpA_Set9_Item3"), "Prop"] <- ((ei_04_1_3 + ei_04_4_2) / 2)
```


```{r Save clean, interpolated verb time sequence}
#save the cleaned file to run in verb/object window analysis scripts. Or plot/analyse overall sentence below.
#saveRDS(end_time_seq, file = "007_end_data_clean_interpolated.rds")
saveRDS(end_time_seq, file = "007_end_data_clean_interpolated_alltrials.rds")
```
Saved 10/21

```{r load clean, interpolated verb time sequence}
#end_time_seq <- readRDS("007_end_data_clean_interpolated.rds", refhook = NULL)
#OR
end_time_seq <- readRDS("007_end_data_clean_interpolated_alltrials.rds", refhook = NULL)

#load required packages
library("Matrix")
library("plyr")
library("tidyr")
library("lme4")
library(lmerTest)
library("ggplot2")
library("eyetrackingR")
library("dplyr")
library(sjstats)
```

##final prep for modeling
```{r final prep for modeling}
end_time_seq$Group <- as.factor(end_time_seq$Group)
end_time_seq$AOI <- as.factor(end_time_seq$AOI)

#create column for natural quadratic time polynomial, if using
end_time_seq %>% 
  mutate(Time2 = Time^2) %>% 
#create column for binomial data instead of raw proportions
  mutate(BiProp = ifelse(Prop > .5, 1, 0)) -> end_time_seq

contrasts(end_time_seq$Group) <- c(-.5, .5)

#subset the data into target AOIs only
`action_end_time_seq` <- subset(end_time_seq, AOI == "ActionDis", 
                    select = c("Prop", "Subject", "Group", "ot1", "ot2", "trial_name", "Time", "BiProp"))

```

##model glmer

1) run null model to find variance components
```{r}
M0S.end <- glmer(BiProp ~ (1|Subject), data=action_end_time_seq, family = binomial, control = glmerControl(optimizer = "bobyqa"))
summary(M0S.end)
performance::icc(M0S.end)

M0T.end <- glmer(BiProp ~ (1|trial_name), data=action_end_time_seq, family = binomial, control = glmerControl(optimizer = "bobyqa"))
summary(M0T.end)
performance::icc(M0T.end)

M0.end <- glmer(BiProp ~ (1|Subject) + (1|trial_name), data=action_end_time_seq, family = binomial, control = glmerControl(optimizer = "bobyqa"))
summary(M0.end)
performance::icc(M0.end)
```
ICC (SUBJECT): 0.099. (.098 without anticipated trials)
ICC (TRIAL): 0.045  (.057 without anticipated trials)
ICC both components: .141  (.150 without anticipated trials)

##Then the full model, below. 
Use either random intercepts only or random slopes only depending on ANOVA results from prior step.
in this case, use random slopes for ot1 (p < .0001) and ot2 (p < .0001), based on correlated random effects model
```{r}

end.target.binom.model.alltrials <- glmer(BiProp ~ (ot1+ot2)*Group + (1+ot1+ot2 | Subject) + (1+ot1+ot2 | trial_name), data=action_end_time_seq, family = binomial, glmerControl(optimizer = c("bobyqa")))
summary(end.target.binom.model.alltrials)

saveRDS(end.target.binom.model.alltrials, "./end.target.binom.model.alltrials.rds")

#or:
# end.target.binom.model <- glmer(BiProp ~ (ot1+ot2)*Group + (1+ot1+ot2 | Subject) + (1+ot1+ot2 | trial_name), data=action_end_time_seq, family = binomial, glmerControl(optimizer = c("bobyqa")))
# summary(end.target.binom.model)
# 
# saveRDS(end.target.binom.model, "./end.target.binom.model.rds")

```

```{r reload models}
#re-load model(s)
end.target.binom.model <- readRDS("./end.target.binom.model.rds")
summary(end.target.binom.model)
#OR
end.target.binom.model.alltrials <- readRDS("./end.target.binom.model.alltrials.rds")
summary(end.target.binom.model.alltrials)
```