PREPARATION

Load Packages

rm(list = ls())

library(pwr) # power analysis
library(afex) # package for running ANOVAs 
library(lsr) # useful functions 
library(MASS) # useful functions
library(lme4) # regressions 
library(xtable) # create latex tables 
library(png) # load pngs 
library(grid) # plotting
library(gridExtra) # plotting
library(stringr) # string handling  
library(magick) # image manipulation
library(Hmisc) # bootstrapping functions 
library(tidyverse) # dplyr, tidyr, ggplot2

Helper functions

# function to make percentages sum up to 100%
roundAndSum = function(x){
  rounded = floor(x)
  difference = 100-sum(rounded)
  decimals = x-floor(x)
  if (difference != 0){
    for (i in 1:difference){
      rounded[which.max(decimals)] = rounded[which.max(decimals)]+1
      decimals[which.max(decimals)] = NA
    }
  }
  return(rounded)
}

# root mean squared error 
rmse = function(x,y){
  return(sqrt(mean((x-y)^2)))
}

# judgment bar graphs 

judgmentPlot = function(df.plot,question,name){
  p = ggplot(df.plot,aes(x=index,y=value,fill=outcome,group=model))+
    geom_bar(stat = "identity",position = position_dodge(0.9), width=0.9, color = "black")+
    geom_linerange(aes(ymin = rating.cl.low, ymax = rating.cl.high),position = position_dodge(0.9), size = 1)+
    geom_text(data = df.plot[1:18,], aes(y=-14,x=index,label = as.character(clip)),
              size = 6,position = position_dodge(width = .7),hjust=0.5)+
    facet_grid(outcome.actual ~ outcome.counterfactual)+
    scale_fill_manual(values=c("red","green","black"))+
    theme_bw()+
    ylab(question)+
    coord_cartesian(xlim = c(0.5, 2.5),ylim=c(0,100))+
    theme_bw()+
    theme(legend.position="none",
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          panel.spacing=unit(c(.8),"cm"),
          plot.margin=unit(c(0.2,0.2,.5,.2),"cm"),
          axis.title.x = element_blank(),
          panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          strip.background = element_blank(),panel.border = element_rect(colour = "black"),
          text=element_text(size=20)
    )
  gt = ggplot_gtable(ggplot_build(p))
  gt$layout$clip = "off"
  grid.draw(gt)
  # print(p)
}

# print ANOVA results 

printAnova = function(fit,correction,row = 1){
  cat(
    "F(", fit$anova_table$`num Df`[row],",",fit$anova_table$`den Df`[row],") = ", fit$anova_table$F[row] %>% round(2),
    ", p = ", fit$anova_table$`Pr(>F)`[row] %>% round(3) * correction,
    ", \\eta_p^2 = ", fit$anova_table$pes[row] %>% round(2),
    sep = "")
}

DATA PROCESSING

Read in data file

load("../../data/summary/trackingDataFrames.RData")
# load("trackingDataFrames.RData")

Exclude participants based on track loss

  • Removes practice clips, and excludes participants based on track loss
  • Displays average track loss for each condition
df.trackloss = df.eyes %>% 
  group_by(condition,participant) %>% 
  summarise(n = n(),
            nas = sum(is.na(x) | is.na(y))
  ) %>% 
  mutate(freq = nas/n) %>% 
  arrange(freq) %>%
  mutate(index = 1:length(condition)) %>%
  group_by(condition) %>% 
  mutate(exclude = ifelse(index > 10,1,0)) %>% 
  ungroup

df.eyes = df.eyes %>% 
  merge(.,df.trackloss %>% select(condition,participant,exclude)) %>% 
  filter(exclude == 0) %>% 
  select(-exclude)

df.judgments = df.judgments %>% 
  merge(.,df.trackloss %>% select(condition,participant,exclude)) %>% 
  filter(exclude == 0) %>% 
  select(-exclude)

# remove practice trials 
df.eyes = df.eyes %>% 
  filter(clip < 19)

df.judgments = df.judgments %>% 
  filter(clip < 19)

# trackloss per condition
df.eyes %>% 
  group_by(condition,participant) %>% 
  summarise(trackloss = mean(is.na(x) | is.na(y))) %>% 
  group_by(condition) %>% 
  summarise(trackloss.mean = mean(trackloss),
            trackloss.sd = sd(trackloss)) %>% 
  mutate_at(vars(contains("track")),funs(round(.*100,2)))
## # A tibble: 3 x 3
##        condition trackloss.mean trackloss.sd
##           <fctr>          <dbl>        <dbl>
## 1 counterfactual           9.39         6.65
## 2         causal           7.31         3.49
## 3        outcome           8.27         3.35
rm("df.trackloss")

Demographics

  • Displays demographic information
df.demographics = df.judgments %>% 
  group_by(condition,participant) %>% 
  filter(row_number() == 1) %>% 
  ungroup() %>% 
  summarise(age.mean = round(mean(age),2),
            age.sd = round(sd(age),2),
            nfemale = sum(sex == "female"),
            time.mean = round(mean(time),2),
            time.sd = round(sd(time),2)) %>% 
  print()
## # A tibble: 1 x 5
##   age.mean age.sd nfemale time.mean time.sd
##      <dbl>  <dbl>   <int>     <dbl>   <dbl>
## 1    36.23  15.85       9     13.58    0.91

ANALYSES & PLOTS

Behavioral judgments: Counterfactual condition

  • ANOVA results and group means
# ANOVA
fit = aov_ez(id = "participant",
             dv = "rating",
             data = df.judgments %>% filter(condition == "counterfactual"),
             between = NULL, 
             anova_table = list(es = "pes",correction="none"),
             within = c("outcome.actual","outcome.counterfactual"))
fit$anova_table
## Anova Table (Type 3 tests)
## 
## Response: rating
##                                       num Df den Df    MSE       F     pes
## outcome.actual                             2     18 137.47  2.7428 0.23357
## outcome.counterfactual                     2     18 506.27 79.5097 0.89832
## outcome.actual:outcome.counterfactual      4     36 139.13  4.3615 0.32642
##                                          Pr(>F)    
## outcome.actual                         0.091254 .  
## outcome.counterfactual                1.162e-09 ***
## outcome.actual:outcome.counterfactual  0.005601 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
printAnova(fit,correction = 1, row = 1)
## F(2,18) = 2.74, p = 0.091, \eta_p^2 = 0.23
printAnova(fit,correction = 1, row = 2)
## F(2,18) = 79.51, p = 0, \eta_p^2 = 0.9
# Means and SDs 
df.judgments %>% filter(condition == "counterfactual") %>% 
  group_by(outcome.counterfactual) %>% 
  summarise(mean = mean(rating),
            sd = sd(rating))
## # A tibble: 3 x 3
##   outcome.counterfactual     mean       sd
##                   <fctr>    <dbl>    <dbl>
## 1                   miss 14.70904 18.25289
## 2                  close 26.77709 23.04242
## 3                    hit 83.32196 24.38691

