---
title: "Is depression associated with reduced optimistic belief updating?"
output: html_notebook
---

```{r}
library(tidyverse)#General cleaning
library(expss) #Descriptive Tables
expss_output_rnotebook()
expss_digits(digits = 2)
library(lmerTest) #MLMs
library(jtools) #Summary Tables for Regression
library(ggeffects) #Regression Graphs
library(see) #Multiple panel figures
library(effectsize) #Standardising parameters
library(gridExtra)
library(emmeans)#Group contrasts
```

```{r}
#Importing event level data (from p.s that passed all attention checks)
load("task_event_level_qs.RData")
```

# Sample Characteristics

```{r}
sample_characteristics <- task_event_level_qs %>%
  select(subject, Group, age, gender, ethnicity, employment_aggregated, education_aggregated, relationship, living, screening_PHQ_tot, PHQ_tot, BDI_tot, GAD_tot, BFNE_tot, LOT_tot, PANAS_pos_change, PANAS_neg_change, depression, psych_ther_curr, antid_curr, depress_family_hist) %>%
  group_by(subject) %>%
  distinct() %>% #Currently at trial level so multiple rows per participant - taking one row per participant for characteristics
  ungroup()


#Sample Characteristics 
sample_characteristics %>%
  tab_cols(Group) %>%
  tab_cells(age) %>%
  tab_stat_mean_sd_n() %>%
  tab_cells(gender, ethnicity, employment_aggregated, education_aggregated, relationship, living) %>%
  tab_stat_cases(label = "N", total_row_position = "none") %>%
  tab_stat_cpct(label = "%", total_row_position = "none") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Demographics by Group")


#Clinical Characteristics 
sample_characteristics %>%
  tab_cols(Group) %>%  
  tab_cells(screening_PHQ_tot, PHQ_tot, BDI_tot, GAD_tot, BFNE_tot, LOT_tot, PANAS_pos_change, PANAS_neg_change) %>%
  tab_stat_mean_sd_n() %>%
  tab_cells(depression, psych_ther_curr, antid_curr, depress_family_hist) %>%
  tab_stat_cases(label = "N", total_row_position = "none") %>%
  tab_stat_cpct(label = "%", total_row_position = "none") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Clinical Characteristics by Group")
```


# Hypothesis 1
Healthy controls will show a greater change in beliefs following desirable versus undesirable information for negative life events, whereas individuals with depression will show a smaller difference in change in beliefs following desirable versus undesirable information for negative life events. 

A mixed-effects linear regression model will be created using only trials for negative life events. Absolute mean update scores will be entered as the outcome, with group (healthy versus depressed), desirability (desirable versus undesirable), and an interaction term between group and desirability entered as predictors. Subject will be entered as a random effect to account for the repeated measures design. 

```{r}
#Aggregating scores from event level data
negative_eligible_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group, desirability) %>%
  summarise(update = mean(update, na.rm = TRUE)) %>%
  ungroup()

#Descriptive Statistics
negative_eligible_aggregate_data %>%
  tab_cols(Group, desirability, Group %nest% desirability) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores for Negative Life Events by Group and Desirability")

#Bar Graph by Group
negative_eligible_aggregate_data %>%
  group_by(Group) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = Group)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Group") +
  ggtitle("Update Scores by Group (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

negative_eligible_aggregate_data %>%
  tab_cols(Group, desirability, Group %nest% desirability) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores for Negative Life Events by Group and Desirability")#Bar Graph by Desirability
negative_eligible_aggregate_data %>%
  group_by(desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = desirability, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Desirability") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Group and Desirability
negative_eligible_aggregate_data %>%
  group_by(Group, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Group & Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

negative_eligible_aggregate_data %>%
  group_by(Group, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp1_model <- lmer(update ~ Group*desirability + (1|subject), data = negative_eligible_aggregate_data)
summ(hyp1_model, confint = TRUE, digits = 3)
hyp1_model_effects <- ggpredict(hyp1_model, terms = c("Group", "desirability"))
plot(hyp1_model_effects) +
  labs(
    title = "Predicted Values of Update Scores for Negative Life Events",
    y = "Update Scores",
    x = "Group"
  )


standardize_parameters(hyp1_model, method = "refit", ci = 0.95, include_response = TRUE)

#Comparing Desirability by group
contrasts_hyp1_model_desirability_by_group <- emmeans(hyp1_model, pairwise ~ desirability | Group, adjust = "tukey")
contrasts_hyp1_model_desirability_by_group
confint(contrasts_hyp1_model_desirability_by_group)

```

# Hypothesis 2
Hypothesis 2: The effect of optimism bias will be consistent across both positive and negative life events. Healthy controls will show a greater change in beliefs following desirable versus undesirable information irrespective of the valence of the life event. Individuals with depression will show smaller differences in change in beliefs following desirable versus undesirable information for both positive and negative life events.

A mixed-effects linear regression model will be conducted using trials for both positive and negative life events. Absolute mean update scores will be entered as the outcome, with group (healthy versus depressed), desirability (desirable versus undesirable) and valence (positive versus negative) entered as predictors. Additionally, interaction terms between group and desirability, group and valence, desirability and valence, and group, desirability and valence will be entered as additional predictors. Subject will be entered as a random effect to account for the repeated measures design. 

```{r}
#Aggregating Event Level Data
positive_negative_eligible_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability) %>%
  summarise(update = mean(update, na.rm = TRUE)) %>%
  ungroup()

#Descriptives
positive_negative_eligible_aggregate_data %>%
  tab_cols(Group, valence, desirability) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores Group, Valence and Desirability (Main Effects)")


positive_negative_eligible_aggregate_data %>%
  tab_cols(Group %nest% valence, Group %nest% desirability, valence %nest% desirability) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores Group, Valence and Desirability (2 Way Interaction Effects)")

positive_negative_eligible_aggregate_data %>%
  tab_cols(Group %nest% valence %nest% desirability) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores Group, Valence and Desirability (3 way Interaction Effects)")

## Switching columns for table 1
positive_negative_eligible_aggregate_data %>%
  tab_cols(valence %nest% Group %nest% desirability) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores Group, Valence and Desirability (3 way Interaction Effects)")


#Bar Graph by Group
positive_negative_eligible_aggregate_data %>%
  group_by(Group) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = Group)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Group") +
  ggtitle("Update Scores by Group (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Valence
positive_negative_eligible_aggregate_data %>%
  group_by(valence) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = valence, y = mean, fill = valence)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Valence") + 
  labs(fill = "Valence") +
  ggtitle("Update Scores by Valence (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Desirability
positive_negative_eligible_aggregate_data %>%
  group_by(desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = desirability, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Desirability") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

##Bar Graph by Group & Valence
positive_negative_eligible_aggregate_data %>%
  group_by(Group, valence) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = valence)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Valence") +
  ggtitle("Update Scores by Group and Valence (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

##Bar Graph by Group & Desirability
positive_negative_eligible_aggregate_data %>%
  group_by(Group, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Group and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

##Bar Graph by Valence & Desirability
positive_negative_eligible_aggregate_data %>%
  group_by(valence, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = valence, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Valence") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Valence and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

##Bar Graph by Group, Valence & Desirability
positive_negative_eligible_aggregate_data %>%
  mutate(valence_labelled = ifelse(valence == "negative", "Negative", "Positive")) %>%
  group_by(Group, valence_labelled, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  facet_wrap(~ valence_labelled) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
    ggtitle("Update Scores by Group, Valence and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

tiff(file="positive_negative_update_group_desirability_int.tiff", units="cm", width=15.91, height=9.82, res=300)
positive_negative_eligible_aggregate_data %>%
  mutate(valence_labelled = ifelse(valence == "negative", "Negative", "Positive")) %>%
  group_by(Group, valence_labelled, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  facet_wrap(~ valence_labelled) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()
dev.off()

#Regression Model
hyp2_model <- lmer(update ~ Group*valence*desirability + (1|subject), data = positive_negative_eligible_aggregate_data)
summ(hyp2_model, confint = TRUE, digits = 3)
hyp2_model_effects <- ggpredict(hyp2_model, terms = c("valence", "desirability", "Group"))
plot(hyp2_model_effects) +
  labs(
    title = "Predicted Values of Update Scores",
    y = "Update Scores",
    x = "Group"
  )

standardize_parameters(hyp2_model, method = "refit", ci = 0.95, include_response = TRUE)

#Group comparisons by desirability across valences
contrasts_hyp2_model_group_desirability <- emmeans(hyp2_model, pairwise ~ desirability | Group, adjust = "tukey")
contrasts_hyp2_model_group_desirability
confint(contrasts_hyp2_model_group_desirability)

#Group comparisons by desirability & valence
contrasts_hyp2_model_group_desirability_valence <- emmeans(hyp2_model, pairwise ~ desirability | Group + valence, adjust = "tukey")
contrasts_hyp2_model_group_desirability_valence
confint(contrasts_hyp2_model_group_desirability_valence)

```

# Hypothesis 3
Hypothesis 3: Initial estimates of negative life events will be associated with depression, such that depressed individuals will on average initially rate their chances of experiencing a negative life event as higher compared to healthy controls.

This hypothesis will be tested using a linear regression using only trials for negative life events. Initial scores will be entered as the outcome, and group (healthy versus depressed) as the predictor. 

```{r}
#Aggregating scores from event level data
negative_initial_estimates <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()


#Descriptive Statistics
negative_initial_estimates %>%
  tab_cols(Group) %>%
  tab_cells(p_estimate1) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Initial Estimates for Negative Life Events by Group")

#Bar Graph by Group
negative_initial_estimates %>%
  group_by(Group) %>%
  summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = Group)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  labs(fill = "Group") +
  ggtitle("Initial Estimates by Group (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

tiff(file="negative_initial_estimate_group.tiff", units="cm", width=15.91, height=9.82, res=300)
negative_initial_estimates %>%
  group_by(Group) %>%
  summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = Group)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  labs(fill = "Group") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()
dev.off()

#Regression Model
hyp3_model <- lm(p_estimate1 ~ Group, data = negative_initial_estimates)
summ(hyp3_model, confint = TRUE, digits = 3)
hyp3_model_effects <- ggpredict(hyp3_model, terms = c("Group"))
plot(hyp3_model_effects) +
  labs(
    title = "Predicted Values of Initial Estimates for Negative Life Events",
    y = "Initial Estimates",
    x = "Group"
  )


standardize_parameters(hyp3_model, method = "refit", ci = 0.95, include_response = TRUE)

```

