PREPARATION

  • Load packages
  • Read in raw data
  • Structure data
  • Calculate demographics

Structure data

goalie_spinner_situations = c("wrong_20","wrong_80",
                      "correct_80","correct_20",
                      "wrong_40","wrong_60",
                      "correct_40","correct_60",
                      "wrong_20","wrong_40",
                      "correct_40","correct_20",
                      "wrong_60","wrong_80",
                      "correct_80","correct_60")

situation_labels = paste0(rep(c("wrong_","correct_"),each=4), seq(20,80,20))

situation_labels2 = paste0(
  rep(c("negative_","positive_"),each=4),
  rep(rep(c("suboptimal_","optimal_"),each=2),2),
  rep(c("nonpivotal","pivotal"),2))

agent_labels = c("bad_agent","average_agent","good_agent")

#GOALIE EXPERIMENT

df.goalie_original = df.data %>%
  filter(experiment == "goalie_original") %>%
  mutate(data_tmp = data %>% strsplit('\n') %>% unlist %>% strsplit(','),
         extra_tmp = extra %>% strsplit('\n') %>% unlist %>% strsplit(','))

df.goalie_original.wide = matrix(unlist(df.goalie_original$data_tmp),nrow=nrow(df.goalie_original),byrow=T) %>%
  as.data.frame(stringsAsFactors=F) %>%
  mutate_all(funs(as.numeric)) %>%
  setNames(paste0(c("round_","ratingLeft_","ratingRight_","time_"),rep(1:8,each=4)))

df.goalie_original.long = df.goalie_original.wide %>%
  mutate(id = 1:nrow(.)) %>%
  wideToLong %>%
  select(-within,-time) %>%
  mutate(round = round + 1) %>%
  gather(key = situation, value = rating,c(ratingLeft,ratingRight)) %>%
  arrange(id,round,situation) %>%
  mutate(situation = rep(goalie_spinner_situations,length(unique(id))),
         situation = factor(situation,levels = situation_labels),
         outcome = str_extract_all(situation,c("wrong|correct"),simplify=T) %>% as.vector,
         probability = str_extract_all(situation,c("20|40|60|80"),simplify=T) %>% as.vector) %>%
  select(id,round,situation,outcome,probability,rating)


#SPINNER EXPERIMENT

df.spinner_original = df.data %>%
  filter(experiment == "spinner_original") %>%
  mutate(data_tmp = data %>% strsplit('\n') %>% unlist %>% strsplit(','),
         extra_tmp = extra %>% strsplit('\n') %>% unlist %>% strsplit(','))

df.spinner_original.wide = matrix(unlist(df.spinner_original$data_tmp),nrow=nrow(df.spinner_original),byrow=T) %>%
  as.data.frame(stringsAsFactors=F) %>%
  mutate_all(funs(as.numeric)) %>%
  setNames(paste0(c("round_","ratingLeft_","ratingRight_","time_"),rep(1:8,each=4)))

df.spinner_original.long = df.spinner_original.wide %>%
  mutate(id = 1:nrow(.)) %>%
  wideToLong %>%
  select(-within,-time) %>%
  mutate(round = round + 1) %>%
  gather(key = situation, value = rating,c(ratingLeft,ratingRight)) %>%
  arrange(id,round,situation) %>%
  mutate(situation = rep(goalie_spinner_situations,length(unique(id))),
         situation = factor(situation,levels = paste0(rep(c("wrong_","correct_"),each=4), seq(20,80,20))),
         outcome = str_extract_all(situation,c("wrong|correct"),simplify=T) %>% as.vector,
         probability = str_extract_all(situation,c("20|40|60|80"),simplify=T) %>% as.vector) %>%
  select(id,round,situation,outcome,probability,rating)

# GOALIE AGENTS EXPERIMENT

df.goalie_skill = df.data %>%
  filter(experiment == "goalie_skill") %>%
  mutate(data = str_replace_all(data,"penalty_miss","wrong"),
         data = str_replace_all(data,"penalty_save","correct"),
         data = str_replace_all(data,"bad goalkeeper","bad_agent"),
         data = str_replace_all(data,"average goalkeeper","average_agent"),
         data = str_replace_all(data,"skilled goalkeeper","good_agent")) %>%
  mutate(data_tmp = data %>% strsplit('\n') %>% unlist %>% strsplit(','),
         extra_tmp = extra %>% strsplit('\n') %>% unlist %>% strsplit(','))

df.goalie_skill.wide = matrix(unlist(df.goalie_skill$data_tmp),nrow=nrow(df.goalie_skill),byrow=T) %>%
  as.data.frame(stringsAsFactors=F) %>%
  setNames(c(paste(c("judgmentRound_","situationLeft_","ratingLeft_","situationRight_","ratingRight_","time_"),rep(1:8,each=6),sep=""),
             paste(c("inferenceRound_","agentTop_","ratingTop_","agentMiddle_","ratingMiddle_","agentBottom_","ratingBottom_"),rep(1:8,each=7),sep="")))

df.goalie_skill.long = df.goalie_skill.wide %>%
  mutate(id = 1:nrow(.)) %>%
  wideToLong() %>%
  select(-within) %>%
  mutate_at(vars(contains("rating"),judgmentRound),funs(as.numeric)) %>%
  mutate(judgmentRound = judgmentRound + 1) %>%
  arrange(id,judgmentRound)

# responsibility judgments
df.goalie_skill.long.responsibility = df.goalie_skill.long %>%
  select(id,judgmentRound,situationLeft,ratingLeft,situationRight,ratingRight) %>%
  gather(key = situation, value = rating,c(ratingLeft,ratingRight)) %>%
  select(-situationLeft,-situationRight) %>%
  rename(round = judgmentRound) %>%
  arrange(id,round,situation) %>%
  mutate(situation = rep(goalie_spinner_situations,length(unique(id))),
         situation = factor(situation,levels = situation_labels),
         outcome = str_extract_all(situation,c("wrong|correct"),simplify=T) %>% as.vector,
         probability = str_extract_all(situation,c("20|40|60|80"),simplify=T) %>% as.vector) %>%
  select(id,round,situation,outcome,probability,rating)

# agent inferences
df.goalie_skill.long.inference = df.goalie_skill.long %>%
  select(id,inferenceRound,contains("Top"),contains("Middle"),contains("Bottom")) %>%
  rename(agent_1 = agentTop, agent_2 = agentMiddle, agent_3 = agentBottom,
         rating_1 = ratingTop, rating_2 = ratingMiddle, rating_3 = ratingBottom,
         situation = inferenceRound) %>%
  wideToLong() %>%
  select(-within) %>%
  mutate(situation = factor(situation,levels = situation_labels),
         agent = factor(agent,levels = agent_labels),
         outcome = str_extract_all(situation,c("wrong|correct"),simplify=T) %>% as.vector,
         probability = str_extract_all(situation,c("20|40|60|80"),simplify=T) %>% as.vector) %>%
  arrange(id,situation,agent) %>%
  select(id,situation,outcome,probability,agent,rating)

# SPINNER AGENTS EXPERIMENT

df.spinner_skill = df.data %>%
  filter(experiment == "spinner_skill") %>%
  mutate(data = str_replace_all(data,"prediction_wrong","wrong"),
         data = str_replace_all(data,"prediction_correct","correct"),
         data = str_replace_all(data,"bad player","bad_agent"),
         data = str_replace_all(data,"average player","average_agent"),
         data = str_replace_all(data,"skilled player","good_agent")) %>%
  mutate(data_tmp = data %>% strsplit('\n') %>% unlist %>% strsplit(','),
         extra_tmp = extra %>% strsplit('\n') %>% unlist %>% strsplit(','))

