---
title: The effect of acute citalopram on self-referential emotional processing and
  social cognition in healthy volunteers
author: "Catherine Hobbs (c.hobbs@bath.ac.uk)"
output: html_document
---

```{r eval=FALSE, include=FALSE}
library(tidyverse)
library(expss)
library(tadaatoolbox)
library(lmerTest)
library(jtools)
library(ggeffects)
library(data.table)
library(emmeans)
library(lmtest)
```

# Loading clean dataframes
```{r}
setwd("/Volumes/files/RDSSD/Katie Studies/Citalopram/8. Data/Analysis/")
load(file = "Questionnaires/demo_qs")
load(file = "Associative Learning/associative_long")
load(file = "GNAT/GNAT_agg_long_excluded")
load(file = "GNAT/GNAT_agg_long")
load(file = "Recall/cat_recall_long")
load(file = "SEL/SEL_clean_aggregate")
load(file = "SEL/SEL_bias")
load(file = "Prisoner/prisoner_clean_aggregate")
```
# Unblinding Data

This data was initially analysed blind with drug group categorised as "A" or "B". To aid comprehension when writing up drug group was unblinded.

```{r}
#Unblinding
demo_qs$group_unblind <- ifelse(demo_qs$group == "A", 1, 2)

demo_qs$group_unblind <- as.factor(demo_qs$group_unblind)

demo_qs$group_unblind <- factor(demo_qs$group_unblind,
                  levels = c(1,2),
                  labels = c("Citalopram", "Placebo"))

associative_long$group_unblind <- ifelse(associative_long$group == "A", 1, 2)

associative_long$group_unblind <- as.factor(associative_long$group_unblind)

associative_long$group_unblind <- factor(associative_long$group_unblind,
                  levels = c(1,2),
                  labels = c("Citalopram", "Placebo"))

GNAT_agg_long$group_unblind <- ifelse(GNAT_agg_long$group == "A", 1, 2)

GNAT_agg_long$group_unblind <- as.factor(GNAT_agg_long$group_unblind)

GNAT_agg_long$group_unblind <- factor(GNAT_agg_long$group_unblind,
                  levels = c(1,2),
                  labels = c("Citalopram", "Placebo"))


GNAT_agg_long_excluded$group_unblind <- ifelse(GNAT_agg_long_excluded$group == "A", 1, 2)

GNAT_agg_long_excluded$group_unblind <- as.factor(GNAT_agg_long_excluded$group_unblind)

GNAT_agg_long_excluded$group_unblind <- factor(GNAT_agg_long_excluded$group_unblind,
                  levels = c(1,2),
                  labels = c("Citalopram", "Placebo"))


cat_recall_long$group_unblind <- ifelse(cat_recall_long$group == "A", 1, 2)

cat_recall_long$group_unblind <- as.factor(cat_recall_long$group_unblind)

cat_recall_long$group_unblind <- factor(cat_recall_long$group_unblind,
                  levels = c(1,2),
                  labels = c("Citalopram", "Placebo"))

SEL_aggregate$group_unblind <- ifelse(SEL_aggregate$group == "A", 1, 2)

SEL_aggregate$group_unblind <- as.factor(SEL_aggregate$group_unblind)

SEL_aggregate$group_unblind <- factor(SEL_aggregate$group_unblind,
                  levels = c(1,2),
                  labels = c("Citalopram", "Placebo"))

SEL_bias$group_unblind <- ifelse(SEL_bias$group == "A", 1, 2)

SEL_bias$group_unblind <- as.factor(SEL_bias$group_unblind)

SEL_bias$group_unblind <- factor(SEL_bias$group_unblind,
                  levels = c(1,2),
                  labels = c("Citalopram", "Placebo"))

prisoner_clean_aggregate$group_unblind <- ifelse(prisoner_clean_aggregate$group == "A", 1, 2)

prisoner_clean_aggregate$group_unblind <- as.factor(prisoner_clean_aggregate$group_unblind)

prisoner_clean_aggregate$group_unblind <- factor(prisoner_clean_aggregate$group_unblind,
                  levels = c(1,2),
                  labels = c("Citalopram", "Placebo"))

```

# Adding labels for tables

Note, a dataframe was kept unlabelled for each task for the regression analyses.