Plot: Counterfactual condition

  • Plots the model predictions and results of the counterfactual condition
  • Shows r and RMSE of the model
# read noisy simulation results
load("../../data/summary/simple_stepNoise_predictions.RData")

df.plot = df.judgments %>% 
  filter(condition == "counterfactual") %>% 
  group_by(condition,clip,outcome,outcome.actual,outcome.counterfactual,index) %>%
  summarise(rating.mean = smean.cl.boot(rating)[1],
            rating.cl.low = smean.cl.boot(rating)[2],
            rating.cl.high = smean.cl.boot(rating)[3]) %>%
  ungroup %>% 
  mutate(outcome.actual = factor(outcome.actual,labels = paste("actual",c("miss","close","hit"))),
         outcome.counterfactual = factor(outcome.counterfactual,labels = paste("counterfactual",c("miss","close","hit"))))

# find noise values that predicts participants' counterfactual judgments best
tmp = numeric(length=length(simple_noise_counterfactual.means))
for (i in 1:length(simple_noise_counterfactual.means)){
  tmp[i] = cor(simple_noise_counterfactual.means[[i]],df.plot$rating.mean)
}

#counterfactual outcomes
df.plot$outcome = c(0, 0, 0, 1, 1, 1,
                    0, 0, 1, 0, 1, 1,
                    0, 0, 1, 0, 1, 1)
df.plot$outcome = as.factor(df.plot$outcome)

df.plot$model = simple_noise_counterfactual.means[[which.max(tmp)]]
df.plot$model = as.numeric(lm(rating.mean~model,data = df.plot)$fitted.values)

df.plot = df.plot %>% 
  gather(model,value,c(rating.mean,model)) %>% 
  mutate(outcome = as.factor(ifelse(model == "model",3,outcome)),
         model = factor(model,levels = c("rating.mean","model"))) %>% 
  mutate_at(vars(contains("cl.")),funs(ifelse(model == "model",NA,.)))

#plot 
judgmentPlot(df.plot,expression(paste('Agreement with ', bold(counterfactual), ' statement')),'Figure 2a')

# r and RMSE 
cor(df.plot$value[df.plot$model == "rating.mean"],
    df.plot$value[df.plot$model == "model"]) %>% round(2)
## [1] 0.92
rmse(df.plot$value[df.plot$model == "rating.mean"],
     df.plot$value[df.plot$model == "model"]) %>% round(2)
## [1] 12.57

Behavioral judgments: Causal condition

  • ANOVA results and group means (separate for cause and prevention judgments)
# ANOVA on prevention judgments
fit = aov_ez(id = "participant",
             dv = "rating",
             data = df.judgments %>% filter(condition == "causal", outcome == 0),
             between = NULL, 
             anova_table = list(es = "pes",correction="none"),
             within = c("outcome.counterfactual","outcome.actual"))
printAnova(fit,correction = 1, row = 1)
## F(2,18) = 21.86, p = 0, \eta_p^2 = 0.71
printAnova(fit,correction = 1, row = 2)
## F(1,9) = 0.02, p = 0.881, \eta_p^2 = 0
# Means and SDs 
df.judgments %>% 
  filter(condition == "causal",outcome == 0) %>% 
  group_by(outcome.counterfactual) %>% 
  summarise(mean = mean(rating),
            sd = sd(rating))
## # A tibble: 3 x 3
##   outcome.counterfactual     mean       sd
##                   <fctr>    <dbl>    <dbl>
## 1                   miss 19.37836 28.29847
## 2                  close 48.18286 35.12835
## 3                    hit 86.33149 24.94903
# ANOVA on causation judgments
fit = aov_ez(id = "participant",
             dv = "rating",
             data = df.judgments %>% filter(condition == "causal", outcome == 1),
             between = NULL, 
             anova_table = list(es = "pes",correction="none"),
             within = c("outcome.counterfactual","outcome.actual"))
fit$anova_table
## Anova Table (Type 3 tests)
## 
## Response: rating
##                                       num Df den Df     MSE       F
## outcome.counterfactual                     2     18 1166.54 15.7817
## outcome.actual                             1      9  245.50  0.2059
## outcome.counterfactual:outcome.actual      2     18  195.66  0.6097
##                                           pes    Pr(>F)    
## outcome.counterfactual                0.63683 0.0001099 ***
## outcome.actual                        0.02237 0.6607451    
## outcome.counterfactual:outcome.actual 0.06344 0.5543753    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
printAnova(fit,correction = 1, row = 1)
## F(2,18) = 15.78, p = 0, \eta_p^2 = 0.64
printAnova(fit,correction = 1, row = 2)
## F(1,9) = 0.21, p = 0.661, \eta_p^2 = 0.02
# Means and SDs 
df.judgments %>% 
  filter(condition == "causal",outcome == 1) %>% 
  group_by(outcome.counterfactual) %>% 
  summarise(mean = mean(rating),
            sd = sd(rating))
## # A tibble: 3 x 3
##   outcome.counterfactual     mean       sd
##                   <fctr>    <dbl>    <dbl>
## 1                   miss 91.31172 16.28326
## 2                  close 72.82504 29.10346
## 3                    hit 31.64303 40.12222

Plot: Causal condition

  • Plots the model predictions and results of the causal condition, shows r and RMSE of the model
df.data = df.judgments %>% 
  filter(condition == "causal") %>% 
  group_by(clip,outcome,outcome.actual,outcome.counterfactual,index) %>%
  summarise(rating.mean = smean.cl.boot(rating)[1],
            rating.cl.low = smean.cl.boot(rating)[2],
            rating.cl.high = smean.cl.boot(rating)[3]) %>%
  ungroup

# use participants' counterfactual judgments as prediction
df.model = df.judgments %>% 
  filter(condition == "counterfactual") %>% 
  group_by(clip,outcome,outcome.actual,outcome.counterfactual,index) %>%
  summarise(model = mean(rating)) %>% 
  mutate(model = ifelse(outcome == 0,model,100-model)) %>% 
  ungroup

df.plot = merge(df.data,df.model)

regression = function(x){
  tmp = as.numeric(lm(rating.mean~model,data=x)$fitted.values)
  return(tmp)
}

df.plot = df.plot %>%
  mutate(model = regression(.),
         outcome.actual = factor(outcome.actual,labels = paste("actual",c("miss","close","hit"))),
         outcome.counterfactual = factor(outcome.counterfactual,labels = paste("counterfactual",c("miss","close","hit"))))

df.plot = df.plot %>% 
  gather(model,value,c(rating.mean,model)) %>% 
  mutate(outcome = as.factor(ifelse(model == "model",3,outcome)),
         model = factor(model,levels = c("rating.mean","model"))) %>% 
  mutate_at(vars(contains("cl.")),funs(ifelse(model == "model",NA,.)))

