Load required packages

library(reshape2)
library(tables)
library(ez)
library(knitr)
library(plyr)
sessionInfo()
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## locale:
## [1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] plyr_1.8.4      knitr_1.15.1    ez_4.3          tables_0.7.92  
##  [5] Hmisc_3.17-4    ggplot2_2.2.1   Formula_1.2-1   survival_2.39-5
##  [9] lattice_0.20-34 reshape2_1.4.1 
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.6         nloptr_1.0.4        RColorBrewer_1.1-2 
##  [4] tools_3.3.2         lme4_1.1-12         rpart_4.1-10       
##  [7] digest_0.6.10       nlme_3.1-128        evaluate_0.10      
## [10] tibble_1.2          gtable_0.2.0        mgcv_1.8-15        
## [13] Matrix_1.2-7.1      parallel_3.3.2      yaml_2.1.13        
## [16] SparseM_1.7         gridExtra_2.2.1     stringr_1.0.0      
## [19] cluster_2.0.5       MatrixModels_0.4-1  rprojroot_1.1      
## [22] grid_3.3.2          nnet_7.3-12         data.table_1.9.6   
## [25] foreign_0.8-67      rmarkdown_1.3       minqa_1.2.4        
## [28] latticeExtra_0.6-28 car_2.1-3           magrittr_1.5       
## [31] MASS_7.3-45         backports_1.0.4     scales_0.4.1       
## [34] htmltools_0.3.5     splines_3.3.2       pbkrtest_0.4-6     
## [37] assertthat_0.1      colorspace_1.2-6    quantreg_5.26      
## [40] stringi_1.1.1       acepack_1.3-3.3     lazyeval_0.2.0     
## [43] munsell_0.4.3       chron_2.3-47

Functions

get_Boxplot

get_Boxplot <- function (Data, DV, yMin=0, yMax=1, yBreaks=0.2, yLabel="Test") {
  ggplot(Data, aes_string(x="factor(Group)", y=DV, fill="Condition")) +
    geom_boxplot(outlier.colour=NA) +
    scale_fill_brewer(palette="BuGn") + 
    theme_minimal(base_size = 16, base_family = "") +
    scale_x_discrete("") +
    geom_point(position=position_jitterdodge(dodge.width=0.75, jitter.width=0.3), alpha=0.2, size=3) +
    scale_y_continuous(yLabel,limits=c(yMin,yMax), breaks=seq(yMin,yMax,yBreaks)) 
}

Data wrangling

Extract accuracy data from trial report.

  AccData <- read.delim("data/Trial_Report_Accuracy.txt", stringsAsFactors=FALSE) # Read in data for premature saccades
  Data <- AccData[,c("RECORDING_SESSION_LABEL", "trial_id_num", "INDEX", "ERROR_TYPE")] # Select relevant columns
  Data$SubNum <- as.numeric(lapply(Data$RECORDING_SESSION_LABEL, function(x) substr(x, 1, 3))) #Extract SubNum from first three characters of recording session label (and convert string to numeric)
  Data$UTI <- Data$SubNum*1000 + Data$trial_id_num # Calculate unique trial identifier

Extract saccadic reaction times from interest area report and add to trial data.

  SRTData <- read.delim("data/IA_Report_RJA_SRT.txt", stringsAsFactors=FALSE) # Read in data for SRTs
  SRTData <- SRTData [(SRTData$IA_LABEL=="IA_BURGLAR "),] # Filter out all interest areas apart from Burglar
  SRTData$IPStartTime <- SRTData$IP_START_TIME - SRTData$TRIAL_START_TIME # Calculate the IP start time relative to the trial start time
  SRTData$IA_FIRST_SACCADE_START_TIME[SRTData$IA_FIRST_SACCADE_START_TIME=="."] <- NA # Convert periods (missing data) to NA
  SRTData$IA_FIRST_SACCADE_START_TIME <- as.numeric(SRTData$IA_FIRST_SACCADE_START_TIME) # Convert to numeric (because missing data was represented by "." all numbers in the column were treated as strings)
  SRTData$SRT <- SRTData$IA_FIRST_SACCADE_START_TIME - (SRTData$IP_START_TIME - SRTData$TRIAL_START_TIME) # Subtract IP Start Time from Saccade time to calculate SRT
  SRTData$SubNum <- as.numeric(lapply(SRTData$RECORDING_SESSION_LABEL, function(x) substr(x, 1, 3))) #Extract SubNum from first three characters of recording session label (and convert string to numeric)
  SRTData$UTI <- SRTData$SubNum*1000 + SRTData$trial_id_num # Calculate unique trial identifier
  Data$FirstSaccadeSRT <- SRTData$SRT[match(Data$UTI, SRTData$UTI)] # Add SRT data to main dataframe using UTI

