1 Load packages

library("RSQLite")
library("tidyjson")
library("knitr")
library("png")
library("grid")
library("kableExtra")
library("DT")
library("lme4")
library("broom.mixed")
library("brms")
library("xtable")
library("emmeans")
library("tidyverse")

2 Set options

theme_set(theme_classic() + 
    theme(text = element_text(size = 24)))

opts_chunk$set(comment = "",
               fig.show = "hold")

# suppress grouping warning 
options(dplyr.summarise.inform = F)

# set default color scheme in ggplot 
options(ggplot2.discrete.color = RColorBrewer::brewer.pal(9,"Set1"))
options(ggplot2.discrete.fill = RColorBrewer::brewer.pal(9,"Set1"))

3 Helper functions

# function for printing out html or latex tables 
print_table = function(data, format = "html", digits = 2){
  if(format == "html"){
    data %>% 
      kable(digits = digits) %>% 
      kable_styling()
  }else if(format == "latex"){
    data %>% 
      xtable(digits = digits,
             caption = "Caption",
             label = "tab:table") %>%
      print(include.rownames = F,
            booktabs = T,
            sanitize.colnames.function = identity,
            caption.placement = "top")
  }
}

# function for RMSE 
rmse = function(x, y){
  sqrt(mean((x - y)^2))
}

# extract means and credible intervals from brms model fit
func_brms_params = function(fit, name){
  fit %>% 
    tidy(effects = "fixed",
         fix.intercept = F,
         conf.method = "HPDinterval") %>% 
    select(term, estimate, contains("conf")) %>% 
    mutate(across(-term, ~ round(., 2)),
           estimate = str_c(estimate, " [", conf.low, ", ", conf.high, "]")) %>% 
    select(term, estimate) %>% 
    pivot_wider(names_from = term,
                values_from = estimate) %>% 
    mutate(name = name) %>% 
    relocate(name)
}

4 Data

4.1 Read in data

Restuls from experiment 1 (causal judgment), experiment 2 (counterfactual judgment) and experiment 3 (hypothetical judgment).

# read in data 
con = dbConnect(SQLite(), dbname = "../../data/experiments_anonymized.db");
df.data = dbReadTable(con, "hypothetical_counterfactual")
dbDisconnect(con)

# filter out incompletes 
df.data = df.data %>% 
  filter(status %in% 3:5) %>%
  filter(codeversion %in% c("experiment_1", "experiment_2", "experiment_3"))

# demographics 
df.demographics = df.data$datastring %>% 
  spread_values(age = jstring("questiondata", "age"),
                gender = jstring("questiondata", "gender"),
                race = jstring("questiondata", "race"),
                ethnicity = jstring("questiondata", "ethnicity"),
                feedback = jstring("questiondata", "feedback")) %>% 
  bind_cols(df.data %>% select(codeversion)) %>% 
  as_tibble() %>% 
  rename(participant = document.id,
         condition = codeversion) %>% 
  mutate(time = difftime(df.data$endhit,
                         df.data$beginhit,
                         units = "mins"),
         condition = factor(condition,
                            levels = str_c("experiment_", 1:3),
                            labels = c("causal", "counterfactual", "hypothetical")),
         age = as.numeric(age))

# main data 
df.long = df.data$datastring %>% 
  as.tbl_json() %>% 
  enter_object("data") %>%
  gather_array("order") %>% 
  enter_object("trialdata") %>% 
  gather_object("index") %>% 
  append_values_string("value") %>% 
  as_tibble() %>% 
  rename(participant = document.id) %>% 
  pivot_wider(names_from = index,
              values_from = value) %>% 
  rename(clip = id,
         outcome = gate_pass,
         rating = response) %>% 
  mutate(across(.cols = c(outcome, rating, replay_times, time),
                .fns = ~ as.numeric(.))) %>% 
  mutate(time = time / 1000) %>% 
  select(participant, clip, order, outcome, replay_times, time, rating) %>% 
  filter(!clip %in% c("p1", "p2")) %>% 
  mutate(clip = factor(clip, levels = 1:8),
         outcome = factor(outcome,
                          levels = c(1, 0),
                          labels = c("ball B went in", "ball B missed"))) %>% 
  arrange(participant, clip) %>% 
  left_join(df.demographics %>% 
              select(participant, condition),
            by = "participant") %>% 
  mutate(condition = factor(condition,
                            levels = c("causal", "hypothetical", "counterfactual")))

