---
title: "Clinical Study Analysis"
output: html_notebook
---


```{r}
library(tidyverse)#General cleaning
library(expss) #Descriptive Tables
expss_output_rnotebook()
expss_digits(digits = 2)
library(lmerTest) #MLMs
library(jtools) #Summary Tables for Regression
library(ggeffects) #Regression Graphs
library(effectsize) #Standardising parameters
library(ggpubr)#Merging plots
library(grid)
```


```{r}
# Importing Files
load("SEL/SEL_bias.RData")
load("Questionnaires/All_Questionnaires.RData")

#Merging SEL & Q Data
SEL_bias_Qs <- merge(SEL_bias, All_Questionnaires, by = c("subject", "session"), all = TRUE)
```

# Number of participants at each timepoint
```{r}
# Identifying number of datapoints completed
key_data_long <- SEL_bias_Qs %>%
  select(subject, session, Condition, Positive, Negative, PHQ_tot, BDI_tot) %>% #Pulling out key data
  group_by(subject, session) %>%
  summarise(Positive = mean(Positive, na.rm = TRUE), Negative = mean(Negative, na.rm = TRUE), PHQ_tot = mean(PHQ_tot, na.rm = TRUE), BDI_tot = mean(BDI_tot, na.rm = TRUE)) %>% #collapsing down by condition so that if there was any data in the SEL it will be provided in the summary statistics
  mutate(SEL_missing = ifelse(is.na(Positive) & is.na(Negative), 1, 0)) %>%
  mutate(Qs_missing = ifelse(is.na(PHQ_tot) & is.na (BDI_tot), 1, 0)) %>%
  mutate(data_available = ifelse(SEL_missing == 0 | Qs_missing == 0, 1, 0)) %>%
  group_by(subject) %>%
  mutate(total_datapoints = sum(data_available)) %>%
  mutate(baseline_data_available = ifelse(session == "1" & data_available == 1, 1, 0)) %>%
  mutate(drop_out = ifelse(baseline_data_available == 1 & total_datapoints == 1, 1, 0)) %>%
  ungroup()

#Identifying p.s with only baseline data
drop_out_identifiers <- key_data_long %>%
  select(subject, session, drop_out) %>%
  filter(session == "1") %>%
  select(-session)

#Calculating Ns & %s of p.s at each timepoint
participant_n <- key_data_long %>%
  select(subject, session, data_available) %>%
  pivot_wider(names_from = session, names_prefix = "data_available_session_", values_from = "data_available")

participant_n_summarised <- participant_n %>%
  summarise_at(c("data_available_session_1", "data_available_session_2", "data_available_session_3", "data_available_session_4", "data_available_session_5"), sum, na.rm = TRUE) %>%
  pivot_longer(cols = starts_with("data_available_session_"),
                                  names_to = "session",
                                  values_to = "N") %>%
  mutate(session = gsub(pattern = "data_available_session_", replacement = "", x = session)) %>%
  mutate(percentage = round((N/max(N)*100), digits = 2))
  
print(participant_n_summarised)

```

# Demographics & Clinical Characteristics
```{r}
demographics_clinical <- All_Questionnaires %>%
  filter(session == "1") %>%
  select(subject, session, age, gender, ethnicity, employment,relationship, education, housing, depress_dur, antid_curr, antid_curr_type, antid_length_baseline, psych_curr, other_meds, antid_past, psych_past)

#Sample Characteristics 
demographics_clinical %>%
  tab_cells(age) %>%
  tab_stat_mean_sd_n() %>%
  tab_stat_min() %>%
  tab_stat_max() %>%
  tab_cells(gender, ethnicity, employment, education, relationship, housing, depress_dur, antid_curr, antid_curr_type) %>%
  tab_stat_cases(label = "N", total_row_position = "none") %>%
  tab_stat_cpct(label = "%", total_row_position = "none") %>%
  tab_cells(antid_length_baseline) %>%
  tab_stat_mean_sd_n() %>%
  tab_cells(other_meds, antid_past, psych_past) %>%
    tab_stat_cases(label = "N", total_row_position = "none") %>%
  tab_stat_cpct(label = "%", total_row_position = "none") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Demographics")


CIS <- All_Questionnaires %>%
  filter(session == "1") %>%
  select(CIS_primary, CIS_secondary)

CIS %>%
  tab_cells(CIS_primary, CIS_secondary) %>%
  tab_stat_cases(label = "N", total_row_position = "none") %>%
  tab_stat_cpct(label = "%", total_row_position = "none") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Demographics")
  

```
# Time varying treatment characteristics
```{r}
treatment_by_time <- All_Questionnaires %>%
  select(subject, session, antid_curr, antid_curr_change, antid_curr_type, antid_curr_type_change, antid_curr_dose, antid_curr_dose_change, antid_any_change, antid_curr_adherence, antid_side_eff, antid_side_eff_freq, antid_side_eff_sev, other_meds, other_meds_type, psych_curr)

dose_change <- treatment_by_time %>%
  select(subject, session, antid_curr_type, antid_curr_dose, antid_curr_dose_change) %>%
  group_by(subject) %>%
  filter(any(antid_curr_dose_change== "Change")) %>%
  ungroup()

meds_change <- treatment_by_time %>%
  select(subject, session, antid_curr_type, antid_curr_type_change) %>%
  group_by(subject) %>%
  filter(any(antid_curr_type_change== "Change")) %>%
  ungroup()

All_Questionnaires %>%
  tab_cols(session) %>%
  tab_cells(antid_curr, antid_curr_type, antid_curr_change, antid_curr_type_change, antid_curr_dose_change, antid_curr_adherence, antid_side_eff, antid_side_eff_freq, antid_side_eff_sev, psych_curr, other_meds, data_collection) %>%
  tab_stat_cases(label = "N", total_row_position = "none") %>%
  tab_stat_cpct(label = "%", total_row_position = "none") %>%
  tab_pivot(stat_position = "inside_rows")%>%
  set_caption("Treatment Characteristics by timepoint")
```


