--- title: "007_AgentAnalysis_11_12" author: "Katharine Aveni" date: "November 12, 2021" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ```{r Load Libraries & Data} #pick ONE of these: #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") ``` ##Agent window ### zoom in to agent window (already offset by 200ms) extended earlier into "the" window by 24ms ```{r} agent_window <- subset_by_window(sentence_window_clean, window_start_time = 58, window_end_time = 658, rezero = TRUE, remove = TRUE) ``` ###create agent window time series ```{r Create verb window time series} agent_time_seq <-make_time_sequence_data(agent_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} agent_time_seq <- as.data.frame(agent_time_seq) #SUBJECT 007_05 #interpolate trial 9_3 `ai_05_3_2` <- subset(agent_time_seq, Subject == "007_05" & trial_name == "ExpA_Set3_Item2", select = c("Prop")) `ai_05_7_3` <- subset(agent_time_seq, Subject == "007_05" & trial_name == "ExpA_Set7_Item3", select = c("Prop")) agent_time_seq[which(agent_time_seq$Subject == "007_05" & agent_time_seq$trial_name == "ExpA_Set9_Item3"), "Prop"] <- ((`ai_05_3_2` + `ai_05_7_3`)/2) #interpolate trial 9_2 "ai_05_7_2" <- subset(agent_time_seq, Subject == "007_05" & trial_name == "ExpA_Set7_Item2", select = c("Prop")) "ai_05_12_1" <- subset(agent_time_seq, Subject == "007_05" & trial_name == "ExpA_Set12_Item1", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_05" & agent_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((ai_05_7_2 + ai_05_12_1) / 2) #interpolate trial 9_4 "ai_05_11_4" <- subset(agent_time_seq, Subject == "007_05" & trial_name == "ExpA_Set11_Item4", select = c("Prop")) "ai_05_12_2" <- subset(agent_time_seq, Subject == "007_05" & trial_name == "ExpA_Set12_Item2", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_05" & agent_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((ai_05_11_4 + ai_05_12_2) / 2) #SUBJECT 007_02 #interpolate trial 9_4 "ai_02_10_2" <- subset(agent_time_seq, Subject == "007_02" & trial_name == "ExpA_Set10_Item2", select = c("Prop")) "ai_02_3_2" <- subset(agent_time_seq, Subject == "007_02" & trial_name == "ExpA_Set3_Item2", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_02" & agent_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((ai_02_10_2 + ai_02_3_2) / 2) #interpolate trial 9_1 "ai_02_3_2" <- subset(agent_time_seq, Subject == "007_02" & trial_name == "ExpA_Set3_Item2", select = c("Prop")) "ai_02_5_3" <- subset(agent_time_seq, Subject == "007_02" & trial_name == "ExpA_Set5_Item3", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_02" & agent_time_seq$trial_name == "ExpA_Set9_Item1"), "Prop"] <- ((ai_02_3_2 + ai_02_5_3) / 2) #interpolate trial 9_2 "ai_02_3_4" <- subset(agent_time_seq, Subject == "007_02" & trial_name == "ExpA_Set3_Item4", select = c("Prop")) "ai_02_5_1" <- subset(agent_time_seq, Subject == "007_02" & trial_name == "ExpA_Set5_Item1", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_02" & agent_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((ai_02_3_4 + ai_02_5_1) / 2) #SUBJECT 007_06 #interpolate trial 9_4 "ai_06_2_1" <- subset(agent_time_seq, Subject == "007_06" & trial_name == "ExpA_Set2_Item1", select = c("Prop")) "ai_06_7_2" <- subset(agent_time_seq, Subject == "007_06" & trial_name == "ExpA_Set7_Item2", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_06" & agent_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((ai_06_2_1 + ai_06_7_2) / 2) #interpolate trial 9_2 "ai_06_7_4" <- subset(agent_time_seq, Subject == "007_06" & trial_name == "ExpA_Set7_Item4", select = c("Prop")) "ai_06_4_2" <- subset(agent_time_seq, Subject == "007_06" & trial_name == "ExpA_Set4_Item2", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_06" & agent_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((ai_06_7_4 + ai_06_4_2) / 2) #SUBJECT 007_03 #interpolate trial 9_4 "ai_03_11_2" <- subset(agent_time_seq, Subject == "007_03" & trial_name == "ExpA_Set11_Item2", select = c("Prop")) "ai_03_12_1" <- subset(agent_time_seq, Subject == "007_03" & trial_name == "ExpA_Set12_Item1", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_03" & agent_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((ai_03_11_2 + ai_03_12_1) / 2) #interpolate trial 9_2 "ai_03_6_2" <- subset(agent_time_seq, Subject == "007_03" & trial_name == "ExpA_Set6_Item2", select = c("Prop")) "ai_03_5_1" <- subset(agent_time_seq, Subject == "007_03" & trial_name == "ExpA_Set5_Item1", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_03" & agent_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((ai_03_6_2 + ai_03_5_1) / 2) #interpolate trial 9_3 "ai_03_7_4" <- subset(agent_time_seq, Subject == "007_03" & trial_name == "ExpA_Set7_Item4", select = c("Prop")) "ai_03_4_2" <- subset(agent_time_seq, Subject == "007_03" & trial_name == "ExpA_Set4_Item2", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_03" & agent_time_seq$trial_name == "ExpA_Set9_Item3"), "Prop"] <- ((ai_03_7_4 + ai_03_4_2) / 2) #SUBJECT 007_04 #interpolate trial 9_4 "ai_04_11_2" <- subset(agent_time_seq, Subject == "007_04" & trial_name == "ExpA_Set11_Item2", select = c("Prop")) "ai_04_7_3" <- subset(agent_time_seq, Subject == "007_04" & trial_name == "ExpA_Set7_Item3", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_04" & agent_time_seq$trial_name == "ExpA_Set9_Item4"), "Prop"] <- ((ai_04_11_2 + ai_04_7_3) / 2) #interpolate trial 9_1 "ai_04_10_4" <- subset(agent_time_seq, Subject == "007_04" & trial_name == "ExpA_Set10_Item4", select = c("Prop")) "ai_04_1_2" <- subset(agent_time_seq, Subject == "007_04" & trial_name == "ExpA_Set1_Item2", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_04" & agent_time_seq$trial_name == "ExpA_Set9_Item1"), "Prop"] <- ((ai_04_10_4 + ai_04_1_2) / 2) #interpolate trial 9_2 "ai_04_6_2" <- subset(agent_time_seq, Subject == "007_04" & trial_name == "ExpA_Set6_Item2", select = c("Prop")) "ai_04_5_1" <- subset(agent_time_seq, Subject == "007_04" & trial_name == "ExpA_Set5_Item1", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_04" & agent_time_seq$trial_name == "ExpA_Set9_Item2"), "Prop"] <- ((ai_04_6_2 + ai_04_5_1) / 2) #interpolate trial 9_3 "ai_04_1_3" <- subset(agent_time_seq, Subject == "007_04" & trial_name == "ExpA_Set1_Item3", select = c("Prop")) "ai_04_4_2" <- subset(agent_time_seq, Subject == "007_04" & trial_name == "ExpA_Set4_Item2", select = c("Prop")) agent_time_seq[ which(agent_time_seq$Subject == "007_04" & agent_time_seq$trial_name == "ExpA_Set9_Item3"), "Prop"] <- ((ai_04_1_3 + ai_04_4_2) / 2) ``` ```{r Save clean, interpolated verb time sequence} #save the cleaned file to run in verb/agent window analysis scripts. Or plot/analyse overall sentence below. #saveRDS(agent_time_seq, file = "007_agent_data_clean_interpolated.rds") saveRDS(agent_time_seq, file = "007_agent_data_interpolated_with_anticipated.rds") ``` Saved 10/1/21 ```{r load clean, interpolated verb time sequence} #agent_time_seq <- readRDS("007_agent_data_clean_interpolated.rds", refhook = NULL) agent_time_seq <- readRDS("007_agent_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") #library(sjstats) #library(corpora) ``` ###Final prep for modeling ```{r final prep for modeling} agent_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")) -> agent_time_seq agent_time_seq$mt_content <- as.factor(agent_time_seq$mt_content) agent_time_seq$Group <- as.factor(agent_time_seq$Group) #create column for natural quadratic time polynomial, if using agent_time_seq %>% mutate(Time2 = Time^2) %>% #create column for binomial data instead of raw proportions mutate(BiProp = ifelse(Prop > .5, 1, 0)) -> agent_time_seq #center factors. mt: low -.5, hi .5 group: control -.5, pd .5 contrasts(agent_time_seq$mt_content) <- c(.5, -.5) contrasts(agent_time_seq$Group) <- c(-.5, .5) #create dataframe that excludes set 9 trials for early participants agent_time_seq_9rm <- agent_time_seq[!((agent_time_seq$Subject=="007_02" | agent_time_seq$Subject=="007_03" | agent_time_seq$Subject=="007_04" | agent_time_seq$Subject=="007_05" | agent_time_seq$Subject=="007_06") & (agent_time_seq$trial_name=="ExpA_Set9_Item1" | agent_time_seq$trial_name=="ExpA_Set9_Item2" | agent_time_seq$trial_name=="ExpA_Set9_Item3" | agent_time_seq$trial_name=="ExpA_Set9_Item4")),] #subset the data into target AOIs only `target_agent_time_seq` <- subset(agent_time_seq, AOI == "Target", select = c("Prop", "Subject", "Group", "ot1", "ot2", "trial_name", "Time", "Time2", "BiProp", "mt_content")) `target_agent_time_seq_9rm` <- subset(agent_time_seq_9rm, AOI == "Target", select = c("Prop", "Subject", "Group", "ot1", "ot2", "trial_name", "Time", "Time2", "BiProp", "mt_content")) #create dataframe with average proportions from each subject, splits AOIs into different columns target_agent_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_agent_time_seq_avg #add new column that converts probabilities to logits target_agent_time_seq_avg$Logit <- qlogis(target_agent_time_seq_avg$MeanProp) ``` #model glmer run null model to find variance components ```{r} M0S.agent <- glmer(BiProp ~ (1|Subject), data=target_agent_time_seq, family = binomial, control = glmerControl(optimizer = "bobyqa")) summary(M0S.agent) performance::icc(M0S.agent) M0T.agent <- glmer(BiProp ~ (1|trial_name), data=target_agent_time_seq, family = binomial, control = glmerControl(optimizer = "bobyqa")) summary(M0T.agent) performance::icc(M0T.agent) M0.agent <- glmer(BiProp ~ (1|Subject) + (1|trial_name), data=target_agent_time_seq, family = binomial, control = glmerControl(optimizer = "bobyqa")) summary(M0.agent) performance::icc(M0.agent) ``` ICC (SUBJECT): 0.030 ICC (TRIAL): .056 ICC total: 0.088 ##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 both ```{r} #RENAME AS NEEDED '.ALLTRIALS' OR '__'. model code. group*time1 &2 interactions, etc. agent.target.binom.model.alltrials <- glmer(BiProp ~ (ot1+ot2)*Group + (1+ot1| Subject) + (1+ot1| trial_name), data=target_agent_time_seq, family = binomial, glmerControl(optimizer = c("bobyqa"))) summary(agent.target.binom.model.alltrials) #save model results saveRDS(agent.target.binom.model.alltrials, "./agent.target.binom.model.alltrials.rds") # re-load the model agent.target.binom.model.alltrials <- readRDS("./agent.target.binom.model.alltrials.rds") print(summary(agent.target.binom.model.alltrials)) #OR anticipated trials out, in which case quadratic random effects in due to different target curve shape # agent.target.binom.model <- glmer(BiProp ~ (ot1+ot2)*Group + (1+ot1+ot2| Subject) + (1+ot1+ot2| trial_name), data=target_agent_time_seq, family = binomial, glmerControl(optimizer = c("bobyqa"))) # summary(agent.target.binom.model) #save model results # saveRDS(agent.target.binom.model, "./agent.target.binom.model.rds") # re-load the model # agent.target.binom.model <- readRDS("./agent.target.binom.model.rds") # print(summary(agent.target.binom.model)) ``` ```{r motion content model} MT.agent.target.binom.model.alltrials.9rm <- glmer(BiProp ~ (ot1+ot2) * Group * mt_content + (1+ot1| Subject) + (1+ot1| trial_name), data=target_agent_time_seq_9rm, family = binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000))) #summary(MT.agent.target.binom.model.alltrials) summary(MT.agent.target.binom.model.alltrials.9rm) #note: same pattern of significance in smaller set of trials (9rm) #saveRDS(MT.agent.target.binom.model.alltrials, "./MT.agent.target.binom.model.alltrials.rds") saveRDS(MT.agent.target.binom.model.alltrials.9rm, "./MT.agent.target.binom.model.alltrials.9rm.rds") #reload MT.agent.target.binom.model.alltrials.9rm <- readRDS("./MT.agent.target.binom.model.alltrials.9rm.rds") print(summary(MT.agent.target.binom.model.alltrials.9rm)) ``` ```{r reload motion content models} #re-load model(s) MT.agent.target.binom.model <- readRDS("./MT.agent.target.binom.model.rds") print(summary(MT.agent.target.binom.model)) #OR MT.agent.target.binom.model.9rm <- readRDS("./MT.agent.target.binom.model.9rm.rds") print(summary(MT.agent.target.binom.model.9rm)) #OR MT.agent.target.binom.model.alltrials <- readRDS("./MT.agent.target.binom.model.alltrials.rds") print(summary(MT.agent.target.binom.model.alltrials)) #OR MT.agent.target.binom.model.alltrials.9rm <- readRDS("./MT.agent.target.binom.model.alltrials.9rm.rds") print(summary(MT.agent.target.binom.model.alltrials.9rm)) ``` ###test for sig. above chance level looks to target and agent distractor ```{r} #test <25% chance looks to target in INITIAL time bin-- pre-anticipation effect? of interest only if running analyses with anticipated trials in target_agent_time_seq_end <- target_agent_time_seq[which(target_agent_time_seq$Time==0),] t.test(target_agent_time_seq_end$BiProp, mu=.25) #test <25% chance looks to target in END time bin-- to show sig. linear increase meaningful target_agent_time_seq_end <- target_agent_time_seq[which(target_agent_time_seq$Time==550),] table(target_agent_time_seq_end$BiProp) prop.test(565, 1565, p = .25, alternative = c("two.sided"), conf.level = 0.95, correct = FALSE) #test <25% chance looks to agent-distractor in initial time bin agent_agent_time_seq <- subset(agent_time_seq, AOI == "AgentDis", select = c("Prop", "Subject", "Group", "ot1", "ot2", "trial_name", "Time", "Time2", "BiProp")) agent_agent_time_seq_end <- agent_agent_time_seq[which(agent_agent_time_seq$Time==0),] t.test(agent_agent_time_seq_end$BiProp, mu=.25) ``` ###ggplot agent window ```{r} #change order of levels in AOI factor in order to change plot default levels(as.factor(agent_time_seq$AOI)) agent_time_seq$AOI <- factor(agent_time_seq$AOI, levels = c("Target", "AgentDis", "ActionDis", "Unrelated")) #tiff("Z:\\LCAN Internal\\LCAN_007_LangPrediction_QRP\\Analyses\\ET_Data\\Results\\Agent_plot_10_24.tiff", units="in", width=7.5, height=5, res=300) ggplot(agent_time_seq_9rm, aes(x = Time, y = BiProp, col= AOI)) + #geom_smooth() + #stat_summary(fun.data = mean_se, geom = "errorbar", aes(color=paste("mean_se", AOI))) + stat_summary(geom="line", size = 2) + theme_minimal() + theme(legend.text=element_text(size=rel(1.2))) + theme(legend.position="bottom", legend.box = "horizontal") + facet_grid(Group ~ .) + #scale_color_manual(labels = c("Target", "AgentDistractor", "ActionDistractor", "Unrelated"), values=c('darkorchid2','turquoise3', 'red2', 'chartreuse3')) + labs(title="Predictive Sentences-- Agent Window", x ="Time (ms)", y = "Proportion of Looks") #dev.off() ``` ```{r target plot data vs model} #reload model if needed agent.target.binom.model.alltrials <- readRDS("./agent.target.binom.model.alltrials.rds") new.agent <- target_agent_time_seq new.agent$predicted.agent <- predict(agent.target.binom.model.alltrials, newdata=target_agent_time_seq) tiff("agent_plot_modelfit_12_2.tiff", units="in", width=3, height=5, res=600) ggplot() + geom_smooth(data = target_agent_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.agent, fun.data = mean_se, geom="line", size = 1, aes(x = Time, y = predicted.agent, 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), labels = c(0, 100, 200, 300, 400, 500))+ labs(title="Agent 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() ``` ```{r agent window mean looks to target by group} table(target_agent_time_seq_9rm$Group, target_agent_time_seq_9rm$mt_content, target_agent_time_seq_9rm$BiProp) #control looks to target: 29.5% in hi, 29.7% in low #PD looks to target: 25.23% in hi, 33.28% in low ```