# filter participants for whom more than 10 clips were incorrectly shown 
# (5 participants total) 
df.long = df.long %>% 
  group_by(participant) %>% 
  filter(max(order) <= 10) %>% 
  ungroup()

4.2 Trial information

df.trialinfo = tibble(clip = 1:8,
                      blocked_initial = c(1, 1, 0, 0, 1, 1, 0, 0),
                      blocked_final = c(0, 1, 1, 0, 0, 1, 1, 0)) %>% 
  mutate(block_moved = (blocked_initial != blocked_final)*1,
         across(everything(), ~ as.factor(.)))

4.3 Model predictions

# noise simulations 
df.noise_simulations = read.csv("data/grid_search.csv") %>% 
  rename(ball_noise = unoise,
         brick_noise = bnoise) %>% 
  select(-X)

# predictions for best fitting parameters
df.model = read.csv("data/model_predictions.csv") %>% 
  rename(clip = X) %>% 
  mutate(clip = clip + 1,
         clip = as.character(clip),
         across(c(hypothetical, counterfactual),
                ~ . * 100))

df.model.ball_only = read.csv("data/model_predictions_ball_only.csv") %>% 
  rename(clip = X) %>% 
  mutate(clip = clip + 1,
         clip = as.character(clip),
         across(c(hypothetical, counterfactual),
                ~ . * 100))

df.model.brick_only = read.csv("data/model_predictions_brick_only.csv") %>% 
  rename(clip = X) %>% 
  mutate(clip = clip + 1,
         clip = as.character(clip),
         across(c(hypothetical, counterfactual),
                ~ . * 100))

5 Experiment 1: Hypothetical and counterfactual judgments

5.1 Stats

5.1.1 Demographics

df.tmp = df.long %>% 
  filter(condition != "causal") %>% 
  distinct(participant, condition) %>% 
  left_join(df.demographics,
            by = c("participant", "condition"))
# age & time
df.tmp %>% 
  summarize(age_mean = mean(age),
            age_sd = sd(age),
            n = n(),
            time_mean = mean(time),
            time_sd = sd(time)) %>% 
  mutate(across(contains("age"), ~ round(.)),
         across(contains("time"), ~ round(., 1))) %>% 
  pivot_longer(cols = everything(),
               values_transform = list(value = as.character))

# A tibble: 5 × 2
  name      value
  <chr>     <chr>
1 age_mean  35   
2 age_sd    10   
3 n         110  
4 time_mean 7.9  
5 time_sd   3.6

# condition
df.tmp %>% 
  count(condition)

# A tibble: 2 × 2
  condition          n
  <fct>          <int>
1 hypothetical      50
2 counterfactual    60

# gender
df.tmp %>% 
  count(gender)

# A tibble: 4 × 2
  gender         n
  <chr>      <int>
1 female        35
2 male          72
3 non-binary     1
4 optout         2

# race
df.tmp %>% 
  count(race)

# A tibble: 5 × 2
  race        n
  <chr>   <int>
1 African    11
2 Asian      10
3 NA          1
4 Native      1
5 White      87

# ethnicity
df.tmp %>% 
  count(ethnicity)

# A tibble: 3 × 2
  ethnicity        n
  <chr>        <int>
1 hispanic        13
2 non-hispanic    96
3 optout           1

5.1.2 Feedback

df.long %>% 
  filter(condition != "causal") %>% 
  distinct(participant, condition) %>% 
  left_join(df.demographics,
            by = c("participant", "condition")) %>% 
  select(participant, condition, feedback) %>% 
  datatable()

5.1.3 Noise parameter search

5.1.3.1 Save means from hypothetical and counterfactual condition

df.means = df.long %>% 
  group_by(condition, clip) %>% 
  summarize(rating = mean(rating)) %>% 
  ungroup()

# write csv with means from hypothetical and counterfactual condition 
df.means %>% 
  filter(condition != "causal") %>% 
  pivot_wider(names_from = condition, 
              values_from = rating) %>% 
  rename(human_cf = counterfactual, 
         human_hp = hypothetical) %>% 
  write.csv(file = str_c("../python/data/", "hpcf_means.csv"))

5.1.3.2 Best parameters

# model with both parameters
df.param.both = df.noise_simulations %>% 
  filter(loss == min(loss)) %>% 
  mutate(name = "both")