df.spinner_skill.wide = matrix(unlist(df.spinner_skill$data_tmp),nrow=nrow(df.spinner_skill),byrow=T) %>%
  as.data.frame(stringsAsFactors=F) %>%
  setNames(c(paste(c("judgmentRound_","situationLeft_","ratingLeft_","situationRight_","ratingRight_","time_"),rep(1:8,each=6),sep=""),
             paste(c("inferenceRound_","agentTop_","ratingTop_","agentMiddle_","ratingMiddle_","agentBottom_","ratingBottom_"),rep(1:8,each=7),sep="")))

df.spinner_skill.long = df.spinner_skill.wide %>%
  mutate(id = 1:nrow(.)) %>%
  wideToLong() %>%
  select(-within) %>%
  mutate_at(vars(contains("rating"),judgmentRound),funs(as.numeric)) %>%
  mutate(judgmentRound = judgmentRound + 1) %>%
  arrange(id,judgmentRound)

# responsibility judgments
df.spinner_skill.long.responsibility = df.spinner_skill.long %>%
  select(id,judgmentRound,situationLeft,ratingLeft,situationRight,ratingRight) %>%
  gather(key = situation, value = rating,c(ratingLeft,ratingRight)) %>%
  select(-situationLeft,-situationRight) %>%
  rename(round = judgmentRound) %>%
  arrange(id,round,situation) %>%
  mutate(situation = rep(goalie_spinner_situations,length(unique(id))),
         situation = factor(situation,levels = situation_labels),
         outcome = str_extract_all(situation,c("wrong|correct"),simplify=T) %>% as.vector,
         probability = str_extract_all(situation,c("20|40|60|80"),simplify=T) %>% as.vector) %>%
  select(id,round,situation,outcome,probability,rating)

# agent inferences
df.spinner_skill.long.inference = df.spinner_skill.long %>%
  select(id,inferenceRound,contains("Top"),contains("Middle"),contains("Bottom")) %>%
  rename(agent_1 = agentTop, agent_2 = agentMiddle, agent_3 = agentBottom,
         rating_1 = ratingTop, rating_2 = ratingMiddle, rating_3 = ratingBottom,
         situation = inferenceRound) %>%
  wideToLong() %>%
  select(-within) %>%
  mutate(situation = factor(situation,levels = situation_labels),
         agent = factor(agent,levels = agent_labels),
         outcome = str_extract_all(situation,c("wrong|correct"),simplify=T) %>% as.vector,
         probability = str_extract_all(situation,c("20|40|60|80"),simplify=T) %>% as.vector) %>%
  arrange(id,situation,agent) %>%
  select(id,situation,outcome,probability,agent,rating)

# GARDENER AGENTS EXPERIMENT

df.gardener_skill = df.data %>%
  filter(experiment == "gardener_skill") %>%
  mutate(data = str_replace_all(data,"bad gardener","bad_agent"),
         data = str_replace_all(data,"average gardener","average_agent"),
         data = str_replace_all(data,"skilled gardener","good_agent")) %>%
  mutate(data_tmp = data %>% strsplit('\n') %>% unlist %>% strsplit(','),
         extra_tmp = extra %>% strsplit('\n') %>% unlist %>% strsplit(','))

df.gardener_skill.wide = matrix(unlist(df.gardener_skill$data_tmp),nrow=nrow(df.gardener_skill),byrow=T) %>%
  as.data.frame(stringsAsFactors=F) %>%
  select(-1) %>%
  setNames(c(
    paste(c("situation_","rating_"),rep(1:8,each=2),sep=""),
    paste(c("inferenceRound_","agentTop_","ratingTop_","agentMiddle_","ratingMiddle_","agentBottom_","ratingBottom_"),rep(1:8,each=7),sep="")))

df.gardener_skill.long = df.gardener_skill.wide %>%
  mutate(id = 1:nrow(.)) %>%
  wideToLong() %>%
  # select(-within) %>%
  mutate_at(vars(contains("rating")),funs(as.numeric)) %>%
  arrange(id)

# responsibility judgments
df.gardener_skill.long.responsibility = df.gardener_skill.long %>%
  select(id,situation,rating) %>%
  mutate(situation = factor(situation,levels = situation_labels2),
         outcome = str_extract_all(situation,c("negative|positive"),simplify=T) %>% as.vector,
         probability = str_extract_all(situation,c("suboptimal|optimal"),simplify=T) %>% as.vector,
         pivotality = str_extract_all(situation,c("nonpivotal|pivotal"),simplify=T) %>% as.vector) %>%
  select(id,situation,outcome,probability,pivotality,rating)

# agent inferences
df.gardener_skill.long.inference = df.gardener_skill.long %>%
  select(id,inferenceRound,contains("Top"),contains("Middle"),contains("Bottom")) %>%
  rename(situation = inferenceRound) %>%
  rename(agent_1 = agentTop, agent_2 = agentMiddle, agent_3 = agentBottom,
         rating_1 = ratingTop, rating_2 = ratingMiddle, rating_3 = ratingBottom) %>%
  wideToLong() %>%
  select(-within) %>%
  mutate(situation = factor(situation,levels = situation_labels2),
         agent = factor(agent,levels = agent_labels),
         outcome = str_extract_all(situation,c("negative|positive"),simplify=T) %>% as.vector,
         probability = str_extract_all(situation,c("suboptimal|optimal"),simplify=T) %>% as.vector,
         pivotality = str_extract_all(situation,c("nonpivotal|pivotal"),simplify=T) %>% as.vector) %>%
  arrange(id,situation,outcome,probability,pivotality,agent,rating) %>%
  select(id,situation,outcome,probability,pivotality,agent,rating)

# save data frames in lists

responsibility.list = list(goalie_original = df.goalie_original.long,
                       spinner_original = df.spinner_original.long,
                       goalie = df.goalie_skill.long.responsibility,
                       spinner = df.spinner_skill.long.responsibility,
                       # manager = df.manager_skill.long.responsibility,
                       gardener = df.gardener_skill.long.responsibility)
                       # gardener_noskill = df.gardener_noskill.long)

inference.list = list(goalie = df.goalie_skill.long.inference,
                   spinner = df.spinner_skill.long.inference,
                   # manager = df.manager_skill.long.inference,
                   gardener = df.gardener_skill.long.inference)

# remove unnecessary variables
rm(list = setdiff(ls(),c("df.data","responsibility.list","inference.list")))

Demographics

demographics.list = list()

# Experiment 1

var.names = c("difficulty","gender","age","time",paste0("decision_",1:10))

demographics.list[["goalie_spinner_original"]] = df.data %>%
  filter(experiment %in% c("goalie_original","spinner_original")) %>%
  select(extra) %>%
  unlist %>%
  strsplit("\n") %>%
  unlist %>%
  strsplit(",") %>%
  unlist %>%
  matrix(ncol = length(var.names),byrow=T) %>%
  as.data.frame(stringsAsFactors = F) %>%
  setNames(var.names)

# Experiment 2

demographics.list[["goalie_spinner"]] = df.data %>%
  filter(experiment %in% c("goalie_skill","spinner_skill")) %>%
  select(extra) %>%
  unlist %>%
  strsplit("\n") %>%
  unlist %>%
  strsplit(",") %>%
  unlist %>%
  matrix(ncol = length(var.names),byrow=T) %>%
  as.data.frame(stringsAsFactors = F) %>%
  setNames(var.names)

