---
title: "E1_Analysis"
output: pdf_document
date: "2023-06-29"
---

In this file, we complete the analysis for Experiment 1, which primarily investigates differences between modes of communication.

```{r setup, include=FALSE}

library(tidyverse)
library(ggdist)
library(lme4)
library(lmerTest)
library(MuMIn)

# Set random seed
set.seed(6)

#### define visualization attributes ####
theme_set(theme_minimal())
rational_colors <- c("Conservative" = "#EFD287", "Rational" = "#EFBA87", "Risky" = "#EF9387") 
mode_fill_colors <- c("Speech" = "#C1A6CE", "Text" = "#AAB2D0", "Visualization" = "#9DC3C6") 
mode_line_colors <- c("Speech" = "#9E78AF", "Text" = "#7E89B3", "Visualization" = "#6D9DA2") 

#### define factor levels ####
ageLevels <-  c("18 - 24", "25 - 34", "35 - 44", "45 - 54", "55 - 64", "65 - 74", "75 - 84", "85 or older", "Prefer not to answer")
educationLevels <-  c("Less than high school", "High school graduate", "Some college", "2 year degree", "4 year degree", "Professional degree", "Doctorate", "Prefer not to answer")
forecastLevels = c("quantile.20", "density.solid", "concrete.dist.text", "fuzzy.dist.text", "concrete.dist.audio", "fuzzy.dist.audio")
modalityLevels <-  c("Speech", "Text", "Visualization")
precisionLevels <- c("Concrete", "Fuzzy")
# further condition information
precision_dict <- list(
  'quantile.20' = 'Concrete',
  'density.solid' = 'Fuzzy',
  'concrete.dist.text' = 'Concrete',
  'fuzzy.dist.text' = 'Fuzzy',
  'concrete.dist.audio' = 'Concrete',
  'fuzzy.dist.audio' = 'Fuzzy'
)

#### add data for stimuli information ####
dists <- read.csv("Experiment 1 Stimuli/E1_stimuli_data.csv")

```

# Data Cleaning

First, we clean the data with the code below. The original data import contains Qualtrics information that is not necessary for our purposes. We also calculate several new columns, such as the graphical literacy. Since we are using mixed effects and linear models in the analysis, we also create factors for our variables of interest and demographics. 