#plot 
judgmentPlot(df.plot,expression(paste('Agreement with ', bold(causal), ' statement')), 'Figure 2b')

# r and RMSE 
cor(df.plot$value[df.plot$model == "rating.mean"],
    df.plot$value[df.plot$model == "model"]) %>% round(2)
## [1] 0.92
rmse(df.plot$value[df.plot$model == "rating.mean"],
     df.plot$value[df.plot$model == "model"]) %>% round(2)
## [1] 10.8

Behavioral judgments: Outcome condition

  • ANOVA results and group means (separate for positive and negative outcomes)
# ANOVA 
fit = aov_ez(id = "participant",
             dv = "rating",
             data = df.judgments %>% filter(condition == "outcome", outcome == 0),
             between = NULL, 
             anova_table = list(es = "pes",correction="none"),
             within = c("outcome.actual","outcome.counterfactual"))
fit$anova_table
## Anova Table (Type 3 tests)
## 
## Response: rating
##                                       num Df den Df     MSE       F
## outcome.actual                             1      9 1270.14 40.8861
## outcome.counterfactual                     2     18  246.37  0.2571
## outcome.actual:outcome.counterfactual      2     18  138.67  9.0754
##                                           pes    Pr(>F)    
## outcome.actual                        0.81959 0.0001261 ***
## outcome.counterfactual                0.02778 0.7760542    
## outcome.actual:outcome.counterfactual 0.50209 0.0018810 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
printAnova(fit,correction = 1, row = 1)
## F(1,9) = 40.89, p = 0, \eta_p^2 = 0.82
printAnova(fit,correction = 1, row = 2)
## F(2,18) = 0.26, p = 0.776, \eta_p^2 = 0.03
# Means and SDs 
df.judgments %>% 
  filter(condition == "outcome",outcome == 0) %>% 
  group_by(outcome.actual) %>% 
  summarise(mean = mean(rating),
            sd = sd(rating))
## # A tibble: 2 x 3
##   outcome.actual     mean       sd
##           <fctr>    <dbl>    <dbl>
## 1           miss 87.95569 17.40096
## 2          close 29.11643 34.34083
# ANOVA 
fit = aov_ez(id = "participant",
             dv = "rating",
             data = df.judgments %>% filter(condition == "outcome", outcome == 1),
             between = NULL, 
             anova_table = list(es = "pes",correction="none"),
             within = c("outcome.actual","outcome.counterfactual"))
fit$anova_table
## Anova Table (Type 3 tests)
## 
## Response: rating
##                                       num Df den Df    MSE       F     pes
## outcome.actual                             1      9 784.88 17.9201 0.66568
## outcome.counterfactual                     2     18 385.69  7.5283 0.45548
## outcome.actual:outcome.counterfactual      2     18 230.12  0.0981 0.01078
##                                         Pr(>F)   
## outcome.actual                        0.002196 **
## outcome.counterfactual                0.004209 **
## outcome.actual:outcome.counterfactual 0.907080   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
printAnova(fit,correction = 1, row = 1)
## F(1,9) = 17.92, p = 0.002, \eta_p^2 = 0.67
printAnova(fit,correction = 1, row = 2)
## F(2,18) = 7.53, p = 0.004, \eta_p^2 = 0.46
df.judgments %>% 
  filter(condition == "outcome",outcome == 1) %>% 
  group_by(outcome.actual) %>% 
  summarise(mean = mean(rating),
            sd = sd(rating))
## # A tibble: 2 x 3
##   outcome.actual     mean       sd
##           <fctr>    <dbl>    <dbl>
## 1          close 48.95892 32.48492
## 2            hit 79.58041 27.09174

Plot: Outcome condition

  • Plots the model predictions and results of the outcome condition, shows r and RMSE of the model
# calculate distance from center of the goal 
df.model = df.eyes %>% 
  filter(participant == 1, instance == 1) %>% 
  group_by(clip) %>%
  filter(frame == t.outcome.actual) %>%
  mutate(distance = abs(both.B.y-384)) %>% 
  ungroup() %>% 
  select(clip,distance)

regression = function(x){
  tmp = as.numeric(lm(rating.mean~distance*outcome,data=x)$fitted.values)
  return(tmp)
}

df.plot = df.judgments %>% 
  filter(condition == "outcome") %>% 
  group_by(clip,outcome,outcome.actual,outcome.counterfactual,index) %>%
  summarise(rating.mean = smean.cl.boot(rating)[1],
            rating.cl.low = smean.cl.boot(rating)[2],
            rating.cl.high = smean.cl.boot(rating)[3]) %>%
  ungroup %>% 
  left_join(df.model)

lm(rating.mean~distance*outcome,data=df.plot) %>% summary()
## 
## Call:
## lm(formula = rating.mean ~ distance * outcome, data = df.plot)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.249  -3.418   0.947   6.867  20.501 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -3.83171   13.38057  -0.286 0.778794    
## distance          0.34864    0.06023   5.788 4.70e-05 ***
## outcome          87.85603   15.05820   5.834 4.33e-05 ***
## distance:outcome -0.96765    0.21615  -4.477 0.000522 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.56 on 14 degrees of freedom
## Multiple R-squared:  0.7519, Adjusted R-squared:  0.6987 
## F-statistic: 14.14 on 3 and 14 DF,  p-value: 0.000161
df.plot = df.plot %>%
  mutate(model = regression(.),
         outcome.actual = factor(outcome.actual,labels = paste("actual",c("miss","close","hit"))),
         outcome.counterfactual = factor(outcome.counterfactual,labels = paste("counterfactual",c("miss","close","hit"))))

df.plot = df.plot %>% 
  gather(model,value,c(rating.mean,model)) %>% 
  mutate(outcome = as.factor(ifelse(model == "model",3,outcome)),
         model = factor(model,levels = c("rating.mean","model"))) %>% 
  mutate_at(vars(contains("cl.")),funs(ifelse(model == "model",NA,.)))

#plot 
judgmentPlot(df.plot,expression(paste('Agreement with ', bold(outcome), ' statement')), 'Figure 2c')

# r and RMSE 
cor(df.plot$value[df.plot$model == "rating.mean"],
    df.plot$value[df.plot$model == "model"]) %>% round(2)
## [1] 0.87
rmse(df.plot$value[df.plot$model == "rating.mean"],
     df.plot$value[df.plot$model == "model"]) %>% round(2)
## [1] 12.84

Eye-tracking: Static analysis

  • Categorizes participants’ looks prior to the two balls colliding into counterfactual looks and other looks
  • Displays results table separated by condition, and separated by condition x clip
  • Reports statistical tests on whether looks differ between conditions
threshold.distance = 100 #distance in pixels from the counterfactual path 
threshold.x = 50 #minimal distance from collision point

