---
title: "007_LookTargetAnalysis_eR_3_3_20"
author: "Katharine Aveni"
date: "March 3, 2020"
output: html_document
---

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

```{r Load previously cleaned data}
look_data_clean <- readRDS("007_Look_data_clean_alltrials.rds", refhook = NULL) #or remove _alltrials, to run data without anticipated trials

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

## Make and Plot Time Sequence for Target  window only
###Target window
```{r}
look_target_window_clean <- subset_by_window(look_data_clean,
                                window_start_time = 391,
                                window_end_time = 991,
                                rezero = TRUE, remove = TRUE)
```
extended by 25ms to create equal bin size. shorten by 25, or lengthen into "the" instead?

```{r}
#create time sequence
look_target_time_seq <-make_time_sequence_data(look_target_window_clean, time_bin_size = 50,
                                   predictor_columns = c("Group", "trial_name", "list")
                                   )
```


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

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

# store vectors for: original trial info; props for each bin in surrounding 2 trials; averaged new trial.

#SUBJECT 007_05
#interpolate trial Q/17
`it_05_16` <- subset(look_target_time_seq, Subject == "007_05" & trial_name == "ExpB_Item 16", 
                    select = c("Prop"))
`it_05_15` <- subset(look_target_time_seq, Subject == "007_05" & trial_name == "ExpB_Item 15", 
                    select = c("Prop"))
look_target_time_seq[which(look_target_time_seq$Subject == "007_05" & look_target_time_seq$trial_name == "ExpB_Item 17"), "Prop"] <- ((`it_05_16` + `it_05_15`)/2)

#interpolate trial J/10
`it_05_8` <- subset(look_target_time_seq, Subject == "007_05" & trial_name == "ExpB_Item 8", 
                    select = c("Prop"))
`it_05_14` <- subset(look_target_time_seq, Subject == "007_05" & trial_name == "ExpB_Item 14", 
                    select = c("Prop"))
look_target_time_seq[which(look_target_time_seq$Subject == "007_05" & look_target_time_seq$trial_name == "ExpB_Item 10"), "Prop"] <- ((`it_05_8` + `it_05_14`)/2)


#SUBJECT 02
#interpolate trial J/10
`it_02_4` <- subset(look_target_time_seq, Subject == "007_02" & trial_name == "ExpB_Item 4", 
                    select = c("Prop"))
`it_02_16` <- subset(look_target_time_seq, Subject == "007_02" & trial_name == "ExpB_Item 16", 
                    select = c("Prop"))
look_target_time_seq[which(look_target_time_seq$Subject == "007_02" & look_target_time_seq$trial_name == "ExpB_Item 10"), "Prop"] <- ((`it_02_4` + `it_02_16`)/2)


#interpolate trial Q/17
`it_02_15` <- subset(look_target_time_seq, Subject == "007_02" & trial_name == "ExpB_Item 15", 
                    select = c("Prop"))
`it_02_5` <- subset(look_target_time_seq, Subject == "007_02" & trial_name == "ExpB_Item 5", 
                    select = c("Prop"))
look_target_time_seq[which(look_target_time_seq$Subject == "007_02" & look_target_time_seq$trial_name == "ExpB_Item 17"), "Prop"] <- ((`it_02_15` + `it_02_5`)/2)


#SUBJECT 06
#interpolate trial Q/17
`it_06_15` <- subset(look_target_time_seq, Subject == "007_06" & trial_name == "ExpB_Item 15", 
                    select = c("Prop"))
`it_06_11` <- subset(look_target_time_seq, Subject == "007_06" & trial_name == "ExpB_Item 11", 
                    select = c("Prop"))
look_target_time_seq[which(look_target_time_seq$Subject == "007_06" & look_target_time_seq$trial_name == "ExpB_Item 17"), "Prop"] <- ((`it_06_15` + `it_06_11`)/2)

#SUBJECT 04
#interpolate trial J/10
`it_04_12` <- subset(look_target_time_seq, Subject == "007_04" & trial_name == "ExpB_Item 12", 
                    select = c("Prop"))
`it_04_6` <- subset(look_target_time_seq, Subject == "007_04" & trial_name == "ExpB_Item 6", 
                    select = c("Prop"))
look_target_time_seq[which(look_target_time_seq$Subject == "007_04" & look_target_time_seq$trial_name == "ExpB_Item 10"), "Prop"] <- ((`it_04_12` + `it_04_6`)/2)
```


