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)

Quick recap

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.

The GAMM stuff

Setting up the data frame

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

Logic of the model

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:

Non-linear predictors:

Random effects:

Potential issues:

Setting the $k$ basis functions

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 for k as upper bound. By default, the value of k for s() is around 9, and for te() and ti() 5 per dimension. Importantly, the value of k 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, using s or te 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.

Including linear and non-linear versions of variables

In the current model, I've included both linear and non-linear versions of BlockContinuous and Condition.

Finnicky details in GAM models

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.

The model

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)

GAMM plots

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")

Visualization of predicted values

```r"}

make sure you don't plot old plot

g1<-NULL

axis defaults

rt_coord_lim <- c(-150, 1000) rt_y_label <- "Adjusted RTs (ms)\n(experiment RTs - baseline RT)"

aesthetic defaults

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

exp2 defaults

palette_exp2 <- c("red", "grey30")

df containing shape info for background shading of blocks

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

RTs by trial and condition

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)

change grob widths so exposure phase is 3x wider than test phase

(i.e., since there are 3x more exposure trials)

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

Overall thoughts/feelings/questions

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

Results from older experiments

Exp1

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

Exp1b: Different wording

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)

Results with a different trial and block predictor

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

Similar results from previous experiments

Exp1

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

Exp1b: Different wording

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)


burchill/pupilr documentation built on May 22, 2019, 2:27 p.m.