```{r}
#Adding Labels for Tables

##Questionnaires & Demographics
demo_qs_no_labels <- demo_qs
demo_qs = apply_labels(demo_qs,
                       group = "Drug Group",
                       gender = "Gender", 
                       age_years = "Age (years)",
                       employment = "Employment",
                       smoke = "Smoker",
                       smoke_amount = "Amount Smoked per day",
                       alcohol_week = "Alcoholic units per week",
                       alcohol_max = "Maximum amount of Alcohol consumed in one session",
                       caffeine_drinks = "Number of caffeinated drinks per day",
                       ethnicity = "Ethnicity",
                       education = "Highest Educational Attainment",
                       english_first_language = "English Spoken as First Language",
                       PHQ_tot = "PHQ-9",
                       BDI_tot = "BDI-II",
                       BFNE_tot = "BFNE",
                       GAD_tot = "GAD-7",
                       STAIT_tot = "Trait Anxiety",
                       STAIS_tot = "State Anxiety Baseline",
                       PANAS_pos = "Positive Affect Baseline",
                       PANAS_neg = "Negative Affect Baseline",
                       VAS_1 = "VAS Happy Baseline",
                       VAS_2 = "VAS Sad Baseline",
                       VAS_3 = "VAS Disgust Baseline",
                       VAS_4 = "VAS Anger Baseline",
                       VAS_5 = "VAS Fear Baseline",
                       VAS_6 = "VAS Anxiety Baseline",
                       VAS_7 = "VAS Alert Baseline",
                       SideEffects_Q1 = "Nausea Baseline",
                       SideEffects_Q2 = "Dizziness Baseline",
                       SideEffects_Q3 = "Dry Mouth Baseline",
                       SideEffects_Q4 = "Headache Baseline",
                       SideEffects_Q5 = "Alert Baseline",
                       SideEffects_Q6 = "Agitation Baseline",
                       STAIS_tot_pd = "State Anxiety Post Drug",
                       PANAS_pos_pd = "Positive Affect Post Drug",
                       PANAS_neg_pd = "Negative Affect Post Drug",
                       VAS_pd_1 = "VAS Happy Post Drug",
                       VAS_pd_2 = "VAS Sad Post Drug",
                       VAS_pd_3 = "VAS Disgust Post Drug",
                       VAS_pd_4 = "VAS Anger Post Drug",
                       VAS_pd_5 = "VAS Fear Post Drug",
                       VAS_pd_6 = "VAS Anxiety Post Drug",
                       VAS_pd_7 = "VAS Alert Post Drug",
                       SideEffects_Q1_pd = "Nausea Post Drug",
                       SideEffects_Q2_pd = "Dizziness Post Drug",
                       SideEffects_Q3_pd = "Dry Mouth Post Drug",
                       SideEffects_Q4_pd = "Headache Post Drug",
                       SideEffects_Q5_pd = "Alert Post Drug",
                       SideEffects_Q6_pd = "Agitation Post Drug",
                       STAIS_tot_pc = "State Anxiety Post Cognitive Testing",
                       PANAS_pos_pc = "Positive Affect Post Cognitive Testing",
                       PANAS_neg_pc = "Negative Affect Post Cognitive Testing",
                       VAS_pc_1 = "VAS Happy Post Cognitive Testing",
                       VAS_pc_2 = "VAS Sad Post Cognitive Testing",
                       VAS_pc_3 = "VAS Disgust Post Cognitive Testing",
                       VAS_pc_4 = "VAS Anger Post Cognitive Testing",
                       VAS_pc_5 = "VAS Fear Post Cognitive Testing",
                       VAS_pc_6 = "VAS Anxiety Post Cognitive Testing",
                       VAS_pc_7 = "VAS Alert Post Cognitive Testing",
                       SideEffects_Q1_pc = "Nausea Post Cognitive Testing",
                       SideEffects_Q2_pc = "Dizziness Post Cognitive Testing",
                       SideEffects_Q3_pc = "Dry Mouth Post Cognitive Testing",
                       SideEffects_Q4_pc = "Headache Post Cognitive Testing",
                       SideEffects_Q5_pc = "Alert Post Cognitive Testing",
                       SideEffects_Q6_pc = "Agitation Post Cognitive Testing",
                       STAIS_base_pd = "State Anxiety Change between Baseline and Post-Drug",
                       PANAS_pos_base_pd = "Positive Affect Change between Baseline and Post-Drug",
                       PANAS_neg_base_pd = "Negative Affect Change between Baseline and Post-Drug",
                       VAS_1_base_pd = "VAS Happy Change between Baseline and Post-Drug",
                       VAS_2_base_pd = "VAS Sad Change between Baseline and Post-Drug",
                       VAS_3_base_pd = "VAS Disgust Change between Baseline and Post-Drug",
                       VAS_4_base_pd = "VAS Anger Change between Baseline and Post-Drug",
                       VAS_5_base_pd = "VAS Fear Change between Baseline and Post-Drug",
                       VAS_6_base_pd = "VAS Anxiety Change between Baseline and Post-Drug",
                       VAS_7_base_pd = "VAS Alert Change between Baseline and Post-Drug",
                       SideEffects_Q1_base_pd = "Nausea Change between Baseline and Post-Drug",
                       SideEffects_Q2_base_pd = "Dizziness Change between Baseline and Post-Drug",
                       SideEffects_Q3_base_pd = "Dry Mouth Change between Baseline and Post-Drug",
                       SideEffects_Q4_base_pd = "Headache Change between Baseline and Post-Drug",
                       SideEffects_Q5_base_pd = "Alert Change between Baseline and Post-Drug",
                       SideEffects_Q6_base_pd = "Agitation Post Change between Baseline and Post-Drug",
                       STAIS_pd_pc = "State Anxiety Change between Post-Drug and Post-Cognitive Testing",
                       PANAS_pos_pd_pc = "Positive Affect Change between Post-Drug and Post-Cognitive Testing",
                       PANAS_neg_pd_pc = "Negative Affect Change between Post-Drug and Post-Cognitive Testing",
                       VAS_1_pd_pc = "VAS Happy Change between Post-Drug and Post-Cognitive Testing",
                       VAS_2_pd_pc = "VAS Sad Change between Post-Drug and Post-Cognitive Testing",
                       VAS_3_pd_pc = "VAS Disgust Change between Post-Drug and Post-Cognitive Testing",
                       VAS_4_pd_pc = "VAS Anger Change between Post-Drug and Post-Cognitive Testing",
                       VAS_5_pd_pc = "VAS Fear Change between Post-Drug and Post-Cognitive Testing",
                       VAS_6_pd_pc = "VAS Anxiety Change between Post-Drug and Post-Cognitive Testing",
                       VAS_7_pd_pc = "VAS Alert Change between Post-Drug and Post-Cognitive Testing",
                       SideEffects_Q1_pd_pc = "Nausea Change between Post-Drug and Post-Cognitive Testing",
                       SideEffects_Q2_pd_pc = "Dizziness Change between Post-Drug and Post-Cognitive Testing",
                       SideEffects_Q3_pd_pc = "Dry Mouth Change between Post-Drug and Post-Cognitive Testing",
                       SideEffects_Q4_pd_pc = "Headache Change between Post-Drug and Post-Cognitive Testing",
                       SideEffects_Q5_pd_pc = "Alert Change between Post-Drug and Post-Cognitive Testing",
                       SideEffects_Q6_pd_pc = "Agitation Post Change between Post-Drug and Post-Cognitive Testing",
                       group_guess_participant = "Participant Group Guess",
                       group_guess_researcher = "Researcher Group Guess",
                       placebo_guess_certainty_participant = "Participant Certainty of being in Placebo Group",
                       drug_guess_certainty_participant = "Participant Certainty of being in Citalopram Group",
                       placebo_guess_certainty_researcher = "Researcher Certainty of being in Placebo Group",
                       drug_guess_certainty_researcher  = "Researcher Certainty of being in Citalopram Group",
                       EPQR_Extraversion = "EPQR Extraversion",
                       EPQR_Neuroticism = "EPQR Neuroticism",
                       EPQR_Lie = "EPQR Lie",
                       EPQR_Psychoticism = "EPQR Psychoticism",
                       group_unblind = "Drug Group")


associative_long_no_labels <- associative_long
associative_long = apply_labels(associative_long,
                       group = "Drug Group (blinded)",
                       group_unblind = "Drug Group",
                       prop_acc = "Accuracy (%)",
                       mean_rt_mult = "Reaction Times (ms)",
                       mean_rt_unmult = "Reaction Times (seconds)",
                       ShapeNo = "Stimuli Type")

GNAT_agg_long_no_labels <- GNAT_agg_long
GNAT_agg_long = apply_labels(GNAT_agg_long,
                             group = "Drug Group (blind)",
                             group_unblind = "Drug Group",
                             latency = "Reaction Times (ms)",
                             correct_prop = "Accuracy (%)",
                             hit_prop = "Hits (%)",
                             miss_prop = "Misses (%)",
                             correct_reject_prop = "Correct Rejections (%)",
                             fa_prop = "False Alarms (%)",
                             correct_sum = "Accuracy (total no.)",
                             hit_sum = "Hits (total no.)",
                             miss_sum = "Misses (total no.)",
                             fa_sum = "False alarms ((total no.)",
                             correct_reject_sum = "Correct Rejections (total no.)",
                             z_Hits = "Hits (z scores)",
                             z_Misses = "Misses (z scores)",
                             z_Correct_Reject = "Correct Rejections (z scores)",
                             z_FA = "False Alarms (z scores)",
                             Other_bias = "Other Bias Score",
                             Self_bias = "Self Biase Score",
                             Pos_Self_other_bias = "Positive Self-Other Bias",
                             Neg_Self_Other_bias = "Negative Self-Other Bias",
                             Overall_bias = "Overall Self-Esteem Bias",
                             ref_condition = "Condition",
                             valence = "Emotion")

GNAT_agg_long_excluded_no_labels <- GNAT_agg_long_excluded
GNAT_agg_long_excluded = apply_labels(GNAT_agg_long_excluded,
                             group = "Drug Group (blind)",
                             group_unblind = "Drug Group",
                             latency = "Reaction Times (ms)",
                             correct_prop = "Accuracy (%)",
                             hit_prop = "Hits (%)",
                             miss_prop = "Misses (%)",
                             correct_reject_prop = "Correct Rejections (%)",
                             fa_prop = "False Alarms (%)",
                             correct_sum = "Accuracy (total no.)",
                             hit_sum = "Hits (total no.)",
                             miss_sum = "Misses (total no.)",
                             fa_sum = "False alarms ((total no.)",
                             correct_reject_sum = "Correct Rejections (total no.)",
                             z_Hits = "Hits (z scores)",
                             z_Misses = "Misses (z scores)",
                             z_Correct_Reject = "Correct Rejections (z scores)",
                             z_FA = "False Alarms (z scores)",
                             Other_bias = "Other Bias Score",
                             Self_bias = "Self Biase Score",
                             Pos_Self_other_bias = "Positive Self-Other Bias",
                             Neg_Self_Other_bias = "Negative Self-Other Bias",
                             Overall_bias = "Overall Self-Esteem Bias",
                             ref_condition = "Condition",
                             valence = "Emotion")    


cat_recall_long_no_labels <- cat_recall_long
cat_recall_long = apply_labels(cat_recall_long,
                       group = "Drug Group (blind)",
                       group_unblind = "Drug Group",
                       condition = "Referential Condition",
                       stimuli_assignment = "Word List Used",
                       overall_number_recalled = "Total Words Correctly Recalled (irrespective of valence)",
                       number_recall_words_entered = "Number of Words submitted in Recall",
                       recall_tot_incorrect_other_condition = "Number of Words Incorrectly Recalled due to being in other Referential Condition",
                       recall_tot_fa = "Number of Words Incorrectly Recalled",
                       valence = "Valence",
                       cat_tot_yes = "Number of words categorised as referential (yes)",
                       cat_tot_no = "Number of words not categorised as referential (no)",
                       cat_rt = "Categorisation Reaction Times",
                       recall_tot_correct = "Number of Words Correctly Recalled")

SEL_aggregate_no_labels <- SEL_aggregate
SEL_aggregate = apply_labels(SEL_aggregate,
                             group = "Drug Group (blind)",
                             group_unblind = "Drug Group",
                             Condition = "Referential Condition",
                             baseline = "Baseline Perception Rating",
                             finalrating = "Final Perception Rating",
                             globalrating = "Global Rating",
                             perception = "Perception Rating",
                             ACC = "Accuracy (%)",
                             RT = "Reaction Times (ms)",
                             ERR = "Errors (%)",
                             ACC_tot = "Accuracy (total)",
                             ERR_tot = "Errors (total)",
                             ERR_to_crit = "Errors to Criterion",
                             POS_resp = "Positive Response Rate",
                             Overall_rule = "Rule")

SEL_bias_no_labels <- SEL_bias
SEL_bias = apply_labels(SEL_bias,
                        group = "Drug Group (blind)",
                        group_unblind = "Drug Group",
                        bias = "Errors to Criterion Bias Scores",
                        Condition = "Referential Condition")

prisoner_clean_aggregate_no_labels <- prisoner_clean_aggregate
prisoner_clean_aggregate = apply_labels(prisoner_clean_aggregate,
                        group = "Drug Group (blind)",
                        group_unblind = "Drug Group",
                        initial_other_choice = "Condition",
                        other_points_total = "Total Other Points Won",
                        participant_points_total = "Total Participant Points Won",
                        resp.time = "Reaction Times",
                        cooperation_proportion = "Proportion Cooperative Behaviours",
                        cooperation_total = "Total Number Cooperative Behaviours",
                        prop_previous_both_cooperate_choice = "Proportion Cooperative Behaviours After Both Players Cooperate",
                        prop_previous_p_coop_o_defect_choice = "Proportion Cooperative Behaviours After the Participant Cooperates and the Other Defects",
                        prop_previous_p_defect_o_coop_choice = "Proportion Cooperative Behaviours After the Participant Defects and the Other Cooperates",
                        prop_previous_both_defect_choice = "Proportion Cooperative Behaviours After Both Players Defect",
                        resp.time_coop = "Reaction Times for CooperateChoices",
                        resp.time_defect = "Reaction Times for Defect Choices")
                       
```




# Sample Demographics

```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

demo_qs %>%
  tab_cols(group_unblind) %>%
  tab_cells(gender, employment, ethnicity, education, english_first_language) %>%
  tab_stat_cases(total_row_position = c("none"))%>%
  tab_cells(age_years) %>%
  tab_stat_mean() %>%
  tab_stat_sd() %>%
  tab_pivot() %>%
  set_caption("Demographics by Group")

```

# Baseline Trait Mood

```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

demo_qs %>%
  tab_cols(group_unblind) %>%
  tab_cells(PHQ_tot, BDI_tot, GAD_tot, BFNE_tot, STAIT_tot, EPQR_Extraversion, EPQR_Neuroticism, EPQR_Lie, EPQR_Psychoticism) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot() %>%
  set_caption("Baseline Trait Mood Measures by Group")

```

# Change in Self-Reported State Mood

