--- 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() ```