---
title: Self processing in relation to emotion and reward processing in depression
  - Data Analysis
author: "Catherine Hobbs (c.hobbs@bath.ac.uk)"
output:
  html_document:
    df_print: paged
---

```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(expss)
library(data.table)
library(irr)
library(lmerTest)
library(jtools)
library(lmerTest)
library(ggeffects)
library(epiDisplay)
library(broom.mixed)
library(nlme)
library(emmeans)#Contrasts
library(effectsize) #Standardising parameters
expss_output_rnotebook()
expss_digits(digits = 2)
library(psych)#ICC calculation
```

```{r message=FALSE, warning=FALSE}
#Questionnaires
load("qs_anon.Rda")

#Associative Learning
load("associative_long_matching_anon.Rda")
load("associative_long_matching_collapsed_anon.Rda")
load("associative_wide_matching_anon.Rda")
load("associative_wide_matching_collapsed_anon.Rda")
load("associative_df_trial_anon.Rda")

#GNAT
load("GNAT_long_Qs_anon.Rda")
load("GNAT_wide_Qs_anon.Rda")
load("GNAT_trial_clean_Qs_anon.Rda")

#SEL
load("SEL_long_anon.Rda")
load("SEL_long_sep_rules_anon.Rda")
load("SEL_wide_anon.Rda")
load("SEL_trial_anon.Rda")
load("SEL_bias_anon.Rda")

```

# Sample Characteristics

```{r}

#Sample Characteristics 
qs %>%
  tab_cols(session %nest% PHQ_dep_group_3) %>%
  tab_cells(age) %>%
  tab_stat_mean_sd_n() %>%
  tab_cells(gender, ethnicity, employment_status, primary_depression_CIS, psychother_curr_, antids_curr_) %>%
  tab_stat_cases(label = "N", total_row_position = "none") %>%
  tab_stat_cpct(label = "%", total_row_position = "none") %>%
  tab_cells(PHQ_tot, BDI_tot, BFNE_tot, GAD_tot, FPE_tot, DAS24_tot, RSES_tot, PANAS_pos_change, PANAS_neg_change) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Participant Characteristics by Session and PHQ-9 Depression Severity")

#Test-Retest Reliability of PHQ and BDI
PHQ_ICC_df <- qs %>%
  dplyr:::select(subject, session, PHQ_tot) %>%
  spread(key = session, value = PHQ_tot) %>%
  rename(PHQ_1 = `1`, PHQ_2 = `2`) %>%
  dplyr::select(PHQ_1, PHQ_2)
  
icc(PHQ_ICC_df, model = "twoway", unit = "average", type = "agreement")

BDI_ICC_df <- qs %>%
  dplyr::select(subject, session, BDI_tot) %>%
  spread(key = session, value = BDI_tot) %>%
  rename(BDI_1 = `1`, BDI_2 = `2`) %>%
  dplyr::select(BDI_1, BDI_2)

icc(BDI_ICC_df, model = "twoway", unit = "average", type = "agreement")

#Correlation between BDI and PHQ
cor.test(qs$PHQ_tot, qs$BDI_tot)


```
# Associative Learning

## Do Accuracy/Reaction Times differ according to stimuli?

### Descriptive Statistics for all associative learning tasks
```{r}
associative_long_matching_collapsed %>%
  tab_cols(stimuli) %>%
  tab_cells(prop_acc) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("Associative Learning Self Accuracy by Stimuli")

associative_long_matching_collapsed %>%
  group_by(Task, stimuli) %>%
  summarize(Mean = mean(prop_acc),  SD = sd(prop_acc), Minimum = min(prop_acc), Maximum = max(prop_acc)) %>%
  ggplot(., aes(x = stimuli, y = Mean, fill = Task)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  coord_cartesian(ylim=c(50, 100)) +
  ylab("Accuracy (%)") + 
  xlab ("Stimuli") + 
  labs(fill = "Stimuli") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))+
  scale_fill_brewer()

associative_long_matching_collapsed %>%
  tab_cols(stimuli) %>%
  tab_cells(mean_rt_mult) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("Associative Learning Self Reaction Times by Stimuli")

associative_long_matching_collapsed %>%
  group_by(Task, stimuli) %>%
  summarize(Mean = mean(mean_rt_mult),  SD = sd(mean_rt_mult), Minimum = min(mean_rt_mult), Maximum = max(mean_rt_mult)) %>%
  ggplot(., aes(x = stimuli, y = Mean, fill = Task)) + 
  coord_cartesian(ylim=c(600, 900)) +
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Reaction Times (ms)") + 
  xlab ("Stimuli") + 
  labs(fill = "Stimuli") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))+
  scale_fill_brewer()
```

### Main effect analysis
```{r}
#Self
Assoc_self <- associative_long_matching_collapsed %>%
  filter(Task == "Self") %>%
  mutate(prop_acc_standard = (prop_acc-mean(prop_acc))/sd(prop_acc))

##Accuracy
assoc_self_acc_main_model <- lmer(prop_acc ~ stimuli + (1|subject), data = Assoc_self)
jtools::summ(assoc_self_acc_main_model, confint = TRUE, digits = 3)

###Standardised
assoc_self_acc_main_model_standard <- standardize_parameters(assoc_self_acc_main_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_self_acc_main_model_standard

##Reaction Times
assoc_self_rt_main_model <- lmer(mean_rt_mult ~ stimuli + (1|subject), data = Assoc_self)
jtools::summ(assoc_self_rt_main_model, confint = TRUE, digits = 3)

###Standardised
assoc_self_rt_main_model_standard <- standardize_parameters(assoc_self_rt_main_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_self_rt_main_model_standard

#Reward
Assoc_reward <- associative_long_matching_collapsed %>%
  filter(Task == "Reward")

##Accuracy
assoc_reward_acc_main_model <- lmer(prop_acc ~ stimuli + (1|subject), data = Assoc_reward)
jtools::summ(assoc_reward_acc_main_model, confint = TRUE, digits = 3)

###Standardised
assoc_reward_acc_main_model_standard <- standardize_parameters(assoc_reward_acc_main_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_reward_acc_main_model_standard

##Reaction Times
assoc_reward_rt_main_model <- lmer(mean_rt_mult ~ stimuli + (1|subject), data = Assoc_reward)
jtools::summ(assoc_reward_rt_main_model, confint = TRUE, digits = 3)

###Standardised
assoc_reward_rt_main_model_standard <- standardize_parameters(assoc_reward_rt_main_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_reward_rt_main_model_standard

#Valence
Assoc_valence <- associative_long_matching_collapsed %>%
  filter(Task == "Valence")

##Accuracy
assoc_valence_acc_main_model <- lmer(prop_acc ~ stimuli + (1|subject), data = Assoc_valence)
jtools::summ(assoc_valence_acc_main_model, confint = TRUE, digits = 3)

###Standardised
assoc_valence_acc_main_model_standard <- standardize_parameters(assoc_valence_acc_main_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_valence_acc_main_model_standard

##Reaction Times
assoc_valence_rt_main_model <- lmer(mean_rt_mult ~ stimuli + (1|subject), data = Assoc_valence)
jtools::summ(assoc_valence_rt_main_model, confint = TRUE, digits = 3)

###Standardised
assoc_valence_rt_main_model_standard <- standardize_parameters(assoc_valence_rt_main_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_valence_rt_main_model_standard
```



## Are accuracy/reaction times for each stimuli associated with depression severity as measured by the PHQ-9?

### Descriptive Statistics for all associative tasks

```{r}
associative_long_matching_collapsed %>%
  tab_cols(stimuli) %>%
  tab_rows(PHQ_dep_group_3) %>%
  tab_cells(prop_acc) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot() %>%
  set_caption("Associative Learning Self Accuracy (%) by Stimuli and PHQ-9 Severity")

ggplot(data = associative_long_matching_collapsed, aes(x = prop_acc, y = PHQ_tot, color = stimuli)) + facet_wrap(~ Task) + geom_smooth(method = lm, fullrange = TRUE) + ylab("PHQ-9 Total Scores") + xlab("Accuracy (%)") + ggtitle("Relationship between PHQ-9 Scores and Accuracy")

associative_long_matching_collapsed %>%
  tab_cols(stimuli) %>%
  tab_rows(PHQ_dep_group_3) %>%
  tab_cells(mean_rt_mult) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot() %>%
  set_caption("Associative Learning Self Reaction Times (ms) by Stimuli and PHQ-9 Severity")

ggplot(data = associative_long_matching_collapsed, aes(x = mean_rt_mult, y = PHQ_tot, color = stimuli)) + facet_wrap(~ Task) + geom_smooth(method = lm, fullrange = TRUE) + ylab("PHQ-9 Total Scores") + xlab("Reaction Times (ms)") + ggtitle("Relationship between PHQ-9 Scores and Reaction Times")


```

### Self

#### Accuracy
```{r}
associative_long_matching_collapsed_self <- associative_long_matching_collapsed %>%
  filter(Task == "Self")

#Descriptives Graphs
ggplot(data = associative_long_matching_collapsed_self, aes(x = prop_acc, y = PHQ_tot, color = stimuli)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("PHQ-9 Total Scores") + xlab("Accuracy (%)") + ggtitle("Self Task Linear Relationship between PHQ-9 Scores and Accuracy")

#Regression Model
associative_wide_matching_collapsed_self <- associative_wide_matching_collapsed %>%
  filter(Task == "Self") #Restricting Task to "Self"

assoc_self_acc_PHQ_model <- lm(PHQ_tot ~ prop_acc_High + prop_acc_Medium + prop_acc_Low, data = associative_wide_matching_collapsed_self)
jtools::summ(assoc_self_acc_PHQ_model, confint = TRUE, digits = 3)


#Standardised coefficients
assoc_self_acc_PHQ_model_standard <- standardize_parameters(assoc_self_acc_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_self_acc_PHQ_model_standard

#Regression Graphs
assoc_self_acc_PHQ_model_effects <- ggpredict(assoc_self_acc_PHQ_model)

prop_acc_High_effects <- as.data.frame(assoc_self_acc_PHQ_model_effects$prop_acc_High)
prop_acc_Medium_effects <- as.data.frame(assoc_self_acc_PHQ_model_effects$prop_acc_Medium)
prop_acc_Low_effects <- as.data.frame(assoc_self_acc_PHQ_model_effects$prop_acc_Low)

associative_wide_matching_collapsed_self_effects_graph <- bind_rows(prop_acc_High_effects, prop_acc_Medium_effects, prop_acc_Low_effects)

associative_wide_matching_collapsed_self_effects_graph <- associative_wide_matching_collapsed_self_effects_graph %>%
  mutate(Stimuli = ifelse(group == "prop_acc_High", 1, ifelse(group == "prop_acc_Medium", 2, 3)))

associative_wide_matching_collapsed_self_effects_graph$Stimuli <- factor(associative_wide_matching_collapsed_self_effects_graph$Stimuli,
                                       levels = c(1, 2, 3),
                                       labels = c("Self", "Friend", "Stranger"))

ggplot(associative_wide_matching_collapsed_self_effects_graph, aes(x = x, y = predicted, color = Stimuli)) + geom_line(size = 1) + geom_ribbon(data = associative_wide_matching_collapsed_self_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Accuracy (%)") + ylab("PHQ-9 Total") + ggtitle("Linear Regression Self Task Predicted Effects of Accuracy (%) on PHQ-9 Scores")
```
#### Reaction Times
```{r}
#Descriptives Graphs
ggplot(data = associative_long_matching_collapsed_self, aes(x = mean_rt_mult, y = PHQ_tot, color = stimuli)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("PHQ-9 Total Scores") + xlab("Reaction Times (ms)") + ggtitle("Self Task Linear Relationship between PHQ-9 Scores and Accuracy")

#Regression Model
associative_wide_matching_collapsed_self <- associative_wide_matching_collapsed %>%
  filter(Task == "Self") #Restricting Task to "Self"

assoc_self_rt_PHQ_model <- lm(PHQ_tot ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low, data = associative_wide_matching_collapsed_self)
jtools::summ(assoc_self_rt_PHQ_model, confint = TRUE, digits = 3)

#Standardised coefficients
assoc_self_rt_PHQ_model_standard <- standardize_parameters(assoc_self_rt_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_self_rt_PHQ_model_standard

#Regression Graphs
assoc_self_rt_PHQ_model_effects <- ggpredict(assoc_self_rt_PHQ_model)

mean_rt_mult_High_effects <- as.data.frame(assoc_self_rt_PHQ_model_effects$mean_rt_mult_High)
mean_rt_mult_Medium_effects <- as.data.frame(assoc_self_rt_PHQ_model_effects$mean_rt_mult_Medium)
mean_rt_mult_Low_effects <- as.data.frame(assoc_self_rt_PHQ_model_effects$mean_rt_mult_Low)

associative_wide_matching_collapsed_self_effects_graph <- bind_rows(mean_rt_mult_High_effects, mean_rt_mult_Medium_effects, mean_rt_mult_Low_effects)

associative_wide_matching_collapsed_self_effects_graph <- associative_wide_matching_collapsed_self_effects_graph %>%
  mutate(Stimuli = ifelse(group == "mean_rt_mult_High", 1, ifelse(group == "mean_rt_mult_Medium", 2, 3)))

associative_wide_matching_collapsed_self_effects_graph$Stimuli <- factor(associative_wide_matching_collapsed_self_effects_graph$Stimuli,
                                       levels = c(1, 2, 3),
                                       labels = c("Self", "Friend", "Stranger"))

ggplot(associative_wide_matching_collapsed_self_effects_graph, aes(x = x, y = predicted, color = Stimuli)) + geom_line(size = 1) + geom_ribbon(data = associative_wide_matching_collapsed_self_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Reaction Times (ms)") + ylab("PHQ-9 Total") + ggtitle("Linear Regression Self Task Predicted Effects of Reaction Times (ms) on PHQ-9 Scores")
```
### Reward

