1 Background

In this study, we compared the eye-movements of participants as they completed a joint attention game with a virtual partner or avatar. Half of the participants believed that the avatar was controlled by another human. The other half were correctly informed that the avatar was controlled by a computer algorithm.

The code imports trial and interest area reports and combines them into a single dataframe describing each trial. Accuracy and eyetracking data are presented via boxplots. Eyetracking data are analysed via ANOVAs and t-tests. We use non-parmetric stats to explore subject ratings of the task.

1.1 Reference

Caruana, N., Spirou, D., & Brock, J. (2017). Human agency beliefs influence joint attention behaviour. Manuscript submitted for publication.

2 Load required packages

library(reshape2)
library(tables)
library(ez)
library(knitr)
library(plyr)
library(gridExtra)
library(ReporteRs)
library(stringr)

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] stringr_1.2.0       ReporteRs_0.8.8     ReporteRsjars_0.0.2
##  [4] gridExtra_2.2.1     plyr_1.8.4          knitr_1.15.1       
##  [7] ez_4.4-0            tables_0.8          Hmisc_4.0-2        
## [10] ggplot2_2.2.1       Formula_1.2-1       survival_2.40-1    
## [13] lattice_0.20-34     reshape2_1.4.2     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.9         png_0.1-7           assertthat_0.1     
##  [4] rprojroot_1.2       digest_0.6.12       mime_0.5           
##  [7] R6_2.2.0            backports_1.0.5     acepack_1.4.1      
## [10] MatrixModels_0.4-1  evaluate_0.10       gdtools_0.1.3      
## [13] lazyeval_0.2.0      minqa_1.2.4         data.table_1.10.4  
## [16] SparseM_1.74        car_2.1-4           nloptr_1.0.4       
## [19] R.utils_2.5.0       R.oo_1.21.0         rpart_4.1-10       
## [22] Matrix_1.2-7.1      checkmate_1.8.2     rmarkdown_1.3      
## [25] splines_3.3.2       lme4_1.1-12         foreign_0.8-67     
## [28] htmlwidgets_0.8     munsell_0.4.3       shiny_1.0.0        
## [31] httpuv_1.3.3        base64enc_0.1-3     mgcv_1.8-15        
## [34] rvg_0.1.2           htmltools_0.3.5     nnet_7.3-12        
## [37] tibble_1.2          htmlTable_1.9       MASS_7.3-45        
## [40] R.methodsS3_1.7.1   grid_3.3.2          xtable_1.8-2       
## [43] nlme_3.1-128        gtable_0.2.0        magrittr_1.5       
## [46] scales_0.4.1        stringi_1.1.2       latticeExtra_0.6-28
## [49] xml2_1.1.1          RColorBrewer_1.1-2  tools_3.3.2        
## [52] parallel_3.3.2      pbkrtest_0.4-6      yaml_2.1.14        
## [55] colorspace_1.3-2    cluster_2.0.5       rJava_0.9-8        
## [58] quantreg_5.29

3 Functions

3.1 get_Boxplot

Draws boxplot with overlaid datapoints

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

3.2 get_legend

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

4 Data wrangling

We used Eyelink DataViewer software to export trial and interest area reports as text files.