Extract dwell times from interest area report and add to trial data.

  DwellData <- read.delim("data/IA_Report_IJA_Dwell.txt", stringsAsFactors=FALSE) # Read in data for burglar dwelltime
  DwellData <- DwellData [(DwellData$IA_LABEL=="IA_BURGLAR "),] # Filter out all interest areas apart from Burglar
  DwellData$SubNum <- as.numeric(lapply(DwellData$RECORDING_SESSION_LABEL, function(x) substr(x, 1, 3))) #Extract SubNum from first three characters of recording session label (and convert string to numeric)
  DwellData$UTI <- DwellData$SubNum*1000 + DwellData$trial_id_num # Calculate unique trial identifier
  Data$BurglarDwellTime <- DwellData$IA_DWELL_TIME[match(Data$UTI, DwellData$UTI)] # Add Dwelltime data to main dataframe using UTI

Extract premature saccade data from interest area report and add to trial data.

  PremiData <- read.delim("data/IA_Report_IJA_Premi.txt", stringsAsFactors=FALSE) # Read in data for premature saccades
  PremiData <- PremiData [(PremiData$IA_LABEL=="IA_BURGLAR "),] # Filter out all interest areas apart from Burglar
  PremiData$PrematureSaccade <- lapply(PremiData$IA_FIXATION_COUNT, function(x) if(x==0){0}else{1}) # Identify premature saccades (Fixation count is >1)
  PremiData$SubNum <- as.numeric(lapply(PremiData$RECORDING_SESSION_LABEL, function(x) substr(x, 1, 3)))
  PremiData$UTI <- PremiData$SubNum*1000 + PremiData$trial_id_num
  Data$PrematureSaccade <- PremiData$PrematureSaccade[match(Data$UTI, PremiData$UTI)]
  
  Data$PrematureSaccade <- vapply(Data$PrematureSaccade, paste, collapse = ", ", character(1L)) # For some reason PrematureSaccade is identified as a list. This command "flattens" the list, allowing it to be written to csv.

Decode correct response / error type

  Data$CorrectResponse <- ifelse(Data$ERROR_TYPE==0,1,0) 
  Data$SearchError <- ifelse(Data$ERROR_TYPE==1,1,0) 
  Data$TimeOut <- ifelse(Data$ERROR_TYPE==2,1,0) 
  Data$LocationError <- ifelse(Data$ERROR_TYPE==3,1,0) 
  Data$RecalibrationError <- ifelse(Data$ERROR_TYPE==4,1,0)

Add information about trials.

  TrialInfo <- read.csv("data/TrialInfo.csv", stringsAsFactors=FALSE)
  Data$Condition <- TrialInfo$Condition[match(Data$trial_id_num, TrialInfo$trial_id_num)] # Add Condition (RJA RJAc IJA IJAc)
  Data$Interface <- TrialInfo$Interface[match(Data$trial_id_num, TrialInfo$trial_id_num)] # Add Interface (Test or Control)
  Data$SubjectRole <- TrialInfo$SubjectRole[match(Data$trial_id_num, TrialInfo$trial_id_num)] # Add SubjectRole (Respond or Initiate)

Add information about subjects.

  SubjectData <- read.csv("data/SubjectInfo.csv", stringsAsFactors=FALSE)
  Data$Group <- SubjectData$Group[match(Data$SubNum, SubjectData$SubNum)] # Add Group membership
  Data$Block1 <- SubjectData$Block.1[match(Data$SubNum, SubjectData$SubNum)] # Add Block1 Session Name
  Data$Block2 <- SubjectData$Block.2[match(Data$SubNum, SubjectData$SubNum)] # Add Block2 Session Name
  Data$Block <- ifelse(Data$RECORDING_SESSION_LABEL==Data$Block1,"1",ifelse(Data$RECORDING_SESSION_LABEL==Data$Block2,"2",NA)) # Determine whether current Session matches Block1 or Block2