# model with only ball noise
df.param.ball = df.noise_simulations %>% 
  filter(brick_noise == 0) %>% 
  filter(loss == min(loss)) %>% 
  mutate(name = "ball only")

# model with only brick noise
df.param.brick = df.noise_simulations %>% 
  filter(ball_noise == 0) %>% 
  filter(loss == min(loss)) %>% 
  mutate(name = "brick only")

bind_rows(df.param.ball,
          df.param.brick,
          df.param.both) %>% 
  relocate(name) %>% 
  print_table()

name	ball_noise	brick_noise	loss
ball only	0.9	0	0.29
brick only	0.0	175	0.17
both	0.6	175	0.13

5.1.4 Regression

df.regression = df.long %>% 
  filter(condition != "causal") %>% 
  left_join(df.trialinfo,
            by = "clip") %>% 
  mutate(outcome = factor(outcome,
                          levels = c("ball B missed", "ball B went in"),
                          labels = c("miss", "hit")),
         across(contains("blocked"), ~ as.factor(.)),
         across(c(contains("blocked"), outcome), ~ C(., sum))) # sum contrasts
 
fit_hypothetical = brm(formula = rating ~ 1 + blocked_initial + 
                           blocked_final +
                           outcome + 
                           (1 | participant),
     data = df.regression %>% 
       filter(condition == "hypothetical"),
     file = "cache/fit_hypothetical",
     seed = 1)

fit_counterfactual = brm(formula = rating ~ 1 + blocked_initial + 
                           blocked_final +
                           outcome + 
                           (1 | participant),
     data = df.regression %>% 
       filter(condition == "counterfactual"),
     file = "cache/fit_counterfactual",
     seed = 1)

5.1.5 Model comparison

5.1.5.1 Overall

df.long %>% 
  filter(condition != "causal") %>% 
  group_by(condition,
           outcome,
           clip) %>% 
  summarize(rating = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.model %>% 
              pivot_longer(cols = -clip, 
                           names_to = "condition",
                           values_to = "prediction"),
            by = c("clip", "condition")) %>% 
  summarize(r = cor(rating, prediction),
            rs = cor(rating, prediction, method = "spearman"),
            rmse = rmse(rating, prediction)) %>% 
  mutate(across(.fns = ~ round(., 2))) %>% 
  print_table()

r	rs	rmse
0.97	0.97	8.97

5.1.5.2 Separated by condition

df.long %>% 
  filter(condition != "causal") %>% 
  group_by(condition,
           outcome,
           clip) %>% 
  summarize(rating = mean(rating)) %>% 
  ungroup() %>% 
  pivot_wider(names_from = condition,
              values_from = rating) %>% 
  rename_with(.cols = c(hypothetical, counterfactual),
              .fn = ~ str_c(., "_rating")) %>% 
  left_join(df.model %>% 
              rename_with(.cols = c(hypothetical, counterfactual),
              .fn = ~ str_c(., "_model")),
            by = "clip") %>% 
  summarize(r.hypothetical = cor(hypothetical_rating, hypothetical_model),
            r.counterfactual = cor(counterfactual_rating, counterfactual_model),
            rs.hypothetical = cor(hypothetical_rating, hypothetical_model,
                                  method = "spearman"),
            rs.counterfactual = cor(counterfactual_rating, counterfactual_model,
                                    method = "spearman"),
            rmse.hypothetical = rmse(hypothetical_rating, hypothetical_model),
            rmse.counterfactual = rmse(counterfactual_rating, counterfactual_model)) %>% 
  mutate(across(.fns = ~ round(., 2))) %>% 
  pivot_longer(cols = everything(), 
               names_to = c("measure", "comparison"),
               names_sep = "\\.") %>% 
  pivot_wider(names_from = measure,
              values_from = value) %>% 
  print_table()

comparison	r	rs	rmse
hypothetical	0.89	0.76	5.19
counterfactual	0.99	0.98	11.57

5.2 Plots

5.2.1 Hypothetical and counterfactual judgments vs. model predictions

set.seed(1)

df.plot = df.long %>% 
  filter(condition %in% c("counterfactual", "hypothetical")) %>% 
  mutate(outcome = factor(outcome, levels = c("ball B went in", "ball B missed")))