Extract accuracy data from trial report.

  AccData <- read.delim("data/Trial_Report_Accuracy.txt", stringsAsFactors=FALSE)
  Data <- AccData[,c("RECORDING_SESSION_LABEL", "trial_id_num", "INDEX", "ERROR_TYPE")]
  Data$SubNum <- as.numeric(lapply(Data$RECORDING_SESSION_LABEL, function(x) substr(x, 1, 3)))
  Data$UTI <- Data$SubNum*1000 + Data$trial_id_num

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

  SRTData <- read.delim("data/IA_Report_RJA_SRT.txt", stringsAsFactors=FALSE)
  SRTData <- SRTData [(SRTData$IA_LABEL=="IA_BURGLAR "),]
  SRTData$IPStartTime <- SRTData$IP_START_TIME - SRTData$TRIAL_START_TIME
  SRTData$IA_FIRST_SACCADE_START_TIME[SRTData$IA_FIRST_SACCADE_START_TIME=="."] <- NA
  SRTData$IA_FIRST_SACCADE_START_TIME <- as.numeric(SRTData$IA_FIRST_SACCADE_START_TIME)
  SRTData$SRT <- SRTData$IA_FIRST_SACCADE_START_TIME - (SRTData$IP_START_TIME - SRTData$TRIAL_START_TIME)
  SRTData$SubNum <- as.numeric(lapply(SRTData$RECORDING_SESSION_LABEL, function(x) substr(x, 1, 3)))
  SRTData$UTI <- SRTData$SubNum*1000 + SRTData$trial_id_num
  Data$FirstSaccadeSRT <- SRTData$SRT[match(Data$UTI, SRTData$UTI)]

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

  DwellData <- read.delim("data/IA_Report_IJA_Dwell.txt", stringsAsFactors=FALSE) 
  DwellData <- DwellData [(DwellData$IA_LABEL=="IA_BURGLAR "),] 
  DwellData$SubNum <- as.numeric(lapply(DwellData$RECORDING_SESSION_LABEL, function(x) substr(x, 1, 3)))
  DwellData$UTI <- DwellData$SubNum*1000 + DwellData$trial_id_num 
  Data$BurglarDwellTime <- DwellData$IA_DWELL_TIME[match(Data$UTI, DwellData$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)
  PremiData <- PremiData [(PremiData$IA_LABEL=="IA_BURGLAR "),]
  PremiData$PrematureSaccade <- lapply(PremiData$IA_FIXATION_COUNT, function(x) if(x==0){0}else{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)] 
  Data$Interface <- TrialInfo$Interface[match(Data$trial_id_num, TrialInfo$trial_id_num)]
  Data$SubjectRole <- TrialInfo$SubjectRole[match(Data$trial_id_num, TrialInfo$trial_id_num)]

Add information about subjects.

  SubjectData <- read.csv("data/SubjectData.csv", stringsAsFactors=FALSE)
  Data$Group <- SubjectData$Group[match(Data$SubNum, SubjectData$SubNum)]

Exclude two participants who weren’t deceived.

  Data$Exclude <- SubjectData$Exclude[match(Data$SubNum, SubjectData$SubNum)] 
  Data <- Data [(Data$Exclude==0),]

Rename dataframe columns for consistency with other datasets.

  names(Data)[names(Data)=="trial_id_num"] <- "TrialID" 
  names(Data)[names(Data)=="INDEX"] <- "TrialNum" 

Remove columns that are now redundant.

Data <- subset(Data, select = -c(Exclude, 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", "Interface", "SubjectRole", "Group", "Condition")
Data[,FactorNames] <- colwise(as.factor)(Data[,FactorNames])
str(Data)
## 'data.frame':    10368 obs. of  16 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 ...

5 Respond condition

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

5.1 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

5.2 Saccadic 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

6 Initiating condition

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

6.1 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

6.2 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

6.3 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

7 Ratings data

RatingsData <- subset(SubjectData[(SubjectData$Exclude==0),], select = -c(Exclude, Sex, Age, EdHAnd, Block.1, Block.2, ET_Accuracy) )

RatLong <- melt(RatingsData, id=c("Group", "SubNum"), variable.name="Measure", value.name="Rating")

RatLong$Question <- as.factor(str_split_fixed(RatLong$Measure, "_", 2)[,1])
RatLong$Condition <- as.factor(str_split_fixed(RatLong$Measure, "_", 2)[,2])

RatLong$Question <- mapvalues(RatLong$Question, from = c("PreferVirtual", "PreferAlone", "Feel", "Behave", "Appear"), to = c("Prefer Virtual", "Prefer Alone", "Feel Humanlike", "Behave Humanlike", "Appear Humanlike"))

QuestionList <- c("Prefer Virtual", "Prefer Alone", "Behave Humanlike", "Appear Humanlike", "Feel Humanlike", "Cooperative", "Pleasant","Intuitive","Natural","Difficult")

RatingsBoxplot <- ggplot(RatLong[(RatLong$Condition!="Arrow"),], aes(x=factor(Question), y=Rating, fill=Group)) +
    geom_boxplot(outlier.colour=NA) +
    theme_minimal(base_size = 16, base_family = "") +
    xlim(QuestionList) +
    xlab("Question") +
    scale_y_continuous("Rating",limits=c(1,10), breaks=seq(1,10,1)) +
    scale_fill_brewer(palette="Blues")+
    coord_flip()
RatingsBoxplot

Wilcoxon.DF <- data.frame(Question=rev(QuestionList),
                          W=numeric(length(QuestionList)),
                          p=numeric(length(QuestionList)))

for (i in 1:length(QuestionList)) {
  QuestionData <- RatLong[(RatLong$Question==Wilcoxon.DF$Question[i]) & (RatLong$Condition=="Alan" | RatLong$Condition=="Both"),]
  Wilcox <- wilcox.test(Rating ~ Group, data = QuestionData, paired=FALSE, exact=FALSE)
  Wilcoxon.DF$W[i] <- Wilcox$statistic
  Wilcoxon.DF$p[i] <- Wilcox$p.value
}

kable(Wilcoxon.DF, row.names=FALSE, caption="Wilcoxon rank sum test with continuity correction", digits=3)
Wilcoxon rank sum test with continuity correction
Question W p
Difficult 346.5 0.209
Natural 216.5 0.139
Intuitive 384.5 0.045
Pleasant 188.5 0.038
Cooperative 193.0 0.039
Feel Humanlike 206.5 0.090
Appear Humanlike 284.5 0.950
Behave Humanlike 214.0 0.124
Prefer Alone 358.5 0.142
Prefer Virtual 299.5 0.817

8 Output

8.1 Tables

TTest.DF <- rbind(TTest.Dwell, TTest.Premi)
doc <- docx()
doc <- addTitle(doc, "Table 2: ", level=2)
doc <- addFlexTable( doc, vanilla.table(TTest.DF))
writeDoc(doc, file = "output/doc/Tables.docx")

8.2 Figures

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
png(filename="output/figures/Ratings.png", width=600, height=600)
RatingsBoxplot
dev.off()
## quartz_off_screen 
##                 2

```