var.names = c("index","difficulty","gender","age","time")

# Experiment 3
demographics.list[["gardener"]] = df.data %>%
  filter(experiment %in% c("gardener_skill")) %>%
  select(extra) %>%
  unlist %>%
  strsplit("\n") %>%
  unlist %>%
  strsplit(",") %>%
  unlist %>%
  matrix(ncol = length(var.names),byrow=T) %>%
  as.data.frame(stringsAsFactors = F) %>%
  setNames(var.names)

# descriptive information summary

df.descriptive_summary = as.data.frame(matrix(NA,ncol=6,nrow=0)) %>%
  setNames(c("experiment","age.mean","age.sd","n.female","time.mean","time.sd"))

for(i in 1:length(demographics.list)){
  df.tmp = demographics.list[[names(demographics.list)[i]]] %>%
    mutate_all(funs(as.numeric)) %>%
    summarise(experiment = names(demographics.list)[i],
              age.mean = mean(age),
              age.sd = sd(age),
              n.female = sum(gender == 1),
              time.mean = mean(time),
              time.sd = sd(time)
              ) %>%
    mutate_at(vars(-experiment),funs(round(.,2)))
  df.descriptive_summary = rbind(df.descriptive_summary,df.tmp)
}

kable(df.descriptive_summary, caption = "Demographics")
Demographics
experiment age.mean age.sd n.female time.mean time.sd
goalie_spinner_original 33.70 11.14 39 6.19 1.75
goalie_spinner 34.80 11.26 27 11.44 4.46
gardener 33.54 10.61 17 10.38 4.86

MODEL

  • Setting up the model
  • Defining the agent decision functions
  • Fitting beta in softmax and priors over the agent types
  • Calculate difference between prior expected reward and posterior expected reward
  • Run different regression models

Fit beta and priors to posteriors

df.posteriors = rbind(
  inference.list[["goalie"]] %>% group_by(probability,outcome,agent) %>%summarise(rating = mean(rating)) %>% mutate(experiment = "goalie", pivotality = "pivotal"),
  inference.list[["spinner"]] %>% group_by(probability,outcome,agent) %>%summarise(rating = mean(rating)) %>% mutate(experiment = "spinner", pivotality = "pivotal"),
  inference.list[["gardener"]] %>% group_by(probability,pivotality,outcome,agent) %>%summarise(rating = mean(rating)) %>% mutate(experiment = "gardener")
)

#TODO: change the probabilities around here
df.posteriors = df.posteriors %>%
  spread(agent,rating) %>%
  ungroup %>%
  mutate(probability = ifelse(experiment == "gardener", str_replace_all(probability,"suboptimal","30"),probability),
         probability = ifelse(experiment == "gardener", str_replace_all(probability,"optimal","70"),probability),
         probability = as.numeric(probability),
         outcome = str_replace_all(outcome,"wrong","negative"),
         outcome = str_replace_all(outcome,"correct","positive"),
         outcome = as.numeric(as.character(factor(outcome,levels=c("negative","positive"),labels = c("0","1")))),
         outcome_counterfactual = ifelse(pivotality == "pivotal",1-outcome,outcome),
         probability_counterfactual = 100-probability,
         experiment = factor(experiment,levels = c("goalie", "spinner", "gardener")))%>%
  mutate_at(vars(probability,probability_counterfactual),funs(./100)) %>%
  select(experiment,outcome,probability,everything()) %>%
  arrange(experiment,outcome,probability)

# fit the parameters
fit_prior = function(parameters){
  prior = matrix(c(rep(c(parameters[2],parameters[3],1-(parameters[2]+parameters[3])),8),
                   rep(c(parameters[4],parameters[5],1-(parameters[4]+parameters[5])),8),
                   rep(c(parameters[6],parameters[7],1-(parameters[6]+parameters[7])),8)),
                   ncol=3,byrow=T)

  likelihood = cbind(
    agent.bad(parameters[1],df.posteriors,reward),
    agent.average(parameters[1],df.posteriors,reward),
    agent.good(parameters[1],df.posteriors,reward)
  )
  posterior_fit = prior*likelihood/rowSums(prior*likelihood)
  posterior_data = select(df.posteriors,contains("agent"))/100
  return(sum((posterior_data-posterior_fit)^2))
}

# output the predicted posteriors, based on the fit priors and beta
predicted_posterior = function(parameters){
  prior = matrix(c(rep(c(parameters[2],parameters[3],1-(parameters[2]+parameters[3])),8),
                   rep(c(parameters[4],parameters[5],1-(parameters[4]+parameters[5])),8),
                   rep(c(parameters[6],parameters[7],1-(parameters[6]+parameters[7])),8)),
                   ncol=3,byrow=T)

  likelihood = cbind(
    agent.bad(parameters[1],df.posteriors,reward),
    agent.average(parameters[1],df.posteriors,reward),
    agent.good(parameters[1],df.posteriors,reward)
  )
  posterior_fit = prior*likelihood/rowSums(prior*likelihood)
  return(posterior_fit)
}

# reward
reward = 1

# optimize the function
fit = optim(par = runif(7,0,1),fit_prior,method = "L-BFGS-B",
                  lower = rep(0, 7), upper = c(5,rep(1,6)))

# fitted parameters
beta = fit$par[1]
prior.list = list()
prior.list[["goalie"]] = c(fit$par[2],fit$par[3],1-sum(fit$par[2],fit$par[3]))
prior.list[["spinner"]] = c(fit$par[4],fit$par[5],1-sum(fit$par[4],fit$par[5]))
prior.list[["gardener"]] = c(fit$par[6],fit$par[7],1-sum(fit$par[6],fit$par[7]))

# add predicted posteriors to dataframe
tmp = predicted_posterior(fit$par) %>%
  as.data.frame() %>%
  setNames(c("p.bad_agent","p.average_agent","p.good_agent"))

df.posteriors = df.posteriors %>%
  cbind(.,tmp)

#round beta
beta = beta %>% round(1)

Expected reward prior

situations = list()
exp_r_prior = list()
exp_r_prior_scaled = list()

# SITUATIONS: GOALIE & SPINNER
situations[["goalie_spinner"]] = expand.grid(probability = seq(0.1,0.9,0.05), outcome = c(0,1)) %>%
  mutate(probability_counterfactual = 1-probability,
         outcome_counterfactual = 1-outcome,
         probability.overall = ifelse(outcome == 1, probability, probability_counterfactual),
         probability.overall = probability.overall/sum(probability.overall))

# SITUATIONS: GARDENER
situations[["gardener"]] = expand.grid(probability = seq(0.1,0.9,0.05), outcome = c(0,1),
                                       probability_counterfactual = seq(0.1,0.9,0.05), outcome_counterfactual = c(0,1)) %>%
  mutate(probability.overall = (1-abs(probability-outcome))*(1-abs(probability_counterfactual-outcome_counterfactual)),
         probability.overall = probability.overall/sum(probability.overall))

situation.names = c("goalie_spinner","gardener")