# Hypothesis 4
Hypothesis 4: Initial estimates of positive life events will be associated with depression, such that depressed individuals will on average initially rate their chances of experiencing a positive life event as lower compared to healthy controls.

```{r}
#Aggregating scores from event level data
positive_initial_estimates <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "positive") %>%
  group_by(subject, Group) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()

#Descriptive Statistics
positive_initial_estimates %>%
  tab_cols(Group) %>%
  tab_cells(p_estimate1) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Initial Estimates for Positive Life Events by Group")

#Bar Graph by Group
positive_initial_estimates %>%
  group_by(Group) %>%
  summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = Group)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  labs(fill = "Group") +
  ggtitle("Initial Estimates by Group (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp4_model <- lm(p_estimate1 ~ Group, data = positive_initial_estimates)
summ(hyp4_model, confint = TRUE, digits = 3)
hyp4_model_effects <- ggpredict(hyp4_model, terms = c("Group"))
plot(hyp4_model_effects) +
  labs(
    title = "Predicted Values of Initial Estimates for Positive Life Events",
    y = "Initial Estimates",
    x = "Group"
  )

standardize_parameters(hyp4_model, method = "refit", ci = 0.95, include_response = TRUE)

```

# Initial Estimates Descriptives and Graph for Positive and Negative Life Events

```{r}
#Aggregating scores from event level data
negative_positive_initial_estimates <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()

#Descriptive Statistics
negative_positive_initial_estimates %>%
  tab_cols(valence %nest% Group) %>%
  tab_cells(p_estimate1) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Initial Estimates for Negative and Positive Life Events by Group")

negative_positive_initial_estimates %>%
  mutate(valence_labelled = ifelse(valence == "negative", "Negative", "Positive")) %>%
  group_by(Group, valence_labelled) %>%
    summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean)) + 
  facet_grid(~valence_labelled) +
  geom_bar(stat = 'identity', position = position_dodge(), fill = "#C7E9C0") +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  theme_minimal()

#Faceted by Valence, Group on X axis, no desirability
tiff(file="positive_negative_initial_estimates_group.tiff", units="cm", width=15.91, height=9.82, res=300)
negative_positive_initial_estimates %>%
  mutate(valence_labelled = ifelse(valence == "negative", "Negative", "Positive")) %>%
  group_by(Group, valence_labelled) %>%
    summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean)) + 
  facet_grid(~valence_labelled) +
  geom_bar(stat = 'identity', position = position_dodge(), fill = "#C7E9C0") +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  theme_minimal()
dev.off()
```

# Descriptives for estimation errors

```{r}
positive_negative_estimation_errors_aggregate <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability) %>%
  mutate(abs_estimation_error = abs(estimation_error)) %>% # Creating absolute values
  summarise(update = mean(update), estimation_error = mean(estimation_error), abs_estimation_error = mean(abs_estimation_error)) %>%
  ungroup()

# Descriptives 
positive_negative_estimation_errors_aggregate %>%
  tab_cols(valence %nest% Group %nest% desirability) %>%
  tab_cells(update, estimation_error, abs_estimation_error) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot(stat_position = "inside_columns")%>%
  set_caption("Update Scores, Raw Estimation Errors and Absolute Estimation Errors by Group, Valence and Desirability")

positive_negative_estimation_errors_aggregate %>%
  mutate(valence_labelled = ifelse(valence == "negative", "Negative", "Positive")) %>%
  group_by(Group, valence_labelled, desirability) %>%
  summarise(mean = mean(abs_estimation_error, na.rm = TRUE), n = n(), sd = sd(abs_estimation_error, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  facet_wrap(~ valence_labelled) +
  ylab("Absolute Estimation Errors") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  ggtitle("Absolute Estimation Errors (Mean & SE)") + 
  scale_fill_brewer(palette = "Purples") + theme_minimal()

```



# Additional Analyses

## Scaled Update Scores

We will investigate whether any observed differences in update scores between depressed individuals and healthy controls are attributable to estimation errors (e.g. the difference between initial estimates and presented average probability). 

Scaled absolute mean update scores will be calculated by dividing absolute mean update scores by absolute mean estimation errors. A mixed-effects linear regression model will be created using only trials for negative life events. Scaled absolute mean update scores will be entered as the outcome, with group (healthy versus depressed), desirability (desirable versus undesirable), and an interaction term between group and desirability entered as predictors. Subject will be entered as a random effect to account for the repeated measures design. 

### Negative Events

```{r}
negative_scaled_update_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group, valence, desirability) %>%
  mutate(abs_estimation_error = abs(estimation_error)) %>% # Creating absolute values
  select(Group, valence, desirability, estimation_error, abs_estimation_error, update) %>%
    summarise(estimation_error = mean(estimation_error), abs_estimation_error = mean(abs_estimation_error), update = mean(update)) %>%
    mutate(scaled_update = update/abs_estimation_error) %>%
  ungroup()
  
# Descriptives 
negative_scaled_update_aggregate_data %>%
  tab_cols(Group, desirability, Group %nest% desirability) %>%
  tab_cells(scaled_update, estimation_error) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot(stat_position = "inside_columns")%>%
  set_caption("Scaled Update Scores and Absolute Estimation Errorsfor Negative Life Events by Group and Desirability")

#Bar Graph by Group
negative_scaled_update_aggregate_data %>%
  group_by(Group) %>%
  summarise(mean = mean(scaled_update, na.rm = TRUE), n = n(), sd = sd(scaled_update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = Group)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Scaled Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Group") +
  ggtitle("Scaled Update Scores for Negative Life Events by Group (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Desirability
negative_scaled_update_aggregate_data %>%
  group_by(desirability) %>%
  summarise(mean = mean(scaled_update, na.rm = TRUE), n = n(), sd = sd(scaled_update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = desirability, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Scaled Update Scores") + 
  xlab ("Desirability") + 
  labs(fill = "Desirability") +
  ggtitle("Scaled Update Scores for Negative Life Events by Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Group & Desirability
negative_scaled_update_aggregate_data %>%
  group_by(Group, desirability) %>%
  summarise(mean = mean(scaled_update, na.rm = TRUE), n = n(), sd = sd(scaled_update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Scaled Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  ggtitle("Scaled Update Scores for Negative Life Events by Group and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
scaled_update_negative_model <- lmer(scaled_update ~ Group*desirability + (1|subject), data = negative_scaled_update_aggregate_data)
summ(scaled_update_negative_model, confint = TRUE, digits = 3)
scaled_update_negative_model_effects <- ggpredict(scaled_update_negative_model, terms = c("Group", "desirability"))
plot(scaled_update_negative_model_effects) +
  labs(
    title = "Predicted Values of Scaled Update Scores for Negative Life Events",
    y = "Scaled Update Scores",
    x = "Group"
  )

standardize_parameters(scaled_update_negative_model, method = "refit", ci = 0.95, include_response = TRUE)

#Comparing Desirability by group
contrasts_scaled_update_negative_model_desirability_by_group <- emmeans(scaled_update_negative_model, pairwise ~ desirability | Group, adjust = "tukey")
contrasts_scaled_update_negative_model_desirability_by_group
confint(contrasts_scaled_update_negative_model_desirability_by_group)
```

This model will also be repeated using both trials for positive and negative events, with valence, and an interaction term between valence and desirability, valence and group, and valence, desirability and group added as additional predictors.