```{r Save interpolated data file}
#save the cleaned file to run in verb/object window analysis scripts. Or plot/analyse overall sentence below.
saveRDS(look_target_time_seq, file = "007_Look_target_data_clean_interpolated_alltrials.rds")
#OR
#saveRDS(look_target_time_seq, file = "007_Look_target_data_clean_interpolated.rds")

```

```{r Load previously cleaned and interpolated data}
#pick one of the below to load
#look_target_time_seq <- readRDS("007_Look_target_data_clean_interpolated.rds", refhook = NULL)
look_target_time_seq <- readRDS("007_Look_target_data_clean_interpolated_alltrials.rds", refhook = NULL)

look_target_time_seq$Group <- as.factor(look_target_time_seq$Group)
look_target_time_seq$AOI <- as.factor(look_target_time_seq$AOI)
look_target_time_seq$list <- as.factor(look_target_time_seq$list)

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

###final prep for modeling
```{r final prep for modeling}

#create column for natural quadratic time polynomial, if using
look_target_time_seq %>% 
  mutate(Time2 = Time^2) %>% 
#create column of proportions rounded to 0 or 1
  mutate(BiProp = ifelse(Prop > .5, 1, 0)) -> look_target_time_seq
  
#set group sum contrasts  
  contrasts(look_target_time_seq$Group) <- c(-.5, .5)
#center groupvariable
#  mutate(GroupCent= ifelse(Group.x== "PD", .5, -.5))-> look_target_time_seq

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

#create dataframe with average proportions from each subject, splits AOIs into different columns
look_target_time_seq %>%
  filter(!is.na(Prop)) %>%
  group_by(Subject, Time, Group) %>% 
  summarise(MeanProp = mean(Prop), y=sum(Prop), N=length(Prop), Time=min(Time)) -> look_target_time_seq_avg

#add new column that converts probabilities to logits
look_target_time_seq_avg$Logit <- qlogis(look_target_time_seq_avg$MeanProp)

```


#model glmer

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

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

M0.look <- glmer(BiProp ~ (1|Subject) + (1|trial_name), data=look_target_time_seq, family = binomial, control = glmerControl(optimizer = "bobyqa"))
summary(M0.look)
performance::icc(M0.look) #0.077

```

##Then the full model, below. 
```{r}

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

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

```{r plot target line vs model fit}
#reload model if needed
look.target.binom.model.alltrials <- readRDS("./look.target.binom.model.alltrials.rds")

new.target <- look_target_time_seq
new.target$predicted.target <- predict(look.target.binom.model.alltrials, newdata=look_target_time_seq)

tiff("baseline_plot_modelfit_12_2.tiff", units="in", width=6, height=4.5, res=600)
ggplot() +
  geom_smooth(data = look_target_time_seq_avg, fun.data = mean_se, geom="line", size = 2, aes(x = Time, y = Logit, color='black'), fill='darkgrey') + 
  #stat_summary(fun.data = mean_se, geom = "errorbar", aes(color='darkorchid')) +
  stat_summary(data = new.target, fun.data = mean_se, geom="line", size = 2, aes(x = Time, y = predicted.target, color='red')) +
  theme_minimal() +
  theme(legend.text=element_text(size=rel(1.2))) +
  theme(legend.position="none") +
  facet_grid(. ~ Group) +
  theme(legend.position = "bottom") +
  scale_color_manual(
    name=NULL,
    values=c('black', 'red'),
    labels = c('Data', 'Model prediction'))+
   scale_fill_manual(
    name = NULL,
    values=c('grey', 'grey'),
    labels = NULL) +
  scale_x_continuous(breaks = c(0, 100, 200, 300, 400, 500, 600), labels = c(0, 100, 200, 300, 400, 500, 600)) +
  labs(title="Baseline sentence target window: Model fit",
        x ="Time (ms)", y = "Log odds of fixations on target AOI")+
  guides(color=guide_legend(override.aes=list(fill=c('grey', 'white'))))
dev.off()
```

``` {r reload and print models}
look.target.binom.model <- readRDS("./look.target.binom.model.rds")
summary(look.target.binom.model)

look.target.binom.model.alltrials <- readRDS("./look.target.binom.model.alltrials.rds")
summary(look.target.binom.model.alltrials)

```