df.plot2 = df.model %>%
# df.plot2 = df.model.ball_only %>% 
# df.plot2 = df.model.brick_only %>% 
  pivot_longer(cols = -clip,
               names_to = "condition",
               values_to = "rating") %>% 
  mutate(outcome = ifelse(clip <= 4, "ball B went in", "ball B missed"),
         outcome = factor(outcome, levels = c("ball B went in","ball B missed")),
         condition = factor(condition, levels = c("hypothetical", "counterfactual"))) %>% 
  mutate(clip = factor(clip))

p = ggplot(data = df.plot,
           mapping = aes(x = clip, 
                         y = rating,
                         fill = condition)) +
  stat_summary(fun = "mean",
               geom = "bar",
               position = position_dodge(width = 0.9),
               color = "black",
               alpha = 0.5) + 
  stat_summary(fun.data = "mean_cl_boot",
               geom = "linerange",
               position = position_dodge(width = 0.9),
               size = 1.5) + 
  geom_point(data = df.plot2,
             shape = 21,
             position = position_dodge(0.9),
             size = 4,
             show.legend = F) + 
  facet_grid(cols = vars(outcome),
             scales = "free") +
  ylab("mean rating") + 
  scale_y_continuous(breaks = seq(0, 100, 25),
                     labels = seq(0, 100, 25),
                     expand = c(0, 0)) + 
  scale_fill_brewer(palette = "Set1",
                    direction = -1) + 
  coord_cartesian(clip = "off",
                  ylim = c(0, 100)) + 
  theme(legend.position = "bottom",
        legend.margin = margin(t = 2.8, unit = "cm"),
        legend.background = element_blank(),
        axis.title.x = element_blank()) 
    
  
# add images of bricks as x-axis labels  
g = list()
image_names = c("down_up", "down", "up_down", "up")

for(i in 1:length(image_names)){
  img = readPNG(str_c("../../figures/diagrams/blocks/", image_names[i], ".png"))
  g[[i]] = rasterGrob(img, interpolate = TRUE)
  
  p = p +
    annotation_custom(grob = g[[i]],
                      xmin = i - 1,
                      xmax = i + 1,
                      ymin = -40,
                      ymax = -10)
}
print(p)

ggsave(filename = str_c("../../figures/plots/", 
                        "hypothetical_counterfactual_judgments.pdf"),
       plot = p,
       width = 8,
       height = 6)

5.2.2 Noise parameters

ggplot(data = df.noise_simulations,
       mappin = aes(x = ball_noise, y = brick_noise)) +
  geom_tile(aes(fill = loss),
            color = "black") +
  scale_fill_gradient(low = "white", high = "black") + 
  coord_cartesian(expand = F) +
  scale_x_continuous(breaks = seq(0, 1.2, 0.2)) + 
  labs(x = expression("ball motion noise"~sigma["ball"]),
       y = expression("block movement noise"~~sigma["block"]))
  
ggsave(filename = "../../figures/plots/parameter_search.pdf",
       width = 8,
       height = 6)

5.3 Tables

func_brms_params(fit_hypothetical, "hypothetical") %>% 
  bind_rows(func_brms_params(fit_counterfactual, "counterfactual")) %>% 
  rename(intercept = Intercept,
         block_initial = blocked_initial1,
         block_final = blocked_final1,
         outcome = outcome1) %>% 
  print_table()

name	intercept	block_initial	block_final	outcome
hypothetical	49.26 [44.95, 53.61]	9.67 [6.64, 12.66]	1.84 [-1.08, 4.94]	-2.08 [-5.3, 1.03]
counterfactual	51.2 [48.02, 54.46]	3.05 [0.48, 5.75]	28.51 [25.91, 31.11]	0.59 [-2.03, 3.29]

6 Experiment 2: Causal judgments

6.1 Stats

6.1.1 Demographics

df.tmp = df.long %>% 
  filter(condition == "causal") %>% 
  distinct(participant, condition) %>% 
  left_join(df.demographics,
            by = c("participant", "condition"))

# condition
df.tmp %>% 
  count(condition)

# A tibble: 1 × 2
  condition     n
  <fct>     <int>
1 causal       67

# age & time
df.tmp %>% 
  summarize(age_mean = mean(age),
            age_sd = sd(age),
            n = n(),
            time_mean = mean(time),
            time_sd = sd(time)) %>% 
  mutate(across(contains("age"), ~ round(.)),
         across(contains("time"), ~ round(., 1))) %>% 
  pivot_longer(cols = everything(),
               values_transform = list(value = as.character))

# A tibble: 5 × 2
  name      value
  <chr>     <chr>