#### Accuracy
```{r}
associative_long_matching_collapsed_reward <- associative_long_matching_collapsed %>%
  filter(Task == "Reward")

#Descriptives Graphs
ggplot(data = associative_long_matching_collapsed_reward, aes(x = prop_acc, y = PHQ_tot, color = stimuli)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("PHQ-9 Total Scores") + xlab("Accuracy (%)") + ggtitle("Reward Task Linear Relationship between PHQ-9 Scores and Accuracy")

#Regression Model
associative_wide_matching_collapsed_reward <- associative_wide_matching_collapsed %>%
  filter(Task == "Reward") #Restricting Task to "Reward"

assoc_reward_acc_PHQ_model <- lm(PHQ_tot ~ prop_acc_High + prop_acc_Medium + prop_acc_Low, data = associative_wide_matching_collapsed_reward)
jtools::summ(assoc_reward_acc_PHQ_model, confint = TRUE, digits = 3)

#Standardised coefficients
assoc_reward_acc_PHQ_model_standard <- standardize_parameters(assoc_reward_acc_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_reward_acc_PHQ_model_standard

#Regression Graphs
assoc_reward_acc_PHQ_model_effects <- ggpredict(assoc_reward_acc_PHQ_model)

prop_acc_High_effects <- as.data.frame(assoc_reward_acc_PHQ_model_effects$prop_acc_High)
prop_acc_Medium_effects <- as.data.frame(assoc_reward_acc_PHQ_model_effects$prop_acc_Medium)
prop_acc_Low_effects <- as.data.frame(assoc_reward_acc_PHQ_model_effects$prop_acc_Low)

associative_wide_matching_collapsed_reward_effects_graph <- bind_rows(prop_acc_High_effects, prop_acc_Medium_effects, prop_acc_Low_effects)

associative_wide_matching_collapsed_reward_effects_graph <- associative_wide_matching_collapsed_reward_effects_graph %>%
  mutate(Stimuli = ifelse(group == "prop_acc_High", 1, ifelse(group == "prop_acc_Medium", 2, 3)))

associative_wide_matching_collapsed_reward_effects_graph$Stimuli <- factor(associative_wide_matching_collapsed_reward_effects_graph$Stimuli,
                                       levels = c(1, 2, 3),
                                       labels = c("High (£9)", "Medium (£3)", "Low (£1)"))

ggplot(associative_wide_matching_collapsed_reward_effects_graph, aes(x = x, y = predicted, color = Stimuli)) + geom_line(size = 1) + geom_ribbon(data = associative_wide_matching_collapsed_reward_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Accuracy (%)") + ylab("PHQ-9 Total") + ggtitle("Linear Regression Reward Task Predicted Effects of Accuracy (%) on PHQ-9 Scores")
```
#### Reaction Times
```{r}
#Descriptives Graphs
ggplot(data = associative_long_matching_collapsed_reward, aes(x = mean_rt_mult, y = PHQ_tot, color = stimuli)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("PHQ-9 Total Scores") + xlab("Reaction Times (ms)") + ggtitle("Reward Task Linear Relationship between PHQ-9 Scores and Accuracy")

#Regression Model
associative_wide_matching_collapsed_reward <- associative_wide_matching_collapsed %>%
  filter(Task == "Reward") #Restricting Task to "Reward"

assoc_reward_rt_PHQ_model <- lm(PHQ_tot ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low, data = associative_wide_matching_collapsed_reward)
jtools::summ(assoc_reward_rt_PHQ_model, confint = TRUE, digits = 3)

#Standardised coefficients
assoc_reward_rt_PHQ_model_standard <- standardize_parameters(assoc_reward_rt_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_reward_rt_PHQ_model_standard

#Regression Graphs
assoc_reward_rt_PHQ_model_effects <- ggpredict(assoc_reward_rt_PHQ_model)

mean_rt_mult_High_effects <- as.data.frame(assoc_reward_rt_PHQ_model_effects$mean_rt_mult_High)
mean_rt_mult_Medium_effects <- as.data.frame(assoc_reward_rt_PHQ_model_effects$mean_rt_mult_Medium)
mean_rt_mult_Low_effects <- as.data.frame(assoc_reward_rt_PHQ_model_effects$mean_rt_mult_Low)

associative_wide_matching_collapsed_reward_effects_graph <- bind_rows(mean_rt_mult_High_effects, mean_rt_mult_Medium_effects, mean_rt_mult_Low_effects)

associative_wide_matching_collapsed_reward_effects_graph <- associative_wide_matching_collapsed_reward_effects_graph %>%
  mutate(Stimuli = ifelse(group == "mean_rt_mult_High", 1, ifelse(group == "mean_rt_mult_Medium", 2, 3)))

associative_wide_matching_collapsed_reward_effects_graph$Stimuli <- factor(associative_wide_matching_collapsed_reward_effects_graph$Stimuli,
                                       levels = c(1, 2, 3),
                                       labels = c("High (£9)", "Medium (£3)", "Low (£1)"))

ggplot(associative_wide_matching_collapsed_reward_effects_graph, aes(x = x, y = predicted, color = Stimuli)) + geom_line(size = 1) + geom_ribbon(data = associative_wide_matching_collapsed_reward_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Reaction Times (ms)") + ylab("PHQ-9 Total") + ggtitle("Linear Regression Reward Task Predicted Effects of Reaction Times (ms) on PHQ-9 Scores")
```

### Valence

#### Accuracy
```{r}
associative_long_matching_collapsed_valence <- associative_long_matching_collapsed %>%
  filter(Task == "Valence")

#Descriptives Graphs
ggplot(data = associative_long_matching_collapsed_valence, aes(x = prop_acc, y = PHQ_tot, color = stimuli)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("PHQ-9 Total Scores") + xlab("Accuracy (%)") + ggtitle("Valence Task Linear Relationship between PHQ-9 Scores and Accuracy")

#Regression Model
associative_wide_matching_collapsed_valence <- associative_wide_matching_collapsed %>%
  filter(Task == "Valence") #Restricting Task to "Valence"

assoc_valence_acc_PHQ_model <- lm(PHQ_tot ~ prop_acc_High + prop_acc_Medium + prop_acc_Low, data = associative_wide_matching_collapsed_valence)
jtools::summ(assoc_valence_acc_PHQ_model, confint = TRUE, digits = 3)

#Standardised coefficients
assoc_valence_acc_PHQ_model_standard <- standardize_parameters(assoc_valence_acc_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_valence_acc_PHQ_model_standard

#Regression Graphs
assoc_valence_acc_PHQ_model_effects <- ggpredict(assoc_valence_acc_PHQ_model)

prop_acc_High_effects <- as.data.frame(assoc_valence_acc_PHQ_model_effects$prop_acc_High)
prop_acc_Medium_effects <- as.data.frame(assoc_valence_acc_PHQ_model_effects$prop_acc_Medium)
prop_acc_Low_effects <- as.data.frame(assoc_valence_acc_PHQ_model_effects$prop_acc_Low)

associative_wide_matching_collapsed_valence_effects_graph <- bind_rows(prop_acc_High_effects, prop_acc_Medium_effects, prop_acc_Low_effects)

associative_wide_matching_collapsed_valence_effects_graph <- associative_wide_matching_collapsed_valence_effects_graph %>%
  mutate(Stimuli = ifelse(group == "prop_acc_High", 1, ifelse(group == "prop_acc_Medium", 2, 3)))

associative_wide_matching_collapsed_valence_effects_graph$Stimuli <- factor(associative_wide_matching_collapsed_valence_effects_graph$Stimuli,
                                       levels = c(1, 2, 3),
                                       labels = c("Happy", "Neutral", "Sad"))

ggplot(associative_wide_matching_collapsed_valence_effects_graph, aes(x = x, y = predicted, color = Stimuli)) + geom_line(size = 1) + geom_ribbon(data = associative_wide_matching_collapsed_valence_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Accuracy (%)") + ylab("PHQ-9 Total") + ggtitle("Linear Regression Valence Task Predicted Effects of Accuracy (%) on PHQ-9 Scores")
```
#### Reaction Times
```{r}
#Descriptives Graphs
ggplot(data = associative_long_matching_collapsed_valence, aes(x = mean_rt_mult, y = PHQ_tot, color = stimuli)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("PHQ-9 Total Scores") + xlab("Reaction Times (ms)") + ggtitle("Valence Task Linear Relationship between PHQ-9 Scores and Accuracy")

#Regression Model
associative_wide_matching_collapsed_valence <- associative_wide_matching_collapsed %>%
  filter(Task == "Valence") #Restricting Task to "Valence"

assoc_valence_rt_PHQ_model <- lm(PHQ_tot ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low, data = associative_wide_matching_collapsed_valence)
jtools::summ(assoc_valence_rt_PHQ_model, confint = TRUE, digits = 3)

#Standardised coefficients
assoc_valence_rt_PHQ_model_standard <- standardize_parameters(assoc_valence_rt_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_valence_rt_PHQ_model_standard

#Regression Graphs
assoc_valence_rt_PHQ_model_effects <- ggpredict(assoc_valence_rt_PHQ_model)

mean_rt_mult_High_effects <- as.data.frame(assoc_valence_rt_PHQ_model_effects$mean_rt_mult_High)
mean_rt_mult_Medium_effects <- as.data.frame(assoc_valence_rt_PHQ_model_effects$mean_rt_mult_Medium)
mean_rt_mult_Low_effects <- as.data.frame(assoc_valence_rt_PHQ_model_effects$mean_rt_mult_Low)

associative_wide_matching_collapsed_valence_effects_graph <- bind_rows(mean_rt_mult_High_effects, mean_rt_mult_Medium_effects, mean_rt_mult_Low_effects)

associative_wide_matching_collapsed_valence_effects_graph <- associative_wide_matching_collapsed_valence_effects_graph %>%
  mutate(Stimuli = ifelse(group == "mean_rt_mult_High", 1, ifelse(group == "mean_rt_mult_Medium", 2, 3)))

associative_wide_matching_collapsed_valence_effects_graph$Stimuli <- factor(associative_wide_matching_collapsed_valence_effects_graph$Stimuli,
                                       levels = c(1, 2, 3),
                                       labels = c("Happy", "Neutral", "Sad"))

ggplot(associative_wide_matching_collapsed_valence_effects_graph, aes(x = x, y = predicted, color = Stimuli)) + geom_line(size = 1) + geom_ribbon(data = associative_wide_matching_collapsed_valence_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Reaction Times (ms)") + ylab("PHQ-9 Total") + ggtitle("Linear Regression Valence Task Predicted Effects of Reaction Times (ms) on PHQ-9 Scores")
```

## Are accuracy/reaction times for each stimuli associated with depression severity as measured by the BDI-II?

### Descriptive Statistics for all associative tasks

```{r}
associative_long_matching_collapsed %>%
  tab_cols(stimuli) %>%
  tab_rows(BDI_dep_group_3) %>%
  tab_cells(prop_acc) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot() %>%
  set_caption("Associative Learning Self Accuracy (%) by Stimuli and BDI-I2I Severity")

ggplot(data = associative_long_matching_collapsed, aes(x = prop_acc, y = BDI_tot, color = stimuli)) + facet_wrap(~ Task) + geom_smooth(method = lm, fullrange = TRUE) + ylab("BDI-II Total Scores") + xlab("Accuracy (%)") + ggtitle("Relationship between BDI-II Scores and Accuracy")

associative_long_matching_collapsed %>%
  tab_cols(stimuli) %>%
  tab_rows(BDI_dep_group_3) %>%
  tab_cells(mean_rt_mult) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot() %>%
  set_caption("Associative Learning Self Reaction Times (ms) by Stimuli and BDI-II Severity")

ggplot(data = associative_long_matching_collapsed, aes(x = mean_rt_mult, y = BDI_tot, color = stimuli)) + facet_wrap(~ Task) + geom_smooth(method = lm, fullrange = TRUE) + ylab("BDI-II Total Scores") + xlab("Reaction Times (ms)") + ggtitle("Relationship between BDI-II Scores and Reaction Times")


```

### Self

#### Accuracy
```{r}
#Descriptives Graphs
ggplot(data = associative_long_matching_collapsed_self, aes(x = prop_acc, y = BDI_tot, color = stimuli)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("BDI-II Total Scores") + xlab("Accuracy (%)") + ggtitle("Self Task Linear Relationship between BDI-II Scores and Accuracy")

#Regression Model
associative_wide_matching_collapsed_self <- associative_wide_matching_collapsed %>%
  filter(Task == "Self") #Restricting Task to "Self"

assoc_self_acc_BDI_model <- lm(BDI_tot ~ prop_acc_High + prop_acc_Medium + prop_acc_Low, data = associative_wide_matching_collapsed_self)
jtools::summ(assoc_self_acc_BDI_model, confint = TRUE, digits = 3)


#Standardised coefficients
assoc_self_acc_BDI_model_standard <- standardize_parameters(assoc_self_acc_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_self_acc_BDI_model_standard

#Regression Graphs
assoc_self_acc_BDI_model_effects <- ggpredict(assoc_self_acc_BDI_model)

prop_acc_High_effects <- as.data.frame(assoc_self_acc_BDI_model_effects$prop_acc_High)
prop_acc_Medium_effects <- as.data.frame(assoc_self_acc_BDI_model_effects$prop_acc_Medium)
prop_acc_Low_effects <- as.data.frame(assoc_self_acc_BDI_model_effects$prop_acc_Low)

associative_wide_matching_collapsed_self_effects_graph <- bind_rows(prop_acc_High_effects, prop_acc_Medium_effects, prop_acc_Low_effects)

associative_wide_matching_collapsed_self_effects_graph <- associative_wide_matching_collapsed_self_effects_graph %>%
  mutate(Stimuli = ifelse(group == "prop_acc_High", 1, ifelse(group == "prop_acc_Medium", 2, 3)))

associative_wide_matching_collapsed_self_effects_graph$Stimuli <- factor(associative_wide_matching_collapsed_self_effects_graph$Stimuli,
                                       levels = c(1, 2, 3),
                                       labels = c("Self", "Friend", "Stranger"))

ggplot(associative_wide_matching_collapsed_self_effects_graph, aes(x = x, y = predicted, color = Stimuli)) + geom_line(size = 1) + geom_ribbon(data = associative_wide_matching_collapsed_self_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Accuracy (%)") + ylab("BDI-II Total") + ggtitle("Linear Regression Self Task Predicted Effects of Accuracy (%) on BDI-II Scores")
```
#### Reaction Times
```{r}
#Descriptives Graphs
ggplot(data = associative_long_matching_collapsed_self, aes(x = mean_rt_mult, y = BDI_tot, color = stimuli)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("BDI-II Total Scores") + xlab("Reaction Times (ms)") + ggtitle("Self Task Linear Relationship between BDI-II Scores and Accuracy")

#Regression Model
associative_wide_matching_collapsed_self <- associative_wide_matching_collapsed %>%
  filter(Task == "Self") #Restricting Task to "Self"

assoc_self_rt_BDI_model <- lm(BDI_tot ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low, data = associative_wide_matching_collapsed_self)
jtools::summ(assoc_self_rt_BDI_model, confint = TRUE, digits = 3)

#Standardised coefficients
assoc_self_rt_BDI_model_standard <- standardize_parameters(assoc_self_rt_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_self_rt_BDI_model_standard

#Regression Graphs
assoc_self_rt_BDI_model_effects <- ggpredict(assoc_self_rt_BDI_model)

mean_rt_mult_High_effects <- as.data.frame(assoc_self_rt_BDI_model_effects$mean_rt_mult_High)
mean_rt_mult_Medium_effects <- as.data.frame(assoc_self_rt_BDI_model_effects$mean_rt_mult_Medium)
mean_rt_mult_Low_effects <- as.data.frame(assoc_self_rt_BDI_model_effects$mean_rt_mult_Low)

associative_wide_matching_collapsed_self_effects_graph <- bind_rows(mean_rt_mult_High_effects, mean_rt_mult_Medium_effects, mean_rt_mult_Low_effects)

associative_wide_matching_collapsed_self_effects_graph <- associative_wide_matching_collapsed_self_effects_graph %>%
  mutate(Stimuli = ifelse(group == "mean_rt_mult_High", 1, ifelse(group == "mean_rt_mult_Medium", 2, 3)))

associative_wide_matching_collapsed_self_effects_graph$Stimuli <- factor(associative_wide_matching_collapsed_self_effects_graph$Stimuli,
                                       levels = c(1, 2, 3),
                                       labels = c("Self", "Friend", "Stranger"))

ggplot(associative_wide_matching_collapsed_self_effects_graph, aes(x = x, y = predicted, color = Stimuli)) + geom_line(size = 1) + geom_ribbon(data = associative_wide_matching_collapsed_self_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Reaction Times (ms)") + ylab("BDI-II Total") + ggtitle("Linear Regression Self Task Predicted Effects of Reaction Times (ms) on BDI-II Scores")
```
### Reward