```{r}

positive_negative_scaled_update_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability) %>%
  mutate(abs_estimation_error = abs(estimation_error)) %>% # Creating absolute values
  select(Group, valence, desirability, estimation_error, abs_estimation_error, update) %>%
    summarise(estimation_error = mean(estimation_error), abs_estimation_error = mean(abs_estimation_error), update = mean(update)) %>%
    mutate(scaled_update = update/abs_estimation_error) %>%
  ungroup()
  
# Descriptives 
positive_negative_scaled_update_aggregate_data %>%
  tab_cols(Group, valence, desirability) %>%
  tab_cells(scaled_update, abs_estimation_error) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot(stat_position = "inside_columns")%>%
  set_caption("Scaled Update Scores and Absolute Estimation Errors by Group Valence, and Desirability (Main effects)")

positive_negative_scaled_update_aggregate_data %>%
  tab_cols(Group %nest% valence, Group %nest% desirability, valence %nest% desirability) %>%
  tab_cells(scaled_update, abs_estimation_error) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot(stat_position = "inside_columns")%>%
  set_caption("Scaled Update Scores and Absolute Estimation Errors by Group Valence, and Desirability (2 way Interaction effects)")


positive_negative_scaled_update_aggregate_data %>%
  tab_cols(Group %nest% valence %nest% desirability) %>%
  tab_cells(scaled_update, abs_estimation_error) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot(stat_position = "inside_columns")%>%
  set_caption("Scaled Update Scores and Absolute Estimation Errors by Group Valence, and Desirability (3 way Interaction effects)")

# Switching columns for table 3
positive_negative_scaled_update_aggregate_data %>%
  tab_cols(valence %nest% Group %nest% desirability) %>%
  tab_cells(scaled_update, abs_estimation_error) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot(stat_position = "inside_columns")%>%
  set_caption("Scaled Update Scores and Absolute Estimation Errors by Group Valence, and Desirability (3 way Interaction effects)")


#Bar Graph by Group
positive_negative_scaled_update_aggregate_data %>%
  group_by(Group) %>%
  summarise(mean = mean(scaled_update, na.rm = TRUE), n = n(), sd = sd(scaled_update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = Group)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Scaled Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Group") +
  ggtitle("Scaled Update Scores by Group (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Valence
positive_negative_scaled_update_aggregate_data %>%
  group_by(valence) %>%
  summarise(mean = mean(scaled_update, na.rm = TRUE), n = n(), sd = sd(scaled_update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = valence, y = mean, fill = valence)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Scaled Update Scores") + 
  xlab ("Valence") + 
  labs(fill = "Valence") +
  ggtitle("Scaled Update Scores by Valence (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()


#Bar Graph by Desirability
positive_negative_scaled_update_aggregate_data %>%
  group_by(desirability) %>%
  summarise(mean = mean(scaled_update, na.rm = TRUE), n = n(), sd = sd(scaled_update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = desirability, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Scaled Update Scores") + 
  xlab ("Desirability") + 
  labs(fill = "Desirability") +
  ggtitle("Scaled Update Scores by Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Group & Valence
positive_negative_scaled_update_aggregate_data %>%
  group_by(Group, valence) %>%
  summarise(mean = mean(scaled_update, na.rm = TRUE), n = n(), sd = sd(scaled_update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = valence)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Scaled Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Valence") +
  ggtitle("Scaled Update Scores by Group and Valence (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()


#Bar Graph by Group & Desirability
positive_negative_scaled_update_aggregate_data %>%
  group_by(Group, desirability) %>%
  summarise(mean = mean(scaled_update, na.rm = TRUE), n = n(), sd = sd(scaled_update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Scaled Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  ggtitle("Scaled Update Scores by Group and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Valence & Desirability
positive_negative_scaled_update_aggregate_data %>%
  group_by(valence, desirability) %>%
  summarise(mean = mean(scaled_update, na.rm = TRUE), n = n(), sd = sd(scaled_update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = valence, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Scaled Update Scores") + 
  xlab ("Valence") + 
  labs(fill = "Desirability") +
  ggtitle("Scaled Update Scores by Valence and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Group, Valence & Desirability
positive_negative_scaled_update_aggregate_data %>%
  mutate(valence_labelled = factor(ifelse(valence == "negative", "Negative", "Positive"))) %>%
  group_by(Group, valence_labelled, desirability) %>%
  summarise(mean = mean(scaled_update, na.rm = TRUE), n = n(), sd = sd(scaled_update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  facet_wrap(~valence_labelled) +
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Scaled Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
    ggtitle("Scaled Update Scores by Group, Valence and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Faceted by Valence, Group on X axis
tiff(file="positive_negative_scaled_update_group_desirability_int_v2.tiff", units="cm", width=15.91, height=9.82, res=300)
positive_negative_scaled_update_aggregate_data %>%
  mutate(valence_labelled = factor(ifelse(valence == "negative", "Negative", "Positive"))) %>%
  group_by(Group, valence_labelled, desirability) %>%
  summarise(mean = mean(scaled_update, na.rm = TRUE), n = n(), sd = sd(scaled_update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  facet_wrap(~valence_labelled) +
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Scaled Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()
dev.off()

#Regression Model
scaled_update_positive_negative_model <- lmer(scaled_update ~ Group*valence*desirability + (1|subject), data = positive_negative_scaled_update_aggregate_data)
summ(scaled_update_positive_negative_model, confint = TRUE, digits = 3)
scaled_update_positive_negative_model_effects <- ggpredict(scaled_update_positive_negative_model, terms = c("valence", "desirability", "Group"))
plot(scaled_update_positive_negative_model_effects) +
  labs(
    title = "Predicted Values of Scaled Update Scores",
    y = "Scaled Update Scores",
    x = "Group"
  )

standardize_parameters(scaled_update_positive_negative_model, method = "refit", ci = 0.95, include_response = TRUE)

#Group comparisons by desirability & valence
contrasts_scaled_update_positive_negative_model_group_desirability_valence <- emmeans(scaled_update_positive_negative_model, pairwise ~ desirability | Group + valence, adjust = "tukey")
contrasts_scaled_update_positive_negative_model_group_desirability_valence
confint(contrasts_scaled_update_positive_negative_model_group_desirability_valence)

```

## Dose Response Depression

### Hypothesis 1

#### BDI-II

##### Update Scores

```{r}
negative_eligible_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group, desirability) %>%
  summarise(update = mean(update), BDI_tot = mean(BDI_tot))
  

#Regression Model
hyp1_dose_model <- lmer(update ~ BDI_tot*desirability + Group + (1|subject), data = negative_eligible_aggregate_data)
summ(hyp1_dose_model, confint = TRUE, digits = 3)
hyp1_dose_model_effects <- ggpredict(hyp1_dose_model, terms = c("BDI_tot", "desirability"))
plot(hyp1_dose_model_effects) 

standardize_parameters(hyp1_dose_model, method = "refit", ci = 0.95, include_response = TRUE)

```

#### PHQ

```{r}
negative_eligible_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group, desirability) %>%
  summarise(update = mean(update), PHQ_tot = mean(PHQ_tot)) %>%
  ungroup(0)

#Regression Model
hyp1_dose_phq_model <- lmer(update ~ PHQ_tot*desirability + Group + (1|subject), data = negative_eligible_aggregate_data)
summ(hyp1_dose_phq_model, confint = TRUE, digits = 3)
hyp1_dose_phq_model_effects <- ggpredict(hyp1_dose_phq_model, terms = c("PHQ_tot", "desirability"))
plot(hyp1_dose_phq_model_effects)  +
  labs(
    y = "Update Scores",
    x = "PHQ-9 Scores",
    colour = "Desirability",
    title = ""
  )


standardize_parameters(hyp1_dose_phq_model, method = "refit", ci = 0.95, include_response = TRUE)

tiff(file="neg_PHQ_dose_regression.tiff", units="cm", width=17.51, height=10.82, res=300)
ggplot(data = negative_eligible_aggregate_data, aes(x = PHQ_tot, y = update, colour = desirability)) + geom_point() + geom_smooth(method = "lm") +  ylab("Update Scores") + xlab ("PHQ-9 Scores") +  labs(colour = "Desirability") + theme_minimal()
dev.off()

```

### Hypothesis 2

#### BDI-II

##### Update Scores

```{r}
eligible_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability) %>%
  summarise(update = mean(update), BDI_tot = mean(BDI_tot))
  

#Regression Model
hyp2_dose_BDI_model <- lmer(update ~ BDI_tot*valence*desirability + Group + (1|subject), data = eligible_aggregate_data)
summ(hyp2_dose_BDI_model, confint = TRUE, digits = 3)
hyp2_dose_BDI_model_effects <- ggpredict(hyp2_dose_BDI_model, terms = c("BDI_tot", "desirability", "valence"))
plot(hyp2_dose_BDI_model_effects)

standardize_parameters(hyp2_dose_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
```

#### PHQ

##### Update Scores

```{r}
eligible_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability) %>%
  summarise(update = mean(update), PHQ_tot = mean(PHQ_tot))
  

#Regression Model
hyp2_dose_PHQ_model <- lmer(update ~ PHQ_tot*valence*desirability + Group + (1|subject), data = eligible_aggregate_data)
summ(hyp2_dose_PHQ_model, confint = TRUE, digits = 3)
hyp2_dose_PHQ_model_effects <- ggpredict(hyp2_dose_PHQ_model, terms = c("PHQ_tot", "desirability", "valence"))
plot(hyp2_dose_PHQ_model_effects)

standardize_parameters(hyp2_dose_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)

```


## Adjusting for LOT-R

To assess whether any observed differences in updates according to group are attributable to trait levels of optimism the regression models outlined for hypothesis 1 and 2 will be repeated with LOT-R scores as an additional predictor.

### Hypothesis 1

#### Update Scores

```{r}
negative_eligible_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group, desirability) %>%
  summarise(update = mean(update), LOT_tot = mean(LOT_tot))

#Regression Model
hyp1_LOTR_model <- lmer(update ~ Group*desirability + LOT_tot + (1|subject), data = negative_eligible_aggregate_data)
summ(hyp1_LOTR_model, confint = TRUE, digits = 3)
hyp1_LOTR_model_effects <- ggpredict(hyp1_LOTR_model, terms = c("LOT_tot", "desirability", "Group"))
plot(hyp1_LOTR_model_effects)

standardize_parameters(hyp1_LOTR_model, method = "refit", ci = 0.95, include_response = TRUE)

```

### Hypothesis 2

#### Update Scores

```{r}
positive_negative_eligible_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability) %>%
  summarise(update = mean(update), LOT_tot = mean(LOT_tot))

#Regression Model
hyp2_LOTR_model <- lmer(update ~ Group*valence*desirability + LOT_tot + (1|subject), data = positive_negative_eligible_aggregate_data)
summ(hyp2_LOTR_model, confint = TRUE, digits = 3)
hyp2_LOTR_model_effects <- ggpredict(hyp2_LOTR_model, terms = c("LOT_tot", "desirability", "valence", "Group"))
plot(hyp2_LOTR_model_effects)

standardize_parameters(hyp2_LOTR_model, method = "refit", ci = 0.95, include_response = TRUE)

```


## Memory Errors

### Memory Errors by Group, Valence & Desirability