```{r load data}

#### read in raw e1 data ####
e1_raw <- e1_raw <- read.csv("e1_data_qualtrics.csv")

#### select and calculate columns ####
# column definitions included here
e1 <- e1_raw %>% 
  select(
    timer.0_Page.Submit, # time on the introduction to forecasts
    timer.0_Click.Count, # click count on the introduction to forecasts
    
    # trial 1
    timer.1_Page.Submit, # time on the decision page, trial 1
    timer.1_Click.Count, # click count on the decision page, trial 1
    salt.1, # decision to salt, trial 1
    numeric_temp.1, # numeric estimation, trial 1
    confidence.1_1, # confidence rating, trial 1
    pct_likelihood.1_1, # likelihood of freezing, trial 1
    pct_tail.1_1, # likelihood of tail value (correct: 0), trial 1
    
    # trial 2
    timer.2_Page.Submit, # time on the decision page, trial 2
    timer.2_Click.Count, # click count on the decision page, trial 2
    salt.2, # decision to salt, trial 2
    numeric_temp.2, # numeric estimation, trial 2
    confidence.2_1, # confidence rating, trial 2
    pct_likelihood.2_1, # likelihood of freezing, trial 2
    pct_tail.2_1, # likelihood of tail value (correct: 0), trial 2
    
    # trial 3
    timer.3_Page.Submit, # time on the decision page, trial 3
    timer.3_Click.Count, # click count on the decision page, trial 3
    salt.3, # decision to salt, trial 3
    numeric_temp.3, # numeric estimation, trial 3
    confidence.3_1, # confidence rating, trial 3
    pct_likelihood.3_1, # likelihood of freezing, trial 3
    pct_tail.3_1, # likelihood of tail value (correct: 0), trial 3
    
    # trial 4
    timer.4_Page.Submit, # time on the decision page, trial 4
    timer.4_Click.Count, # click count on the decision page, trial 4
    salt.4, # decision to salt, trial 4
    numeric_temp.4, # numeric estimation, trial 4
    confidence.4_1, # confidence rating, trial 4
    pct_likelihood.4_1, # likelihood of freezing, trial 4
    pct_tail.4_1, # likelihood of tail value (correct: 0), trial 4
    
    # trial 5
    timer.5_Page.Submit, # time on the decision page, trial 5
    timer.5_Click.Count, # click count on the decision page, trial 5
    salt.5, # decision to salt, trial 5
    numeric_temp.5, # numeric estimation, trial 5
    confidence.5_1, # confidence rating, trial 5
    pct_likelihood.5_1, # likelihood of freezing, trial 5
    pct_tail.5_1, # likelihood of tail value (correct: 0), trial 5
    
    # trial 6
    timer.6_Page.Submit, # time on the decision page, trial 6
    timer.6_Click.Count, # click count on the decision page, trial 6
    salt.6, # decision to salt, trial 6
    numeric_temp.6, # numeric estimation, trial 6
    confidence.6_1, # confidence rating, trial 6
    pct_likelihood.6_1, # likelihood of freezing, trial 6
    pct_tail.6_1, # likelihood of tail value (correct: 0), trial 6
    
    # trial 7
    timer.7_Page.Submit, # time on the decision page, trial 7
    timer.7_Click.Count, # click count on the decision page, trial 7
    salt.7, # decision to salt, trial 7
    numeric_temp.7, # numeric estimation, trial 7
    confidence.7_1, # confidence rating, trial 7
    pct_likelihood.7_1, # likelihood of freezing, trial 7
    pct_tail.7_1, # likelihood of tail value (correct: 0), trial 7
    
    # trial 8
    timer.8_Page.Submit, # time on the decision page, trial 8
    timer.8_Click.Count, # click count on the decision page, trial 8
    salt.8, # decision to salt, trial 8
    numeric_temp.8, # numeric estimation, trial 8
    confidence.8_1, # confidence rating, trial 8
    pct_likelihood.8_1, # likelihood of freezing, trial 8
    pct_tail.8_1, # likelihood of tail value (correct: 0), trial 8
    
    # trial 9
    timer.9_Page.Submit, # time on the decision page, trial 9
    timer.9_Click.Count, # click count on the decision page, trial 9
    salt.9, # decision to salt, trial 9
    numeric_temp.9, # numeric estimation, trial 9
    confidence.9_1, # confidence rating, trial 9
    pct_likelihood.9_1, # likelihood of freezing, trial 9
    pct_tail.9_1, # likelihood of tail value (correct: 0), trial 9
    
    # trial 10
    timer.10_Page.Submit, # time on the decision page, trial 10
    timer.10_Click.Count, # click count on the decision page, trial 10
    salt.10, # decision to salt, trial 10
    numeric_temp.10, # numeric estimation, trial 10
    confidence.10_1, # confidence rating, trial 10
    pct_likelihood.10_1, # likelihood of freezing, trial 10
    pct_tail.10_1, # likelihood of tail value (correct: 0), trial 10
    
    # trial 11
    timer.11_Page.Submit, # time on the decision page, trial 11
    timer.11_Click.Count, # click count on the decision page, trial 11
    salt.11, # decision to salt, trial 11
    numeric_temp.11, # numeric estimation, trial 11
    confidence.11_1, # confidence rating, trial 11
    pct_likelihood.11_1, # likelihood of freezing, trial 11
    pct_tail.11_1, # likelihood of tail value (correct: 0), trial 11
    
    # trial 12
    timer.12_Page.Submit, # time on the decision page, trial 12
    timer.12_Click.Count, # click count on the decision page, trial 12
    salt.12, # decision to salt, trial 12
    numeric_temp.12, # numeric estimation, trial 12
    confidence.12_1, # confidence rating, trial 12
    pct_likelihood.12_1, # likelihood of freezing, trial 12
    pct_tail.12_1, # likelihood of tail value (correct: 0), trial 12
    
    # trust
    clear_1, # clarity rating
    accurate_1, # accuracy rating
    useful_1, # usefulness rating
    
    # graphical literacy
    lit_item1,
    lit_item2,
    lit_item3,
    lit_item4,
    
    # free response questions
    strategy, # strategy used by the pt for decision
    liked, # what pt liked about forecast
    disliked, # what pt disliked about forecast
    comments, # pt comments on the survey
    
    # demographics
    age, # pt age range
    education, # pt education level
    info_pref_1, # preference ranking for visuals
    info_pref_2, # preference ranking for text
    info_pref_3, # preference ranking for audio
    # embedded data
    Budget, # pt remaining budget
    Forecast, # pt forecast condition
    Cost, # the amount lost when salting roads
    Penalty, # the amount lost when roads freeze
    # audio specific
    PlayCount.0, # play count for introduction of forecast
    PlayCount.1, # play count for trial 1
    PlayCount.2, # play count for trial 2
    PlayCount.3, # play count for trial 3
    PlayCount.4, # play count for trial 4
    PlayCount.5, # play count for trial 5
    PlayCount.6, # play count for trial 6
    PlayCount.7, # play count for trial 7
    PlayCount.8, # play count for trial 8
    PlayCount.9, # play count for trial 9
    PlayCount.10, # play count for trial 10
    PlayCount.11, # play count for trial 11
    PlayCount.12 # play count for trial 12
  ) %>% 
  # calculate columns here
  mutate(
    # create anonymous participant id
    pt_id = 1:nrow(e1_raw), 
    # calculate graphical literacy
    graph_lit = if_else(lit_item1 < 27 & lit_item1 > 23, 1, 0) +
      if_else(lit_item2 == "They are equal", 1, 0) +
      if_else(lit_item3 == 20, 1, 0) + 
      if_else(lit_item4 == "Can't say", 1, 0),
    ageFactor = factor(age, levels = ageLevels),
    # create column for mode
    mode = case_when(
      str_detect(Forecast, "audio") ~ "Speech",
      str_detect(Forecast, "text") ~ "Text",
      TRUE ~ "Visualization"),
    # add in concrete and fuzzy information
    precision_level = case_when(
      Forecast %in% names(precision_dict)[precision_dict == 'Concrete'] ~ 'Concrete',
      Forecast %in% names(precision_dict)[precision_dict == 'Fuzzy'] ~ 'Fuzzy',
      TRUE ~ 'unknown',
    ),
    # create factor level columns for plotting and modeling
    forecastFactor = factor(Forecast, levels = forecastLevels), 
    modalityFactor = factor(mode, levels = modalityLevels),
    educationFactor = factor(education, levels = educationLevels),
    # rename the preference columns
    pref_visual = info_pref_1,
    pref_text = info_pref_2,
    pref_audio = info_pref_3,
    # create column for the most preferred
    pref1 = case_when(
      pref_visual == 1 ~ 'visual',
      pref_text == 1 ~ 'text',
      pref_audio == 1 ~ 'audio',
      TRUE ~ 'unknown',
    ), 
  ) %>% 
  # rename columns
  rename_with(~ sub("_Page.Submit", "", .), ends_with("_Page.Submit")) %>% # remove page submit from timer values
  rename_with(~ str_replace(., "_1$", ""), everything()) %>% # take out the _1 from the scale responses 
  rename_with(~ gsub("timer.(\\d+)_Click.Count", "click.\\1", .)) %>%  # change the click count variable names
  # calculate average trust
  rowwise() %>% 
  mutate(trust_average = mean(c(clear, accurate, useful))) %>% 
  ungroup() 

#### clean data into long format, perform additional exclusions ####
# exclusion based on tail percentage estimation, only easily accessible after the data has been transformed
e1_long <- e1 %>% 
  # move to long
  pivot_longer(
    cols = matches("\\.\\d+$"),  # select columns which have .digit in their names
    names_to = c(".value", "trial"),  # put the first part into value, the digit into trial
    names_pattern = "(.*)\\.(\\d+)$",  # old column names
    names_sep = NULL # no separation
  ) %>% 
  mutate(
    trial = as.integer(trial) # make trial an integer for joining purposes
  ) %>% 
  # add stimuli information
  left_join(dists[, c("trial", "freeze.pct", "mean", "crossover_distance")], by = "trial") %>% 
  # update format and calculate columns
  mutate(
    salt_code = ifelse(salt == "Yes, apply salt", 1, 0), # binary salt code for models
    EL = if_else(salt_code == 1, Cost, Penalty * freeze.pct), # the penalty multiplied by the probability of freezing for each trial they decided not to salt.cost if they decided to salt
    rational_EL = if_else(freeze.pct > Cost / Penalty, Cost, Penalty * freeze.pct), # if the percentage of freeze is greater than the ratio of cost:penalty, it is rational to pay the cost. otherwise, it is more rational overall to risk the freeze. cost:penalty is the point at which the cost of salting = the expected loss from penalty
    avg_temp = mean, # rename mean to average temp to better indicate what it represents
    numeric_estimate_diff = numeric_temp - mean, # calculate error in temperature estimate
    pct_estimate_diff = pct_likelihood - (freeze.pct *100), # calculate error in probability estimate
    # were choices rational?
    is_rational = ifelse(rational_EL == EL, 1,0),
    rational_code = ifelse(is_rational == 1, "Rational",
                           ifelse(salt_code == 1, "Conservative", "Risky")),
    # using expected loss
    EL_difference = EL - rational_EL,
    # opt_temp = avg_temp - optimal_temperature, # difference for each trial from optimal temperature
    # account for possible errors in playcount by checking with overall click count from timer variable
    PlayCount = ifelse(PlayCount == 0, ifelse(click - 2 > 0, click -2, 1), PlayCount)
  ) 

e1_long <- e1_long[!(e1_long$pt_id %in% e1_long$pt_id[e1_long$pct_tail >= 50]), ] # exclude responses which gave > 50% likelihood for tail value

#### select columns for analysis ####
e1_main <- e1_long %>% 
  select(
    pt_id,
    forecastFactor, # pt forecast condition
    modalityFactor, # mode of condition
    precision_level, # precision level of the stimuli
    trial,
    timer, # time to decision 
    click, # number of clicks on page
    PlayCount, # count of audio plays
    salt_code, # decision to salt
    numeric_temp, # numeric estimation
    confidence, # confidence rating
    pct_likelihood, # likelihood of freezing
    pct_tail, # likelihood of tail value (correct: 0)
    clear, # clarity rating
    accurate, # accuracy rating
    useful, # usefulness rating
    trust_average, # average trust rating
    graph_lit, 
    freeze.pct, # actual likelihood of freezing
    # calculated columns
    EL, 
    rational_EL, 
    EL_difference,
    avg_temp,
    crossover_distance, # distance from the mean of the distribution to the optimal crossover
    is_rational,
    rational_code,
    numeric_estimate_diff,
    pct_estimate_diff,
    # demographics
    ageFactor, # pt age range
    educationFactor, # pt education level
    pref1, # mode rated first
    # embedded data
    Budget # pt remaining budget
  ) %>% 
  mutate(
    trial = as.factor(trial) # make trial a factor for modeling purposes
  )

#### create data frame for participant information ####
pt_e1 <- e1 %>% 
  select(
    pt_id,
    ageFactor, # pt age range
    educationFactor, # pt education level
    graph_lit, # graph_lit score
    forecastFactor, # specific forecast viewed
    modalityFactor, # mode of condition
    precision_level, # precision level of the stimuli,
    trust_average, # trust in the forecast (model DV)
    pref1 # top preference for information
  ) %>% 
  left_join(
    e1_long %>% group_by(pt_id) %>%  summarise(rationality = mean(is_rational, na.rm = TRUE))
  ) %>% 
  filter(pt_id %in% e1_long$pt_id)

```