1 age_mean  34   
2 age_sd    10   
3 n         67   
4 time_mean 7.3  
5 time_sd   2.1

# gender
df.tmp %>% 
  count(gender)

# A tibble: 3 × 2
  gender         n
  <chr>      <int>
1 female        20
2 male          46
3 non-binary     1

# ethnicity
df.tmp %>% 
  count(ethnicity)

# A tibble: 3 × 2
  ethnicity        n
  <chr>        <int>
1 hispanic         9
2 non-hispanic    57
3 optout           1

6.1.2 Feedback

df.long %>% 
  filter(condition == "causal") %>% 
  distinct(participant, condition) %>% 
  left_join(df.demographics,
            by = c("participant", "condition")) %>% 
  select(participant, condition, feedback) %>% 
  datatable()

6.1.3 Regression

df.regression = df.long %>% 
  filter(condition == "causal") %>% 
  left_join(df.trialinfo,
            by = "clip") %>% 
  mutate(outcome = factor(outcome,
                          levels = c("ball B missed", "ball B went in"),
                          labels = c("miss", "hit"))) %>% 
  mutate(across(contains("blocked"), ~ as.factor(.)),
         across(c(contains("blocked"), outcome), ~ C(., sum))) # setting sum contrasts
 
fit_causal = brm(formula = rating ~ 1 + (blocked_initial + 
                           blocked_final) *
                           outcome + 
                           (1 | participant),
     data = df.regression,
     file = "cache/fit_causal",
     seed = 1)

fit_causal

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ 1 + (blocked_initial + blocked_final) * outcome + (1 | participant) 
   Data: df.regression (Number of observations: 536) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~participant (Number of levels: 67) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     9.19      1.76     5.78    12.65 1.00     1417     2231

Population-Level Effects: 
                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                    61.46      1.60    58.19    64.57 1.00     3016
blocked_initial1              0.11      1.15    -2.15     2.41 1.00     5853
blocked_final1                0.71      1.12    -1.48     2.87 1.00     5673
outcome1                     -1.02      1.12    -3.15     1.17 1.00     5054
blocked_initial1:outcome1     1.91      1.16    -0.41     4.18 1.01     5639
blocked_final1:outcome1      23.97      1.13    21.75    26.17 1.00     5293
                          Tail_ESS
Intercept                     2285
blocked_initial1              2892
blocked_final1                3019
outcome1                      2939
blocked_initial1:outcome1     3118
blocked_final1:outcome1       2657

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    26.23      0.88    24.61    28.04 1.00     3474     2961

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

6.1.4 Model comparison

df.means = df.long %>% 
  group_by(condition,
           outcome,
           clip) %>% 
  summarize(rating = mean(rating)) %>% 
  ungroup() %>% 
  mutate(rating = ifelse(condition != "causal" & 
                             outcome == "ball B went in",
                         100 - rating,
                         rating)) %>% 
  pivot_wider(names_from = condition,
              values_from = rating)

df.means %>% 
  summarize(r.causal_hypothetical = cor(causal, hypothetical),
            r.causal_counterfactual = cor(causal, counterfactual),
            rs.causal_hypothetical = cor(causal, hypothetical, method = "spearman"),
            rs.causal_counterfactual = cor(causal, counterfactual, method = "spearman"),
            rmse.causal_hypothetical = rmse(causal, hypothetical),
            rmse.causal_counterfactual = rmse(causal, counterfactual)) %>% 
  mutate(across(.fns = ~ round(., 2))) %>% 
  pivot_longer(cols = everything(), 
               names_to = c("measure", "comparison"),
               names_sep = "\\.") %>% 
  pivot_wider(names_from = measure,
              values_from = value) %>% 
  print_table()

comparison	r	rs	rmse
causal_hypothetical	0.24	0.40	27.31
causal_counterfactual	0.99	0.79	12.50

6.2 Plots

6.2.1 Causal judgments

set.seed(1)

df.plot = df.long %>% 
  filter(condition == "causal") %>% 
  mutate(outcome = factor(outcome, levels = c("ball B went in", "ball B missed")))