#### Accuracy
```{r}
#Descriptives Graphs
ggplot(data = associative_long_matching_collapsed_reward, aes(x = prop_acc, y = BDI_tot, color = stimuli)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("BDI-II Total Scores") + xlab("Accuracy (%)") + ggtitle("Reward Task Linear Relationship between BDI-II Scores and Accuracy")

#Regression Model
associative_wide_matching_collapsed_reward <- associative_wide_matching_collapsed %>%
  filter(Task == "Reward") #Restricting Task to "Reward"

assoc_reward_acc_BDI_model <- lm(BDI_tot ~ prop_acc_High + prop_acc_Medium + prop_acc_Low, data = associative_wide_matching_collapsed_reward)
jtools::summ(assoc_reward_acc_BDI_model, confint = TRUE, digits = 3)

#Standardised coefficients
assoc_reward_acc_BDI_model_standard <- standardize_parameters(assoc_reward_acc_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_reward_acc_BDI_model_standard

#Regression Graphs
assoc_reward_acc_BDI_model_effects <- ggpredict(assoc_reward_acc_BDI_model)

prop_acc_High_effects <- as.data.frame(assoc_reward_acc_BDI_model_effects$prop_acc_High)
prop_acc_Medium_effects <- as.data.frame(assoc_reward_acc_BDI_model_effects$prop_acc_Medium)
prop_acc_Low_effects <- as.data.frame(assoc_reward_acc_BDI_model_effects$prop_acc_Low)

associative_wide_matching_collapsed_reward_effects_graph <- bind_rows(prop_acc_High_effects, prop_acc_Medium_effects, prop_acc_Low_effects)

associative_wide_matching_collapsed_reward_effects_graph <- associative_wide_matching_collapsed_reward_effects_graph %>%
  mutate(Stimuli = ifelse(group == "prop_acc_High", 1, ifelse(group == "prop_acc_Medium", 2, 3)))

associative_wide_matching_collapsed_reward_effects_graph$Stimuli <- factor(associative_wide_matching_collapsed_reward_effects_graph$Stimuli,
                                       levels = c(1, 2, 3),
                                       labels = c("High (£9)", "Medium (£3)", "Low (£1)"))

ggplot(associative_wide_matching_collapsed_reward_effects_graph, aes(x = x, y = predicted, color = Stimuli)) + geom_line(size = 1) + geom_ribbon(data = associative_wide_matching_collapsed_reward_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Accuracy (%)") + ylab("BDI-II Total") + ggtitle("Linear Regression Reward Task Predicted Effects of Accuracy (%) on BDI-II Scores")
```
#### Reaction Times
```{r}
#Descriptives Graphs
ggplot(data = associative_long_matching_collapsed_reward, aes(x = mean_rt_mult, y = BDI_tot, color = stimuli)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("BDI-II Total Scores") + xlab("Reaction Times (ms)") + ggtitle("Reward Task Linear Relationship between BDI-II Scores and Accuracy")

#Regression Model
associative_wide_matching_collapsed_reward <- associative_wide_matching_collapsed %>%
  filter(Task == "Reward") #Restricting Task to "Reward"

assoc_reward_rt_BDI_model <- lm(BDI_tot ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low, data = associative_wide_matching_collapsed_reward)
jtools::summ(assoc_reward_rt_BDI_model, confint = TRUE, digits = 3)

#Standardised coefficients
assoc_reward_rt_BDI_model_standard <- standardize_parameters(assoc_reward_rt_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_reward_rt_BDI_model_standard

#Regression Graphs
assoc_reward_rt_BDI_model_effects <- ggpredict(assoc_reward_rt_BDI_model)

mean_rt_mult_High_effects <- as.data.frame(assoc_reward_rt_BDI_model_effects$mean_rt_mult_High)
mean_rt_mult_Medium_effects <- as.data.frame(assoc_reward_rt_BDI_model_effects$mean_rt_mult_Medium)
mean_rt_mult_Low_effects <- as.data.frame(assoc_reward_rt_BDI_model_effects$mean_rt_mult_Low)

associative_wide_matching_collapsed_reward_effects_graph <- bind_rows(mean_rt_mult_High_effects, mean_rt_mult_Medium_effects, mean_rt_mult_Low_effects)

associative_wide_matching_collapsed_reward_effects_graph <- associative_wide_matching_collapsed_reward_effects_graph %>%
  mutate(Stimuli = ifelse(group == "mean_rt_mult_High", 1, ifelse(group == "mean_rt_mult_Medium", 2, 3)))

associative_wide_matching_collapsed_reward_effects_graph$Stimuli <- factor(associative_wide_matching_collapsed_reward_effects_graph$Stimuli,
                                       levels = c(1, 2, 3),
                                       labels = c("High (£9)", "Medium (£3)", "Low (£1)"))

ggplot(associative_wide_matching_collapsed_reward_effects_graph, aes(x = x, y = predicted, color = Stimuli)) + geom_line(size = 1) + geom_ribbon(data = associative_wide_matching_collapsed_reward_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Reaction Times (ms)") + ylab("BDI-II Total") + ggtitle("Linear Regression Reward Task Predicted Effects of Reaction Times (ms) on BDI-II Scores")
```

### Valence

#### Accuracy
```{r}
#Descriptives Graphs
ggplot(data = associative_long_matching_collapsed_valence, aes(x = prop_acc, y = BDI_tot, color = stimuli)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("BDI-II Total Scores") + xlab("Accuracy (%)") + ggtitle("Valence Task Linear Relationship between BDI-II Scores and Accuracy")

#Regression Model
associative_wide_matching_collapsed_valence <- associative_wide_matching_collapsed %>%
  filter(Task == "Valence") #Restricting Task to "Valence"

assoc_valence_acc_BDI_model <- lm(BDI_tot ~ prop_acc_High + prop_acc_Medium + prop_acc_Low, data = associative_wide_matching_collapsed_valence)
jtools::summ(assoc_valence_acc_BDI_model, confint = TRUE, digits = 3)

#Standardised coefficients
assoc_valence_acc_BDI_model_standard <- standardize_parameters(assoc_valence_acc_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_valence_acc_BDI_model_standard

#Regression Graphs
assoc_valence_acc_BDI_model_effects <- ggpredict(assoc_valence_acc_BDI_model)

prop_acc_High_effects <- as.data.frame(assoc_valence_acc_BDI_model_effects$prop_acc_High)
prop_acc_Medium_effects <- as.data.frame(assoc_valence_acc_BDI_model_effects$prop_acc_Medium)
prop_acc_Low_effects <- as.data.frame(assoc_valence_acc_BDI_model_effects$prop_acc_Low)

associative_wide_matching_collapsed_valence_effects_graph <- bind_rows(prop_acc_High_effects, prop_acc_Medium_effects, prop_acc_Low_effects)

associative_wide_matching_collapsed_valence_effects_graph <- associative_wide_matching_collapsed_valence_effects_graph %>%
  mutate(Stimuli = ifelse(group == "prop_acc_High", 1, ifelse(group == "prop_acc_Medium", 2, 3)))

associative_wide_matching_collapsed_valence_effects_graph$Stimuli <- factor(associative_wide_matching_collapsed_valence_effects_graph$Stimuli,
                                       levels = c(1, 2, 3),
                                       labels = c("Happy", "Neutral", "Sad"))

ggplot(associative_wide_matching_collapsed_valence_effects_graph, aes(x = x, y = predicted, color = Stimuli)) + geom_line(size = 1) + geom_ribbon(data = associative_wide_matching_collapsed_valence_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Accuracy (%)") + ylab("BDI-II Total") + ggtitle("Linear Regression Valence Task Predicted Effects of Accuracy (%) on BDI-II Scores")
```
#### Reaction Times
```{r}
#Descriptives Graphs
ggplot(data = associative_long_matching_collapsed_valence, aes(x = mean_rt_mult, y = BDI_tot, color = stimuli)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("BDI-II Total Scores") + xlab("Reaction Times (ms)") + ggtitle("Valence Task Linear Relationship between BDI-II Scores and Accuracy")

#Regression Model
associative_wide_matching_collapsed_valence <- associative_wide_matching_collapsed %>%
  filter(Task == "Valence") #Restricting Task to "Valence"

assoc_valence_rt_BDI_model <- lm(BDI_tot ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low, data = associative_wide_matching_collapsed_valence)
jtools::summ(assoc_valence_rt_BDI_model, confint = TRUE, digits = 3)

#Standardised coefficients
assoc_valence_rt_BDI_model_standard <- standardize_parameters(assoc_valence_rt_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_valence_rt_BDI_model_standard

#Regression Graphs
assoc_valence_rt_BDI_model_effects <- ggpredict(assoc_valence_rt_BDI_model)

mean_rt_mult_High_effects <- as.data.frame(assoc_valence_rt_BDI_model_effects$mean_rt_mult_High)
mean_rt_mult_Medium_effects <- as.data.frame(assoc_valence_rt_BDI_model_effects$mean_rt_mult_Medium)
mean_rt_mult_Low_effects <- as.data.frame(assoc_valence_rt_BDI_model_effects$mean_rt_mult_Low)

associative_wide_matching_collapsed_valence_effects_graph <- bind_rows(mean_rt_mult_High_effects, mean_rt_mult_Medium_effects, mean_rt_mult_Low_effects)

associative_wide_matching_collapsed_valence_effects_graph <- associative_wide_matching_collapsed_valence_effects_graph %>%
  mutate(Stimuli = ifelse(group == "mean_rt_mult_High", 1, ifelse(group == "mean_rt_mult_Medium", 2, 3)))

associative_wide_matching_collapsed_valence_effects_graph$Stimuli <- factor(associative_wide_matching_collapsed_valence_effects_graph$Stimuli,
                                       levels = c(1, 2, 3),
                                       labels = c("Happy", "Neutral", "Sad"))

ggplot(associative_wide_matching_collapsed_valence_effects_graph, aes(x = x, y = predicted, color = Stimuli)) + geom_line(size = 1) + geom_ribbon(data = associative_wide_matching_collapsed_valence_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Reaction Times (ms)") + ylab("BDI-II Total") + ggtitle("Linear Regression Valence Task Predicted Effects of Reaction Times (ms) on BDI-II Scores")
```


# GNAT

```{r}
#Excluding participants according to a priori criteria
GNAT_long_Qs_eligible <- GNAT_long_Qs %>%
  filter(ineligible == "Include")

GNAT_wide_Qs_eligible <- GNAT_wide_Qs %>%
  filter(ineligible == "Include")

GNAT_trial_clean_Qs_eligible <- GNAT_trial_clean_Qs %>%
  filter(ineligible == "Include")

```


## Overall Descriptives
```{r}

#Discriminative Accuracy
GNAT_long_Qs_eligible %>%
  tab_cols(Condition) %>%
  tab_cells(d) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("GNAT Discriminative Accuracy by Condition")

GNAT_long_Qs_eligible %>%
  group_by(Condition) %>%
  summarize(Mean = mean(d),  SD = sd(d), Minimum = min(d), Maximum = max(d)) %>%
  ggplot(., aes(x = Condition, y = Mean, fill = Condition)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Discriminative Accuracy") + 
  xlab ("Condition") + 
  labs(fill = "Condition") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))+
  scale_fill_brewer()
```

## Main Effects
```{r}
GNAT_long_Qs_eligible <- GNAT_long_Qs_eligible %>%
  dplyr::mutate(Referential_Condition = ifelse(Condition == "Self Positive" | Condition == "Self Negative", 1, 2)) %>%
  dplyr::mutate(Valence = ifelse(Condition == "Self Positive" | Condition == "Other Positive", 1, 2)) %>%
  dplyr::select(subject, Condition, Referential_Condition, Valence, everything())


GNAT_long_Qs_eligible$Referential_Condition <- factor(GNAT_long_Qs_eligible$Referential_Condition,
                                             levels = c(1, 2),
                                             labels = c("Self", "Other"))


GNAT_long_Qs_eligible$Valence <- factor(GNAT_long_Qs_eligible$Valence,
                                             levels = c(1, 2),
                                             labels = c("Positive", "Negative"))


gnat_d_main_model <- lmer(d ~ Referential_Condition*Valence + (1|subject), data = GNAT_long_Qs_eligible)
jtools::summ(gnat_d_main_model, confint = TRUE, digits = 3)
gnat_d_main_model_effects <- ggpredict(gnat_d_main_model, terms = c("Referential_Condition", "Valence"))
plot(gnat_d_main_model_effects)

contrasts_gnat_d_main_model <- emmeans(gnat_d_main_model, pairwise ~ Valence | Referential_Condition, adjust = "tukey")
contrasts_gnat_d_main_model
confint(contrasts_gnat_d_main_model)
 
 
#Standardised coefficients
gnat_d_main_model_standard <- standardize_parameters(gnat_d_main_model, method = "refit", ci = 0.95, include_response = TRUE)
gnat_d_main_model_standard
 

```


## Is discriminative accuracy associated with depression severity?

### PHQ-9

#### Descriptive Data
```{r}
GNAT_long_Qs_eligible %>%
  tab_cols(Condition) %>%
  tab_rows(PHQ_dep_group_3) %>%
  tab_cells(d) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("GNAT Discriminative Accuracy by Condition and PHQ-9 Severity")

ggplot(data = GNAT_long_Qs_eligible, aes(x = d, y = PHQ_tot, color = Condition)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("PHQ-9 Total Scores") + xlab("Discriminative Accuracy") + ggtitle("Association between Discriminative Accuracy and PHQ-9 Severity by Condition")

```
#### Regression Models
```{r}
GNAT_discrim_PHQ_model <- lm(PHQ_tot ~ d_Self_Pos + d_Self_Neg + d_Other_Pos + d_Other_Neg, data = GNAT_wide_Qs_eligible)
jtools::summ(GNAT_discrim_PHQ_model, confint = TRUE, digits = 3)

#Standardised coefficients
GNAT_discrim_PHQ_model_standard <- standardize_parameters(GNAT_discrim_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
GNAT_discrim_PHQ_model_standard

#Regression Graphs
GNAT_discrim_PHQ_model_effects <- ggpredict(GNAT_discrim_PHQ_model)

GNAT_d_Self_Pos_PHQ_effects <- as.data.frame(GNAT_discrim_PHQ_model_effects$d_Self_Pos)
GNAT_d_Self_Neg_PHQ_effects <- as.data.frame(GNAT_discrim_PHQ_model_effects$d_Self_Neg)
GNAT_d_Other_Pos_PHQ_effects <- as.data.frame(GNAT_discrim_PHQ_model_effects$d_Other_Pos)
GNAT_d_Other_Neg_PHQ_effects <- as.data.frame(GNAT_discrim_PHQ_model_effects$d_Other_Neg)


GNAT_discrim_PHQ_model_effects_graph <- bind_rows(GNAT_d_Self_Pos_PHQ_effects, GNAT_d_Self_Neg_PHQ_effects, GNAT_d_Other_Pos_PHQ_effects, GNAT_d_Other_Neg_PHQ_effects)

GNAT_discrim_PHQ_model_effects_graph <- GNAT_discrim_PHQ_model_effects_graph %>%
  mutate(Condition = ifelse(group == "d_Self_Pos", 1, ifelse(group == "d_Self_Neg", 2, ifelse(group == "d_Other_Pos", 3, 4))))

GNAT_discrim_PHQ_model_effects_graph$Condition <- factor(GNAT_discrim_PHQ_model_effects_graph$Condition,
                                       levels = c(1, 2, 3, 4),
                                       labels = c("Self Positive", "Self Negative", "Other Positive", "Other Negative"))

ggplot(GNAT_discrim_PHQ_model_effects_graph, aes(x = x, y = predicted, color = Condition)) + geom_line(size = 1) + geom_ribbon(data = GNAT_discrim_PHQ_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Discriminative Accuracy") + ylab("PHQ-9 Total") + ggtitle("Linear Regression GNAT Predicted Effects of Discriminative Accuracy on PHQ-9 Scores")
```

### BDI-II

#### Descriptive Data
```{r}
GNAT_long_Qs_eligible %>%
  tab_cols(Condition) %>%
  tab_rows(BDI_dep_group_3) %>%
  tab_cells(d) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("GNAT Discriminative Accuracy by Condition and BDI-II Severity")

ggplot(data = GNAT_long_Qs_eligible, aes(x = d, y = BDI_tot, color = Condition)) + geom_point() + geom_smooth(method = lm, fullrange = TRUE) + ylab("BDI-II Total Scores") + xlab("Discriminative Accuracy") + ggtitle("Association between Discriminative Accuracy and BDI-II Severity by Condition")

```
#### Regression Models
```{r}
GNAT_discrim_BDI_model <- lm(BDI_tot ~ d_Self_Pos + d_Self_Neg + d_Other_Pos + d_Other_Neg, data = GNAT_wide_Qs_eligible)
jtools::summ(GNAT_discrim_BDI_model, confint = TRUE, digits = 3)

#Standardised coefficients
GNAT_discrim_BDI_model_standard <- standardize_parameters(GNAT_discrim_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
GNAT_discrim_BDI_model_standard

#Regression Graphs
GNAT_discrim_BDI_model_effects <- ggpredict(GNAT_discrim_BDI_model)

GNAT_d_Self_Pos_BDI_effects <- as.data.frame(GNAT_discrim_BDI_model_effects$d_Self_Pos)
GNAT_d_Self_Neg_BDI_effects <- as.data.frame(GNAT_discrim_BDI_model_effects$d_Self_Neg)
GNAT_d_Other_Pos_BDI_effects <- as.data.frame(GNAT_discrim_BDI_model_effects$d_Other_Pos)
GNAT_d_Other_Neg_BDI_effects <- as.data.frame(GNAT_discrim_BDI_model_effects$d_Other_Neg)


GNAT_discrim_BDI_model_effects_graph <- bind_rows(GNAT_d_Self_Pos_BDI_effects, GNAT_d_Self_Neg_BDI_effects, GNAT_d_Other_Pos_BDI_effects, GNAT_d_Other_Neg_BDI_effects)

GNAT_discrim_BDI_model_effects_graph <- GNAT_discrim_BDI_model_effects_graph %>%
  mutate(Condition = ifelse(group == "d_Self_Pos", 1, ifelse(group == "d_Self_Neg", 2, ifelse(group == "d_Other_Pos", 3, 4))))

GNAT_discrim_BDI_model_effects_graph$Condition <- factor(GNAT_discrim_BDI_model_effects_graph$Condition,
                                       levels = c(1, 2, 3, 4),
                                       labels = c("Self Positive", "Self Negative", "Other Positive", "Other Negative"))

ggplot(GNAT_discrim_BDI_model_effects_graph, aes(x = x, y = predicted, color = Condition)) + geom_line(size = 1) + geom_ribbon(data = GNAT_discrim_BDI_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Discriminative Accuracy") + ylab("BDI-II Total") + ggtitle("Linear Regression GNAT Predicted Effects of Discriminative Accuracy on BDI-II Scores")
```