df.counterfactual.path = df.eyes %>% 
  filter(participant == 1, instance == 1, condition == "causal") %>% 
  group_by(clip) %>% 
  filter(frame == t.collision | frame == t.outcome.counterfactual) %>% 
  ungroup() %>% 
  arrange(clip,frame) %>% 
  mutate(path.marker = rep(c("start","end"),(nrow(.)/2))) %>% 
  select(clip,onlyB.x,onlyB.y,path.marker)

df.counterfactual.path = cbind(df.counterfactual.path %>% filter(path.marker == "start") %>% select(-path.marker),
                               df.counterfactual.path %>% filter(path.marker == "end") %>% select(-path.marker,-clip)) %>% 
  setNames(c("clip","path.start.x","path.start.y","path.end.x","path.end.y"))

df.eyes = df.eyes %>% 
  left_join(df.counterfactual.path)

# distance of each fixation to the counterfactual path 
fun.distance_to_line = function(x,y,path.start.x,path.start.y,path.end.x,path.end.y){
  distance = abs((path.end.y-path.start.y)*x-(path.end.x-path.start.x)*y+path.end.x*path.start.y-path.end.y*path.start.x)/
    sqrt((path.end.y-path.start.y)^2+(path.end.x-path.start.x)^2) #distance of point to line 
  return(distance)
}

df.eyes = df.eyes %>% 
  rowwise %>% 
  mutate(counterfactual.distance = fun.distance_to_line(x,y,path.start.x,path.start.y,path.end.x,path.end.y)) %>% 
  ungroup() %>% 
  mutate(counterfactual.look = ifelse(counterfactual.distance < threshold.distance & 
                                        !is.na(counterfactual.distance) & 
                                        x < (path.start.x-threshold.x), "counterfactual", "other")) %>% 
  mutate(counterfactual.look = factor(counterfactual.look,levels = c("counterfactual","other"))) %>% 
  arrange(condition,participant,clip,instance,frame)

# result table 
df.tmp = df.eyes %>% 
  na.omit() %>% 
  group_by(clip) %>% 
  filter(saccade == 'sacc.end') %>% 
  filter(frame > 15,frame < t.collision) %>% 
  group_by(condition,participant,counterfactual.look) %>%
  summarise (n = n()) %>%
  mutate(freq = n / sum(n),
         freq_rounded = roundAndSum(freq*100)) %>% 
  group_by(condition,counterfactual.look) %>% 
  summarise(mfreq = mean(freq),sdfreq = round(sd(freq)*100,2)) %>% 
  mutate(mfreq = roundAndSum(mfreq*100)) %>% 
  filter(counterfactual.look == "counterfactual") %>% 
  ungroup %>% 
  mutate(data = paste0(mfreq,"% (", sdfreq,")")) %>% 
  select(condition,data) %>%
  spread(condition,data) %>% 
  print()
## # A tibble: 1 x 3
##   counterfactual      causal   outcome
## *          <chr>       <chr>     <chr>
## 1     22% (9.75) 23% (11.24) 6% (4.78)
# results per clip
df.tmp = df.eyes %>%
  na.omit() %>%
  filter(saccade == "sacc.end", frame < t.collision, frame > 15) %>%
  group_by(clip,condition,participant,counterfactual.look) %>%
  summarise (n = n()) %>%
  mutate(freq = n / sum(n),
         freq_rounded = roundAndSum(freq*100)) %>%
  complete(counterfactual.look,fill=list(n=0,freq=0,freq_rounded=0)) %>% 
  group_by(condition,clip,counterfactual.look) %>% 
  summarise(mfreq = mean(freq),sdfreq = round(sd(freq)*100,2)) %>% 
  ungroup() %>% 
  mutate(mfreq = roundAndSum(mfreq*100)) %>% 
  filter(counterfactual.look == "counterfactual") %>% 
  ungroup %>% 
  mutate(data = paste0(mfreq,"% (", sdfreq,")")) %>% 
  select(condition,clip,data) %>%
  spread(condition,data) %>% 
  print()
## # A tibble: 18 x 4
##     clip counterfactual      causal     outcome
##  * <int>          <chr>       <chr>       <chr>
##  1     1    12% (13.73) 28% (20.92)   2% (3.83)
##  2     2    17% (18.62) 24% (12.69)     5% (11)
##  3     3    29% (26.29) 25% (21.03)  9% (10.38)
##  4     4    36% (19.93) 31% (19.84)   1% (2.87)
##  5     5    37% (17.18)  20% (18.5)  10% (12.2)
##  6     6     10% (9.96)   10% (9.2)   6% (5.21)
##  7     7    17% (20.88) 20% (18.51)      1% (0)
##  8     8    19% (15.76) 20% (16.73)   3% (4.22)
##  9     9    28% (15.33) 26% (16.92)      1% (0)
## 10    10    17% (16.71) 28% (11.51)   5% (7.85)
## 11    11    45% (19.62) 26% (20.85)  7% (11.74)
## 12    12    28% (22.12) 33% (20.96)   7% (7.29)
## 13    13     10% (8.55)  26% (24.9)   4% (4.01)
## 14    14    22% (23.83) 28% (26.62)  7% (14.13)
## 15    15    20% (16.36)  26% (13.3)   7% (8.58)
## 16    16    25% (13.54) 20% (14.29)  5% (10.84)
## 17    17     16% (9.63)  15% (8.85) 11% (10.92)
## 18    18    15% (10.74) 15% (12.98)   6% (5.55)
# Chi-square test of types of looks per condition 

df.chi = df.eyes %>% 
  filter(saccade == "sacc.end") %>% 
  filter(frame < t.collision, frame > 15)
  
df.chi = table(df.chi$counterfactual.look,df.chi$condition)
chisq.test(df.chi)
## 
##  Pearson's Chi-squared test
## 
## data:  df.chi
## X-squared = 302.49, df = 2, p-value < 2.2e-16

Plot: Fixation map prior to collision

  • Plots fixation maps for clip 4 in the causal condition (code can be adapted to show plots for other conditions and clips)
# condition.name = "counterfactual"
condition.name = "causal"
# condition.name = "outcome"

# clips = 1:18
clips = 4

df.plot = df.eyes %>% 
  filter(saccade == "sacc.end") %>%
  filter(condition == condition.name) %>% 
  filter(frame < t.collision, frame > 15) %>% 
  na.omit()

df.legend_labels = df.plot %>% 
  group_by(clip,counterfactual.look) %>% 
  summarise (n = n()) %>%
  complete(counterfactual.look,fill = list(n=0)) %>% 
  group_by(clip) %>% 
  mutate(freq = n / sum(n),
         freq_rounded = roundAndSum(freq*100)) %>% 
  ungroup() %>% 
  mutate(counterfactual.look = factor(counterfactual.look, levels = c("counterfactual", "other"), 
                                      labels = c('Counterfactual Saccade', 'Other Saccade')),
         counterfactual.look = as.character(counterfactual.look),
         label = paste0(counterfactual.look," (",freq_rounded,"%)")) %>%
  ungroup