```{r}
positive_negative_eligible_memory_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability) %>%
  mutate(memory_errors = abs(memory_errors)) %>% # Creating absolute values
  summarise(memory_errors = mean(memory_errors, na.rm = TRUE))


#By Group
positive_negative_eligible_memory_aggregate_data %>%
  tab_cols(Group) %>%
  tab_cells(memory_errors) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Absolute Memory Errors by Group")

positive_negative_eligible_memory_aggregate_data %>%
  group_by(Group) %>%
  summarise(mean = mean(memory_errors, na.rm = TRUE), n = n(), sd = sd(memory_errors, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = Group)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Absolute Memory Errors") + 
  xlab ("Group") + 
  labs(fill = "Group") +
  ggtitle("Absolute Memory Errors by Group (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#By Valence
positive_negative_eligible_memory_aggregate_data %>%
  tab_cols(valence) %>%
  tab_cells(memory_errors) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Absolute Memory Errors by Valence")

positive_negative_eligible_memory_aggregate_data %>%
  group_by(valence) %>%
  summarise(mean = mean(memory_errors, na.rm = TRUE), n = n(), sd = sd(memory_errors, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = valence, y = mean, fill = valence)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Absolute Memory Errors") + 
  xlab ("Valence") + 
  labs(fill = "Valence") +
  ggtitle("Absolute Memory Errors by Valence (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#By Desirability
positive_negative_eligible_memory_aggregate_data %>%
  tab_cols(desirability) %>%
  tab_cells(memory_errors) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Absolute Memory Errors by Desirability")

positive_negative_eligible_memory_aggregate_data %>%
  group_by(desirability) %>%
  summarise(mean = mean(memory_errors, na.rm = TRUE), n = n(), sd = sd(memory_errors, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = desirability, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Absolute Memory Errors") + 
  xlab ("Desirability") + 
  labs(fill = "Desirability") +
  ggtitle("Absolute Memory Errors by Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

##By Group & Valence
positive_negative_eligible_memory_aggregate_data %>%
  tab_cols(Group %nest% valence) %>%
  tab_cells(memory_errors) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Absolute Memory Errors by Group and Valence")

positive_negative_eligible_memory_aggregate_data %>%
  group_by(Group, valence) %>%
  summarise(mean = mean(memory_errors, na.rm = TRUE), n = n(), sd = sd(memory_errors, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = valence)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Absolute Memory Errors") + 
  xlab ("Group") + 
  labs(fill = "Valence") +
  ggtitle("Absolute Memory Errors by Group and Valence (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

##By Group & Desirability
positive_negative_eligible_memory_aggregate_data %>%
  tab_cols(Group %nest% desirability) %>%
  tab_cells(memory_errors) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Absolute Memory Errors by Group and Desirability")

positive_negative_eligible_memory_aggregate_data %>%
  group_by(Group, desirability) %>%
  summarise(mean = mean(memory_errors, na.rm = TRUE), n = n(), sd = sd(memory_errors, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Absolute Memory Errors") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  ggtitle("Absolute Memory Errors by Group and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

## By Valence & Desirability
positive_negative_eligible_memory_aggregate_data %>%
  tab_cols(valence %nest% desirability) %>%
  tab_cells(memory_errors) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Absolute Memory Errors by Valence and Desirability")

positive_negative_eligible_memory_aggregate_data %>%
  group_by(valence, desirability) %>%
  summarise(mean = mean(memory_errors, na.rm = TRUE), n = n(), sd = sd(memory_errors, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = valence, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Absolute Memory Errors") + 
  xlab ("Valence") + 
  labs(fill = "Desirability") +
  ggtitle("Absolute Memory Errors by Valence and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

### By Group, Valence & Desirability
positive_negative_eligible_memory_aggregate_data %>%
  tab_cols(Group %nest% valence %nest% desirability) %>%
  tab_cells(memory_errors) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Absolute Memory Errors by Group, Valence and Desirability")

#Switching columns for table 3
positive_negative_eligible_memory_aggregate_data %>%
  tab_cols(valence %nest% Group %nest% desirability) %>%
  tab_cells(memory_errors) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Absolute Memory Errors by Group, Valence and Desirability")

positive_negative_eligible_memory_aggregate_data %>%
  mutate(valence_labelled = ifelse(valence == "negative", "Negative", "Positive")) %>%
  group_by(Group, valence_labelled, desirability) %>%
  summarise(mean = mean(memory_errors, na.rm = TRUE), n = n(), sd = sd(memory_errors, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  facet_wrap(~ valence_labelled) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  ggtitle("Absolute Memory Errors by Group, Valence and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
memory_model <- lmer(memory_errors ~ Group*valence*desirability + (1|subject), data = positive_negative_eligible_memory_aggregate_data)
summ(memory_model, confint = TRUE, digits = 3)
memory_model_effects <- ggpredict(memory_model, terms = c("valence", "desirability", "Group"))
plot(memory_model_effects)

standardize_parameters(memory_model, method = "refit", ci = 0.95, include_response = TRUE)

```

### Adjusting Hypothesis 1 for Memory Errors

```{r}
negative_eligible_memory_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group, desirability) %>%
  mutate(memory_errors = abs(memory_errors)) %>% # Creating absolute values
  summarise(memory_errors = mean(memory_errors), update = mean(update))

#Regression Model
hyp1_memory_model <- lmer(update ~ Group*desirability + memory_errors + (1|subject), data = negative_eligible_memory_aggregate_data)
summ(hyp1_memory_model, confint = TRUE, digits = 3)
hyp1_memory_model_effects <- ggpredict(hyp1_memory_model, terms = c("Group", "desirability"))
plot(hyp1_memory_model_effects)

standardize_parameters(hyp1_memory_model, method = "refit", ci = 0.95, include_response = TRUE)
```

### Adjusting Hypothesis 2 for Memory Errors

```{r}
eligible_memory_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability) %>%
  mutate(memory_errors = abs(memory_errors)) %>% # Creating absolute values
  summarise(memory_errors = mean(memory_errors), update = mean(update)) %>%
  ungroup()


hyp2_memory_model <- lmer(update ~ Group*valence*desirability + memory_errors + (1|subject), data = eligible_memory_aggregate_data)
summ(hyp2_memory_model, confint = TRUE, digits = 3)
hyp2_memory_model_effects <- ggpredict(hyp2_memory_model, terms = c("valence", "desirability", "Group"))
plot(hyp2_memory_model_effects)

standardize_parameters(hyp2_memory_model, method = "refit", ci = 0.95, include_response = TRUE)
```

## Framing

### Hypothesis 1

```{r}
negative_eligible_aggregate_framing_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group, desirability, happenCondition) %>%
  summarise(update = mean(update))

negative_eligible_aggregate_framing_data %>%
  tab_cols(Group %nest% desirability %nest% happenCondition) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores for Negative Life Events by Group, Desirability and Framing")

negative_eligible_aggregate_framing_data %>%
  group_by(Group, desirability, happenCondition) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = desirability, y = mean, fill = happenCondition)) + 
  facet_wrap(~ Group) +
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Desirability") + 
  labs(fill = "Framing") +
  ggtitle("Update Scores for Negative Life Events by Group, Desirability and Framing (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp1_framing_model <- lmer(update ~ Group*desirability + happenCondition + (1|subject), data = negative_eligible_aggregate_framing_data)
summ(hyp1_framing_model, confint = TRUE, digits = 3)
hyp1_framing_model_effects <- ggpredict(hyp1_framing_model, terms = c("desirability", "happenCondition","Group"))
plot(hyp1_framing_model_effects)

standardize_parameters(hyp1_framing_model, method = "refit", ci = 0.95, include_response = TRUE)
```

### Hypothesis 2

```{r}
positive_negative_framing_eligible_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability, happenCondition) %>%
  summarise(update = mean(update)) %>%
  ungroup()

positive_negative_framing_eligible_aggregate_data %>%
  tab_cols(Group %nest% valence %nest% desirability %nest% happenCondition) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores by Group, Valence, Desirability and Framing")

positive_negative_framing_eligible_aggregate_data %>%
    mutate(valence_labelled = ifelse(valence == "negative", "Negative", "Positive")) %>%
  group_by(Group, valence_labelled, desirability, happenCondition) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = desirability, y = mean, fill = happenCondition)) + 
  facet_wrap(~ valence_labelled + Group) +
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Desirability") + 
  labs(fill = "Framing") +
  ggtitle("Update Scores by Group, Valence, Desirability and Framing (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp2_framing_model <- lmer(update ~ Group*valence*desirability + happenCondition + (1|subject), data = positive_negative_framing_eligible_aggregate_data)
summ(hyp2_framing_model, confint = TRUE, digits = 3)
hyp2_framing_model_effects <- ggpredict(hyp2_framing_model, terms = c("desirability", "happenCondition", "valence", "Group"))
plot(hyp2_framing_model_effects)

standardize_parameters(hyp2_framing_model, method = "refit", ci = 0.95, include_response = TRUE)
```

### Hypothesis 3

```{r}

#Aggregating scores from event level data
negative_initial_estimates_participant_framing <- task_event_level_qs %>%
  filter(valence == "negative" & eligible_trial == "Eligible") %>%
  group_by(subject, Group, happenCondition) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()

#Descriptive Statistics
negative_initial_estimates_participant_framing %>%
  tab_cols(Group %nest% happenCondition) %>%
  tab_cells(p_estimate1) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Initial Estimates for Negative Life Events by Group and Framing")

#Bar Graph by Group
negative_initial_estimates_participant_framing %>%
  group_by(Group, happenCondition) %>%
  summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = happenCondition)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  labs(fill = "Happening Condition") +
  ggtitle("Initial Estimates for Negative Life Events by Group and Framing (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp3_framing_model <- lmer(p_estimate1 ~ Group + happenCondition + (1|subject), data = negative_initial_estimates_participant_framing)
summ(hyp3_framing_model, confint = TRUE, digits = 3)
hyp3_framing_model_effects <- ggpredict(hyp3_framing_model, terms = c("Group"))
plot(hyp3_framing_model_effects) +
  labs(
    title = "Predicted Values of Initial Estimates for Negative Life Events",
    y = "Initial Estimates",
    x = "Group"
  )
standardize_parameters(hyp3_framing_model, method = "refit", ci = 0.95, include_response = TRUE)

```

### Hypothesis 4

```{r}
#Aggregating scores from event level data
positive_initial_estimates_participant_framing <- task_event_level_qs %>%
  filter(valence == "positive" & eligible_trial == "Eligible") %>%
  group_by(subject, Group, happenCondition) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()

#Descriptive Statistics
positive_initial_estimates_participant_framing %>%
  tab_cols(Group %nest% happenCondition) %>%
  tab_cells(p_estimate1) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Initial Estimates for positive Life Events by Group and Framing")

#Bar Graph by Group
positive_initial_estimates_participant_framing %>%
  group_by(Group, happenCondition) %>%
  summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = happenCondition)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  labs(fill = "Happening Condition") +
  ggtitle("Initial Estimates for positive Life Events by Group and Framing (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp4_framing_model <- lmer(p_estimate1 ~ Group + happenCondition + (1|subject), data = positive_initial_estimates_participant_framing)
summ(hyp4_framing_model, confint = TRUE, digits = 3)
hyp4_framing_model_effects <- ggpredict(hyp4_framing_model, terms = c("Group"))
plot(hyp4_framing_model_effects) +
  labs(
    title = "Predicted Values of Initial Estimates for positive Life Events",
    y = "Initial Estimates",
    x = "Group"
  )
standardize_parameters(hyp4_framing_model, method = "refit", ci = 0.95, include_response = TRUE)
```

## Valence Order Effects

### Hypothesis 1