Exclude two participants who weren’t deceived.

  Data$SubInclusion <- SubjectData$SubjectInclude[match(Data$SubNum, SubjectData$SubNum)] # Add info about excluding participants
  Data <- Data [(Data$SubInclusion==1),] # Exclude participants

Rename dataframe columns for consistency with other datasets.

  names(Data)[names(Data)=="trial_id_num"] <- "TrialID" # Rename TrialID column
  names(Data)[names(Data)=="INDEX"] <- "TrialNum" # Rename TrialNum column
  #Data$ScreenIn <- 1 # Add column to be used for screening data (default is 1)

Remove columns that are now redundant.

Data <- subset(Data, select = -c(Block1, Block2, SubInclusion, RECORDING_SESSION_LABEL, ERROR_TYPE))

Write trial data to CSV file.

  write.csv(Data, file="output/data/Data.csv", row.names=FALSE)

Identify factors

FactorNames <- c("SubNum", "TrialID", "Block", "Interface", "SubjectRole", "Group", "Condition")
Data[,FactorNames] <- colwise(as.factor)(Data[,FactorNames])
str(Data)
## 'data.frame':    10368 obs. of  17 variables:
##  $ TrialID           : Factor w/ 432 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ TrialNum          : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ SubNum            : Factor w/ 48 levels "101","102","103",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ UTI               : num  103001 103002 103003 103004 103005 ...
##  $ FirstSaccadeSRT   : num  NA NA 486 NA NA 504 627 NA NA NA ...
##  $ BurglarDwellTime  : int  0 2956 0 1256 698 0 0 636 624 1938 ...
##  $ PrematureSaccade  : chr  "0" "1" "0" "1" ...
##  $ CorrectResponse   : num  1 1 1 1 1 1 1 0 1 1 ...
##  $ SearchError       : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ TimeOut           : num  0 0 0 0 0 0 0 1 0 0 ...
##  $ LocationError     : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ RecalibrationError: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Condition         : Factor w/ 4 levels "IJA","IJAc","RJA",..: 4 2 4 2 2 4 3 1 1 1 ...
##  $ Interface         : Factor w/ 2 levels "Control","Test": 1 1 1 1 1 1 2 2 2 2 ...
##  $ SubjectRole       : Factor w/ 2 levels "Initiate","Respond": 2 1 2 1 1 2 2 1 1 1 ...
##  $ Group             : Factor w/ 2 levels "Computer","Human": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Block             : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...

Respond condition

CurrentData <- Data [(Data$SubjectRole=="Respond"), ]
CurrentData$Condition <- factor(CurrentData$Condition)

Correct responses

CurrentData$ScreenIn <-ifelse((CurrentData$RecalibrationError==0 & CurrentData$SearchError==0 & !is.na(CurrentData$CorrectResponse)), 1, 0)
mean(CurrentData$ScreenIn)
## [1] 0.9876543
ScreenedData <- CurrentData [(CurrentData$ScreenIn==1), ]
BySubjectsData <- ddply(ScreenedData, .(Group, SubNum, Condition), summarise, ProportionCorrect = mean(CorrectResponse))
RJA_Acc_Boxplot <- get_Boxplot(BySubjectsData, "ProportionCorrect", yMin=0, yMax=1, yBreaks=0.2, yLabel="Proportion Correct")
RJA_Acc_Boxplot

Reaction times in Responding condition

Mean reaction times

CurrentData$ScreenIn <-ifelse((CurrentData$CorrectResponse==1 & CurrentData$FirstSaccadeSRT>=150 & CurrentData$FirstSaccadeSRT<=3000 & !is.na(CurrentData$FirstSaccadeSRT)), 1, 0)
mean(CurrentData$ScreenIn)
## [1] 0.8661265
ScreenedData <- CurrentData [(CurrentData$ScreenIn==1), ]

BySubjectsData <- ddply(ScreenedData, .(Group, SubNum, Condition), summarise, FirstSaccadeSRT = mean(FirstSaccadeSRT))
RJA_SRT_Boxplot <- get_Boxplot(BySubjectsData, "FirstSaccadeSRT", yMin=0, yMax=1200, yBreaks=200, yLabel="Reaction Time")
RJA_SRT_Boxplot