for(i in 1:length(situation.names)){
  #expected reward prior
  situations[[situation.names[i]]]$agent.bad = agent.bad(beta,situations[[situation.names[i]]],reward)
  situations[[situation.names[i]]]$agent.average = agent.average(beta,situations[[situation.names[i]]],reward)
  situations[[situation.names[i]]]$agent.good = agent.good(beta,situations[[situation.names[i]]],reward)

  exp_r_prior[[situation.names[i]]] = situations[[situation.names[i]]] %>%
    summarise(exp_r_prior.bad = sum((agent.bad*outcome + ((1-agent.bad)*outcome_counterfactual))*probability.overall),
              exp_r_prior.average = sum((agent.average*outcome + ((1-agent.average)*outcome_counterfactual))*probability.overall),
              exp_r_prior.good = sum((agent.good*outcome + ((1-agent.good)*outcome_counterfactual))*probability.overall))

  #rescale expected prior rewards
  exp_r_prior_scaled[[situation.names[i]]] = rescale(exp_r_prior[[situation.names[i]]] %>% as.numeric(),to=c(0,1),from = range(exp_r_prior[[situation.names[i]]]))
}

Regression models

df.regression.goalie = responsibility.list[["goalie"]] %>%
  group_by(situation,outcome,probability) %>%
  summarise(mean = mean(rating),
            low = smean.cl.boot(rating)[2],
            high = smean.cl.boot(rating)[3]
  ) %>%
  ungroup() %>%
  mutate(pivotality = 1) %>%
  select(situation,outcome,probability,pivotality,mean,low,high)

df.regression.spinner = responsibility.list[["spinner"]] %>%
  group_by(situation,outcome,probability) %>%
  summarise(mean = mean(rating),
            low = smean.cl.boot(rating)[2],
            high = smean.cl.boot(rating)[3]
  ) %>%
  ungroup() %>%
  mutate(pivotality = 1) %>%
  select(situation,outcome,probability,pivotality,mean,low,high)

df.regression.gardener = responsibility.list[["gardener"]] %>%
  group_by(situation,outcome,probability,pivotality) %>%
  summarise(mean = mean(rating),
            low = smean.cl.boot(rating)[2],
            high = smean.cl.boot(rating)[3]
  ) %>%
  ungroup() %>%
  select(situation,outcome,probability,pivotality,mean,low,high)


df.difference = rbind(exp_r_difference.list[["goalie"]],
                      exp_r_difference.list[["spinner"]],
                      exp_r_difference.list[["gardener"]]
)

df.regression = rbind(df.regression.goalie,
                      df.regression.spinner,
                      df.regression.gardener
) %>%
  mutate(outcome = str_replace_all(outcome,"wrong|negative","0"),
         outcome = str_replace_all(outcome,"correct|positive","1"),
         outcome = as.numeric(outcome),
         pivotality = str_replace_all(pivotality,"nonpivotal","0.5"), #can set different values here for pivotality
         pivotality = str_replace_all(pivotality,"pivotal","1"),
         pivotality = as.numeric(pivotality),
         mean = ifelse(outcome == 0,-mean,mean),
         pivotality = ifelse(outcome == 0,-pivotality,pivotality),
         difference = as.vector(df.difference$exp_r_difference),
         difference.model = as.vector(df.difference$exp_r_difference.model),
         probability = str_replace_all(probability,"suboptimal","30"),
         probability = str_replace_all(probability,"optimal","70"),
         probability = as.numeric(probability)
  ) %>%
  mutate(experiment = rep(c("goalie","spinner","gardener"),each=8)) %>%
  select(experiment,situation,outcome,difference,difference.model,pivotality,mean,low,high,probability)

# list for regression fits
fit.regression = list()

# difference and pivotality
fit.regression[["difference_pivotality"]] = lm(mean~scale(difference)+scale(pivotality),data=df.regression)

# round coefficients
fit.regression[["difference_pivotality"]]$coefficients = fit.regression[["difference_pivotality"]]$coefficients %>% round(0)

df.regression$p.difference_pivotality = fit.regression[["difference_pivotality"]]$coefficients[1] +
  fit.regression[["difference_pivotality"]]$coefficients[2] * scale(df.regression$difference)[,] +
  fit.regression[["difference_pivotality"]]$coefficients[3] * scale(df.regression$pivotality)[,]

# difference only
fit.regression[["difference"]] = lm(mean~difference,data=df.regression)
df.regression$p.difference = fit.regression[["difference"]]$fitted.values

# pivotality only
fit.regression[["pivotality"]] = lm(mean~pivotality,data=df.regression)
df.regression$p.pivotality = fit.regression[["pivotality"]]$fitted.values

# optimality model
df.regression = df.regression %>%
  mutate(optimal = 1,
         optimal = ifelse(str_detect(situation,"20|40|suboptimal"),0,optimal))

fit.regression[["optimality"]] = lm(mean~optimal*outcome,data=df.regression)
df.regression$p.optimality = fit.regression[["optimality"]]$fitted.values

TABLES

  • Table with average agent inferences across all experiments
  • Table with average responsibility judgments across all experiments

Agent inferences longtable

Agent inferences
experiment probability outcome pivotality p.bad_agent p.average_agent p.good_agent bad_agent average_agent good_agent
goalie 20 0 pivotal 53.53 31.63 14.85 61.44 24.05 14.51
goalie 40 0 pivotal 46.87 40.13 13.00 47.07 35.37 17.56
goalie 60 0 pivotal 41.29 47.25 11.45 29.27 53.29 17.44
goalie 80 0 pivotal 37.21 52.46 10.32 33.41 47.49 19.10
goalie 20 1 pivotal 11.69 29.43 58.88 12.32 19.51 68.17
goalie 40 1 pivotal 10.32 37.67 52.00 12.39 32.20 55.41
goalie 60 1 pivotal 9.16 44.68 46.15 11.10 54.02 34.88
goalie 80 1 pivotal 8.30 49.88 41.82 12.56 54.88 32.56
spinner 20 0 pivotal 53.90 33.46 12.63 59.90 29.66 10.44
spinner 40 0 pivotal 46.86 42.16 10.98 36.66 49.49 13.85
spinner 60 0 pivotal 41.04 49.34 9.62 32.39 47.05 20.56
spinner 80 0 pivotal 36.82 54.55 8.63 32.15 40.71 27.15
spinner 20 1 pivotal 12.66 33.48 53.86 23.20 19.54 57.27
spinner 40 1 pivotal 11.00 42.18 46.82 21.27 33.46 45.27
spinner 60 1 pivotal 9.63 49.36 41.00 9.10 60.22 30.68
spinner 80 1 pivotal 8.64 54.57 36.79 8.59 60.17 31.24
gardener 30 0 nonpivotal 33.36 30.05 36.60 54.10 25.12 20.78
gardener 30 0 pivotal 55.14 30.66 14.20 60.98 25.44 13.59
gardener 70 0 nonpivotal 26.99 43.41 29.61 31.56 46.15 22.29
gardener 70 0 pivotal 44.44 44.12 11.44 33.66 49.27 17.07
gardener 30 1 nonpivotal 33.36 30.05 36.60 23.88 26.49 49.63
gardener 30 1 pivotal 12.43 29.46 58.11 12.44 17.80 69.76
gardener 70 1 nonpivotal 26.99 43.41 29.61 14.66 43.85 41.49
gardener 70 1 pivotal 10.10 42.71 47.19 13.90 52.32 33.78

Responsibility judgments longtable

