1 Setup

2 Helper functions

2.1 Demographics

This classification is based on the US Census Bureau, which counts people of MENA (Middle Eastern and North African) descent as White. Hispanic/Latinx is not a race according to the Census (it’s an ethnicity).

fem = c('female', 'feminine', 'f', 'women')
masc = c('male', 'man', 'm', 'trans man', 'masculine', 'masculino')
nonb = c('non-binary')

white = c('caucasian', 'caucasian/white', 'caucasic', 'caucasoid', 'european', 'greek - white',
          'italian', 'polish', 'middle eastern', 'egyptian', 'mediterranean white', 'mena',
          'south european caucasian', 'white', 'white british', 'white/caucasian', 'whitr')
black = c('black', 'black/african american', 'north african')
asian = c('asian', 'indian')
ind_nat = c('american indian/alaska native')
multi = c('mixed', 'asian and white', 'asian/white', 'black/aa + asian', 'multi-racial',
          'white + asian')
other = c('brown', 'latinx/hispanic/caribbean', 'prefer not to say')

print_demographics = function(data) {
  data = data %>%
    mutate(gender = tolower(gender),
           race = tolower(race),
           gender = case_when(gender %in% fem ~ 'Female',
                              gender %in% masc ~ 'Male',
                              gender %in% nonb ~ 'Non-binary',
                              TRUE ~ NA_character_),
           race = case_when(race %in% white ~ 'White',
                            race %in% black ~ 'Black',
                            race %in% asian ~ 'Asian',
                            race %in% ind_nat ~ 'American Indian/Alaska Native',
                            race %in% multi ~ 'Multiracial',
                            race %in% other ~ 'Undisclosed/other',
                            TRUE ~ NA_character_))
  # gender
  data %>%
    group_by(gender) %>%
    summarise(n = n()) %>%
    print()
  
  # age
  print('age:')
  mean(data$age, na.rm = T) %>% print()
  sd(data$age, na.rm = T) %>% print()
  
  # race
  data %>%
    group_by(race) %>%
    summarise(n = n()) %>%
    print()
  
  # ethnicity
  data %>%
    group_by(ethnicity) %>%
    summarise(n = n()) %>%
    print()

  # time taken
  print('time taken (min):')
  print(mean(data$time, na.rm = T)/60/1000)
  print(sd(data$time, na.rm = T)/60/1000)
}

2.2 Parameter search

2.2.1 Uniform probability

update_search_one = function(exp, condition, p) {
  type = ifelse(condition == 'agent', 'craft', 'color')
  
  data = get(paste0('data.exp', exp, '.', condition)) %>%
    mutate(prob = 1 - (1-p)^num)
  info = get(paste0('info.exp', exp, '.', condition))
  model = lmer(resp ~ 1 + prob + (1 | id),
               data = data)
  data_means = data %>%
    mutate(pred = model %>%
            predict(newdata = data,
                    re.form = NA)) %>%
    group_by_at(c('trial', type, 'pred')) %>%
    summarise(resp = mean(resp),
              pred = mean(pred)) %>%
    merge(info %>%
            select(trial, type),
          by = c('trial', type))
  rmse = rmse(data_means$resp, data_means$pred)
  data_rmse = get(paste0('rmse.exp', exp, '.', condition, '_uniform')) %>%
    rbind(c(p, rmse))
  
  return(data_rmse)
}

plot_search_one = function(exp, condition) {
  data_rmse =  get(paste0('rmse.exp', exp, '.', condition, '_uniform'))
  
  g = ggplot(data_rmse,
             aes(x = p,
                 y = rmse)) +
    geom_point() +
    # best parametrization
    geom_point(data = data_rmse %>%
                 filter(rmse == min(rmse)),
               color = 'red') +
    ggtitle(paste(condition, 'condition')) +
    theme(text = element_text(size = 12),
          axis.ticks.x = element_blank()) +
    scale_x_continuous(name = 'p',
                       limits = c(0, 1),
                       breaks = seq(0, 1, by = 0.2)) +
    scale_y_continuous(name = 'RMSE',
                       limits = c(0, 5.5),
                       breaks = seq(0, 5))
  return(g)
}

2.2.2 Varying probabilities

update_search_two = function(exp, condition, p_low, p_high) {
  type = ifelse(condition == 'agent', 'craft', 'color')
  
  data = get(paste0('data.exp', exp, '.', condition)) %>%
    mutate(prob = 1 - (1-p_low)^low * (1-p_high)^high)
  info = get(paste0('info.exp', exp, '.', condition))
  model = lmer(resp ~ 1 + prob + (1 | id),
               data = data)
  data_means = data %>%
    mutate(pred = model %>%
            predict(newdata = data,
                    re.form = NA)) %>%
    group_by_at(c('trial', type, 'pred')) %>%
    summarise(resp = mean(resp),
              pred = mean(pred)) %>%
    merge(info %>%
            select(trial, type),
          by = c('trial', type))
  rmse = rmse(data_means$resp, data_means$pred)
  data_rmse = get(paste0('rmse.exp', exp, '.', condition, '_full')) %>%
    rbind(c(p_low, p_high, rmse))
  
  return(data_rmse)
}

plot_search_two = function(exp, condition) {
  data_rmse =  get(paste0('rmse.exp', exp, '.', condition, '_full')) %>%
    mutate(p_low = as.factor(p_low),
           p_high = as.factor(p_high))
  
  g = ggplot(data_rmse,
             aes(x = p_low,
                 y = p_high,
                 fill = rmse)) +
    geom_tile(color = 'black') +
    coord_fixed() + # square tiles
    # best parametrization
    geom_point(data = data_rmse %>%
                 filter(rmse == min(rmse)),
               color = 'red',
               size = 2) +
    scale_fill_gradient(low = "gray20", high = "white") +
    guides(fill = guide_colourbar(barwidth = 0.5,
                                  barheight = 5,
                                  title = 'RMSE')) +
    ggtitle(paste(condition, 'condition')) +
    theme(text = element_text(size = 12),
          legend.position = c(0.9, 0.3),
          axis.ticks.x = element_blank()) +
    scale_x_discrete(name = 'low availability p',
                     breaks = seq(0, 1, by = 0.2)) +
    scale_y_discrete(name = 'high availability p',
                     breaks = seq(0, 1, by = 0.2))
  return(g)
}

2.3 Make scatter plot

make_scatter = function(data_means, model_name, xlabel, breaks) {
  col_name = paste('pred_', model_name, sep = '')
  
  if (model_name == 'baseline') {
    r = 'N/A'
  } else {
    r = cor(data_means$resp_mean, data_means[[col_name]]) %>%
      round(2)
  }
  
  rmse = rmse(data_means$resp_mean, data_means[[col_name]])
  
  g = ggplot(data = data_means,
             aes(x = get(col_name),
                 y = resp_mean)) +
    geom_abline(intercept = 0, slope = 1, linetype = 2) +
    geom_linerange(size = 0.5, 
                   mapping = aes(ymin = resp_low,
                                 ymax = resp_high),
                   color = "lightgray") +
    geom_point(size = 2, shape = 21, stroke = 0.5,
               fill = 'gray', color = 'black') +
    theme(text = element_text(size = 12),
          legend.position = "none") +
    annotate('text',
             label = sprintf('r = %s\nRMSE = %.2f', r, rmse),
             hjust = 0,   # left align
             x = breaks[1], y = breaks[length(breaks)] - 2) +
    scale_x_continuous(name = xlabel,
                       breaks = breaks,
                       limits = c(breaks[1], breaks[length(breaks)])) +
    scale_y_continuous(name = 'responsibility judgment',
                       breaks = breaks,
                       limits = c(breaks[1], breaks[length(breaks)])) +
    coord_fixed(ratio = 1)

  return(g)
}

agent_icon = readPNG(paste0(figures_dir, 'icons/agent.png')) %>%
  rasterGrob(interpolate = T)
object_icon = readPNG(paste0(figures_dir, 'icons/object.png')) %>%
  rasterGrob(interpolate = T)

3 Experiment 1

3.1 Data

3.1.1 Agent condition

3.1.1.1 Trial info

info.exp1.agent = read_csv(paste0(data_dir, 'experiment1_agent/trial_info.csv'),
                           show_col_types = F) %>%
  pivot_longer(cols = !trial,
               names_to = c('.value', 'craft'),
               names_pattern = '(.+)_(.)')

3.1.1.2 Response data

data.exp1.agent = read_csv(paste0(data_dir, 'experiment1_agent/trials.csv'),
                           show_col_types = F) %>%
  select(-error, -proliferate.condition) %>%
  rename(id = workerid,
         b = resp_b, c = resp_c, t = resp_t) %>%
  # filter by exclusion trial X1
  pivot_wider(names_from = trial, values_from = c(c, b, t)) %>%
  rowwise(id) %>%
  filter(max(c_X1, b_X1, t_X1) - min(c_X1, b_X1, t_X1) < 30) %>%
  # pivot back to long form
  pivot_longer(cols = !id,
               names_to = 'trial',
               values_to = 'resp') %>%
  separate(col = trial,
           into = c('craft', 'trial'),
           sep = '_') %>%
  # ignore exclusion trial and practice trials
  filter(!(trial %in% c("X1", 1, 2))) %>%
  # merge with trial info
  merge(info.exp1.agent,
        by = c('trial', 'craft')) %>%
  # factor variables for plotting
  mutate(id = factor(id),
         plot_id = as.factor(as.integer(id)),
         craft = factor(craft))

Number of participants:

length(unique(data.exp1.agent$id))

## [1] 50

3.1.1.3 Participant data

Participants were asked: “What factors influenced how you decided to respond? Do you have any questions or comments regarding the experiment?”

participants.exp1.agent = read_csv(paste0(data_dir, 'experiment1_agent/participants.csv'),
                                   show_col_types = F) %>%
  select(-error, -proliferate.condition) %>%
  rename(id = workerid) %>%
  filter(id %in% data.exp1.agent$id) %>%
  mutate(id = factor(id)) %>%
  left_join(data.exp1.agent %>%
              distinct(id, plot_id),
            by = 'id')

participants.exp1.agent %>%
  select(plot_id, feedback) %>% 
  datatable()

3.1.2 Object condition

3.1.2.1 Trial info

info.exp1.object = read_csv(paste0(data_dir, 'experiment1_object/trial_info.csv'),
                            show_col_types = F) %>%
  pivot_longer(cols = !trial,
               names_to = c('.value', 'color'),
               names_pattern = '(.+)_(.)')

3.1.2.2 Response data