```{r}
#Creating a variable to show order effects (see supplementary table 5 for full details)
task_event_level_qs <- task_event_level_qs %>%
  mutate(valence_order = ifelse(task_group == 1 | task_group == 2 | task_group == 3 | task_group == 4, 1, 2))

task_event_level_qs$valence_order <- factor(task_event_level_qs$valence_order,
                                            levels = c(1, 2),
                                            labels = c("Negative First", "Positive First"))


negative_eligible_aggregate_valence_order_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group, desirability, valence_order) %>%
  summarise(update = mean(update))

negative_eligible_aggregate_valence_order_data %>%
  tab_cols(Group %nest% desirability %nest% valence_order) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores for Negative Life Events by Group, Desirability and Valence Order")

negative_eligible_aggregate_valence_order_data %>%
  group_by(Group, desirability, valence_order) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = desirability, y = mean, fill = valence_order)) + 
  facet_wrap(~ Group) +
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Desirability") + 
  labs(fill = "Valence Order") +
  ggtitle("Update Scores for Negative Life Events by Group, Desirability and Valence Order (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()


#Regression Model
hyp1_valence_order_model <- lmer(update ~ Group*desirability + valence_order + (1|subject), data = negative_eligible_aggregate_valence_order_data)
summ(hyp1_valence_order_model, confint = TRUE, digits = 3)
hyp1_valence_order_model_effects <- ggpredict(hyp1_valence_order_model, terms = c("desirability", "valence_order","Group"))
plot(hyp1_valence_order_model_effects)

standardize_parameters(hyp1_valence_order_model, method = "refit", ci = 0.95, include_response = TRUE)

```


### Hypothesis 2

```{r}
positive_negative_valence_order_eligible_aggregate_data <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability, valence_order) %>%
  summarise(update = mean(update))

positive_negative_valence_order_eligible_aggregate_data %>%
  tab_cols(Group %nest% valence %nest% desirability %nest% valence_order) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores by Group, Valence, Desirability and Valence Order")

positive_negative_valence_order_eligible_aggregate_data %>%
  group_by(Group, valence, desirability, valence_order) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = desirability, y = mean, fill = valence_order)) + 
  facet_wrap(~ valence + Group) +
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Desirability") + 
  labs(fill = "Valence Order") +
  ggtitle("Update Scores by Group, Valence, Desirability and Valence Order (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp2_valence_order_model <- lmer(update ~ Group*valence*desirability + valence_order + (1|subject), data = positive_negative_valence_order_eligible_aggregate_data)
summ(hyp2_valence_order_model, confint = TRUE, digits = 3)
hyp2_valence_order_model_effects <- ggpredict(hyp2_valence_order_model, terms = c("desirability", "valence_order", "valence", "Group"))
plot(hyp2_valence_order_model_effects)

standardize_parameters(hyp2_valence_order_model, method = "refit", ci = 0.95, include_response = TRUE)


```

### Hypothesis 3

```{r}

#Aggregating scores from event level data
negative_initial_estimates_valence_order <- task_event_level_qs %>%
  filter(valence == "negative" & eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence_order) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()

#Descriptive Statistics
negative_initial_estimates_valence_order %>%
  tab_cols(Group %nest% valence_order) %>%
  tab_cells(p_estimate1) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Initial Estimates for Negative Life Events by Group and Valence Order")

#Bar Graph by Group
negative_initial_estimates_valence_order %>%
  group_by(Group, valence_order) %>%
  summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = valence_order)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  labs(fill = "Valence Order") +
  ggtitle("Initial Estimates for Negative Life Events by Group and Valence Order (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp3_valenceorder_model <- lm(p_estimate1 ~ Group + valence_order, data = negative_initial_estimates_valence_order)
summ(hyp3_valenceorder_model, confint = TRUE, digits = 3)
hyp3_valenceorder_model_effects <- ggpredict(hyp3_valenceorder_model, terms = c("Group"))
plot(hyp3_valenceorder_model_effects) +
  labs(
    title = "Predicted Values of Initial Estimates for Negative Life Events",
    y = "Initial Estimates",
    x = "Group"
  )
standardize_parameters(hyp3_valenceorder_model, method = "refit", ci = 0.95, include_response = TRUE)

```

### Hypothesis 4

```{r}
#Aggregating scores from event level data
positive_initial_estimates_valence_order <- task_event_level_qs %>%
  filter(valence == "positive" & eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence_order) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()

#Descriptive Statistics
positive_initial_estimates_valence_order %>%
  tab_cols(Group %nest% valence_order) %>%
  tab_cells(p_estimate1) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Initial Estimates for positive Life Events by Group and Valence Order")

#Bar Graph by Group
positive_initial_estimates_valence_order %>%
  group_by(Group, valence_order) %>%
  summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = valence_order)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  labs(fill = "Valence Order") +
  ggtitle("Initial Estimates for positive Life Events by Group and Valence Order (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp4_valenceorder_model <- lm(p_estimate1 ~ Group + valence_order, data = positive_initial_estimates_valence_order)
summ(hyp4_valenceorder_model, confint = TRUE, digits = 3)
hyp4_valenceorder_model_effects <- ggpredict(hyp4_valenceorder_model, terms = c("Group"))
plot(hyp4_valenceorder_model_effects) +
  labs(
    title = "Predicted Values of Initial Estimates for positive Life Events",
    y = "Initial Estimates",
    x = "Group"
  )
standardize_parameters(hyp4_valenceorder_model, method = "refit", ci = 0.95, include_response = TRUE)
```

## Ratings

### Descriptive Statistics for Ratings

```{r}
ratings_wide <- task_event_level_qs %>%
  group_by(subject, Group, valence, desirability) %>%
  summarise(rating_controllability = mean(rating_controllability, na.rm = TRUE),
            rating_emotionalarousal = mean(rating_emotionalarousal, na.rm = TRUE),
            rating_familiarity = mean(rating_familiarity, na.rm = TRUE),
            rating_negativity = mean(rating_negativity, na.rm = TRUE),
            rating_positivity = mean(rating_positivity, na.rm = TRUE),
            rating_priorexperience = mean(rating_priorexperience, na.rm = TRUE),
            rating_vividness = mean(rating_vividness, na.rm = TRUE)) %>%
  ungroup()

ratings_long <- ratings_wide %>%
  pivot_longer(cols = starts_with("rating_"), names_prefix = "rating_", names_to = "rating_type", values_to = "rating") %>%
  mutate(rating_type = str_to_title(rating_type)) %>%
  mutate(rating_type = ifelse(rating_type == "Emotionalarousal", "Emotional Arousal",
                              ifelse(rating_type == "Priorexperience", "Prior Experience", rating_type))) #Renaming for graphs

  
ratings_long %>%
  tab_cols(Group %nest% valence) %>%
  tab_rows(rating_type) %>%
  tab_cells(rating) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_columns")%>%
  set_caption("Event Ratings by Group and Valence")


ratings_long %>%
  tab_cols(Group %nest% valence %nest% desirability) %>%
  tab_rows(rating_type) %>%
  tab_cells(rating) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_columns")%>%
  set_caption("Event Ratings by Group, Valence, Desirability")

#Switching columns for table 3
ratings_long %>%
  tab_cols(valence %nest% Group %nest% desirability) %>%
  tab_rows(rating_type) %>%
  tab_cells(rating) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_columns")%>%
  set_caption("Event Ratings by Group, Valence, Desirability")

tiff(file="ratings_group_valence.tiff", units="cm", width=27, height=12, res=300)
ratings_long %>%
  filter(!is.na(desirability)) %>%
    mutate(valence_labelled = ifelse(valence == "negative", "Negative", "Positive")) %>%
  group_by(Group, valence_labelled, rating_type) %>%
  summarise(mean = mean(rating, na.rm = TRUE), n = n(), sd = sd(rating, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = valence_labelled)) + 
  facet_grid(~ rating_type) +
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Rating") + 
  xlab ("Group") + 
  labs(fill = "Valence") +
  scale_fill_brewer(palette = "Reds") + 
  theme_minimal() + 
  theme(strip.background = element_rect(color="grey", fill="#FFFFFF", size=1, linetype="solid"))
dev.off()


ratings_long %>%
  filter(!is.na(desirability)) %>%
  group_by(Group, valence, desirability, rating_type) %>%
  summarise(mean = mean(rating, na.rm = TRUE), n = n(), sd = sd(rating, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  facet_wrap(~ rating_type + valence) +
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Rating") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  ggtitle("Ratings by Group, Valence and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()



```



### Adjusting Hypothesis 1 for Ratings

```{r}
#Aggregating scores from event level data
negative_eligible_aggregate_data_ratings <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group, desirability) %>%
  summarise(update = mean(update, na.rm = TRUE), rating_controllability = mean(rating_controllability, na.rm = TRUE), rating_emotionalarousal = mean(rating_emotionalarousal, na.rm = TRUE), rating_familiarity = mean(rating_familiarity, na.rm = TRUE), rating_negativity = mean(rating_negativity, na.rm = TRUE), rating_positivity = mean(rating_positivity, na.rm = TRUE), rating_priorexperience = mean(rating_priorexperience, na.rm = TRUE), rating_vividness = mean(rating_vividness, na.rm = TRUE)) %>%
  ungroup()

#Regression Model
hyp1_adj_ratings_model <- lmer(update ~ Group*desirability + rating_controllability + rating_emotionalarousal + rating_familiarity + rating_negativity + rating_positivity + rating_priorexperience + rating_vividness + (1|subject), data = negative_eligible_aggregate_data_ratings)
summ(hyp1_adj_ratings_model, confint = TRUE, digits = 3)
hyp1_adj_ratings_model_effects <- ggpredict(hyp1_adj_ratings_model, terms = c("Group", "desirability"))
plot(hyp1_adj_ratings_model_effects) +
  labs(
    title = "Predicted Values of Update Scores for Negative Life Events",
    y = "Update Scores",
    x = "Group"
  )

standardize_parameters(hyp1_adj_ratings_model, method = "refit", ci = 0.95, include_response = TRUE)

```

### Adjusting Hypothesis 2 for ratings

```{r}
#Aggregating Event Level Data
positive_negative_eligible_aggregate_data_ratings <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability) %>%
  summarise(update = mean(update, na.rm = TRUE), rating_controllability = mean(rating_controllability, na.rm = TRUE), rating_emotionalarousal = mean(rating_emotionalarousal, na.rm = TRUE), rating_familiarity = mean(rating_familiarity, na.rm = TRUE), rating_negativity = mean(rating_negativity, na.rm = TRUE), rating_positivity = mean(rating_positivity, na.rm = TRUE), rating_priorexperience = mean(rating_priorexperience, na.rm = TRUE), rating_vividness = mean(rating_vividness, na.rm = TRUE))

#Regression Model
hyp2_adj_ratings_model <- lmer(update ~ Group*valence*desirability + rating_controllability + rating_emotionalarousal + rating_familiarity + rating_negativity + rating_positivity + rating_priorexperience + rating_vividness + (1|subject), data = positive_negative_eligible_aggregate_data_ratings)
summ(hyp2_adj_ratings_model, confint = TRUE, digits = 3)
hyp2_adj_ratings_model_effects <- ggpredict(hyp2_adj_ratings_model, terms = c("valence", "desirability", "Group"))
plot(hyp2_adj_ratings_model_effects) +
  labs(
    title = "Predicted Values of Update Scores",
    y = "Update Scores",
    x = "Group"
  )

standardize_parameters(hyp2_adj_ratings_model, method = "refit", ci = 0.95, include_response = TRUE)
```