for (cli in clips){
  ins = c(1,2)
  df = df.plot %>% 
    filter(clip == cli,
           instance %in% ins)
  legend_label = df.legend_labels %>% filter(clip == cli) %>% pull(label)
  
  m = readPNG(paste0("../../figures/diagrams/png/clip_",cli,".png"), FALSE)
  w = matrix(rgb(m[,,1],m[,,2],m[,,3]), nrow=dim(m)[1])
  
p = ggplot(df,aes(x=x,y=y))+
    annotation_custom(xmin=-Inf, ymin=-Inf, xmax=Inf, ymax=Inf,
                      rasterGrob(w)) +
    geom_point(data = df,size=5, alpha = 0.6, aes(color = counterfactual.look, shape = counterfactual.look, fill = counterfactual.look),stroke=2) +
    scale_shape_manual("",values = c(21,16),
                       labels = legend_label, limits = c("counterfactual","other")
                         )+
    scale_color_manual("",values = c('black','gray20'),
                       labels = legend_label, limits = c("counterfactual","other"))+
    scale_fill_manual("",values = c('white','gray20'),
                      labels = legend_label, limits = c("counterfactual","other"))+
    expand_limits(x=c(0,1024),y=c(0,768))+
    scale_y_continuous(expand = c(0,0)) + 
    scale_x_continuous(expand = c(0,0)) +
    labs(x=NULL, y=NULL)+
    theme_void()+
    theme(
      text = element_text(size=30),
      axis.ticks.length = unit(0,"null"),
      legend.position = c(0.5,0), 
      legend.justification = c(0.5,-0.2),
      # legend.position = c(0.5,1), 
      # legend.justification = c(0.5,1.2),
      legend.title = element_blank(),
      legend.key = element_rect(fill=alpha('white', 0)),
      legend.background = element_rect(fill='white',color="black"),
      legend.key.height = unit(1.2,"cm"),
      legend.key.width = unit(1,"cm"),
      legend.margin = margin(t=0.1,r=0.2,b=0.1,l=0.1,unit="cm"),
      legend.direction = 'vertical'
    )
}
print(p)

Eye-tracking: Analysis of HMM results

  • Runs ANOVA on the probability that participants made a counterfactual look to where ball B would have gone
  • Reports the ANOVA results, as well as the results of post-hoc comparisons between conditions
# Mixed ANOVAs

# variable = "post.A.look"
# variable = "post.B.look"
# variable = "post.A.predict"
# variable = "post.B.predict"
# variable = "post.A.counterfactual"
variable = "post.B.counterfactual"
# variable = "post.other"

df.analysis = df.eyes %>% 
  filter(saccade == "sacc.end") %>% 
  filter(frame < t.outcome.actual, frame > 15) %>% 
  select_("participant","condition",variable) %>% 
  group_by(participant,condition) %>%
  summarize_(value = paste0("mean(",variable,",na.rm=T)")) %>% 
  ungroup()

fit = aov_ez(id = "participant",
             dv = "value",
             data = df.analysis,
             between = "condition", 
             anova_table = list(es = "pes",correction="none"),
             within = NULL)
summary(fit)
## Anova Table (Type 3 tests)
## 
## Response: value
##           num Df den Df      MSE      F     pes    Pr(>F)    
## condition      2     27 0.002543 17.632 0.56636 1.262e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#post-hoc tests 
tmp = lsmeans(fit, "condition", contr = "pairwise") 
ttest = summary(tmp,adjust="none")$contrasts

# effect sizes 
cohen = c(cohensD(df.analysis %>% filter(condition == "counterfactual") %>% select(value) %>% unlist %>% as.numeric,
                  df.analysis %>% filter(condition == "causal") %>% select(value) %>% unlist %>% as.numeric),
            cohensD(df.analysis %>% filter(condition == "counterfactual") %>% select(value) %>% unlist %>% as.numeric,
                    df.analysis %>% filter(condition == "outcome") %>% select(value) %>% unlist %>% as.numeric),
            cohensD(df.analysis %>% filter(condition == "causal") %>% select(value) %>% unlist %>% as.numeric,
                    df.analysis %>% filter(condition == "outcome") %>% select(value) %>% unlist %>% as.numeric))


# print results 
correction = 7 # bonferroni correction for multiple comparisons 

cat(
  "F(", fit$anova_table$`num Df`,",",fit$anova_table$`den Df`,") = ", fit$anova_table$F %>% round(2),
  ", p = ", fit$anova_table$`Pr(>F)` %>% round(3),
  ", \\eta_p^2 = ", fit$anova_table$pes %>% round(2),
  sep = "")
## F(2,27) = 17.63, p = 0, \eta_p^2 = 0.57
for (i in 1:nrow(ttest)){
  cat("\n")
  cat(ttest$contrast[i] %>% as.character,": ",
    "t(",ttest$df[i] ,") = ", ttest$t.ratio[i] %>% round(2),
  ", p = ", ttest$p.value[i] %>% round(3),
  ", d = ", cohen[i] %>% round(2),
  sep = "")
}
## 
## counterfactual - causal: t(27) = 0.39, p = 0.7, d = 0.15
## counterfactual - outcome: t(27) = 5.33, p = 0, d = 2.85
## causal - outcome: t(27) = 4.94, p = 0, d = 2.26

Plot: HMM results

  • Plots the average probability of participants being in different states separated by condition (only taking into account end points of participants’ saccades)
df.plot = df.eyes %>% 
  filter(saccade == "sacc.end") %>%
  filter(frame < t.outcome.actual, frame > 15) %>%
  group_by(participant,condition) %>% 
  select(contains("post")) %>% 
  summarise_all(funs(mean(.,na.rm=T))) %>% 
  gather(look,percentage,-c(condition,participant)) %>% 
  ungroup %>%
  mutate(look = factor(look,levels = c("post.A.look", "post.B.look", "post.A.predict", "post.B.predict",
                                       "post.A.counterfactual", "post.B.counterfactual","post.other"),
                       labels=c("A look", "B look", "A predict ", "B predict ", "A counterfactual ",
                                "B counterfactual ", "other")),
         condition = factor(condition,levels = c("counterfactual","causal","outcome"),
                            labels=c(expression("Counterfactual\ncondition"),
                                     expression("Causal\ncondition"),
                                     expression("Outcome\ncondition")
                            ))) %>% 
  arrange(condition,participant,look)