# Visualize

Now, we visualize preliminary results & information, starting with demographic and other related information.

```{r demographics}

# participant age
ggplot(pt_e1, aes(x = ageFactor))+
  geom_bar(fill = "#aaaaaa")+
  scale_x_discrete(limits = ageLevels, labels = c("18-24", "25-34", "35-44", "45-54", "55-64", "65-74", "75-84", ">85", "No answer\n"))+
  labs(
    x = "Age Range",
    y = "Count",
    title = "Distribution of Participant Ages"
  )

# participant education
ggplot(pt_e1, aes(x = educationFactor))+
  geom_bar(fill = "#aaaaaa")+
  scale_x_discrete(limits = educationLevels, labels = c("Less than\nhigh school\n", "High school\ngraduate", "Some\ncollege", "2 year\ndegree", "4 year\ndegree", "Professional\ndegree", "Doctorate", "No answer"))+
  labs(
    x = "Education Level",
    y = "Count",
    title = "Distribution of Participant Education"
  )

# literacy scores 
ggplot(pt_e1, aes(x = as.factor(graph_lit), fill = as.factor(graph_lit)))+
  geom_bar()+
  scale_fill_brewer()+
  labs(
    x = "Graphical Literacy Score",
    y = NULL,
    title = "Distribution of Literacy Scores",
    fill = "Graphical Literacy Score"
  )

# preference for each Mode
prefs_e1 <- e1  %>% 
  filter(pt_id %in% e1_long$pt_id) %>%
  pivot_longer(
    cols = starts_with("pref_"),
    names_to = "source_mode",
    values_to = "rank"
  ) %>%
  filter(!is.na(rank)) %>% 
  select(pt_id, forecastFactor, modalityFactor, source_mode, rank)
ggplot(prefs_e1, aes(x = rank, fill = as.factor(source_mode)))+
  geom_bar(stat = "count", position = "dodge")+
  scale_fill_manual(values = mode_fill_colors) + # Use manually defined colors
  scale_color_manual(values = mode_line_colors) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nRank",
    title = "Mode Rankings\n",
    fill = "Mode"
  )

# play counts between audio variants
ggplot(e1_main %>% filter(!is.na(PlayCount)), aes(x = PlayCount))+
  geom_bar()+
  labs(x = "Play count", y = "Count", title = "Play Count for Speech Forecasts")

```