## Descriptive Statistics
```{r}
#Creating a long dataset with state mood measures
state_mood <- demo_qs_no_labels %>%
  select(participant, group, group_unblind, STAIS_tot, PANAS_pos, PANAS_neg, VAS_1, VAS_2, VAS_3, VAS_4, VAS_5, VAS_6, VAS_7, SideEffects_Q1, SideEffects_Q2, SideEffects_Q3, SideEffects_Q4, SideEffects_Q5, SideEffects_Q6,STAIS_tot_pd, PANAS_pos_pd, PANAS_neg_pd, VAS_pd_1, VAS_pd_2, VAS_pd_3, VAS_pd_4, VAS_pd_5, VAS_pd_6, VAS_pd_7, SideEffects_Q1_pd, SideEffects_Q2_pd, SideEffects_Q3_pd, SideEffects_Q4_pd, SideEffects_Q5_pd, SideEffects_Q6_pd, STAIS_tot_pc, PANAS_pos_pc, PANAS_neg_pc, VAS_pc_1, VAS_pc_2, VAS_pc_3, VAS_pc_4, VAS_pc_5, VAS_pc_6, VAS_pc_7, SideEffects_Q1_pc, SideEffects_Q2_pc, SideEffects_Q3_pc, SideEffects_Q4_pc, SideEffects_Q5_pc, SideEffects_Q6_pc) %>%
  rename(STAIS_baseline = STAIS_tot, PANAS_pos_baseline = PANAS_pos, PANAS_neg_baseline = PANAS_neg, VAS_1_baseline = VAS_1, VAS_2_baseline = VAS_2, VAS_3_baseline = VAS_3, VAS_4_baseline = VAS_4, VAS_5_baseline = VAS_5, VAS_6_baseline = VAS_6, VAS_7_baseline = VAS_7, SideEffects_Q1_baseline = SideEffects_Q1, SideEffects_Q2_baseline = SideEffects_Q2, SideEffects_Q3_baseline = SideEffects_Q3, SideEffects_Q4_baseline = SideEffects_Q4, SideEffects_Q5_baseline = SideEffects_Q5, SideEffects_Q6_baseline = SideEffects_Q6, VAS_1_pd = VAS_pd_1, VAS_2_pd = VAS_pd_2, VAS_3_pd = VAS_pd_3, VAS_4_pd = VAS_pd_4, VAS_5_pd = VAS_pd_5, VAS_6_pd = VAS_pd_6, VAS_7_pd = VAS_pd_7, VAS_1_pc = VAS_pc_1, VAS_2_pc = VAS_pc_2, VAS_3_pc = VAS_pc_3, VAS_4_pc = VAS_pc_4, VAS_5_pc = VAS_pc_5, VAS_6_pc = VAS_pc_6, VAS_7_pc = VAS_pc_7)
state_mood_long <- melt(setDT(state_mood), measure = patterns('^STAIS', '^PANAS_pos', '^PANAS_neg', '^VAS_1', '^VAS_2', '^VAS_3', '^VAS_4', '^VAS_5', '^VAS_6', '^VAS_7', '^SideEffects_Q1', '^SideEffects_Q2', '^SideEffects_Q3', '^SideEffects_Q4', '^SideEffects_Q5', '^SideEffects_Q6'), value.name = c('STAIS', 'PANAS_pos', 'PANAS_neg', 'VAS_1', 'VAS_2', 'VAS_3', 'VAS_4', 'VAS_5', 'VAS_6', 'VAS_7', 'SideEffects_Q1', 'SideEffects_Q2', 'SideEffects_Q3', 'SideEffects_Q4', 'SideEffects_Q5', 'SideEffects_Q6'), variable.name = 'Time')

state_mood_long$Time <- as.factor(state_mood_long$Time)

state_mood_long$Time <- factor(state_mood_long$Time,
                  levels = c(1,2,3),
                  labels = c("Baseline", "Post-Drug", "Post-Testing"))

state_mood_long_no_labels <- state_mood_long
state_mood_long = apply_labels(state_mood_long,
                       group = "Drug Group (blinded)",
                       group_unblind = "Drug Group",
                       STAIS = "State Anxiety (STAI-S)",
                       PANAS_pos = "PANAS Positive",
                       PANAS_neg = "PANAS Negative",
                       VAS_1 = "VAS Happy",
                       VAS_2 = "VAS Sad",
                       VAS_3 = "VAS Disgust",
                       VAS_4 = "VAS Anger",
                       VAS_5 = "VAS Fear",
                       VAS_6 = "VAS Anxious",
                       VAS_7 = "VAS Alert",
                       SideEffects_Q1 = "Nausea",
                       SideEffects_Q2 = "Dizziness",
                       SideEffects_Q3 = "Dry Mouth",
                       SideEffects_Q4 = "Headache",
                       SideEffects_Q5 = "Alertness",
                       SideEffects_Q6 = "Agitation")
                       
#Overall descriptive statistics
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

state_mood_long %>%
  tab_cols(group_unblind %nest% Time) %>%
  tab_cells(STAIS, PANAS_pos, PANAS_neg, VAS_1, VAS_2, VAS_3, VAS_4, VAS_5, VAS_6, VAS_7) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot() %>%
  set_caption("State Mood Measures by Group and Time")
                       
```

## Regression Analyses examining change over time by drug group

### STAI-S
```{r}
STAIS_group_time <- lmer(STAIS ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(STAIS_group_time, confint = TRUE, digits = 3)
STAIS_group_time_effects <- ggpredict(STAIS_group_time, terms = c("Time", "group_unblind"))
plot(STAIS_group_time_effects) +
  labs(
    title = "Predicted Values of State Anxiety (STAIS) from Mixed Regression Model",
    y = "State Anxiety (STAIS)",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of time 
STAIS_no_time <- lmer(STAIS ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(STAIS_group_time, STAIS_no_time)

#Getting overall p-value for interaction effect of time
STAIS_main_time <- lmer(STAIS ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(STAIS_group_time, STAIS_main_time)
```

### PANAS Positive
```{r}
PANAS_pos_group_time <- lmer(PANAS_pos ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(PANAS_pos_group_time, confint = TRUE, digits = 3)
PANAS_pos_group_time_effects <- ggpredict(PANAS_pos_group_time, terms = c("Time", "group_unblind"))
plot(PANAS_pos_group_time_effects) +
  labs(
    title = "Predicted Values of PANAS Positive from Mixed Regression Model",
    y = "PANAS Positive",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of time 
PANAS_pos_no_time <- lmer(PANAS_pos ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(PANAS_pos_group_time, PANAS_pos_no_time)

#Getting overall p-value for interaction effect of time
PANAS_pos_main_time <- lmer(PANAS_pos ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(PANAS_pos_group_time, PANAS_pos_main_time)
```

### PANAS Negative
```{r}
PANAS_neg_group_time <- lmer(PANAS_neg ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(PANAS_neg_group_time, confint = TRUE, digits = 3)
PANAS_neg_group_time_effects <- ggpredict(PANAS_neg_group_time, terms = c("Time", "group_unblind"))
plot(PANAS_neg_group_time_effects) +
  labs(
    title = "Predicted Values of PANAS Negative from Mixed Regression Model",
    y = "PANAS Negative",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of time 
PANAS_neg_no_time <- lmer(PANAS_neg ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(PANAS_neg_group_time, PANAS_neg_no_time)

#Getting overall p-value for interaction effect of time
PANAS_neg_main_time <- lmer(PANAS_neg ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(PANAS_neg_group_time, PANAS_neg_main_time)
```

### VAS Happy
```{r}
VAS_1_group_time <- lmer(VAS_1 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(VAS_1_group_time, confint = TRUE, digits = 3)
VAS_1_group_time_effects <- ggpredict(VAS_1_group_time, terms = c("Time", "group_unblind"))
plot(VAS_1_group_time_effects) +
  labs(
    title = "Predicted Values of VAS Happy from Mixed Regression Model",
    y = "VAS Happy",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Contrasts for drug effect at post-drug and post-testing
contrast_VAS_1 <- emmeans(VAS_1_group_time, ~ group_unblind*Time)
contrast(contrast_VAS_1, "revpairwise", by = "Time", adjust = "bonferroni")

state_mood_long %>%
  group_by(group_unblind, Time) %>%
  summarize(Mean = mean(VAS_1, na.rm = TRUE),  SD = sd(VAS_1, na.rm = TRUE), Median = median(VAS_1, na.rm = TRUE), Minimum = min(VAS_1, na.rm = TRUE), Maximum = max(VAS_1, na.rm = TRUE)) %>%
  ggplot(., aes(x = Time, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("VAS Happy according to Drug Group over Time") + 
  ylab("VAS Happy") + 
  xlab("Timepoint") +
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

#Getting overall p-value for main effect of time 
VAS_1_no_time <- lmer(VAS_1 ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(VAS_1_group_time, VAS_1_no_time)

#Getting overall p-value for interaction effect of time
VAS_1_main_time <- lmer(VAS_1 ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(VAS_1_group_time, VAS_1_main_time)
```