# Descriptives
```{r}
# Creating a binary GCR - better vs. same/worse
SEL_bias_Qs <- SEL_bias_Qs %>%
  mutate(GRC_binary = ifelse(GRC < 3, 2, 1)) %>%
  mutate(GRC_binary_factor = factor(GRC_binary, levels = c(1, 2), labels = c("Same/Worse", "Better")))

# Table
descriptives_table <- SEL_bias_Qs %>%
  select(subject, session, Condition, bias, PHQ_tot, BDI_tot, GAD_tot, GRC_binary_factor) %>%
  pivot_wider(values_from = bias, names_from = Condition) %>% mutate(session_FU = as.numeric(session)) %>%
  mutate(session_FU = factor(session, levels = c(1, 2, 3, 4, 5), labels = c("Baseline", "2 Weeks", "6 Weeks", "8 Weeks", "6 Months")))

descriptives_table %>%
  tab_cols(session_FU) %>%
  tab_cells(PHQ_tot, BDI_tot, GAD_tot, Self, Friend, Stranger) %>%
  tab_stat_mean() %>%
  tab_stat_sd() %>%
  tab_stat_valid_n() %>%
  tab_pivot(stat_position = "inside_columns")%>%
  set_caption("Mean outcomes by timepoint")

descriptives_table %>%
  tab_cols(session_FU) %>%
  tab_cells(GRC_binary_factor) %>%
  tab_stat_cases() %>%
  tab_stat_cpct() %>%
    tab_pivot(stat_position = "inside_columns")%>%
  set_caption("GRC Responses by timepoint")

#Graphs
SEL_bias_Qs_short <- SEL_bias_Qs %>%
  filter(session != "5") #Not including 6 month FU

### PHQ Scores
PHQ_total_graph <- SEL_bias_Qs_short %>%
  select(subject, session, PHQ_tot) %>%
  mutate(session_FU = as.numeric(session)) %>%
  mutate(session_FU = factor(session, levels = c(1, 2, 3, 4), labels = c("Baseline", "2 Weeks", "6 Weeks", "8 Weeks"))) %>%
  group_by(session_FU) %>%
  summarise( 
    n=n(),
    mean=mean(PHQ_tot, na.rm = TRUE),
    sd=sd(PHQ_tot, na.rm = TRUE)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
  ungroup()
PHQ_total_graph %>%
  ggplot(., aes(x = session_FU, y = mean, group = 1)) + 
  geom_point(stat = 'identity', size = 2) +
  geom_line(stat = 'identity', linetype = 2) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2, position=position_dodge(.2), size = 0.8) +
  ylab("PHQ Scores") + 
  xlab ("Timepoint") + 
  theme_minimal()

### BDI Scores
BDI_total_graph <- SEL_bias_Qs_short %>%
  select(subject, session, BDI_tot) %>%
  mutate(session_FU = as.numeric(session)) %>%
  mutate(session_FU = factor(session, levels = c(1, 2, 3, 4), labels = c("Baseline", "2 Weeks", "6 Weeks", "8 Weeks"))) %>%
  group_by(session_FU) %>%
  summarise( 
    n=n(),
    mean=mean(BDI_tot, na.rm = TRUE),
    sd=sd(BDI_tot, na.rm = TRUE)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
  ungroup()
BDI_total_graph %>%
  ggplot(., aes(x = session_FU, y = mean, group = 1)) + 
  geom_point(stat = 'identity', size = 2) +
  geom_line(stat = 'identity', linetype = 2) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2, position=position_dodge(.2), size = 0.8) +
  ylab("BDI Scores") + 
  xlab ("Timepoint") + 
  theme_minimal()

#GAD
GAD_total_graph <- SEL_bias_Qs_short %>%
  select(subject, session, GAD_tot) %>%
  mutate(session_FU = as.numeric(session)) %>%
  mutate(session_FU = factor(session, levels = c(1, 2, 3, 4), labels = c("Baseline", "2 Weeks", "6 Weeks", "8 Weeks"))) %>%
  group_by(session_FU) %>%
  summarise( 
    n=n(),
    mean=mean(GAD_tot, na.rm = TRUE),
    sd=sd(GAD_tot, na.rm = TRUE)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
  ungroup()
GAD_total_graph %>%
  ggplot(., aes(x = session_FU, y = mean, group = 1)) + 
  geom_point(stat = 'identity', size = 2) +
  geom_line(stat = 'identity', linetype = 2) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2, position=position_dodge(.2), size = 0.8) +
  ylab("GAD Scores") + 
  xlab ("Timepoint") + 
  theme_minimal()

##SEL Bias Scores by Condition

SEL_bias_graph <- SEL_bias_Qs_short %>%
  select(subject, session, Condition, bias) %>%
  mutate(session_FU = as.numeric(session)) %>%
  mutate(session_FU = factor(session, levels = c(1, 2, 3, 4), labels = c("Baseline", "2 Weeks", "6 Weeks", "8 Weeks"))) %>%
  group_by(session_FU, Condition) %>%
  summarise( 
    n=n(),
    mean=mean(bias, na.rm = TRUE),
    sd=sd(bias, na.rm = TRUE)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
  ungroup()


SEL_bias_graph %>%
  ggplot(., aes(x = session_FU, y = mean, colour = Condition, group = Condition)) + 
  geom_point(stat = 'identity', position = position_dodge(.2), size = 2) +
  geom_line(stat = 'identity', position = position_dodge(.2), linetype = 2) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2, position=position_dodge(.2), size = 0.8) +
  ylab("Bias Scores") + 
  xlab ("Timepoint") + 
  labs(colour = "Condition") +
  theme_minimal(base_size = 10)


## Self Bias Scores
Self_SEL_bias_graph <- SEL_bias_graph %>%
  filter(Condition == "Self")
Self_SEL_bias_graph %>%
  ggplot(., aes(x = session_FU, y = mean, group = 1)) + 
  geom_point(stat = 'identity', size = 2) +
  geom_line(stat = 'identity', linetype = 2) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2, position=position_dodge(.2), size = 0.8) +
  ylab("Self Bias Scores") + 
  xlab ("Timepoint") + 
  theme_minimal()
```