Responsibility judgments
experiment probability outcome difference pivotality p.difference_pivotality rating
goalie 20 0 -0.28 -1.0 -70.35 -76.84
goalie 40 0 -0.18 -1.0 -54.55 -62.29
goalie 60 0 -0.07 -1.0 -37.21 -42.50
goalie 80 0 -0.09 -1.0 -40.31 -32.24
goalie 20 1 0.23 1.0 75.47 76.50
goalie 40 1 0.18 1.0 68.08 73.00
goalie 60 1 0.11 1.0 57.55 64.60
goalie 80 1 0.09 1.0 54.79 60.95
spinner 20 0 -0.27 -1.0 -68.13 -60.85
spinner 40 0 -0.11 -1.0 -43.44 -50.46
spinner 60 0 -0.06 -1.0 -35.41 -34.10
spinner 80 0 -0.03 -1.0 -31.39 -29.43
spinner 20 1 0.14 1.0 61.63 54.62
spinner 40 1 0.10 1.0 56.63 56.27
spinner 60 1 0.13 1.0 60.16 62.98
spinner 80 1 0.13 1.0 60.98 63.49
gardener 30 0 -0.20 -0.5 -40.37 -35.56
gardener 30 0 -0.27 -1.0 -67.97 -73.68
gardener 70 0 -0.05 -0.5 -17.37 -20.27
gardener 70 0 -0.08 -1.0 -39.15 -32.88
gardener 30 1 0.10 0.5 39.25 41.85
gardener 30 1 0.25 1.0 78.69 77.29
gardener 70 1 0.13 0.5 43.66 39.61
gardener 70 1 0.10 1.0 56.76 52.27

STATS

Test for context effects in Experiments 1 and 2

## Anova Table (Type 3 tests)
## 
## Response: rating
##                        Effect           df     MSE         F   ges p.value
## 1                     outcome        1, 40  644.13 42.37 ***   .06  <.0001
## 2                 probability  1.66, 66.52 1964.02 37.37 ***   .21  <.0001
## 3                     context        1, 40  272.03      2.12  .001     .15
## 4         outcome:probability  1.39, 55.73 2388.68   8.39 **   .06    .002
## 5             outcome:context        1, 40  361.19      0.63 .0005     .43
## 6         probability:context  2.14, 85.49  394.62      0.34 .0006     .72
## 7 outcome:probability:context 2.52, 100.81  260.65      1.24  .002     .30
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
## 
## Sphericity correction method: GG

ANOVAs for responsibility judgments (all experiments)

## $F(1,81) = 9.55, p = 0.003, \eta^2 = 0.11$ $F(3,243) = 71.58, p = 0, \eta^2 = 0.47$ $F(3,243) = 0.33, p = 0.806, \eta^2 = 0$
## $F(1,40) = 45.29, p = 0, \eta^2 = 0.53$ $F(1,40) = 32.19, p = 0, \eta^2 = 0.45$ $F(1,40) = 19.91, p = 0, \eta^2 = 0.33$

Model comparisons

## Likelihood ratio test
## 
## Model 1: mean ~ scale(difference) + scale(pivotality)
## Model 2: mean ~ difference
##   #Df   LogLik Df Chisq Pr(>Chisq)    
## 1   4  -72.939                        
## 2   3 -100.104 -1 54.33  1.695e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Variance Table
## 
## Response: mean
##                   Df Sum Sq Mean Sq F value  Pr(>F)    
## scale(difference)  1  67203   67203    2302 < 2e-16 ***
## scale(pivotality)  1   5284    5284     181 8.6e-12 ***
## Residuals         21    613      29                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Model comparison
Model r rmse BIC
difference_pivotality 0.9531723 5.116225 158.5906
difference 0.8379497 15.674916 209.7422
pivotality 0.6592709 12.671361 199.5317
optimality 0.7461116 11.084670 199.4663

Ordinal logistic regression on agent inferences

##                 Wald Statistics          Response: mostlikely 
## 
##  Factor                                                           
##  experiment  (Factor+Higher Order Factors)                        
##   All Interactions                                                
##  probability  (Factor+Higher Order Factors)                       
##   All Interactions                                                
##  outcome  (Factor+Higher Order Factors)                           
##   All Interactions                                                
##  experiment * probability  (Factor+Higher Order Factors)          
##  experiment * outcome  (Factor+Higher Order Factors)              
##  probability * outcome  (Factor+Higher Order Factors)             
##  experiment * probability * outcome  (Factor+Higher Order Factors)
##  TOTAL INTERACTION                                                
##  TOTAL                                                            
##  Chi-Square d.f. P     
##    5.66     4    0.2258
##    5.65     3    0.1297
##   69.44     4    <.0001
##   69.43     3    <.0001
##  199.81     4    <.0001
##   71.76     3    <.0001
##    1.27     2    0.5299
##    4.96     2    0.0839
##   68.75     2    <.0001
##    0.37     1    0.5416
##   72.30     4    <.0001
##  200.92     7    <.0001
##                 Wald Statistics          Response: mostlikely 
## 
##  Factor     Chi-Square d.f. P     
##  outcome    144.78     1    <.0001
##  TOTAL      144.78     1    <.0001
##                 Wald Statistics          Response: mostlikely 
## 
##  Factor      Chi-Square d.f. P     
##  probability 0.41       1    0.5243
##  TOTAL       0.41       1    0.5243
##                 Wald Statistics          Response: mostlikely 
## 
##  Factor                                               Chi-Square d.f.
##  probability  (Factor+Higher Order Factors)            68.18     2   
##   All Interactions                                     68.18     1   
##  outcome  (Factor+Higher Order Factors)               198.42     2   
##   All Interactions                                     68.18     1   
##  probability * outcome  (Factor+Higher Order Factors)  68.18     1   
##  TOTAL                                                198.66     3   
##  P     
##  <.0001
##  <.0001
##  <.0001
##  <.0001
##  <.0001
##  <.0001
## # A tibble: 6 x 3
## # Groups:   experiment [?]
##   experiment agent         rating
##   <chr>      <fct>          <dbl>
## 1 goalie     bad_agent        12.
## 2 goalie     average_agent    26.
## 3 goalie     good_agent       62.
## 4 spinner    bad_agent        22.
## 5 spinner    average_agent    26.
## 6 spinner    good_agent       51.
##                 Wald Statistics          Response: mostlikely 
## 
##  Factor     Chi-Square d.f. P     
##  outcome    89.23      1    <.0001
##  TOTAL      89.23      1    <.0001
##                 Wald Statistics          Response: mostlikely 
## 
##  Factor      Chi-Square d.f. P     
##  probability 0.67       1    0.4137
##  TOTAL       0.67       1    0.4137
##                 Wald Statistics          Response: mostlikely 
## 
##  Factor     Chi-Square d.f. P     
##  pivotality 0.13       1    0.7235
##  TOTAL      0.13       1    0.7235
##                 Wald Statistics          Response: mostlikely 
## 
##  Factor                                               Chi-Square d.f.
##  probability  (Factor+Higher Order Factors)            68.21     2   
##   All Interactions                                     68.10     1   
##  outcome  (Factor+Higher Order Factors)               139.82     2   
##   All Interactions                                     68.10     1   
##  probability * outcome  (Factor+Higher Order Factors)  68.10     1   
##  TOTAL                                                139.89     3   
##  P     
##  <.0001
##  <.0001
##  <.0001
##  <.0001
##  <.0001
##  <.0001
##                 Wald Statistics          Response: mostlikely 
## 
##  Factor                                              Chi-Square d.f.
##  outcome  (Factor+Higher Order Factors)              92.48      2   
##   All Interactions                                    4.29      1   
##  pivotality  (Factor+Higher Order Factors)            4.42      2   
##   All Interactions                                    4.29      1   
##  outcome * pivotality  (Factor+Higher Order Factors)  4.29      1   
##  TOTAL                                               92.58      3   
##  P     
##  <.0001
##  0.0384
##  0.1096
##  0.0384
##  0.0384
##  <.0001
##                 Wald Statistics          Response: mostlikely 
## 
##  Factor                                                           
##  probability  (Factor+Higher Order Factors)                       
##   All Interactions                                                
##  outcome  (Factor+Higher Order Factors)                           
##   All Interactions                                                
##  pivotality  (Factor+Higher Order Factors)                        
##   All Interactions                                                
##  probability * outcome  (Factor+Higher Order Factors)             
##  probability * pivotality  (Factor+Higher Order Factors)          
##  outcome * pivotality  (Factor+Higher Order Factors)              
##  probability * outcome * pivotality  (Factor+Higher Order Factors)
##  TOTAL INTERACTION                                                
##  TOTAL                                                            
##  Chi-Square d.f. P     
##   70.97     4    <.0001
##   70.00     3    <.0001
##  134.86     4    <.0001
##   69.40     3    <.0001
##   12.69     4    0.0129
##   12.68     3    0.0054
##   68.27     2    <.0001
##    9.45     2    0.0089
##    8.42     2    0.0149
##    4.79     1    0.0286
##   71.01     4    <.0001
##  136.60     7    <.0001