### VAS Sad
```{r}
VAS_2_group_time <- lmer(VAS_2 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(VAS_2_group_time, confint = TRUE, digits = 3)
VAS_2_group_time_effects <- ggpredict(VAS_2_group_time, terms = c("Time", "group_unblind"))
plot(VAS_2_group_time_effects) +
  labs(
    title = "Predicted Values of VAS Sad from Mixed Regression Model",
    y = "VAS Sad",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Contrasts for drug effect at post-drug and post-testing
contrast_VAS_2 <- emmeans(VAS_2_group_time, ~ group_unblind*Time)
contrast(contrast_VAS_2, "revpairwise", by = "Time", adjust = "bonferroni")

state_mood_long %>%
  group_by(group_unblind, Time) %>%
  summarize(Mean = mean(VAS_2, na.rm = TRUE),  SD = sd(VAS_1, na.rm = TRUE), Median = median(VAS_1, na.rm = TRUE), Minimum = min(VAS_1, na.rm = TRUE), Maximum = max(VAS_1, na.rm = TRUE)) %>%
  ggplot(., aes(x = Time, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("VAS Sad according to Drug Group over Time") + 
  ylab("VAS Sad") + 
  xlab("Timepoint") +
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

#Getting overall p-value for main effect of time 
VAS_2_no_time <- lmer(VAS_2 ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(VAS_2_group_time, VAS_2_no_time)

#Getting overall p-value for interaction effect of time
VAS_2_main_time <- lmer(VAS_2 ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(VAS_2_group_time, VAS_2_main_time)

```
### VAS Disgust
```{r}
VAS_3_group_time <- lmer(VAS_3 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(VAS_3_group_time, confint = TRUE, digits = 3)
VAS_3_group_time_effects <- ggpredict(VAS_3_group_time, terms = c("Time", "group_unblind"))
plot(VAS_3_group_time_effects) +
  labs(
    title = "Predicted Values of VAS Disgust from Mixed Regression Model",
    y = "VAS Disgust",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of time 
VAS_3_no_time <- lmer(VAS_3 ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(VAS_3_group_time, VAS_3_no_time)

#Getting overall p-value for interaction effect of time
VAS_3_main_time <- lmer(VAS_3 ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(VAS_3_group_time, VAS_3_main_time)
```
### VAS Anger
```{r}
VAS_4_group_time <- lmer(VAS_4 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(VAS_4_group_time, confint = TRUE, digits = 3)
VAS_4_group_time_effects <- ggpredict(VAS_4_group_time, terms = c("Time", "group_unblind"))
plot(VAS_4_group_time_effects) +
  labs(
    title = "Predicted Values of VAS Anger from Mixed Regression Model",
    y = "VAS Anger",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of time 
VAS_4_no_time <- lmer(VAS_4 ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(VAS_4_group_time, VAS_4_no_time)

#Getting overall p-value for interaction effect of time
VAS_4_main_time <- lmer(VAS_4 ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(VAS_4_group_time, VAS_4_main_time)
```
### VAS Fear
```{r}
VAS_5_group_time <- lmer(VAS_5 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(VAS_5_group_time, confint = TRUE, digits = 3)
VAS_5_group_time_effects <- ggpredict(VAS_5_group_time, terms = c("Time", "group_unblind"))
plot(VAS_5_group_time_effects) +
  labs(
    title = "Predicted Values of VAS Fear from Mixed Regression Model",
    y = "VAS Fear",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of time 
VAS_5_no_time <- lmer(VAS_5 ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(VAS_5_group_time, VAS_5_no_time)

#Getting overall p-value for interaction effect of time
VAS_5_main_time <- lmer(VAS_5 ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(VAS_5_group_time, VAS_5_main_time)
```
### VAS Anxiety
```{r}
VAS_6_group_time <- lmer(VAS_6 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(VAS_6_group_time, confint = TRUE, digits = 3)
VAS_6_group_time_effects <- ggpredict(VAS_6_group_time, terms = c("Time", "group_unblind"))
plot(VAS_6_group_time_effects) +
  labs(
    title = "Predicted Values of VAS Anxiety from Mixed Regression Model",
    y = "VAS Anxiety",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of time 
VAS_6_no_time <- lmer(VAS_6 ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(VAS_6_group_time, VAS_6_no_time)

#Getting overall p-value for interaction effect of time
VAS_6_main_time <- lmer(VAS_6 ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(VAS_6_group_time, VAS_6_main_time)
```
### VAS Alertness
```{r}
VAS_7_group_time <- lmer(VAS_7 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(VAS_7_group_time, confint = TRUE, digits = 3)
VAS_7_group_time_effects <- ggpredict(VAS_7_group_time, terms = c("Time", "group_unblind"))
plot(VAS_7_group_time_effects) +
  labs(
    title = "Predicted Values of VAS Alertness from Mixed Regression Model",
    y = "VAS Alertness",
    x = "Timepoint",
    colour = "Drug Group"
  )
```


# Prisoner's Dilemma

## Proportion of Cooperative Behaviours
```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

prisoner_clean_aggregate %>%
  tab_cells(cooperation_proportion) %>%
  tab_cols(group_unblind %nest% initial_other_choice) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()

prisoner_coop_prop_graph <- prisoner_clean_aggregate %>%
  group_by(group_unblind, initial_other_choice) %>%
  summarize(Mean = mean(cooperation_proportion),  SD = sd(cooperation_proportion), Median = median(cooperation_proportion), Minimum = min(cooperation_proportion), Maximum = max(cooperation_proportion)) %>%
  ggplot(., aes(x = initial_other_choice, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Cooperative Behaviours (%)") + 
  xlab ("Social Context") + 
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))+
  scale_fill_brewer()

print(prisoner_coop_prop_graph)
ggsave("prisoner_coop_prop_graph.jpeg")

coop_prop_model <- lmer(cooperation_proportion ~ group_unblind*initial_other_choice + (1|participant), data = prisoner_clean_aggregate_no_labels)
summ(coop_prop_model, confint = TRUE, digits = 3)
coop_prop_model_effects <- ggpredict(coop_prop_model, terms = c("initial_other_choice", "group_unblind"))
plot(coop_prop_model_effects) + 
  labs(
    x = "Social Context", 
    y = "Cooperative Behaviours (%)", 
    title = "Predicted Values of Cooperative Behaviours (%)",
    colour = "Drug Group"
  )
```


# Social Evaluation Learning

## Bias Scores

```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

SEL_bias %>%
  tab_cells(bias) %>%
  tab_cols(group_unblind, group_unblind %nest% Condition) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()

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

print(SEL_bias_graph)
ggsave("SEL_bias_graph.jpeg")

bias_model <- lmer(bias ~ group_unblind*Condition + (1|participant), data = SEL_bias_no_labels)
summ(bias_model, confint = TRUE, digits = 3)
bias_model_effects <- ggpredict(bias_model, terms = c("Condition", "group_unblind"))
plot(bias_model_effects) + 
  labs(
    x = "Referential Condition", 
    y = "Bias Scores", 
    title = "Predicted Values of Bias Scores",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of referential condition
bias_model_no_referential_condition <- lmer(bias ~ group_unblind + (1|participant), data = SEL_bias_no_labels)
anova(bias_model, bias_model_no_referential_condition)

#Getting overall p-value for interaction effect of referential condition
bias_model_main_referential_condition <- lmer(bias ~ group_unblind + Condition + (1|participant), data = SEL_bias_no_labels)
anova(bias_model, bias_model_main_referential_condition)


```

### Exploratory Analyses - Each Referential Condition Separately with main effect of Drug

#### Self
```{r}
SEL_bias_self <- SEL_bias_no_labels %>%
  filter(Condition == "Self")

bias_self_model <- lm(bias ~ group_unblind, data = SEL_bias_self)
summ(bias_self_model, confint = TRUE, digits = 3)
bias_self_model_effects <- ggpredict(bias_self_model, terms = c("group_unblind"))
plot(bias_self_model_effects)
```

#### Friend

```{r}
SEL_bias_friend <- SEL_bias_no_labels %>%
  filter(Condition == "Friend")

bias_friend_model <- lm(bias ~ group_unblind, data = SEL_bias_friend)
summ(bias_friend_model, confint = TRUE, digits = 3)
bias_friend_model_effects <- ggpredict(bias_friend_model, terms = c("group_unblind"))
plot(bias_friend_model_effects)
```

#### Stranger

```{r}
SEL_bias_stranger <- SEL_bias_no_labels %>%
  filter(Condition == "Stranger")

bias_stranger_model <- lm(bias ~ group_unblind, data = SEL_bias_stranger)
summ(bias_stranger_model, confint = TRUE, digits = 3)
bias_stranger_model_effects <- ggpredict(bias_stranger_model, terms = c("group_unblind"))
plot(bias_stranger_model_effects)
```

## Errors to Criterion

```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

SEL_aggregate %>%
  tab_cells(ERR_to_crit) %>%
  tab_cols(group_unblind %nest% Condition %nest% Overall_rule) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()

SEL_aggregate %>%
  group_by(group_unblind, Condition, Overall_rule) %>%
  summarize(Mean = mean(ERR_to_crit),  SD = sd(ERR_to_crit), Median = median(ERR_to_crit), Minimum = min(ERR_to_crit), Maximum = max(ERR_to_crit)) %>%
  ggplot(., aes(x = Overall_rule, y = Mean, fill = group_unblind)) + 
  facet_wrap(~Condition) +
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("Errors to Criterion According to Drug Group, Referential Condition and Rule") + 
  ylab("Errors to Criterion") + 
  xlab ("Rule") + 
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

err_to_crit_model <- lmer(ERR_to_crit ~ group_unblind*Condition*Overall_rule + (1|participant), data = SEL_aggregate_no_labels)
summ(err_to_crit_model, confint = TRUE, digits = 3)
err_to_crit_model_effects <- ggpredict(err_to_crit_model, terms = c("Overall_rule", "group_unblind", "Condition"))
plot(err_to_crit_model_effects) + 
  labs(
    x = "Referential Condition & Rule", 
    y = "Errors to Criterion", 
    title = "Predicted Values of Errors to Criterion",
    colour = "Drug Group"
  )

#Overall p value referential condition ME
err_to_crit_model_no_referential_condition <- lmer(ERR_to_crit ~ group_unblind*Overall_rule + (1|participant), data = SEL_aggregate_no_labels)

anova(err_to_crit_model, err_to_crit_model_no_referential_condition)

#Overall p value referential condition interaction effect with group
err_to_crit_model_referential_condition_ME_drug <- lmer(ERR_to_crit ~ (group_unblind + Condition)*Overall_rule + (1|participant), data = SEL_aggregate_no_labels)

anova(err_to_crit_model, err_to_crit_model_referential_condition_ME_drug)

#Overall p value referential condition interaction effect with rule
err_to_crit_model_referential_condition_ME_rule <- lmer(ERR_to_crit ~ (Overall_rule + Condition)*group_unblind + (1|participant), data = SEL_aggregate_no_labels)

anova(err_to_crit_model, err_to_crit_model_referential_condition_ME_rule)

#Overall p value referential condition interaction effect with rule and group
err_to_crit_model_referential_condition_int_drug_rule <- lmer(ERR_to_crit ~ group_unblind*Overall_rule + Condition + (1|participant), data = SEL_aggregate_no_labels)

anova(err_to_crit_model, err_to_crit_model_referential_condition_int_drug_rule)


```