# Hypothesis 1
Social evaluation learning regarding the self will become more positively biased during antidepressant treatment, as measured by a better learning of positive evaluations relative to negative evaluations towards the self. 

A multi-level linear regression model will be used to evaluate changes in social evaluation learning. Bias scores, calculated through errors to criterion, will be entered as a continuous outcome, and time as a categorical fixed effect predictor with baseline as the reference category. To examine whether changes in learning are specific to the self, referential condition and the interaction of this with time, will be entered as additional categorical predictors. Participant ID will be entered as a random effect to account for clustering within repeated measures

```{r}
hyp1 <- lmer(bias ~ session*Condition + (1|subject), data = SEL_bias_Qs_short)
summ(hyp1, confint = TRUE, digits = 3)
hyp1_model_effects <- ggpredict(hyp1, terms = c("session", "Condition"))
plot(hyp1_model_effects)

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

#Getting overall p-value for main effect of session
hyp1_no_session <- lmer(bias ~ Condition + (1|subject), data = SEL_bias_Qs_short)
anova(hyp1, hyp1_no_session)

#Getting overall p-value for main effect of condition
hyp1_no_condition <- lmer(bias ~ session + (1|subject), data = SEL_bias_Qs_short)
anova(hyp1, hyp1_no_condition)

#Getting overall p-value for interaction effect of referential condition
hyp1_no_int <- lmer(bias ~ session + Condition + (1|subject), data = SEL_bias_Qs_short)
anova(hyp1, hyp1_no_int)
```


# Hypothesis 2
Changes in social evaluation learning will predict a reduction in depressive symptoms, as indicated by a decrease in PHQ-9 scores

A multi-level linear regression will be used to evaluate the relationship between change in social evaluation learning biases towards the self and change in depression severity over the first 8-weeks of antidepressant treatment. Change in PHQ-9 scores will be entered as the outcome and change in bias based on errors to criterion in the self and other conditions will be entered as exposures in wide format. Time-point will be entered as a fixed effect, and participant ID as a random effect, to allow for clustering over time-points. Depression severity will be adjusted for by entering baseline PHQ-9 scores as a fixed effect. This model will be repeated with change in BDI-II scores as the outcome.