kable(ezANOVA(BySubjectsData, dv = FirstSaccadeSRT, wid = SubNum, within = .(Condition), between = Group, type = 3), row.names=FALSE, caption="First saccade reaction time", digits=3)
First saccade reaction time
Effect DFn DFd F p p<.05 ges
Group 1 46 5.710 0.021 * 0.075
Condition 1 46 264.629 0.000 * 0.664
Group:Condition 1 46 2.340 0.133 0.017

Median reaction times

CurrentData$ScreenIn <-ifelse((CurrentData$CorrectResponse==1 & CurrentData$FirstSaccadeSRT>0 & !is.na(CurrentData$FirstSaccadeSRT)), 1, 0)
BySubjectsData <- ddply(ScreenedData, .(Group, SubNum, Condition), summarise, FirstSaccadeSRT = median(FirstSaccadeSRT))

get_Boxplot(BySubjectsData, "FirstSaccadeSRT", yMin=0, yMax=1200, yBreaks=200, yLabel="Reaction Time")

kable(ezANOVA(BySubjectsData, dv = FirstSaccadeSRT, wid = SubNum, within = .(Condition), between = Group, type = 3), row.names=FALSE, caption="First saccade reaction time", digits=3)
First saccade reaction time
Effect DFn DFd F p p<.05 ges
Group 1 46 7.307 0.010 * 0.093
Condition 1 46 127.747 0.000 * 0.496
Group:Condition 1 46 2.185 0.146 0.017

Initiating condition

CurrentData <- Data [(Data$SubjectRole=="Initiate"), ]
CurrentData$Condition <- factor(CurrentData$Condition)

Correct responses

CurrentData$ScreenIn <-ifelse((CurrentData$RecalibrationError==0 & CurrentData$SearchError==0 & !is.na(CurrentData$CorrectResponse)), 1, 0)
mean(CurrentData$ScreenIn)
## [1] 0.9907407
ScreenedData <- CurrentData [(CurrentData$ScreenIn==1), ]
BySubjectsData <- ddply(ScreenedData, .(Group, SubNum, Condition), summarise, ProportionCorrect = mean(CorrectResponse))
IJA_Acc_Boxplot <- get_Boxplot(BySubjectsData, "ProportionCorrect", yMin=0, yMax=1, yBreaks=0.2, yLabel="Proportion Correct")
IJA_Acc_Boxplot

Dwell times

CurrentData$ScreenIn <-ifelse((CurrentData$CorrectResponse==1 & CurrentData$BurglarDwellTime>=150 & CurrentData$BurglarDwellTime<=3000 & !is.na(CurrentData$BurglarDwellTime)), 1, 0)
mean(CurrentData$ScreenIn)
## [1] 0.9371142
ScreenedData <- CurrentData [(CurrentData$ScreenIn==1), ]
BySubjectsData <- ddply(ScreenedData, .(Group, SubNum, Condition), summarise, BurglarDwellTime = mean(BurglarDwellTime))
IJA_Dwell_Boxplot <- get_Boxplot (BySubjectsData, "BurglarDwellTime", yMin=0, yMax=1800, yBreaks=200, yLabel="Dwell Time / ms")
IJA_Dwell_Boxplot

kable(ezANOVA(BySubjectsData, dv = BurglarDwellTime, wid = SubNum, within = .(Condition), between = Group, type = 3), row.names=FALSE, caption="ANOVA: Burglar Dwell Time", digits=3)
ANOVA: Burglar Dwell Time
Effect DFn DFd F p p<.05 ges
Group 1 46 0.055 0.816 0.001
Condition 1 46 24.361 0.000 * 0.043
Group:Condition 1 46 14.723 0.000 * 0.026 
TTest.Dwell <- data.frame(Measure=c("Dwell Time", "Dwell Time"),
                       Group=c("Computer", "Human"),
                       df=numeric(2),
                       t=numeric(2),
                       p=numeric(2))

for (i in 1:2) {
  Group <- TTest.Dwell$Group[i]
  TTest <- t.test(BurglarDwellTime ~ Condition, BySubjectsData[(BySubjectsData$Group==Group),], paired=TRUE)
  TTest.Dwell$df[i] <- TTest$parameter
  TTest.Dwell$t[i] <- TTest$statistic
  TTest.Dwell$p[i] <- TTest$p.value
}
kable(TTest.Dwell, row.names=FALSE, caption="T-tests: Effect of Condition on Dwell Times", digits=3)
T-tests: Effect of Condition on Dwell Times
Measure Group df t p
Dwell Time Computer 23 0.888 0.383
Dwell Time Human 23 5.581 0.000

