library(knitr) opts_chunk$set(comment="", autodep = TRUE, cache=FALSE, echo=FALSE, warning=FALSE, message=FALSE)
# Kodi's libraries library(dplyr) # uses pre-0.5 version library(magrittr) library(stringr) library(ggplot2) library(ggthemes) require(lme4) library(lmerTest) library(tidyr) library(grid) library(gridExtra) library(mgcv) library(itsadug) library(pupilr) library(zplyr) source("/Users/zburchill/clarkegarrett2004/Exp1/Analysis/helper_functions/create_lmer_table.R")
# Given a ggplot object, a dataframe, a list of vectors of conditions, and a list of names, builds a plot build_facet_plot <-function(p,df,l_of_v,l_of_titles,rect_df=FALSE) { combined_ls<-purrr::map2(l_of_v,l_of_titles,list) v <- purrr::reduce(combined_ls, function(a,b) { if (rect_df==FALSE) { t <- a } else { t <- a + geom_rect(data = transform(rect_df,Comparison=b[[2]]), aes(xmin = xstart_rect, xmax = xend_rect, ymin = -Inf, ymax = Inf), fill = rect_df$colors,inherit.aes = FALSE) } t <- t + stat_summary(fun.data="mean_cl_boot", geom="errorbar", width=0.25, data=transform(subset(df,Condition %in% b[[1]]), Comparison=b[[2]])) t <- t + stat_summary(fun.y=mean,geom="point", data=transform(subset(df,Condition %in% b[[1]]), Comparison=b[[2]])) t <- t + stat_summary(fun.y=mean,geom="line", data=transform(subset(df,Condition %in% b[[1]]), Comparison=b[[2]])) return(t)}, .init=p) v <- v + facet_wrap(~Comparison) return(v) } ## ------------------------------------------ # Change font color with proper rendering across latex and html outputs # credit to Nicholas Hamilton (http://stackoverflow.com/questions/29067541/rmarkdown-how-to-change-the-font-color) font_color = function(x, color){ outputFormat = opts_knit$get("rmarkdown.pandoc.to") if(outputFormat == 'latex') paste("\\textcolor{",color,"}{",x,"}",sep="") else if(outputFormat == 'html') paste("<font color='",color,"'>",x,"</font>",sep="") else x }
#color scheme for accented stuff (dark to light): #e6550d #fdae6b #fee6ce #color scheme for unaccented stuff (''): #3182bd #9ecae1 #deebf7 color_values <- c("accented" = "#e6550d", "unaccented" = "#3182bd", "accented new" = "#fdae6b","unaccented new" = "#9ecae1")
d.all <- d.ex2# %>% filter(Run!="Exp2v2") d.all$ExpRun <- ifelse(d.all$Run=="Exp2v2", "Run2","Run1") d.noexp1 <- d.ex1 # print("Does the data pass the sanity check?") # sanity_check_data(d.all,pilot=FALSE) # sanity_check_data(d.noexp1,pilot=FALSE) d.all$PartOfExp <- factor(d.all$PartOfExp,levels=c("practice","main","baseline")) d.noexp1$PartOfExp <- factor(d.noexp1$PartOfExp,levels=c("practice","main","baseline")) # Load data for wording test d.wording <- d.ex1_wording_test d.wording$PartOfExp <- factor(d.wording$PartOfExp,levels=c("practice","main","newtest","baseline")) #d.wording$PartOfExp <- ifelse(d.wording$PartOfExp == "newtest","baseline",d.wording$PartOfExp) #d.wording$Condition <- "New" #d.wording$LoadTime <- NULL # Load data for the exposure blocks d.exposure_trials <- d.ex2.exposure d.catch_trials <- d.exposure_trials %>% group_by(WorkerId) %>% summarise(Misses=sum(ifelse(isCatch==TRUE & pressedSpace==FALSE,1,0)), FalsePositives=sum(ifelse(isCatch==FALSE & pressedSpace==TRUE,1,0)), AvgCorrectRT=mean(RT[isCatch==TRUE & pressedSpace==TRUE])) # Load in language background free response questions d.LgFreeNoExp <- read.csv("/Users/zburchill/clarkegarrett2004/Exp1/Analysis/no_exposure_pupillometry_language_free_response.csv") d.LgFreeExp <- read.csv("/Users/zburchill/clarkegarrett2004/Exp2/Analysis/exposure_pupil_lgbackground_free_response.csv") %>% select(WorkerId,LgRating) suppressWarnings( d.all <- left_join(d.all,d.LgFreeExp,by="WorkerId") ) suppressWarnings( d.noexp1 <- left_join(d.noexp1, d.LgFreeNoExp,by="WorkerId") ) suppressWarnings( d.wording <-left_join(d.wording,d.LgFreeNoExp,by="WorkerId") ) # calculate subject-wise mean RT for each Block d.all <- initialize_scores(d.all) d.noexp1 <- initialize_scores(d.noexp1) d.wording <- initialize_scores(d.wording) # Prepare some new survey questions d.all <- d.all %>% mutate_each( funs(factor(., levels=c("","dont_remember","never","hardly","sometimes","weekly","daily_plus"))), matches("last_week|last_month|childhood"))
duplicate_turkers <- d.all %>% select(WorkerId,UniqueID) %>% rbind(d.noexp1 %>% select(WorkerId,UniqueID), d.wording %>% select(WorkerId,UniqueID)) %>% distinct(WorkerId,UniqueID,.keep_all=TRUE) %>% group_by(WorkerId) %>% tally() %>% filter(n > 1) if (nrow(duplicate_turkers) > 0) { print("WARNING! DUPLICATE TURKERS! DO SOMETHING ABOUT THIS BETTER THAN WHAT WE HAVE RIGHT NOW!") } ## ---------------------------------------- # collect all non eligible subjets into a single df non_eligible_subjs_all <- pupilr::non_eligible_participants_f( d.all, duplicate_turkers = duplicate_turkers, newsurvey = TRUE) non_eligible_subjs_word <- pupilr::non_eligible_participants_f( d.wording, duplicate_turkers = duplicate_turkers, newsurvey = FALSE) non_eligible_subjs_noexp <- pupilr::non_eligible_participants_f( d.noexp1, duplicate_turkers = duplicate_turkers, newsurvey = FALSE) ## ---------------------------------------- # remove those inelligible subjects dat_out1 <- d.all %>% filter( !(WorkerId %in% non_eligible_subjs_all$WorkerId) ) dat_noexp_out1 <- d.noexp1 %>% filter( !(WorkerId %in% non_eligible_subjs_noexp$WorkerId) ) d.wording_out <- d.wording %>% filter( !(WorkerId %in% non_eligible_subjs_word$WorkerId) )
dat_out2 <- pupilr::exclude_extreme_rts_f(dat_out1) d.wording_out2 <- pupilr::exclude_extreme_rts_f(d.wording_out) dat_noexp_out2 <- pupilr::exclude_extreme_rts_f(dat_noexp_out1)
# calculate each subject's mean Baseline RT # and subtract that value from experimental RTs dat_out2 <- pupilr::calculate_adjustedRT_f(dat_out2) dat_noexp_out2 <- pupilr::calculate_adjustedRT_f(dat_noexp_out2) d.wording_out2 <- pupilr::calculate_adjustedRT_f(d.wording_out2)
In this somewhat replication of C&G'04, we had participants listen to chunks of a story read by their exposure talkers between each block of the RT task, with the intent that this would give them even more exposure to adapt to, letting us study accent adaptation over a longer scale. We had planned to look at how pupil dilation corresponds to accent adaptation after this behavioral study turned out ok.
Although the experiments without the exposure passages seemed encouraging, the results from the experiments with exposure passages interspersed with the C&G task were somewhat confusing--each time the participants came back to the task, it seemed that their response times reset, especially for the first one or two trials. These non-linearities prompted us to approach the problem with GAM models. We discussed using the three exposure blocks for each condition to make predictions about the test block, but I can't quite remember the comparisons we had planned.
Regardless, I've implemented a few different sorts of predictions that the GAM model could make about the test block---namely how the model would predict the responses if the participants heard the test block with the same talker as exposure (i.e. participants in the unaccented condition heard the same unaccented talker for test) and if all the participants heard the accented talker in test (as was done in the actual experiment). Predicting the data in the test block given only the exposure (during which each participant only heard one talker) lets us establish a baseline to compare the actual results.
Just the nitty-gritty about what I used in the GAM model. Feel free to skip.
d.gamm <- dat_out2 %>% filter(Condition %in% c("accented","unaccented") & Run %in% c("Exp2v1","ExpPilotv1","Exp2v2")) %>% within(., { Word <- factor(Word) WorkerId <- factor(WorkerId) # sum coding for accent condition Condition <- factor(Condition, levels = c("unaccented","accented")) contrasts(Condition) <- cbind("accented" = c(-1,1)) # Make the unordered version of Condition Condition.unordered <- Condition # Make the ordered version Condition <- ordered(Condition) Trial <- TrialInExperiment-max(Trial[Block=="0"]) Trial.c <- Trial-12.5 BlockTrial <- Trial-floor((Trial-1)/8)*8 BlockTrial.c <- BlockTrial-mean(BlockTrial) }) %>% filter(Block %in% c("1","2","3")) %>% within(., { # Make Block an ordered factor with helmert coding Block <- factor(Block) Block <- as.ordered(Block) contrasts(Block) <- contr.helmert(3) colnames(contrasts(Block)) <- c("Block:2v1","Block:3v12") # Make a continuous variable `Block` BlockContinuous <- ifelse(Block=="1",1, ifelse(Block=="2",2, ifelse(Block=="3",3,NA))) # Make a factor of the combination of `Block` and `Condition` BlockXCondition <- factor(paste0(Condition,"-",Block)) # sum contrast code List (counterbalancing nuisance factor) List <- factor(List) contrasts(List) <- contr.sum(nlevels(List)) colnames(contrasts(List)) <- rownames(contrasts(List))[1:(length(levels(List))-1)] # sum code ListID # Currently considers the new runs completely unrelated ListID <- sapply(List,function(x){as.numeric(strsplit(as.character(x),"_")[[1]][3])}, USE.NAMES = FALSE) ListID <- factor(ifelse(Run=="Exp2v2", paste0(ListID,"_new"), ListID)) contrasts(ListID) <- contr.sum(nlevels(ListID)) #sum code ListOrder ListOrder <- factor(sapply(List,function(x){strsplit(as.character(x),"_")[[1]][4]}, USE.NAMES = FALSE)) contrasts(ListOrder) <- contr.sum(nlevels(ListOrder)) }) %>% as.data.frame()
d.gamm.ex1 <- dat_noexp_out2 %>% filter(Condition %in% c("accented","unaccented") & Run %in% c("Exp1V1", "PilotV3", "PilotV2")) %>% within(., { Word <- factor(Word) WorkerId <- factor(WorkerId) # sum coding for accent condition Condition <- factor(Condition, levels = c("unaccented","accented")) contrasts(Condition) <- cbind("accented" = c(-1,1)) # Make the unordered version of Condition Condition.unordered <- Condition # Make the ordered version Condition <- ordered(Condition) Trial <- TrialInExperiment-max(Trial[Block=="0"]) Trial.c <- Trial-mean(Trial) BlockTrial <- Trial-floor((Trial-1)/8)*8 BlockTrial.c <- BlockTrial-mean(BlockTrial) }) %>% filter(Block %in% c("1","2","3")) %>% within(., { # Make Block an ordered factor with helmert coding Block <- factor(Block) Block <- as.ordered(Block) contrasts(Block) <- contr.helmert(3) colnames(contrasts(Block)) <- c("Block:2v1","Block:3v12") # Make a continuous variable `Block` BlockContinuous <- ifelse(Block=="1",1, ifelse(Block=="2",2, ifelse(Block=="3",3,NA))) # Make a factor of the combination of `Block` and `Condition` BlockXCondition <- factor(paste0(Condition,"-",Block)) # sum contrast code List (counterbalancing nuisance factor) List <- factor(List) contrasts(List) <- contr.sum(nlevels(List)) colnames(contrasts(List)) <- rownames(contrasts(List))[1:(length(levels(List))-1)] # sum code ListID # Currently considers the new runs completely unrelated # ListID <- sapply(List,function(x){as.numeric(strsplit(as.character(x),"_")[[1]][3])}, # USE.NAMES = FALSE) # ListID <- factor(ifelse(Run=="Exp2v2", # paste0(ListID,"_new"), # ListID)) # contrasts(ListID) <- contr.sum(nlevels(ListID)) # # #sum code ListOrder # ListOrder <- factor(sapply(List,function(x){strsplit(as.character(x),"_")[[1]][4]}, # USE.NAMES = FALSE)) # contrasts(ListOrder) <- contr.sum(nlevels(ListOrder)) }) %>% as.data.frame() d.gamm.wording <- d.wording_out2 %>% # filter(Condition %in% c("accented","unaccented") & # Run %in% c("Exp1V1", "PilotV3", "PilotV2")) %>% within(., { Word <- factor(Word) WorkerId <- factor(WorkerId) # sum coding for accent condition Condition <- factor(Condition, levels = c("unaccented","accented")) contrasts(Condition) <- cbind("accented" = c(-1,1)) # Make the unordered version of Condition Condition.unordered <- Condition # Make the ordered version Condition <- ordered(Condition) Trial <- TrialInExperiment-max(Trial[Block=="0"]) Trial.c <- Trial-mean(Trial) BlockTrial <- Trial-floor((Trial-1)/8)*8 BlockTrial.c <- BlockTrial-mean(BlockTrial) }) %>% filter(Block %in% c("1","2","3")) %>% within(., { # Make Block an ordered factor with helmert coding Block <- factor(Block) Block <- as.ordered(Block) contrasts(Block) <- contr.helmert(3) colnames(contrasts(Block)) <- c("Block:2v1","Block:3v12") # Make a continuous variable `Block` BlockContinuous <- ifelse(Block=="1",1, ifelse(Block=="2",2, ifelse(Block=="3",3,NA))) # Make a factor of the combination of `Block` and `Condition` BlockXCondition <- factor(paste0(Condition,"-",Block)) # sum contrast code List (counterbalancing nuisance factor) List <- factor(List) contrasts(List) <- contr.sum(nlevels(List)) colnames(contrasts(List)) <- rownames(contrasts(List))[1:(length(levels(List))-1)] # sum code ListID # Currently considers the new runs completely unrelated # ListID <- sapply(List,function(x){as.numeric(strsplit(as.character(x),"_")[[1]][3])}, # USE.NAMES = FALSE) # ListID <- factor(ifelse(Run=="Exp2v2", # paste0(ListID,"_new"), # ListID)) # contrasts(ListID) <- contr.sum(nlevels(ListID)) # # #sum code ListOrder # ListOrder <- factor(sapply(List,function(x){strsplit(as.character(x),"_")[[1]][4]}, # USE.NAMES = FALSE)) # contrasts(ListOrder) <- contr.sum(nlevels(ListOrder)) }) %>% as.data.frame()
I begin by making my assumptions clear (or as clear as I can). The model will predict AdjustedRT
(the adjusted reaction times) with the following linear predictors:
Linear predictors:
BlockContinuous
), 1-3, so we can extend our predictions to Block 4, as a linear predictor~~~Condition
), a binary centered ordered factor (1 and -1), as a linear predictor (so that we can more easily test differences between conditions)Block
and Condition
, as a linear predictor~~~Non-linear predictors:
r font_color("Block number (BlockContinuous), 1-3, so we can extend our predictions to Block 4, as a linear predictor.", color="red")
I decided to not have Block number as a linear and non-linear predictor.BlockTrial
), 1-8, then centered around the mean (i.e. 4.5), as a continuous variable, as a non-linear predictorBlock
and BlockTrial
, so that different smooths for BlockTrial
are fitted to each level of Block
, as a non-linear predictor. Because Block
is an ordered factor, the smooths for each level of the factor will be a "reference" smooth for the first level and "difference" smooths for all the other levels, representing the difference from the reference. (I'm not actually sure how/if the coding of the factor actually does anything if each level gets its own smooth. In fact, I'm pretty sure the coding doesn't matter at all in the interaction, based on examples from the package)Condition
and BlockTrial
, so that different smooths for BlockTrial
are fitted to each level of Condition
BlockTrial
and the combination of Condition
and Block
, so that different smooths for BlockTrial
are fitted to each combination of Condition
and Block
Random effects:
BlockTrial
each participant (WorkerId
), so that each participant has their own smooth fitted for BlockTrial
(also maybe: the same smoothing parameter for each participant)BlockTrial
only has 8 values (1-8) and BlockContinuous
only has 3 (1-3). Jacolien van Rij says the following:
"Note that the model will base the number of base functions (reflected in the
edf
of the summary) on the data with the setting fork
as upper bound. By default, the value ofk
fors()
is around 9, and forte()
andti()
5 per dimension. Importantly, the value ofk
should be at most one less than the number of unique data points, otherwise it will fit the density of that predictor.
The mgcv
documentation says the following:
"When setting up models in the
mgcv
package, usings
orte
terms in a model formula,k
must be chosen: the defaults are essentially arbitrary.
To me, this means that since our variables have so few different values, we should manually set k
so that it's one less than the number of unique data points for all non-linear predictors, and completely ignore defaults.
In the current model, I've included both linear and non-linear versions of BlockContinuous
and Condition
.
There's also the possibility that something about the way I fit the model is unoptimal. There are a bunch of finnicky things that I tried to account for, but I'm no expert, so I might have gotten them wrong. For example, I tried to set the k basis functions as well as I understand them, but I'm not 100% confident in them. Likewise, choices such as whether to center certain variables, whether to use factors vs. continuous predictors vs. ordered factors, etc.
The aspect I'm least sure is the d
parameter in mgcv::ti
. The predictors in my interactions are not on the same scale, and I think the default for mgcv::ti
assumes that, but I can't really make heads or tails of the documentation.
This model is only fit to correct responses.
m.full.block_contin <- gam(AdjustedRT ~ # Linear predictors Condition + #BlockContinuous + BlockContinuous:Condition + # Non linear predictors ti(BlockContinuous, k=3) + ti(BlockTrial, k=7) + ti(BlockContinuous, by=Condition, k=3) + ti(BlockTrial, BlockContinuous, k=3) + ti(BlockTrial, by=Condition, k=7) + ti(BlockTrial, BlockContinuous, by=Condition, k=3) + # Random effects s(WorkerId, bs="re") + #Random intercept s(WorkerId, BlockTrial, bs="re"), # Random slope data = d.gamm %>% filter(BinaryCorrect==1) ) summary(m.full.block_contin)
These are plots of the components of the model. I can't really decipher the interaction plots, but the splines seem relatively interpretable.
I made the Condition
variable an ordered factor in the model---from what I understand, this means that the splines that have :Conditionaccented
in the y axis are "difference" smooths---how the accented
splines differ from the unaccented
ones.
It seems that Condition
doesn't really change the slope of the effect of Block number (BlockContinuous
), although it seems to be doing something with the effect of the trial number within the block (BlockTrial
).
plot(m.full.block_contin,pages=2)
# # # The model: # gam(AdjustedRT ~ # # Linear predictors # Block + Condition + Block:Condition + # # Non linear predictors # ti(BlockTrial, k=7) + # ti(BlockTrial, by=Block, k=7) + # ti(BlockTrial, by=Condition, k=7) + # ti(BlockTrial, by=BlockXCondition, k=7) + # # Random effects # s(WorkerId, bs="re") + #Random intercept # s(WorkerId, BlockTrial, bs="re") + # Random slope # # Possible!!?!? # s(BlockTrial, WorkerId, bs="fs", m=1), # data=d.gamm) # # m.full <- gam(AdjustedRT ~ # # Linear predictors # Block*Condition + # # Non linear predictors # ti(BlockTrial, k=7) + # ti(BlockTrial, by=Block, k=7) + # ti(BlockTrial, by=Condition, k=7) + # ti(BlockTrial, by=BlockXCondition, k=7) + # # Random effects # s(WorkerId, bs="re") + #Random intercept # s(WorkerId, BlockTrial, bs="re"), # data=d.gamm) # summary(m.full) # plot(m.full,pages=2) # # # m.full.block_contin <- gam(AdjustedRT ~ # # Linear predictors # Condition*BlockContinuous + # # Non linear predictors # ti(BlockContinuous, k=3) + # ti(BlockTrial, k=7) + # ti(BlockContinuous, by=Condition, k=3) + # ti(BlockTrial, BlockContinuous, k=3) + # ti(BlockTrial, by=Condition, k=7) + # ti(BlockTrial, BlockContinuous, by=Condition, k=3) + # # Random effects # s(WorkerId, bs="re") + #Random intercept # s(WorkerId, BlockTrial, bs="re"), # Random slope # data=d.gamm #%>% # #mutate(Condition=ifelse(Condition=="accented",1,0)) # ) # summary(m.full.block_contin) # # m.simple <- gam(AdjustedRT ~ # # Linear predictors # Block + Condition + # # Non linear predictors # ti(BlockTrial, k=7) + # ti(BlockTrial, by=Block, k=7) + # ti(BlockTrial, by=Condition, k=7) + # # Random effects # s(WorkerId, bs="re") + #Random intercept # s(WorkerId, BlockTrial, bs="re"), # data=d.gamm) # summary(m.simple) # plot(m.simple,pages=1)
# Have participants hear a test block with the same accent d.gamm.fake_extended.block4 <- d.gamm %>% filter(BlockContinuous==3) %>% mutate(Trial=Trial+8, BlockContinuous=4, Block=factor("4",levels=c("1","2","3","4")), Phase = factor("Test phase",levels=c("Exposure phase","Test phase"))) # Have participants hear a test block like in the experiment d.gamm.fake_experiment.block4 <- d.gamm %>% filter(BlockContinuous==3) %>% mutate(Trial=Trial+8, BlockContinuous=4, Block=factor("4",levels=c("1","2","3","4")), Phase = factor("Test phase",levels=c("Exposure phase","Test phase")), OrigCondition=Condition, Condition=factor("accented",levels=c("unaccented","accented"))) # Getting the actual data d.gamm.real.block4 <- dat_out2 %>% filter(Condition %in% c("accented","unaccented") & Run %in% c("Exp2v1","ExpPilotv1","Exp2v2")) %>% within(., { WorkerId <- factor(WorkerId) # sum coding for accent condition Condition <- factor(Condition, levels = c("unaccented","accented")) contrasts(Condition) <- cbind("accented" = c(-1,1)) # Make the unordered version of Condition Condition.unordered <- Condition # Make the ordered version Condition <- ordered(Condition) Trial <- TrialInExperiment-max(Trial[Block=="0"]) BlockTrial <- Trial-floor((Trial-1)/8)*8 BlockTrial.c <- BlockTrial-mean(BlockTrial) }) %>% filter(Block=="4") %>% mutate(BlockContinuous=4, Block=factor("4",levels=c("1","2","3","4")), OrigCondition=Condition)
d.gamm$Predict <- m.full.block_contin %>% predict.gam(d.gamm) d.gamm.fake_extended.block4$Predict <- m.full.block_contin %>% predict.gam(d.gamm.fake_extended.block4) d.gamm.fake_experiment.block4$Predict <- m.full.block_contin %>% predict.gam(d.gamm.fake_experiment.block4) d.gamm.fake_experiment.block4$Condition<-d.gamm.fake_experiment.block4$OrigCondition d.gamm.compare <- bind_rows( "Actual" = bind_rows(d.gamm, d.gamm.real.block4) %>% mutate(Data=AdjustedRT), "PredSameAcc" = bind_rows(d.gamm, d.gamm.fake_extended.block4) %>% mutate(Data=Predict), "PredExperiment" = bind_rows(d.gamm, d.gamm.fake_experiment.block4) %>% mutate(Data=Predict), .id="DataType")
```r"}
g1<-NULL
rt_coord_lim <- c(-150, 1000) rt_y_label <- "Adjusted RTs (ms)\n(experiment RTs - baseline RT)"
block_fill_palette <- rep(c("#f1f1f1", "white"),100)[1:(length(unique(d.gamm.compare$BlockContinuous))* length(unique(d.gamm.compare$DataType)))] fatten_val <- 5 base_font_size <- 15 strip_text_size <- 14
palette_exp2 <- c("red", "grey30")
rects <- d.gamm.compare %>% filter(!is.na(Phase)) %>% group_by(Block) %>% summarise(Phase=first(Phase), xstart_rect = min(Trial)-0.5, xend_rect = max(Trial)+0.5, x_seg = mean(Trial), y_seg_RTs = rt_coord_lim[2], y_seg_errors = 0.5) %>% mutate(Block=paste("Block",Block))
p1 <- d.gamm.compare %>% #filter(!(Trial %in% c(1,9,17,25))) %>% filter(BinaryCorrect == 1) %>% filter(Block %in% c("1","2","3","4")) %>% ggplot(aes(x = Trial, y = Data, colour = Condition)) + facet_grid(DataType~Phase, scales = "free_x") + geom_rect(data = rects, aes(xmin = xstart_rect, xmax = xend_rect, ymin = -Inf, ymax = Inf), fill = block_fill_palette, inherit.aes = FALSE) + geom_text(data = rects, aes(x = x_seg, y = y_seg_RTs, label = Block), inherit.aes = FALSE) + stat_summary(fun.data = "mean_cl_boot", geom = "pointrange", position = position_dodge(.15), fatten = fatten_val, alpha = .7) + stat_smooth(method = "lm", show.legend = FALSE) + geom_hline(yintercept = 0, linetype = "dashed") + scale_x_continuous("Trial", breaks = seq(1,32,1), expand = c(0,0)) + scale_y_continuous(rt_y_label) + scale_colour_manual("Exposure condition", labels = c(accented="Indian-accented", unaccented="Control (unaccented)"), values = color_values) + scale_fill_manual("Exposure condition", labels = c(accented="Indian-accented", unaccented="Control (unaccented)"), values = color_values) + labs(y = rt_y_label) + guides(colour = guide_legend(override.aes = list(alpha = 1))) + theme_bw(base_size = base_font_size) + theme(panel.border = element_blank(), strip.text = element_text(size = strip_text_size), legend.position = "top") + geom_rangeframe(colour = "black") + coord_cartesian(ylim=rt_coord_lim)
g1 <- ggplotGrob(p1) g1$widths[[4]] <- unit(3, "null")
library(grid) grid.newpage() grid.draw(g1)
```r"} rull_story<-NULL # axis defaults rt_y_label <- "Unadjusted RTs (ms)\n(experiment RTs - baseline RT)" # aesthetic defaults dodge_amt <- 0.9 dodge_amt2 <- 0.0 base_font_size <- 13 strip_text_size <- 11 rull_story <- d.gamm.compare %>% filter(BinaryCorrect == 1) %>% filter(Block %in% c("1","2","3","4")) %>% #filter(!(Trial %in% c(1,9,17,25))) %>% ggplot(aes(x = Trial, y = Data, colour = Condition, group=paste0(Block,Condition,DataType))) + facet_wrap(~ DataType, scales = "free_x",nrow=3) + geom_rect(data = rects, aes(xmin = xstart_rect, xmax = xend_rect, ymin = -Inf, ymax = Inf), fill = rep(block_fill_palette,1), inherit.aes = FALSE) + geom_text(data = rects, aes(x = x_seg, y = y_seg_RTs-300, label = Block), inherit.aes = FALSE) + stat_smooth(method = "lm", show.legend = FALSE) + geom_hline(yintercept = 0, linetype = "dashed") + scale_x_continuous("", breaks = NULL, expand = c(0,0)) + scale_y_continuous(rt_y_label) + scale_colour_manual("Exposure condition", labels = c(accented="Indian-accented", unaccented="Control (unaccented)"), values = color_values) + scale_fill_manual("Exposure condition", labels = c(accented="Indian-accented", unaccented="Control (unaccented)"), values = color_values) + labs(y = rt_y_label) + guides(colour = guide_legend(override.aes = list(alpha = 1))) + theme_bw(base_size = base_font_size) + theme(panel.border = element_blank(), strip.text = element_text(size = strip_text_size), legend.position = "top") + geom_rangeframe(colour = "black") + coord_cartesian(ylim=c(-150,750)) rull_story
These results don't strike me as super encouraging, at least in the sense that they don't seem to suggest that there's adaptation in the exposure phase, which might matter for the pupillometry. However, it might be interesting to compare these results (which come from the experiments where there were passive listening sections between blocks) to the results where there were no passive listening blocks (which seemed to have much stronger results).
Like I said earlier, I'm less certain of how to use the actual test block data to measure adaptation in the context of a GAM model. Would it be fair to compare the actual results with the predicted values?
\FloatBarrier
m.ex1.full.block_contin <- gam(AdjustedRT ~ # Linear predictors Condition + #BlockContinuous + BlockContinuous:Condition + # Non linear predictors ti(BlockContinuous, k=3) + ti(BlockTrial, k=7) + ti(BlockContinuous, by=Condition, k=3) + ti(BlockTrial, BlockContinuous, k=3) + ti(BlockTrial, by=Condition, k=7) + ti(BlockTrial, BlockContinuous, by=Condition, k=3) + # Random effects s(WorkerId, bs="re") + #Random intercept s(WorkerId, BlockTrial, bs="re"), # Random slope data = d.gamm.ex1 %>% filter(BinaryCorrect==1) ) summary(m.ex1.full.block_contin) #summary(m..full.block_contin) plot(m.ex1.full.block_contin,pages=2)
\FloatBarrier
m.wording.full.block_contin <- gam(AdjustedRT ~ # Linear predictors Condition + #BlockContinuous + BlockContinuous:Condition + # Non linear predictors ti(BlockContinuous, k=3) + ti(BlockTrial, k=7) + ti(BlockContinuous, by=Condition, k=3) + ti(BlockTrial, BlockContinuous, k=3) + ti(BlockTrial, by=Condition, k=7) + ti(BlockTrial, BlockContinuous, by=Condition, k=3) + # Random effects s(WorkerId, bs="re") + #Random intercept s(WorkerId, BlockTrial, bs="re"), # Random slope data = d.gamm.wording %>% filter(BinaryCorrect==1) ) summary(m.wording.full.block_contin)
m.full.no_block <- gam(AdjustedRT ~ # Linear predictors Condition + #BlockContinuous + BlockContinuous:Condition + # Non linear predictors ti(Trial.c) + ti(BlockTrial, k=7) + ti(Trial.c, by=Condition) + ti(BlockTrial, Trial.c, k=7) + ti(BlockTrial, by=Condition, k=7) + ti(BlockTrial, Trial.c, by=Condition, k=7) + # Random effects s(WorkerId, bs="re") + #Random intercept s(WorkerId, BlockTrial, bs="re"), # Random slope data = d.gamm %>% filter(BinaryCorrect==1) ) summary(m.full.no_block) plot(m.full.no_block,pages=2)
\FloatBarrier
m.ex1.full.no_block <- gam(AdjustedRT ~ # Linear predictors Condition + #BlockContinuous + BlockContinuous:Condition + # Non linear predictors ti(Trial.c) + ti(BlockTrial, k=7) + ti(Trial.c, by=Condition) + ti(BlockTrial, Trial.c, k=7) + ti(BlockTrial, by=Condition, k=7) + ti(BlockTrial, Trial.c, by=Condition, k=7) + # Random effects s(WorkerId, bs="re") + #Random intercept s(WorkerId, BlockTrial, bs="re"), # Random slope data = d.gamm.ex1 %>% filter(BinaryCorrect==1) ) summary(m.ex1.full.no_block) plot(m.ex1.full.no_block,pages=2)
\FloatBarrier
m.wording.full.no_block <- gam(AdjustedRT ~ # Linear predictors Condition + #BlockContinuous + BlockContinuous:Condition + # Non linear predictors ti(Trial.c) + ti(BlockTrial, k=7) + ti(Trial.c, by=Condition) + ti(BlockTrial, Trial.c, k=7) + ti(BlockTrial, by=Condition, k=7) + ti(BlockTrial, Trial.c, by=Condition, k=7) + # Random effects s(WorkerId, bs="re") + #Random intercept s(WorkerId, BlockTrial, bs="re"), # Random slope data = d.gamm.ex1 %>% filter(BinaryCorrect==1) ) summary(m.wording.full.no_block)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.