```{r}
## Creating change score dataframe
bias_change <- SEL_bias_Qs %>%
  filter(session != "5") %>% #Removing 6 month FU
  select(subject, session, Condition, bias, PHQ_tot, BDI_tot, GAD_tot, GRC_binary_factor) %>%
  pivot_wider(names_from = "Condition", values_from = "bias") %>%
  arrange(subject, session) %>%
  group_by(subject) %>%
  mutate(Self_bias_change = Self - lag(Self),  Friend_bias_change = Friend - lag(Friend), Stranger_bias_change = Stranger - lag(Stranger), PHQ_change =  PHQ_tot - lag(PHQ_tot), BDI_change = BDI_tot - lag(BDI_tot), GAD_change = GAD_tot - lag(GAD_tot)) %>%
  ungroup()
  
baseline <- bias_change %>%
  filter(session == 1) %>%
  select(subject, PHQ_tot, BDI_tot, GAD_tot, Self, Friend, Stranger) %>%
  rename(PHQ_tot_baseline = PHQ_tot, BDI_tot_baseline = BDI_tot, GAD_tot_baseline = GAD_tot, Self_bias_baseline = Self, Friend_bias_baseline = Friend, Stranger_bias_baseline = Stranger)

bias_change <- bias_change %>%
  filter(session != "1") %>%
  select(subject, session, Self_bias_change, Friend_bias_change, Stranger_bias_change, PHQ_change, BDI_change, GAD_change, GRC_binary_factor) %>%
  mutate(change_timepoint = ifelse(session == "2", 1,
                                   ifelse(session == "3", 2, 3))) %>%
  mutate(change_timepoint = factor(change_timepoint, levels = c(1, 2, 3),
                                   labels = c("Baseline-2 weeks", "2-6 weeks", "6-8 weeks")))

bias_change <- merge(bias_change, baseline, by = "subject")


```

## Visualising Data

```{r}
#Prepping dataframe
bias_change_wide <- bias_change %>%
  pivot_longer(cols = c(Self_bias_change, Friend_bias_change, Stranger_bias_change), names_to = "condition", values_to = "bias_change") %>%
  mutate(Condition = ifelse(condition == "Self_bias_change", 1,
                            ifelse(condition == "Friend_bias_change", 2, 3))) %>%
           mutate(Condition = factor(Condition, levels = c(1, 2, 3), labels = c("Self", "Friend", "Stranger"))) %>%
  select(-condition)

# PHQ
PHQ_change_scatter <- ggplot(bias_change_wide, aes(x = bias_change, y = PHQ_change, colour = Condition)) +
  geom_point() +
  geom_smooth(method = "lm") +
  ylim(-17, 10) +
  facet_grid(. ~ Condition + . ~ change_timepoint) +
  xlab("Change in Bias Scores (Positive - Negative)") +
  ylab("Change in PHQ-9") +
  theme(strip.text.y = element_text(size = 8, colour = "black", angle = 270),
        strip.background = element_rect(colour = "black", fill = "white"),
        panel.border = element_rect(colour = "black", fill=NA, size=1))

# BDI
BDI_change_scatter <- ggplot(bias_change_wide, aes(x = bias_change, y = BDI_change, colour = Condition)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(. ~ Condition + . ~ change_timepoint) +
    xlab("Change in Bias Scores (Positive - Negative)") +
  ylab("Change in BDI-II") +
  theme(strip.text.y = element_text(size = 8, colour = "black", angle = 270),
        strip.background = element_rect(colour = "black", fill = "white"),
        panel.border = element_rect(colour = "black", fill=NA, size=1))


#GAD
GAD_change_scatter <- ggplot(bias_change_wide, aes(x = bias_change, y = GAD_change, colour = Condition)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(. ~ Condition + . ~ change_timepoint) +
  xlab("Change in Bias Scores (Positive - Negative)") +
  ylab("Change in GAD-7") +
  theme(strip.text.y = element_text(size = 8, colour = "black", angle = 270),
        strip.background = element_rect(colour = "black", fill = "white"),
        panel.border = element_rect(colour = "black", fill=NA, size=1))




#Bringing together plots
all_change_scatter <- ggarrange(PHQ_change_scatter, NULL, BDI_change_scatter, NULL, GAD_change_scatter, 
          labels = c("PHQ-9", "", "BDI-II", "", "GAD-7"), 
          font.label = list(size = 12, face = "bold"),
          ncol = 1, nrow = 5, heights = c(1, 0.1, 1, 0.1, 1),
          label.y = 1.1,
          common.legend = TRUE,
          legend = "top")


all_change_scatter

             
```


## Statistical Analysis

### PHQ - Primary Measure for Depression
Change in PHQ-9 scores will be entered as the outcome and change in bias based on errors to criterion in the self and other conditions will be entered as exposures in wide format. Time-point will be entered as a fixed effect, and participant ID as a random effect, to allow for clustering over time-points. Depression severity will be adjusted for by entering baseline PHQ-9 scores as a fixed effect. 