Premature saccades

CurrentData$ScreenIn <-ifelse((CurrentData$RecalibrationError==0 & CurrentData$SearchError==0 & !is.na(CurrentData$PrematureSaccade)), 1, 0)
mean(CurrentData$ScreenIn)
## [1] 0.9907407
ScreenedData <- CurrentData [(CurrentData$ScreenIn==1), ]
ScreenedData$PrematureSaccade <- as.numeric(ScreenedData$PrematureSaccade) 
BySubjectsData <- ddply(ScreenedData, .(Group, SubNum, Condition), summarise, PrematureSaccades = mean(PrematureSaccade))
IJA_Premi_Boxplot <- get_Boxplot (BySubjectsData, "PrematureSaccades", yMin=0, yMax=1, yBreaks=0.2, yLabel="Proportion of Premature Saccades")
IJA_Premi_Boxplot

kable(ezANOVA(BySubjectsData, dv = PrematureSaccades, wid = SubNum, within = .(Condition), between = Group, type = 3), row.names=FALSE, caption="PrematureSaccade", digits=3)
PrematureSaccade
Effect DFn DFd F p p<.05 ges
Group 1 46 3.780 0.058 0.062
Condition 1 46 38.494 0.000 * 0.141
Group:Condition 1 46 13.788 0.001 * 0.055 
TTest.Premi <- data.frame(Measure=c("Premature Saccades", "Premature Saccades"),
                       Group=c("Computer", "Human"),
                       df=numeric(2),
                       t=numeric(2),
                       p=numeric(2))

for (i in 1:2) {
  Group <- TTest.Premi$Group[i]
  TTest <- t.test(PrematureSaccades ~ Condition, BySubjectsData[(BySubjectsData$Group==Group),], paired=TRUE)
  TTest.Premi$df[i] <- TTest$parameter
  TTest.Premi$t[i] <- TTest$statistic
  TTest.Premi$p[i] <- TTest$p.value
}
kable(TTest.Premi, row.names=FALSE, caption="T-tests: Effect of Condition on Premature Saccades", digits=3)
T-tests: Effect of Condition on Premature Saccades
Measure Group df t p
Premature Saccades Computer 23 2.109 0.046
Premature Saccades Human 23 6.145 0.000 
TTest.DF <- rbind(TTest.Dwell, TTest.Premi)
kable(TTest.DF, row.names=FALSE, caption="T-tests: Effect of Condition on DwellTime and Premature Saccades", digits=3)
T-tests: Effect of Condition on DwellTime and Premature Saccades
Measure Group df t p
Dwell Time Computer 23 0.888 0.383
Dwell Time Human 23 5.581 0.000
Premature Saccades Computer 23 2.109 0.046
Premature Saccades Human 23 6.145 0.000 
library(ReporteRs)
## Loading required package: ReporteRsjars
doc <- docx()
doc <- addTitle(doc, "Table 2: ", level=2)
doc <- addFlexTable( doc, vanilla.table(TTest.DF))
writeDoc(doc, file = "output/doc/Tables.docx")
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:Hmisc':
## 
##     combine
# Grabs the legend from one of the figures so it can be placed in a multi-panel figure
# Source: http://www.sthda.com/english/wiki/print.php?id=177
get_legend<-function(myggplot){
  tmp <- ggplot_gtable(ggplot_build(myggplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)
}

legend <- get_legend(RJA_Acc_Boxplot)
png(filename="output/figures/Respond.png", width=650, height=400)
grid.arrange (RJA_Acc_Boxplot + theme(legend.position="none"), RJA_SRT_Boxplot + theme(legend.position="none"), legend,
              ncol=3, nrow =1, widths=c(2.3, 2.3, 0.8), heights=c(2.3))
dev.off()
## quartz_off_screen 
##                 2
legend <- get_legend(IJA_Acc_Boxplot)
png(filename="output/figures/Initiate.png", width=950, height=400)
grid.arrange (IJA_Acc_Boxplot + theme(legend.position="none"), IJA_Dwell_Boxplot + theme(legend.position="none"), IJA_Premi_Boxplot + theme(legend.position="none"), legend,
  ncol=4, nrow =1, widths=c(2.3, 2.3, 2.3, 0.8), heights=c(2.3))
dev.off()
## quartz_off_screen 
##                 2