### If we include all participants are results consistent?

#### PHQ
```{r}
GNAT_discrim_PHQ_model_sensitivity <- lm(PHQ_tot ~ d_Self_Pos + d_Self_Neg + d_Other_Pos + d_Other_Neg, data = GNAT_wide_Qs)
jtools::summ(GNAT_discrim_PHQ_model_sensitivity, confint = TRUE, digits = 3)

#Standardised coefficients
GNAT_discrim_PHQ_model_sensitivity_standard <- standardize_parameters(GNAT_discrim_PHQ_model_sensitivity, method = "refit", ci = 0.95, include_response = TRUE)
GNAT_discrim_PHQ_model_sensitivity_standard
```

#### BDI
```{r}
GNAT_discrim_BDI_model_sensitivity <- lm(BDI_tot ~ d_Self_Pos + d_Self_Neg + d_Other_Pos + d_Other_Neg, data = GNAT_wide_Qs)
jtools::summ(GNAT_discrim_BDI_model_sensitivity, confint = TRUE, digits = 3)

#Standardised coefficients
GNAT_discrim_BDI_model_sensitivity_standard <- standardize_parameters(GNAT_discrim_BDI_model_sensitivity, method = "refit", ci = 0.95, include_response = TRUE)
GNAT_discrim_BDI_model_sensitivity_standard
```


# Social Evaluation Learning

## Learning Curves

### PHQ-9

```{r}
SEL_trial <- SEL_trial %>%
  group_by(subject, session, condition, rule) %>%
  mutate(cum_acc = cumsum(acc)) %>%
  mutate(prop_acc = cum_acc/SubTrial) %>%
  ungroup()

SEL_trial$rule <- as.factor(SEL_trial$rule)

SEL_trial <- SEL_trial %>%
  mutate(condition_rule = ifelse(condition == 1 & rule == 60 | condition == 1 & rule == 80, 1, ifelse(condition == 1 & rule == 20 | condition == 1 & rule == 40, 2, ifelse(condition == 2 & rule == 60 | condition == 2 & rule == 80, 3, ifelse(condition == 2 & rule == 20 | condition == 2 & rule == 40, 4, ifelse(condition == 3 & rule == 60 | condition == 3 & rule == 80, 5, 6))))))

SEL_trial$condition_rule <- factor(SEL_trial$condition_rule,
                                       levels = c(1, 2, 3, 4, 5, 6),
                                       labels = c("Self-Like", "Self-Dislike", "Friend-Like", "Friend-Dislike", "Stranger-Like", "Stranger-Dislike"))

SEL_trial %>%
  group_by(PHQ_dep_group_3, condition_rule, SubTrial) %>%
  summarize(Mean = mean(prop_acc)) %>%
  ggplot(., aes(x = SubTrial, y = Mean, color = condition_rule)) + facet_wrap(~PHQ_dep_group_3)  + geom_line() + ylab("Accuracy (%)") + xlab("Trial Number") + labs(color = "Rule")
      
SEL_trial %>%
  filter(condition_rule == "Self-Like" | condition_rule == "Self-Dislike") %>%
  group_by(PHQ_dep_group_3, condition_rule, SubTrial) %>%
  summarize(Mean = mean(prop_acc)) %>%
  ggplot(., aes(x = SubTrial, y = Mean, color = condition_rule)) + facet_wrap(~PHQ_dep_group_3)  + geom_line() + ylab("Accuracy (%)") + xlab("Trial Number") + labs(color = "Rule")


SEL_trial$Condition <- factor(SEL_trial$condition,
                              levels = c(1, 2, 3),
                              labels = c("Self", "Friend", "Stranger"))

SEL_trial <- SEL_trial %>%
  mutate(Overall_rule = ifelse(rule == 60 | rule == 80, 1, 2))

SEL_trial$Overall_rule <- factor(SEL_trial$Overall_rule,
                              levels = c(1, 2),
                              labels = c("Positive", "Negative"))

SEL_trial %>%
  group_by(PHQ_dep_group_3, Condition, Overall_rule, SubTrial) %>%
  summarize(Mean = mean(prop_acc)) %>%
  ggplot(., aes(x = SubTrial, y = Mean, color = Overall_rule)) + facet_wrap(~ Condition + PHQ_dep_group_3)  + geom_line() + ylab("Accuracy (%)") + xlab("Trial Number") + labs(color = "Rule")

# New facet label names for BDI
PHQ_supp.labs <- c("None \n (PHQ-9 \u003C 5)", "Mild \n (PHQ-9 5-9)", "Moderate to Severe \n (PHQ-9 \u2265 10)")
names(PHQ_supp.labs) <- c("None", "Mild", "Moderate to Severe")


# Self learning rates by BDI
SEL_learning_rates_PHQ <- SEL_trial %>%
  filter(condition_rule == "Self-Like" | condition_rule == "Self-Dislike") %>%
  group_by(PHQ_dep_group_3, condition_rule, SubTrial) %>%
  summarize(Mean = mean(prop_acc)) %>%
  ggplot(., aes(x = SubTrial, y = Mean)) + facet_wrap(~PHQ_dep_group_3, labeller = labeller(PHQ_dep_group_3 = PHQ_supp.labs))  + geom_line(aes(linetype = condition_rule), size = 0.8) + ylab("Accuracy (%)") + xlab("Trial Number") + labs(linetype = "Rule") + theme(panel.spacing = unit(1.5, "lines"),    panel.background = element_blank(),
                                                                                                                                                                                            panel.grid.major = element_line(size = 0.1, linetype = 'solid',
                                                                                                                                                                                                                            colour = "grey60"), 
                                                                                                                                                                                            panel.grid.minor = element_line(size = 0.1, linetype = 'solid',
                                                                                                                                                                                                                            colour = "grey60"))

SEL_learning_rates_PHQ

ggsave("SEL_learning_rates_PHQ.TIFF")


```

### BDI-II
```{r}
SEL_trial <- SEL_trial %>%
  group_by(subject, session, condition, rule) %>%
  mutate(cum_acc = cumsum(acc)) %>%
  mutate(prop_acc = cum_acc/SubTrial) %>%
  ungroup()

SEL_trial$rule <- as.factor(SEL_trial$rule)

SEL_trial <- SEL_trial %>%
  mutate(condition_rule = ifelse(condition == 1 & rule == 60 | condition == 1 & rule == 80, 1, ifelse(condition == 1 & rule == 20 | condition == 1 & rule == 40, 2, ifelse(condition == 2 & rule == 60 | condition == 2 & rule == 80, 3, ifelse(condition == 2 & rule == 20 | condition == 2 & rule == 40, 4, ifelse(condition == 3 & rule == 60 | condition == 3 & rule == 80, 5, 6))))))

SEL_trial$condition_rule <- factor(SEL_trial$condition_rule,
                                       levels = c(1, 2, 3, 4, 5, 6),
                                       labels = c("Self-Like", "Self-Dislike", "Friend-Like", "Friend-Dislike", "Stranger-Like", "Stranger-Dislike"))

SEL_trial %>%
  group_by(BDI_dep_group_3, condition_rule, SubTrial) %>%
  summarize(Mean = mean(prop_acc)) %>%
  ggplot(., aes(x = SubTrial, y = Mean, color = condition_rule)) + facet_wrap(~BDI_dep_group_3)  + geom_line() + ylab("Accuracy (%)") + xlab("Trial Number") + labs(color = "Rule")
      
SEL_trial %>%
  filter(condition_rule == "Self-Like" | condition_rule == "Self-Dislike") %>%
  group_by(BDI_dep_group_3, condition_rule, SubTrial) %>%
  summarize(Mean = mean(prop_acc)) %>%
  ggplot(., aes(x = SubTrial, y = Mean, color = condition_rule)) + facet_wrap(~BDI_dep_group_3)  + geom_line() + ylab("Accuracy (%)") + xlab("Trial Number") + labs(color = "Rule")

# New facet label names for BDI
BDI_supp.labs <- c("None \n (BDI-II \u2264 13)", "Mild \n (BDI-II 14-19)", "Moderate to Severe \n (BDI-II \u2265 20)")
names(BDI_supp.labs) <- c("None", "Mild", "Moderate to Severe")

# Self learning rates by BDI
SEL_learning_rates_BDI <- SEL_trial %>%
  filter(condition_rule == "Self-Like" | condition_rule == "Self-Dislike") %>%
  group_by(BDI_dep_group_3, condition_rule, SubTrial) %>%
  summarize(Mean = mean(prop_acc)) %>%
  ggplot(., aes(x = SubTrial, y = Mean)) + facet_wrap(~BDI_dep_group_3, labeller = labeller(BDI_dep_group_3 = BDI_supp.labs))  + geom_line(aes(linetype = condition_rule), size = 0.8) + ylab("Accuracy (%)") + xlab("Trial Number") + labs(linetype = "Rule") + theme(panel.spacing = unit(1.5, "lines"),    panel.background = element_blank(),
                                                                                                                                                                                            panel.grid.major = element_line(size = 0.1, linetype = 'solid',
                                                                                                                                                                                                                            colour = "grey60"), 
                                                                                                                                                                                            panel.grid.minor = element_line(size = 0.1, linetype = 'solid',
                                                                                                                                                                                                                            colour = "grey60"))
SEL_learning_rates_BDI

ggsave("SEL_learning_rates_BDI.TIFF")

```


## Bias Scores

### Overall Descriptives
```{r}
SEL_bias %>%
  tab_cols(Condition) %>%
  tab_cells(bias) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Bias Scores by Referential Condition")

SEL_bias %>%
  group_by(Condition) %>%
  summarize(Mean = mean(bias),  SD = sd(bias), Minimum = min(bias), Maximum = max(bias)) %>%
  ggplot(., aes(x = Condition, y = Mean, fill = Condition)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Bias Scores (Positive - Negative)") + 
  xlab ("Condition") + 
  labs(fill = "Condition") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))+
  scale_fill_brewer()
```

## ICC
```{r}
overall_bias_ICC_df <- SEL_bias %>%
  dplyr::select(subject, session, bias) %>%
  group_by(subject, session) %>%
  summarise(bias = mean(bias)) %>%
  pivot_wider(values_from = bias, names_from = session) %>%
  ungroup() %>%
  dplyr::select(- subject)

ICC(overall_bias_ICC_df)
```

## Main Effects
```{r}
SEL_bias_main_model <- lmer(bias ~ Condition + (1|subject), data = SEL_bias)
jtools::summ(SEL_bias_main_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_bias_main_model_standard <- standardize_parameters(SEL_bias_main_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_bias_main_model_standard
```

### Do Bias Scores differ by PHQ-9 Severity?

```{r}

#Descriptives
SEL_bias %>%
  tab_cols(Condition) %>%
  tab_rows(PHQ_dep_group_3) %>%
  tab_cells(bias) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Bias Scores by Referential Condition and PHQ-9 Severity")


ggplot(data = SEL_bias, aes(x = bias, y = PHQ_tot, color = Condition)) + geom_smooth(method = lm, fullrange = TRUE) + ylab("PHQ-9 Total Scores") + xlab("Bias Scores (Positive - Negative)") + ggtitle("Association between Bias Scores and PHQ-9 Severity by Condition")

SEL_bias_PHQ_data <- ggplot(SEL_bias, aes(x = bias, y = PHQ_tot)) + geom_point(alpha=0.7) + geom_smooth(method = lm, fullrange = TRUE, colour = "black") + facet_wrap(~ Condition) + ylab("PHQ-9 Total Scores") + xlab("Bias Scores (Positive - Negative)") + theme_minimal() + scale_colour_grey()
ggsave("SEL_bias_PHQ_data.TIFF")
SEL_bias_PHQ_data

#Reshaping to wide for regression models
SEL_bias_wide <- SEL_bias %>%
  dplyr::select(subject, session, Condition, bias, PHQ_tot, BDI_tot, BFNE_tot, primary_depression_CIS) %>%
  spread(key = Condition, value = bias) 
  
#Regression Models
SEL_bias_PHQ_model <- lmer(PHQ_tot ~ Self + Friend + Stranger + session + (1|subject), data = SEL_bias_wide) 
jtools::summ(SEL_bias_PHQ_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_bias_PHQ_model_standard <- standardize_parameters(SEL_bias_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_bias_PHQ_model_standard

#Regression Graphs
SEL_bias_PHQ_model_effects <- ggpredict(SEL_bias_PHQ_model)

Self_effects <- as.data.frame(SEL_bias_PHQ_model_effects$Self)
Friend_effects <- as.data.frame(SEL_bias_PHQ_model_effects$Friend)
Stranger_effects <- as.data.frame(SEL_bias_PHQ_model_effects$Stranger)

SEL_bias_PHQ_model_effects_graph <- bind_rows(Self_effects, Friend_effects, Stranger_effects)

SEL_bias_PHQ_model_effects_graph <- SEL_bias_PHQ_model_effects_graph %>%
  mutate(Condition = ifelse(group == "Self", 1, ifelse(group == "Friend", 2, 3)))

SEL_bias_PHQ_model_effects_graph$Condition <- factor(SEL_bias_PHQ_model_effects_graph$Condition,
                                       levels = c(1, 2, 3),
                                       labels = c("Self", "Friend", "Stranger"))

ggplot(SEL_bias_PHQ_model_effects_graph, aes(x = x, y = predicted, color = Condition)) + geom_line(size = 1) + geom_ribbon(data = SEL_bias_PHQ_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Bias Scores (Positive - Negative)") + ylab("PHQ-9 Total") + ggtitle("Linear Regression SEL Bias Scores on PHQ-9 Scores")

PHQ_SEL_bias <- ggplot(SEL_bias_PHQ_model_effects_graph, aes(x = x, y = predicted, color = Condition)) + geom_line(size = 1) + geom_ribbon(data = SEL_bias_PHQ_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Bias Scores (Positive - Negative)") + ylab("PHQ-9 Total")

ggsave("PHQ_SEL_bias.jpeg")

```
### Do Bias Scores differ by BDI-II Severity?

```{r}

#Descriptives
SEL_bias %>%
  tab_cols(Condition) %>%
  tab_rows(BDI_dep_group_3) %>%
  tab_cells(bias) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Bias Scores by Referential Condition and BDI-II Severity")

ggplot(data = SEL_bias, aes(x = bias, y = BDI_tot, color = Condition)) + geom_smooth(method = lm, fullrange = TRUE) + ylab("BDI-II Total Scores") + xlab("Bias Scores (Positive - Negative)") + ggtitle("Association between Bias Scores and BDI-II Severity by Condition")

SEL_bias_BDI_data <- ggplot(SEL_bias, aes(x = bias, y = BDI_tot)) + geom_point(alpha=0.7) + geom_smooth(method = lm, fullrange = TRUE, colour = "black") + facet_wrap(~ Condition) + ylab("BDI-II Total Scores") + xlab("Bias Scores (Positive - Negative)")+ theme_minimal() + scale_colour_grey()
ggsave("SEL_bias_BDI_data.TIFF")
SEL_bias_BDI_data

#Reshaping to wide for regression models
SEL_bias_wide <- SEL_bias %>%
  dplyr::select(subject, session, Condition, bias, BDI_tot, BDI_tot, BFNE_tot, primary_depression_CIS) %>%
  spread(key = Condition, value = bias) 
  
#Regression Models
SEL_bias_BDI_model <- lmer(BDI_tot ~ Self + Friend + Stranger + session + (1|subject), data = SEL_bias_wide) 
jtools::summ(SEL_bias_BDI_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_bias_BDI_model_standard <- standardize_parameters(SEL_bias_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_bias_BDI_model_standard

#Regression Graphs
SEL_bias_BDI_model_effects <- ggpredict(SEL_bias_BDI_model)

Self_effects <- as.data.frame(SEL_bias_BDI_model_effects$Self)
Friend_effects <- as.data.frame(SEL_bias_BDI_model_effects$Friend)
Stranger_effects <- as.data.frame(SEL_bias_BDI_model_effects$Stranger)

SEL_bias_BDI_model_effects_graph <- bind_rows(Self_effects, Friend_effects, Stranger_effects)

SEL_bias_BDI_model_effects_graph <- SEL_bias_BDI_model_effects_graph %>%
  mutate(Condition = ifelse(group == "Self", 1, ifelse(group == "Friend", 2, 3)))

SEL_bias_BDI_model_effects_graph$Condition <- factor(SEL_bias_BDI_model_effects_graph$Condition,
                                       levels = c(1, 2, 3),
                                       labels = c("Self", "Friend", "Stranger"))

ggplot(SEL_bias_BDI_model_effects_graph, aes(x = x, y = predicted, color = Condition)) + geom_line(size = 1) + geom_ribbon(data = SEL_bias_BDI_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Bias Scores (Positive - Negative)") + ylab("BDI-II Total") + ggtitle("Linear Regression SEL Bias Scores on BDI-II Scores")

BDI_SEL_bias <- ggplot(SEL_bias_BDI_model_effects_graph, aes(x = x, y = predicted, color = Condition)) + geom_line(size = 1) + geom_ribbon(data = SEL_bias_BDI_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Bias Scores (Positive - Negative)") + ylab("BDI-II Total") 

ggsave("BDI_SEL_bias.jpeg")

```