data.exp1.object = read_csv(paste0(data_dir, 'experiment1_object/trials.csv'),
                            show_col_types = F) %>%
  select(-error, -proliferate.condition) %>%
  rename(id = workerid,
         y = resp_y, g = resp_g, b = resp_b) %>%
  # filter by exclusion trial X1
  pivot_wider(names_from = trial, values_from = c(y, g, b)) %>%
  rowwise(id) %>%
  filter(max(y_X1, g_X1, b_X1) - min(y_X1, g_X1, b_X1) < 30) %>%
  # pivot back to long form
  pivot_longer(cols = !id,
               names_to = 'trial',
               values_to = 'resp') %>%
  separate(col = trial,
           into = c('color', 'trial'),
           sep = '_') %>%
  # ignore exclusion trial and practice trials
  filter(!(trial %in% c("X1", 1, 2))) %>%
  # merge with trial info
  merge(info.exp1.object,
        by = c('trial', 'color')) %>%
  # factor variables for plotting
  mutate(id = factor(id),
         plot_id = as.factor(as.integer(id)),
         color = factor(color))

Number of participants:

length(unique(data.exp1.object$id))

## [1] 50

3.1.2.3 Participant data

Participants were asked: “What factors influenced how you decided to respond? Do you have any questions or comments regarding the experiment?”

participants.exp1.object = read_csv(paste0(data_dir, 'experiment1_object/participants.csv'),
                                    show_col_types = F) %>%
  select(-error, -proliferate.condition) %>%
  rename(id = workerid) %>%
  filter(id %in% data.exp1.object$id) %>%
  mutate(id = factor(id)) %>%
  left_join(data.exp1.object %>%
              distinct(id, plot_id),
            by = 'id')

participants.exp1.object %>%
  select(plot_id, feedback) %>% 
  datatable()

3.1.3 Demographics

participants.exp1.agent %>%
  rbind(participants.exp1.object) %>%
  print_demographics()

## # A tibble: 4 × 2
##   gender         n
##   <chr>      <int>
## 1 Female        34
## 2 Male          63
## 3 Non-binary     1
## 4 <NA>           2
## [1] "age:"
## [1] 24.84
## [1] 5.814854
## # A tibble: 5 × 2
##   race            n
##   <chr>       <int>
## 1 Asian           7
## 2 Black           7
## 3 Multiracial     3
## 4 White          64
## 5 <NA>           19
## # A tibble: 3 × 2
##   ethnicity        n
##   <chr>        <int>
## 1 Hispanic        29
## 2 Non-Hispanic    70
## 3 <NA>             1
## [1] "time taken (min):"
## [1] 9.846081
## [1] 5.648796

3.2 Models

3.2.1 Agent condition

3.2.1.1 Parameter search

rmse.exp1.agent_uniform = data.frame()
for (p in seq(0, 1, by = 0.05)) {
  rmse.exp1.agent_uniform = update_search_one(1, 'agent', p)
}
colnames(rmse.exp1.agent_uniform) = c('p', 'rmse')

# get optimal values
p.exp1.agent = rmse.exp1.agent_uniform[which.min(rmse.exp1.agent_uniform$rmse),]$p

# create plot
plot_rmse.exp1.agent = plot_search_one(1, 'agent')

# update response data
data.exp1.agent = data.exp1.agent %>%
  mutate(prob = 1 - (1 - p.exp1.agent)^num)

3.2.1.2 Overall models

Fitting baseline model:

fit.exp1.agent_baseline = brm(
  formula = resp ~ 1 + (1 | id),
  data = data.exp1.agent,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp1.agent_baseline")

Fitting counterfactual replacement model (CRM):

fit.exp1.agent_crm = brm(
  formula = resp ~ 1 + prob +  (1 + prob | id),
  data = data.exp1.agent,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp1.agent_crm")

Comparing models with leave-one-out cross-validation:

fit.exp1.agent_baseline = add_criterion(
  fit.exp1.agent_baseline,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp1.agent_baseline")

fit.exp1.agent_crm = add_criterion(
  fit.exp1.agent_crm,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp1.agent_crm")

loo_compare(fit.exp1.agent_baseline,
            fit.exp1.agent_crm)

##                         elpd_diff se_diff
## fit.exp1.agent_crm          0.0       0.0
## fit.exp1.agent_baseline -1591.1      74.2

3.2.1.3 Individual models

prior = c(prior(normal(0, 20), class = 'b', ub = 0))

# initialize the models 
fit.exp1.agent_baseline_individual = brm(
  formula = resp ~ 1,
  data = data.exp1.agent %>% 
    filter(plot_id == 1),
  seed = 1,
  file = "cache/fit.exp1.agent_baseline_individual")

fit.exp1.agent_crm_individual = brm(
  formula = resp ~ 1 + prob,
  data = data.exp1.agent %>% 
    filter(plot_id == 1),
  prior = prior,
  seed = 1,
  file = "cache/fit.exp1.agent_crm_individual")

# update model fits for each participant 
data.exp1.agent_model_fits = data.exp1.agent %>% 
  group_by(plot_id) %>% 
  nest() %>% 
  ungroup() %>% 
  mutate(fit_baseline = map(.x = data,
                            .f = ~ update(fit.exp1.agent_baseline_individual,
                                          newdata = .x,
                                          seed = 1)),
         fit_crm = map(.x = data,
                       .f = ~ update(fit.exp1.agent_crm_individual,
                                     newdata = .x,
                                     seed = 1))) %>%
  mutate(fit_baseline = map(.x = fit_baseline,
                            .f = ~ add_criterion(.x, criterion = "loo")),
         fit_crm = map(.x = fit_crm,
                       .f = ~ add_criterion(.x, criterion = "loo")),
         model_comparison = pmap(.l = list(baseline = fit_baseline,
                                           crm = fit_crm),
                                 .f = ~ loo_compare(..1, ..2)),
         best_model = map_chr(.x = model_comparison,
                              .f = ~ rownames(.) %>% 
                                .[1]),
         best_model = factor(best_model,
                             levels = c("..1", "..2"),
                             labels = c("baseline",
                                        "crm")))

save(list = c("data.exp1.agent_model_fits"),
     file = paste0(data_dir, 'experiment1_agent/model-fits.RData'))

load(file = paste0(data_dir, 'experiment1_agent/model-fits.RData'))

data.exp1.agent_model_fits %>% 
  count(best_model) %>%
  arrange(desc(n))

## # A tibble: 2 × 2
##   best_model     n
##   <fct>      <int>
## 1 crm           32
## 2 baseline      18

3.2.2 Object condition

3.2.2.1 Parameter search

rmse.exp1.object_uniform = data.frame()
for (p in seq(0, 1, by = 0.05)) {
  rmse.exp1.object_uniform = update_search_one(1, 'object', p)
}
colnames(rmse.exp1.object_uniform) = c('p', 'rmse')

# get optimal values
p.exp1.object = rmse.exp1.object_uniform[which.min(rmse.exp1.object_uniform$rmse),]$p

# create plot
plot_rmse.exp1.object = plot_search_one(1, 'object')

# update response data
data.exp1.object = data.exp1.object %>%
  mutate(prob = 1 - (1 - p.exp1.object)^num)

3.2.2.2 Overall models

Fitting baseline model:

fit.exp1.object_baseline = brm(
  formula = resp ~ 1 + (1 | id),
  data = data.exp1.object,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp1.object_baseline")

Fitting counterfactual replacement model (CRM):

fit.exp1.object_crm = brm(
  formula = resp ~ 1 + prob +  (1 + prob | id),
  data = data.exp1.object,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp1.object_crm")

Comparing models with leave-one-out cross-validation:

fit.exp1.object_baseline = add_criterion(
  fit.exp1.object_baseline,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp1.object_baseline")

fit.exp1.object_crm = add_criterion(
  fit.exp1.object_crm,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp1.object_crm")

loo_compare(fit.exp1.object_baseline,
            fit.exp1.object_crm)

##                          elpd_diff se_diff
## fit.exp1.object_crm          0.0       0.0
## fit.exp1.object_baseline -1778.2      80.2

3.2.2.3 Individual models

prior = c(prior(normal(0, 20), class = 'b', ub = 0))

# initialize the models 
fit.exp1.object_baseline_individual = brm(
  formula = resp ~ 1,
  data = data.exp1.object %>% 
    filter(plot_id == 1),
  seed = 1,
  file = "cache/fit.exp1.object_baseline_individual")

fit.exp1.object_crm_individual = brm(
  formula = resp ~ 1 + prob,
  data = data.exp1.object %>% 
    filter(plot_id == 1),
  prior = prior,
  seed = 1,
  file = "cache/fit.exp1.object_crm_individual")

# update model fits for each participant 
data.exp1.object_model_fits = data.exp1.object %>% 
  group_by(plot_id) %>% 
  nest() %>% 
  ungroup() %>% 
  mutate(fit_baseline = map(.x = data,
                            .f = ~ update(fit.exp1.object_baseline_individual,
                                          newdata = .x,
                                          seed = 1)),
         fit_crm = map(.x = data,
                       .f = ~ update(fit.exp1.object_crm_individual,
                                     newdata = .x,
                                     seed = 1))) %>%
  mutate(fit_baseline = map(.x = fit_baseline,
                            .f = ~ add_criterion(.x, criterion = "loo")),
         fit_crm = map(.x = fit_crm,
                       .f = ~ add_criterion(.x, criterion = "loo")),
         model_comparison = pmap(.l = list(baseline = fit_baseline,
                                           crm = fit_crm),
                                 .f = ~ loo_compare(..1, ..2)),
         best_model = map_chr(.x = model_comparison,
                              .f = ~ rownames(.) %>% 
                                .[1]),
         best_model = factor(best_model,
                             levels = c("..1", "..2"),
                             labels = c("baseline",
                                        "crm")))

save(list = c("data.exp1.object_model_fits"),
     file = paste0(data_dir, 'experiment1_object/model-fits.RData'))

load(file = paste0(data_dir, 'experiment1_object/model-fits.RData'))

data.exp1.object_model_fits %>% 
  count(best_model) %>%
  arrange(desc(n))

## # A tibble: 2 × 2
##   best_model     n
##   <fct>      <int>
## 1 crm           37
## 2 baseline      13

3.3 Plots

3.3.1 Parameter search

plot_rmse.exp1.agent + plot_spacer() + plot_rmse.exp1.object +
  plot_annotation(tag_levels = 'A') +
  plot_layout(widths = c(4, 1, 4)) &
  theme(plot.tag = element_text(size = 24))

# ggsave(paste0(figures_dir, '/experiment1/p_search.pdf'), width = 9, height = 4)

3.3.2 Setting up data

data.exp1.agent = data.exp1.agent %>%
  cbind(fit.exp1.agent_baseline %>%
          fitted(newdata = data.exp1.agent,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_baseline = Estimate) %>%
          select(pred_baseline)) %>%
  cbind(fit.exp1.agent_crm %>%
          fitted(newdata = data.exp1.agent,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_crm = Estimate) %>%
          select(pred_crm))

data.exp1.agent_means = data.exp1.agent %>%
  group_by(trial, craft) %>%
  summarise(resp_mean = mean(resp),
            resp_low = smean.cl.boot(resp)[2],
            resp_high = smean.cl.boot(resp)[3]) %>%
  merge(info.exp1.agent, by = c('trial', 'craft')) %>%
  left_join(data.exp1.agent %>%
              distinct(trial, craft, pred_baseline, pred_crm),
            by = c('trial', 'craft'))

data.exp1.object = data.exp1.object  %>%
  cbind(fit.exp1.object_baseline %>%
          fitted(newdata = data.exp1.object ,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_baseline = Estimate) %>%
          select(pred_baseline)) %>%
  cbind(fit.exp1.object_crm %>%
          fitted(newdata = data.exp1.object ,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_crm = Estimate) %>%
          select(pred_crm))

data.exp1.object_means = data.exp1.object  %>%
  group_by(trial, color) %>%
  summarise(resp_mean = mean(resp),
            resp_low = smean.cl.boot(resp)[2],
            resp_high = smean.cl.boot(resp)[3]) %>%
  merge(info.exp1.object, by = c('trial', 'color')) %>%
  left_join(data.exp1.object  %>%
              distinct(trial, color, pred_baseline, pred_crm),
            by = c('trial', 'color'))

model_labels.exp1 = c('baseline' = 'Contribution model',
                      'crm' = 'CRM')

model_order.exp1 = c('baseline', 'crm')

3.3.3 Scatter plots

scatter.exp1.agent_crm = make_scatter(
  data.exp1.agent_means, 'crm', 'CRM', seq(0, 100, 20))
scatter.exp1.object_crm = make_scatter(
  data.exp1.object_means, 'crm', 'CRM', seq(0, 100, 20))

scatter.exp1.agent_crm + inset_element(agent_icon, 0.85, 0.01, 1, 0.25) +
  scatter.exp1.object_crm + inset_element(object_icon, 0.85, 0.01, 1, 0.15) +
  plot_layout(ncol = 2) +
  plot_annotation(tag_levels = list(c('A', '', 'B', ''))) &
  theme(plot.title = element_text(hjust = 0.5),
        plot.tag = element_text(size = 24))

# ggsave(paste0(figures_dir, 'experiment1/scatter.pdf'), width = 8, height = 4)

3.3.4 Bar plot

barplot.exp1 = function(data, fill) {
  g = ggplot(data = data,
             mapping = aes(x = num,
                           y = resp)) +
    stat_summary(fun = "mean",
                 geom = "bar",
                 color = 'black',
                 fill = fill,
                 width = 0.8) +
    stat_summary(fun.data = "mean_cl_boot",
                 geom = "linerange",
                 size = 1) +
    stat_summary(mapping = aes(y = pred,
                               fill = model,
                               shape = model,
                               size = model),
                 fun = 'mean',
                 geom = 'point',
                 color = 'black',
                 position = position_dodge(width = 0.5)) +
    theme(text = element_text(size = 16),
          legend.position = 'bottom',
          legend.direction = 'horizontal',
          legend.title = element_blank(),
          #axis.text.x = element_blank(),
          axis.line.x = element_blank(),
          axis.ticks.x = element_blank(),
          strip.background = element_blank(),
          strip.text = element_blank(),
          panel.grid.major.y = element_line()) +
    scale_fill_manual(values = c('baseline' = 'black',
                                 'crm' = 'white'),
                      labels = model_labels.exp1,
                      name = 'model') +
    scale_size_manual(values = c('baseline' = 2.5,
                                 'crm' = 3),
                      labels = model_labels.exp1,
                      name = 'model') +
    scale_shape_manual(values = c('baseline' = 21,  # circle
                                  'crm' = 23), # diamond
                       labels = model_labels.exp1,
                       name = 'model') +
    scale_y_continuous(name = 'responsibility judgment',
                       breaks = seq(0, 100, 20),
                       expand = expansion(mult = 0.01)) +
    xlab('number of replacements') +
    coord_cartesian(clip = 'off', ylim = c(0, 100))
  
  return(g)
}

barplot.exp1.agent = barplot.exp1(
  data = data.exp1.agent %>%
    pivot_longer(cols = starts_with('pred'),
                 names_to = 'model',
                 values_to = 'pred') %>%
    mutate(model = factor(sub('^pred_', '', model),
                          levels = model_order.exp1)),
  fill = 'gray90'
)

barplot.exp1.object = barplot.exp1(
  data = data.exp1.object %>%
    pivot_longer(cols = starts_with('pred'),
                 names_to = 'model',
                 values_to = 'pred') %>%
    mutate(model = factor(sub('^pred_', '', model),
                          levels = model_order.exp1)),
  fill = 'gray75'
)

barplot.exp1.agent + inset_element(agent_icon, 0.85, 0.8, 1, 0.95,
                                   align_to = 'full') +
  barplot.exp1.object + inset_element(object_icon, 0.85, 0.9, 1, 1) +
  plot_layout(nrow = 1,
              guides = 'collect') +
  plot_annotation(tag_levels = list(c('A', '', 'B', ''))) &
  theme(plot.tag = element_text(size = 24),
        legend.position = 'bottom')

# ggsave(paste0(figures_dir, 'experiment1/bar.pdf'), width = 10, height = 4.5)

4 Experiment 2

4.1 Data

4.1.1 Agent condition

4.1.1.1 Trial info

info.exp2.agent = read_csv(paste0(data_dir, 'experiment2_agent/trial_info.csv'),
                           show_col_types = F) %>%
  pivot_longer(cols = !trial,
               names_to = c('.value', 'craft'),
               names_pattern = '(.+)_(.)') %>%
  mutate(num = low + high)

4.1.1.2 Response data

data.exp2.agent = read_csv(paste0(data_dir, 'experiment2_agent/trials.csv'),
                           show_col_types = F) %>%
  select(-error, -proliferate.condition) %>%
  rename(id = workerid,
         b = resp_b, c = resp_c, t = resp_t) %>%
  # filter by exclusion trial X1
  pivot_wider(names_from = trial, values_from = c(c, b, t)) %>%
  rowwise(id) %>%
  # pivot back to long form
  pivot_longer(cols = !id,
               names_to = 'trial',
               values_to = 'resp') %>%
  separate(col = trial,
           into = c('craft', 'trial'),
           sep = '_') %>%
  # ignore exclusion trial and practice trials
  filter(!(trial %in% c("X1", 1, 2))) %>%
  # merge with trial info
  merge(info.exp2.agent,
        by = c('trial', 'craft')) %>%
  # factor variables for plotting
  mutate(id = factor(id),
         plot_id = as.factor(as.integer(id)),
         craft = factor(craft))

Number of participants:

length(unique(data.exp2.agent$id))

## [1] 50

4.1.1.3 Participant data

Participants were asked: “What factors influenced how you decided to respond? Do you have any questions or comments regarding the experiment?”

participants.exp2.agent = read_csv(paste0(data_dir, 'experiment2_agent/participants.csv'),
                                   show_col_types = F) %>%
  select(-error, -proliferate.condition) %>%
  rename(id = workerid) %>%
  filter(id %in% data.exp2.agent$id) %>%
  mutate(id = factor(id)) %>%
  left_join(data.exp2.agent %>%
              distinct(id, plot_id),
            by = 'id')

participants.exp2.agent %>%
  select(plot_id, feedback) %>% 
  datatable()

4.1.2 Object condition

4.1.2.1 Trial info

info.exp2.object = read_csv(paste0(data_dir, 'experiment2_object/trial_info.csv'),
                            show_col_types = F) %>%
  pivot_longer(cols = !trial,
               names_to = c('.value', 'color'),
               names_pattern = '(.+)_(.)') %>%
  mutate(num = low + high)

4.1.2.2 Response data

data.exp2.object = read_csv(paste0(data_dir, 'experiment2_object/trials.csv'),
                            show_col_types = F) %>%
  select(-error, -proliferate.condition) %>%
  rename(id = workerid,
         y = resp_y, g = resp_g, b = resp_b) %>%
  # filter by exclusion trial X1
  pivot_wider(names_from = trial, values_from = c(y, g, b)) %>%
  rowwise(id) %>%
  filter(max(y_X1, g_X1, b_X1) - min(y_X1, g_X1, b_X1) < 30) %>%
  # pivot back to long form
  pivot_longer(cols = !id,
               names_to = 'trial',
               values_to = 'resp') %>%
  separate(col = trial,
           into = c('color', 'trial'),
           sep = '_') %>%
  # ignore exclusion trial and practice trials
  filter(!(trial %in% c("X1", 1, 2))) %>%
  # merge with trial info
  merge(info.exp2.object,
        by = c('trial', 'color')) %>%
  # factor variables for plotting
  mutate(id = factor(id),
         plot_id = as.factor(as.integer(id)),
         color = factor(color))

Number of participants:

length(unique(data.exp2.object$id))

## [1] 50

4.1.2.3 Participant data

Participants were asked: “What factors influenced how you decided to respond? Do you have any questions or comments regarding the experiment?”

participants.exp2.object = read_csv(paste0(data_dir, 'experiment2_object/participants.csv'),
                                    show_col_types = F) %>%
  select(-error, -proliferate.condition) %>%
  rename(id = workerid) %>%
  filter(id %in% data.exp2.object$id) %>%
  mutate(id = factor(id)) %>%
  left_join(data.exp2.object %>%
              distinct(id, plot_id),
            by = 'id')

participants.exp2.object %>%
  select(plot_id, feedback) %>% 
  datatable()

4.1.3 Demographics

participants.exp2.agent %>%
  rbind(participants.exp2.object) %>%
  print_demographics()

## # A tibble: 3 × 2
##   gender         n
##   <chr>      <int>
## 1 Female        40
## 2 Male          58
## 3 Non-binary     2
## [1] "age:"
## [1] 24.86
## [1] 6.446391
## # A tibble: 6 × 2
##   race                              n
##   <chr>                         <int>
## 1 American Indian/Alaska Native     2
## 2 Asian                             5
## 3 Black                             7
## 4 Multiracial                       3
## 5 White                            58
## 6 <NA>                             25
## # A tibble: 3 × 2
##   ethnicity        n
##   <chr>        <int>
## 1 Hispanic        40
## 2 Non-Hispanic    58
## 3 <NA>             2
## [1] "time taken (min):"
## [1] 12.33368
## [1] 6.544082

4.2 Models

4.2.1 Agent condition

4.2.1.1 Parameter search

Uniform probabilities:

rmse.exp2.agent_uniform = data.frame()
for (p in seq(0, 1, by = 0.05)) {
  rmse.exp2.agent_uniform = update_search_one(2, 'agent', p)
}
colnames(rmse.exp2.agent_uniform) = c('p', 'rmse')

# get optimal values
p.exp2.agent = rmse.exp2.agent_uniform[which.min(rmse.exp2.agent_uniform$rmse),]$p

# create plot
plot_rmse.exp2.agent_uniform = plot_search_one(2, 'agent')

# update response data
data.exp2.agent = data.exp2.agent %>%
  mutate(prob_uniform = 1 - (1 - p.exp2.agent)^num)

Varying probabilities:

rmse.exp2.agent_full = data.frame()
for (p_low in seq(0, 1, by = 0.05)) {
  for (p_high in seq(p_low, 1, by = 0.05)) {
    rmse.exp2.agent_full = update_search_two(2, 'agent', p_low, p_high)
  }
}
colnames(rmse.exp2.agent_full) = c('p_low', 'p_high', 'rmse')

# get optimal values
p.exp2.agent_low = rmse.exp2.agent_full[which.min(rmse.exp2.agent_full$rmse),]$p_low
p.exp2.agent_high = rmse.exp2.agent_full[which.min(rmse.exp2.agent_full$rmse),]$p_high

# create plot
plot_rmse.exp2.agent_full = plot_search_two(2, 'agent')

# update response data
data.exp2.agent = data.exp2.agent %>%
  mutate(prob_full = 1 - (1 - p.exp2.agent_low)^low * (1 - p.exp2.agent_high)^high)

4.2.1.2 Overall models

Fitting baseline model:

fit.exp2.agent_baseline = brm(
  formula = resp ~ 1 + (1 | id),
  data = data.exp2.agent,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp2.agent_baseline")

Fitting CRM (uniform):

fit.exp2.agent_crm_uniform = brm(
  formula = resp ~ 1 + prob_uniform + (1 + prob_uniform | id),
  data = data.exp2.agent,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp2.agent_crm_uniform")

Fitting CRM (full):

fit.exp2.agent_crm_full = brm(
  formula = resp ~ 1 + prob_full + (1 + prob_full | id),
  data = data.exp2.agent,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp2.agent_crm_full")

Comparing models with leave-one-out cross-validation:

fit.exp2.agent_baseline = add_criterion(
  fit.exp2.agent_baseline,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp2.agent_baseline")

fit.exp2.agent_crm_uniform = add_criterion(
  fit.exp2.agent_crm_uniform,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp2.agent_crm_uniform")

fit.exp2.agent_crm_full = add_criterion(
  fit.exp2.agent_crm_full,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp2.agent_crm_full")

loo_compare(fit.exp2.agent_baseline,
            fit.exp2.agent_crm_uniform,
            fit.exp2.agent_crm_full)

##                            elpd_diff se_diff
## fit.exp2.agent_crm_full       0.0       0.0 
## fit.exp2.agent_crm_uniform -411.7      35.9 
## fit.exp2.agent_baseline    -586.4      40.2

4.2.1.3 Individual models

prior = c(prior(normal(0, 20), class = 'b', ub = 0))

fit.exp2.agent_baseline_individual = brm(
  formula = resp ~ 1,
  data = data.exp2.agent %>% 
    filter(plot_id == 1),
  seed = 1,
  file = "cache/fit.exp2.agent_baseline_individual")

fit.exp2.agent_crm_uniform_individual = brm(
  formula = resp ~ 1 + prob_uniform,
  data = data.exp2.agent %>% 
    filter(plot_id == 1),
  prior = prior,
  seed = 1,
  file = "cache/fit.exp2.agent_crm_uniform_individual")

fit.exp2.agent_crm_full_individual = brm(
  formula = resp ~ 1 + prob_full,
  data = data.exp2.agent %>% 
    filter(plot_id == 1),
  prior = prior,
  seed = 1,
  file = "cache/fit.exp2.agent_crm_full_individual")

# update model fits for each participant 
data.exp2.agent_model_fits = data.exp2.agent %>% 
  group_by(plot_id) %>% 
  nest() %>% 
  ungroup() %>% 
  mutate(fit_baseline = map(.x = data,
                            .f = ~ update(fit.exp2.agent_baseline_individual,
                                          newdata = .x,
                                          seed = 1)),
         fit_crm_uniform = map(.x = data,
                               .f = ~ update(fit.exp2.agent_crm_uniform_individual,
                                             newdata = .x,
                                             seed = 1)),
         fit_crm_full = map(.x = data,
                            .f = ~ update(fit.exp2.agent_crm_full_individual,
                                          newdata = .x,
                                          seed = 1))) %>%
  mutate(fit_baseline = map(.x = fit_baseline,
                            .f = ~ add_criterion(.x, criterion = "loo")),
         fit_crm_uniform = map(.x = fit_crm_uniform,
                               .f = ~ add_criterion(.x, criterion = "loo")),
         fit_crm_full = map(.x = fit_crm_full,
                            .f = ~ add_criterion(.x, criterion = "loo")),
         model_comparison = pmap(.l = list(baseline = fit_baseline,
                                           crm_uniform = fit_crm_uniform,
                                           crm_full = fit_crm_full),
                                 .f = ~ loo_compare(..1, ..2, ..3)),
         best_model = map_chr(.x = model_comparison,
                              .f = ~ rownames(.) %>% 
                                .[1]),
         best_model = factor(best_model,
                             levels = c("..1", "..2", "..3"),
                             labels = c("baseline",
                                        "crm_uniform",
                                        "crm_full")))

save(list = c("data.exp2.agent_model_fits"),
     file = paste0(data_dir, 'experiment2_agent/model-fits.RData'))

load(file = paste0(data_dir, 'experiment2_agent/model-fits.RData'))

data.exp2.agent_model_fits %>% 
  count(best_model) %>%
  arrange(desc(n))

## # A tibble: 3 × 2
##   best_model      n
##   <fct>       <int>
## 1 crm_full       26
## 2 baseline       17
## 3 crm_uniform     7

4.2.2 Object condition

4.2.2.1 Parameter search

Uniform probabilities:

rmse.exp2.object_uniform = data.frame()
for (p in seq(0, 1, by = 0.05)) {
  rmse.exp2.object_uniform = update_search_one(2, 'object', p)
}
colnames(rmse.exp2.object_uniform) = c('p', 'rmse')

# get optimal values
p.exp2.object = rmse.exp2.object_uniform[which.min(rmse.exp2.object_uniform$rmse),]$p

# create plot
plot_rmse.exp2.object_uniform = plot_search_one(2, 'object')

# update response data
data.exp2.object = data.exp2.object %>%
  mutate(prob_uniform = 1 - (1 - p.exp2.object)^num)

Varying probabilities:

rmse.exp2.object_full = data.frame()
for (p_low in seq(0, 1, by = 0.05)) {
  for (p_high in seq(p_low, 1, by = 0.05)) {
    rmse.exp2.object_full = update_search_two(2, 'object', p_low, p_high)
  }
}
colnames(rmse.exp2.object_full) = c('p_low', 'p_high', 'rmse')

# get optimal values
p.exp2.object_low = rmse.exp2.object_full[which.min(rmse.exp2.object_full$rmse),]$p_low
p.exp2.object_high = rmse.exp2.object_full[which.min(rmse.exp2.object_full$rmse),]$p_high

# create plot
plot_rmse.exp2.object_full = plot_search_two(2, 'object')

# update response data
data.exp2.object = data.exp2.object %>%
  mutate(prob_full = 1 - (1 - p.exp2.object_low)^low * (1 - p.exp2.object_high)^high)

4.2.2.2 Overall models

Fitting baseline model:

fit.exp2.object_baseline = brm(
  formula = resp ~ 1 + (1 | id),
  data = data.exp2.object,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp2.object_baseline")

Fitting CRM (uniform):

fit.exp2.object_crm_uniform = brm(
  formula = resp ~ 1 + prob_uniform + (1 + prob_uniform | id),
  data = data.exp2.object,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp2.object_crm_uniform")

Fitting CRM (full):

fit.exp2.object_crm_full = brm(
  formula = resp ~ 1 + prob_full + (1 + prob_full | id),
  data = data.exp2.object,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp2.object_crm_full")

Comparing models with leave-one-out cross-validation:

fit.exp2.object_baseline = add_criterion(
  fit.exp2.object_baseline,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp2.object_baseline")

fit.exp2.object_crm_uniform = add_criterion(
  fit.exp2.object_crm_uniform,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp2.object_crm_uniform")

fit.exp2.object_crm_full = add_criterion(
  fit.exp2.object_crm_full,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp2.object_crm_full")

loo_compare(fit.exp2.object_baseline,
            fit.exp2.object_crm_uniform,
            fit.exp2.object_crm_full)

##                             elpd_diff se_diff
## fit.exp2.object_crm_full       0.0       0.0 
## fit.exp2.object_crm_uniform  -61.8       6.8 
## fit.exp2.object_baseline    -307.3      24.8

4.2.2.3 Individual models

prior = c(prior(normal(0, 20), class = 'b', ub = 0))

fit.exp2.object_baseline_individual = brm(
  formula = resp ~ 1,
  data = data.exp2.object %>% 
    filter(plot_id == 1),
  seed = 1,
  file = "cache/fit.exp2.object_baseline_individual")

fit.exp2.object_crm_uniform_individual = brm(
  formula = resp ~ 1 + prob_uniform,
  data = data.exp2.object %>% 
    filter(plot_id == 1),
  prior = prior,
  seed = 1,
  file = "cache/fit.exp2.object_crm_uniform_individual")

fit.exp2.object_crm_full_individual = brm(
  formula = resp ~ 1 + prob_full,
  data = data.exp2.object %>% 
    filter(plot_id == 1),
  prior = prior,
  seed = 1,
  file = "cache/fit.exp2.object_crm_full_individual")

# update model fits for each participant 
data.exp2.object_model_fits = data.exp2.object %>% 
  group_by(plot_id) %>% 
  nest() %>% 
  ungroup() %>% 
  mutate(fit_baseline = map(.x = data,
                            .f = ~ update(fit.exp2.object_baseline_individual,
                                          newdata = .x,
                                          seed = 1)),
         fit_crm_uniform = map(.x = data,
                               .f = ~ update(fit.exp2.object_crm_uniform_individual,
                                             newdata = .x,
                                             seed = 1)),
         fit_crm_full = map(.x = data,
                            .f = ~ update(fit.exp2.object_crm_full_individual,
                                          newdata = .x,
                                          seed = 1))) %>%
  mutate(fit_baseline = map(.x = fit_baseline,
                            .f = ~ add_criterion(.x, criterion = "loo")),
         fit_crm_uniform = map(.x = fit_crm_uniform,
                               .f = ~ add_criterion(.x, criterion = "loo")),
         fit_crm_full = map(.x = fit_crm_full,
                            .f = ~ add_criterion(.x, criterion = "loo")),
         model_comparison = pmap(.l = list(baseline = fit_baseline,
                                           crm_uniform = fit_crm_uniform,
                                           crm_full = fit_crm_full),
                                 .f = ~ loo_compare(..1, ..2, ..3)),
         best_model = map_chr(.x = model_comparison,
                              .f = ~ rownames(.) %>% 
                                .[1]),
         best_model = factor(best_model,
                             levels = c("..1", "..2", "..3"),
                             labels = c("baseline",
                                        "crm_uniform",
                                        "crm_full")))

save(list = c("data.exp2.object_model_fits"),
     file = paste0(data_dir, 'experiment2_object/model-fits.RData'))

load(file = paste0(data_dir, 'experiment2_object/model-fits.RData'))

data.exp2.object_model_fits %>% 
  count(best_model) %>%
  arrange(desc(n))

## # A tibble: 3 × 2
##   best_model      n
##   <fct>       <int>
## 1 crm_full       21
## 2 baseline       16
## 3 crm_uniform    13

4.3 Plots

4.3.1 Parameter search

plot_rmse.exp2.agent_uniform + plot_spacer() + plot_rmse.exp2.object_uniform +
  plot_annotation(tag_levels = 'A') +
  plot_layout(widths = c(4, 1, 4)) &
  theme(plot.tag = element_text(size = 24))

# ggsave(paste0(figures_dir, 'experiment2/p_search_uniform.pdf'), width = 9, height = 4)

plot_rmse.exp2.agent_full + plot_spacer() + plot_rmse.exp2.object_full +
  plot_annotation(tag_levels = 'A') +
  plot_layout(widths = c(4, 1, 4)) &
  theme(plot.tag = element_text(size = 24))

# ggsave(paste0(figures_dir, 'experiment2/p_search_full.pdf'), width = 9, height = 4)

4.3.2 Setting up data

data.exp2.agent = data.exp2.agent %>%
  cbind(fit.exp2.agent_baseline %>%
          fitted(newdata = data.exp2.agent,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_baseline = Estimate) %>%
          select(pred_baseline)) %>%
  cbind(fit.exp2.agent_crm_uniform %>%
          fitted(newdata = data.exp2.agent,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_crm_uniform = Estimate) %>%
          select(pred_crm_uniform)) %>%
  cbind(fit.exp2.agent_crm_full %>%
          fitted(newdata = data.exp2.agent,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_crm_full = Estimate) %>%
          select(pred_crm_full))

data.exp2.agent_means = data.exp2.agent %>%
  group_by(trial, craft) %>%
  summarise(resp_mean = mean(resp),
            resp_low = smean.cl.boot(resp)[2],
            resp_high = smean.cl.boot(resp)[3]) %>%
  merge(info.exp2.agent, by = c('trial', 'craft')) %>%
  left_join(data.exp2.agent %>%
              distinct(trial, craft, pred_baseline,
                       pred_crm_uniform, pred_crm_full),
            by = c('trial', 'craft'))

data.exp2.object = data.exp2.object  %>%
  cbind(fit.exp2.object_baseline %>%
          fitted(newdata = data.exp2.object ,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_baseline = Estimate) %>%
          select(pred_baseline)) %>%
  cbind(fit.exp2.object_crm_uniform %>%
          fitted(newdata = data.exp2.object ,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_crm_uniform = Estimate) %>%
          select(pred_crm_uniform)) %>%
  cbind(fit.exp2.object_crm_full %>%
          fitted(newdata = data.exp2.object ,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_crm_full = Estimate) %>%
          select(pred_crm_full))

data.exp2.object_means = data.exp2.object  %>%
  group_by(trial, color) %>%
  summarise(resp_mean = mean(resp),
            resp_low = smean.cl.boot(resp)[2],
            resp_high = smean.cl.boot(resp)[3]) %>%
  merge(info.exp2.object, by = c('trial', 'color')) %>%
  left_join(data.exp2.object  %>%
              distinct(trial, color, pred_baseline,
                       pred_crm_uniform, pred_crm_full),
            by = c('trial', 'color'))

labels.exp2 = info.exp2.agent %>% 
  distinct(low, high, num) %>% 
  arrange(num, high) %>%
  mutate(condition = NA,
         index = 1:n(),
         grob = map2(.x = low,
                     .y = high,
                     .f = ~ readPNG(paste0(figures_dir, 'experiment2/sets/',
                                           .x, .y, '.png'))))

model_labels.exp2 = c('baseline' = 'Contribution model',
                      'crm_uniform' = 'CRM (uniform)',
                      'crm_full' = 'CRM (full)')

model_order.exp2 = c('baseline', 'crm_uniform', 'crm_full')

4.3.3 Scatter plots

scatter.exp2.agent_uniform = make_scatter(
  data.exp2.agent_means, 'crm_uniform', 'CRM (uniform)', seq(40, 100, 20))
scatter.exp2.agent_full = make_scatter(
  data.exp2.agent_means, 'crm_full', 'CRM (full)', seq(40, 100, 20))

scatter.exp2.object_uniform = make_scatter(
  data.exp2.object_means, 'crm_uniform', 'CRM (uniform)', seq(40, 100, 20))
scatter.exp2.object_full = make_scatter(
  data.exp2.object_means, 'crm_full', 'CRM (full)', seq(40, 100, 20))

scatter.exp2.agent_uniform + inset_element(agent_icon, 0.85, 0.01, 1, 0.25) +
  scatter.exp2.agent_full + inset_element(agent_icon, 0.85, 0.01, 1, 0.25) +
  scatter.exp2.object_uniform + inset_element(object_icon, 0.85, 0.01, 1, 0.15) +
  scatter.exp2.object_full + inset_element(object_icon, 0.85, 0.01, 1, 0.15) +
  plot_layout(ncol = 2) +
  plot_annotation(tag_levels = list(c('A', '', '', '',
                                      'B', '', '', ''))) &
  theme(plot.title = element_text(hjust = 0.5),
        plot.tag = element_text(size = 24))

# ggsave(paste0(figures_dir, 'experiment2/scatter.pdf'), width = 8, height = 8)

4.3.4 Bar plots

make_bar.exp2 = function(data, fill) {
  g = ggplot(data = data,
             mapping = aes(x = index,
                           y = resp)) +
    stat_summary(fun = "mean",
                 geom = "bar",
                 color = 'black',
                 fill = fill,
                 width = 0.8) +
    stat_summary(fun.data = "mean_cl_boot",
                 geom = "linerange",
                 size = 1) +
    stat_summary(mapping = aes(y = pred,
                               fill = model,
                               shape = model,
                               size = model),
                 fun = 'mean',
                 geom = 'point',
                 color = 'black',
                 position = position_dodge(width = 0.5)) +
    facet_wrap(~ index, scales = 'free_x', nrow = 1) +
    geom_rect(aes(xmin = index - 0.5, xmax = index + 0.5,
                  ymin = 0, ymax = 39.9),
              fill = 'white') +
    geom_custom(data = labels.exp2,
                mapping = aes(data = grob,
                              x = -Inf, y = 36),
                grob_fun = function(x) rasterGrob(x,
                                                  interpolate = T,
                                                  vjust = 0,
                                                  hjust = 0)) +
    geom_text(aes(label = index,
                  y = 33),
              show.legend = F,
              check_overlap = T) +
    theme(text = element_text(size = 16),
          legend.position = c(0.5, 0.87),
          legend.direction = 'horizontal',
          legend.title = element_blank(),
          axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.line.x = element_blank(),
          axis.ticks.x = element_blank(),
          strip.background = element_blank(),
          strip.text = element_blank(),
          panel.grid.major.y = element_line(),
          panel.spacing.x = unit(0.1, "cm"),
          plot.margin = margin(t = 0.2, l = 0.2, r = 0.2, b = 0.9, unit = "cm")) +
    scale_fill_manual(values = c('baseline' = 'black',
                                 'crm_uniform' = 'gray50',
                                 'crm_full' = 'white'),
                      labels = model_labels.exp2,
                      name = 'model') +
    scale_size_manual(values = c('baseline' = 1.5,
                                 'crm_uniform' = 1.5,
                                 'crm_full' = 2),
                      labels = model_labels.exp2,
                      name = 'model') +
    scale_shape_manual(values = c('baseline' = 21,  # circle
                                  'crm_uniform' = 24, # triangle
                                  'crm_full' = 23), # diamond
                       labels = model_labels.exp2,
                       name = 'model') +
    scale_y_continuous(name = 'responsibility judgment',
                       breaks = seq(40, 100, 20)) +
    coord_cartesian(clip = 'off', ylim = c(40, 100))
  
  return(g)
}

bar.exp2.agent = make_bar.exp2(
  data = data.exp2.agent %>%
    left_join(labels.exp2 %>%
                select(low, high, index),
              by = c('low', 'high')) %>%
    pivot_longer(cols = starts_with('pred'),
                 names_to = 'model',
                 values_to = 'pred') %>%
    mutate(model = factor(sub('^pred_', '', model),
                          levels = model_order.exp2)),
  fill = 'gray90' 
)

bar.exp2.object = make_bar.exp2(
  data = data.exp2.object %>%
    left_join(labels.exp2 %>%
                select(low, high, index),
              by = c('low', 'high')) %>%
    pivot_longer(cols = starts_with('pred'),
                 names_to = 'model',
                 values_to =  'pred') %>%
    mutate(model = factor(sub('^pred_', '', model),
                          levels = model_order.exp2)),
  fill = 'gray75'
)

bar.exp2.agent + inset_element(agent_icon, 0.85, 0.7, 1, 1) +
  bar.exp2.object + inset_element(object_icon, 0.85, 0.8, 1, 1) +
  plot_layout(ncol = 1) +
  plot_annotation(tag_levels = list(c('A', '', 'B', ''))) &
  theme(plot.tag = element_text(size = 24))

# ggsave(paste0(figures_dir, 'experiment2/bar.pdf'), width = 10, height = 8)

5 Experiment 3

5.1 Data

5.1.1 Agent condition

5.1.1.1 Trial info

info.exp3.agent = read_csv(paste0(data_dir, 'experiment3_agent/trial_info.csv'),
                           show_col_types = F) %>%
  pivot_longer(cols = !trial,
               names_to = c('.value', 'craft'),
               names_pattern = '(.+)_(.)') %>%
  mutate(num = low + high)

5.1.1.2 Response data

data.exp3.agent = read_csv(paste0(data_dir, 'experiment3_agent/trials.csv'),
                           show_col_types = F)

# filter participants who accidentally responded more than once
rm_ids_agent = data.exp3.agent %>%
  group_by(id) %>%
  summarise(n = n()/length(unique(info.exp3.agent$trial))) %>%
  filter(n > 1) %>%
  .$id

data.exp3.agent = data.exp3.agent %>%
  filter(!(id %in% rm_ids_agent)) %>%
  pivot_longer(cols = c(c, t),
               names_to = 'craft',
               values_to = 'resp') %>%
  # ignore exclusion trial and practice trials
  filter(!(trial %in% c("X1", 1, 2))) %>%
  # merge with trial info
  merge(info.exp3.agent,
        by = c('trial', 'craft')) %>%
  # factor variables for plotting
  mutate(id = factor(id),
         plot_id = as.factor(as.integer(id)),
         craft = factor(craft))

Number of participants:

length(unique(data.exp3.agent$id))

## [1] 50

5.1.1.3 Participant data

Participants were asked: “What factors influenced how you decided to respond? Do you have any questions or comments regarding the experiment?”

participants.exp3.agent = read_csv(paste0(data_dir, 'experiment3_agent/participants.csv'),
                                   show_col_types = F) %>%
  filter(id %in% data.exp3.agent$id) %>%
  mutate(id = factor(id)) %>%
  left_join(data.exp3.agent %>%
              distinct(id, plot_id),
            by = 'id')

participants.exp3.agent %>%
  select(plot_id, feedback) %>% 
  datatable()

5.1.2 Object condition

5.1.2.1 Trial info

info.exp3.object = read_csv(paste0(data_dir, 'experiment3_object/trial_info.csv'),
                            show_col_types = F) %>%
  pivot_longer(cols = !trial,
               names_to = c('.value', 'color'),
               names_pattern = '(.+)_(.)') %>%
  mutate(num = low + high)

5.1.2.2 Response data

data.exp3.object = read_csv(paste0(data_dir, 'experiment3_object/trials.csv'),
                            show_col_types = F)

# filter participants who accidentally responded more than once
rm_ids_object = data.exp3.object %>%
  group_by(id) %>%
  summarise(n = n()/length(unique(info.exp3.object$trial))) %>%
  filter(n > 1) %>%
  .$id

data.exp3.object = data.exp3.object %>%
  filter(!(id %in% rm_ids_object)) %>%
  pivot_longer(cols = c(b, y),
               names_to = 'color',
               values_to = 'resp') %>%
  # ignore exclusion trial and practice trials
  filter(!(trial %in% c("X1", 1, 2))) %>%
  # merge with trial info
  merge(info.exp3.object,
        by = c('trial', 'color')) %>%
  # factor variables for plotting
  mutate(id = factor(id),
         plot_id = as.factor(as.integer(id)),
         color = factor(color))

Number of participants:

length(unique(data.exp3.object$id))

## [1] 50

5.1.2.3 Participant data

Participants were asked: “What factors influenced how you decided to respond? Do you have any questions or comments regarding the experiment?”

participants.exp3.object = read_csv(paste0(data_dir, 'experiment3_object/participants.csv'),
                                    show_col_types = F) %>%
  filter(id %in% data.exp3.object$id) %>%
  mutate(id = factor(id)) %>%
  left_join(data.exp3.object %>%
              distinct(id, plot_id),
            by = 'id')

participants.exp3.object %>%
  select(plot_id, feedback) %>% 
  datatable()

5.1.3 Demographics

participants.exp3.agent %>%
  rbind(participants.exp3.object) %>%
  print_demographics()

## # A tibble: 3 × 2
##   gender     n
##   <chr>  <int>
## 1 Female    56
## 2 Male      43
## 3 <NA>       1
## [1] "age:"
## [1] 19.74
## [1] 1.060017
## # A tibble: 7 × 2
##   race                              n
##   <chr>                         <int>
## 1 American Indian/Alaska Native     1
## 2 Asian                            44
## 3 Black                             7
## 4 Multiracial                       6
## 5 Undisclosed/other                 3
## 6 White                            36
## 7 <NA>                              3
## # A tibble: 3 × 2
##   ethnicity        n
##   <chr>        <int>
## 1 Hispanic        21
## 2 Non-Hispanic    74
## 3 <NA>             5
## [1] "time taken (min):"
## [1] 14.64606
## [1] 65.07773

5.2 Models

5.2.1 Agent condition

5.2.1.1 Parameter search

rmse.exp3.agent_full = data.frame()
for (p_low in seq(0, 1, by = 0.05)) {
  for (p_high in seq(p_low, 1, by = 0.05)) {
    rmse.exp3.agent_full = update_search_two(3, 'agent', p_low, p_high)
  }
}
colnames(rmse.exp3.agent_full) = c('p_low', 'p_high', 'rmse')

# get optimal values
p.exp3.agent_low = rmse.exp3.agent_full[which.min(rmse.exp3.agent_full$rmse),]$p_low
p.exp3.agent_high = rmse.exp3.agent_full[which.min(rmse.exp3.agent_full$rmse),]$p_high

# create plot
plot_rmse.exp3.agent = plot_search_two(3, 'agent')

# update response data
data.exp3.agent = data.exp3.agent %>%
  mutate(prob = 1 - (1 - p.exp3.agent_low)^low * (1 - p.exp3.agent_high)^high,
         contributor = ifelse(availability == 'high', p.exp3.agent_high, p.exp3.agent_low),
         CP = (1 - contributor) * (1 - prob))

5.2.1.2 Overall models

Fitting baseline model:

fit.exp3.agent_baseline = brm(
  formula = resp ~ 1 + (1 | id),
  data = data.exp3.agent,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp3.agent_baseline")

Fitting CP (multiplicative) model:

fit.exp3.agent_cp = brm(
  formula = resp ~ 1 + CP + (1 + CP | id),
  data = data.exp3.agent,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp3.agent_cp")

Fitting CRM-based account (additive):

fit.exp3.agent_crm = brm(
  formula = resp ~ 1 + contributor + prob + (1 + contributor + prob | id),
  data = data.exp3.agent,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp3.agent_crm")

Comparing models with leave-one-out cross-validation:

fit.exp3.agent_baseline = add_criterion(
  fit.exp3.agent_baseline,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp3.agent_baseline")

fit.exp3.agent_cp = add_criterion(
  fit.exp3.agent_cp,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp3.agent_cp")

fit.exp3.agent_crm = add_criterion(
  fit.exp3.agent_crm,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp3.agent_crm")

loo_compare(fit.exp3.agent_baseline,
            fit.exp3.agent_crm,
            fit.exp3.agent_cp)

##                         elpd_diff se_diff
## fit.exp3.agent_crm         0.0       0.0 
## fit.exp3.agent_cp       -579.9      55.8 
## fit.exp3.agent_baseline -793.0      61.8

5.2.1.3 Individual models

fit.exp3.agent_baseline_individual = brm(
  formula = resp ~ 1,
  data = data.exp3.agent %>% 
    filter(plot_id == 1),
  seed = 1,
  file = "cache/fit.exp3.agent_baseline_individual")

fit.exp3.agent_cp_individual = brm(
  formula = resp ~ 1 + CP,
  data = data.exp3.agent %>% 
    filter(plot_id == 1),
  seed = 1,
  file = "cache/fit.exp3.agent_cp_individual")

fit.exp3.agent_crm_individual = brm(
  formula = resp ~ 1 + contributor + prob,
  data = data.exp3.agent %>% 
    filter(plot_id == 1),
  seed = 1,
  file = "cache/fit.exp3.agent_crm_individual")

# update model fits for each participant 
data.exp3.agent_model_fits = data.exp3.agent %>% 
  group_by(plot_id) %>% 
  nest() %>% 
  ungroup() %>% 
  mutate(fit_baseline = map(.x = data,
                            .f = ~ update(fit.exp3.agent_baseline_individual,
                                          newdata = .x,
                                          seed = 1)),
         fit_cp = map(.x = data,
                      .f = ~ update(fit.exp3.agent_cp_individual,
                                    newdata = .x,
                                    seed = 1)),
         fit_crm = map(.x = data,
                       .f = ~ update(fit.exp3.agent_crm_individual,
                                     newdata = .x,
                                     seed = 1))) %>%
  mutate(fit_baseline = map(.x = fit_baseline,
                            .f = ~ add_criterion(.x, criterion = "loo")),
         fit_cp = map(.x = fit_cp,
                      .f = ~ add_criterion(.x, criterion = "loo")),
         fit_crm = map(.x = fit_crm,
                       .f = ~ add_criterion(.x, criterion = "loo")),
         model_comparison = pmap(.l = list(baseline = fit_baseline,
                                           cp = fit_cp,
                                           crm = fit_crm),
                                 .f = ~ loo_compare(..1, ..2, ..3)),
         best_model = map_chr(.x = model_comparison,
                              .f = ~ rownames(.) %>% 
                                .[1]),
         best_model = factor(best_model,
                             levels = c("..1", "..2", "..3"),
                             labels = c("baseline",
                                        "cp",
                                        "crm")))

save(list = c("data.exp3.agent_model_fits"),
     file = paste0(data_dir, 'experiment3_agent/model-fits.RData'))

load(file = paste0(data_dir, 'experiment3_agent/model-fits.RData'))

data.exp3.agent_model_fits %>% 
  count(best_model) %>%
  arrange(desc(n))

## # A tibble: 3 × 2
##   best_model     n
##   <fct>      <int>
## 1 crm           30
## 2 baseline      10
## 3 cp            10

# count number of positive and negative posteriors on replaceability predictor
data.exp3.agent_model_fits %>%
  filter(best_model == 'crm') %>%
  .$fit_crm %>%
  sapply(function(x) fixef(x)[['prob', 'Estimate']]) %>%
  sign() %>%
  as_tibble() %>%
  count(value)

## # A tibble: 2 × 2
##   value     n
##   <dbl> <int>
## 1    -1    25
## 2     1     5

# count number of positive and negative posteriors on contributor predictor
data.exp3.agent_model_fits %>%
  filter(best_model == 'crm') %>%
  .$fit_crm %>%
  sapply(function(x) fixef(x)[['contributor', 'Estimate']]) %>%
  sign() %>%
  as_tibble() %>%
  count(value)

## # A tibble: 2 × 2
##   value     n
##   <dbl> <int>
## 1    -1    17
## 2     1    13

5.2.2 Object condition

5.2.2.1 Parameter search

rmse.exp3.object_full = data.frame()
for (p_low in seq(0, 1, by = 0.05)) {
  for (p_high in seq(p_low, 1, by = 0.05)) {
    rmse.exp3.object_full = update_search_two(3, 'object', p_low, p_high)
  }
}
colnames(rmse.exp3.object_full) = c('p_low', 'p_high', 'rmse')

# get optimal values
p.exp3.object_low = rmse.exp3.object_full[which.min(rmse.exp3.object_full$rmse),]$p_low
p.exp3.object_high = rmse.exp3.object_full[which.min(rmse.exp3.object_full$rmse),]$p_high

# create plot
plot_rmse.exp3.object = plot_search_two(3, 'object')

# update response data
data.exp3.object = data.exp3.object %>%
  mutate(prob = 1 - (1 - p.exp3.object_low)^low * (1 - p.exp3.object_high)^high,
         contributor = ifelse(availability == 'high', p.exp3.object_high, p.exp3.object_low),
         CP = (1 - contributor) * (1 - prob))

5.2.2.2 Overall models

Fitting baseline model:

fit.exp3.object_baseline = brm(
  formula = resp ~ 1 + (1 | id),
  data = data.exp3.object,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp3.object_baseline")

Fitting CP (multiplicative) model:

fit.exp3.object_cp = brm(
  formula = resp ~ 1 + CP + (1 + CP | id),
  data = data.exp3.object,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp3.object_cp")

Fitting CRM-based model (additive):

fit.exp3.object_crm = brm(
  formula = resp ~ 1 + contributor + prob + (1 + contributor + prob | id),
  data = data.exp3.object,
  iter = 8000,
  seed = 1,
  file = "cache/fit.exp3.object_crm")

Comparing models with leave-one-out cross-validation:

fit.exp3.object_baseline = add_criterion(
  fit.exp3.object_baseline,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp3.object_baseline")

fit.exp3.object_cp = add_criterion(
  fit.exp3.object_cp,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp3.object_cp")

fit.exp3.object_crm = add_criterion(
  fit.exp3.object_crm,
  criterion = "loo",
  reloo = T,
  file = "cache/fit.exp3.object_crm")

loo_compare(fit.exp3.object_baseline,
            fit.exp3.object_crm,
            fit.exp3.object_cp)

##                          elpd_diff se_diff
## fit.exp3.object_crm         0.0       0.0 
## fit.exp3.object_cp       -611.6      40.5 
## fit.exp3.object_baseline -710.7      42.5

5.2.2.3 Individual models

fit.exp3.object_baseline_individual = brm(
  formula = resp ~ 1,
  data = data.exp3.object %>% 
    filter(plot_id == 1),
  seed = 1,
  file = "cache/fit.exp3.object_baseline_individual")

fit.exp3.object_cp_individual = brm(
  formula = resp ~ 1 + CP,
  data = data.exp3.object %>% 
    filter(plot_id == 1),
  seed = 1,
  file = "cache/fit.exp3.object_cp_individual")

fit.exp3.object_crm_individual = brm(
  formula = resp ~ 1 + contributor + prob,
  data = data.exp3.object %>% 
    filter(plot_id == 1),
  seed = 1,
  file = "cache/fit.exp3.object_crm_individual")

# update model fits for each participant 
data.exp3.object_model_fits = data.exp3.object %>% 
  group_by(plot_id) %>% 
  nest() %>% 
  ungroup() %>% 
  mutate(fit_baseline = map(.x = data,
                            .f = ~ update(fit.exp3.object_baseline_individual,
                                          newdata = .x,
                                          seed = 1)),
         fit_cp = map(.x = data,
                      .f = ~ update(fit.exp3.object_cp_individual,
                                    newdata = .x,
                                    seed = 1)),
         fit_crm = map(.x = data,
                       .f = ~ update(fit.exp3.object_crm_individual,
                                     newdata = .x,
                                     seed = 1))) %>%
  mutate(fit_baseline = map(.x = fit_baseline,
                            .f = ~ add_criterion(.x, criterion = "loo")),
         fit_cp = map(.x = fit_cp,
                      .f = ~ add_criterion(.x, criterion = "loo")),
         fit_crm = map(.x = fit_crm,
                       .f = ~ add_criterion(.x, criterion = "loo")),
         model_comparison = pmap(.l = list(baseline = fit_baseline,
                                           cp = fit_cp,
                                           crm = fit_crm),
                                 .f = ~ loo_compare(..1, ..2, ..3)),
         best_model = map_chr(.x = model_comparison,
                              .f = ~ rownames(.) %>% 
                                .[1]),
         best_model = factor(best_model,
                             levels = c("..1", "..2", "..3"),
                             labels = c("baseline",
                                        "cp",
                                        "crm")))

save(list = c("data.exp3.object_model_fits"),
     file = paste0(data_dir, 'experiment3_object/model-fits.RData'))

load(file = paste0(data_dir, 'experiment3_object/model-fits.RData'))

data.exp3.object_model_fits %>% 
  count(best_model) %>%
  arrange(desc(n))

## # A tibble: 3 × 2
##   best_model     n
##   <fct>      <int>
## 1 crm           37
## 2 baseline      11
## 3 cp             2

# count number of positive and negative posteriors on replaceability predictor
data.exp3.object_model_fits %>%
  filter(best_model == 'crm') %>%
  .$fit_crm %>%
  sapply(function(x) fixef(x)[['prob', 'Estimate']]) %>%
  sign() %>%
  as_tibble() %>%
  count(value)

## # A tibble: 2 × 2
##   value     n
##   <dbl> <int>
## 1    -1    29
## 2     1     8

# count number of positive and negative posteriors on contributor predictor
data.exp3.object_model_fits %>%
  filter(best_model == 'crm') %>%
  .$fit_crm %>%
  sapply(function(x) fixef(x)[['contributor', 'Estimate']]) %>%
  sign() %>%
  as_tibble() %>%
  count(value)

## # A tibble: 2 × 2
##   value     n
##   <dbl> <int>
## 1    -1     5
## 2     1    32

5.3 Plots

5.3.1 Parameter search

plot_rmse.exp3.agent + plot_spacer() + plot_rmse.exp3.object +
  plot_annotation(tag_levels = 'A') +
  plot_layout(widths = c(4, 1, 4)) &
  theme(plot.tag = element_text(size = 24))

# ggsave(paste0(figures_dir, 'experiment3/p_search.pdf'), width = 9, height = 4)

5.3.2 Setting up data

data.exp3.agent = data.exp3.agent %>%
  cbind(fit.exp3.agent_baseline %>%
          fitted(newdata = data.exp3.agent,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_baseline = Estimate) %>%
          select(pred_baseline)) %>%
  cbind(fit.exp3.agent_cp %>%
          fitted(newdata = data.exp3.agent,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_cp = Estimate) %>%
          select(pred_cp)) %>%
  cbind(fit.exp3.agent_crm %>%
          fitted(newdata = data.exp3.agent,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_crm = Estimate) %>%
          select(pred_crm))

data.exp3.agent_means = data.exp3.agent %>%
  group_by(trial, craft) %>%
  summarise(resp_mean = mean(resp),
            resp_low = smean.cl.boot(resp)[2],
            resp_high = smean.cl.boot(resp)[3]) %>%
  merge(info.exp3.agent, by = c('trial', 'craft')) %>%
  left_join(data.exp3.agent %>%
              distinct(trial, craft, pred_baseline,
                       pred_cp, pred_crm),
            by = c('trial', 'craft'))

data.exp3.object = data.exp3.object  %>%
  cbind(fit.exp3.object_baseline %>%
          fitted(newdata = data.exp3.object ,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_baseline = Estimate) %>%
          select(pred_baseline)) %>%
  cbind(fit.exp3.object_cp %>%
          fitted(newdata = data.exp3.object ,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_cp = Estimate) %>%
          select(pred_cp)) %>%
  cbind(fit.exp3.object_crm %>%
          fitted(newdata = data.exp3.object ,
                 re_formula = NA) %>%
          as_tibble() %>%
          mutate(pred_crm = Estimate) %>%
          select(pred_crm))

data.exp3.object_means = data.exp3.object  %>%
  group_by(trial, color) %>%
  summarise(resp_mean = mean(resp),
            resp_low = smean.cl.boot(resp)[2],
            resp_high = smean.cl.boot(resp)[3]) %>%
  merge(info.exp3.object, by = c('trial', 'color')) %>%
  left_join(data.exp3.object  %>%
              distinct(trial, color, pred_baseline,
                       pred_cp, pred_crm),
            by = c('trial', 'color'))

labels.exp3 = info.exp3.agent %>% 
  distinct(low, high, availability) %>% 
  arrange(low + high, high, desc(availability)) %>%
  mutate(condition = NA,
         index = 1:n(),
         index_prob = rep(1:(n()/2), each = 2),
         grob = pmap(.l = list(availability, low, high),
                     .f = ~ readPNG(paste0(figures_dir, 'experiment3/sets/',
                                          ..1, ..2, ..3, '.png'))))

model_labels.exp3 = c('baseline' = 'Contribution model',
                      'cp' = 'CP model',
                      'crm' = 'CRM-based model')

model_order.exp3 = c('baseline', 'cp', 'crm')

5.3.3 Scatter plots

scatter.exp3.agent_cp = make_scatter(
  data.exp3.agent_means, 'cp', 'CP model', seq(20, 100, 20))
scatter.exp3.agent_crm = make_scatter(
  data.exp3.agent_means, 'crm', 'CRM-based model', seq(20, 100, 20))

scatter.exp3.object_cp = make_scatter(
  data.exp3.object_means, 'cp', 'CP model', seq(20, 100, 20))
scatter.exp3.object_crm = make_scatter(
  data.exp3.object_means, 'crm', 'CRM-based model', seq(20, 100, 20))

scatter.exp3.agent_cp + inset_element(agent_icon, 0.85, 0.01, 1, 0.25) +
  scatter.exp3.agent_crm + inset_element(agent_icon, 0.85, 0.01, 1, 0.25) +
  scatter.exp3.object_cp + inset_element(object_icon, 0.85, 0.01, 1, 0.15) +
  scatter.exp3.object_crm + inset_element(object_icon, 0.85, 0.01, 1, 0.15) +
  plot_layout(ncol = 2) +
  plot_annotation(tag_levels = list(c('A', '', '', '',
                                      'B', '', '', ''))) &
  theme(plot.title = element_text(hjust = 0.5),
        plot.tag = element_text(size = 24))

# ggsave(paste0(figures_dir, 'experiment3/scatter.pdf'), width = 8, height = 8)

5.3.4 Bar plots

make_bar.exp3 = function(data, fill) {
  g = ggplot(data = data,
             mapping = aes(x = index,
                           y = resp)) +
    stat_summary(fun = "mean",
                 geom = "bar",
                 color = 'black',
                 fill = fill,
                 width = 0.8) +
    stat_summary(fun.data = "mean_cl_boot",
                 geom = "linerange",
                 size = 1) +
    stat_summary(mapping = aes(y = pred,
                               fill = model,
                               shape = model,
                               size = model),
                 fun = 'mean',
                 geom = 'point',
                 color = 'black',
                 position = position_dodge(width = 0.5)) +
    facet_wrap(~ index, scales = 'free_x', nrow = 2) +
    geom_rect(aes(xmin = index - 0.5, xmax = index + 0.5,
                  ymin = -1, ymax = 19.9),
              fill = 'white') +
    geom_custom(data = labels.exp3,
                mapping = aes(data = grob,
                              x = -Inf, y = 9),
                grob_fun = function(x) rasterGrob(x,
                                                  interpolate = T,
                                                  vjust = 0,
                                                  hjust = 0)) +
    geom_text(aes(label = index,
                  y = 5),
              show.legend = F,
              check_overlap = T) +
    theme(text = element_text(size = 16),
          legend.position = c(0.5, 0.94),
          legend.direction = 'horizontal',
          legend.title = element_blank(),
          axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.line.x = element_blank(),
          axis.ticks.x = element_blank(),
          strip.background = element_blank(),
          strip.text = element_blank(),
          panel.grid.major.y = element_line(),
          panel.spacing.y = unit(1.3, "cm"),
          panel.spacing.x = unit(0.1, "cm"),
          plot.margin = margin(t = 0.2, l = 0.2, r = 0.2, b = 1.3, unit = "cm")) +
    scale_fill_manual(values = c('baseline' = 'black',
                                 'cp' = 'gray50',
                                 'crm' = 'white'),
                      labels = model_labels.exp3,
                      name = 'model') +
    scale_size_manual(values = c('baseline' = 2,
                                 'cp' = 2,
                                 'crm' = 2.5),
                      labels = model_labels.exp3,
                      name = 'model') +
    scale_shape_manual(values = c('baseline' = 21,  # circle
                                  'cp' = 24, # triangle
                                  'crm' = 23), # diamond
                       labels = model_labels.exp3,
                       name = 'model') +
    scale_y_continuous(name = 'responsibility judgment',
                       breaks = seq(20, 100, 20)) +
    coord_cartesian(clip = 'off', ylim = c(20, 100))
  
  return (g)
}

barplot.exp3.agent = make_bar.exp3(
  data = data.exp3.agent %>%
    left_join(labels.exp3 %>%
                select(low, high, availability, index, index_prob),
              by = c('low', 'high', 'availability')) %>%
    pivot_longer(cols = starts_with('pred'),
                 names_to = 'model',
                 values_to = 'pred') %>%
    mutate(model = factor(sub('^pred_', '', model),
                          levels = model_order.exp3)),
  fill = 'gray90'
)

barplot.exp3.object = make_bar.exp3(
  data = data.exp3.object %>%
    left_join(labels.exp3 %>%
                select(low, high, availability, index, index_prob),
              by = c('low', 'high', 'availability')) %>%
    pivot_longer(cols = starts_with('pred'),
                 names_to = 'model',
                 values_to = 'pred') %>%
    mutate(model = factor(sub('^pred_', '', model),
                          levels = model_order.exp3)),
  fill = 'gray75'
)

barplot.exp3.agent + inset_element(agent_icon, 0.9, 0.85, 1, 1, align_to = 'full') +
  barplot.exp3.object + inset_element(object_icon, 0.9, 0.9, 1, 1) +
  plot_layout(ncol = 1) +
  plot_annotation(tag_levels = list(c('A', '', 'B', ''))) &
  theme(plot.tag = element_text(size = 24))

# ggsave(paste0(figures_dir, 'experiment3/bar.pdf'), width = 10, height = 12)

6 Session info

sessionInfo()

## R version 4.1.3 (2022-03-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 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.8        
##  [4] purrr_0.3.4         readr_2.1.2         tidyr_1.2.0        
##  [7] tibble_3.1.6        tidyverse_1.3.2     patchwork_1.1.1    
## [10] egg_0.4.5           gridExtra_2.3       png_0.1-7          
## [13] broom.mixed_0.2.9.3 showtext_0.9-5      showtextdb_3.0     
## [16] sysfonts_0.8.8      ggplotify_0.1.0     lemon_0.4.5        
## [19] scales_1.2.0        xtable_1.8-4        Hmisc_4.6-0        
## [22] ggplot2_3.3.5       Formula_1.2-4       survival_3.2-13    
## [25] lattice_0.20-45     Metrics_0.1.4       extraDistr_1.9.1   
## [28] brms_2.17.0         Rcpp_1.0.8.3        DT_0.22            
## [31] lme4_1.1-29         Matrix_1.5-1        ggstatsplot_0.9.1  
## [34] knitr_1.38         
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2             tidyselect_1.1.2       htmlwidgets_1.5.4     
##   [4] munsell_0.5.0          codetools_0.2-18       future_1.24.0         
##   [7] miniUI_0.1.1.1         withr_2.5.0            Brobdingnag_1.2-9     
##  [10] colorspace_2.0-3       highr_0.9              rstudioapi_0.13       
##  [13] stats4_4.1.3           bayesplot_1.9.0        listenv_0.8.0         
##  [16] labeling_0.4.2         emmeans_1.8.2          rstan_2.21.3          
##  [19] bit64_4.0.5            farver_2.1.0           datawizard_0.4.0      
##  [22] bridgesampling_1.1-2   coda_0.19-4            parallelly_1.31.0     
##  [25] vctrs_0.4.0            generics_0.1.2         TH.data_1.1-0         
##  [28] xfun_0.30              R6_2.5.1               markdown_1.1          
##  [31] reshape_0.8.8          gridGraphics_0.5-1     assertthat_0.2.1      
##  [34] vroom_1.5.7            promises_1.2.0.1       multcomp_1.4-18       
##  [37] nnet_7.3-17            googlesheets4_1.0.0    gtable_0.3.0          
##  [40] globals_0.14.0         processx_3.5.3         sandwich_3.0-1        
##  [43] rlang_1.1.1            zeallot_0.1.0          splines_4.1.3         
##  [46] gargle_1.2.0           broom_0.8.0            checkmate_2.0.0       
##  [49] inline_0.3.19          modelr_0.1.8           yaml_2.3.5            
##  [52] reshape2_1.4.4         abind_1.4-5            threejs_0.3.3         
##  [55] crosstalk_1.2.0        backports_1.4.1        httpuv_1.6.5          
##  [58] tensorA_0.36.2         tools_4.1.3            bookdown_0.26         
##  [61] ellipsis_0.3.2         jquerylib_0.1.4        posterior_1.2.1       
##  [64] RColorBrewer_1.1-3     WRS2_1.1-3             ggridges_0.5.3        
##  [67] plyr_1.8.7             base64enc_0.1-3        ps_1.6.0              
##  [70] prettyunits_1.1.1      rpart_4.1.16           correlation_0.8.0     
##  [73] zoo_1.8-9              haven_2.4.3            cluster_2.1.2         
##  [76] fs_1.5.2               furrr_0.2.3            magrittr_2.0.3        
##  [79] data.table_1.14.2      reprex_2.0.1           colourpicker_1.1.1    
##  [82] googledrive_2.0.0      mvtnorm_1.1-3          matrixStats_0.61.0    
##  [85] hms_1.1.1              shinyjs_2.1.0          mime_0.12             
##  [88] evaluate_0.15          shinystan_2.6.0        jpeg_0.1-9            
##  [91] readxl_1.4.0           rstantools_2.2.0       compiler_4.1.3        
##  [94] crayon_1.5.1           minqa_1.2.4            StanHeaders_2.21.0-7  
##  [97] htmltools_0.5.2        tzdb_0.3.0             mc2d_0.1-21           
## [100] later_1.3.0            RcppParallel_5.1.5     lubridate_1.8.0       
## [103] DBI_1.1.2              dbplyr_2.1.1           MASS_7.3-55           
## [106] boot_1.3-28            cli_3.2.0              parallel_4.1.3        
## [109] insight_0.17.0         igraph_1.3.0           pkgconfig_2.0.3       
## [112] statsExpressions_1.3.1 foreign_0.8-82         xml2_1.3.3            
## [115] paletteer_1.4.0        dygraphs_1.1.1.6       bslib_0.3.1           
## [118] estimability_1.4.1     rvest_1.0.2            yulab.utils_0.0.4     
## [121] distributional_0.3.0   callr_3.7.0            digest_0.6.29         
## [124] parameters_0.17.0      rmarkdown_2.13         cellranger_1.1.0      
## [127] htmlTable_2.4.0        shiny_1.7.1            gtools_3.9.2          
## [130] nloptr_2.0.0           lifecycle_1.0.1        nlme_3.1-155          
## [133] jsonlite_1.8.0         fansi_1.0.3            pillar_1.7.0          
## [136] loo_2.5.1              httr_1.4.2             fastmap_1.1.0         
## [139] pkgbuild_1.3.1         glue_1.6.2             xts_0.12.1            
## [142] bayestestR_0.11.5      shinythemes_1.2.0      bit_4.0.4             
## [145] stringi_1.7.6          sass_0.4.1             performance_0.9.0     
## [148] rematch2_2.1.2         latticeExtra_0.6-29

If not me, then who? Responsibility and replacement

Sarah A. Wu and Tobias Gerstenberg

August 24, 2023

1 Setup

2 Helper functions

2.1 Demographics

2.2 Parameter search

2.2.1 Uniform probability

2.2.2 Varying probabilities

2.3 Make scatter plot

3 Experiment 1

3.1 Data

3.1.1 Agent condition

3.1.1.1 Trial info

3.1.1.2 Response data

3.1.1.3 Participant data

3.1.2 Object condition

3.1.2.1 Trial info

3.1.2.2 Response data

3.1.2.3 Participant data

3.1.3 Demographics

3.2 Models

3.2.1 Agent condition

3.2.1.1 Parameter search

3.2.1.2 Overall models

3.2.1.3 Individual models

3.2.2 Object condition

3.2.2.1 Parameter search

3.2.2.2 Overall models

3.2.2.3 Individual models

3.3 Plots

3.3.1 Parameter search

3.3.2 Setting up data

3.3.3 Scatter plots

3.3.4 Bar plot

4 Experiment 2

4.1 Data

4.1.1 Agent condition

4.1.1.1 Trial info

4.1.1.2 Response data

4.1.1.3 Participant data

4.1.2 Object condition

4.1.2.1 Trial info

4.1.2.2 Response data

4.1.2.3 Participant data

4.1.3 Demographics

4.2 Models

4.2.1 Agent condition

4.2.1.1 Parameter search

4.2.1.2 Overall models

4.2.1.3 Individual models

4.2.2 Object condition

4.2.2.1 Parameter search

4.2.2.2 Overall models

4.2.2.3 Individual models

4.3 Plots

4.3.1 Parameter search

4.3.2 Setting up data

4.3.3 Scatter plots

4.3.4 Bar plots

5 Experiment 3

5.1 Data

5.1.1 Agent condition

5.1.1.1 Trial info

5.1.1.2 Response data

5.1.1.3 Participant data

5.1.2 Object condition

5.1.2.1 Trial info

5.1.2.2 Response data

5.1.2.3 Participant data

5.1.3 Demographics

5.2 Models

5.2.1 Agent condition

5.2.1.1 Parameter search

5.2.1.2 Overall models

5.2.1.3 Individual models

5.2.2 Object condition

5.2.2.1 Parameter search

5.2.2.2 Overall models

5.2.2.3 Individual models