ggplot(df.plot,aes(x=look,y=percentage,fill = look))+
  stat_summary(fun.y = mean, geom = "bar", color="black", 
               position = position_dodge(0.8), width=0.8)+
  stat_summary(fun.data = mean_cl_boot, geom = "linerange",size = 1, 
               position = position_dodge(0.8))+
  facet_grid(~condition)+
  labs(y = 'probability of each type of look', fill = '')+
  scale_y_continuous(breaks = seq(0,0.4,0.1),labels = paste0(seq(0,40,10),"%"),
                     expand=c(0,0))+
  scale_fill_brewer(type = 'qual', palette = 3)+
  coord_cartesian(ylim=c(0,0.45))+
  theme_bw()+
  theme(text = element_text(size=30),
        legend.position = "bottom",
        panel.grid.major.x = element_blank(), panel.grid.minor = element_blank(),
        panel.grid.major.y = element_line(color = "gray60",linetype=2),
        strip.background = element_blank(),panel.border = element_rect(colour = "black"),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        strip.text = element_text(size=36),
        legend.text = element_text(size=30),
        panel.spacing.y=unit(c(.8),"cm"),
        plot.margin=unit(c(0.2,0.2,.2,.2),"cm"),
        axis.title.y = element_text(margin = margin(0,0.5,0,0,unit="cm")),
        legend.key = element_blank(),
        legend.key.size = unit(1.2,"cm"),
        legend.key.height = unit(1.2,"cm")
  )+
  guides(fill=guide_legend(nrow=2,byrow=TRUE))

ggsave('test.pdf',width=20,height=8)

Eyes & Judgments: Relationship between counterfactual looks and outcome certainty

  • Checks whether there is a relationship between participants’ eye-movements and their uncertainty in the causal judgment
# check whether participants make more counterfactual looks for the clips in which the 
# counterfactual outcome was close 
df.tmp = df.eyes %>%
  na.omit() %>%
  filter(condition == "causal") %>% 
  filter(saccade == "sacc.end") %>%
  group_by(participant,outcome.counterfactual) %>%
  summarise_at(vars(contains("post")),funs(mean)) %>% 
  select(participant,outcome.counterfactual,post.B.counterfactual) %>% 
  ungroup()

fit = aov_ez(id = "participant",
             dv = "post.B.counterfactual",
             data = df.tmp,
             between = NULL, 
             anova_table = list(es = "pes",correction="none"),
             within = c("outcome.counterfactual"))
fit$anova_table
## Anova Table (Type 3 tests)
## 
## Response: post.B.counterfactual
##                        num Df den Df        MSE     F     pes   Pr(>F)   
## outcome.counterfactual      2     18 0.00087717 9.503 0.51359 0.001524 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tmp = lsmeans(fit, "outcome.counterfactual", contr = "pairwise")
ttest = summary(tmp,adjust="none")$contrasts

cohen = c(cohensD(df.tmp %>% filter(outcome.counterfactual == "miss") %>% select(post.B.counterfactual) %>% unlist %>% as.numeric,
                  df.tmp %>% filter(outcome.counterfactual == "close") %>% select(post.B.counterfactual) %>% unlist %>% as.numeric),
          cohensD(df.tmp %>% filter(outcome.counterfactual == "miss") %>% select(post.B.counterfactual) %>% unlist %>% as.numeric,
                  df.tmp %>% filter(outcome.counterfactual == "hit") %>% select(post.B.counterfactual) %>% unlist %>% as.numeric),
          cohensD(df.tmp %>% filter(outcome.counterfactual == "close") %>% select(post.B.counterfactual) %>% unlist %>% as.numeric,
                  df.tmp %>% filter(outcome.counterfactual == "hit") %>% select(post.B.counterfactual) %>% unlist %>% as.numeric))

printAnova(fit,1)
## F(2,18) = 9.5, p = 0.002, \eta_p^2 = 0.51
for (i in 1:nrow(ttest)){
  cat("\n")
  cat(ttest$contrast[i] %>% as.character,": ",
      "t(",ttest$df[i] ,") = ", ttest$t.ratio[i] %>% round(2),
      ", p = ", ttest$p.value[i] %>% round(3),
      ", d = ", cohen[i] %>% round(2),
      sep = "")
}
## 
## miss - close: t(18) = -3.37, p = 0.003, d = 0.82
## miss - hit: t(18) = 0.71, p = 0.49, d = 0.21
## close - hit: t(18) = 4.08, p = 0.001, d = 1.05
# Plot 
df.plot = df.eyes %>% 
  filter(condition == "causal") %>%
  filter(saccade == "sacc.end") %>%
  filter(frame < t.outcome.actual, frame > 15) %>%
  group_by(participant,clip,outcome) %>%
  mutate(rating = abs(50-rating)) %>%
  summarize(counterfactual = mean(post.B.counterfactual,na.rm=T),rating = mean(rating))

ggplot(df.plot,aes(x=counterfactual,y=rating))+
  geom_smooth(method = lm,color = "black")+
  geom_point()+
  labs(y = "certainty in causal jugment", x = "p(look = B counterfactual)")+
  theme_bw()+
  theme(panel.grid = element_blank(),
        text = element_text(size = 24),
        axis.text.y = element_text(angle = 0),
        legend.position = "none"
  )

cor.test(df.plot$counterfactual,df.plot$rating,test = T)
## 
##  Pearson's product-moment correlation
## 
## data:  df.plot$counterfactual and df.plot$rating
## t = -4.3032, df = 178, p-value = 2.774e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4337526 -0.1682546
## sample estimates:
##        cor 
## -0.3069637

Eyes & Judgments: Relationship between counterfactual looks and model predictions

  • Checks whether there is a relationship between participants’ eye-movements and how well the counterfactual simulation model explains their judgments
# model predictions based on participants' counterfactual judgments 
df.model = df.judgments %>% 
  filter(condition == "counterfactual") %>% 
  select(outcome,participant,clip,rating) %>% 
  group_by(outcome,clip) %>% 
  summarise(model = mean(rating)) %>% 
  mutate(model = ifelse(outcome == 0,model,100-model))

# merge counterfactual looks and model predictions, aggregate per participant  
df.plot = df.eyes %>% 
  filter(condition == "causal") %>%
  filter(saccade == "sacc.end") %>%
  filter(frame < t.outcome.actual, frame > 15) %>%
  group_by(participant,clip,outcome) %>%
  summarize(counterfactual = mean(post.B.counterfactual,na.rm=T),rating = mean(rating)) %>% 
  merge(.,df.model) %>% 
  group_by(participant) %>% 
  summarize(counterfactual = mean(counterfactual), fit = cor(rating,model)) %>% 
  arrange(counterfactual) %>% 
  print()
## # A tibble: 10 x 3
##    participant counterfactual        fit
##          <dbl>          <dbl>      <dbl>
##  1          49     0.05123800 -0.3092899
##  2          40     0.08786687  0.6908162
##  3          50     0.10901883  0.5668944
##  4          11     0.11267575  0.2020171
##  5           7     0.16560696  0.9658276
##  6           8     0.18987878  0.7729559
##  7           5     0.20302528  0.8470754
##  8          41     0.21292387  0.7891971
##  9           4     0.22678369  0.8675713
## 10           1     0.23295458  0.8757335
# plot 
ggplot(df.plot,aes(x=counterfactual,y=fit))+
  geom_smooth(method = lm,color = "black")+
  geom_point(size=4)+
  labs(y = "model fit", x = "p(look = B counterfactual)")+
  theme_bw()+
  theme(panel.grid = element_blank(),
        text = element_text(size = 24),
        axis.text.y = element_text(angle = 0)
  )