### Exploratory Analyses - Each Referential Condition Seperately, drug * valence effect

#### Self
```{r}
SEL_self <- SEL_aggregate_no_labels %>%
  filter(Condition == "Self")

self_model <- lmer(ERR_to_crit ~ group_unblind*Overall_rule + (1|participant), data = SEL_self)
summ(self_model, confint = TRUE, digits = 3)
self_model_effects <- ggpredict(self_model, terms = c("Overall_rule", "group_unblind"))
plot(self_model_effects)
```

#### Friend
```{r}
SEL_friend <- SEL_aggregate_no_labels %>%
  filter(Condition == "Friend")

friend_model <- lmer(ERR_to_crit ~ group_unblind*Overall_rule + (1|participant), data = SEL_friend)
summ(friend_model, confint = TRUE, digits = 3)

friend_model_effects <- ggpredict(friend_model, terms = c("Overall_rule", "group_unblind"))
plot(friend_model_effects)+ 
  labs(
    x = "Rule", 
    y = "Errors to Criterion", 
    colour = "Drug Group"
  )


SEL_friend_errors_graph <- SEL_friend %>%
  group_by(group_unblind, Overall_rule) %>%
  summarize(Mean = mean(ERR_to_crit),  SD = sd(ERR_to_crit), Median = median(ERR_to_crit), Minimum = min(ERR_to_crit), Maximum = max(ERR_to_crit)) %>%
  ggplot(., aes(x = Overall_rule, y = Mean, fill = group_unblind)) +
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Errors to Criterion") + 
  xlab ("Rule") + 
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9)) +
  scale_fill_brewer()

print(SEL_friend_errors_graph)

ggsave("SEL_friend_errors_graph.jpeg")

```


#### Stranger
```{r}
SEL_stranger <- SEL_aggregate_no_labels %>%
  filter(Condition == "Stranger")

stranger_model <- lmer(ERR_to_crit ~ group_unblind*Overall_rule + (1|participant), data = SEL_stranger)
summ(stranger_model, confint = TRUE, digits = 3)
stranger_model_effects <- ggpredict(stranger_model, terms = c("Overall_rule", "group_unblind"))
plot(stranger_model_effects)
```

# Categorisation and Recall

## Categorisation

```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

cat_recall_long %>%
  tab_cells(cat_tot_yes) %>%
  tab_cols(group_unblind %nest% condition %nest% valence) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()

cat_recall_long %>%
  group_by(group_unblind, condition, valence) %>%
  summarize(Mean = mean(cat_tot_yes),  SD = sd(cat_tot_yes), Median = median(cat_tot_yes), Minimum = min(cat_tot_yes), Maximum = max(cat_tot_yes)) %>%
  ggplot(., aes(x = valence, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  facet_wrap(~condition) + 
  ggtitle("Number of Words Categorised as Descriptive of the Referential Person According to Drug Group, Condition and Valence") + 
  ylab("Words Categorised as Descriptive (yes)") + 
  xlab ("Referential Condition") + 
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

cat_model <- lmer(cat_tot_yes ~ group_unblind*condition*valence + (1|participant), data = cat_recall_long_no_labels)
summ(cat_model, confint = TRUE, digits = 3)
cat_model_effects <- ggpredict(cat_model, terms = c("valence", "group_unblind", "condition"))
plot(cat_model_effects) + 
  labs(
    x = "Valence", 
    y = "Words Categorised as Descriptive (yes)", 
    title = "Predicted Values of Categorised Words",
    colour = "Drug Group"
  )

```

## Recall
```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

cat_recall_long %>%
  tab_cells(recall_tot_correct) %>%
  tab_cols(group_unblind %nest% condition %nest% valence) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()

cat_recall_long %>%
  group_by(group_unblind, condition, valence) %>%
  summarize(Mean = mean(recall_tot_correct),  SD = sd(recall_tot_correct), Median = median(recall_tot_correct), Minimum = min(recall_tot_correct), Maximum = max(recall_tot_correct)) %>%
  ggplot(., aes(x = valence, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  facet_wrap(~condition) + 
  ggtitle("Correctly Recalled Words According to Drug Group, Condition and Valence") + 
  ylab("Total Number Correctly Recalled Words") + 
  xlab ("Referential Condition") + 
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

recall_model <- lmer(recall_tot_correct ~ group_unblind*condition*valence + (1|participant), data = cat_recall_long_no_labels)
summ(recall_model, confint = TRUE, digits = 3)
recall_model_effects <- ggpredict(recall_model, terms = c("valence", "group_unblind", "condition"))
plot(recall_model_effects) + 
  labs(
    x = "Valence", 
    y = "Number of Correctly Recalled Words", 
    title = "Predicted Values of Correctly Recalled Words",
    colour = "Drug Group"
  )

```


#### Exploratory Analyses - Recall by Referential Condition, Drug Group * Valence effect

##### Self
```{r}
recall_self <- cat_recall_long_no_labels %>%
  filter(condition == "Self")

recall_self_model <- lmer(recall_tot_correct ~ group_unblind*valence + (1|participant), data = recall_self)
summ(recall_self_model, confint = TRUE, digits = 3)
recall_self_model_effects <- ggpredict(recall_self_model, terms = c("valence", "group_unblind"))
plot(recall_self_model_effects) + 
  labs(
    x = "Valence", 
    y = "Number of Correctly Recalled Words", 
    title = "Predicted Values of Correctly Recalled Words",
    colour = "Drug Group"
  )
```

##### Other
```{r}
recall_other <- cat_recall_long_no_labels %>%
  filter(condition == "Other")

recall_other_model <- lmer(recall_tot_correct ~ group_unblind*valence + (1|participant), data = recall_other)
summ(recall_other_model, confint = TRUE, digits = 3)
recall_other_model_effects <- ggpredict(recall_other_model, terms = c("valence", "group_unblind"))
plot(recall_other_model_effects) + 
  labs(
    x = "Valence", 
    y = "Number of Correctly Recalled Words", 
    title = "Predicted Values of Correctly Recalled Words",
    colour = "Drug Group"
  )

recall_other_graph <- recall_other %>%
  group_by(group_unblind, valence) %>%
  summarize(Mean = mean(recall_tot_correct),  SD = sd(recall_tot_correct), Median = median(recall_tot_correct), Minimum = min(recall_tot_correct), Maximum = max(recall_tot_correct)) %>%
  ggplot(., aes(x = valence, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ylab("Total Correct Recalled Words") + 
  xlab ("Valence") + 
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))+
  scale_fill_brewer()

print(recall_other_graph)
ggsave("recall_other_graph.jpeg")

```



### Assessing the effect of English as a first language on Recall


```{r}
#Merging required demographic data with task data
english_language <- demo_qs_no_labels %>%
  select(participant, english_first_language)

cat_recall_long_no_labels <- merge(cat_recall_long_no_labels, english_language, by = "participant")


recall_adj_model <- lmer(recall_tot_correct ~ group_unblind*condition*valence + english_first_language + (1|participant), data = cat_recall_long_no_labels)
summ(recall_adj_model, confint = TRUE, digits = 3)
recall_adj_model_effects <- ggpredict(recall_adj_model, terms = c("valence", "group_unblind", "condition"))
plot(recall_adj_model_effects) + 
  labs(
    x = "Valence", 
    y = "Number of Correctly Recalled Words", 
    title = "Predicted Values of Correctly Recalled Words",
    colour = "Drug Group"
  )

```

# Go No-Go Associative Learning Task

Initially conducting this with excluded participants as a priori exclusion criteria specified. However as a large number of participants are excluded due to this (final N = 31) I have also repeated the analysis with participants excluded on the basis of bias and discrimination scores included. One participant was excluded from all analyses as they did not press the spacebar throughout and therefore clearly did not engage in the task. This was done in the cleaning script for this task. 

## D'
```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

GNAT_agg_long_excluded %>%
  tab_cells(d) %>%
  tab_cols(group_unblind %nest% ref_condition %nest% valence) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()

GNAT_agg_long_excluded %>%
  group_by(group_unblind, ref_condition, valence) %>%
  summarize(Mean = mean(d),  SD = sd(d), Median = median(d), Minimum = min(d), Maximum = max(d)) %>%
  ggplot(., aes(x = valence, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  facet_wrap(~ref_condition) + 
  ggtitle("d' According to Drug Group, Condition and Valence") + 
  ylab("d’") + 
  xlab ("Referential Condition") + 
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

d_model <- lmer(d ~ group_unblind*ref_condition*valence + (1|participant), data = GNAT_agg_long_excluded_no_labels)
summ(d_model, confint = TRUE, digits = 3)
d_model_effects <- ggpredict(d_model, terms = c("valence", "group_unblind", "ref_condition"))
plot(d_model_effects) + 
  labs(
    x = "Valence", 
    y = "D'", 
    title = "Predicted Values of D'",
    colour = "Drug Group"
  )


GNAT_agg_long_excluded %>%
  group_by(ref_condition, valence) %>%
  summarize(Mean = mean(d),  SD = sd(d), Median = median(d), Minimum = min(d), Maximum = max(d)) %>%
  ggplot(., aes(x = valence, y = Mean, fill = ref_condition)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("d' According to Condition and Valence") + 
  ylab("d’") + 
  xlab ("Valence") + 
  labs(fill = "Condition") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))
```

## Hits