```{r}
hyp2_phq <- lmer(PHQ_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + PHQ_tot_baseline + change_timepoint + (1|subject), data = bias_change)
summ(hyp2_phq, confint = TRUE, digits = 3)
hyp2_phq_model_effects <- ggpredict(hyp2_phq, terms = c("Self_bias_change"))
plot(hyp2_phq_model_effects)

#Getting overall p-value for main effect of session
hyp2_phq_no_session <- lmer(PHQ_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + PHQ_tot_baseline + (1|subject), data = bias_change)
anova(hyp2_phq, hyp2_phq_no_session)


#Regression Graphs
SEL_bias_PHQ_model_effects <- ggpredict(hyp2_phq)

Self_effects <- as.data.frame(SEL_bias_PHQ_model_effects$Self_bias_change)
Friend_effects <- as.data.frame(SEL_bias_PHQ_model_effects$Friend_bias_change)
Stranger_effects <- as.data.frame(SEL_bias_PHQ_model_effects$Stranger_bias_change)

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_bias_change", 1, ifelse(group == "Friend_bias_change", 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("Change in SEL Bias Scores (Positive - Negative)") + ylab("Change in PHQ-9 Scores") + theme_minimal()

```
### BDI - Secondary Measure of Depression
```{r}
hyp2_BDI <- lmer(BDI_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + BDI_tot_baseline + change_timepoint + (1|subject), data = bias_change)
summ(hyp2_BDI, confint = TRUE, digits = 3)

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

#Getting overall p-value for main effect of session
hyp2_BDI_no_session <- lmer(BDI_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + BDI_tot_baseline + (1|subject), data = bias_change)
anova(hyp2_BDI, hyp2_BDI_no_session)

#Regression Graphs
SEL_bias_BDI_model_effects <- ggpredict(hyp2_BDI)

Self_effects <- as.data.frame(SEL_bias_BDI_model_effects$Self_bias_change)
Friend_effects <- as.data.frame(SEL_bias_BDI_model_effects$Friend_bias_change)
Stranger_effects <- as.data.frame(SEL_bias_BDI_model_effects$Stranger_bias_change)

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_bias_change", 1, ifelse(group == "Friend_bias_change", 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("Change in SEL Bias Scores (Positive - Negative)") + ylab("Change in BDI-9 Scores") + theme_minimal()


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("Change in SEL Bias Scores (Positive - Negative)") + ylab("Change in BDI-9 Scores") + theme_minimal()
```


### Exploratory Analysis with GAD-7 Scores

```{r}
#Exploratory analysis with anxiety
hyp2_GAD <- lmer(GAD_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + GAD_tot_baseline + change_timepoint + (1|subject), data = bias_change)
summ(hyp2_GAD, confint = TRUE, digits = 3)
standardize_parameters(hyp2_GAD, method = "refit", ci = 0.95, include_response = TRUE)

#Getting overall p-value for main effect of session
hyp2_GAD_no_session <- lmer(GAD_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + GAD_tot_baseline + (1|subject), data = bias_change)
anova(hyp2_GAD, hyp2_GAD_no_session)


#Regression Graphs
SEL_bias_GAD_model_effects <- ggpredict(hyp2_GAD)

Self_effects <- as.data.frame(SEL_bias_GAD_model_effects$Self_bias_change)
Friend_effects <- as.data.frame(SEL_bias_GAD_model_effects$Friend_bias_change)
Stranger_effects <- as.data.frame(SEL_bias_GAD_model_effects$Stranger_bias_change)

SEL_bias_GAD_model_effects_graph <- bind_rows(Self_effects, Friend_effects, Stranger_effects)

SEL_bias_GAD_model_effects_graph <- SEL_bias_GAD_model_effects_graph %>%
  mutate(Condition = ifelse(group == "Self_bias_change", 1, ifelse(group == "Friend_bias_change", 2, 3)))

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


ggplot(SEL_bias_GAD_model_effects_graph, aes(x = x, y = predicted, color = Condition)) + geom_line(size = 1) + geom_ribbon(data = SEL_bias_GAD_model_effects_graph, aes(ymin = conf.low, ymax = conf.high), alpha = .1, linetype=0) + xlab("Change in SEL Bias Scores (Positive - Negative)") + ylab("Change in GAD-7 Scores") + theme_minimal(base_size = 10)


```

#### Adjusting for depression

```{r}
#Adjusting for depression

##PHQ
hyp2_GAD_adj_PHQ <- lmer(GAD_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + GAD_tot_baseline +  change_timepoint + PHQ_change + (1|subject), data = bias_change)
summ(hyp2_GAD_adj_PHQ, confint = TRUE, digits = 3)
standardize_parameters(hyp2_GAD_adj_PHQ, method = "refit", ci = 0.95, include_response = TRUE)
#Getting overall p-value for main effect of session
hyp2_GAD_adj_PHQ_no_session <- lmer(GAD_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + GAD_tot_baseline +  PHQ_change + (1|subject), data = bias_change)
anova(hyp2_GAD_adj_PHQ, hyp2_GAD_adj_PHQ_no_session)

##BDI
hyp2_GAD_adj_BDI <- lmer(GAD_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + GAD_tot_baseline + change_timepoint + BDI_change + (1|subject), data = bias_change)
summ(hyp2_GAD_adj_BDI, confint = TRUE, digits = 3)
standardize_parameters(hyp2_GAD_adj_BDI, method = "refit", ci = 0.95, include_response = TRUE)
#Getting overall p-value for main effect of session
hyp2_GAD_adj_BDI_no_session <- lmer(GAD_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + GAD_tot_baseline +  BDI_change + (1|subject), data = bias_change)
anova(hyp2_GAD_adj_BDI, hyp2_GAD_adj_BDI_no_session)
```