cor.test(df.plot$counterfactual,df.plot$fit)
## 
##  Pearson's product-moment correlation
## 
## data:  df.plot$counterfactual and df.plot$fit
## t = 3.5849, df = 8, p-value = 0.007137
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3073868 0.9467283
## sample estimates:
##       cor 
## 0.7850713
rmse(df.plot$counterfactual,df.plot$fit)
## [1] 0.571153

SUPPLEMENTARY PLOTS

Change in r and RMSE for counterfactual judgments

  • Plots the change in r and RMSE as a function of how much noise is added to ball B’s motion vector in the counterfactual simulation
df.plot = df.judgments %>% 
  filter(condition == "counterfactual") %>% 
  select(condition,outcome,participant,clip,rating) %>% 
  group_by(clip,outcome) %>% 
  summarise(rating.mean = mean(rating)) %>% 
  arrange(clip) %>% 
  ungroup

# find noise values that predicts participants' counterfactual judgments best
tmp = as.data.frame(matrix(nrow= length(simple_noise_counterfactual.means),ncol = 3))
colnames(tmp) = c("noise", "r", "RMSE")

for (i in 1:length(simple_noise_counterfactual.means)){
  tmp$noise[i] = names(simple_noise_counterfactual.means)[i]
  df = data.frame(dv = df.plot$rating.mean, iv = simple_noise_counterfactual.means[[i]]*100)
  tmp$r[i] = cor(simple_noise_counterfactual.means[[i]],df.plot$rating.mean)
  tmp$RMSE[i] = sqrt(mean((lm(dv~iv,data=df)$fitted.values-df.plot$rating.mean)^2))
}

df.plot = tmp %>% 
  mutate(noise = as.numeric(str_replace(noise,"noise","")))

# Correlation 
ggplot(df.plot,aes(x = noise, y = r))+
  geom_line(size=2)+
  geom_vline(xintercept = 0.7,linetype = 2)+
  theme_bw()+
  labs(y = "Pearson's r", x = "degree of noise")+
  coord_cartesian(ylim = c(0.8,0.925))+
  scale_y_continuous(breaks = seq(0.8,0.925,0.025))+
  scale_x_continuous(breaks = seq(0,2,0.2))+
  theme(panel.grid = element_blank(),
        text = element_text(size=20))

# RMSE
ggplot(df.plot,aes(x = noise, y = RMSE))+
  geom_line(size=2)+
  geom_vline(xintercept = 0.7,linetype = 2)+
  theme_bw()+
  labs(y = "root mean squared error", x = "degree of noise")+
  coord_cartesian(ylim = c(12,18))+
  scale_y_continuous(breaks = seq(12,18,2))+
  scale_x_continuous(breaks = seq(0,2,0.2))+
  theme(panel.grid = element_blank(),
        text = element_text(size=20))

Eye-tracking: Static analysis: Results for each individual participant

  • Classification of the different types of looks shown for each participant
df.plot = df.eyes %>% 
  na.omit() %>% 
  filter(frame > 15, frame < t.collision, saccade == "sacc.end") %>% 
  count(condition,participant,counterfactual.look) %>% 
  group_by(condition,participant) %>%
  mutate(freq = n/sum(n)) %>% 
  ungroup() %>% 
  left_join(df.eyes %>% expand(participant,counterfactual.look) %>% na.omit(),.) %>% 
  mutate(condition = ifelse(is.na(condition),3,condition),
         condition = factor(condition,levels=1:3,labels = c("counterfactual","causal","outcome"))) %>% 
  mutate_at(vars(n,freq),funs(ifelse(is.na(.),0,.))) %>% 
  arrange(condition) %>%
  mutate(index = rep(rep(1:10,each=2),3))

ggplot(df.plot,aes(x=index,y=freq,fill = counterfactual.look)) +
  geom_bar(stat = "identity",color = "black")+
  facet_wrap(~condition,scales = "free_x")+
  theme_bw()+
  scale_fill_manual(values= c("green","gray"))+
  labs(y = "percentage", x = "participant index", fill = "type of look")+
  scale_x_continuous(breaks = 1:10)+
  scale_y_continuous(breaks = seq(0,1,.25), labels = paste(seq(0,100,25),'%',sep=""))+
  theme(text = element_text(size=20),
        legend.position = "bottom",
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        strip.background = element_blank(),panel.border = element_rect(colour = "black"),
        strip.text = element_text(size=24),
        legend.text = element_text(size=16),
        axis.title.y = element_text(vjust = 1.5)
  )

Eye-tracking: Static analysis: Counterfacutal looks as a function of the distance parameter

  • Classification of the different types of looks as a function of the distance parameter which determines how close the end point of a saccade needs to be to B’s counterfactual path in order for it to qualify as a ‘counterfactual saccade’
df.fixations = df.eyes %>%
  na.omit() %>%
  expand(condition,counterfactual.look)

threshold.x = 50

for (i in seq(0,300,10)){
# for (i in seq(0,300,50)){
  threshold.distance = i
  df.tmp = df.eyes %>% 
    na.omit() %>% 
    filter(saccade == "sacc.end") %>% 
    group_by(clip) %>% 
    filter(frame < t.collision, frame > 15) %>% 
    mutate(counterfactual.look = ifelse(counterfactual.distance < threshold.distance & 
                                          !is.na(counterfactual.distance) & 
                                          x < (path.start.x-threshold.x), "counterfactual", "other")) %>% 
    mutate(counterfactual.look = factor(counterfactual.look,levels = c("counterfactual","other"))) %>% 
    ungroup %>% 
    group_by(condition,counterfactual.look) %>% 
    summarise(n = n()) %>%
    complete(counterfactual.look,fill = list(n = 0)) %>% 
    mutate(freq = n/sum(n)) %>% 
    select(condition,counterfactual.look,freq) %>%
    left_join(df.fixations,.) %>%  #hack to make sure that all levels are being used
    mutate(freq = ifelse(is.na(freq),0,freq)) %>%
    ungroup
  df.fixations[[paste("threshold_",i,sep="")]] = df.tmp$freq
}

df.fixations = df.fixations %>% 
  gather(threshold,percentage,threshold_0:threshold_300) %>% 
  mutate(threshold = str_replace_all(threshold,"threshold_","") %>% as.numeric())

df.plot = df.fixations %>% 
  mutate(condition = factor(condition, levels = c("counterfactual","causal","outcome"),
                            labels=c(expression("counterfactual\ncondition"),
                                     expression("causal\ncondition"),
                                     expression("outcome\ncondition"))))