### Does the relationship between depression and self bias differ between sessions?

#### PHQ
```{r}

SEL_bias_wide <- SEL_bias %>%
  dplyr::select(subject, session, Condition, bias, PHQ_tot, BDI_tot, BFNE_tot, primary_depression_CIS) %>%
  spread(key = Condition, value = bias) 
  
#Regression Models
SEL_bias_PHQ_session_model <- lmer(PHQ_tot ~ Self*session + (1|subject), data = SEL_bias_wide) 
jtools::summ(SEL_bias_PHQ_session_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_bias_PHQ_session_model_standard <- standardize_parameters(SEL_bias_PHQ_session_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_bias_PHQ_session_model_standard
```
#### BDI
```{r}
#Regression Models
SEL_bias_BDI_session_model <- lmer(BDI_tot ~ Self*session + (1|subject), data = SEL_bias_wide) 
jtools::summ(SEL_bias_BDI_session_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_bias_BDI_session_model_standard <- standardize_parameters(SEL_bias_BDI_session_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_bias_BDI_session_model_standard
```

## Errors to Criterion

### Overall Descriptives
```{r}
SEL_long %>%
  tab_cols(Condition %nest% Overall_Rule) %>%
  tab_cells(ERR_to_crit) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Errors to Criterion by Referential Condition and Rule")

SEL_long %>%
  group_by(Condition, Overall_Rule) %>%
  summarize(Mean = mean(ERR_to_crit),  SD = sd(ERR_to_crit), Minimum = min(ERR_to_crit), Maximum = max(ERR_to_crit)) %>%
  ggplot(., aes(x = Condition, y = Mean, fill = Overall_Rule)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Errors to Criterion") + 
  xlab ("Condition") + 
  labs(fill = "Rule") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))+
  scale_fill_brewer()

SEL_long %>%
  tab_cols(Condition %nest% Overall_Rule) %>%
  tab_cells(POS_resp) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Positive Response Rate by Referential Condition and Rule")

SEL_long %>%
  group_by(Condition, Overall_Rule) %>%
  summarize(Mean = mean(POS_resp),  SD = sd(POS_resp), Minimum = min(POS_resp), Maximum = max(POS_resp)) %>%
  ggplot(., aes(x = Condition, y = Mean, fill = Overall_Rule)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Positive Response Rate") + 
  xlab ("Condition") + 
  labs(fill = "Rule") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))+
  scale_fill_brewer()


```


## Main Effects
```{r}
SEL_err_main_model <- lmer(ERR_to_crit ~ Condition*Overall_Rule + (1|subject), data = SEL_long)
jtools::summ(SEL_err_main_model, confint = TRUE, digits = 3)

#Overall p value referential condition
SEL_err_main_model_no_referential_condition <- lmer(ERR_to_crit ~ Overall_Rule + (1|subject), data = SEL_long)

anova(SEL_err_main_model, SEL_err_main_model_no_referential_condition)

#Overall p value rule
SEL_err_main_model_no_rule <- lmer(ERR_to_crit ~ Condition + (1|subject), data = SEL_long)

anova(SEL_err_main_model, SEL_err_main_model_no_rule)

#Overall p value referential condition interaction effect with rule
SEL_err_main_model_no_interaction <- lmer(ERR_to_crit ~ Overall_Rule + Condition + (1|subject), data = SEL_long)

anova(SEL_err_main_model, SEL_err_main_model_no_interaction)

#Standardised coefficients
SEL_err_main_model_standard <- standardize_parameters(SEL_err_main_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_err_main_model_standard
```




### Do Errors to Criterion differ by PHQ-9 Severity?
```{r}

#Descriptives
SEL_long %>%
  tab_cols(Condition %nest% Overall_Rule) %>%
  tab_rows(PHQ_dep_group_3) %>%
  tab_cells(ERR_to_crit) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Errors to Criterion by Referential Condition, Rule and PHQ-9 Severity")

ggplot(data = SEL_long, aes(x = ERR_to_crit, y = PHQ_tot, color = Overall_Rule)) + facet_wrap(~ Condition) + geom_smooth(method = lm, fullrange = TRUE) + ylab("PHQ-9 Total Scores") + xlab("Errors to Criterion") + ggtitle("Association between Errors to Criterion and PHQ-9 Severity by Condition and Rule") + labs(color = "Rule")



#Regression Models
SEL_err_PHQ_model <- lmer(PHQ_tot ~ ERR_to_crit_Self_Positive + ERR_to_crit_Self_Negative + ERR_to_crit_Friend_Positive + ERR_to_crit_Friend_Negative + ERR_to_crit_Stranger_Positive + ERR_to_crit_Stranger_Negative + session + (1|subject), data = SEL_wide) 
jtools::summ(SEL_err_PHQ_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_err_PHQ_model_standard <- standardize_parameters(SEL_err_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_err_PHQ_model_standard

#Regression Graphs
SEL_err_PHQ_model_effects <- ggpredict(SEL_err_PHQ_model)

Self_Positive_effects <- as.data.frame(SEL_err_PHQ_model_effects$ERR_to_crit_Self_Positive)
Self_Negative_effects <- as.data.frame(SEL_err_PHQ_model_effects$ERR_to_crit_Self_Negative)
Friend_Positive_effects <- as.data.frame(SEL_err_PHQ_model_effects$ERR_to_crit_Friend_Positive)
Friend_Negative_effects <- as.data.frame(SEL_err_PHQ_model_effects$ERR_to_crit_Friend_Negative)
Stranger_Positive_effects <- as.data.frame(SEL_err_PHQ_model_effects$ERR_to_crit_Stranger_Positive)
Stranger_Negative_effects <- as.data.frame(SEL_err_PHQ_model_effects$ERR_to_crit_Stranger_Negative)


SEL_err_PHQ_model_effects_graph <- bind_rows(Self_Positive_effects, Self_Negative_effects, Friend_Positive_effects, Friend_Negative_effects, Stranger_Positive_effects, Stranger_Negative_effects)

SEL_err_PHQ_model_effects_graph <- SEL_err_PHQ_model_effects_graph %>%
  mutate(ConditionRule = ifelse(group == "ERR_to_crit_Self_Positive", 1, ifelse(group == "ERR_to_crit_Self_Negative", 2, ifelse(group == "ERR_to_crit_Friend_Positive", 3, ifelse(group == "ERR_to_crit_Friend_Negative", 4, ifelse(group == "ERR_to_crit_Stranger_Positive", 5, 6)))))) %>%
  mutate(Condition = ifelse(ConditionRule <=2, 1, ifelse(ConditionRule > 2 & ConditionRule <= 4, 2, 3))) %>%
  mutate(Rule = ifelse(ConditionRule == 1 | ConditionRule == 3 | ConditionRule == 5, 1, 2))
                                 
SEL_err_PHQ_model_effects_graph$Condition <- factor(SEL_err_PHQ_model_effects_graph$Condition,
                                       levels = c(1, 2, 3),
                                       labels = c("Self", "Friend", "Stranger"))

SEL_err_PHQ_model_effects_graph$Rule <- factor(SEL_err_PHQ_model_effects_graph$Rule,
                                       levels = c(1, 2),
                                       labels = c("Positive", "Negative"))


ggplot(SEL_err_PHQ_model_effects_graph, aes(x = x, y = predicted, color = Rule)) + facet_wrap(~ Condition) + geom_line(size = 1) + geom_ribbon(data = SEL_err_PHQ_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Errors to Criterion") + ylab("PHQ-9 Total") + ggtitle("Linear Regression SEL Errors to Criterion on PHQ-9 Scores")

PHQ_SEL_err <- ggplot(SEL_err_PHQ_model_effects_graph, aes(x = x, y = predicted, color = Rule)) + facet_wrap(~ Condition) + geom_line(size = 1) + geom_ribbon(data = SEL_err_PHQ_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Errors to Criterion") + ylab("PHQ-9 Total")

ggsave("PHQ_SEL_err.jpeg")

```

### Do Errors to Criterion differ by BDI-II Severity?
```{r}

#Descriptives
SEL_long %>%
  tab_cols(Condition %nest% Overall_Rule) %>%
  tab_rows(BDI_dep_group_3) %>%
  tab_cells(ERR_to_crit) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Errors to Criterion by Referential Condition, Rule and BDI-II Severity")

ggplot(data = SEL_long, aes(x = ERR_to_crit, y = BDI_tot, color = Overall_Rule)) + facet_wrap(~ Condition) + geom_smooth(method = lm, fullrange = TRUE) + ylab("BDI-II Total Scores") + xlab("Errors to Criterion") + ggtitle("Association between Errors to Criterion and BDI-II Severity by Condition and Rule") + labs(color = "Rule")
  
#Regression Models
SEL_err_BDI_model <- lmer(BDI_tot ~ ERR_to_crit_Self_Positive + ERR_to_crit_Self_Negative + ERR_to_crit_Friend_Positive + ERR_to_crit_Friend_Negative + ERR_to_crit_Stranger_Positive + ERR_to_crit_Stranger_Negative + (1|subject), data = SEL_wide) 
jtools::summ(SEL_err_BDI_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_err_BDI_model_standard <- standardize_parameters(SEL_err_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_err_BDI_model_standard

#Regression Graphs
SEL_err_BDI_model_effects <- ggpredict(SEL_err_BDI_model)

Self_Positive_effects <- as.data.frame(SEL_err_BDI_model_effects$ERR_to_crit_Self_Positive)
Self_Negative_effects <- as.data.frame(SEL_err_BDI_model_effects$ERR_to_crit_Self_Negative)
Friend_Positive_effects <- as.data.frame(SEL_err_BDI_model_effects$ERR_to_crit_Friend_Positive)
Friend_Negative_effects <- as.data.frame(SEL_err_BDI_model_effects$ERR_to_crit_Friend_Negative)
Stranger_Positive_effects <- as.data.frame(SEL_err_BDI_model_effects$ERR_to_crit_Stranger_Positive)
Stranger_Negative_effects <- as.data.frame(SEL_err_BDI_model_effects$ERR_to_crit_Stranger_Negative)


SEL_err_BDI_model_effects_graph <- bind_rows(Self_Positive_effects, Self_Negative_effects, Friend_Positive_effects, Friend_Negative_effects, Stranger_Positive_effects, Stranger_Negative_effects)

SEL_err_BDI_model_effects_graph <- SEL_err_BDI_model_effects_graph %>%
  mutate(ConditionRule = ifelse(group == "ERR_to_crit_Self_Positive", 1, ifelse(group == "ERR_to_crit_Self_Negative", 2, ifelse(group == "ERR_to_crit_Friend_Positive", 3, ifelse(group == "ERR_to_crit_Friend_Negative", 4, ifelse(group == "ERR_to_crit_Stranger_Positive", 5, 6)))))) %>%
  mutate(Condition = ifelse(ConditionRule <=2, 1, ifelse(ConditionRule > 2 & ConditionRule <= 4, 2, 3))) %>%
  mutate(Rule = ifelse(ConditionRule == 1 | ConditionRule == 3 | ConditionRule == 5, 1, 2))
                                 
SEL_err_BDI_model_effects_graph$Condition <- factor(SEL_err_BDI_model_effects_graph$Condition,
                                       levels = c(1, 2, 3),
                                       labels = c("Self", "Friend", "Stranger"))

SEL_err_BDI_model_effects_graph$Rule <- factor(SEL_err_BDI_model_effects_graph$Rule,
                                       levels = c(1, 2),
                                       labels = c("Positive", "Negative"))


ggplot(SEL_err_BDI_model_effects_graph, aes(x = x, y = predicted, color = Rule)) + facet_wrap(~ Condition) + geom_line(size = 1) + geom_ribbon(data = SEL_err_BDI_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Errors to Criterion") + ylab("BDI-II Total") + ggtitle("Linear Regression SEL Errors to Criterion on BDI-II Scores")

BDI_SEL_err <- ggplot(SEL_err_BDI_model_effects_graph, aes(x = x, y = predicted, color = Rule)) + facet_wrap(~ Condition) + geom_line(size = 1) + geom_ribbon(data = SEL_err_BDI_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Errors to Criterion") + ylab("BDI-II Total")

ggsave("BDI_SEL_err.jpeg")
```





## Global Ratings

### Overall Descriptives
```{r}
SEL_long %>%
  tab_cols(Condition %nest% Overall_Rule) %>%
  tab_cells(globalrating) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Global Ratings by Referential Condition and Rule")

SEL_long %>%
  group_by(Condition, Overall_Rule) %>%
  summarize(Mean = mean(globalrating),  SD = sd(globalrating), Minimum = min(globalrating), Maximum = max(globalrating)) %>%
  ggplot(., aes(x = Condition, y = Mean, fill = Overall_Rule)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Global Ratings") + 
  xlab ("Condition") + 
  labs(fill = "Rule") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))+
  scale_fill_brewer()
```

### Main Effects

```{r}
SEL_global_main_model <- lmer(globalrating ~ Condition*Overall_Rule + (1|subject), data = SEL_long)
jtools::summ(SEL_global_main_model, confint = TRUE, digits = 3)

#Overall p value referential condition
SEL_global_main_model_no_referential_condition <- lmer(globalrating ~ Overall_Rule + (1|subject), data = SEL_long)

anova(SEL_global_main_model, SEL_global_main_model_no_referential_condition)

#Overall p value rule
SEL_global_main_model_no_rule <- lmer(globalrating ~ Condition + (1|subject), data = SEL_long)

anova(SEL_global_main_model, SEL_global_main_model_no_rule)

#Overall p value referential condition interaction effect with rule
SEL_global_main_model_no_interaction <- lmer(globalrating ~ Overall_Rule + Condition + (1|subject), data = SEL_long)

anova(SEL_global_main_model, SEL_global_main_model_no_interaction)

#Standardised coefficients
SEL_global_main_model_standard <- standardize_parameters(SEL_global_main_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_global_main_model_standard
```