PLOTS

Plotting functions

# responsiblity judgments without model predictions
plot.responsibility.no_model = function(x,bottom_space){
  ggplot(x,aes(x=probability,y=rating,fill=outcome))+
    stat_summary(fun.y = mean, geom = "bar",color="black") +
    stat_summary(fun.data = mean_cl_boot, geom = "errorbar",width=0, size = 2)+
    facet_wrap(~outcome)+
    scale_y_continuous(breaks = seq(0,100,25), expand = c(0,0))+
    coord_cartesian(ylim=c(0,100))+
    labs(x = "situation", y = "responsibility rating")+
    scale_fill_manual(values=c("red3","green3"))+
    theme_bw()+
    theme(legend.position = "none",
          text = element_text(size=24),
          strip.text = element_text(size=24),
          axis.text.x = element_text(size = 20),
          strip.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.margin = unit(c(0,0.2,bottom_space,0.2),"cm")
    )
}

# responsibility judgments with model predictions
plot.responsibility.with_model = function(x,bottom_space){
  ggplot(x,aes(x=probability,y=value,group=index,fill=color_index))+
    stat_summary(fun.y = mean, geom = "bar",color="black",position = position_dodge(0.9)) +
    stat_summary(fun.data = mean_cl_boot, geom = "errorbar",position = position_dodge(0.9),width=0, size = 2)+
    facet_wrap(~outcome)+
    scale_y_continuous(breaks = seq(0,100,25), expand = c(0,0))+
    coord_cartesian(ylim=c(0,100))+
    labs(x = "situation", y = "responsibility rating")+
    scale_fill_manual(values=c("red3","green3","black"))+
    theme_bw()+
    theme(legend.position = "none",
          text = element_text(size=24),
          strip.text = element_text(size=24),
          axis.text.x = element_text(size = 20),
          strip.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.margin = unit(c(0,0.2,bottom_space,0.2),"cm")
    )
}

plot.inferences = function(x,experiment.name, title.name){
  if (experiment.name %in% c("spinner","goalie")){
    agentlabels = c("unskilled  ","average  ","skilled  ")
  }else{
    agentlabels = c("bad  ","average  ","skilled  ")
  }

  ggplot(x,aes(x=probability,y=rating,group=rev(agent),fill=rev(agent)))+
    stat_summary(fun.y = mean, geom = "bar",color="black",position=position_fill(),width=0.75) +
    facet_wrap(~outcome,ncol=2)+
    labs(y = "agent inference", fill = "type of agent:  ", x = "situation")+
    scale_fill_manual(values = rev(c("peru", rgb(0.7,0.7,0.7), "gold")),labels = rev(agentlabels))+
    scale_alpha_manual(values = c(1,0.5),guide='none')+
    scale_y_continuous(breaks = seq(0,1,0.25), labels = paste0(seq(0,100,25),"%"))+
    coord_cartesian(xlim = c(0.5, 4.5), ylim = c(0,1.01), expand=F)+
    theme_bw()+
    theme(
      axis.title.x=element_text(margin = margin(t = 15,unit="pt")),
      legend.position="bottom",
      text = element_text(size=24),
      strip.text = element_text(size=24),
      axis.text.x = element_text(size = 20),
      axis.ticks.x = element_blank(),
      legend.key = element_rect(color="black"),
      panel.grid.minor = element_blank(),
      panel.grid.major = element_blank(),
      strip.background = element_blank()
    )+
    guides(fill = guide_legend(reverse = T, title.vjust = 1.5))
}

Responsibility plots with model predictions

model.name = "difference_pivotality"
# model.name = "pivotality"
# model.name = "difference"

experiment.name = "goalie"
# experiment.name = "spinner"
# experiment.name = "gardener"

df.prediction = df.regression %>%
  filter(experiment == experiment.name)

if (experiment.name %in% c("goalie","spinner")){
  df.plot = responsibility.list[[experiment.name]] %>%
    select(id,situation,outcome,probability,rating) %>%
    merge(.,select(df.prediction,situation,p.difference,p.pivotality,p.difference_pivotality)) %>%
    mutate_at(vars(p.difference,p.pivotality,p.difference_pivotality),funs(abs)) %>%
    mutate(probability = paste0(probability,"%"),
           outcome = str_replace_all(outcome,"wrong","negative outcome"),
           outcome = str_replace_all(outcome,"correct","positive outcome")
    )
}else{
  df.plot = responsibility.list[[experiment.name]] %>%
    select(id,situation,outcome,probability,pivotality,rating) %>%
    rename(pivotality.index = pivotality) %>%
    merge(.,select(df.prediction,situation,p.difference,p.pivotality,p.difference_pivotality)) %>%
    mutate_at(vars(p.difference,p.pivotality,p.difference_pivotality),funs(abs)) %>%
    mutate(pivotality.index = str_replace_all(pivotality.index,"nonpivotal","not\npivotal"),
           probability = str_replace_all(probability,"suboptimal","30%\n"),
           probability = str_replace_all(probability,"optimal","70%\n"),
           probability = paste0(probability,pivotality.index),
           outcome = str_replace_all(outcome,"negative","negative outcome"),
           outcome = str_replace_all(outcome,"positive","positive outcome")
    )
}

df.plot = df.plot %>%
  select_("id","outcome","probability",paste0("p.",model.name),"rating") %>%
  gather_("index","value",c("rating",paste0("p.",model.name))) %>%
  mutate(color_index = 1,
         color_index = ifelse(outcome == "negative outcome",0,color_index),
         color_index = ifelse(index == paste0("p.",model.name),2,color_index),
         color_index = as.factor(color_index),
         index = factor(index,levels = c("rating",model.name),labels = c("rating",model.name))) %>%
  arrange(index,id,outcome,probability)

if (experiment.name %in% c("goalie","spinner")){
  bottom_space = 0.2
}else{
  bottom_space = 2
}

# plot results
plot.responsibility.with_model(df.plot,bottom_space)

Scatter plots (responsibility)

model.name = "difference_pivotality"
# model.name = "pivotality"
# model.name = "difference"
# model.name = "optimality"