# 6 Month Follow-Up

## Descriptives

```{r}

descriptives_table_LTFU <- SEL_bias_Qs %>%
  select(subject, session, Condition, bias, PHQ_tot, BDI_tot, GAD_tot, GRC_binary_factor) %>%
  pivot_wider(values_from = bias, names_from = Condition) %>% mutate(session_FU = as.numeric(session)) %>%
  mutate(session_FU = factor(session, levels = c(1, 2, 3, 4, 5), labels = c("Baseline", "2 Weeks", "6 Weeks", "8 Weeks", "6 Months")))

descriptives_table_LTFU %>%
  tab_cols(session_FU) %>%
  tab_cells(PHQ_tot, BDI_tot, GAD_tot, Self, Friend, Stranger) %>%
  tab_stat_mean() %>%
  tab_stat_sd() %>%
  tab_stat_valid_n() %>%
  tab_pivot(stat_position = "inside_columns")%>%
  set_caption("Mean self-report measures of mood by timepoint")

descriptives_table_LTFU %>%
  tab_cols(session_FU) %>%
  tab_cells(GRC_binary_factor) %>%
  tab_stat_cases() %>%
  tab_stat_cpct() %>%
    tab_pivot(stat_position = "inside_columns")%>%
  set_caption("GRC Responses by timepoint")

### PHQ Scores
PHQ_total_LTFU_graph <- SEL_bias_Qs %>%
  select(subject, session, PHQ_tot) %>%
  mutate(session_FU = as.numeric(session)) %>%
  mutate(session_FU = factor(session, levels = c(1, 2, 3, 4, 5), labels = c("Baseline", "2 Weeks", "6 Weeks", "8 Weeks", "6 Months"))) %>%
  group_by(session_FU) %>%
  summarise( 
    n=n(),
    mean=mean(PHQ_tot, na.rm = TRUE),
    sd=sd(PHQ_tot, na.rm = TRUE)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
  ungroup()
PHQ_total_LTFU_graph %>%
  ggplot(., aes(x = session_FU, y = mean, group = 1)) + 
  geom_point(stat = 'identity', size = 2) +
  geom_line(stat = 'identity', linetype = 2) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2, position=position_dodge(.2), size = 0.8) +
  ylab("PHQ Scores") + 
  xlab ("Timepoint") + 
  theme_minimal()




### BDI Scores
BDI_total_LTFU_graph <- SEL_bias_Qs %>%
  select(subject, session, BDI_tot) %>%
  mutate(session_FU = as.numeric(session)) %>%
  mutate(session_FU = factor(session, levels = c(1, 2, 3, 4, 5), labels = c("Baseline", "2 Weeks", "6 Weeks", "8 Weeks", "6 Months"))) %>%
  group_by(session_FU) %>%
  summarise( 
    n=n(),
    mean=mean(BDI_tot, na.rm = TRUE),
    sd=sd(BDI_tot, na.rm = TRUE)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
  ungroup()
BDI_total_LTFU_graph %>%
  ggplot(., aes(x = session_FU, y = mean, group = 1)) + 
  geom_point(stat = 'identity', size = 2) +
  geom_line(stat = 'identity', linetype = 2) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2, position=position_dodge(.2), size = 0.8) +
  ylab("BDI Scores") + 
  xlab ("Timepoint") + 
  theme_minimal()



#GAD
GAD_total_LTFU_graph <- SEL_bias_Qs %>%
  select(subject, session, GAD_tot) %>%
  mutate(session_FU = as.numeric(session)) %>%
  mutate(session_FU = factor(session, levels = c(1, 2, 3, 4, 5), labels = c("Baseline", "2 Weeks", "6 Weeks", "8 Weeks", "6 Months"))) %>%
  group_by(session_FU) %>%
  summarise( 
    n=n(),
    mean=mean(GAD_tot, na.rm = TRUE),
    sd=sd(GAD_tot, na.rm = TRUE)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
  ungroup()
GAD_total_LTFU_graph %>%
  ggplot(., aes(x = session_FU, y = mean, group = 1)) + 
  geom_point(stat = 'identity', size = 2) +
  geom_line(stat = 'identity', linetype = 2) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2, position=position_dodge(.2), size = 0.8) +
  ylab("GAD Scores") + 
  xlab ("Timepoint") + 
  theme_minimal()

##SEL Bias Scores by Condition

SEL_bias_LTFU_graph <- SEL_bias_Qs %>%
  select(subject, session, Condition, bias) %>%
  mutate(session_FU = as.numeric(session)) %>%
  mutate(session_FU = factor(session, levels = c(1, 2, 3, 4, 5), labels = c("Baseline", "2 Weeks", "6 Weeks", "8 Weeks", "6 Months"))) %>%
  group_by(session_FU, Condition) %>%
  summarise( 
    n=n(),
    mean=mean(bias, na.rm = TRUE),
    sd=sd(bias, na.rm = TRUE)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
  ungroup()

SEL_bias_LTFU_graph %>%
  ggplot(., aes(x = session_FU, y = mean, colour = Condition, group = Condition)) + 
  geom_point(stat = 'identity', position = position_dodge(.2), size = 2) +
  geom_line(stat = 'identity', position = position_dodge(.2), linetype = 2) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2, position=position_dodge(.2), size = 0.8) +
  ylab("Bias Scores") + 
  xlab ("Timepoint") + 
  labs(colour = "Condition") +
  theme_minimal(base_size = 10)


## Self Bias Scores
Self_SEL_bias_LTFU_graph <- SEL_bias_LTFU_graph %>%
  filter(Condition == "Self")
Self_SEL_bias_LTFU_graph %>%
  ggplot(., aes(x = session_FU, y = mean, group = 1)) + 
  geom_point(stat = 'identity', size = 2) +
  geom_line(stat = 'identity', linetype = 2) +
 geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2, position=position_dodge(.2), size = 0.8) +
  ylab("Self Bias Scores") + 
  xlab ("Timepoint") + 
  theme_minimal()


```

