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

In this file, we complete the analysis for each mode of information presentation, along with several exploratory analyses. 

```{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") 
precision_fill_colors <- c("Most Concrete" = "#8DC66C", "Somewhat Concrete" = "#AFD272", "Mixed" = "#C8DB77", "Somewhat Fuzzy" = "#DDE17A", "Most Fuzzy" = "#E3DA7B")
precision_line_colors <- c("Most Concrete" = "#6DA53F", "Somewhat Concrete" = "#84AB41", "Mixed" = "#A0B344", "Somewhat Fuzzy" = "#B4B846",  "Most Fuzzy" = "#B9AF46")

#### 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")
# further condition information
precision_dict <- list(
  'forecast_1' = 'Most Concrete',
  'forecast_2' = 'Somewhat Concrete',
  'forecast_3' = 'Mixed',
  'forecast_4' = 'Somewhat Fuzzy',
  'forecast_5' = 'Most Fuzzy'
)
# convert the list to a named vector for easy subsetting
precision_vector <- unlist(precision_dict)
precision_levels <- c('Most Concrete', 'Somewhat Concrete', 'Mixed','Somewhat Fuzzy', 'Most Fuzzy')

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

```

# Speech

We'll start by analyzing the speech survey. Here, we compared variations of speech forecasts, including different delays, pauses, and inflections.

## Speech 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 speech load data}

#### read in raw e2_speech data ####
e2_speech_raw <- e2_speech_raw <- read.csv("e2_speech_qualtrics.csv") 

#### select and calculate columns ####
# column definitions included here
e2_speech <- e2_speech_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 speech
    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(e2_speech_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),
    # create column for mode
    mode = "Speech",
    # add in concrete and fuzzy information
    precision_level =  precision_vector[Forecast],
    # create factor level columns for plotting and modeling
    educationFactor = factor(education, levels = educationLevels),
    ageFactor = factor(age, levels = ageLevels), 
    # rename the preference columns to make sense
    pref_visual = info_pref_1,
    pref_text = info_pref_2,
    pref_audio = info_pref_3,
    pref1 = case_when(
      pref_visual == 1 ~ 'visual',
      pref_text == 1 ~ 'text',
      pref_audio == 1 ~ 'audio',
      TRUE ~ 'unknown',
    ), # make single column for preference - what was the most preferred? 
  ) %>% 
  # 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 exclusion ####
# exclusion based on tail percentage estimation, only easily accessible after the data has been transformed
e2_speech_long <- e2_speech %>% 
  # 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 = round(pct_likelihood - (freeze.pct *100), 0), # 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")),
    EL_difference = EL - rational_EL,
    # crossover_distance = 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)
  ) 

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

#### select columns for analysis ####
e2_speech_main <- e2_speech_long %>% 
  select(
    pt_id,
    mode, # 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
    precision_level = factor(precision_level, levels = c("Most Concrete", "Somewhat Concrete", "Mixed", "Somewhat Fuzzy", "Most Fuzzy"))
  )

#### create data frame for participant information ####
pt_e2_speech <- e2_speech %>% 
  select(
    pt_id,
    Budget,
    ageFactor, # pt age range
    educationFactor, # pt education level
    graph_lit, # graph_lit score
    mode, # 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(
    e2_speech_long %>% group_by(pt_id) %>%  summarise(rationality = mean(is_rational, na.rm = TRUE))
  ) %>% 
  filter(pt_id %in% e2_speech_long$pt_id)

```

## Speech Visualize

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

```{r speech demographics}

# participant age
ggplot(pt_e2_speech, 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 [Speech]"
  )

# participant education
ggplot(pt_e2_speech, 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 [Speech]"
  )

# literacy scores 
ggplot(pt_e2_speech, 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 [Speech]",
    fill = "Graphical Literacy Score"
  )

# preference for each Mode
prefs_e2_speech <- e2_speech %>%
  pivot_longer(
    cols = starts_with("pref_"),
    names_to = "source_mode",
    values_to = "rank"
  ) %>%
  filter(!is.na(rank)) %>% 
  select(pt_id, precision_level, mode, source_mode, rank)
ggplot(prefs_e2_speech, aes(x = rank, fill = as.factor(source_mode)))+
  geom_bar(position = "dodge")+
  scale_fill_manual(values = c("pref_audio" = "#C1A6CE", "pref_text" = "#AAB2D0", "pref_visual" = "#9DC3C6") ) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nRank",
    title = "Mode Rankings [Speech]\n",
    fill = "Mode"
  )


# play counts between audio variants
ggplot(e2_speech_main %>% filter(trial != 0), aes(x = PlayCount))+
  geom_bar()+
  scale_x_continuous(breaks = seq(0, 13, 1))+
  facet_grid(precision_level ~ .)+
  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 speech decision}

# decision type by precision
decision_precision_e2_speech <- e2_speech_main %>%
  filter(trial != 0) %>% 
  group_by(precision_level, rational_code) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(precision_level) %>%
  mutate(total = sum(n), 
         pct = n / total)
ggplot(decision_precision_e2_speech, aes(x = rational_code, y = pct, fill = rational_code))+
  geom_col(position = "dodge")+
  scale_fill_manual(values = rational_colors) + # Use manually defined colors
  facet_grid(.~ precision_level)+
  scale_x_discrete(labels = c("Cons.", "Ration.", "Risky"))+
  labs(x = "\nPrecision Level",
       fill = "Decision Type",
       y = "Proportion of Decisions\n",
       title = "Decision Types by Precision Level")+
  ylim(c(0, 1))+
  theme(legend.position = "none") # Remove the legend


# confidence in choice by forecast variant
ggplot(e2_speech_main %>% filter(trial != 0), aes(x = confidence, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =5)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_x_continuous(limits = c(50,100))+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nDecision Confidence",
    title = "Confidence in Decisions by Forecast",
    fill = "Precision Level"
  )+
  theme(legend.position = "none",
        strip.text.y = element_text(angle = 0, hjust = 0)) # Remove the legend

```


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


```{r speech trust ratings}