df.plot = df.regression %>%
  mutate_at(vars(p.difference,p.pivotality,p.difference_pivotality,p.optimality,mean),funs(abs)) %>%
  mutate(experiment = factor(experiment,levels = c("goalie","spinner","gardener","manager")),
         outcome = factor(outcome,levels = c(0,1), labels = c("negative", "positive")))

ggplot(df.plot,aes_string(x=paste0("p.",model.name),y="mean",color="experiment"))+
  geom_abline(intercept = 0,slope = 1, linetype = 2)+
  geom_smooth(aes(group=1), method="lm",color="black")+
  geom_errorbar(aes(ymin=low,ymax=high),alpha = 1,width=0,size=1)+
  geom_point(aes(shape=outcome),size=4)+
  annotate(x = 10, y = 95, geom = 'text', label = paste("R^2 == ", cor(df.plot$mean,df.plot[[paste0("p.",model.name)]])^2 %>% round(2)), size = 8, parse = T)+
  theme_bw()+
  labs(y = "responsibility rating", x="model prediction")+
  scale_x_continuous(breaks = seq(0,100,25))+
  scale_y_continuous(breaks = seq(0,100,25))+
  coord_cartesian(xlim=c(0,100),ylim=c(0,100))+
  theme(text = element_text(size=26),
        legend.position=c(1,0),
        legend.justification=c(1.01,-0.01),
        legend.text = element_text(size=20),
        legend.title = element_text(size=22),
        legend.background = element_blank(),
        legend.key = element_blank(),
        panel.grid = element_blank(),
        panel.border = element_blank(),
        axis.line = element_line(colour = "black", size = 1),
        axis.title.x = element_text(margin = margin(0.2,0,0,0,"cm")),
        axis.title.y = element_text(margin = margin(0,0.1,0,0,"cm"))
  )+
  guides(color = guide_legend(order = 1),
         shapes = guide_legend(order = 2))

CROSS VALIDATION

Data frame

df.cv = rbind(
  inference.list[["goalie"]] %>% group_by(probability,outcome,agent) %>%summarise(rating = mean(rating)) %>% mutate(experiment = "goalie", pivotality = "pivotal"),
  inference.list[["spinner"]] %>% group_by(probability,outcome,agent) %>%summarise(rating = mean(rating)) %>% mutate(experiment = "spinner", pivotality = "pivotal"),
  inference.list[["gardener"]] %>% group_by(probability,pivotality,outcome,agent) %>%summarise(rating = mean(rating)) %>% mutate(experiment = "gardener")
)

df.cv = df.cv %>%
  spread(agent,rating) %>%
  ungroup %>%
  mutate(probability = ifelse(experiment == "gardener", str_replace_all(probability,"suboptimal","30"),probability),
         probability = ifelse(experiment == "gardener", str_replace_all(probability,"optimal","70"),probability),
         outcome = str_replace_all(outcome,"wrong","negative"),
         outcome = str_replace_all(outcome,"correct","positive"),
         outcome = as.numeric(as.character(factor(outcome,levels=c("negative","positive"),labels = c("0","1")))),
         probability = as.numeric(probability),
         pivotality = str_replace_all(pivotality,"nonpivotal","0.5"),
         pivotality = str_replace_all(pivotality,"pivotal","1"),
         pivotality = as.numeric(pivotality),
         pivotality = ifelse(outcome == 0,-pivotality,pivotality),
         outcome_counterfactual = ifelse(pivotality == "pivotal",1-outcome,outcome),
         probability_counterfactual = 100-probability,
         experiment = factor(experiment,levels = c("goalie", "spinner","gardener")))%>%
  mutate_at(vars(probability,probability_counterfactual),funs(./100)) %>%
  select(experiment,outcome,probability,everything()) %>%
  arrange(experiment,outcome,probability)

df.regression.goalie = responsibility.list[["goalie"]] %>%
  group_by(situation,outcome,probability) %>%
  summarise(rating = mean(rating)) %>%
  ungroup() %>%
  mutate(pivotality = 1) %>%
  select(situation,outcome,probability,pivotality,rating)

df.regression.spinner = responsibility.list[["spinner"]] %>%
  group_by(situation,outcome,probability) %>%
  summarise(rating = mean(rating)) %>%
  ungroup() %>%
  mutate(pivotality = 1) %>%
  select(situation,outcome,probability,pivotality,rating)

df.regression.gardener = responsibility.list[["gardener"]] %>%
  group_by(situation,outcome,probability,pivotality) %>%
  summarise(rating = mean(rating)) %>%
  ungroup() %>%
  select(situation,outcome,probability,pivotality,rating)

df.cv2 = rbind(df.regression.goalie,
               df.regression.spinner,
               df.regression.gardener) %>%
  mutate(experiment = rep(c("goalie","spinner","gardener"),each=8),
         probability = ifelse(experiment == "gardener", str_replace_all(probability,"suboptimal","30"),probability),
         probability = ifelse(experiment == "gardener", str_replace_all(probability,"optimal","70"),probability),
         probability = as.numeric(probability),
         probability = probability/100,
         outcome = str_replace_all(outcome,"wrong|negative","0"),
         outcome = str_replace_all(outcome,"correct|positive","1"),
         outcome = as.numeric(outcome),
         pivotality = str_replace_all(pivotality,"nonpivotal","0.5"),
         pivotality = str_replace_all(pivotality,"pivotal","1"),
         pivotality = as.numeric(pivotality),
         rating = ifelse(outcome == 0,-rating,rating),
         pivotality = ifelse(outcome == 0,-pivotality,pivotality),
         experiment = factor(experiment,levels = c("goalie", "spinner", "gardener"))
  ) %>%
  select(experiment,outcome,pivotality,probability,rating)

df.cv = df.cv %>%
  merge(.,df.cv2) %>%
  mutate(situation.index = rep(1:8,3),
         optimality = ifelse(probability > 0.4, 1, 0)) %>%
  arrange(experiment)

rm(list = setdiff(ls(),c("df.cv","df.data","inference.list", "responsibility.list")))

Functions

# softmax decision function
softmax = function(beta,theta,reward){
  exp(beta*theta[,1]*reward)/rowSums(exp(beta*theta*reward))
}

# agent decision functions
agent.bad = function(beta,df,reward){
  df = df %>%
    select(outcome,outcome_counterfactual)
  softmax(beta,1-df,reward)
}

agent.average = function(beta,df,reward){
  df = df %>%
    select(probability,probability_counterfactual)
  softmax(beta,df,reward)
}

agent.good = function(beta,df,reward){
  df = df %>%
    select(outcome,outcome_counterfactual)
  softmax(beta,df,reward)
}

fit_prior_cv = function(parameters){
  ngoalie = sum(df.tmp$experiment == "goalie")
  nspinner = sum(df.tmp$experiment == "spinner")
  ngardener = sum(df.tmp$experiment == "gardener")

  prior = matrix(c(rep(c(parameters[2],parameters[3],1-(parameters[2]+parameters[3])),ngoalie),
                   rep(c(parameters[4],parameters[5],1-(parameters[4]+parameters[5])),nspinner),
                   rep(c(parameters[6],parameters[7],1-(parameters[6]+parameters[7])),ngardener)),ncol=3,byrow=T)

  likelihood = cbind(
    agent.bad(parameters[1],df.tmp,reward),
    agent.average(parameters[1],df.tmp,reward),
    agent.good(parameters[1],df.tmp,reward)
  )
  posterior_fit = prior*likelihood/rowSums(prior*likelihood)
  posterior_data = select(df.tmp,contains("agent"))/100
  return(sum((posterior_data-posterior_fit)^2))
}