# Outlier Senstivity Analysis
```{r}
#Replacing outliers with NA
bias_change_outlier <- bias_change %>%
  mutate(Self_bias_zscore = (Self_bias_change - mean(Self_bias_change, na.rm = TRUE))/sd(Self_bias_change, na.rm = TRUE),
         Friend_bias_zscore = (Friend_bias_change - mean(Friend_bias_change, na.rm = TRUE))/sd(Friend_bias_change, na.rm = TRUE),
         Stranger_bias_zscore = (Stranger_bias_change - mean(Stranger_bias_change, na.rm = TRUE))/sd(Stranger_bias_change, na.rm = TRUE)) %>% #Calculating Z scores
  mutate(Self_bias_change = ifelse(Self_bias_zscore > 3 | Self_bias_zscore < -3, NA, Self_bias_change),
         Friend_bias_change = ifelse(Friend_bias_zscore > 3 | Friend_bias_zscore < -3, NA, Friend_bias_change),
         Stranger_bias_change = ifelse(Stranger_bias_zscore > 3 | Stranger_bias_zscore < -3, NA, Stranger_bias_change)) %>% #Removing values 3 SDs from mean
  mutate(Self_outlier = ifelse(Self_bias_zscore > 3 | Self_bias_zscore < -3, TRUE, FALSE),
         Friend_outlier = ifelse(Friend_bias_zscore > 3 | Friend_bias_zscore < -3, TRUE, FALSE),
         Stranger_outlier = ifelse(Stranger_bias_zscore > 3 | Stranger_bias_zscore < -3, TRUE, FALSE)) #Identifying outliers to report number


##Number of outliers
bias_change_outlier %>%
  group_by(change_timepoint) %>%
  summarise(Self_outlier = sum(Self_outlier, na.rm = TRUE), Friend_outlier = sum(Friend_outlier, na.rm = TRUE), Stranger_outlier = sum(Stranger_outlier, na.rm = TRUE))

```
## Visualising Data
```{r}
#Prepping dataframe
bias_change_outlier_wide <- bias_change_outlier %>%
  pivot_longer(cols = c(Self_bias_change, Friend_bias_change, Stranger_bias_change), names_to = "condition", values_to = "bias_change") %>%
  mutate(Condition = ifelse(condition == "Self_bias_change", 1,
                            ifelse(condition == "Friend_bias_change", 2, 3))) %>%
           mutate(Condition = factor(Condition, levels = c(1, 2, 3), labels = c("Self", "Friend", "Stranger"))) %>%
  select(-condition)

# PHQ
ggplot(bias_change_outlier_wide, aes(x = bias_change, y = PHQ_change, colour = Condition)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(. ~ change_timepoint + . ~ Condition) +
  theme_minimal()+
  ggtitle("Change in Bias scores and PHQ-9")


# BDI
ggplot(bias_change_outlier_wide, aes(x = bias_change, y = BDI_change, colour = Condition)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(. ~ change_timepoint + . ~ Condition) +
  theme_minimal()+
  ggtitle("Change in Bias scores and BDI-II")


#GAD
ggplot(bias_change_outlier_wide, aes(x = bias_change, y = GAD_change, colour = Condition)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(. ~ change_timepoint + . ~ Condition) +
  theme_minimal()+
  ggtitle("Change in Bias scores and GAD-7")

```