Now, we want to examine participant decisions overall - what kinds of decisions did they make and how confident were they? We consider three types of decisions: Conservative (it was rational to not salt and they salted), Rational (it was rational to salt and they salted OR it was rational not to salt and they did not salt), and Risky (it was rational to salt and they did not salt).

```{r decision}

# decision types 
decision_interaction_e1 <- e1_main %>%
  filter(trial != 0) %>% 
  group_by(modalityFactor, precision_level, rational_code) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(modalityFactor, precision_level) %>%
  mutate(total = sum(n), 
         pct = n / total)
ggplot(decision_interaction_e1, aes(x = rational_code, y = pct, fill = rational_code))+
  geom_col(position = "dodge")+
  facet_grid(modalityFactor ~ precision_level)+
  scale_fill_manual(values = rational_colors) + # Use manually defined colors
  labs(x = "\nDecision Type",
       fill = "Decision Type",
       y = "Proportion of Decisions\n",
       title = "Decision Types by Mode and Precision Level")+
  theme_minimal()


# confidence in choice 
ggplot(e1_main %>% filter(trial != 0), aes(x = confidence, fill = modalityFactor, color = modalityFactor))+
  geom_histogram(binwidth =5)+
  stat_slabinterval(size = 1)+
  scale_x_continuous(limits = c(50,100))+
  scale_y_continuous(n.breaks = 3)+
  scale_fill_manual(values = mode_fill_colors) + # Use manually defined colors
  scale_color_manual(values = mode_line_colors) + # Use manually defined colors
  facet_grid(modalityFactor~precision_level)+
  labs(
    y = "Count\n",
    x = "\nDecision Confidence",
    title = "Confidence in Decisions by Mode and Precision",
    fill = "Forecast Mode",
    # subtitle = "No differences reach significance"
  )+
  theme_minimal()+
  theme(legend.position = "none",
        axis.title = element_text(size = 10), # Axis titles
        axis.text = element_text(size = 8), # Axis text
        plot.title = element_text(size = 12),
        strip.text = element_text(size = 8),
        strip.text.y = element_text(angle = 0, hjust = 0)) # Remove the legend

ggsave("forecast_confidence_ratings.png", width = 5, height = 3, dpi = 300)


```