### Hypothesis 3

```{r}
#Aggregating scores from event level data
negative_initial_estimates_ratings <- task_event_level_qs %>%
  filter(valence == "negative" & eligible_trial == "Eligible") %>%
  group_by(subject, Group) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE), rating_controllability = mean(rating_controllability, na.rm = TRUE), rating_emotionalarousal = mean(rating_emotionalarousal, na.rm = TRUE), rating_familiarity = mean(rating_familiarity, na.rm = TRUE), rating_negativity = mean(rating_negativity, na.rm = TRUE), rating_positivity = mean(rating_positivity, na.rm = TRUE), rating_priorexperience = mean(rating_priorexperience, na.rm = TRUE), rating_vividness = mean(rating_vividness, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()


#Regression Model
hyp3_ratings_model <- lm(p_estimate1 ~ Group + rating_controllability + rating_emotionalarousal + rating_familiarity + rating_negativity + rating_positivity + rating_priorexperience + rating_vividness, data = negative_initial_estimates_ratings)
summ(hyp3_ratings_model, confint = TRUE, digits = 3)
hyp3_ratings_model_effects <- ggpredict(hyp3_ratings_model, terms = c("Group"))
plot(hyp3_ratings_model_effects) +
  labs(
    title = "Predicted Values of Initial Estimates for Negative Life Events",
    y = "Initial Estimates",
    x = "Group"
  )
standardize_parameters(hyp3_ratings_model, method = "refit", ci = 0.95, include_response = TRUE)
```

### Hypothesis 4

```{r}
#Aggregating scores from event level data
positive_initial_estimates_ratings <- task_event_level_qs %>%
  filter(valence == "positive" & eligible_trial == "Eligible") %>%
  group_by(subject, Group) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE), rating_controllability = mean(rating_controllability, na.rm = TRUE), rating_emotionalarousal = mean(rating_emotionalarousal, na.rm = TRUE), rating_familiarity = mean(rating_familiarity, na.rm = TRUE), rating_negativity = mean(rating_negativity, na.rm = TRUE), rating_positivity = mean(rating_positivity, na.rm = TRUE), rating_priorexperience = mean(rating_priorexperience, na.rm = TRUE), rating_vividness = mean(rating_vividness, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()


#Regression Model
hyp4_ratings_model <- lm(p_estimate1 ~ Group + rating_controllability + rating_emotionalarousal + rating_familiarity + rating_negativity + rating_positivity + rating_priorexperience + rating_vividness, data = positive_initial_estimates_ratings)
summ(hyp4_ratings_model, confint = TRUE, digits = 3)
hyp4_ratings_model_effects <- ggpredict(hyp4_ratings_model, terms = c("Group"))
plot(hyp4_ratings_model_effects) +
  labs(
    title = "Predicted Values of Initial Estimates for positive Life Events",
    y = "Initial Estimates",
    x = "Group"
  )
standardize_parameters(hyp4_ratings_model, method = "refit", ci = 0.95, include_response = TRUE)
```


## Categorising event valence according to participant ratings

```{r}

# Averaging positivity & negativity ratings
valence_participant_ratings <- task_event_level_qs %>%
  mutate(reverse_code_rating_positivity = dplyr::recode(.$rating_positivity, `1` = 6L, `2` = 5L, `3` = 4L, `4` = 3L, `5` = 2L, `6` = 1L)) %>% #Reverse coding positive ratings
  mutate(mean_rating_valence = (rating_negativity + reverse_code_rating_positivity)/2) %>% #Averaging across negative & reverse coded positivity ratings (higher values = more negative & less positive)
  mutate(mean_rating_valence = ifelse(is.na(rating_negativity) & !is.na(reverse_code_rating_positivity), reverse_code_rating_positivity,
                                      ifelse(is.na(reverse_code_rating_positivity) & !is.na(rating_negativity), rating_negativity, mean_rating_valence))) %>% #If positive rating is missing just use negative and vice versa
  mutate(participant_valence_avg_rating = ifelse(mean_rating_valence > 3.5, 1,
                                                 ifelse(mean_rating_valence < 3.5, 2,
                                                        ifelse(mean_rating_valence == 3.5, 3, NA)))) #Categorising events so greater than 3.5 = negative, less than 3.5 = positive, and == 3.5 = neutral

valence_participant_ratings$participant_valence_avg_rating <- factor(valence_participant_ratings$participant_valence_avg_rating,
                                                          levels = c(1, 2, 3),
                                                          labels = c("Negative", "Positive", "Neutral"))

table(valence_participant_ratings$valence, valence_participant_ratings$participant_valence_avg_rating, useNA = "ifany")


valence_participant_ratings_selected <- valence_participant_ratings %>%
  select(subject, lifeEvent, valence, rating_negativity, rating_positivity, reverse_code_rating_positivity, mean_rating_valence, participant_valence_avg_rating)
```

### Hypothesis 1

```{r}
#Aggregating scores from event level data
negative_eligible_aggregate_data_participant_valence <- valence_participant_ratings %>%
  filter(eligible_trial == "Eligible" & participant_valence_avg_rating == "Negative") %>%
  group_by(subject, Group, desirability) %>%
  summarise(update = mean(update, na.rm = TRUE))

#Descriptive Statistics
negative_eligible_aggregate_data_participant_valence %>%
  tab_cols(Group, desirability, Group %nest% desirability) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores for Negative Life Events (according to participant ratings) by Group and Desirability")

#Bar Graph by Group
negative_eligible_aggregate_data_participant_valence %>%
  group_by(Group) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = Group)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Group") +
  ggtitle("Update Scores by Group (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Desirability
negative_eligible_aggregate_data_participant_valence %>%
  group_by(desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = desirability, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Desirability") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Group and Desirability
negative_eligible_aggregate_data_participant_valence %>%
  group_by(Group, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Group & Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp1_participant_valence_model <- lmer(update ~ Group*desirability + (1|subject), data = negative_eligible_aggregate_data_participant_valence)
summ(hyp1_participant_valence_model, confint = TRUE, digits = 3)
hyp1_participant_valence_model_effects <- ggpredict(hyp1_participant_valence_model, terms = c("Group", "desirability"))
plot(hyp1_participant_valence_model_effects) +
  labs(
    title = "Predicted Values of Update Scores for Negative Life Events (according to participant ratings)",
    y = "Update Scores",
    x = "Group"
  )

standardize_parameters(hyp1_participant_valence_model, method = "refit", ci = 0.95, include_response = TRUE)
```

### Hypothesis 2

```{r}
#Aggregating Event Level Data
positive_negative_eligible_aggregate_data_participant_valence <- valence_participant_ratings %>%
  filter(eligible_trial == "Eligible") %>%
  filter(participant_valence_avg_rating == "Positive" | participant_valence_avg_rating == "Negative") %>%
  group_by(subject, Group, participant_valence_avg_rating, desirability) %>%
  summarise(update = mean(update, na.rm = TRUE))

#Descriptives
positive_negative_eligible_aggregate_data_participant_valence %>%
  tab_cols(Group, participant_valence_avg_rating, desirability) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores Group, Valence (participant ratings) and Desirability (Main Effects)")


positive_negative_eligible_aggregate_data_participant_valence %>%
  tab_cols(Group %nest% participant_valence_avg_rating, Group %nest% desirability, participant_valence_avg_rating %nest% desirability, Group %nest% participant_valence_avg_rating %nest% desirability) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores Group, Valence (participant ratings) and Desirability (Interaction Effects)")


#Bar Graph by Group
positive_negative_eligible_aggregate_data_participant_valence %>%
  group_by(Group) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = Group)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Group") +
  ggtitle("Update Scores by Group (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Valence
positive_negative_eligible_aggregate_data_participant_valence %>%
  group_by(participant_valence_avg_rating) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = participant_valence_avg_rating, y = mean, fill = participant_valence_avg_rating)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Valence (participant ratings)") + 
  labs(fill = "Valence (participant ratings)") +
  ggtitle("Update Scores by Valence (participant ratings) (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Desirability
positive_negative_eligible_aggregate_data_participant_valence %>%
  group_by(desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = desirability, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Desirability") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

##Bar Graph by Group & Valence
positive_negative_eligible_aggregate_data_participant_valence %>%
  group_by(Group, participant_valence_avg_rating) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = participant_valence_avg_rating)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Valence (participant ratings)") +
  ggtitle("Update Scores by Group and Valence (participant ratings) (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

##Bar Graph by Group & Desirability
positive_negative_eligible_aggregate_data_participant_valence %>%
  group_by(Group, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Group and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

##Bar Graph by Valence & Desirability
positive_negative_eligible_aggregate_data_participant_valence %>%
  group_by(participant_valence_avg_rating, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = participant_valence_avg_rating, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Valence (participant ratings)") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Valence and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

##Bar Graph by Group, Valence & Desirability
positive_negative_eligible_aggregate_data_participant_valence %>%
  group_by(Group, participant_valence_avg_rating, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  facet_wrap(~ participant_valence_avg_rating) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Group, Valence (according to participants) and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp2_participant_valence_model <- lmer(update ~ Group*participant_valence_avg_rating*desirability + (1|subject), data = positive_negative_eligible_aggregate_data_participant_valence)
summ(hyp2_participant_valence_model, confint = TRUE, digits = 3)
hyp2_participant_valence_model_effects <- ggpredict(hyp2_participant_valence_model, terms = c("participant_valence_avg_rating", "desirability", "Group"))
plot(hyp2_participant_valence_model_effects)

standardize_parameters(hyp2_participant_valence_model, method = "refit", ci = 0.95, include_response = TRUE)
```

### Hypothesis 3