```{r}

options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

GNAT_agg_long_excluded %>%
  tab_cells(hit_prop) %>%
  tab_cols(group_unblind %nest% ref_condition %nest% valence) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()

GNAT_agg_long_excluded %>%
  group_by(group_unblind, ref_condition, valence) %>%
  summarize(Mean = mean(hit_prop),  SD = sd(hit_prop), Median = median(hit_prop), Minimum = min(hit_prop), Maximum = max(hit_prop)) %>%
  ggplot(., aes(x = valence, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  facet_wrap(~ref_condition) + 
  ggtitle("Hits According to Drug Group, Condition and Valence") + 
  ylab("Hits (proportion)") + 
  xlab ("Referential Condition") + 
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

Hits_model <- lmer(hit_prop ~ group_unblind*ref_condition*valence + (1|participant), data = GNAT_agg_long_excluded_no_labels)
summ(Hits_model, confint = TRUE, digits = 3)
Hits_model_effects <- ggpredict(Hits_model, terms = c("valence", "group_unblind", "ref_condition"))
plot(Hits_model_effects) + 
  labs(
    x = "Valence", 
    y = "Hits (Proportion)", 
    title = "Predicted Values of Hits",
    colour = "Drug Group"
  )
```
## False Alarms
```{r}

options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

GNAT_agg_long_excluded %>%
  tab_cells(fa_prop) %>%
  tab_cols(group_unblind %nest% ref_condition %nest% valence) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()

GNAT_agg_long_excluded %>%
  group_by(group_unblind, ref_condition, valence) %>%
  summarize(Mean = mean(fa_prop),  SD = sd(fa_prop), Median = median(fa_prop), Minimum = min(fa_prop), Maximum = max(fa_prop)) %>%
  ggplot(., aes(x = valence, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  facet_wrap(~ref_condition) + 
  ggtitle("False Alarms According to Drug Group, Condition and Valence") + 
  ylab("False Alarms (proportion)") + 
  xlab ("Referential Condition") + 
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

FA_model <- lmer(fa_prop ~ group_unblind*ref_condition*valence + (1|participant), data = GNAT_agg_long_excluded_no_labels)
summ(FA_model, confint = TRUE, digits = 3)
FA_model_effects <- ggpredict(FA_model, terms = c("valence", "group_unblind", "ref_condition"))
plot(FA_model_effects) + 
  labs(
    x = "Valence", 
    y = "False Alarms (proportion)", 
    title = "Predicted Values of False Alarms",
    colour = "Drug Group"
  )
```

## Repeating Analyses with Excluded Participants Included (Sensitivity Analysis)

In the following analyses participants initially excluded on the basis of bias and discrimination scores are included (overall n = 40). One participant was excluded from all analyses as they did not press the spacebar throughout and therefore clearly did not engage in the task. 

### D' (all participants)
```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

GNAT_agg_long %>%
  tab_cells(d) %>%
  tab_cols(group_unblind %nest% ref_condition %nest% valence) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()

GNAT_agg_long %>%
  group_by(group_unblind, ref_condition, valence) %>%
  summarize(Mean = mean(d),  SD = sd(d), Median = median(d), Minimum = min(d), Maximum = max(d)) %>%
  ggplot(., aes(x = valence, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  facet_wrap(~ref_condition) + 
  ggtitle("d' According to Drug Group, Condition and Valence") + 
  ylab("d’") + 
  xlab ("Referential Condition") + 
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

d_model <- lmer(d ~ group_unblind*ref_condition*valence + (1|participant), data = GNAT_agg_long_no_labels)
summ(d_model, confint = TRUE, digits = 3)
d_model_effects <- ggpredict(d_model, terms = c("valence", "group_unblind", "ref_condition"))
plot(d_model_effects) + 
  labs(
    x = "Valence", 
    y = "D'", 
    title = "Predicted Values of D'",
    colour = "Drug Group"
  )
```

### Hits (all participants)

```{r}

options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

GNAT_agg_long %>%
  tab_cells(hit_prop) %>%
  tab_cols(group_unblind %nest% ref_condition %nest% valence) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()

GNAT_agg_long %>%
  group_by(group_unblind, ref_condition, valence) %>%
  summarize(Mean = mean(hit_prop),  SD = sd(hit_prop), Median = median(hit_prop), Minimum = min(hit_prop), Maximum = max(hit_prop)) %>%
  ggplot(., aes(x = valence, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  facet_wrap(~ref_condition) + 
  ggtitle("Hits According to Drug Group, Condition and Valence") + 
  ylab("Hits (proportion)") + 
  xlab ("Referential Condition") + 
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

Hits_model <- lmer(hit_prop ~ group_unblind*ref_condition*valence + (1|participant), data = GNAT_agg_long_no_labels)
summ(Hits_model, confint = TRUE, digits = 3)
Hits_model_effects <- ggpredict(Hits_model, terms = c("valence", "group_unblind", "ref_condition"))
plot(Hits_model_effects) + 
  labs(
    x = "Valence", 
    y = "Hits (Proportion)", 
    title = "Predicted Values of Hits",
    colour = "Drug Group"
  )
```
## False Alarms (all participants)
```{r}

options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

GNAT_agg_long %>%
  tab_cells(fa_prop) %>%
  tab_cols(group_unblind %nest% ref_condition %nest% valence) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()

GNAT_agg_long %>%
  group_by(group_unblind, ref_condition, valence) %>%
  summarize(Mean = mean(fa_prop),  SD = sd(fa_prop), Median = median(fa_prop), Minimum = min(fa_prop), Maximum = max(fa_prop)) %>%
  ggplot(., aes(x = valence, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  facet_wrap(~ref_condition) + 
  ggtitle("False Alarms According to Drug Group, Condition and Valence") + 
  ylab("False Alarms (proportion)") + 
  xlab ("Referential Condition") + 
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

FA_model <- lmer(fa_prop ~ group_unblind*ref_condition*valence + (1|participant), data = GNAT_agg_long_no_labels)
summ(FA_model, confint = TRUE, digits = 3)
FA_model_effects <- ggpredict(FA_model, terms = c("valence", "group_unblind", "ref_condition"))
plot(FA_model_effects) + 
  labs(
    x = "Valence", 
    y = "False Alarms (proportion)", 
    title = "Predicted Values of False Alarms",
    colour = "Drug Group"
  )
```

# Associative Learning

```{r}
#Collapsing across matching & non-matching trials
associative_long_no_labels_collapsed <- associative_long_no_labels %>%
  group_by(participant, group_unblind, Task, stimuli) %>%
  summarise_at(c("mean_rt_mult", "prop_acc"), mean)

associative_long_collapsed <- associative_long_no_labels_collapsed

associative_long_collapsed = apply_labels(associative_long_collapsed,
                       group_unblind = "Drug Group",
                       prop_acc = "Accuracy (%)",
                       mean_rt_mult = "Reaction Times (ms)",
                       stimuli = "Stimuli")
```


## Self

### Accuracy (%)
```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

associative_long_collapsed_self <- associative_long_collapsed %>%
  filter(Task == "Self")

associative_long_collapsed_self_no_labels <- associative_long_no_labels_collapsed %>%
  filter(Task == "Self")

associative_long_collapsed_self %>%
  tab_cells(prop_acc) %>%
  tab_cols(group_unblind %nest% stimuli) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()


associative_long_collapsed_self %>%
  group_by(group_unblind, stimuli) %>%
  summarize(Mean = mean(prop_acc),  SD = sd(prop_acc), Median = median(prop_acc), Minimum = min(prop_acc), Maximum = max(prop_acc)) %>%
  ggplot(., aes(x = stimuli, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("Proportion of Accurate Trials in the Self Associative Learning Task") + 
  ylab("Average Proportion of Accurate Trials") + 
  xlab("Stimuli") +
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

associative_self_acc_model <- lmer(prop_acc ~ group_unblind*stimuli + (1|participant), data = associative_long_collapsed_self_no_labels)
summ(associative_self_acc_model, confint = TRUE, digits = 3)
associative_self_acc_model_effects <- ggpredict(associative_self_acc_model, terms = c("stimuli", "group_unblind"))
plot(associative_self_acc_model_effects) +
  labs(
    title = "Predicted Values of Accuracy (%) from Mixed Regression Model",
    y = "Accuracy (%)",
    x = "Stimuli",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of stimuli 
associative_self_acc_model_no_stimuli<- lmer(prop_acc ~ group_unblind + (1|participant), data = associative_long_collapsed_self_no_labels)
anova(associative_self_acc_model, associative_self_acc_model_no_stimuli)

#Getting overall p-value for interaction effect of stimuli
associative_self_acc_model_main_stimuli <- lmer(prop_acc ~ group_unblind + stimuli + (1|participant), data = associative_long_collapsed_self_no_labels)
anova(associative_self_acc_model, associative_self_acc_model_main_stimuli)

```
### Reaction Times (ms)
```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

associative_long_collapsed_self <- associative_long_collapsed %>%
  filter(Task == "Self")

associative_long_collapsed_self_no_labels <- associative_long_no_labels_collapsed %>%
  filter(Task == "Self")

associative_long_collapsed_self %>%
  tab_cells(mean_rt_mult) %>%
  tab_cols(group_unblind %nest% stimuli) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()


associative_long_collapsed_self %>%
  group_by(group_unblind, stimuli) %>%
  summarize(Mean = mean(mean_rt_mult),  SD = sd(mean_rt_mult), Median = median(mean_rt_mult), Minimum = min(mean_rt_mult), Maximum = max(mean_rt_mult)) %>%
  ggplot(., aes(x = stimuli, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("Reaction Times (ms) in the Self Associative Learning Task") + 
  ylab("Average Reaction Times (ms)") + 
  xlab("Stimuli") +
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

associative_self_rt_model <- lmer(mean_rt_mult ~ group_unblind*stimuli + (1|participant), data = associative_long_collapsed_self_no_labels)
summ(associative_self_rt_model, confint = TRUE, digits = 3)
associative_self_rt_model_effects <- ggpredict(associative_self_rt_model, terms = c("stimuli", "group_unblind"))
plot(associative_self_rt_model_effects) +
  labs(
    title = "Predicted Values of Reaction Times from Mixed Regression Model",
    y = "Reaction Times (ms)",
    x = "Stimuli",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of stimuli 
associative_self_rt_model_no_stimuli<- lmer(mean_rt_mult ~ group_unblind + (1|participant), data = associative_long_collapsed_self_no_labels)
anova(associative_self_rt_model, associative_self_rt_model_no_stimuli)

#Getting overall p-value for interaction effect of stimuli
associative_self_rt_model_main_stimuli <- lmer(mean_rt_mult ~ group_unblind + stimuli + (1|participant), data = associative_long_collapsed_self_no_labels)
anova(associative_self_rt_model, associative_self_rt_model_main_stimuli)

```