Now, we look at trust ratings overall, also examining the individual components.

```{r trust ratings}


# ratings for clarity, accuracy, usefulness (and calculated overall trust)
trust_e1 <- e1 %>%
  pivot_longer(cols = c(clear, accurate, useful, trust_average),
               names_to = "rating_type",
               values_to = "rating_value") %>% 
  select(rating_type, rating_value, modalityFactor, precision_level, forecastFactor)

# average clarity
ggplot(trust_e1 %>% filter(rating_type == "clear"), aes(x = rating_value, fill = modalityFactor, color = modalityFactor))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  scale_fill_manual(values = mode_fill_colors) + # Use manually defined colors
  scale_color_manual(values = mode_line_colors) + # Use manually defined colors
  facet_grid(modalityFactor~precision_level)+
  labs(
    y = "Count\n",
    x = "\nAverage Forecast Clarity",
    title = "Distribution of Forecast Clarity Ratings",
    fill = "Forecast Mode",
  )+
  theme_minimal()+
  theme(legend.position = "none") # Remove the legend


# average accuracy
ggplot(trust_e1 %>% filter(rating_type == "accurate"), aes(x = rating_value, fill = modalityFactor, color = modalityFactor))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  scale_fill_manual(values = mode_fill_colors) + # Use manually defined colors
  scale_color_manual(values = mode_line_colors) + # Use manually defined colors
  facet_grid(modalityFactor~precision_level)+
  labs(
    y = "Count\n",
    x = "\nAverage Forecast Accuracy",
    title = "Distribution of Forecast Accuracy Ratings",
    fill = "Forecast Mode",
  )+
  theme_minimal()+
  theme(legend.position = "none") # Remove the legend

# average usefulness
ggplot(trust_e1 %>% filter(rating_type == "useful"), aes(x = rating_value, fill = modalityFactor, color = modalityFactor))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  scale_fill_manual(values = mode_fill_colors) + # Use manually defined colors
  scale_color_manual(values = mode_line_colors) + # Use manually defined colors
  facet_grid(modalityFactor~precision_level)+
  labs(
    y = "Count\n",
    x = "\nAverage Forecast Usefulness",
    title = "Distribution of Forecast Usefulness Ratings",
    fill = "Forecast Mode",
  )+
  theme_minimal()+
  theme(legend.position = "none") # Remove the legend

# average trust (mean values of clarity, accurate, and useful)
# plot <- 
ggplot(e1_main, aes(x = trust_average, fill = modalityFactor, color = modalityFactor))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  scale_fill_manual(values = mode_fill_colors) + # Use manually defined colors
  scale_color_manual(values = mode_line_colors) + # Use manually defined colors
  facet_grid(modalityFactor~precision_level)+
  scale_y_continuous(n.breaks = 3)+
  labs(
    y = "Count\n",
    x = "\nAverage Trust in Forecast",
    title = "Distribution of Forecast Trust Ratings",
    fill = "Forecast Mode",
  )+
  theme_minimal()
#   theme(legend.position = "none",
#         axis.title = element_text(size = 10), # Axis titles
#         axis.text = element_text(size = 8), # Axis text
#         plot.title = element_text(size = 12),
#         strip.text = element_text(size = 8)
#   ) # Remove the legend
# # 
# ggsave("forecast_trust_ratings.png", plot, width = 5, height = 3, dpi = 300)


```