# output the predicted posteriors, based on the fit priors and beta
predicted_posterior = function(parameters){
  ngoalie = sum(df.tmp$experiment == "goalie")
  nspinner = sum(df.tmp$experiment == "spinner")
  ngardener = sum(df.tmp$experiment == "gardener")

  prior = matrix(c(rep(c(parameters[2],parameters[3],1-(parameters[2]+parameters[3])),ngoalie),
                   rep(c(parameters[4],parameters[5],1-(parameters[4]+parameters[5])),nspinner),
                   rep(c(parameters[6],parameters[7],1-(parameters[6]+parameters[7])),ngardener)),ncol=3,byrow=T)

  likelihood = cbind(
    agent.bad(parameters[1],df.tmp,reward),
    agent.average(parameters[1],df.tmp,reward),
    agent.good(parameters[1],df.tmp,reward)
  )
  posterior_fit = prior*likelihood/rowSums(prior*likelihood)
  return(posterior_fit)
}

expected_reward = function(){
  situations = list()
  exp_r_prior = list()
  exp_r_prior_scaled = list()

  # SITUATIONS: GOALIE & SPINNER
  situations[["goalie_spinner"]] = expand.grid(probability = seq(0.1,0.9,0.05), outcome = c(0,1)) %>%
    mutate(probability_counterfactual = 1-probability,
           outcome_counterfactual = 1-outcome,
           probability.overall = ifelse(outcome == 1, probability, probability_counterfactual),
           probability.overall = probability.overall/sum(probability.overall))

  # SITUATIONS: GARDENER
  situations[["gardener"]] = expand.grid(probability = seq(0.1,0.9,0.05), outcome = c(0,1),
                                         probability_counterfactual = seq(0.1,0.9,0.05), outcome_counterfactual = c(0,1)) %>%
    mutate(probability.overall = (1-abs(probability-outcome))*(1-abs(probability_counterfactual-outcome_counterfactual)),
           probability.overall = probability.overall/sum(probability.overall))

  situation.names = c("goalie_spinner","gardener")

  for(i in 1:length(situation.names)){
    #expected reward prior
    situations[[situation.names[i]]]$agent.bad = agent.bad(beta,situations[[situation.names[i]]],reward)
    situations[[situation.names[i]]]$agent.average = agent.average(beta,situations[[situation.names[i]]],reward)
    situations[[situation.names[i]]]$agent.good = agent.good(beta,situations[[situation.names[i]]],reward)

    exp_r_prior[[situation.names[i]]] = situations[[situation.names[i]]] %>%
      summarise(exp_r_prior.bad = sum((agent.bad*outcome + ((1-agent.bad)*outcome_counterfactual))*probability.overall),
                exp_r_prior.average = sum((agent.average*outcome + ((1-agent.average)*outcome_counterfactual))*probability.overall),
                exp_r_prior.good = sum((agent.good*outcome + ((1-agent.good)*outcome_counterfactual))*probability.overall))

    #rescale expected prior rewards
    exp_r_prior_scaled[[situation.names[i]]] = rescale(exp_r_prior[[situation.names[i]]] %>% as.numeric(),to=c(0,1),from = range(exp_r_prior[[situation.names[i]]]))
  }
  return(exp_r_prior_scaled)
}

difference_in_expected_reward = function(){
  inference.names = c(rep("goalie_spinner",2),"gardener")
  experiment.names = c("goalie","spinner","gardener")
  exp_r_difference.list = list()
  for(i in 1:length(inference.list)){
    df.tmp = inference.list[[i]] %>%
      group_by(situation,agent) %>%
      summarise(rating = mean(rating)/100) %>%
      ungroup %>%
      spread(agent,rating) %>%
      mutate(exp_r_prior_sum = sum(prior.list[[names(inference.list[i])]]*exp_r_prior_scaled[[inference.names[i]]]), #prior expected reward
             exp_r_posterior = rowSums(matrix(c(bad_agent,average_agent,good_agent),ncol=3) *
                                         matrix(unlist(rep(exp_r_prior_scaled[[inference.names[i]]],8)),byrow=T,ncol=3)), #posterior expected reward
             exp_r_difference = exp_r_posterior - exp_r_prior_sum#expected reward difference
             )

    exp_r_difference.list[[names(inference.list[i])]] = select(df.tmp,situation,exp_r_difference)
  }
  return(exp_r_difference.list)
}

Run (Models with difference in expectation)

reward = 1
# nsims = 100
nsims = 2

set.seed(5)
seeds = sample(nsims,nsims)

experiment.names = c("goalie", "spinner", "gardener")

cv.results = list()

# Models which include difference in expectations

for (i in 1:nsims){

  set.seed(seeds[i])
  cases = sample(1:nrow(df.cv))[1:(nrow(df.cv)/2)] #split half cross-validation

  df.training = df.cv %>%
    filter(row_number() %in% cases) %>%
    arrange(experiment)

  df.test = df.cv %>%
    filter(!row_number() %in% cases) %>%
    arrange(experiment)

  df.tmp = df.training
  fit.prior = optim(par = runif(7,0,1),fit_prior_cv,method = "L-BFGS-B",
                    lower = rep(0, 7), upper = c(100,rep(1,6)))

  beta = fit.prior$par[1]
  prior.list = list()
  prior.list[["goalie"]] = c(fit.prior$par[2],fit.prior$par[3],1-sum(fit.prior$par[2],fit.prior$par[3]))
  prior.list[["spinner"]] = c(fit.prior$par[4],fit.prior$par[5],1-sum(fit.prior$par[4],fit.prior$par[5]))
  prior.list[["gardener"]] = c(fit.prior$par[6],fit.prior$par[7],1-sum(fit.prior$par[6],fit.prior$par[7]))

  exp_r_prior_scaled = expected_reward()
  difference = difference_in_expected_reward()

  differences = data.frame()

  for (j in 1:length(difference)){
    tmp = difference[[j]]
    differences = rbind(differences,tmp)
  }

  differences = differences %>%
    mutate(experiment = rep(experiment.names,each=8),
           # situation.index = rep(1:8,4))
           situation.index = rep(1:8,3))

  df.tmp = df.tmp %>%
    merge(.,differences %>% select(experiment,situation.index,exp_r_difference))

  # CHOOSE WHICH MODEL TO RUN HERE
  # model = "rating~pivotality+exp_r_difference"
  model = "rating~exp_r_difference"
  fit.regression = lm(formula = paste0(model),data=df.tmp)

  model.fit = df.test %>%
    merge(.,differences %>% select(experiment,situation.index,exp_r_difference)) %>%
    mutate(prediction = predict(fit.regression,.)) %>%
    mutate_at(vars(rating,prediction),funs(abs)) %>%
    summarise(r = cor(rating,prediction),
              rmse = sqrt(mean((rating-prediction)^2))
              )

  cv.results[['parameters']][[i]] = as.data.frame(matrix(fit.prior$par,nrow=1)) %>%
  setNames(c('beta','goalie_bad','goalie_average','spinner_bad','spinner_average','gardener_bad','gardener_average'))
  cv.results[['regression']][[i]] = as.data.frame(matrix(fit.regression$coefficients,nrow=1))
  cv.results[['model']][[i]] = model.fit
  cv.results[['training']][[i]] = cases
}