### Do Global Ratings differ by PHQ-9 Severity?
```{r}

#Descriptives
SEL_long %>%
  tab_cols(Condition %nest% Overall_Rule) %>%
  tab_rows(PHQ_dep_group_3) %>%
  tab_cells(globalrating) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Global Ratings by Referential Condition, Rule and PHQ-9 Severity")

ggplot(data = SEL_long, aes(x = globalrating, y = PHQ_tot, color = Overall_Rule)) + facet_wrap(~ Condition) + geom_smooth(method = lm, fullrange = TRUE) + ylab("PHQ-9 Total Scores") + xlab("Global Ratings") + ggtitle("Association between Global Ratings and PHQ-9 Severity by Condition and Rule") + labs(color = "Rule")
  
#Regression Models
SEL_global_PHQ_model <- lmer(PHQ_tot ~ globalrating_Self_Positive + globalrating_Self_Negative + globalrating_Friend_Positive + globalrating_Friend_Negative + globalrating_Stranger_Positive + globalrating_Stranger_Negative + (1|subject), data = SEL_wide) 
jtools::summ(SEL_global_PHQ_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_global_PHQ_model_standard <- standardize_parameters(SEL_global_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_global_PHQ_model_standard

#Regression Graphs
SEL_global_PHQ_model_effects_global <- ggpredict(SEL_global_PHQ_model)

Self_Positive_effects_global <- as.data.frame(SEL_global_PHQ_model_effects_global$globalrating_Self_Positive)
Self_Negative_effects_global <- as.data.frame(SEL_global_PHQ_model_effects_global$globalrating_Self_Negative)
Friend_Positive_effects_global <- as.data.frame(SEL_global_PHQ_model_effects_global$globalrating_Friend_Positive)
Friend_Negative_effects_global <- as.data.frame(SEL_global_PHQ_model_effects_global$globalrating_Friend_Negative)
Stranger_Positive_effects_global <- as.data.frame(SEL_global_PHQ_model_effects_global$globalrating_Stranger_Positive)
Stranger_Negative_effects_global <- as.data.frame(SEL_global_PHQ_model_effects_global$globalrating_Stranger_Negative)


SEL_global_PHQ_model_effects_global_graph <- bind_rows(Self_Positive_effects_global, Self_Negative_effects_global, Friend_Positive_effects_global, Friend_Negative_effects_global, Stranger_Positive_effects_global, Stranger_Negative_effects_global)

SEL_global_PHQ_model_effects_global_graph <- SEL_global_PHQ_model_effects_global_graph %>%
  mutate(ConditionRule = ifelse(group == "globalrating_Self_Positive", 1, ifelse(group == "globalrating_Self_Negative", 2, ifelse(group == "globalrating_Friend_Positive", 3, ifelse(group == "globalrating_Friend_Negative", 4, ifelse(group == "globalrating_Stranger_Positive", 5, 6)))))) %>%
  mutate(Condition = ifelse(ConditionRule <=2, 1, ifelse(ConditionRule > 2 & ConditionRule <= 4, 2, 3))) %>%
  mutate(Rule = ifelse(ConditionRule == 1 | ConditionRule == 3 | ConditionRule == 5, 1, 2))
                                 
SEL_global_PHQ_model_effects_global_graph$Condition <- factor(SEL_global_PHQ_model_effects_global_graph$Condition,
                                       levels = c(1, 2, 3),
                                       labels = c("Self", "Friend", "Stranger"))

SEL_global_PHQ_model_effects_global_graph$Rule <- factor(SEL_global_PHQ_model_effects_global_graph$Rule,
                                       levels = c(1, 2),
                                       labels = c("Positive", "Negative"))


ggplot(SEL_global_PHQ_model_effects_global_graph, aes(x = x, y = predicted, color = Rule)) + facet_wrap(~ Condition) + geom_line(size = 1) + geom_ribbon(data = SEL_global_PHQ_model_effects_global_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Global Ratings") + ylab("PHQ-9 Total") + ggtitle("Linear Regression SEL Global Ratings on PHQ-9 Scores")
```

### Do Global Ratings differ by BDI-II Severity?
```{r}

#Descriptives
SEL_long %>%
  tab_cols(Condition %nest% Overall_Rule) %>%
  tab_rows(BDI_dep_group_3) %>%
  tab_cells(globalrating) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Global Ratings by Referential Condition, Rule and BDI-II Severity")

ggplot(data = SEL_long, aes(x = globalrating, y = BDI_tot, color = Overall_Rule)) + facet_wrap(~ Condition) + geom_smooth(method = lm, fullrange = TRUE) + ylab("BDI-II Total Scores") + xlab("Global Ratings") + ggtitle("Association between Global Ratings and BDI-II Severity by Condition and Rule") + labs(color = "Rule")
  
#Regression Models
SEL_global_BDI_model <- lmer(BDI_tot ~ globalrating_Self_Positive + globalrating_Self_Negative + globalrating_Friend_Positive + globalrating_Friend_Negative + globalrating_Stranger_Positive + globalrating_Stranger_Negative + (1|subject), data = SEL_wide) 
jtools::summ(SEL_global_BDI_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_global_BDI_model_standard <- standardize_parameters(SEL_global_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_global_BDI_model_standard

#Regression Graphs
SEL_global_BDI_model_effects_global <- ggpredict(SEL_global_BDI_model)

Self_Positive_effects_global <- as.data.frame(SEL_global_BDI_model_effects_global$globalrating_Self_Positive)
Self_Negative_effects_global <- as.data.frame(SEL_global_BDI_model_effects_global$globalrating_Self_Negative)
Friend_Positive_effects_global <- as.data.frame(SEL_global_BDI_model_effects_global$globalrating_Friend_Positive)
Friend_Negative_effects_global <- as.data.frame(SEL_global_BDI_model_effects_global$globalrating_Friend_Negative)
Stranger_Positive_effects_global <- as.data.frame(SEL_global_BDI_model_effects_global$globalrating_Stranger_Positive)
Stranger_Negative_effects_global <- as.data.frame(SEL_global_BDI_model_effects_global$globalrating_Stranger_Negative)


SEL_global_BDI_model_effects_global_graph <- bind_rows(Self_Positive_effects_global, Self_Negative_effects_global, Friend_Positive_effects_global, Friend_Negative_effects_global, Stranger_Positive_effects_global, Stranger_Negative_effects_global)

SEL_global_BDI_model_effects_global_graph <- SEL_global_BDI_model_effects_global_graph %>%
  mutate(ConditionRule = ifelse(group == "globalrating_Self_Positive", 1, ifelse(group == "globalrating_Self_Negative", 2, ifelse(group == "globalrating_Friend_Positive", 3, ifelse(group == "globalrating_Friend_Negative", 4, ifelse(group == "globalrating_Stranger_Positive", 5, 6)))))) %>%
  mutate(Condition = ifelse(ConditionRule <=2, 1, ifelse(ConditionRule > 2 & ConditionRule <= 4, 2, 3))) %>%
  mutate(Rule = ifelse(ConditionRule == 1 | ConditionRule == 3 | ConditionRule == 5, 1, 2))
                                 
SEL_global_BDI_model_effects_global_graph$Condition <- factor(SEL_global_BDI_model_effects_global_graph$Condition,
                                       levels = c(1, 2, 3),
                                       labels = c("Self", "Friend", "Stranger"))

SEL_global_BDI_model_effects_global_graph$Rule <- factor(SEL_global_BDI_model_effects_global_graph$Rule,
                                       levels = c(1, 2),
                                       labels = c("Positive", "Negative"))


ggplot(SEL_global_BDI_model_effects_global_graph, aes(x = x, y = predicted, color = Rule)) + facet_wrap(~ Condition) + geom_line(size = 1) + geom_ribbon(data = SEL_global_BDI_model_effects_global_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Global Ratings") + ylab("BDI-II Total") + ggtitle("Linear Regression SEL Global Ratings on BDI-II Scores")
```

## Do effects persist when social anxiety is taken into account?

### Bias Scores

#### PHQ
```{r}
SEL_bias_wide <- SEL_bias %>%
  dplyr::select(subject, session, Condition, bias, BDI_tot, PHQ_tot, BFNE_tot, primary_depression_CIS) %>%
  spread(key = Condition, value = bias) 

SEL_bias_PHQ_BFNE_model <- lmer(PHQ_tot ~ Self + Friend + Stranger + BFNE_tot + (1|subject), data = SEL_bias_wide) 
jtools::summ(SEL_bias_PHQ_BFNE_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_bias_PHQ_BFNE_model_standard <- standardize_parameters(SEL_bias_PHQ_BFNE_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_bias_PHQ_BFNE_model_standard
```


#### BDI
```{r}
SEL_bias_wide <- SEL_bias %>%
  dplyr::select(subject, session, Condition, bias, BDI_tot, PHQ_tot, BFNE_tot, primary_depression_CIS) %>%
  spread(key = Condition, value = bias) 

SEL_bias_BDI_BFNE_model <- lmer(BDI_tot ~ Self + Friend + Stranger + BFNE_tot + (1|subject), data = SEL_bias_wide) 
jtools::summ(SEL_bias_BDI_BFNE_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_bias_BDI_BFNE_model_standard <- standardize_parameters(SEL_bias_BDI_BFNE_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_bias_BDI_BFNE_model_standard
```


### Errors to Criterion

#### PHQ
```{r}
SEL_err_PHQ_BFNE_model <- lmer(PHQ_tot ~ ERR_to_crit_Self_Positive + ERR_to_crit_Self_Negative + ERR_to_crit_Friend_Positive + ERR_to_crit_Friend_Negative + ERR_to_crit_Stranger_Positive + ERR_to_crit_Stranger_Negative + BFNE_tot + (1|subject), data = SEL_wide) 
jtools::summ(SEL_err_PHQ_BFNE_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_err_PHQ_BFNE_model_standard <- standardize_parameters(SEL_err_PHQ_BFNE_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_err_PHQ_BFNE_model_standard
```


#### BDI
```{r}
SEL_err_BDI_BFNE_model <- lmer(BDI_tot ~ ERR_to_crit_Self_Positive + ERR_to_crit_Self_Negative + ERR_to_crit_Friend_Positive + ERR_to_crit_Friend_Negative + ERR_to_crit_Stranger_Positive + ERR_to_crit_Stranger_Negative + BFNE_tot + (1|subject), data = SEL_wide) 
jtools::summ(SEL_err_BDI_BFNE_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_err_BDI_BFNE_model_standard <- standardize_parameters(SEL_err_BDI_BFNE_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_err_BDI_BFNE_model_standard
```


### Global Ratings

#### PHQ
```{r}
SEL_global_PHQ_BFNE_model <- lmer(PHQ_tot ~ globalrating_Self_Positive + globalrating_Self_Negative + globalrating_Friend_Positive + globalrating_Friend_Negative + globalrating_Stranger_Positive + globalrating_Stranger_Negative + BFNE_tot + (1|subject), data = SEL_wide) 
jtools::summ(SEL_global_PHQ_BFNE_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_global_PHQ_BFNE_model_standard <- standardize_parameters(SEL_global_PHQ_BFNE_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_global_PHQ_BFNE_model_standard
```

#### BDI
```{r}
SEL_global_BDI_BFNE_model <- lmer(BDI_tot ~ globalrating_Self_Positive + globalrating_Self_Negative + globalrating_Friend_Positive + globalrating_Friend_Negative + globalrating_Stranger_Positive + globalrating_Stranger_Negative + BFNE_tot + (1|subject), data = SEL_wide) 
jtools::summ(SEL_global_BDI_BFNE_model, confint = TRUE, digits = 3)

#Standardised coefficients
SEL_global_BDI_BFNE_model_standard <- standardize_parameters(SEL_global_BDI_BFNE_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_global_BDI_BFNE_model_standard
```




# CIS-R Primary Diagnosis MDE

## Associative Learning

### Overall Descriptives
```{r}
associative_long_matching_collapsed %>%
  tab_cols(stimuli) %>%
  tab_rows(primary_depression_CIS) %>%
  tab_cells(prop_acc) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot() %>%
  set_caption("Associative Learning Self Accuracy (%) by Stimuli and MDE Diagnosis")



associative_long_matching_collapsed %>%
  tab_cols(stimuli) %>%
  tab_rows(primary_depression_CIS) %>%
  tab_cells(mean_rt_mult) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot() %>%
  set_caption("Associative Learning Self Reaction Times (ms) by Stimuli and MDE Diagnosis")


```

### Self

#### Accuracy
```{r}
#Regression Model
associative_wide_matching_collapsed_self <- associative_wide_matching_collapsed %>%
  filter(Task == "Self") %>%  #Restricting Task to "Self"
mutate(primary_depression_CIS = as.factor(primary_depression_CIS))

assoc_self_acc_CIS_model <- glm(primary_depression_CIS ~ prop_acc_High + prop_acc_Medium + prop_acc_Low, data = associative_wide_matching_collapsed_self, family = "binomial")
jtools::summ(assoc_self_acc_CIS_model, confint = TRUE, digits = 3)
logistic.display(assoc_self_acc_CIS_model, crude = FALSE)
assoc_self_acc_CIS_model_effects <- ggpredict(assoc_self_acc_CIS_model)
plot(assoc_self_acc_CIS_model_effects)
```
#### Reaction Times
```{r}
#Regression Model
associative_wide_matching_collapsed_self <- associative_wide_matching_collapsed %>%
  filter(Task == "Self") %>%  #Restricting Task to "Self"
mutate(primary_depression_CIS = as.factor(primary_depression_CIS))

assoc_self_rt_CIS_model <- glm(primary_depression_CIS ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low, data = associative_wide_matching_collapsed_self, family = "binomial")
jtools::summ(assoc_self_rt_CIS_model, confint = TRUE, digits = 3)
logistic.display(assoc_self_rt_CIS_model, crude = FALSE)
assoc_self_rt_CIS_model_effects <- ggpredict(assoc_self_rt_CIS_model)
plot(assoc_self_rt_CIS_model_effects)
```

### Reward

#### Accuracy
```{r}
#Regression Model
associative_wide_matching_collapsed_Reward <- associative_wide_matching_collapsed %>%
  filter(Task == "Reward") %>%  #Restricting Task to "Reward"
mutate(primary_depression_CIS = as.factor(primary_depression_CIS))

assoc_Reward_acc_CIS_model <- glm(primary_depression_CIS ~ prop_acc_High + prop_acc_Medium + prop_acc_Low, data = associative_wide_matching_collapsed_Reward, family = "binomial")
jtools::summ(assoc_Reward_acc_CIS_model, confint = TRUE, digits = 3)
logistic.display(assoc_Reward_acc_CIS_model, crude = FALSE)
assoc_Reward_acc_CIS_model_effects <- ggpredict(assoc_Reward_acc_CIS_model)
plot(assoc_Reward_acc_CIS_model_effects)
```
#### Reaction Times
```{r}
#Regression Model
associative_wide_matching_collapsed_Reward <- associative_wide_matching_collapsed %>%
  filter(Task == "Reward") %>%  #Restricting Task to "Reward"
mutate(primary_depression_CIS = as.factor(primary_depression_CIS))

assoc_Reward_rt_CIS_model <- glm(primary_depression_CIS ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low, data = associative_wide_matching_collapsed_Reward, family = "binomial")
jtools::summ(assoc_Reward_rt_CIS_model, confint = TRUE, digits = 3)
test <- logistic.display(assoc_Reward_rt_CIS_model, crude = FALSE)
assoc_Reward_rt_CIS_model_effects <- ggpredict(assoc_Reward_rt_CIS_model)
plot(assoc_Reward_rt_CIS_model_effects)



```
### Valence

#### Accuracy
```{r}
#Regression Model
associative_wide_matching_collapsed_Valence <- associative_wide_matching_collapsed %>%
  filter(Task == "Valence") %>%  #Restricting Task to "Valence"
mutate(primary_depression_CIS = as.factor(primary_depression_CIS))

assoc_Valence_acc_CIS_model <- glm(primary_depression_CIS ~ prop_acc_High + prop_acc_Medium + prop_acc_Low, data = associative_wide_matching_collapsed_Valence, family = "binomial")
jtools::summ(assoc_Valence_acc_CIS_model, confint = TRUE, digits = 3)
logistic.display(assoc_Valence_acc_CIS_model, crude = FALSE)
assoc_Valence_acc_CIS_model_effects <- ggpredict(assoc_Valence_acc_CIS_model)
plot(assoc_Valence_acc_CIS_model_effects)
```
#### Reaction Times
```{r}
#Regression Model
associative_wide_matching_collapsed_Valence <- associative_wide_matching_collapsed %>%
  filter(Task == "Valence") %>%  #Restricting Task to "Valence"
mutate(primary_depression_CIS = as.factor(primary_depression_CIS))

assoc_Valence_rt_CIS_model <- glm(primary_depression_CIS ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low, data = associative_wide_matching_collapsed_Valence, family = "binomial")
jtools::summ(assoc_Valence_rt_CIS_model, confint = TRUE, digits = 3)
logistic.display(assoc_Valence_rt_CIS_model, crude = FALSE)
assoc_Valence_rt_CIS_model_effects <- ggpredict(assoc_Valence_rt_CIS_model)
plot(assoc_Valence_rt_CIS_model_effects)
```

## Go/No-Go Task

## Descriptives