## Emotion
### Accuracy (%)
```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

associative_long_collapsed_valence <- associative_long_collapsed %>%
  filter(Task == "Valence")

associative_long_collapsed_valence_no_labels <- associative_long_no_labels_collapsed %>%
  filter(Task == "Valence")

associative_long_collapsed_valence %>%
  tab_cells(prop_acc) %>%
  tab_cols(group_unblind %nest% stimuli) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()


associative_long_collapsed_valence %>%
  group_by(group_unblind, stimuli) %>%
  summarize(Mean = mean(prop_acc),  SD = sd(prop_acc), Median = median(prop_acc), Minimum = min(prop_acc), Maximum = max(prop_acc)) %>%
  ggplot(., aes(x = stimuli, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("Proportion of Accurate Trials in the Valence Associative Learning Task") + 
  ylab("Average Proportion of Accurate Trials") + 
  xlab("Stimuli") +
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

associative_valence_acc_model <- lmer(prop_acc ~ group_unblind*stimuli + (1|participant), data = associative_long_collapsed_valence_no_labels)
summ(associative_valence_acc_model, confint = TRUE, digits = 3)
associative_valence_acc_model_effects <- ggpredict(associative_valence_acc_model, terms = c("stimuli", "group_unblind"))
plot(associative_valence_acc_model_effects) +
  labs(
    title = "Predicted Values of Accuracy (%) from Mixed Regression Model",
    y = "Accuracy (%)",
    x = "Stimuli",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of stimuli 
associative_valence_acc_model_no_stimuli<- lmer(prop_acc ~ group_unblind + (1|participant), data = associative_long_collapsed_valence_no_labels)
anova(associative_valence_acc_model, associative_valence_acc_model_no_stimuli)

#Getting overall p-value for interaction effect of stimuli
associative_valence_acc_model_main_stimuli <- lmer(prop_acc ~ group_unblind + stimuli + (1|participant), data = associative_long_collapsed_valence_no_labels)
anova(associative_valence_acc_model, associative_valence_acc_model_main_stimuli)

```
### Reaction Times (ms)
```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

associative_long_collapsed_valence <- associative_long_collapsed %>%
  filter(Task == "Valence")

associative_long_collapsed_valence_no_labels <- associative_long_no_labels_collapsed %>%
  filter(Task == "Valence")

associative_long_collapsed_valence %>%
  tab_cells(mean_rt_mult) %>%
  tab_cols(group_unblind %nest% stimuli) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()


associative_long_collapsed_valence %>%
  group_by(group_unblind, stimuli) %>%
  summarize(Mean = mean(mean_rt_mult),  SD = sd(mean_rt_mult), Median = median(mean_rt_mult), Minimum = min(mean_rt_mult), Maximum = max(mean_rt_mult)) %>%
  ggplot(., aes(x = stimuli, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("Reaction Times (ms) in the Valence Associative Learning Task") + 
  ylab("Average Reaction Times (ms)") + 
  xlab("Stimuli") +
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

associative_valence_rt_model <- lmer(mean_rt_mult ~ group_unblind*stimuli + (1|participant), data = associative_long_collapsed_valence_no_labels)
summ(associative_valence_rt_model, confint = TRUE, digits = 3)
associative_valence_rt_model_effects <- ggpredict(associative_valence_rt_model, terms = c("stimuli", "group_unblind"))
plot(associative_valence_rt_model_effects) +
  labs(
    title = "Predicted Values of Reaction Times from Mixed Regression Model",
    y = "Reaction Times (ms)",
    x = "Stimuli",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of stimuli 
associative_valence_rt_model_no_stimuli<- lmer(mean_rt_mult ~ group_unblind + (1|participant), data = associative_long_collapsed_valence_no_labels)
anova(associative_valence_rt_model, associative_valence_rt_model_no_stimuli)

#Getting overall p-value for interaction effect of stimuli
associative_valence_rt_model_main_stimuli <- lmer(mean_rt_mult ~ group_unblind + stimuli + (1|participant), data = associative_long_collapsed_valence_no_labels)
anova(associative_valence_rt_model, associative_valence_rt_model_main_stimuli)
```

## Reward
### Accuracy (%)
```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

associative_long_collapsed_reward <- associative_long_collapsed %>%
  filter(Task == "Reward")

associative_long_collapsed_reward_no_labels <- associative_long_no_labels_collapsed %>%
  filter(Task == "Reward")

associative_long_collapsed_reward %>%
  tab_cells(prop_acc) %>%
  tab_cols(group_unblind %nest% stimuli) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()


associative_long_collapsed_reward %>%
  group_by(group_unblind, stimuli) %>%
  summarize(Mean = mean(prop_acc),  SD = sd(prop_acc), Median = median(prop_acc), Minimum = min(prop_acc), Maximum = max(prop_acc)) %>%
  ggplot(., aes(x = stimuli, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("Proportion of Accurate Trials in the Reward Associative Learning Task") + 
  ylab("Average Proportion of Accurate Trials") + 
  xlab("Stimuli") +
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

associative_reward_acc_model <- lmer(prop_acc ~ group_unblind*stimuli + (1|participant), data = associative_long_collapsed_reward_no_labels)
summ(associative_reward_acc_model, confint = TRUE, digits = 3)
associative_reward_acc_model_effects <- ggpredict(associative_reward_acc_model, terms = c("stimuli", "group_unblind"))
plot(associative_reward_acc_model_effects) +
  labs(
    title = "Predicted Values of Accuracy (%) from Mixed Regression Model",
    y = "Accuracy (%)",
    x = "Stimuli",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of stimuli 
associative_reward_acc_model_no_stimuli<- lmer(prop_acc ~ group_unblind + (1|participant), data = associative_long_collapsed_reward_no_labels)
anova(associative_reward_acc_model, associative_reward_acc_model_no_stimuli)

#Getting overall p-value for interaction effect of stimuli
associative_reward_acc_model_main_stimuli <- lmer(prop_acc ~ group_unblind + stimuli + (1|participant), data = associative_long_collapsed_reward_no_labels)
anova(associative_reward_acc_model, associative_reward_acc_model_main_stimuli)
```
### Reaction Times (ms)
```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

associative_long_collapsed_reward <- associative_long_collapsed %>%
  filter(Task == "Reward")

associative_long_collapsed_reward_no_labels <- associative_long_no_labels_collapsed %>%
  filter(Task == "Reward")

associative_long_collapsed_reward %>%
  tab_cells(mean_rt_mult) %>%
  tab_cols(group_unblind %nest% stimuli) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot()%>%
  drop_empty_columns()


associative_long_collapsed_reward %>%
  group_by(group_unblind, stimuli) %>%
  summarize(Mean = mean(mean_rt_mult),  SD = sd(mean_rt_mult), Median = median(mean_rt_mult), Minimum = min(mean_rt_mult), Maximum = max(mean_rt_mult)) %>%
  ggplot(., aes(x = stimuli, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("Reaction Times (ms) in the Reward Associative Learning Task") + 
  ylab("Average Reaction Times (ms)") + 
  xlab("Stimuli") +
  labs(fill = "Drug Group") +
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9))

associative_reward_rt_model <- lmer(mean_rt_mult ~ group_unblind*stimuli + (1|participant), data = associative_long_collapsed_reward_no_labels)
summ(associative_reward_rt_model, confint = TRUE, digits = 3)
associative_reward_rt_model_effects <- ggpredict(associative_reward_rt_model, terms = c("stimuli", "group_unblind"))
plot(associative_reward_rt_model_effects) +
  labs(
    title = "Predicted Values of Reaction Times from Mixed Regression Model",
    y = "Reaction Times (ms)",
    x = "Stimuli",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of stimuli 
associative_reward_rt_model_no_stimuli<- lmer(mean_rt_mult ~ group_unblind + (1|participant), data = associative_long_collapsed_reward_no_labels)
anova(associative_reward_rt_model, associative_reward_rt_model_no_stimuli)

#Getting overall p-value for interaction effect of stimuli
associative_reward_rt_model_main_stimuli <- lmer(mean_rt_mult ~ group_unblind + stimuli + (1|participant), data = associative_long_collapsed_reward_no_labels)
anova(associative_reward_rt_model, associative_reward_rt_model_main_stimuli)
```



# Blinding

## Group Guesses