p = ggplot(data = df.plot,
           mapping = aes(x = clip,
                         y = rating)) +
  stat_summary(fun = "mean",
               geom = "bar",
               color = "black",
               fill = "gray90") +
  stat_summary(fun.data = "mean_cl_boot",
               geom = "linerange",
               color = "black",
               size = 1.5) +
  stat_summary(data = df.long %>% 
                 filter(condition == "hypothetical") %>% 
                 mutate(rating = ifelse(outcome == "ball B went in", 100 - rating, rating)),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               aes(fill = "hypothetical"),
               position = position_nudge(x = -0.2),
               shape = 21,
               size = 1,
               fatten = 5,
               stroke = 1) + 
  stat_summary(data = df.long %>%
                 filter(condition == "counterfactual") %>%
                 mutate(rating = ifelse(outcome == "ball B went in", 100 - rating, rating)),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               aes(fill = "counterfactual"),
               position = position_nudge(x = 0.2),
               shape = 21,
               size = 1,
               fatten = 5,
               stroke = 1) +
  facet_grid(cols = vars(outcome),
             scales = "free") +
  labs(y = "mean causal rating",
       fill = "model prediction") + 
  scale_y_continuous(breaks = seq(0, 100, 25),
                     labels = seq(0, 100, 25),
                     expand = c(0, 0)) + 
  scale_fill_manual(values = c("hypothetical" = "#377eb8",
                               "counterfactual" = "#e41a1c")) +
  coord_cartesian(clip = "off",
                  ylim = c(0, 100)) + 
  theme(legend.position = "bottom",
        legend.margin = margin(t = 2.8, unit = "cm"),
        legend.background = element_blank(),
        axis.title.x = element_blank()) 

# add images of bricks as x-axis labels  
g = list()
image_names = c("down_up", "down", "up_down", "up")

for(i in 1:length(image_names)){
  img = readPNG(str_c("../../figures/diagrams/blocks/", image_names[i], ".png"))
  g[[i]] = rasterGrob(img, interpolate = TRUE)
  
  p = p +
    annotation_custom(grob = g[[i]],
                      xmin = i - 1,
                      xmax = i + 1,
                      ymin = -40,
                      ymax = -10)
}

print(p)

ggsave(filename = "../../figures/plots/causal_ratings_all.pdf",
       plot = p,
       width = 8,
       height = 6)

6.3 Tables

6.3.1 brms

func_brms_params(fit_causal, "causal") %>% 
  rename(intercept = Intercept,
         block_initial = blocked_initial1,
         block_final = blocked_final1,
         outcome = outcome1,
         `block_initial:outcome` = `blocked_initial1:outcome1`,
         `block_final:outcome` = `blocked_final1:outcome1`) %>% 
  print_table()

name	intercept	block_initial	block_final	outcome	block_initial:outcome	block_final:outcome
causal	61.46 [58.1, 64.42]	0.11 [-2.31, 2.22]	0.71 [-1.58, 2.78]	-1.02 [-3.22, 1.07]	1.91 [-0.3, 4.25]	23.97 [21.72, 26.13]

7 Session info

R version 4.1.2 (2021-11-01)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] forcats_0.5.1     stringr_1.4.0     dplyr_1.0.7       purrr_0.3.4      
 [5] readr_2.1.0       tidyr_1.1.4       tibble_3.1.6      ggplot2_3.3.5    
 [9] tidyverse_1.3.1   emmeans_1.7.1-1   xtable_1.8-4      brms_2.16.3      
[13] Rcpp_1.0.7        broom.mixed_0.2.7 lme4_1.1-27.1     Matrix_1.3-4     
[17] DT_0.20           kableExtra_1.3.4  png_0.1-7         knitr_1.36       
[21] tidyjson_0.3.1    RSQLite_2.2.8    