We also quickly examine how accurate participants were in their estimations overall. 

```{r accuracy in estimations}


# numeric temp estimate error by forecast variant
ggplot(e1_main %>% filter(trial != 0), aes(x = numeric_estimate_diff, fill = modalityFactor))+
  geom_bar()+
  facet_grid(modalityFactor ~ .)+
  scale_fill_manual(values = mode_fill_colors) + # Use manually defined colors
  scale_color_manual(values = mode_line_colors) + # Use manually defined colors
  xlim(-10, 10)+
  labs(
    y = "Forecast Variant\n",
    x = "\nDifference between Actual\nand Estimated Numeric Temperature",
    title = "Numeric Temperature Estimation Accuracy by Mode",
    fill = "Forecast Mode"
  )+
  theme_minimal()+
  theme(legend.position = "none") # Remove the legend

# difference between pct_likelihood and actual likelihood by forecast variant
ggplot(e1_main %>% filter(trial != 0), aes(y = modalityFactor, x = pct_estimate_diff, fill = modalityFactor, color = modalityFactor))+
  stat_halfeye(size =  1)+
  scale_fill_manual(values = mode_fill_colors) + # Use manually defined colors
  scale_color_manual(values = mode_line_colors) + # Use manually defined colors
  labs(
    y = "Forecast Variant\n",
    x = "\nDifference between Actual\n and Estimated Freeze Likelihood (%)",
    title = "Freeze Likelihood Estimation Accuracy by Mode",
    fill = "Forecast Mode"
  )+
  theme_minimal()+
  theme(legend.position = "none") # Remove the legend

```

# Analyze

The following sections address each attribute of decision making in turn.

## Decision to Salt

So here, we're seeing how the forecasts differ in their crossover points. At what temperature do participants make the switch to not salt or salt? How does it differ from the rational crossover point?