# ratings for clarity, accuracy, usefulness (and calculated overall trust)
trust_e2_speech <- e2_speech %>%
  pivot_longer(cols = c(clear, accurate, useful, trust_average),
               names_to = "rating_type",
               values_to = "rating_value") %>% 
  select(rating_type, rating_value, precision_level) %>% 
  mutate(
    precision_level = factor(precision_level, levels = c("Most Concrete", "Somewhat Concrete", "Mixed", "Somewhat Fuzzy", "Most Fuzzy"))
  )


# clarity by forecast variant
ggplot(trust_e2_speech %>% filter(rating_type == "clear"), aes( x = rating_value, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Precision Level\n",
    x = "\nForecast Attribute",
    title = "Distribution of Forecast Clarity Ratings",
    fill = "Precision Level",
    color = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend

# accuracy by forecast variant
ggplot(trust_e2_speech %>% filter(rating_type == "accurate"), aes( x = rating_value, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Precision Level\n",
    x = "\nForecast Attribute",
    title = "Distribution of Forecast Accuracy Ratings",
    fill = "Precision Level",
    color = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend

# useful by forecast variant
ggplot(trust_e2_speech %>% filter(rating_type == "useful"), aes( x = rating_value, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Precision Level\n",
    x = "\nForecast Attribute",
    title = "Distribution of Forecast Usefulness Ratings",
    fill = "Precision Level",
    color = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend

# average trust in forecast variant (mean values of clarity, accurate, and useful)
ggplot(e2_speech_main, aes(x = trust_average, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Precision Level\n",
    x = "\nForecast Attribute",
    title = "Distribution of Forecast Trust Ratings",
    fill = "Precision Level",
    color = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend


```

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


```{r speech accuracy in estimations}

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

# difference between pct_likelihood and actual likelihood by forecast variant
ggplot(e2_speech_main %>% filter(trial != 0), aes(x = pct_estimate_diff, fill = precision_level))+
  geom_bar()+
  facet_grid(precision_level ~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Precision Level\n",
    x = "\nDifference between Actual\n and Estimated Freeze Likelihood (%)",
    title = "Freeze Likelihood Estimation Accuracy by Forecast",
    fill = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend


```


## Speech Analyze

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

If there's an indication that the level of precision may have had an effect, we run additional models with varying referents, as our main factor of interest is categorical.

### Speech Crossover Temperature

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.


```{r speech salt choice model}

#### stepwise models ####

# baseline 
salt0_speech <- glmer(salt_code ~ crossover_distance + (1|pt_id) , 
                      data = e2_speech_main,
                      family = binomial("logit")) 

# add in precision 
salt1_speech <- glmer(salt_code ~ crossover_distance + precision_level + (1|pt_id) , 
                      data = e2_speech_main,
                      family = binomial("logit")) 

# most complex model - not included in comparisons but was pre-registered to run
salt2_speech <- glmer(salt_code ~ crossover_distance + precision_level + ageFactor + educationFactor + pref1 + (1|pt_id) , 
                      data = e2_speech_main,
                      nAGQ = 0,
                      family = binomial("logit"))



#### select and run model(s) ####

# find model with best fit
anova(salt0_speech, salt1_speech)


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


# get model coefs and errors
coef_salt_speech <- fixef(salt1_speech)

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

```


### Speech 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 speech decision rationality}

precision_table_speech <- table(e2_speech_main$precision_level, 
                                e2_speech_main$rational_code)
precision_table_speech
speech_chisq <- chisq.test(precision_table_speech)
speech_chisq

```


### Speech Confidence

Next, we want to look at how confidence participants were in the decisions they made. Here, we look to see if confidence reliably varied between different levels of fuzziness. 

```{r speech confidence model}

#### stepwise models ####

e2_speech_main$precision_level <-  factor(e2_speech_main$precision_level, levels = c('Most Concrete', 'Somewhat Concrete', 'Mixed','Somewhat Fuzzy', 'Most Fuzzy'))

# baseline
confidence0_speech <- lmer(confidence ~ crossover_distance + (1|pt_id) , 
                           data = e2_speech_main) 

# add in rationality
confidence1_speech <- lmer(confidence ~ crossover_distance + is_rational + (1|pt_id) , 
                           data = e2_speech_main) 

# add in precision 
confidence2_speech <- lmer(confidence ~ crossover_distance + is_rational + precision_level + (1|pt_id) , 
                           data = e2_speech_main) 

# most complex model - not included in comparisons but was pre-registered to run
confidence3_speech <- lmer(confidence ~ crossover_distance + precision_level + is_rational + ageFactor + educationFactor + pref1 + (1|pt_id) , 
                           data = e2_speech_main)

#### select and run model(s) ####

# find model with best fit
anova(confidence0_speech, confidence1_speech, confidence2_speech)

# conduct model
summary(confidence1_speech)
r.squaredGLMM(confidence1_speech)

confidence_coef_speech <- fixef(confidence2_speech)

for(precision in precision_levels) {
  
  # Get the coefficients relevant for this combination of levels
  confidence <- sum(c(confidence_coef_speech[1],
                      confidence_coef_speech[paste0("precision_level", precision, sep="")]),
                    na.rm= TRUE)
  
  cat("For", precision, "the confidence is", confidence, "\n")
}


```


### Speech 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 speech select trust model}

#### stepwise models ####

# baseline 
trust0_speech <- lm(trust_average ~ rationality, 
                    data = pt_e2_speech)
# add in precision 
trust1_speech <- lm(trust_average ~ rationality + precision_level, 
                    data = pt_e2_speech)

# most complex model - not included in comparisons but was pre-registered to run
trust2_speech <- lm(trust_average ~ rationality + precision_level +  ageFactor + educationFactor + pref1, 
                    data = pt_e2_speech)

#### select and run model(s)

# find model with best fit
anova(trust0_speech, trust1_speech)

# conduct model
summary(trust0_speech)

```


# Text

Now, we move on to complete the same analysis with the text survey.

## Text 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 text load data}

#### read in raw e2_text data ####
e2_text_raw <- e2_text_raw <- read.csv("e2_text_qualtrics.csv")

#### select and calculate columns ####
# column definitions included here
e2_text <- e2_text_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
  ) %>% 
  # calculate columns here
  mutate(
    # create anonymous participant id
    pt_id = 1:nrow(e2_text_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),
    # create column for mode
    mode = "Text",
    # add in concrete and fuzzy information
    precision_level =  precision_vector[Forecast],
    # create factor level columns for plotting and modeling
    educationFactor = factor(education, levels = educationLevels),
    ageFactor = factor(age, levels = ageLevels), 
    # rename the preference columns to make sense
    pref_visual = info_pref_1,
    pref_text = info_pref_2,
    pref_audio = info_pref_3,
    pref1 = case_when(
      pref_visual == 1 ~ 'visual',
      pref_text == 1 ~ 'text',
      pref_audio == 1 ~ 'audio',
      TRUE ~ 'unknown',
    ), # make single column for preference - what was the most preferred? 
  ) %>% 
  # 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 exclusion ####
# exclusion based on tail percentage estimation, only easily accessible after the data has been transformed
e2_text_long <- e2_text %>% 
  # 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 = round(pct_likelihood - (freeze.pct *100), 0), # 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")),
    EL_difference = EL - rational_EL,
    # crossover_distance = avg_temp - optimal_temperature, # difference for each trial from optimal temperature
    PlayCount = NA
  ) 

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


#### select columns for analysis ####
e2_text_main <- e2_text_long %>% 
  select(
    pt_id,
    mode, # 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
    precision_level = factor(precision_level, levels = c("Most Concrete", "Somewhat Concrete", "Mixed", "Somewhat Fuzzy", "Most Fuzzy"))
  )

#### create data frame for participant information ####
pt_e2_text <- e2_text %>% 
  select(
    pt_id,
    Budget,
    ageFactor, # pt age range
    educationFactor, # pt education level
    graph_lit, # graph_lit score
    mode, # 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(
    e2_text_long %>% group_by(pt_id) %>%  summarise(rationality = mean(is_rational, na.rm = TRUE))
  ) %>% 
  filter(pt_id %in% e2_text_long$pt_id)


```

## Text Visualize

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

```{r text demographics}

# participant age
ggplot(pt_e2_text, 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 [Text]"
  )

# participant education
ggplot(pt_e2_text, 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 [Text]"
  )

# literacy scores 
ggplot(pt_e2_text, 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 [Text]",
    fill = "Graphical Literacy Score"
  )

# preference for each Mode
prefs_e2_text <- e2_text %>%
  pivot_longer(
    cols = starts_with("pref_"),
    names_to = "source_mode",
    values_to = "rank"
  ) %>%
  filter(!is.na(rank)) %>% 
  select(pt_id, precision_level, mode, source_mode, rank)
ggplot(prefs_e2_text, aes(x = rank, fill = as.factor(source_mode)))+
  geom_bar(position = "dodge")+
  scale_fill_manual(values = c("pref_audio" = "#C1A6CE", "pref_text" = "#AAB2D0", "pref_visual" = "#9DC3C6") ) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nRank",
    title = "Mode Rankings [Text]\n",
    fill = "Mode"
  )

```

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 text decision}

# decision type by precision
decision_precision_e2_text <- e2_text_main %>%
  filter(trial != 0) %>% 
  group_by(precision_level, rational_code) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(precision_level) %>%
  mutate(total = sum(n), 
         pct = n / total)
ggplot(decision_precision_e2_text, aes(x = rational_code, y = pct, fill = rational_code))+
  geom_col(position = "dodge")+
  scale_fill_manual(values = rational_colors) + # Use manually defined colors
  facet_grid(.~ precision_level)+
  scale_x_discrete(labels = c("Cons.", "Ration.", "Risky"))+
  labs(x = "\nPrecision Level",
       fill = "Decision Type",
       y = "Proportion of Decisions\n",
       title = "Decision Types by Precision Level")+
  ylim(c(0, 1))+
  theme(legend.position = "none") # Remove the legend


# confidence in choice by forecast variant
ggplot(e2_text_main %>% filter(trial != 0), aes(x = confidence, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =5)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_x_continuous(limits = c(50,100))+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nDecision Confidence",
    title = "Confidence in Decisions by Forecast",
    fill = "Precision Level"
  )+
  theme(legend.position = "none",
        strip.text.y = element_text(angle = 0, hjust = 0)) # Remove the legend


```

```{r text trust ratings}


# ratings for clarity, accuracy, usefulness (and calculated overall trust)
trust_e2_text <- e2_text %>%
  pivot_longer(cols = c(clear, accurate, useful, trust_average),
               names_to = "rating_type",
               values_to = "rating_value") %>% 
  select(rating_type, rating_value, precision_level)%>% 
  mutate(
    precision_level = factor(precision_level, levels = c("Most Concrete", "Somewhat Concrete", "Mixed", "Somewhat Fuzzy", "Most Fuzzy"))
  )


# clarity by forecast variant
ggplot(trust_e2_text %>% filter(rating_type == "clear"), aes( x = rating_value, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nForecast Attribute",
    title = "Distribution of Forecast Clarity Ratings",
    fill = "Precision Level",
    color = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend

# accuracy by forecast variant
ggplot(trust_e2_text %>% filter(rating_type == "accurate"), aes( x = rating_value, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nForecast Attribute",
    title = "Distribution of Forecast Accuracy Ratings",
    fill = "Precision Level",
    color = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend

# useful by forecast variant
ggplot(trust_e2_text %>% filter(rating_type == "useful"), aes( x = rating_value, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nForecast Attribute",
    title = "Distribution of Forecast Usefulness Ratings",
    fill = "Precision Level",
    color = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend

# average trust in forecast variant (mean values of clarity, accurate, and useful)
ggplot(e2_text_main, aes( x = trust_average, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nForecast Attribute",
    title = "Distribution of Forecast Trust Ratings",
    fill = "Precision Level",
    color = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend


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

```{r text accuracy in estimations}

# numeric temp estimate error by forecast variant
ggplot(e2_text_main, aes(x = numeric_estimate_diff, fill = precision_level))+
  geom_bar()+
  facet_grid(precision_level ~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  xlim(-10, 10)+
  labs(
    y = "Precision Level\n",
    x = "\nDifference between Actual\nand Estimated Numeric Temperature",
    title = "Numeric Temperature Estimation Accuracy by Forecast",
    fill = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend

# difference between pct_likelihood and actual likelihood by forecast variant
ggplot(e2_text_main, aes(x = pct_estimate_diff, fill = precision_level))+
  geom_bar()+
  facet_grid(precision_level ~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  labs(
    y = "Precision Level\n",
    x = "\nDifference between Actual\n and Estimated Freeze Likelihood (%)",
    title = "Freeze Likelihood Estimation Accuracy by Forecast",
    fill = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend


```


## Text Analyze

The following sections address each hypothesis in turn.

If there's an indication that the level of precision may have had an effect, we run additional models with varying referents, as our main factor of interest is categorical.

### Text Crossover Temperature

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.


```{r text salt choice model}

#### stepwise models ####

e2_text_main$precision_level <-  factor(e2_text_main$precision_level, levels = c('Most Concrete', 'Somewhat Concrete', 'Mixed','Somewhat Fuzzy', 'Most Fuzzy'))

# baseline - does salt likelihood change with temperature? it should. 
salt0_text <- glmer(salt_code ~ crossover_distance + (1|pt_id) , 
                    data = e2_text_main,
                    family = binomial("logit"),  control = glmerControl(tol = 0.0001))

# add in precision - is salt likelihood different for precision level
salt1_text <- glmer(salt_code ~ crossover_distance + precision_level + (1|pt_id) , 
                    data = e2_text_main,
                    family = binomial("logit"),  control = glmerControl(tol = 0.0001)) 

# most complex model - not included in comparisons but was pre-registered to run
salt2_text <- glmer(salt_code ~ crossover_distance + precision_level + ageFactor + educationFactor + pref1 + (1|pt_id) , 
                    data = e2_text_main,
                    nAGQ = 0,
                    family = binomial("logit"))


#### select and run model(s) ####

# find model with best fit
anova(salt0_text, salt1_text)

summary(salt0_text)
r.squaredGLMM(salt0_text)

# calculate the coefs and errors for model interpretation
coef_salt_text <- fixef(salt1_text)

# For each combination of levels of precision_level, solve for temperature when log odds = 0
for(precision in precision_levels) {
  
  # Get the coefficients relevant for this combination of levels
  intercept <- coef_salt_text[1]
  condition <- sum(c(coef_salt_text[1],
                     coef_salt_text[paste0("precision_level", precision, sep="")]),
                   na.rm= TRUE)
  temp_coef <- coef_salt_text[2]
  
  # Solve for temperature: logit(0.5) = 0 = intercept + temp_coef * temperature
  temperature <- ((-1 * intercept) - condition)  / temp_coef
  
  cat("For precision =", precision, ", the crossover temperature is", temperature, "from optimal.\n")
  
}

```


### Text 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 text decision rationality}

precision_table_text <- table(e2_text_main$precision_level, 
                              e2_text_main$rational_code)
precision_table_text
text_chisq <- chisq.test(precision_table_text)
text_chisq

```


### Text Confidence

Next, we want to look at how confidence participants were in the decisions they made. Here, we look to see if confidence reliably varied between different levels of fuzziness. 

```{r text confidence model}

#### stepwise models ####

e2_text_main$precision_level <-  factor(e2_text_main$precision_level, levels = c('Most Concrete', 'Somewhat Concrete', 'Mixed','Somewhat Fuzzy', 'Most Fuzzy'))

# baseline 
confidence0_text <- lmer(confidence ~ crossover_distance + (1|pt_id) , 
                         data = e2_text_main) 

# add in rationality
confidence1_text <- lmer(confidence ~ crossover_distance + is_rational + (1|pt_id) , 
                         data = e2_text_main) 

# add in precision 
confidence2_text <- lmer(confidence ~ crossover_distance + is_rational + precision_level + (1|pt_id) , 
                         data = e2_text_main) 

# most complex model - not included in comparisons but was pre-registered to run
confidence3_text <- lmer(confidence ~ crossover_distance + precision_level + is_rational + ageFactor + educationFactor + pref1 + (1|pt_id) , 
                         data = e2_text_main)

#### select and run model(s) ####

# find model with best fit
anova(confidence0_text, confidence1_text, confidence2_text)

# conduct model
summary(confidence1_text)
r.squaredGLMM(confidence1_text)

confidence_coef_text <- fixef(confidence2_text)

for(precision in precision_levels) {
  
  # Get the coefficients relevant for this combination of levels
  confidence <- sum(c(confidence_coef_text[1],
                      confidence_coef_text[paste0("precision_level", precision, sep="")]),
                    na.rm= TRUE)
  
  cat("For", precision, "the confidence is", confidence, "\n")
}

```


### Text 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 text select trust model}

#### stepwise models ####
pt_e2_text$precision_level <-  factor(pt_e2_text$precision_level, levels = c('Most Concrete', 'Somewhat Concrete', 'Mixed','Somewhat Fuzzy', 'Most Fuzzy'))

# baseline - does trust change with how rational participants were?
trust0_text <- lm(trust_average ~ rationality, 
                  data = pt_e2_text)
# add in mode - is trust different for the different modalities
trust1_text <- lm(trust_average ~ rationality + precision_level, 
                  data = pt_e2_text)
pt_e2_text$precision_level <-  factor(pt_e2_text$precision_level, levels = c('Somewhat Concrete', 'Most Concrete', 'Mixed','Somewhat Fuzzy', 'Most Fuzzy'))
trust1_text_2 <- lm(trust_average ~ rationality + precision_level, 
                    data = pt_e2_text)
pt_e2_text$precision_level <-  factor(pt_e2_text$precision_level, levels = c('Mixed', 'Somewhat Concrete', 'Most Concrete', 'Somewhat Fuzzy', 'Most Fuzzy'))
trust1_text_3 <- lm(trust_average ~ rationality + precision_level, 
                    data = pt_e2_text)
pt_e2_text$precision_level <-  factor(pt_e2_text$precision_level, levels = c('Somewhat Fuzzy', 'Mixed', 'Somewhat Concrete', 'Most Concrete', 'Most Fuzzy'))
trust1_text_4 <- lm(trust_average ~ rationality + precision_level, 
                    data = pt_e2_text)
pt_e2_text$precision_level <-  factor(pt_e2_text$precision_level, levels = c( 'Most Fuzzy', 'Somewhat Fuzzy', 'Mixed', 'Somewhat Concrete', 'Most Concrete'))
trust1_text_5 <- lm(trust_average ~ rationality + precision_level, 
                    data = pt_e2_text)
pt_e2_text$precision_level <-  factor(pt_e2_text$precision_level, levels = c('Most Concrete', 'Somewhat Concrete', 'Mixed','Somewhat Fuzzy', 'Most Fuzzy'))

# most complex model - not included in comparisons but was pre-registered to run
trust2_text <- lm(trust_average ~ rationality + precision_level +  ageFactor + educationFactor + pref1, 
                  data = pt_e2_text)

#### select and run model(s)

# find model with best fit
anova(trust0_text, trust1_text)

# compare different referents
summary(trust1_text)
summary(trust1_text_2)
summary(trust1_text_3)
summary(trust1_text_4)
summary(trust1_text_5)

# conduct model
summary(trust1_text)

# calculate coefs for optimal model
trust_coef_text <- coef(trust1_text)

for(precision in precision_levels) {
  
  # Get the coefficients relevant for this combination of levels
  trust <- sum(c(trust_coef_text[1],
                 trust_coef_text[paste0("precision_level", precision, sep="")]),
               na.rm= TRUE)
  
  cat("For", precision, "the trust, on average, is", trust, "out of 10 \n")
  
}

```


# Visual

Finally, we analyze the visualization survey. Here, we compared variations of visualization as well as visualization techniques, such as a gradient fill.

## Visual 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 visual load data}


#### read in raw e2_visual data ####
e2_visual_raw <- e2_visual_raw <- read.csv("e2_visual_qualtrics.csv") 

#### select and calculate columns ####
# column definitions included here
e2_visual <- e2_visual_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
  ) %>% 
  # calculate columns here
  mutate(
    # create anonymous participant id
    pt_id = 1:nrow(e2_visual_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),
    # create column for mode
    mode = "Visualization",
    # add in concrete and fuzzy information
    precision_level =  precision_vector[Forecast],
    # create factor level columns for plotting and modeling
    educationFactor = factor(education, levels = educationLevels),
    ageFactor = factor(age, levels = ageLevels), 
    precision_level = factor(precision_level, levels = precision_levels),
    # rename the preference columns to make sense
    pref_visual = info_pref_1,
    pref_text = info_pref_2,
    pref_audio = info_pref_3,
    pref1 = case_when(
      pref_visual == 1 ~ 'visual',
      pref_text == 1 ~ 'text',
      pref_audio == 1 ~ 'audio',
      TRUE ~ 'unknown',
    ), # make single column for preference - what was the most preferred? 
  ) %>% 
  # 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 exclusion ####
# exclusion based on tail percentage estimation, only easily accessible after the data has been transformed
e2_visual_long <- e2_visual %>% 
  # 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 = round(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")),
    EL_difference = EL - rational_EL,
    # crossover_distance = avg_temp - optimal_temperature, # difference for each trial from optimal temperature
    PlayCount = NA
  ) 

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


#### select columns for analysis ####
e2_visual_main <- e2_visual_long %>% 
  select(
    pt_id,
    mode, # 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
    precision_level = factor(precision_level, levels = c("Most Concrete", "Somewhat Concrete", "Mixed", "Somewhat Fuzzy", "Most Fuzzy"))
  )

#### create data frame for participant information ####
pt_e2_visual <- e2_visual %>% 
  select(
    pt_id,
    Budget,
    ageFactor, # pt age range
    educationFactor, # pt education level
    graph_lit, # graph_lit score
    mode, # 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(
    e2_visual_long %>% group_by(pt_id) %>%  summarise(rationality = mean(is_rational, na.rm = TRUE))
  ) %>% 
  filter(pt_id %in% e2_visual_long$pt_id)


```

## Visual Visualize

Now, we visualize preliminary results & information. This helps us get a sense of the data prior to analysis.

```{r visual demographics}


# participant age
ggplot(pt_e2_visual, 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 [Visual]"
  )

# participant education
ggplot(pt_e2_visual, 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 [Visual]"
  )

# literacy scores 
ggplot(pt_e2_visual, 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 [Visual]",
    fill = "Graphical Literacy Score"
  )


# preference for each Mode
prefs_e2_visual <- e2_visual %>%
  pivot_longer(
    cols = starts_with("pref_"),
    names_to = "source_mode",
    values_to = "rank"
  ) %>%
  filter(!is.na(rank)) %>% 
  select(pt_id, precision_level, mode, source_mode, rank)
ggplot(prefs_e2_visual, aes(x = rank, fill = as.factor(source_mode)))+
  geom_bar(position = "dodge")+
  scale_fill_manual(values = c("pref_audio" = "#C1A6CE", "pref_text" = "#AAB2D0", "pref_visual" = "#9DC3C6") ) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nRank",
    title = "Mode Rankings [Visual]\n",
    fill = "Mode"
  )

```

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 visual decision}

# decision type by precision
decision_precision_e2_visual <- e2_visual_main %>%
  filter(trial != 0) %>% 
  group_by(precision_level, rational_code) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(precision_level) %>%
  mutate(total = sum(n), 
         pct = n / total)
ggplot(decision_precision_e2_visual, aes(x = precision_level, y = pct, fill = rational_code))+
  geom_col(position = "dodge")+
  scale_fill_manual(values = rational_colors) + # Use manually defined colors
  labs(x = "\nPrecision Level",
       fill = "Decision Type",
       y = "Proportion of Decisions\n",
       title = "Decision Types by Precision Level")+
  ylim(c(0, 1))+
  theme(legend.position = "none") # Remove the legend


# confidence in choice by forecast variant
ggplot(e2_visual_main %>% filter(trial != 0), aes(x = confidence, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =5)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_x_continuous(limits = c(50,100))+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nDecision Confidence",
    title = "Confidence in Decisions by Forecast",
    fill = "Precision Level"
  )+
  theme(legend.position = "none", 
        strip.text.y = element_text(angle = 0, hjust= 0)) # Remove the legend


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

```{r visual trust ratings}

# ratings for clarity, accuracy, usefulness (and calculated overall trust)
trust_e2_visual <- e2_visual %>%
  pivot_longer(cols = c(clear, accurate, useful, trust_average),
               names_to = "rating_type",
               values_to = "rating_value") %>% 
  select(rating_type, rating_value, precision_level) %>% 
  mutate(
    precision_level = factor(precision_level, levels = c("Most Concrete", "Somewhat Concrete", "Mixed", "Somewhat Fuzzy", "Most Fuzzy"))
  )

# clarity by forecast variant
ggplot(trust_e2_visual %>% filter(rating_type == "clear"), aes( x = rating_value, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nForecast Attribute",
    title = "Distribution of Forecast Clarity Ratings",
    fill = "Precision Level",
    color = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend

# accuracy by forecast variant
ggplot(trust_e2_visual %>% filter(rating_type == "accurate"), aes( x = rating_value, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nForecast Attribute",
    title = "Distribution of Forecast Accuracy Ratings",
    fill = "Precision Level",
    color = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend

# useful by forecast variant
ggplot(trust_e2_visual %>% filter(rating_type == "useful"), aes( x = rating_value, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nForecast Attribute",
    title = "Distribution of Forecast Usefulness Ratings",
    fill = "Precision Level",
    color = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend

# average trust in forecast variant (mean values of clarity, accurate, and useful)
ggplot(e2_visual_main, aes(x = trust_average, fill = precision_level, color = precision_level))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  scale_color_manual(values = precision_line_colors) + # Use manually defined colors
  labs(
    y = "Count\n",
    x = "\nForecast Attribute",
    title = "Distribution of Forecast Trust Ratings",
    fill = "Precision Level",
    color = "Precision Level"
  )+
  theme(legend.position = "none",
        strip.text.y = element_text(angle = 0, hjust= 0)) # Remove the legend


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

```{r visual accuracy in estimations}

# numeric temp estimate error by forecast variant
ggplot(e2_visual_main %>%  filter(trial !=0), aes(x = numeric_estimate_diff, fill = precision_level))+
  geom_bar()+
  facet_grid(precision_level ~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  xlim(-10, 10)+
  labs(
    y = "Precision Level\n",
    x = "\nDifference between Actual\nand Estimated Numeric Temperature",
    title = "Numeric Temperature Estimation Accuracy by Forecast",
    subtitle = "Three outlier responses not shown",
    fill = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend

# difference between pct_likelihood and actual likelihood by forecast variant
ggplot(e2_visual_main %>%  filter(trial !=0), aes(x = pct_estimate_diff, fill = precision_level))+
  geom_bar()+
  facet_grid(precision_level ~.)+
  scale_fill_manual(values = precision_fill_colors) + # Use manually defined colors
  labs(
    y = "Precision Level\n",
    x = "\nDifference between Actual\n and Estimated Freeze Likelihood (%)",
    title = "Freeze Likelihood Estimation Accuracy by Forecast",
    fill = "Precision Level"
  )+
  theme(legend.position = "none") # Remove the legend

```


## Visual Analyze

The following sections address each hypothesis in turn.

If there's an indication that the level of precision may have had an effect, we run additional models with varying referents, as our main factor of interest is categorical.

### Visual Crossover Temperature

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.


```{r visual salt choice model}

#### stepwise models ####

e2_visual_main$precision_level <-  factor(e2_visual_main$precision_level, levels = c('Most Concrete', 'Somewhat Concrete', 'Mixed','Somewhat Fuzzy', 'Most Fuzzy'))

# baseline 
salt0_visual <- glmer(salt_code ~ crossover_distance  + (1|pt_id), 
                      data = e2_visual_main,
                      family = binomial("logit")) 
# add in precision 
salt1_visual <- glmer(salt_code ~ crossover_distance  + precision_level + (1|pt_id), 
                      data = e2_visual_main,
                      family = binomial("logit")) 
e2_visual_main$precision_level <-  factor(e2_visual_main$precision_level, levels = c('Somewhat Concrete', 'Most Concrete', 'Mixed','Somewhat Fuzzy', 'Most Fuzzy'))
salt1_visual_2 <-    glmer(salt_code ~ crossover_distance  + precision_level + (1|pt_id),
                           data = e2_visual_main,
                           family = binomial("logit"))
e2_visual_main$precision_level <-  factor(e2_visual_main$precision_level, levels = c('Mixed', 'Somewhat Concrete', 'Most Concrete', 'Somewhat Fuzzy', 'Most Fuzzy'))
salt1_visual_3 <-    glmer(salt_code ~ crossover_distance  + precision_level + (1|pt_id),
                           data = e2_visual_main,
                           family = binomial("logit"))
e2_visual_main$precision_level <-  factor(e2_visual_main$precision_level, levels = c('Somewhat Fuzzy', 'Mixed', 'Somewhat Concrete', 'Most Concrete', 'Most Fuzzy'))
salt1_visual_4 <-    glmer(salt_code ~ crossover_distance  + precision_level + (1|pt_id),
                           data = e2_visual_main,
                           family = binomial("logit"))
e2_visual_main$precision_level <-  factor(e2_visual_main$precision_level, levels = c( 'Most Fuzzy', 'Somewhat Fuzzy', 'Mixed', 'Somewhat Concrete', 'Most Concrete'))
salt1_visual_5 <-    glmer(salt_code ~ crossover_distance  + precision_level + (1|pt_id),
                           data = e2_visual_main,
                           family = binomial("logit"))
e2_visual_main$precision_level <-  factor(e2_visual_main$precision_level, levels = c('Most Concrete', 'Somewhat Concrete', 'Mixed','Somewhat Fuzzy', 'Most Fuzzy'))

# most complex model - not included in comparisons but was pre-registered to run
salt2_visual <- glmer(salt_code ~ crossover_distance  + precision_level + graph_lit + ageFactor + educationFactor + pref1 + (1|pt_id), 
                      data = e2_visual_main,
                      nAGQ = 0,
                      family = binomial("logit"))

#### select and run model(s) ####

# find model with best fit
anova(salt0_visual, salt1_visual)

summary(salt1_visual)
summary(salt1_visual_2)
summary(salt1_visual_3)
summary(salt1_visual_4)
summary(salt1_visual_5)
# summarize best model
r.squaredGLMM(salt1_visual)

# get model coefs and errors
coef_salt_visual <- fixef(salt1_visual)

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

```

### Visual 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 visual decision rationality}

precision_table_visual <- table(e2_visual_main$precision_level, 
                                e2_visual_main$rational_code)
precision_table_visual
visual_chisq <- chisq.test(precision_table_visual)
visual_chisq
residuals_visual <- visual_chisq$stdres
print(residuals_visual)

```

### Visual Confidence

Next, we want to look at how confidence participants were in the decisions they made. Here, we look to see if confidence reliably varied between different levels of fuzziness. 

```{r visual confidence model}

#### stepwise models ####

e2_visual_main$precision_level <-  factor(e2_visual_main$precision_level, levels = c('Most Concrete', 'Somewhat Concrete', 'Mixed','Somewhat Fuzzy', 'Most Fuzzy'))

# baseline
confidence0_visual <- lmer(confidence ~ crossover_distance + (1|pt_id) , 
                           data = e2_visual_main) 
# add in rationality
confidence1_visual <- lmer(confidence ~ crossover_distance + is_rational + (1|pt_id) , 
                           data = e2_visual_main) 

# add in precision
confidence2_visual <- lmer(confidence ~ crossover_distance + is_rational + precision_level + (1|pt_id) , 
                           data = e2_visual_main) 

# most complex model - not included in comparisons but was pre-registered to run
confidence3_visual <- lmer(confidence ~ crossover_distance + precision_level + is_rational + graph_lit + ageFactor + educationFactor + pref1 + (1|pt_id) , 
                           data = e2_visual_main) 


#### select and run model(s) ####

# find model with best fit
anova(confidence0_visual, confidence1_visual, confidence2_visual)

# summarize model
summary(confidence1_visual)

# calculate the coefs for model interpretation 
confidence_coef_visual <- fixef(confidence2_visual)

for(precision in precision_levels) {
  
  # Get the coefficients relevant for this combination of levels
  confidence <- sum(c(confidence_coef_visual[1],
                      confidence_coef_visual[paste0("precision_level", precision, sep="")]),
                    na.rm= TRUE)
  
  cat("For", precision, "the confidence is", confidence, "\n")
}


```


### Visual 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 visual trust model}

#### stepwise models ####

pt_e2_visual$precision_level <-  factor(pt_e2_visual$precision_level, levels = c('Most Concrete', 'Somewhat Concrete', 'Mixed','Somewhat Fuzzy', 'Most Fuzzy'))

# baseline 
trust0_visual <- lm(trust_average ~ rationality, 
                    data = pt_e2_visual)
# add in precision
trust1_visual <- lm(trust_average ~ rationality + precision_level, 
                    data = pt_e2_visual)

# most complex model - not included in comparisons but was pre-registered to run
trust2_visual <- lm(trust_average ~ rationality + precision_level + graph_lit +  ageFactor + educationFactor + pref1, 
                    data = pt_e2_visual)

#### select and run model(s) ####

# find model with best fit
anova(trust0_visual, trust1_visual)

summary(trust0_visual)

```



# Exploratory Analysis

For the exploratory analyses, we want to work with one large data frame containing all modes of information


```{r combine modes}

e2_combined <- rbind(e2_visual_main, e2_text_main, e2_speech_main)
e2_speech <- e2_speech %>% select(-contains("PlayCount"))
pt_e2_combined <- rbind(pt_e2_visual, pt_e2_text, pt_e2_speech)

pt_e2_combined %>% 
  group_by(pref1) %>% 
  summarise(
    n = n()
  )

# participant age
ggplot(pt_e2_combined, aes(x = ageFactor))+
  geom_bar(fill = "#aaaaaa")+
  scale_x_discrete(limits = rev(ageLevels), labels = rev(c("18-24", "25-34", "35-44", "45-54", "55-64", "65-74", "75-84", ">85", "No answer")))+
  labs(
    x = "Age Range",
    y = "Count",
    title = "Age Distribution [E2]"
  )+
  coord_flip()

# participant education
ggplot(pt_e2_combined, aes(x = educationFactor))+
  geom_bar(fill = "#aaaaaa")+
  scale_x_discrete(limits = rev(educationLevels), labels = rev(c("Less than high school", "High school graduate", "Some college", "2 year degree", "4 year degree", "Professional degree", "Doctorate", "No answer")))+
  labs(
    x = "Education Level",
    y = "Count",
    title = "Education Distribution [E2]"
  )+
  coord_flip()

```

```{r overall visualizations}

# Decision Rationality
decision_interaction_e2 <- e2_combined %>%
  filter(!is.na(rational_code)) %>% 
  group_by(mode, precision_level, rational_code) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(mode, precision_level) %>%
  mutate(total = sum(n), 
         pct = n / total,
         mode = factor(mode, levels = rev(c('Speech', 'Text', 'Visualization'))))

ggplot(decision_interaction_e2, aes(x = rational_code, y = pct, fill = rational_code))+
  geom_col(position = "dodge")+
  facet_grid( fct_rev(factor(mode)) ~ precision_level)+
  scale_fill_manual(values = rational_colors) + 
  labs(x = "\nDecision Type",
       fill = "Decision Type",
       y = "Proportion of Decisions\n",
       title = "Decision Types by Mode and Precision Level"
  )+
  theme(
    legend.position = "none")
ggsave("e2_rationality.pdf", height = 5, width = 12)


# Decision Confidence
ggplot(e2_combined %>% filter(trial != 0), aes(x = confidence, fill = mode, color = mode))+
  geom_histogram(binwidth =5)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~mode)+
  scale_x_continuous(limits = c(50,100))+
  scale_fill_manual(values = mode_fill_colors) + 
  scale_color_manual(values = mode_line_colors) +  
  labs(
    y = "Count\n",
    x = "\nDecision Confidence",
    title = "Confidence in Decisions by Precision and Mode"
  )+
  theme_minimal()+
  theme(legend.position = "none", 
        strip.text.y = element_text(angle = 0, hjust = 0))
ggsave("e2_confidence.pdf", height = 5, width = 10)

# Trust in Forecast

# average trust in forecast Mode (mean values of clarity, accurate, and useful)
ggplot(e2_combined, aes(x = trust_average, fill = mode, color = mode))+
  geom_histogram(binwidth =1)+
  stat_slabinterval(size = 1)+
  facet_grid(precision_level~mode)+
  scale_fill_manual(values = mode_fill_colors) + 
  scale_color_manual(values = mode_line_colors) +  
  labs(
    y = "Forecast Mode\n",
    x = "\nAverage Trust in Forecast",
    title = "Distribution of Forecast Trust Ratings"
  )+
  theme_minimal()+
  theme(legend.position = "none",
        strip.text.y = element_text(angle = 0, hjust = 0))
ggsave("e2_trust.pdf", height = 5, width = 10)


```


## Experiment 1 Replication

First, we want to re-run the selected models from the first experiment.

```{r filter for e1 conditions}

e2_replication <-  e2_combined %>% 
  filter(
    (mode == "Speech" & precision_level %in% c("Most Concrete", "Somewhat Fuzzy")) |
      (mode == "Text" & precision_level %in% c("Most Concrete", "Mixed")) |
      (mode == "Visualization" & precision_level %in% c("Most Concrete", "Somewhat Concrete"))
  ) %>% 
  mutate(
    precision_level = case_when(
      (mode == "Speech" & precision_level == "Most Concrete") ~ "Concrete",
      (mode == "Speech" & precision_level == "Somewhat Fuzzy") ~ "Fuzzy",
      (mode == "Text" & precision_level == "Most Concrete") ~ "Concrete",
      (mode == "Text" & precision_level == "Mixed") ~ "Fuzzy",
      (mode == "Visualization" & precision_level == "Most Concrete") ~ "Concrete",
      (mode == "Visualization" & precision_level == "Somewhat Concrete") ~ "Fuzzy",
      TRUE ~ precision_level
    )
  )
pt_e2_replication <- pt_e2_combined %>% 
  filter(
    (mode == "Speech" & precision_level %in% c("Most Concrete", "Somewhat Fuzzy")) |
      (mode == "Text" & precision_level %in% c("Most Concrete", "Mixed")) |
      (mode == "Visualization" & precision_level %in% c("Most Concrete", "Somewhat Concrete"))
  ) %>% 
  mutate(
    precision_level = case_when(
      (mode == "Speech" & precision_level == "Most Concrete") ~ "Concrete",
      (mode == "Speech" & precision_level == "Somewhat Fuzzy") ~ "Fuzzy",
      (mode == "Text" & precision_level == "Most Concrete") ~ "Concrete",
      (mode == "Text" & precision_level == "Mixed") ~ "Fuzzy",
      (mode == "Visualization" & precision_level == "Most Concrete") ~ "Concrete",
      (mode == "Visualization" & precision_level == "Somewhat Concrete") ~ "Fuzzy",
      TRUE ~ precision_level 
    )
  )

```

```{r e1 replication models}

salt_e1 <- glmer(salt_code ~ crossover_distance + mode + (1|pt_id) , 
                 data = e2_replication,
                 nAGQ = 0,
                 family = binomial("logit"))
summary(salt_e1)

stargazer::stargazer(salt_e1, "latex")


precision_table_replication <- table(e2_replication$precision_level, 
                                     e2_replication$rational_code)
precision_table_replication
precision_replication_chisq <- chisq.test(precision_table_replication)
precision_replication_chisq
precision_replication_chisq$stdres

modality_table_replication <- table(e2_replication$mode, 
                                    e2_replication$rational_code)
modality_table_replication
modality_replication_chisq <- chisq.test(modality_table_replication)
modality_replication_chisq

e2_replication <- e2_replication %>% 
  mutate(
    mode_confidence = factor(mode, levels = c("Visualization", "Text", "Speech"))
  )
confidence_e1 <- lmer(confidence ~ crossover_distance + is_rational + mode_confidence + (1|pt_id) + (1|trial), 
                      data = e2_replication) 
summary(confidence_e1)
confidence_e1.1 <- lmer(confidence ~ crossover_distance + is_rational + mode + (1|pt_id)  + (1|trial), 
                        data = e2_replication) 
summary(confidence_e1.1)

trust_e1 <- lm(trust_average ~ rationality + mode, 
               data = pt_e2_combined)
summary(trust_e1)

```

## Large Overall Models

Then, we expand out models to account for the additional precision conditions studied in E2. While the E1 model for trust did not incorporate precision level, we want to include it in this analysis to account for the additional precision levels. To ensure that this is a valid course of action, we conduct an ANOVA between the two models to ensure that adding the precision levels improves prediction over the original model. 

```{r large combined models}

salt_combined <- glmer(salt_code ~ crossover_distance + mode * precision_level + ageFactor + educationFactor + pref1 + (1|pt_id) , 
                       data = e2_replication,
                       nAGQ = 0,
                       family = binomial("logit"))
summary(salt_combined)


confidence_combined <- lmer(confidence ~ crossover_distance + mode * precision_level + is_rational + ageFactor + educationFactor + (1|pt_id) , 
                            data = e2_replication)
summary(confidence_combined)


trust_combined0 <- lm(trust_average ~ rationality + mode, 
                      data = pt_e2_combined)
trust_combined <- lm(trust_average ~ rationality + mode + precision_level, 
                     data = pt_e2_combined)
anova(trust_combined0, trust_combined)

summary(trust_combined)


```