loaded via a namespace (and not attached):
  [1] readxl_1.3.1         backports_1.3.0      Hmisc_4.6-0         
  [4] systemfonts_1.0.3    plyr_1.8.6           igraph_1.2.9        
  [7] splines_4.1.2        crosstalk_1.2.0      TH.data_1.1-0       
 [10] rstantools_2.1.1     inline_0.3.19        digest_0.6.28       
 [13] htmltools_0.5.2      rsconnect_0.8.25     fansi_0.5.0         
 [16] magrittr_2.0.1       checkmate_2.0.0      memoise_2.0.0       
 [19] cluster_2.1.2        tzdb_0.2.0           modelr_0.1.8        
 [22] RcppParallel_5.1.4   matrixStats_0.61.0   sandwich_3.0-1      
 [25] xts_0.12.1           svglite_2.0.0        prettyunits_1.1.1   
 [28] jpeg_0.1-9           colorspace_2.0-2     blob_1.2.2          
 [31] rvest_1.0.2          haven_2.4.3          xfun_0.28           
 [34] callr_3.7.0          crayon_1.4.2         jsonlite_1.7.2      
 [37] survival_3.2-13      zoo_1.8-9            glue_1.5.0          
 [40] gtable_0.3.0         webshot_0.5.2        V8_3.6.0            
 [43] distributional_0.2.2 pkgbuild_1.2.1       rstan_2.21.2        
 [46] abind_1.4-5          scales_1.1.1         mvtnorm_1.1-3       
 [49] DBI_1.1.1            miniUI_0.1.1.1       htmlTable_2.3.0     
 [52] viridisLite_0.4.0    diffobj_0.3.5        foreign_0.8-81      
 [55] bit_4.0.4            Formula_1.2-4        stats4_4.1.2        
 [58] StanHeaders_2.21.0-7 htmlwidgets_1.5.4    httr_1.4.2          
 [61] threejs_0.3.3        RColorBrewer_1.1-2   posterior_1.1.0     
 [64] ellipsis_0.3.2       pkgconfig_2.0.3      loo_2.4.1           
 [67] farver_2.1.0         nnet_7.3-16          dbplyr_2.1.1        
 [70] sass_0.4.0           utf8_1.2.2           labeling_0.4.2      
 [73] tidyselect_1.1.1     rlang_0.4.12         reshape2_1.4.4      
 [76] later_1.3.0          cellranger_1.1.0     munsell_0.5.0       
 [79] tools_4.1.2          cachem_1.0.6         cli_3.1.0           
 [82] generics_0.1.1       broom_0.7.10         ggridges_0.5.3      
 [85] evaluate_0.14        fastmap_1.1.0        yaml_2.2.1          
 [88] fs_1.5.0             processx_3.5.2       bit64_4.0.5         
 [91] nlme_3.1-153         mime_0.12            rstanarm_2.21.1     
 [94] xml2_1.3.2           compiler_4.1.2       bayesplot_1.8.1     
 [97] shinythemes_1.2.0    rstudioapi_0.13      curl_4.3.2          
[100] reprex_2.0.1         bslib_0.3.1          stringi_1.7.5       
[103] highr_0.9            ps_1.6.0             Brobdingnag_1.2-6   
[106] lattice_0.20-45      nloptr_1.2.2.3       markdown_1.1        
[109] shinyjs_2.0.0        tensorA_0.36.2       vctrs_0.3.8         
[112] pillar_1.6.4         lifecycle_1.0.1      jquerylib_0.1.4     
[115] bridgesampling_1.1-2 estimability_1.3     data.table_1.14.2   
[118] httpuv_1.6.3         latticeExtra_0.6-29  R6_2.5.1            
[121] bookdown_0.24        promises_1.2.0.1     gridExtra_2.3       
[124] codetools_0.2-18     boot_1.3-28          colourpicker_1.1.1  
[127] MASS_7.3-54          gtools_3.9.2         assertthat_0.2.1    
[130] withr_2.4.2          shinystan_2.5.0      multcomp_1.4-17     
[133] hms_1.1.1            parallel_4.1.2       rpart_4.1-15        
[136] coda_0.19-4          minqa_1.2.4          rmarkdown_2.11      
[139] lubridate_1.8.0      shiny_1.7.1          base64enc_0.1-3     
[142] dygraphs_1.1.1.6

Hypotheticals vs. counterfactuals

Tobias Gerstenberg

January 18, 2022

1 Load packages

2 Set options

3 Helper functions

4 Data

4.1 Read in data

4.2 Trial information

4.3 Model predictions

5 Experiment 1: Hypothetical and counterfactual judgments

5.1 Stats

5.1.1 Demographics

5.1.2 Feedback

5.1.3 Noise parameter search

5.1.3.1 Save means from hypothetical and counterfactual condition

5.1.3.2 Best parameters

5.1.4 Regression

5.1.5 Model comparison

5.1.5.1 Overall

5.1.5.2 Separated by condition

5.2 Plots

5.2.1 Hypothetical and counterfactual judgments vs. model predictions

5.2.2 Noise parameters

5.3 Tables

6 Experiment 2: Causal judgments

6.1 Stats

6.1.1 Demographics

6.1.2 Feedback

6.1.3 Regression

6.1.4 Model comparison

6.2 Plots

6.2.1 Causal judgments

6.3 Tables

6.3.1 brms

7 Session info