```{r}
#Aggregating scores from event level data
negative_initial_estimates_participant_valence <- valence_participant_ratings %>%
  filter(participant_valence_avg_rating == "Negative") %>%
  group_by(subject, Group) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()

#Descriptive Statistics
negative_initial_estimates_participant_valence %>%
  tab_cols(Group) %>%
  tab_cells(p_estimate1) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Initial Estimates for Negative Life Events (according to participant ratings) by Group")

#Bar Graph by Group
negative_initial_estimates_participant_valence %>%
  group_by(Group) %>%
  summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = Group)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  labs(fill = "Group") +
  ggtitle("Initial Estimates for Negative Life Events (according to participant ratings) by Group (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp3_participant_valence_model <- lm(p_estimate1 ~ Group, data = negative_initial_estimates_participant_valence)
summ(hyp3_participant_valence_model, confint = TRUE, digits = 3)
hyp3_participant_valence_model_effects <- ggpredict(hyp3_participant_valence_model, terms = c("Group"))
plot(hyp3_participant_valence_model_effects) +
  labs(
    title = "Predicted Values of Initial Estimates for Negative Life Events",
    y = "Initial Estimates",
    x = "Group"
  )
standardize_parameters(hyp3_participant_valence_model, method = "refit", ci = 0.95, include_response = TRUE)

```

### Hypothesis 4

```{r}
#Aggregating scores from event level data
positive_initial_estimates_participant_valence <- valence_participant_ratings %>%
  filter(participant_valence_avg_rating == "Positive") %>%
  group_by(subject, Group) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()

#Descriptive Statistics
positive_initial_estimates_participant_valence %>%
  tab_cols(Group) %>%
  tab_cells(p_estimate1) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Initial Estimates for Positive Life Events (according to participant ratings) by Group")

#Bar Graph by Group
positive_initial_estimates_participant_valence %>%
  group_by(Group) %>%
  summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = Group)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  labs(fill = "Group") +
  ggtitle("Initial Estimates by Group (according to participant ratings) (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Regression Model
hyp4_participant_valence_model <- lm(p_estimate1 ~ Group, data = positive_initial_estimates_participant_valence)
summ(hyp4_participant_valence_model, confint = TRUE, digits = 3)
hyp4_participant_valence_model_effects <- ggpredict(hyp4_participant_valence_model, terms = c("Group"))
plot(hyp4_participant_valence_model_effects) +
  labs(
    title = "Predicted Values of Initial Estimates for Positive Life Events (according to participant ratings)",
    y = "Initial Estimates",
    x = "Group"
  )
standardize_parameters(hyp4_participant_valence_model, method = "refit", ci = 0.95, include_response = TRUE)

```


### Calculating mean ratings per life event

```{r}
mean_event_level_ratings <- task_event_level_qs %>%
  group_by(Group, lifeEvent, valence) %>%
  summarise_at(c("rating_controllability", "rating_emotionalarousal", "rating_familiarity", "rating_negativity", "rating_positivity", "rating_priorexperience", "rating_vividness"), ~ mean(.x, na.rm = TRUE))
colnames(mean_event_level_ratings)[4:10] <- paste("mean", colnames(mean_event_level_ratings)[4:10],  sep = "_") 

sd_event_level_ratings <- task_event_level_qs %>%
  group_by(Group, lifeEvent, valence) %>%
  summarise_at(c("rating_controllability", "rating_emotionalarousal", "rating_familiarity", "rating_negativity", "rating_positivity", "rating_priorexperience", "rating_vividness"), ~ sd(.x, na.rm = TRUE))
colnames(sd_event_level_ratings)[4:10] <- paste("sd", colnames(sd_event_level_ratings)[4:10],  sep = "_") 
 
mean_sd_event_level_ratings <- merge(mean_event_level_ratings, sd_event_level_ratings, by = c("Group", "lifeEvent", "valence"))

mean_sd_event_level_ratings <- mean_sd_event_level_ratings %>%
  pivot_wider(names_from = "Group", values_from = c(mean_rating_controllability, mean_rating_emotionalarousal, mean_rating_familiarity, mean_rating_negativity, mean_rating_positivity, mean_rating_priorexperience, mean_rating_vividness, sd_rating_controllability, sd_rating_emotionalarousal, sd_rating_familiarity, sd_rating_negativity, sd_rating_positivity, sd_rating_priorexperience, sd_rating_vividness), names_sep = "_")


mean_sd_event_level_ratings <- mean_sd_event_level_ratings %>%
  arrange(valence, lifeEvent) %>%
  select(valence, lifeEvent, mean_rating_controllability_Healthy, sd_rating_controllability_Healthy, mean_rating_emotionalarousal_Healthy, sd_rating_emotionalarousal_Healthy, mean_rating_familiarity_Healthy, sd_rating_familiarity_Healthy, mean_rating_negativity_Healthy, sd_rating_negativity_Healthy, mean_rating_positivity_Healthy, sd_rating_positivity_Healthy, mean_rating_priorexperience_Healthy, sd_rating_priorexperience_Healthy, mean_rating_vividness_Healthy, sd_rating_vividness_Healthy, mean_rating_controllability_Depression, sd_rating_controllability_Depression, mean_rating_emotionalarousal_Depression, sd_rating_emotionalarousal_Depression, mean_rating_familiarity_Depression, sd_rating_familiarity_Depression, mean_rating_negativity_Depression, sd_rating_negativity_Depression, mean_rating_positivity_Depression, sd_rating_positivity_Depression, mean_rating_priorexperience_Depression, sd_rating_priorexperience_Depression, mean_rating_vividness_Depression, sd_rating_vividness_Depression)

```


# Scepticism of Life Event Probabilities

```{r}
scepticism <- task_event_level_qs %>%
  select(subject, Group, debrief_factor, debrief) %>%
  group_by(subject) %>%
  distinct() %>% #Currently at trial level so multiple rows per participant - taking one row per participant for characteristics
  ungroup()

scepticism %>%
  tab_cols(Group) %>%
  tab_cells(debrief_factor) %>%
  tab_stat_cases() %>%
  tab_stat_cpct() %>%
  tab_pivot(stat_position = "inside_columns")%>%
  set_caption("The average probabilities in the task were accurate")

#Creating a binary variable of sceptical (somewhat disagree/strongly disagree) and non-sceptical (strongly agree, somewhat agree, neither agree nor disagree)
task_event_level_qs <- task_event_level_qs %>%
  mutate(scepticism_binary = ifelse(debrief < 4, 1, 2)) 

task_event_level_qs$scepticism_binary <- factor(task_event_level_qs$scepticism_binary,
                                       levels = c(1, 2),
                                       labels = c("Sceptical", "Non-Sceptical"))


#Aggregating scores from event level data
negative_eligible_aggregate_data_scepticism <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative" & !is.na(scepticism_binary)) %>%
  group_by(subject, scepticism_binary, Group, desirability) %>%
  summarise(update = mean(update, na.rm = TRUE)) %>%
  ungroup()

#Descriptive Statistics
negative_eligible_aggregate_data_scepticism %>%
  tab_cols(scepticism_binary, Group %nest% scepticism_binary %nest% desirability) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores for Negative Life Events by Scepticism, Group and Desirability")

#Bar Graph by Scepticism
negative_eligible_aggregate_data_scepticism %>%
  group_by(scepticism_binary) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = scepticism_binary, y = mean, fill = scepticism_binary)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Scepticism") + 
  labs(fill = "Scepticism") +
  ggtitle("Update Scores by Scepticism (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Scepticism & Group
negative_eligible_aggregate_data_scepticism %>%
  group_by(scepticism_binary, Group) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = scepticism_binary)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Scepticism") +
  ggtitle("Update Scores by Group and Scepticism (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

#Bar Graph by Scepticism, Group & desirability
negative_eligible_aggregate_data_scepticism %>%
  group_by(scepticism_binary, Group, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = scepticism_binary, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  facet_wrap(~ Group) +
  ylab("Update Scores") + 
  xlab ("Scepticism") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Group, Scepticism and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()

```
## Hypothesis 1

Removing individuals sceptical
```{r}
#Aggregating scores from event level data
negative_eligible_aggregate_data_scepticism_removed <- negative_eligible_aggregate_data_scepticism %>%
  filter(scepticism_binary == "Non-Sceptical")

#Descriptive Statistics
negative_eligible_aggregate_data_scepticism_removed %>%
  tab_cols(Group, desirability, Group %nest% desirability) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores for Negative Life Events by Group and Desirability with Sceptical Participants Removed")


#Regression Model
hyp1_model_sceptical <- lmer(update ~ Group*desirability + (1|subject), data = negative_eligible_aggregate_data_scepticism_removed)
summ(hyp1_model_sceptical, confint = TRUE, digits = 3)
hyp1_model_sceptical_effects <- ggpredict(hyp1_model_sceptical, terms = c("Group", "desirability"))
plot(hyp1_model_sceptical_effects) +
  labs(
    title = "Predicted Values of Update Scores for Negative Life Events",
    y = "Update Scores",
    x = "Group"
  )


standardize_parameters(hyp1_model_sceptical, method = "refit", ci = 0.95, include_response = TRUE)

```
## Hypothesis 2