```{r}
demo_qs %>%
  tab_cols(group_unblind) %>%
  tab_cells(group_guess_participant, group_guess_researcher) %>%
  tab_stat_cases(total_row_position = c("none"))%>%
  tab_cells(drug_guess_certainty_participant, placebo_guess_certainty_participant, drug_guess_certainty_researcher, placebo_guess_certainty_researcher) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot() %>%
  set_caption("Group Guesses")

```
### Participant Group Guess
```{r}
participant_group_guess_table <- as.data.frame(table(demo_qs$group_unblind, demo_qs$group_guess_participant))

ggplot(data = participant_group_guess_table, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat="identity", position = position_dodge()) +
  xlab("Drug Group") +
  ylab("Total Number") +
  labs(fill = "Group Guess")


chisq.test(demo_qs$group_guess_participant, demo_qs$group_unblind)
```
### Participant Citalopram Certainty
```{r}

demo_qs %>%
  group_by(group_unblind) %>%
  summarize(Mean = mean(drug_guess_certainty_participant, na.rm = TRUE),  SD = sd(drug_guess_certainty_participant, na.rm = TRUE), Median = median(drug_guess_certainty_participant, na.rm = TRUE), Minimum = min(drug_guess_certainty_participant, na.rm = TRUE), Maximum = max(drug_guess_certainty_participant, na.rm = TRUE)) %>%
  ggplot(., aes(x = group_unblind, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("Participant Citalopram Group Certainty") + 
  ylab("Certainty") + 
  xlab("Drug Group") + 
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9)) + 
  theme(legend.position = "none")

t.test(demo_qs$drug_guess_certainty_participant~demo_qs$group_unblind)

```
### Participant Placebo Certainty
```{r}
demo_qs %>%
  group_by(group_unblind) %>%
  summarize(Mean = mean(placebo_guess_certainty_participant, na.rm = TRUE),  SD = sd(placebo_guess_certainty_participant, na.rm = TRUE), Median = median(placebo_guess_certainty_participant, na.rm = TRUE), Minimum = min(placebo_guess_certainty_participant, na.rm = TRUE), Maximum = max(placebo_guess_certainty_participant, na.rm = TRUE)) %>%
  ggplot(., aes(x = group_unblind, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("Participant Placebo Group Certainty") + 
  ylab("Certainty") + 
  xlab("Drug Group") + 
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9)) + 
  theme(legend.position = "none")

t.test(demo_qs$placebo_guess_certainty_participant~demo_qs$group_unblind)
```



### Researcher Group Guess
```{r}
researcher_group_guess_table <- as.data.frame(table(demo_qs$group_unblind, demo_qs$group_guess_researcher))

ggplot(data = researcher_group_guess_table, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat="identity", position = position_dodge()) +
  xlab("Drug Group") +
  ylab("Total Number") +
  labs(fill = "Group Guess")

chisq.test(demo_qs$group_guess_researcher, demo_qs$group_unblind)
```

### Researcher Citalopram Certainty
```{r}

demo_qs %>%
  group_by(group_unblind) %>%
  summarize(Mean = mean(drug_guess_certainty_researcher, na.rm = TRUE),  SD = sd(drug_guess_certainty_researcher, na.rm = TRUE), Median = median(drug_guess_certainty_researcher, na.rm = TRUE), Minimum = min(drug_guess_certainty_researcher, na.rm = TRUE), Maximum = max(drug_guess_certainty_researcher, na.rm = TRUE)) %>%
  ggplot(., aes(x = group_unblind, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("researcher Citalopram Group Certainty") + 
  ylab("Certainty") + 
  xlab("Drug Group") + 
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9)) + 
  theme(legend.position = "none")

t.test(demo_qs$drug_guess_certainty_researcher~demo_qs$group_unblind)

```
### Researcher Placebo Certainty
```{r}
demo_qs %>%
  group_by(group_unblind) %>%
  summarize(Mean = mean(placebo_guess_certainty_researcher, na.rm = TRUE),  SD = sd(placebo_guess_certainty_researcher, na.rm = TRUE), Median = median(placebo_guess_certainty_researcher, na.rm = TRUE), Minimum = min(placebo_guess_certainty_researcher, na.rm = TRUE), Maximum = max(placebo_guess_certainty_researcher, na.rm = TRUE)) %>%
  ggplot(., aes(x = group_unblind, y = Mean, fill = group_unblind)) + 
  geom_bar(stat = 'identity', position = position_dodge()) + 
  ggtitle("researcher Placebo Group Certainty") + 
  ylab("Certainty") + 
  xlab("Drug Group") + 
  theme_minimal() + 
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = .2, position=position_dodge(.9)) + 
  theme(legend.position = "none")

t.test(demo_qs$placebo_guess_certainty_researcher~demo_qs$group_unblind)
```

## Side Effects

```{r}
options(expss.digits = 2, expss.output = 'rnotebook') #Setting digits displayed to 2 decimal places

state_mood_long %>%
  tab_cols(group_unblind %nest% Time) %>%
  tab_cells(SideEffects_Q1, SideEffects_Q2, SideEffects_Q3, SideEffects_Q4, SideEffects_Q5, SideEffects_Q6) %>%
  tab_stat_mean_sd_n() %>%
  tab_pivot() %>%
  set_caption("Self-Reported Side Effects by Group and Time")
```

### Nausea
```{r}
SideEffects_Q1_group_time <- lmer(SideEffects_Q1 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(SideEffects_Q1_group_time, confint = TRUE, digits = 3)
SideEffects_Q1_group_time_effects <- ggpredict(SideEffects_Q1_group_time, terms = c("Time", "group_unblind"))
plot(SideEffects_Q1_group_time_effects) +
  labs(
    title = "Predicted Values of Nausea from Mixed Regression Model",
    y = "Nausea",
    x = "Timepoint",
    colour = "Drug Group"
  )


#Contrasts for drug effect at post-drug and post-testing
contrast_SideEffects_Q1_group_time<- emmeans(SideEffects_Q1_group_time, ~ group_unblind*Time)
contrast(contrast_SideEffects_Q1_group_time, "revpairwise", by = "Time", adjust = "bonferroni")

#Getting overall p-value for main effect of time 
SideEffects_Q1_no_time <- lmer(SideEffects_Q1 ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(SideEffects_Q1_group_time, SideEffects_Q1_no_time)

#Getting overall p-value for interaction effect of time
SideEffects_Q1_main_time <- lmer(SideEffects_Q1 ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(SideEffects_Q1_group_time, SideEffects_Q1_main_time)
```


### Dizziness
```{r}
SideEffects_Q2_group_time <- lmer(SideEffects_Q2 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(SideEffects_Q2_group_time, confint = TRUE, digits = 3)
SideEffects_Q2_group_time_effects <- ggpredict(SideEffects_Q2_group_time, terms = c("Time", "group_unblind"))
plot(SideEffects_Q2_group_time_effects) +
  labs(
    title = "Predicted Values of Dizziness from Mixed Regression Model",
    y = "Dizziness",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of time 
SideEffects_Q2_no_time <- lmer(SideEffects_Q2 ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(SideEffects_Q2_group_time, SideEffects_Q2_no_time)

#Getting overall p-value for interaction effect of time
SideEffects_Q2_main_time <- lmer(SideEffects_Q2 ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(SideEffects_Q2_group_time, SideEffects_Q2_main_time)
```


### Dry Mouth
```{r}
SideEffects_Q3_group_time <- lmer(SideEffects_Q3 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(SideEffects_Q3_group_time, confint = TRUE, digits = 3)
SideEffects_Q3_group_time_effects <- ggpredict(SideEffects_Q3_group_time, terms = c("Time", "group_unblind"))
plot(SideEffects_Q3_group_time_effects) +
  labs(
    title = "Predicted Values of Dry Mouth from Mixed Regression Model",
    y = "Dry Mouth",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of time 
SideEffects_Q3_no_time <- lmer(SideEffects_Q3 ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(SideEffects_Q3_group_time, SideEffects_Q3_no_time)

#Getting overall p-value for interaction effect of time
SideEffects_Q3_main_time <- lmer(SideEffects_Q3 ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(SideEffects_Q3_group_time, SideEffects_Q3_main_time)
```


### Headache
```{r}
SideEffects_Q4_group_time <- lmer(SideEffects_Q4 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(SideEffects_Q4_group_time, confint = TRUE, digits = 3)
SideEffects_Q4_group_time_effects <- ggpredict(SideEffects_Q4_group_time, terms = c("Time", "group_unblind"))
plot(SideEffects_Q4_group_time_effects) +
  labs(
    title = "Predicted Values of Dizziness from Mixed Regression Model",
    y = "Headache",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of time 
SideEffects_Q4_no_time <- lmer(SideEffects_Q4 ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(SideEffects_Q4_group_time, SideEffects_Q4_no_time)

#Getting overall p-value for interaction effect of time
SideEffects_Q4_main_time <- lmer(SideEffects_Q4 ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(SideEffects_Q4_group_time, SideEffects_Q4_main_time)
```


### Alertness
```{r}
SideEffects_Q5_group_time <- lmer(SideEffects_Q5 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(SideEffects_Q5_group_time, confint = TRUE, digits = 3)
SideEffects_Q5_group_time_effects <- ggpredict(SideEffects_Q5_group_time, terms = c("Time", "group_unblind"))
plot(SideEffects_Q5_group_time_effects) +
  labs(
    title = "Predicted Values of Alertness from Mixed Regression Model",
    y = "Alertness",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of time 
SideEffects_Q5_no_time <- lmer(SideEffects_Q5 ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(SideEffects_Q5_group_time, SideEffects_Q5_no_time)

#Getting overall p-value for interaction effect of time
SideEffects_Q5_main_time <- lmer(SideEffects_Q5 ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(SideEffects_Q5_group_time, SideEffects_Q5_main_time)
```


### Agitation
```{r}
SideEffects_Q6_group_time <- lmer(SideEffects_Q6 ~ group_unblind*Time + (1|participant), data = state_mood_long_no_labels)
summ(SideEffects_Q6_group_time, confint = TRUE, digits = 3)
SideEffects_Q6_group_time_effects <- ggpredict(SideEffects_Q6_group_time, terms = c("Time", "group_unblind"))
plot(SideEffects_Q6_group_time_effects) +
  labs(
    title = "Predicted Values of Agitation from Mixed Regression Model",
    y = "Agitation",
    x = "Timepoint",
    colour = "Drug Group"
  )

#Getting overall p-value for main effect of time 
SideEffects_Q6_no_time <- lmer(SideEffects_Q6 ~ group_unblind + (1|participant), data = state_mood_long_no_labels)
anova(SideEffects_Q6_group_time, SideEffects_Q6_no_time)

#Getting overall p-value for interaction effect of time
SideEffects_Q6_main_time <- lmer(SideEffects_Q6 ~ group_unblind + Time + (1|participant), data = state_mood_long_no_labels)
anova(SideEffects_Q6_group_time, SideEffects_Q6_main_time)
```