ggplot(df.plot,aes(x=threshold,y=percentage,color=counterfactual.look,group=counterfactual.look)) +
  geom_vline(xintercept = 100, linetype = 2)+
  geom_line(stat = "identity",size=2)+
  facet_wrap(~condition)+
  theme_bw()+
  scale_color_manual(values= c("green","gray"))+
  labs(y = "percentage", x = "distance parameter", color = "type of look")+
  scale_x_continuous(breaks = seq(0,300,50))+
  scale_y_continuous(breaks = seq(0,1,.25), labels = paste(seq(0,100,25),'%',sep=""))+
  coord_cartesian(ylim = c(0,1))+
  theme(text = element_text(size=20),
        legend.position = "bottom",
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        strip.background = element_blank(),panel.border = element_rect(colour = "black"),
        strip.text = element_text(size=24),
        legend.text = element_text(size=16),
        plot.margin=unit(c(0.2,0.2,.2,.2),"cm"),
        axis.title.y = element_text(vjust = 1.5)
  )

Eye-tracking: HMM posterior over time

  • Plots the average probability of participants being in different states separated by condition (taking into account all gaze positions, and not just the end points of saccades)
# condition.name = "counterfactual"
condition.name = "causal"
# condition.name = "outcome"
states = c("A look","B look","A predict","B predict","A counterfactual","B counterfactual","other")

smoothing_func = function(x,y){
  fit = loess(y~x,span=0.1)
  out = fit$fitted
  return(out)
}

# cli = 1:18
cli = 4

for (i in cli){

df.plot = df.eyes %>% 
  group_by(clip) %>% 
  filter(frame > 15) %>% 
  # filter(saccade == "sacc.end") %>%
  ungroup %>% 
  filter(condition == condition.name, clip == i) %>% 
  select(participant,instance,frame,contains("post")) %>% 
  arrange(participant,instance,frame) %>% 
  na.omit() %>% 
  group_by(frame) %>% 
  summarise_at(vars(contains("post")),funs(mean)) %>% 
  gather(state,probability,contains("post")) %>%
  mutate(state = factor(state, levels = paste0("post.",c("A.look","B.look","A.predict","B.predict",
                                                         "A.counterfactual","B.counterfactual",
                                                         "other")),labels = states)) %>% 
  mutate(probability = ifelse(is.nan(probability),0,probability)) %>% 
  group_by(state) %>%
  mutate(smoothed = smoothing_func(frame,probability))

t.collision = unique(df.eyes$t.collision[df.eyes$clip == i])
t.outcome.actual = unique(df.eyes$t.outcome.actual[df.eyes$clip == i])
  
p = ggplot(df.plot,aes(x = frame,y = smoothed, color = state)) +
  geom_vline(xintercept = t.collision, linetype=2)+
  geom_vline(xintercept = t.outcome.actual, linetype=2)+
  geom_line(size=2)+
  geom_rect(colour="white",fill="white",aes(xmin=t.collision-25,xmax=t.collision+25 ,ymin=0.95,ymax=1.1))+
  geom_rect(colour="white",fill="white",aes(xmin=t.outcome.actual-25,xmax=t.outcome.actual+25 ,ymin=0.95,ymax=1.1))+
  annotate("text",x = t.collision,y = Inf, label = "collision", hjust = 0.5, vjust = 1.5,size=6)+
  annotate("text",x = t.outcome.actual,y = Inf, label = "outcome", hjust = 0.5, vjust = 1.5,size=6)+
  scale_x_continuous(limits=c(min(df.plot$frame)-1,max(df.plot$frame)+1),breaks=seq(0,max(df.plot$frame)+1,25))+
  scale_y_continuous(limits=c(-0.02,1.1),breaks=seq(0,1,0.25),labels = paste0(seq(0,100,25),"%"))+
  coord_cartesian(ylim = c(0,1))+
  labs(y = "probability of each look", color = "", x = "time frame")+
  theme_bw()+
  theme(
    legend.position = "top",
    panel.grid.major.x = element_blank(), panel.grid.minor = element_blank(),
    strip.background = element_blank(),panel.border = element_rect(colour = "black"),
    axis.title.x = element_text(margin = margin(0.5,0,0,0,"cm")),
    axis.title.y = element_text(margin = margin(0,0.5,0,0,"cm")),
    text=element_text(size=24),
    legend.text = element_text(size=16),
    legend.title = element_text(size=16)
  )+
  guides(color=guide_legend(nrow=2,byrow=TRUE))

print(p)
}

Eye-tracking: Heatmaps

  • Displays a heat map of clip 4 in the causal condition (the code can be adapted to show heat maps for other clips and conditions)
# condition.name = "counterfactual"
condition.name = "causal"
# condition.name = "outcome"

clips = 4
ins = c(1,2)

for (cli in clips){
  
  df.tmp = df.eyes %>% 
    na.omit() %>% 
    # filter(saccade == "sacc.end") %>%
    # filter(frame < t.collision) %>% 
    filter(condition %in% condition.name,
           clip == cli,
           instance %in% ins) 
    
  fit = kde2d(df.tmp$x, df.tmp$y, h = c(100,100),n=500)

  df.plot = expand.grid(x = fit$x, y = fit$y)
  df.plot$z = as.numeric(matrix(unlist(fit$z),ncol=1))
  
  # adjust grid bounds (to avoid polygon errors ...)
  minValue<-sapply(df.plot,min)
  maxValue<-sapply(df.plot,max)
  arbitaryValue=min(df.plot$z)
  
  padding = 30
  bound1=data.frame(x=minValue[["x"]]-padding,y=unique(df.plot$y),z=arbitaryValue)
  bound2=data.frame(x=unique(df.plot$x),y=minValue[["y"]]-padding,z=arbitaryValue)
  bound3=data.frame(x=maxValue[["x"]]+padding,y=unique(df.plot$y),z=arbitaryValue)
  bound4=data.frame(x=unique(df.plot$x),y=maxValue[["y"]]+padding,z=arbitaryValue)
  
  bound=rbind(bound1,bound2,bound3,bound4)
  df.plot = rbind(df.plot,bound)
  
  m = readPNG(paste0("../../figures/diagrams/png/clip_",cli,".png"), FALSE)
  w = matrix(rgb(m[,,1],m[,,2],m[,,3]), nrow=dim(m)[1])
  p = ggplot(df.plot)+
    annotation_custom(xmin=-Inf, ymin=-Inf, xmax=Inf, ymax=Inf, 
                      rasterGrob(w)) +
    stat_contour(aes(x=x,y=y,z=z, color=z,fill=..level..), geom="polygon",alpha=0.05,bins=50)+
    scale_fill_continuous(low = "green", high = "red") +
    # geom_point(data = df.tmp,aes(x=x,y=y),size=3,color = "black", alpha = 0.1) +
    coord_cartesian(xlim=c(0,1024),
                    ylim=c(0,768))+
    labs(x=NULL, y=NULL)+
    theme_void()+
    theme(
      axis.ticks.length = unit(0,"null"),
      legend.position = "none"
    )
  print(p)
}