```{r}
pos_neg_eligible_aggregate_data_scepticism <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & !is.na(scepticism_binary)) %>%
  group_by(subject, scepticism_binary, Group, valence, desirability) %>%
  summarise(update = mean(update, na.rm = TRUE)) %>%
  ungroup()

#Descriptive Statistics
pos_neg_eligible_aggregate_data_scepticism %>%
  tab_cols(Group %nest% scepticism_binary %nest% valence %nest% desirability) %>%
  tab_cells(update) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Update Scores for by Group, Scepticism, Valence and Desirability")

#Bar Graph by Scepticism, Group, valence & desirability
pos_neg_eligible_aggregate_data_scepticism %>%
  mutate(valence_labelled = ifelse(valence == "negative", "Negative", "Positive")) %>%
  group_by(scepticism_binary, Group, valence_labelled, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  facet_grid(scepticism_binary ~ valence_labelled, scales='free') +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  ggtitle("Update Scores by Group, Valence, Scepticism and Desirability (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal() +  scale_y_continuous(limits=c(0,15.5)) +  scale_x_discrete(labels = c("Healthy", "Depression"))+
 theme(
   strip.background = element_rect(
     color="black", fill="#FFFFFF", size=1.5, linetype="solid"
     )
   )

tiff(file="pos_neg_update_scepticism_group_desirability.tiff", units="cm", width=15.91, height=9.82, res=300)
pos_neg_eligible_aggregate_data_scepticism %>%
  mutate(valence_labelled = ifelse(valence == "negative", "Negative", "Positive")) %>%
  group_by(scepticism_binary, Group, valence_labelled, desirability) %>%
  summarise(mean = mean(update, na.rm = TRUE), n = n(), sd = sd(update, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = desirability)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  facet_grid(scepticism_binary ~ valence_labelled, scales='free') +
  ylab("Update Scores") + 
  xlab ("Group") + 
  labs(fill = "Desirability") +
  scale_fill_brewer(palette = "Blues") + theme_minimal() +  scale_y_continuous(limits=c(0,15.5)) +  scale_x_discrete(labels = c("Healthy", "Depression"))+
 theme(
   strip.background = element_rect(
     color="black", fill="#FFFFFF", size=1.5, linetype="solid"
     )
   )
dev.off()

#Regression Model
##Removing Sceptical
pos_neg_eligible_aggregate_data_scepticism_removed <- pos_neg_eligible_aggregate_data_scepticism %>%
  filter(scepticism_binary == "Non-Sceptical")

hyp2_model_scepticism_removed <- lmer(update ~ Group*valence*desirability + (1|subject), data = pos_neg_eligible_aggregate_data_scepticism_removed)
summ(hyp2_model_scepticism_removed, confint = TRUE, digits = 3)
hyp2_model_scepticism_removed_effects <- ggpredict(hyp2_model_scepticism_removed, terms = c("Group", "desirability", "valence"))
plot(hyp2_model_scepticism_removed_effects) +
  labs(
    title = "Predicted Values of Update Scores",
    y = "Update Scores",
    x = "Group"
  )

standardize_parameters(hyp2_model_scepticism_removed, method = "refit", ci = 0.95, include_response = TRUE)

```
## Hypothesis 3

```{r}
#Aggregating scores from event level data
negative_initial_estimates_scepticism <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group, scepticism_binary) %>%
  filter(!is.na(scepticism_binary)) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()

#Descriptive Statistics
negative_initial_estimates_scepticism %>%
  tab_cols(Group %nest% scepticism_binary) %>%
  tab_cells(p_estimate1) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Initial Estimates for Negative Life Events by Group")

#Bar Graph by Group & Scepticism
negative_initial_estimates_scepticism %>%
  group_by(Group, scepticism_binary) %>%
  summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = scepticism_binary)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  labs(fill = "Scepticism") +
  ggtitle("Initial Estimates by Group & Scepticism (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()



negative_initial_estimates_scepticism_removed <- negative_initial_estimates_scepticism %>%
  filter(scepticism_binary == "Non-Sceptical")
  

#Regression Model
hyp3_model_scepticism_removed <- lm(p_estimate1 ~ Group, data = negative_initial_estimates_scepticism_removed)
summ(hyp3_model_scepticism_removed, confint = TRUE, digits = 3)
hyp3_model_scepticism_removed_effects <- ggpredict(hyp3_model_scepticism_removed, terms = c("Group"))
plot(hyp3_model_scepticism_removed_effects) +
  labs(
    title = "Predicted Values of Initial Estimates for Negative Life Events",
    y = "Initial Estimates",
    x = "Group"
  )


standardize_parameters(hyp3_model_scepticism_removed, method = "refit", ci = 0.95, include_response = TRUE)

```

## Hypothesis 4

```{r}
#Aggregating scores from event level data
positive_initial_estimates_scepticism <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "positive") %>%
  group_by(subject, Group, scepticism_binary) %>%
  filter(!is.na(scepticism_binary)) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()

#Descriptive Statistics
positive_initial_estimates_scepticism %>%
  tab_cols(Group %nest% scepticism_binary) %>%
  tab_cells(p_estimate1) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Initial Estimates for positive Life Events by Group")

#Bar Graph by Group & Scepticism
positive_initial_estimates_scepticism %>%
  group_by(Group, scepticism_binary) %>%
  summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean, fill = scepticism_binary)) + 
  geom_bar(stat = 'identity', position = position_dodge()) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  labs(fill = "Scepticism") +
  ggtitle("Initial Estimates by Group & Scepticism (Mean, SE)") +
  scale_fill_brewer(palette = "Blues") + theme_minimal()


    
    
positive_initial_estimates_scepticism_removed <- positive_initial_estimates_scepticism %>%
  filter(scepticism_binary == "Non-Sceptical")
  

#Regression Model
hyp4_model_scepticism_removed <- lm(p_estimate1 ~ Group, data = positive_initial_estimates_scepticism_removed)
summ(hyp4_model_scepticism_removed, confint = TRUE, digits = 3)
hyp4_model_scepticism_removed_effects <- ggpredict(hyp4_model_scepticism_removed, terms = c("Group"))
plot(hyp4_model_scepticism_removed_effects) +
  labs(
    title = "Predicted Values of Initial Estimates for positive Life Events",
    y = "Initial Estimates",
    x = "Group"
  )


standardize_parameters(hyp4_model_scepticism_removed, method = "refit", ci = 0.95, include_response = TRUE)

```

### Positive & Negative Initial Estimates Graph by Scepticism
```{r}
#Aggregating scores from event level data
positive_negative_initial_estimates_scepticism <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, scepticism_binary, valence) %>%
  filter(!is.na(scepticism_binary)) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE)) %>% #Collapsing across
  ungroup()

#Bar Graph by Group & Scepticism
tiff(file="pos_neg_initial_estimates_scepticism_group_desirability.tiff", units="cm", width=15.91, height=9.82, res=300)
positive_negative_initial_estimates_scepticism %>%
    mutate(valence_labelled = ifelse(valence == "negative", "Negative", "Positive")) %>%
  group_by(Group, scepticism_binary, valence_labelled) %>%
  summarise(mean = mean(p_estimate1, na.rm = TRUE), n = n(), sd = sd(p_estimate1, na.rm = TRUE), se = (sd/sqrt(n))) %>%
  ggplot(., aes(x = Group, y = mean)) + 
  geom_bar(stat = 'identity', position = position_dodge(), fill = "#C7E9C0") +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(.9), width = .2) +
  facet_grid(scepticism_binary ~ valence_labelled) +
  ylab("Initial Estimates") + 
  xlab ("Group") + 
  labs(fill = "Scepticism") +
  scale_fill_brewer() +
 theme_minimal()+
 theme(
   strip.background = element_rect(
     color="black", fill="#FFFFFF", size=1.5, linetype="solid"
     )
   )+  scale_x_discrete(labels = c("Healthy", "Depression"))
dev.off()



```

# Number of Trials Descriptives

## Ineligible Trials

```{r}
number_ineligible_trials <- task_event_level_qs %>%
  filter(eligible_trial == "Missing Participant Estimate" | eligible_trial == "Estimation Error of Zero"  | eligible_trial == "Outside Range") %>%
  group_by(subject, Group, valence, eligible_trial) %>%
  summarise(n_trials = n()) %>%
  pivot_wider(names_from = eligible_trial, values_from = n_trials) %>%
  ungroup()

#Descriptive Statistics
number_ineligible_trials %>%
  tab_cols(valence %nest% Group) %>%
  tab_cells(`Estimation Error of Zero`, `Missing Participant Estimate`, `Outside Range`) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Ineligible Trials by Group and Valence")
```

## Eligible Trials

```{r}
number_eligible_trials <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability) %>%
  summarise(n_eligible_trials = n()) %>%
  ungroup()

#Descriptive Statistics
number_eligible_trials %>%
  tab_cols(valence %nest% Group %nest% desirability) %>%
  tab_cells(n_eligible_trials) %>%
  tab_stat_mean(label = "Mean") %>%
  tab_stat_sd(label = "SD") %>%
  tab_stat_se(label = "SE") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Eligible Trials by Group, Valence and Desirability")



```

# Exploratory analyses adjusting for age

## Hypothesis One

```{r}
##Aggregating scores from event level data
negative_eligible_aggregate_data_age <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group, desirability) %>%
  summarise(update = mean(update, na.rm = TRUE), age = mean(age, na.rm = TRUE)) %>%
  ungroup()

##Adjusted Regression
hyp1_model_age <- lmer(update ~ Group*desirability + age + (1|subject), data = negative_eligible_aggregate_data_age)
summ(hyp1_model_age, confint = TRUE, digits = 3)
standardize_parameters(hyp1_model_age, method = "refit", ci = 0.95, include_response = TRUE)
```

## Hypothesis Two

```{r}
positive_negative_eligible_aggregate_data_age <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible") %>%
  group_by(subject, Group, valence, desirability) %>%
  summarise(update = mean(update, na.rm = TRUE), age = mean(age, na.rm = TRUE)) %>%
  ungroup()

hyp2_model_age <- lmer(update ~ Group*valence*desirability + age + (1|subject), data = positive_negative_eligible_aggregate_data_age)
summ(hyp2_model_age, confint = TRUE, digits = 3)
standardize_parameters(hyp2_model_age, method = "refit", ci = 0.95, include_response = TRUE)
```

## Hypothesis Three

```{r}
negative_initial_estimates_age <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "negative") %>%
  group_by(subject, Group) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE), age = mean(age, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()

hyp3_model_age <- lm(p_estimate1 ~ Group + age, data = negative_initial_estimates_age)
summ(hyp3_model_age, confint = TRUE, digits = 3)
standardize_parameters(hyp3_model_age, method = "refit", ci = 0.95, include_response = TRUE)


```

## Hypothesis Four

```{r}
#Hypothesis Four
positive_initial_estimates_age <- task_event_level_qs %>%
  filter(eligible_trial == "Eligible" & valence == "positive") %>%
  group_by(subject, Group) %>%
  summarise(p_estimate1 = mean(p_estimate1, na.rm = TRUE), age = mean(age, na.rm = TRUE)) %>% #Collapsing across desirability
  ungroup()

hyp4_model_age <- lm(p_estimate1 ~ Group + age, data = positive_initial_estimates_age)
summ(hyp4_model_age, confint = TRUE, digits = 3)
standardize_parameters(hyp4_model_age, method = "refit", ci = 0.95, include_response = TRUE)
```