```{r}
GNAT_long_Qs_eligible$primary_depression_CIS <- as.factor(GNAT_long_Qs_eligible$primary_depression_CIS )

GNAT_long_Qs_eligible %>%
  tab_cols(Condition) %>%
  tab_rows(primary_depression_CIS) %>%
  tab_cells(d) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("GNAT Discriminative Accuracy by Condition and Primary MDE Diagnosis (1 = Yes)")

GNAT_long_Qs_eligible %>%
  group_by(primary_depression_CIS, Condition) %>%
  summarize(Mean = mean(d),  SD = sd(d), Minimum = min(d), Maximum = max(d)) %>%
  ggplot(., aes(x = Condition, y = Mean, fill = primary_depression_CIS)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Discriminative Accuracy") + 
  xlab ("Condition") + 
  scale_fill_discrete(name = "MDE Diagnosis", labels = c("No", "Yes")) +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))
```

## Regression Models
```{r}
GNAT_wide_Qs_eligible$primary_depression_CIS <- as.factor(GNAT_wide_Qs_eligible$primary_depression_CIS)

GNAT_discrim_CIS_model <- glm(primary_depression_CIS ~ d_Self_Pos + d_Self_Neg + d_Other_Pos + d_Other_Neg, data = GNAT_wide_Qs_eligible, family = "binomial")
jtools::summ(GNAT_discrim_CIS_model, confint = TRUE, digits = 3)
GNAT_discrim_CIS_model_odds <- logistic.display(GNAT_discrim_CIS_model, crude = FALSE)
GNAT_discrim_CIS_model_odds
GNAT_discrim_CIS_model_effects <- ggpredict(GNAT_discrim_CIS_model)
plot(GNAT_discrim_CIS_model_effects)

#Plot model with odds ratios
##Creating a variable to show conditions
GNAT_conditions <- c(1, 2, 3, 4)
GNAT_conditions <- factor(GNAT_conditions,
                          levels = c(1, 2, 3, 4),
                          labels = c("Self Positive", "Self Negative", "Other Positive", "Other Negative"))
##Creating a variable for odds
GNAT_d_odds <- c(as.numeric(str_extract(GNAT_discrim_CIS_model_odds$table[1], "[[:digit:]]\\.[[:digit:]][[:digit:]]")), as.numeric(str_extract(GNAT_discrim_CIS_model_odds$table[3], "[[:digit:]]\\.[[:digit:]][[:digit:]]")), as.numeric(str_extract(GNAT_discrim_CIS_model_odds$table[5], "[[:digit:]]\\.[[:digit:]][[:digit:]]")), as.numeric(str_extract(GNAT_discrim_CIS_model_odds$table[7], "[[:digit:]]\\.[[:digit:]][[:digit:]]")))
##Creating a variable for lower CI
GNAT_d_odds_lci_1 <- as.numeric(gsub("\\(", "", str_extract(GNAT_discrim_CIS_model_odds$table[1], "\\([[:digit:]]\\.[[:digit:]][[:digit:]]")))
GNAT_d_odds_lci_2 <- as.numeric(gsub("\\(", "", str_extract(GNAT_discrim_CIS_model_odds$table[3], "\\([[:digit:]]\\.[[:digit:]][[:digit:]]")))
GNAT_d_odds_lci_3 <- as.numeric(gsub("\\(", "", str_extract(GNAT_discrim_CIS_model_odds$table[5], "\\([[:digit:]]\\.[[:digit:]][[:digit:]]")))
GNAT_d_odds_lci_4 <- as.numeric(gsub("\\(", "", str_extract(GNAT_discrim_CIS_model_odds$table[7], "\\([[:digit:]]\\.[[:digit:]][[:digit:]]")))
GNAT_d_odds_lci <- c(GNAT_d_odds_lci_1, GNAT_d_odds_lci_2, GNAT_d_odds_lci_3, GNAT_d_odds_lci_4)
##Creating a variable for upper CI
GNAT_d_odds_uci_1 <- as.numeric(gsub("\\,", "", str_extract(GNAT_discrim_CIS_model_odds$table[1], "\\,[[:digit:]]\\.[[:digit:]][[:digit:]]")))
GNAT_d_odds_uci_2 <- as.numeric(gsub("\\,", "", str_extract(GNAT_discrim_CIS_model_odds$table[3], "\\,[[:digit:]]\\.[[:digit:]][[:digit:]]")))
GNAT_d_odds_uci_3 <- as.numeric(gsub("\\,", "", str_extract(GNAT_discrim_CIS_model_odds$table[5], "\\,[[:digit:]][[:digit:]]\\.[[:digit:]][[:digit:]]")))
GNAT_d_odds_uci_4 <- as.numeric(gsub("\\,", "", str_extract(GNAT_discrim_CIS_model_odds$table[7], "\\,[[:digit:]]\\.[[:digit:]][[:digit:]]")))
GNAT_d_odds_uci <- c(GNAT_d_odds_uci_1, GNAT_d_odds_uci_2, GNAT_d_odds_uci_3, GNAT_d_odds_uci_4)
##Binding
GNAT_d_odds_df <- data.frame(GNAT_conditions, GNAT_d_odds, GNAT_d_odds_lci, GNAT_d_odds_uci)
##Graph
GNAT_d_odds_df %>%
  ggplot(., aes(x = GNAT_conditions, y = GNAT_d_odds)) + 
  geom_point(stat = 'identity', position = position_dodge(.125), size = 2) +
  geom_line(stat = 'identity', position = position_dodge(.125), linetype = 2) +
 geom_errorbar(aes(ymin = GNAT_d_odds_lci, ymax = GNAT_d_odds_uci), width = .2, position=position_dodge(.125), size = 0.8) + coord_flip() +
 geom_hline(yintercept=1, linetype="dashed", color = "red") +
  ylab("Odds") + 
  xlab ("Condition") + 
  theme_minimal()  +
  scale_x_discrete(limits = rev(levels(GNAT_conditions)))

GNAT_odds_CIS_graph <- GNAT_d_odds_df %>%
  ggplot(., aes(x = GNAT_conditions, y = GNAT_d_odds)) + 
  geom_point(stat = 'identity', position = position_dodge(.125), size = 2) +
  geom_line(stat = 'identity', position = position_dodge(.125), linetype = 2) +
  geom_errorbar(aes(ymin = GNAT_d_odds_lci, ymax = GNAT_d_odds_uci), width = .2, position=position_dodge(.125), size = 0.8) + coord_flip() +
  geom_hline(yintercept=1, linetype="dashed", color = "red") +
  ylab("Odds of meeting primary diagnostic criteria for a MDE") + 
  xlab ("Condition") + 
  theme_minimal()  +
  scale_x_discrete(limits = rev(levels(GNAT_conditions)))

ggsave("GNAT_odds_CIS_graph.TIFF")


```

## Social Reinforcement Learning

### Bias Scores

```{r}
SEL_bias$primary_depression_CIS <- as.factor(SEL_bias$primary_depression_CIS)
SEL_bias_wide$primary_depression_CIS <- as.factor(SEL_bias_wide$primary_depression_CIS)

# Using only session 1 data as CIS only measured at time 1
SEL_bias_wide_CIS <- SEL_bias_wide %>%
  filter(session == 1)

SEL_bias_long_CIS <- SEL_bias %>%
  filter(session == 1)

#Descriptives
SEL_bias_long_CIS %>%
  tab_cols(Condition) %>%
  tab_rows(primary_depression_CIS) %>%
  tab_cells(bias) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Bias Scores by Referential Condition and MDE Diagnosis (1 = Yes)")

SEL_bias_long_CIS %>%
  group_by(primary_depression_CIS, Condition) %>%
  summarize(Mean = mean(bias),  SD = sd(bias), Minimum = min(bias), Maximum = max(bias)) %>%
  ggplot(., aes(x = Condition, y = Mean, fill = primary_depression_CIS)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Bias Scores (positive - negative)") + 
  xlab ("Condition") + 
  scale_fill_discrete(name = "MDE Diagnosis", labels = c("No", "Yes")) +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

#Regression Model

SEL_bias_CIS_model <- glm(primary_depression_CIS ~ Self + Friend + Stranger, data = SEL_bias_wide_CIS, family = "binomial")
jtools::summ(SEL_bias_CIS_model, confint = TRUE, digits = 3)
SEL_bias_CIS_model_odds <- logistic.display(SEL_bias_CIS_model, crude = FALSE)
SEL_bias_CIS_model_odds
SEL_bias_CIS_model_effects <- ggpredict(SEL_bias_CIS_model)
plot(SEL_bias_CIS_model_effects)

#Plot model with odds ratios
##Creating a variable to show conditions
SEL_conditions <- c(1, 2, 3)
SEL_conditions <- factor(SEL_conditions,
                          levels = c(1, 2, 3),
                          labels = c("Self", "Friend", "Stranger"))
##Creating a variable for odds
SEL_bias_odds <- c(as.numeric(str_extract(SEL_bias_CIS_model_odds$table[1], "[[:digit:]]\\.[[:digit:]][[:digit:]]")), as.numeric(str_extract(SEL_bias_CIS_model_odds$table[3], "[[:digit:]]\\.[[:digit:]][[:digit:]]")), as.numeric(str_extract(SEL_bias_CIS_model_odds$table[5], "[[:digit:]]\\.[[:digit:]][[:digit:]]")))

##Creating a variable for lower CI
SEL_bias_odds_lci_1 <- as.numeric(gsub("\\(", "", str_extract(SEL_bias_CIS_model_odds$table[1], "\\([[:digit:]]\\.[[:digit:]][[:digit:]]")))
SEL_bias_odds_lci_2 <- as.numeric(gsub("\\(", "", str_extract(SEL_bias_CIS_model_odds$table[3], "\\([[:digit:]]\\.[[:digit:]][[:digit:]]")))
SEL_bias_odds_lci_3 <- as.numeric(gsub("\\(", "", str_extract(SEL_bias_CIS_model_odds$table[5], "\\([[:digit:]]\\.[[:digit:]][[:digit:]]")))
SEL_bias_odds_lci <- c(SEL_bias_odds_lci_1, SEL_bias_odds_lci_2, SEL_bias_odds_lci_3)
##Creating a variable for upper CI
SEL_bias_odds_uci_1 <- as.numeric(gsub("\\,", "", str_extract(SEL_bias_CIS_model_odds$table[1], "\\,[[:digit:]]\\.[[:digit:]][[:digit:]]")))
SEL_bias_odds_uci_2 <- as.numeric(gsub("\\,", "", str_extract(SEL_bias_CIS_model_odds$table[3], "\\,[[:digit:]]\\.[[:digit:]][[:digit:]]")))
SEL_bias_odds_uci_3 <- as.numeric(gsub("\\,", "", str_extract(SEL_bias_CIS_model_odds$table[5], "\\,[[:digit:]]\\.[[:digit:]][[:digit:]]")))
SEL_bias_odds_uci <- c(SEL_bias_odds_uci_1, SEL_bias_odds_uci_2, SEL_bias_odds_uci_3)
##Binding
SEL_bias_odds_df <- data.frame(SEL_conditions, SEL_bias_odds, SEL_bias_odds_lci, SEL_bias_odds_uci)
##Graph
SEL_bias_odds_df %>%
  ggplot(., aes(x = SEL_conditions, y = SEL_bias_odds)) + 
  geom_point(stat = 'identity', position = position_dodge(.125), size = 2) +
  geom_line(stat = 'identity', position = position_dodge(.125), linetype = 2) +
 geom_errorbar(aes(ymin = SEL_bias_odds_lci, ymax = SEL_bias_odds_uci), width = .2, position=position_dodge(.125), size = 0.8) + coord_flip() +
 geom_hline(yintercept=1, linetype="dashed", color = "red") +
  ylab("Odds") + 
  xlab ("Condition") + 
  theme_minimal() +
  scale_x_discrete(limits = rev(levels(SEL_conditions)))


SEL_odds_CIS_graph <- SEL_bias_odds_df %>%
  ggplot(., aes(x = SEL_conditions, y = SEL_bias_odds)) + 
  geom_point(stat = 'identity', position = position_dodge(.125), size = 2) +
  geom_line(stat = 'identity', position = position_dodge(.125), linetype = 2) +
  geom_errorbar(aes(ymin = SEL_bias_odds_lci, ymax = SEL_bias_odds_uci), width = .2, position=position_dodge(.125), size = 0.8) + coord_flip() +
  geom_hline(yintercept=1, linetype="dashed", color = "red") +
  ylab("Odds of meeting primary diagnostic criteria for a MDE") + 
  xlab ("Condition") + 
  theme_minimal() +
  scale_x_discrete(limits = rev(levels(SEL_conditions)))

ggsave("SEL_odds_CIS_graph.TIFF")


```

### Errors to Criterion

```{r}
SEL_long$primary_depression_CIS <- as.factor(SEL_long$primary_depression_CIS)
SEL_wide$primary_depression_CIS <- as.factor(SEL_wide$primary_depression_CIS)


#Keeping only session 1 as CIS only completed at this timepoint
SEL_long_CIS <- SEL_long %>%
  filter(session == 1)

SEL_wide_CIS <- SEL_wide %>%
  filter(session == 1)

#Descriptives
SEL_long_CIS %>%
  tab_cols(Condition %nest% Overall_Rule) %>%
  tab_rows(primary_depression_CIS) %>%
  tab_cells(ERR_to_crit) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Bias Scores by Referential Condition, Rule and MDE Diagnosis (1 = Yes)")

SEL_long_CIS %>%
  group_by(primary_depression_CIS, Condition, Overall_Rule) %>%
  summarize(Mean = mean(ERR_to_crit),  SD = sd(ERR_to_crit), Minimum = min(ERR_to_crit), Maximum = max(ERR_to_crit)) %>%
  ggplot(., aes(x = Overall_Rule, y = Mean, fill = primary_depression_CIS)) + 
  facet_wrap(~Condition) +
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Errors to Criterion") + 
  xlab ("Overall_Rule") + 
  scale_fill_discrete(name = "MDE Diagnosis", labels = c("No", "Yes")) +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

#Regression Model
SEL_err_CIS_model <- glm(primary_depression_CIS ~ ERR_to_crit_Self_Positive + ERR_to_crit_Self_Negative + ERR_to_crit_Friend_Positive + ERR_to_crit_Friend_Negative + ERR_to_crit_Stranger_Positive + ERR_to_crit_Stranger_Negative, data = SEL_wide_CIS, family = "binomial")
jtools::summ(SEL_err_CIS_model, confint = TRUE, digits = 3)
logistic.display(SEL_err_CIS_model, crude = FALSE)
SEL_err_CIS_model_effects <- ggpredict(SEL_err_CIS_model)
plot(SEL_err_CIS_model_effects)

Self_Positive_effects <- as.data.frame(SEL_err_CIS_model_effects$ERR_to_crit_Self_Positive)
Self_Negative_effects <- as.data.frame(SEL_err_CIS_model_effects$ERR_to_crit_Self_Negative)
Friend_Positive_effects <- as.data.frame(SEL_err_CIS_model_effects$ERR_to_crit_Friend_Positive)
Friend_Negative_effects <- as.data.frame(SEL_err_CIS_model_effects$ERR_to_crit_Friend_Negative)
Stranger_Positive_effects <- as.data.frame(SEL_err_CIS_model_effects$ERR_to_crit_Stranger_Positive)
Stranger_Negative_effects <- as.data.frame(SEL_err_CIS_model_effects$ERR_to_crit_Stranger_Negative)

SEL_err_CIS_model_effects_graph <- bind_rows(Self_Positive_effects, Self_Negative_effects, Friend_Positive_effects, Friend_Negative_effects, Stranger_Positive_effects, Stranger_Negative_effects)

SEL_err_CIS_model_effects_graph <- SEL_err_CIS_model_effects_graph %>%
  mutate(ConditionRule = ifelse(group == "ERR_to_crit_Self_Positive", 1, ifelse(group == "ERR_to_crit_Self_Negative", 2, ifelse(group == "ERR_to_crit_Friend_Positive", 3, ifelse(group == "ERR_to_crit_Friend_Negative", 4, ifelse(group == "ERR_to_crit_Stranger_Positive", 5, 6)))))) %>%
  mutate(Condition = ifelse(ConditionRule <=2, 1, ifelse(ConditionRule > 2 & ConditionRule <= 4, 2, 3))) %>%
  mutate(Rule = ifelse(ConditionRule == 1 | ConditionRule == 3 | ConditionRule == 5, 1, 2))
                                 
SEL_err_CIS_model_effects_graph$Condition <- factor(SEL_err_CIS_model_effects_graph$Condition,
                                       levels = c(1, 2, 3),
                                       labels = c("Self", "Friend", "Stranger"))

SEL_err_CIS_model_effects_graph$Rule <- factor(SEL_err_CIS_model_effects_graph$Rule,
                                       levels = c(1, 2),
                                       labels = c("Positive", "Negative"))


ggplot(SEL_err_CIS_model_effects_graph, aes(x = x, y = predicted, color = Rule)) + facet_wrap(~ Condition) + geom_line(size = 1) + geom_ribbon(data = SEL_err_CIS_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Errors to Criterion") + ylab("Probability") + ggtitle("Linear Regression SEL Errors to Criterion on Probability of MDE Diagnosis")

```