Below, we run a generalized linear mixed-effects model which will give us the average difference in crossover point from the rational. This result shows how the variation in fuzziness affect the point at which participants begin to salt the road (e.g., when the temperature is 32º? when it's 33º?). We can compare this to where they *should* begin to salt the road by calculating difference from optimal crossover temperature.

We also test models with different referents, since our main factors are categorical. 


```{r select salt choice model}

#### stepwise models ####

# baseline - does salt likelihood change with temperature? it should. 
salt0 <- glmer(salt_code ~ crossover_distance + (1|pt_id), 
               data = e1_main,
               family = binomial("logit")) 
# add in Mode - is salt likelihood different for the different modalities
salt1 <- glmer(salt_code ~ crossover_distance + modalityFactor + (1|pt_id), 
               data = e1_main,
               family = binomial("logit")) 
# add in precision - is salt likelihood different for the different precision levels
salt2 <- glmer(salt_code ~ crossover_distance + modalityFactor + precision_level +(1|pt_id), 
               data = e1_main,
               family = binomial("logit")) 
# add in interaction - is there an interaction between Mode and precision 
salt3 <- glmer(salt_code ~ crossover_distance + modalityFactor * precision_level + (1|pt_id), 
               data = e1_main,
               family = binomial("logit"))
# most complex model, not included in comparisons but was pre-registered to run
salt4 <- glmer(salt_code ~ crossover_distance + modalityFactor * precision_level + ageFactor + educationFactor + pref1 + (1|pt_id), 
               data = e1_main,
               nAGQ = 0,
               family = binomial("logit"))

#### select and interpret model ####

# find model with best fit
anova(salt0, salt1, salt2, salt3)

summary(salt1) # summarize model
r.squaredGLMM(salt1)

# text baseline
salt1_text <- glmer(salt_code ~ crossover_distance + factor(modalityFactor, levels = c("Text", "Visualization", "Speech")) + (1|pt_id), 
                    data = e1_main,
                    family = binomial("logit")) 
summary(salt1_text)

# visual baseline
salt1_visual <- glmer(salt_code ~ crossover_distance + factor(modalityFactor, levels = c("Visualization", "Text", "Speech")) + (1|pt_id), 
                      data = e1_main,
                      family = binomial("logit")) 
summary(salt1_visual)


# get model coefs 
coef_salt <- fixef(salt2)

### print out crossover temperatures
for(mode in modalityLevels) {
  for(precision in precisionLevels) {
    
    # Get the coefficients relevant for this combination of levels
    intercept <- coef_salt[1]
    condition <- sum(c(coef_salt[paste0("modalityFactor", mode, sep="")],
                       coef_salt[paste0("precision_level", precision, sep="")],
                       coef_salt[paste0("modalityFactor", mode, ":precision_level", precision, sep="")]),
                     na.rm= TRUE)
    temp_coef <- coef_salt[2]
    
    # Solve for temperature: logit(0.5) = 0 = intercept + temp_coef * temperature
    temperature <- ((-1 * intercept) - condition)  / temp_coef
    
    cat("For mode =", mode, "and precision =", precision, ", the crossover temperature is", temperature, "from optimal.\n")
  }
}

```


## Rationality

This analysis pairs with the crossover temperature analysis. Here, rather than look at what temperature participants switched from not salting to salting, we look at the overall count of the kinds of decisions (rational, risky, or conservative). Crossover temperatures for two conditions could be the same if one condition was highly rational and one condition was not frequently rational but had an even spread of risky and conservative decisions. This analysis helps to clear this up and point more clearly to variations in rationality.


```{r decision rationality}

#### precision ####

precision_table <- table(e1_main$precision_level, 
                         e1_main$rational_code)
precision_table
chisq.test(precision_table)


#### mode ####

modality_table <- table(e1_main$modalityFactor, 
                        e1_main$rational_code)
modality_table
modality_chisq <- chisq.test(modality_table)
modality_chisq
modality_chisq$stdres

```


## Confidence

Next, we want to look at how confidence participants were in the decisions they made. Here, we look to if and how confidence reliably varied between forecasts.

```{r select confidence model}

# baseline - does confidence change with temperature difference from crossover? it should. 
confidence0 <- lmer(confidence ~ crossover_distance + (1|pt_id) , 
                    data = e1_main) 

# add in rationality - is confidence different for rational decisions
confidence1 <- lmer(confidence ~ crossover_distance + is_rational + (1|pt_id) , 
                    data = e1_main) 

# add in mode - is confidence different for the different modalities
confidence2 <- lmer(confidence ~ crossover_distance + is_rational + modalityFactor + (1|pt_id) , 
                    data = e1_main) 

# add in precision - is confidence different for the different precision levels
confidence3 <- lmer(confidence ~ crossover_distance + is_rational + modalityFactor + precision_level + (1|pt_id) , 
                    data = e1_main) 

# add in interaction - is there an interaction between mode and precision 
confidence4 <- lmer(confidence ~ crossover_distance + is_rational + modalityFactor * precision_level + is_rational + (1|pt_id) , 
                    data = e1_main) 
# most complex model, not included in comparisons but was pre-registered to run
confidence5 <- lmer(confidence ~ crossover_distance + modalityFactor * precision_level + is_rational + ageFactor + educationFactor + pref1 + (1|pt_id) , 
                    data = e1_main) 

#### select and interpret model #### 

# find model with best fit
anova(confidence0, confidence1, confidence2, confidence3, confidence4)

# run additional models for different pairwise comparisons
# text baseline
confidence2_text <- lmer(confidence ~ crossover_distance + is_rational + factor(modalityFactor, levels = c("Text", "Visualization", "Speech")) + (1|pt_id) , 
                         data = e1_main)
# visual baseline
confidence2_visual <- lmer(confidence ~ crossover_distance + is_rational + factor(modalityFactor, levels = c("Visualization", "Text", "Speech")) + (1|pt_id) , 
                           data = e1_main)


summary(confidence2_text)

# get model coefs
confidence_coef <- fixef(confidence3)

### print out confidence for each
for(mode in modalityLevels) {
  for(precision in precisionLevels) {
    
    # Get the coefficients relevant for this combination of levels
    confidence <- sum(c(confidence_coef[1], confidence_coef[paste0("modalityFactor", mode, sep="")],
                        confidence_coef[paste0("precision_level", precision, sep="")],
                        confidence_coef[paste0("modalityFactor", mode, ":precision_level", precision, sep="")]),
                      na.rm= TRUE)
    
    cat("For", precision, mode, "the confidence is", confidence, "\n")
  }
}
```


## Trust

Finally, we look at measure of overall trust to determine if and how it systematically varies by how fuzzy the forecast is. Trust is an aggregate measure of clarity, usefulness, and accuracy.


```{r select trust model}

# baseline - does trust change with how rational participants were?
trust0 <- lm(trust_average ~ rationality, 
             data = pt_e1)
# add in mode - is trust different for the different modalities
trust1 <- lm(trust_average ~ rationality + modalityFactor, 
             data = pt_e1)

# add in precision - is trust different for the different precision levels
trust2 <- lm(trust_average ~ rationality + modalityFactor + precision_level, 
             data = pt_e1)

# add in interaction - is there an interaction between mode and precision 
trust3 <- lm(trust_average ~ rationality + modalityFactor * precision_level, 
             data = pt_e1)
# most complex model, not included in comparisons but was pre-registered to run
trust4 <- lm(trust_average ~ rationality + modalityFactor * precision_level + graph_lit +  ageFactor + educationFactor + pref1, 
             data = pt_e1)

#### select and interpret model #### 

# find model with best fit
anova(trust0, trust1, trust2, trust3)

# run additional models for different pairwise comparisons
# text baseline
trust1_visual <- lm(trust_average ~ rationality + factor(modalityFactor, levels = c("Visualization", "Text", "Speech")), 
                    data = pt_e1)
trust1_text <- lm(trust_average ~ rationality + factor(modalityFactor, levels = c("Text", "Visualization", "Speech")), 
                  data = pt_e1)

summary(trust1)

trust_coef <- coef(trust2)

# # For each combination of levels of modalityFactor and precision_level, show predicted trust level
for(mode in modalityLevels) {
  for(precision in precisionLevels) {
    
    # Get the coefficients relevant for this combination of levels
    trust <- sum(c(trust_coef[1], trust_coef[paste0("modalityFactor", mode, sep="")],
                   trust_coef[paste0("precision_level", precision, sep="")],
                   trust_coef[paste0("modalityFactor", mode, ":precision_level", precision, sep="")]),
                 na.rm= TRUE)
    
    cat("For", precision, mode, "the trust, on average, is", trust, "out of 10 \n")
  }
}

```