```{r}
#Prepping dataframe
bias_change_outlier_wide <- bias_change_outlier %>%
  pivot_longer(cols = c(Self_bias_change, Friend_bias_change, Stranger_bias_change), names_to = "condition", values_to = "bias_change") %>%
  mutate(Condition = ifelse(condition == "Self_bias_change", 1,
                            ifelse(condition == "Friend_bias_change", 2, 3))) %>%
           mutate(Condition = factor(Condition, levels = c(1, 2, 3), labels = c("Self", "Friend", "Stranger"))) %>%
  select(-condition)

# PHQ
PHQ_change_scatter_outliers <- ggplot(bias_change_outlier_wide, aes(x = bias_change, y = PHQ_change, colour = Condition)) +
  geom_point() +
  geom_smooth(method = "lm") +
  ylim(-17, 10) +
  facet_grid(. ~ Condition + . ~ change_timepoint) +
  xlab("Change in Bias Scores (Positive - Negative)") +
  ylab("Change in PHQ-9") +
  theme(strip.text.y = element_text(size = 8, colour = "black", angle = 270),
        strip.background = element_rect(colour = "black", fill = "white"),
        panel.border = element_rect(colour = "black", fill=NA, size=1))

# BDI
BDI_change_scatter_outliers <- ggplot(bias_change_outlier_wide, aes(x = bias_change, y = BDI_change, colour = Condition)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(. ~ Condition + . ~ change_timepoint) +
    xlab("Change in Bias Scores (Positive - Negative)") +
  ylab("Change in BDI-II") +
  theme(strip.text.y = element_text(size = 8, colour = "black", angle = 270),
        strip.background = element_rect(colour = "black", fill = "white"),
        panel.border = element_rect(colour = "black", fill=NA, size=1))


#GAD
GAD_change_scatter_outliers <- ggplot(bias_change_outlier_wide, aes(x = bias_change, y = GAD_change, colour = Condition)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(. ~ Condition + . ~ change_timepoint) +
  xlab("Change in Bias Scores (Positive - Negative)") +
  ylab("Change in GAD-7") +
  theme(strip.text.y = element_text(size = 8, colour = "black", angle = 270),
        strip.background = element_rect(colour = "black", fill = "white"),
        panel.border = element_rect(colour = "black", fill=NA, size=1))




#Bringing together plots
all_change_scatter_outliers <- ggarrange(PHQ_change_scatter_outliers, NULL, BDI_change_scatter_outliers, NULL, GAD_change_scatter_outliers, 
          labels = c("PHQ-9", "", "BDI-II", "", "GAD-7"), 
          font.label = list(size = 12, face = "bold"),
          ncol = 1, nrow = 5, heights = c(1, 0.1, 1, 0.1, 1),
          label.y = 1.1,
          common.legend = TRUE,
          legend = "top")


all_change_scatter_outliers

```





## Repeating Statistical Analysis For Hypothesis Two

### PHQ
```{r}
hyp2_phq_outlier <- lmer(PHQ_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + PHQ_tot_baseline + change_timepoint + (1|subject), data = bias_change_outlier)
summ(hyp2_phq_outlier, confint = TRUE, digits = 3)
standardize_parameters(hyp2_phq_outlier, method = "refit", ci = 0.95, include_response = TRUE)

#Getting overall p-value for main effect of session
hyp2_phq_outlier_no_session <- lmer(PHQ_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + PHQ_tot_baseline + (1|subject), data = bias_change_outlier)
anova(hyp2_phq_outlier, hyp2_phq_outlier_no_session)
```

### BDI
```{r}
hyp2_BDI_outlier <- lmer(BDI_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + BDI_tot_baseline + change_timepoint + (1|subject), data = bias_change_outlier)
summ(hyp2_BDI_outlier, confint = TRUE, digits = 3)
standardize_parameters(hyp2_BDI_outlier, method = "refit", ci = 0.95, include_response = TRUE)

#Getting overall p-value for main effect of session
hyp2_BDI_outlier_no_session <- lmer(BDI_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + BDI_tot_baseline + (1|subject), data = bias_change_outlier)
anova(hyp2_BDI_outlier, hyp2_BDI_outlier_no_session)
```

### GAD
```{r}

hyp2_GAD_outlier <- lmer(GAD_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + GAD_tot_baseline + change_timepoint + (1|subject), data = bias_change_outlier)
summ(hyp2_GAD_outlier, confint = TRUE, digits = 3)
standardize_parameters(hyp2_GAD_outlier, method = "refit", ci = 0.95, include_response = TRUE)

#Getting overall p-value for main effect of session
hyp2_GAD_outlier_no_session <- lmer(GAD_change ~ Self_bias_change + Friend_bias_change + Stranger_bias_change + GAD_tot_baseline + (1|subject), data = bias_change_outlier)
anova(hyp2_GAD_outlier, hyp2_GAD_outlier_no_session)
```