### Global Ratings

```{r}
SEL_long$primary_depression_CIS <- as.factor(SEL_long$primary_depression_CIS)
SEL_wide$primary_depression_CIS <- as.factor(SEL_wide$primary_depression_CIS)

#Keeping only session 1 as CIS only completed at this timepoint
SEL_long_CIS <- SEL_long %>%
  filter(session == 1)

SEL_wide_CIS <- SEL_wide %>%
  filter(session == 1)

#Descriptives
SEL_long_CIS %>%
  tab_cols(Condition %nest% Overall_Rule) %>%
  tab_rows(primary_depression_CIS) %>%
  tab_cells(globalrating) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  set_caption("SEL Bias Scores by Referential Condition, Rule and MDE Diagnosis (1 = Yes)")

SEL_long_CIS %>%
  group_by(primary_depression_CIS, Condition, Overall_Rule) %>%
  summarize(Mean = mean(globalrating),  SD = sd(globalrating), Minimum = min(globalrating), Maximum = max(globalrating)) %>%
  ggplot(., aes(x = Overall_Rule, y = Mean, fill = primary_depression_CIS)) + 
  facet_wrap(~Condition) +
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Errors to Criterion") + 
  xlab ("Overall_Rule") + 
  scale_fill_discrete(name = "MDE Diagnosis", labels = c("No", "Yes")) +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

#Regression Model
SEL_global_CIS_model <- glm(primary_depression_CIS ~ globalrating_Self_Positive + globalrating_Self_Negative + globalrating_Friend_Positive + globalrating_Friend_Negative + globalrating_Stranger_Positive + globalrating_Stranger_Negative, data = SEL_wide_CIS, family = "binomial")
jtools::summ(SEL_global_CIS_model, confint = TRUE, digits = 3)
logistic.display(SEL_global_CIS_model, crude = FALSE)
SEL_global_CIS_model_effects <- ggpredict(SEL_global_CIS_model)
plot(SEL_global_CIS_model_effects)

Self_Positive_effects <- as.data.frame(SEL_global_CIS_model_effects$globalrating_Self_Positive)
Self_Negative_effects <- as.data.frame(SEL_global_CIS_model_effects$globalrating_Self_Negative)
Friend_Positive_effects <- as.data.frame(SEL_global_CIS_model_effects$globalrating_Friend_Positive)
Friend_Negative_effects <- as.data.frame(SEL_global_CIS_model_effects$globalrating_Friend_Negative)
Stranger_Positive_effects <- as.data.frame(SEL_global_CIS_model_effects$globalrating_Stranger_Positive)
Stranger_Negative_effects <- as.data.frame(SEL_global_CIS_model_effects$globalrating_Stranger_Negative)

SEL_global_CIS_model_effects_graph <- bind_rows(Self_Positive_effects, Self_Negative_effects, Friend_Positive_effects, Friend_Negative_effects, Stranger_Positive_effects, Stranger_Negative_effects)

SEL_global_CIS_model_effects_graph <- SEL_global_CIS_model_effects_graph %>%
  mutate(ConditionRule = ifelse(group == "globalrating_Self_Positive", 1, ifelse(group == "globalrating_Self_Negative", 2, ifelse(group == "globalrating_Friend_Positive", 3, ifelse(group == "globalrating_Friend_Negative", 4, ifelse(group == "globalrating_Stranger_Positive", 5, 6)))))) %>%
  mutate(Condition = ifelse(ConditionRule <=2, 1, ifelse(ConditionRule > 2 & ConditionRule <= 4, 2, 3))) %>%
  mutate(Rule = ifelse(ConditionRule == 1 | ConditionRule == 3 | ConditionRule == 5, 1, 2))
                                 
SEL_global_CIS_model_effects_graph$Condition <- factor(SEL_global_CIS_model_effects_graph$Condition,
                                       levels = c(1, 2, 3),
                                       labels = c("Self", "Friend", "Stranger"))

SEL_global_CIS_model_effects_graph$Rule <- factor(SEL_global_CIS_model_effects_graph$Rule,
                                       levels = c(1, 2),
                                       labels = c("Positive", "Negative"))


ggplot(SEL_global_CIS_model_effects_graph, aes(x = x, y = predicted, color = Rule)) + facet_wrap(~ Condition) + geom_line(size = 1) + geom_ribbon(data = SEL_global_CIS_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Global Rating") + ylab("Probability") + ggtitle("Linear Regression SEL Global Ratings on Probability of MDE Diagnosis")

```

# Adjusting for gender & age

## Associative Learning

### Depression Models

#### PHQ

```{r}
#Self
associative_wide_matching_collapsed_self <- associative_wide_matching_collapsed %>%
  filter(Task == "Self") #Restricting Task to "Self"

##Accuracy
assoc_self_acc_age_gender_PHQ_model <- lm(PHQ_tot ~ prop_acc_High + prop_acc_Medium + prop_acc_Low + gender + age, data = associative_wide_matching_collapsed_self)
jtools::summ(assoc_self_acc_age_gender_PHQ_model, confint = TRUE, digits = 3)
assoc_self_acc_age_gender_PHQ_model_standard <- standardize_parameters(assoc_self_acc_age_gender_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_self_acc_age_gender_PHQ_model_standard



##RT
assoc_self_rt_age_gender_PHQ_model <- lm(PHQ_tot ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low + gender + age, data = associative_wide_matching_collapsed_self)
jtools::summ(assoc_self_rt_age_gender_PHQ_model, confint = TRUE, digits = 3)
assoc_self_rt_age_gender_PHQ_model_standard <- standardize_parameters(assoc_self_rt_age_gender_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_self_rt_age_gender_PHQ_model_standard


#Reward
associative_wide_matching_collapsed_reward <- associative_wide_matching_collapsed %>%
  filter(Task == "Reward") #Restricting Task to "Reward"

##Accuracy
assoc_reward_acc_age_gender_PHQ_model <- lm(PHQ_tot ~ prop_acc_High + prop_acc_Medium + prop_acc_Low + gender + age, data = associative_wide_matching_collapsed_reward)
jtools::summ(assoc_reward_acc_age_gender_PHQ_model, confint = TRUE, digits = 3)
assoc_reward_acc_age_gender_PHQ_model_standard <- standardize_parameters(assoc_reward_acc_age_gender_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_reward_acc_age_gender_PHQ_model_standard

## RT
assoc_reward_rt_age_gender_PHQ_model <- lm(PHQ_tot ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low + gender + age, data = associative_wide_matching_collapsed_reward)
jtools::summ(assoc_reward_rt_age_gender_PHQ_model, confint = TRUE, digits = 3)
assoc_reward_rt_age_gender_PHQ_model_standard <- standardize_parameters(assoc_reward_rt_age_gender_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_reward_rt_age_gender_PHQ_model_standard

#Emotion/Valence
associative_wide_matching_collapsed_valence <- associative_wide_matching_collapsed %>%
  filter(Task == "Valence") #Restricting Task to "Valence"

##Accuracy
assoc_valence_acc_age_gender_PHQ_model <- lm(PHQ_tot ~ prop_acc_High + prop_acc_Medium + prop_acc_Low + gender + age, data = associative_wide_matching_collapsed_valence)
jtools::summ(assoc_valence_acc_age_gender_PHQ_model, confint = TRUE, digits = 3)

assoc_valence_acc_age_gender_PHQ_model_standard <- standardize_parameters(assoc_valence_acc_age_gender_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_valence_acc_age_gender_PHQ_model_standard

##RT
assoc_valence_rt_age_gender_PHQ_model <- lm(PHQ_tot ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low + gender + age, data = associative_wide_matching_collapsed_valence)
jtools::summ(assoc_valence_rt_age_gender_PHQ_model, confint = TRUE, digits = 3)

assoc_valence_rt_age_gender_PHQ_model_standard <- standardize_parameters(assoc_valence_rt_age_gender_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_valence_rt_age_gender_PHQ_model_standard
```
#### BDI

```{r}
#Self
associative_wide_matching_collapsed_self <- associative_wide_matching_collapsed %>%
  filter(Task == "Self") #Restricting Task to "Self"

##Accuracy
assoc_self_acc_age_gender_BDI_model <- lm(BDI_tot ~ prop_acc_High + prop_acc_Medium + prop_acc_Low + gender + age, data = associative_wide_matching_collapsed_self)
jtools::summ(assoc_self_acc_age_gender_BDI_model, confint = TRUE, digits = 3)
assoc_self_acc_age_gender_BDI_model_standard <- standardize_parameters(assoc_self_acc_age_gender_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_self_acc_age_gender_BDI_model_standard



##RT
assoc_self_rt_age_gender_BDI_model <- lm(BDI_tot ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low + gender + age, data = associative_wide_matching_collapsed_self)
jtools::summ(assoc_self_rt_age_gender_BDI_model, confint = TRUE, digits = 3)
assoc_self_rt_age_gender_BDI_model_standard <- standardize_parameters(assoc_self_rt_age_gender_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_self_rt_age_gender_BDI_model_standard


#Reward
associative_wide_matching_collapsed_reward <- associative_wide_matching_collapsed %>%
  filter(Task == "Reward") #Restricting Task to "Reward"

##Accuracy
assoc_reward_acc_age_gender_BDI_model <- lm(BDI_tot ~ prop_acc_High + prop_acc_Medium + prop_acc_Low + gender + age, data = associative_wide_matching_collapsed_reward)
jtools::summ(assoc_reward_acc_age_gender_BDI_model, confint = TRUE, digits = 3)
assoc_reward_acc_age_gender_BDI_model_standard <- standardize_parameters(assoc_reward_acc_age_gender_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_reward_acc_age_gender_BDI_model_standard

## RT
assoc_reward_rt_age_gender_BDI_model <- lm(BDI_tot ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low + gender + age, data = associative_wide_matching_collapsed_reward)
jtools::summ(assoc_reward_rt_age_gender_BDI_model, confint = TRUE, digits = 3)
assoc_reward_rt_age_gender_BDI_model_standard <- standardize_parameters(assoc_reward_rt_age_gender_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_reward_rt_age_gender_BDI_model_standard

#Emotion/Valence
associative_wide_matching_collapsed_valence <- associative_wide_matching_collapsed %>%
  filter(Task == "Valence") #Restricting Task to "Valence"

##Accuracy
assoc_valence_acc_age_gender_BDI_model <- lm(BDI_tot ~ prop_acc_High + prop_acc_Medium + prop_acc_Low + gender + age, data = associative_wide_matching_collapsed_valence)
jtools::summ(assoc_valence_acc_age_gender_BDI_model, confint = TRUE, digits = 3)

assoc_valence_acc_age_gender_BDI_model_standard <- standardize_parameters(assoc_valence_acc_age_gender_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_valence_acc_age_gender_BDI_model_standard

##RT
assoc_valence_rt_age_gender_BDI_model <- lm(BDI_tot ~ mean_rt_mult_High + mean_rt_mult_Medium + mean_rt_mult_Low + gender + age, data = associative_wide_matching_collapsed_valence)
jtools::summ(assoc_valence_rt_age_gender_BDI_model, confint = TRUE, digits = 3)

assoc_valence_rt_age_gender_BDI_model_standard <- standardize_parameters(assoc_valence_rt_age_gender_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
assoc_valence_rt_age_gender_BDI_model_standard
```
## GNAT

### PHQ-9


```{r}

GNAT_discrim_gender_age_PHQ_model <- lm(PHQ_tot ~ d_Self_Pos + d_Self_Neg + d_Other_Pos + d_Other_Neg + gender + age, data = GNAT_wide_Qs_eligible)
jtools::summ(GNAT_discrim_gender_age_PHQ_model, confint = TRUE, digits = 3)

GNAT_discrim_gender_age_PHQ_model_standard <- standardize_parameters(GNAT_discrim_gender_age_PHQ_model, method = "refit", ci = 0.95, include_response = TRUE)
GNAT_discrim_gender_age_PHQ_model_standard

```

### BDI-II

```{r}
GNAT_discrim_gender_age_BDI_model <- lm(BDI_tot ~ d_Self_Pos + d_Self_Neg + d_Other_Pos + d_Other_Neg + gender + age, data = GNAT_wide_Qs_eligible)
jtools::summ(GNAT_discrim_gender_age_BDI_model, confint = TRUE, digits = 3)

GNAT_discrim_gender_age_BDI_model_standard <- standardize_parameters(GNAT_discrim_gender_age_BDI_model, method = "refit", ci = 0.95, include_response = TRUE)
GNAT_discrim_gender_age_BDI_model_standard

```



## Social Reinforcement Learning

### Bias Scores

#### PHQ
```{r}
SEL_bias_wide <- SEL_bias %>%
  dplyr::select(subject, session, Condition, age, gender, bias, BDI_tot, PHQ_tot, BFNE_tot, primary_depression_CIS) %>%
  spread(key = Condition, value = bias) 

SEL_bias_PHQ_age_gender_model <- lmer(PHQ_tot ~ Self + Friend + Stranger + gender + age + (1|subject), data = SEL_bias_wide) 
jtools::summ(SEL_bias_PHQ_age_gender_model, confint = TRUE, digits = 3)

SEL_bias_PHQ_age_gender_model_standard <- standardize_parameters(SEL_bias_PHQ_age_gender_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_bias_PHQ_age_gender_model_standard
```


#### BDI
```{r}
SEL_bias_wide <- SEL_bias %>%
  dplyr::select(subject, session, Condition, bias, BDI_tot, PHQ_tot, BFNE_tot, primary_depression_CIS, gender, age) %>%
  spread(key = Condition, value = bias) 

SEL_bias_BDI_age_gender_model <- lmer(BDI_tot ~ Self + Friend + Stranger + gender + age + (1|subject), data = SEL_bias_wide) 
jtools::summ(SEL_bias_BDI_age_gender_model, confint = TRUE, digits = 3)

SEL_bias_BDI_age_gender_model_standard <- standardize_parameters(SEL_bias_BDI_age_gender_model, method = "refit", ci = 0.95, include_response = TRUE)
SEL_bias_BDI_age_gender_model_standard
```


### Errors to Criterion

#### PHQ
```{r}
SEL_err_PHQ_age_gender_model <- lmer(PHQ_tot ~ ERR_to_crit_Self_Positive + ERR_to_crit_Self_Negative + ERR_to_crit_Friend_Positive + ERR_to_crit_Friend_Negative + ERR_to_crit_Stranger_Positive + ERR_to_crit_Stranger_Negative + age + gender + (1|subject), data = SEL_wide) 
jtools::summ(SEL_err_PHQ_age_gender_model, confint = TRUE, digits = 3)
```


#### BDI
```{r}
SEL_err_BDI_age_gender_model <- lmer(BDI_tot ~ ERR_to_crit_Self_Positive + ERR_to_crit_Self_Negative + ERR_to_crit_Friend_Positive + ERR_to_crit_Friend_Negative + ERR_to_crit_Stranger_Positive + ERR_to_crit_Stranger_Negative + age + gender + (1|subject), data = SEL_wide) 
jtools::summ(SEL_err_BDI_age_gender_model, confint = TRUE, digits = 3)
```


### Global Ratings

#### PHQ
```{r}
SEL_global_PHQ_age_gender_model <- lmer(PHQ_tot ~ globalrating_Self_Positive + globalrating_Self_Negative + globalrating_Friend_Positive + globalrating_Friend_Negative + globalrating_Stranger_Positive + globalrating_Stranger_Negative + age + gender + (1|subject), data = SEL_wide) 
jtools::summ(SEL_global_PHQ_age_gender_model, confint = TRUE, digits = 3)
```

#### BDI
```{r}
SEL_global_BDI_age_gender_model <- lmer(BDI_tot ~ globalrating_Self_Positive + globalrating_Self_Negative + globalrating_Friend_Positive + globalrating_Friend_Negative + globalrating_Stranger_Positive + globalrating_Stranger_Negative + age + gender + (1|subject), data = SEL_wide) 
jtools::summ(SEL_global_BDI_age_gender_model, confint = TRUE, digits = 3)
```


