---
title: "007_ObjAnalysis_3_3_20"
author: "Katharine Aveni"
date: "March 3, 2020"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
##Load Libraries & Data
```{r}
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")
```
##Object window
### zoom in to object window (offset by 200ms)
```{r}
obj_window <- subset_by_window(sentence_window_clean,
window_start_time = 1246,
window_end_time = 1946,
rezero = TRUE, remove = TRUE)
```
Reduced object window by 4 ms in order to avoid final small bin (at this endpoint, stimulus would have ended 196 ms ago).
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}
obj_time_seq <-make_time_sequence_data(obj_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.
007_02: list 9; 007_03 & 007_04: list 10; 007_05: list 1; 007_06: list 2
Each chunk stores vectors for: Proportions for each bin in surrounding 2 trials (unless surrounding trial was a question trial or pre-anticipated trial, in which case the next closest trial was stored).
Stored vectors were averaged and written to new trial.
```{r Create vectors for interpolated trials}
obj_time_seq <- as.data.frame(obj_time_seq)
#SUBJECT 007_05
#interpolate trial 9_3
`oi_05_3_2` <- subset(obj_time_seq, Subject == "007_05" & trial_name == "ExpA_Set3_Item2",
select = c("Prop"))
`oi_05_7_3` <- subset(obj_time_seq, Subject == "007_05" & trial_name == "ExpA_Set7_Item3",
select = c("Prop"))
obj_time_seq[which(obj_time_seq$Subject == "007_05" & obj_time_seq$trial_name == "ExpA_Set9_Item3"), "Prop"] <- ((`oi_05_3_2` + `oi_05_7_3`)/2)
#interpolate trial 9_2
"oi_05_7_2" <- subset(obj_time_seq, Subject == "007_05" & trial_name == "ExpA_Set7_Item2",
select = c("Prop"))
"oi_05_12_1" <- subset(obj_time_seq, Subject == "007_05" & trial_name == "ExpA_Set12_Item1",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_05" & obj_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((oi_05_7_2 + oi_05_12_1) / 2)
#interpolate trial 9_4
"oi_05_11_4" <- subset(obj_time_seq, Subject == "007_05" & trial_name == "ExpA_Set11_Item4",
select = c("Prop"))
"oi_05_12_2" <- subset(obj_time_seq, Subject == "007_05" & trial_name == "ExpA_Set12_Item2",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_05" & obj_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((oi_05_11_4 + oi_05_12_2) / 2)
#SUBJECT 007_02
#interpolate trial 9_4
"oi_02_10_2" <- subset(obj_time_seq, Subject == "007_02" & trial_name == "ExpA_Set10_Item2",
select = c("Prop"))
"oi_02_3_2" <- subset(obj_time_seq, Subject == "007_02" & trial_name == "ExpA_Set3_Item2",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_02" & obj_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((oi_02_10_2 + oi_02_3_2) / 2)
#interpolate trial 9_1
"oi_02_3_2" <- subset(obj_time_seq, Subject == "007_02" & trial_name == "ExpA_Set3_Item2",
select = c("Prop"))
"oi_02_5_3" <- subset(obj_time_seq, Subject == "007_02" & trial_name == "ExpA_Set5_Item3",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_02" & obj_time_seq$trial_name == "ExpA_Set9_Item1"), "Prop"] <- ((oi_02_3_2 + oi_02_5_3) / 2)
#interpolate trial 9_2
"oi_02_3_4" <- subset(obj_time_seq, Subject == "007_02" & trial_name == "ExpA_Set3_Item4",
select = c("Prop"))
"oi_02_5_1" <- subset(obj_time_seq, Subject == "007_02" & trial_name == "ExpA_Set5_Item1",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_02" & obj_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((oi_02_3_4 + oi_02_5_1) / 2)
#SUBJECT 007_06
#interpolate trial 9_4
"oi_06_2_1" <- subset(obj_time_seq, Subject == "007_06" & trial_name == "ExpA_Set2_Item1",
select = c("Prop"))
"oi_06_7_2" <- subset(obj_time_seq, Subject == "007_06" & trial_name == "ExpA_Set7_Item2",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_06" & obj_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((oi_06_2_1 + oi_06_7_2) / 2)
#interpolate trial 9_2
"oi_06_7_4" <- subset(obj_time_seq, Subject == "007_06" & trial_name == "ExpA_Set7_Item4",
select = c("Prop"))
"oi_06_4_2" <- subset(obj_time_seq, Subject == "007_06" & trial_name == "ExpA_Set4_Item2",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_06" & obj_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((oi_06_7_4 + oi_06_4_2) / 2)
#SUBJECT 007_03
#interpolate trial 9_4
"oi_03_11_2" <- subset(obj_time_seq, Subject == "007_03" & trial_name == "ExpA_Set11_Item2",
select = c("Prop"))
"oi_03_12_1" <- subset(obj_time_seq, Subject == "007_03" & trial_name == "ExpA_Set12_Item1",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_03" & obj_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((oi_03_11_2 + oi_03_12_1) / 2)
#interpolate trial 9_2
"oi_03_6_2" <- subset(obj_time_seq, Subject == "007_03" & trial_name == "ExpA_Set6_Item2",
select = c("Prop"))
"oi_03_5_1" <- subset(obj_time_seq, Subject == "007_03" & trial_name == "ExpA_Set5_Item1",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_03" & obj_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((oi_03_6_2 + oi_03_5_1) / 2)
#interpolate trial 9_3
"oi_03_7_4" <- subset(obj_time_seq, Subject == "007_03" & trial_name == "ExpA_Set7_Item4",
select = c("Prop"))
"oi_03_4_2" <- subset(obj_time_seq, Subject == "007_03" & trial_name == "ExpA_Set4_Item2",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_03" & obj_time_seq$trial_name == "ExpA_Set9_Item3"), "Prop"] <- ((oi_03_7_4 + oi_03_4_2) / 2)
#SUBJECT 007_04
#interpolate trial 9_4
"oi_04_11_2" <- subset(obj_time_seq, Subject == "007_04" & trial_name == "ExpA_Set11_Item2",
select = c("Prop"))
"oi_04_7_3" <- subset(obj_time_seq, Subject == "007_04" & trial_name == "ExpA_Set7_Item3",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_04" & obj_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((oi_04_11_2 + oi_04_7_3) / 2)
#interpolate trial 9_1
"oi_04_10_4" <- subset(obj_time_seq, Subject == "007_04" & trial_name == "ExpA_Set10_Item4",
select = c("Prop"))
"oi_04_1_2" <- subset(obj_time_seq, Subject == "007_04" & trial_name == "ExpA_Set1_Item2",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_04" & obj_time_seq$trial_name == "ExpA_Set9_Item1"), "Prop"] <- ((oi_04_10_4 + oi_04_1_2) / 2)
#interpolate trial 9_2
"oi_04_6_2" <- subset(obj_time_seq, Subject == "007_04" & trial_name == "ExpA_Set6_Item2",
select = c("Prop"))
"oi_04_5_1" <- subset(obj_time_seq, Subject == "007_04" & trial_name == "ExpA_Set5_Item1",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_04" & obj_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((oi_04_6_2 + oi_04_5_1) / 2)
#interpolate trial 9_3
"oi_04_1_3" <- subset(obj_time_seq, Subject == "007_04" & trial_name == "ExpA_Set1_Item3",
select = c("Prop"))
"oi_04_4_2" <- subset(obj_time_seq, Subject == "007_04" & trial_name == "ExpA_Set4_Item2",
select = c("Prop"))
obj_time_seq[ which(obj_time_seq$Subject == "007_04" & obj_time_seq$trial_name == "ExpA_Set9_Item3"), "Prop"] <- ((oi_04_1_3 + oi_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(obj_time_seq, file = "007_obj_data_interpolated_with_anticipated.rds")
```
```{r load clean, interpolated verb time sequence}
#pick ONE to read in
#obj_time_seq <- readRDS("007_obj_data_clean_interpolated.rds", refhook = NULL)
obj_time_seq <- readRDS("007_obj_data_interpolated_with_anticipated.rds", refhook = NULL)
#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}
#label high vs. low action verbs
obj_time_seq %>%
mutate(mt_content = ifelse(file == "boxer_hits_punchingbag" | file == "driver_hits_brakes" |
file == "cat_catches_mouse" | file == " child_catches_football" |
file == "dolphin_jumps_waves" | file == "horse_jumps_fence" |
file == "lifeguard_saves_swimmer" | file == "detective_hunts_fugitive" |
file == "lion_hunts_zebra" | file == "lion_crosses_jungle" |
file == "detective_crosses_courtroom" | file == "mechanic_organizes_tools" |
file == "waitress_organizes_dishes", "Hi", "Lo")) -> obj_time_seq
obj_time_seq$mt_content <- as.factor(obj_time_seq$mt_content)
obj_time_seq$Group <- as.factor(obj_time_seq$Group)
#create column for natural quadratic time polynomial, if using
obj_time_seq %>%
mutate(Time2 = Time^2) %>%
#create column for binomial data instead of raw proportions
mutate(BiProp = ifelse(Prop > .5, 1, 0)) -> obj_time_seq
#center factors. mt: low -.5, hi .5 group: control -.5, pd .5
contrasts(obj_time_seq$mt_content) <- c(.5, -.5)
contrasts(obj_time_seq$Group) <- c(-.5, .5)
#create subset with set 9 trials removed for early Ss
obj_time_seq_9rm <- obj_time_seq[!((obj_time_seq$Subject=="007_02" | obj_time_seq$Subject=="007_03" | obj_time_seq$Subject=="007_04" | obj_time_seq$Subject=="007_05" | obj_time_seq$Subject=="007_06") & (obj_time_seq$trial_name=="ExpA_Set9_Item1" | obj_time_seq$trial_name=="ExpA_Set9_Item2" | obj_time_seq$trial_name=="ExpA_Set9_Item3" | obj_time_seq$trial_name=="ExpA_Set9_Item4")),]
#subset the data into target AOIs only
`target_obj_time_seq` <- subset(obj_time_seq, AOI == "Target",
select = c("Prop", "Subject", "Group", "ot1", "ot2", "trial_name", "Time", "BiProp", "mt_content"))
#target AOI only; set 9 removed for early Ss
`target_obj_time_seq_9rm` <- subset(obj_time_seq_9rm, AOI == "Target",
select = c("Prop", "Subject", "Group", "ot1", "ot2", "trial_name", "Time", "Time2", "BiProp", "mt_content"))
#action-based distractor AOI only (all trial sets)
`verb_dis_obj_time_seq` <- subset(obj_time_seq, AOI == "ActionDis",
select = c("Prop", "Subject", "Group", "ot1", "ot2", "trial_name", "Time", "BiProp", "mt_content"))
#PD only for UPDRS modeling
`PD_target_obj_time_seq` <- subset(obj_time_seq, (AOI == "Target" & Group == "PD"),
select = c("Prop", "Subject", "Group", "ot1", "ot2", "trial_name", "Time", "BiProp", "mt_content"))
#create dataframe with average proportions from each subject, splits AOIs into different columns
target_obj_time_seq %>%
filter(!is.na(Prop)) %>%
group_by(Subject, Time, Group) %>%
summarise(MeanProp = mean(Prop), y=sum(Prop), N=length(Prop), Time=min(Time)) -> target_obj_time_seq_avg
#add new column that converts probabilities to logits
target_obj_time_seq_avg$Logit <- qlogis(target_obj_time_seq_avg$MeanProp)
```
#model glmer
1) run null model to find variance components
```{r}
M0S.obj <- glmer(BiProp ~ (1|Subject), data=target_obj_time_seq, family = binomial, control = glmerControl(optimizer = "bobyqa"))
summary(M0S.obj)
performance::icc(M0S.obj)
M0T.obj <- glmer(BiProp ~ (1|trial_name), data=target_obj_time_seq, family = binomial, control = glmerControl(optimizer = "bobyqa"))
summary(M0T.obj)
performance::icc(M0T.obj)
M0.obj <- glmer(BiProp ~ (1|Subject) + (1|trial_name), data=target_obj_time_seq, family = binomial, control = glmerControl(optimizer = "bobyqa"))
summary(M0.obj)
performance::icc(M0.obj)
```
ICC (SUBJECT): 0.205
ICC (TRIAL): 0.048
ICC both components: .256
##Then the full model, below.
in this case, use random slopes for ot1 (p < .0001) and ot2 (p = 0.002151), based on correlated random effects model
```{r}
#including anticipated trials
obj.target.binom.model.alltrials <- glmer(BiProp ~ (ot1+ot2)*Group + (1+ot1+ot2 | Subject) + (1+ot1+ot2 | trial_name), data=target_obj_time_seq, family = binomial, glmerControl(optimizer = c("bobyqa")))
summary(obj.target.binom.model.alltrials)
saveRDS(obj.target.binom.model.alltrials, "./obj.target.binom.model.alltrials.rds")
#reload
obj.target.binom.model.alltrials <- readRDS("./obj.target.binom.model.alltrials.rds")
print(summary(obj.target.binom.model.alltrials))
#only unanticipated trials
# obj.target.binom.model <- glmer(BiProp ~ (ot1+ot2)*Group + (1+ot1+ot2 | Subject) + (1+ot1+ot2 | trial_name), data=target_obj_time_seq, family = binomial, glmerControl(optimizer = c("bobyqa")))
# summary(obj.target.binom.model)
#
# saveRDS(obj.target.binom.model, "./obj.target.binom.model.rds")
```
##Then the motion content model
```{r}
MT.obj.target.binom.model.alltrials.9rm <- glmer(BiProp ~ (ot1+ot2) * Group * mt_content + (1+ot1+ot2 | Subject) + (1+ot1+ot2 | trial_name), data=target_obj_time_seq_9rm, family = binomial, control=glmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=100000)))
summary(MT.obj.target.binom.model.alltrials.9rm)
saveRDS(MT.obj.target.binom.model.alltrials.9rm, "./MT.obj.target.binom.model.alltrials.9rm.rds")
#reload
MT.obj.target.binom.model.alltrials.9rm <- readRDS("./MT.obj.target.binom.model.alltrials.9rm.rds")
print(summary(MT.obj.target.binom.model.alltrials.9rm))
```
```{r reload motion content models}
#re-load model(s)
MT.obj.target.binom.model <- readRDS("./MT.obj.target.binom.model.rds")
print(summary(MT.obj.target.binom.model))
#OR
MT.obj.target.binom.model.9rm <- readRDS("./MT.obj.target.binom.model.9rm.rds")
print(summary(MT.obj.target.binom.model.9rm))
#OR
MT.obj.target.binom.model.alltrials <- readRDS("./MT.obj.target.binom.model.alltrials.rds")
print(summary(MT.obj.target.binom.model.alltrials))
#OR
MT.obj.target.binom.model.alltrials.9rm <- readRDS("./MT.obj.target.binom.model.alltrials.9rm.rds")
print(summary(MT.obj.target.binom.model.alltrials.9rm))
```
##plot
```{r Visualize obj window results}
#change order of levels in AOI factor in order to change plot default
levels(as.factor(obj_time_seq$AOI))
obj_time_seq$AOI <- factor(obj_time_seq$AOI, levels = c("Target", "AgentDis", "ActionDis", "Unrelated"))
obj_time_seq$mt_content <- factor(obj_time_seq$mt_content, levels = c("Lo", "Hi"))
#tiff("obj_plot_motion_content_with_anticipated_2_24_21.tiff", units="in", width=7.5, height=5, res=300)
ggplot(obj_time_seq, aes(x = Time, y = BiProp, col= AOI)) +
geom_smooth() +
theme_minimal() +
theme(legend.text=element_text(size=rel(1.2))) +
theme(legend.position="bottom", legend.box = "horizontal") +
facet_grid(Group ~ mt_content) +
scale_color_manual(labels = c("Target", "AgentDistractor", "ActionDistractor", "Unrelated"), values=c('darkorchid2','turquoise3', 'red2', 'chartreuse3')) +
labs(title="Predictive Sentences-- Target Window",
x ="Time (ms)", y = "Proportion of Looks")
#dev.off()
```
```{r plot target line vs model fit}
#reload model if needed
obj.target.binom.model.alltrials <- readRDS("./obj.target.binom.model.alltrials.rds")
new.obj <- target_obj_time_seq
new.obj$predicted.obj <- predict(obj.target.binom.model.alltrials, newdata=target_obj_time_seq)
tiff("obj_plot_modelfit_12_2.tiff", units="in", width=3, height=5, res=600)
ggplot() +
geom_smooth(data = target_obj_time_seq_avg, fun.data = mean_se, geom="line", alpha = 0.15, size = 1, aes(x = Time, y = Logit, color='darkorchid'), fill= 'darkorchid') +
stat_summary(data = new.obj, fun.data = mean_se, geom="line", size = 1, aes(x = Time, y = predicted.obj, color='black')) +
facet_grid(Group ~ .) +
theme_minimal() +
theme(legend.text=element_text(size=rel(1.2))) +
theme(legend.position = "bottom") +
scale_color_manual(
name=NULL,
values=c('black', 'darkorchid2'),
labels = c('Model prediction', 'Data')) +
scale_fill_manual(
name = NULL,
values=c('darkorchid', 'darkorchid'),
labels = NULL) +
expand_limits(y = c(-1.2, 1.2)) +
scale_x_continuous(breaks = c(0, 100, 200, 300, 400, 500, 600), labels = c(0, 100, 200, 300, 400, 500, 600)) +
labs(title="Target window: Model fit",
x ="Time (ms)", y = "Log odds of fixations on target AOI") +
guides(color=guide_legend(override.aes=list(fill=c('white', 'darkorchid'))))
dev.off()
```