SST_test_intervention_effects.Rmd

---
title: "SST_test_intervention_effects"
author: "Ben Smith"
date: "2023-08-30"
output: html_document

---


```{r setup, include=FALSE}
library(dplyr)
source("SST_processing.R")

knitr::opts_chunk$set(echo = TRUE)
#data_path <- "../../../files/data/"
data_path <- config::get("dev_analysis_data_dir")

full_dataset_with_groups <- readr::read_csv(paste0(data_path,"full_dataset_aim3.csv"))
full_dataset_with_groups$intervention_group<-factor(full_dataset_with_groups$intervention_group,levels=c("willamette","umpqua","mckenzie"))
full_dataset_with_groups$age365_z<-scale(full_dataset_with_groups$age365)


dropbox_file_dir = config::get("dev_analysis_data_dir")
sst_all_data_filepath <- paste0(dropbox_file_dir,"sst_behavioral_data_all.csv")

sst_all_data_raw <- readr::read_csv(sst_all_data_filepath)
sst_all_data_clean <-sst_all_data_raw %>% clean_sst_data

all_ppt_data <- readr::read_csv(paste0(dropbox_file_dir,"data_by_ppt_all_sessions.csv"))
#all_ppt_data <- readr::read_csv(paste0(dropbox_file_dir,"data_by_wave_ppt.csv"))
```

```{r}


healthy_vs_unhealthy <- sst_all_data_clean %>% filter(condition %in% c("CorrectGo","FailedStop")) %>% dplyr::select(subid,waveid,trial_n,condition,reaction_time_clean,health_condition_label) %>%
  group_by(subid,waveid,condition,health_condition_label) %>%
  summarize(mean_rt=mean(reaction_time_clean,na.rm=TRUE), mean_trial_n=mean(trial_n)) %>%
  pivot_wider(names_from="health_condition_label",values_from = c("mean_rt","mean_trial_n"))

healthy_vs_unhealthy$healthy_minus_unhealthy_rt <- healthy_vs_unhealthy$mean_rt_Healthy-healthy_vs_unhealthy$mean_rt_Unhealthy
healthy_vs_unhealthy$healthy_minus_unhealthy_trial_n <- healthy_vs_unhealthy$mean_trial_n_Healthy-healthy_vs_unhealthy$mean_trial_n_Unhealthy

value_cols <- colnames(healthy_vs_unhealthy)[!(colnames(healthy_vs_unhealthy) %in% c("subid", "waveid","condition"))]
results_by_subj<-healthy_vs_unhealthy %>% pivot_wider(id_cols = c("subid","condition"),names_from="waveid",names_prefix="waveid_",values_from = value_cols)

results_by_subj$mean_rt_Healthy_w2_m_w1<-results_by_subj$mean_rt_Healthy_waveid_2-results_by_subj$mean_rt_Healthy_waveid_1

results_by_subj$mean_rt_Unhealthy_w2_m_w1<-results_by_subj$mean_rt_Unhealthy_waveid_2-results_by_subj$mean_rt_Unhealthy_waveid_1

results_by_subj$mean_rt_Healthy_m_Unhealthy_w2_m_w1<-(
  results_by_subj$mean_rt_Healthy_w2_m_w1-results_by_subj$mean_rt_Unhealthy_w2_m_w1
)

```


# Behavioral predictions

## BRT theory:
    - BR Training (vs. Control) decreases RT to healthy foods and increases RT to unhealthy foods; because the SST task calibrates difficulty level, these two will be difficult to pull apart but we should see an increase Unhealthy RT - Healthy RT
    
    

    
```{r}
library(tidyverse)


#all_ppt_data %>% select(SID, session_id, SST_mean_ssrt_0) %>%  


results_by_subj_int <- results_by_subj %>% merge(full_dataset_with_groups %>% dplyr::select(SID,intervention_group),by.x = "subid",by.y="SID",all = FALSE)
results_by_subj_int$intervention_group<-factor(results_by_subj_int$intervention_group,levels=c("willamette","umpqua","mckenzie"))


results_by_subj_brt_test <- results_by_subj_int %>% filter(intervention_group %in% c("willamette", "mckenzie"))
results_by_subj_brt<-results_by_subj_int %>% filter(intervention_group %in% c("mckenzie"))
```
    
### Directional test
First of all, is the prediction _directionally_ right, i.e., does RT to healthy food decrease, unhealthy increase, and Unhealth-Healthy increase even more? We're mainly looking at directions; the hypothesis is best tested with intergroup comparison via a MLM. But we'll look twice: once across the whole dataset; once at BRT specifically


```{r}
library(tidyverse)


#all_ppt_data %>% select(SID, session_id, SST_mean_ssrt_0) %>%  

healthy_vs_unhealthy <- sst_all_data_clean %>% filter(condition %in% c("CorrectGo","FailedStop")) %>% dplyr::select(subid,waveid,trial_n,condition,reaction_time_clean,health_condition_label) %>%
  group_by(subid,waveid,condition,health_condition_label) %>%
  summarize(mean_rt=mean(reaction_time_clean,na.rm=TRUE), mean_trial_n=mean(trial_n)) %>%
  pivot_wider(names_from="health_condition_label",values_from = c("mean_rt","mean_trial_n"))

healthy_vs_unhealthy$healthy_minus_unhealthy_rt <- healthy_vs_unhealthy$mean_rt_Healthy-healthy_vs_unhealthy$mean_rt_Unhealthy
healthy_vs_unhealthy$healthy_minus_unhealthy_trial_n <- healthy_vs_unhealthy$mean_trial_n_Healthy-healthy_vs_unhealthy$mean_trial_n_Unhealthy
  


```


Does RT to healthy food decrease?

```{r}

t.test(results_by_subj$mean_rt_Healthy_w2_m_w1)
t.test(results_by_subj_brt$mean_rt_Healthy_w2_m_w1)
```
Yes! For the groups in genral, and for BRT in particular.

And does Unhealthy _increase_?

```{r}


t.test(results_by_subj$mean_rt_Unhealthy_w2_m_w1)

t.test(results_by_subj_brt$mean_rt_Unhealthy_w2_m_w1)
```

No--there's no change for this; in fact, for Unhealthy food, the directionality is in the REVERSE, but non-significantly so.


It's still going to be important to look at stop signal reaction time. We should be able to calculate that separately for healthy and unhealthy trials.

OK, these results are starting to look at bit odd. But I think--because of the titration design--these are most interpretable when comparing between healthy and unhealthy ratehr htan looking at unhealthy conditions by themselves.

Finally, before we look at conditions, let's summarize all the relevant effects so that we're very clear about the directionality.

```{r}
summary(results_by_subj_int$mean_rt_Healthy_waveid_1)
summary(results_by_subj_int$mean_rt_Healthy_waveid_2)
summary(results_by_subj_int$mean_rt_Unhealthy_waveid_1)
summary(results_by_subj_int$mean_rt_Unhealthy_waveid_2)
```

Specifically for the BRT group:

```{r}
rt_comparison<-results_by_subj_int %>% filter(condition=="CorrectGo") %>%
  select(subid,intervention_group,
        mean_rt_Healthy_waveid_1,
         mean_rt_Healthy_waveid_2,
         mean_rt_Unhealthy_waveid_1,
         mean_rt_Unhealthy_waveid_2) %>% 
  pivot_longer(cols = c(-subid,-intervention_group), names_to = "measure",values_to = "value") %>%
  mutate(measure = sub("waveid_","",sub("mean_rt_","",measure))) %>%
  separate(measure,into=c("Health","Wave"))
  
  
  
ggplot(rt_comparison,aes(x=intervention_group,group=interaction(Health, Wave,intervention_group),color=Wave,y=value))+geom_boxplot()+facet_wrap(~Health)

summary(results_by_subj_brt$mean_rt_Healthy_waveid_1)
summary(results_by_subj_brt$mean_rt_Healthy_waveid_2)
summary(results_by_subj_brt$mean_rt_Unhealthy_waveid_1)
summary(results_by_subj_brt$mean_rt_Unhealthy_waveid_2)
```

### Condition interactions
So--let's now look at change by condition! Because we're laser focused on BRT, just look at that.

```{r}


summary(lm(mean_rt_Healthy_m_Unhealthy_w2_m_w1~1,results_by_subj_brt_test %>% filter(condition=="CorrectGo")))
summary(lm(mean_rt_Healthy_m_Unhealthy_w2_m_w1~intervention_group,results_by_subj_brt_test %>% filter(condition=="CorrectGo")))


```

```{r}
summary(lm(mean_rt_Healthy_w2_m_w1~1,results_by_subj_brt_test %>% filter(condition=="CorrectGo")))

summary(lm(mean_rt_Healthy_w2_m_w1~intervention_group,results_by_subj_brt_test %>% filter(condition=="CorrectGo")))

summary(lm(mean_rt_Unhealthy_w2_m_w1~intervention_group,results_by_subj_brt_test %>% filter(condition=="CorrectGo")))

summary(lm(mean_rt_Healthy_m_Unhealthy_w2_m_w1~intervention_group,results_by_subj_brt_test %>% filter(condition=="CorrectGo")))
```
Let's look closer 

OK--the behavioral intervention worked--in the first month. It seems to fade away after that. What's odd is that a main effect takes over. Perhaps that reflects that the SST itself is practice?

The pre-registered hypothesis is **confirmed**. Behavioral response training (McKenzie) created a difference in unhealthy and healthy food training.




There are further questions:

 - Although we didn't pre-register it, we could look at how this might have affected Stop Signal Reaction Time (but what would be the implications?)
 - How does participants' RT change relate to neural activities? do participants with greater change have greater neural activity changes?
 - Does the magnitude of the change relate to FFQ changes?
 


### Controls

```{r}

results_by_subj_fulldataset <- results_by_subj %>% merge(full_dataset_with_groups,by.x = "subid",by.y="SID",all = FALSE)
summary(lm(healthy_minus_unhealthy_rt_w2_m_w1~intervention_group+age365_z+birthsex_factor,results_by_subj_fulldataset %>% filter(condition=="CorrectGo")))
```

### MLM

#### Contrasting health groups

```{r}
library(lme4)
grt_scores <- sst_all_data_clean %>% dplyr::select(subid,waveid,runid, trial_n,condition,reaction_time_clean,health_condition_label,SSD_recorded) %>%
  #get_more_sst_trial_stats() %>% 
  group_by(subid,waveid,runid) %>% 
  summarize(grt = mean(reaction_time_clean[condition=="CorrectGo"],na.rm=TRUE),
            grt_healthy = mean(reaction_time_clean[condition=="CorrectGo" & health_condition_label=="Healthy"],na.rm=TRUE),
            grt_unhealthy = mean(reaction_time_clean[condition=="CorrectGo" & health_condition_label=="Unhealthy"],na.rm=TRUE),
            ) %>% 
  mutate(grt_healthy_m_unhealthy = grt_healthy - grt_unhealthy) %>%
  merge(full_dataset_with_groups %>% dplyr::select(SID,intervention_group,age365_z,birthsex_factor),by.x = "subid",by.y="SID",all = FALSE)
  # summarize(mean_rt=mean(reaction_time_clean), mean_trial_n=mean(trial_n)) %>%
  # pivot_wider(names_from="health_condition_label",values_from = c("mean_rt","mean_trial_n"))

grt_scores$PostIntervention<-grt_scores$waveid>1
grt_scores_brt_v_cont <- grt_scores %>% filter(intervention_group %in% c("willamette","mckenzie"))

base_model = lmer(
  grt_healthy_m_unhealthy~waveid+intervention_group + (1 |subid),
     grt_scores_brt_v_cont %>% filter(waveid %in% c(1,2)))
summary(base_model)
intervention_model = lmer(
  grt_healthy_m_unhealthy~waveid*intervention_group + (1 |subid),
  grt_scores_brt_v_cont %>% filter(waveid %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)

```
Great. So we can again confirm our BRT hypothesis with Go Response Time.

Now let's do another model which captures the extra value we have over time.
```{r}



base_model = lmer(
  grt_healthy_m_unhealthy~waveid+intervention_group + PostIntervention + (1 |subid),
     grt_scores_brt_v_cont)
summary(base_model)
intervention_model = lmer(
  grt_healthy_m_unhealthy~waveid+PostIntervention*intervention_group + (1 |subid),
  grt_scores_brt_v_cont)

summary(intervention_model)

anova(base_model,intervention_model)
```

#### Separated by health condition label
```{r}


base_model = lmer(
  grt_healthy~waveid+intervention_group + (1 |subid),
     grt_scores_brt_v_cont %>% filter(waveid %in% c(1,2)))
summary(base_model)
intervention_model = lmer(
  grt_healthy~waveid*intervention_group + (1 |subid),
  grt_scores_brt_v_cont %>% filter(waveid %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)


base_model = lmer(
  grt_unhealthy~waveid+intervention_group + (1 |subid),
     grt_scores_brt_v_cont %>% filter(waveid %in% c(1,2)))
summary(base_model)
intervention_model = lmer(
  grt_unhealthy~waveid*intervention_group + (1 |subid),
  grt_scores_brt_v_cont %>% filter(waveid %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)

```

#### Health groups combined
```{r}

base_model = lmer(
  grt~waveid+intervention_group + (1 |subid),
     grt_scores_brt_v_cont %>% filter(waveid %in% c(1,2)))
summary(base_model)
intervention_model = lmer(
  grt~waveid*intervention_group + (1 |subid),
  grt_scores_brt_v_cont %>% filter(waveid %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)
```


## Cognitive reappraisal theory:

cognitive appraisal may lead to similar changes or it might not. CR gives people the ability to habitually or momentarily regulate their response to unhealthy food, but it isn’t clear whether this is going to be applied in very short time periods in an implicit manner as in the SST. Therefore, this is a test of the pathway CR uses. Does the CR modify the initial response to food, or does it take an effortful or at least specific contextual frame for it to be active?
    - If CR elicits a broad and momentary change then we can expect an increase in RT to unhealthy food
        - To confirm this hypothesis we’d have to see a greater increase in RT for CR vs. control
    - If CR requires a level of effort or is context restricted we expect no change in RT to unhealthy food
        - To really confirm this hypothesis, we’d have to see BRT being significantly *more* effective than CR, i.e., A greater increase in RT to unhealthy foods for BRT vs. CR

Let's take a look at that comparison graph again:

```{r}

  
ggplot(rt_comparison,aes(x=intervention_group,group=interaction(Health, Wave,intervention_group),color=Wave,y=value))+geom_boxplot()+facet_wrap(~Health)
```

### MLM

#### Separated by health condition label

What we are mainly hypothesizing is an INCREASE in RT in response to Unhealthy food particularly. Do we see that?
```{r}

grt_scores_cr_v_cont <- grt_scores %>% filter(intervention_group %in% c("willamette","umpqua"))

base_model = lmer(
  grt_unhealthy~waveid+intervention_group + (1 |subid),
     grt_scores_cr_v_cont %>% filter(waveid %in% c(1,2)))
summary(base_model)
intervention_model = lmer(
  grt_unhealthy~waveid*intervention_group + (1 |subid),
  grt_scores_cr_v_cont %>% filter(waveid %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)

```

And that seems to be what happens: for CR, the reaction time for Unhealthy food in particular _increases_. **Hypothesis confirmed**

Let's look across all waves...

```{r}



base_model = lmer(
  grt_unhealthy~waveid+PostIntervention+intervention_group + (1 |subid),
     grt_scores_cr_v_cont)
summary(base_model)
intervention_model = lmer(
  grt_unhealthy~waveid+PostIntervention*intervention_group + (1 |subid),
  grt_scores_cr_v_cont)

summary(intervention_model)

anova(base_model,intervention_model)

```

Yes, it works across all waves, too.


While we're here, let's look at healthy food:

```{r}


base_model = lmer(
  grt_healthy~waveid+intervention_group + (1 |subid),
     grt_scores_cr_v_cont %>% filter(waveid %in% c(1,2)))
summary(base_model)
intervention_model = lmer(
  grt_healthy~waveid*intervention_group + (1 |subid),
  grt_scores_cr_v_cont %>% filter(waveid %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)
```

There's a change in healthy food in the _same_ direction--the CR group seem to be more mindful either way.


#### Contrasting health groups

It will be useful corroborating evidence if the increase for unhealthy food is _greater_ than the increase for healthy food

```{r}
library(lme4)


#flip around to make it easier to interpret.
grt_scores_cr_v_cont$grt_unhealthy_m_healthy <- - grt_scores_cr_v_cont$grt_healthy_m_unhealthy

base_model = lmer(
  grt_unhealthy_m_healthy~waveid+intervention_group + (1 |subid),
     grt_scores_cr_v_cont %>% filter(waveid %in% c(1,2)))
summary(base_model)
intervention_model = lmer(
  grt_unhealthy_m_healthy~waveid*intervention_group + (1 |subid),
  grt_scores_cr_v_cont %>% filter(waveid %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)

```


```{r}



base_model = lmer(
  grt_unhealthy_m_healthy~waveid+intervention_group + PostIntervention + (1 |subid),
     grt_scores_cr_v_cont)
summary(base_model)
intervention_model = lmer(
  grt_unhealthy_m_healthy~waveid+PostIntervention*intervention_group + (1 |subid),
  grt_scores_cr_v_cont)

summary(intervention_model)

anova(base_model,intervention_model)
```
So in this test we do see a INCREASE in reaction time to unhealthy (vs. healthy) food.   This is consistent with the hypothesis for cognitive appraisal.



#### Health groups combined
```{r}

base_model = lmer(
  grt~waveid+intervention_group + (1 |subid),
     grt_scores_cr_v_cont %>% filter(waveid %in% c(1,2)))
summary(base_model)
intervention_model = lmer(
  grt~waveid*intervention_group + (1 |subid),
  grt_scores_cr_v_cont %>% filter(waveid %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)
```



# Relationship with FFQ

Can we detect a relationship between change in the unhealthy food responses, and the FFQ scores?

```{r}



summary(lm(NUTRIENT_RICH_FOODS_INDEX_2wkAverage_FFQ_NutrientDensity_1wkpost~NUTRIENT_RICH_FOODS_INDEX_2wkAverage_FFQ_NutrientDensity_baseline+intervention_group,full_dataset_with_groups))
summary(lm(total_calorie_FFQ_NutrientDensity_1wkpost~total_calorie_FFQ_NutrientDensity_baseline+intervention_group,full_dataset_with_groups))
summary(lm(CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost~CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_baseline+intervention_group,full_dataset_with_groups))
#t.test(full_dataset_with_groups$total_calorie_FFQ_NutrientDensity_1wkpost-full_dataset_with_groups$total_calorie_FFQ_NutrientDensity_baseline)
```

Let's use the third of these--there's the strongest group effect to explain.

NB that there's no change at all for McKenzie (BRT)! So any effect would be hidden in the noise and only applicable for a subset of the subjects.

What I am curious about is: does the magnitude of change in RT predict the 1-week post score?

This is getting pretty messy and is at the point where I should be integrating the mixed effects model but let's run with it so far...

```{r}

full_dataset_with_groups$CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline<-(
  full_dataset_with_groups$CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost - 
    full_dataset_with_groups$CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_baseline)

full_dataset_with_groups$NUTRIENT_RICH_FOODS_INDEX_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline<-(
  full_dataset_with_groups$NUTRIENT_RICH_FOODS_INDEX_2wkAverage_FFQ_NutrientDensity_1wkpost - 
    full_dataset_with_groups$NUTRIENT_RICH_FOODS_INDEX_2wkAverage_FFQ_NutrientDensity_baseline)


summary(lm(CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline~intervention_group+age365_z,full_dataset_with_groups)) #yep, the detected difference persists


results_by_subj_int_ffq <- results_by_subj %>% merge(full_dataset_with_groups %>% dplyr::select(SID,intervention_group,CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline,NUTRIENT_RICH_FOODS_INDEX_2wkAverage_FFQ_NutrientDensity_baseline,NUTRIENT_RICH_FOODS_INDEX_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline,age365_z),by.x = "subid",by.y="SID",all = FALSE)

summary(lm(CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline~intervention_group+age365_z,results_by_subj_int_ffq %>% filter(condition=="CorrectGo"   & intervention_group!="umpqua")))

summary(lm(CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline~healthy_minus_unhealthy_rt_w2_m_w1+age365_z,results_by_subj_int_ffq %>% filter(condition=="CorrectGo"   & intervention_group!="umpqua")))
summary(lm(CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline~healthy_minus_unhealthy_rt_w2_m_w1*intervention_group+age365_z*intervention_group,results_by_subj_int_ffq %>% filter(condition=="CorrectGo" & intervention_group!="umpqua")))


```


## Mediation analysis

```{r}

results_by_subj_fulldataset$CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline <- (
    results_by_subj_fulldataset$CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost - results_by_subj_fulldataset$CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_baseline
) 
selected_rows <- (
  (results_by_subj_fulldataset$intervention_group!="umpqua") & 
  (results_by_subj_fulldataset$condition=="CorrectGo") & 
  (
    rowSums(is.na(results_by_subj_fulldataset %>% dplyr::select(intervention_group,healthy_minus_unhealthy_rt_w2_m_w1,CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline))
    )==0
  ))
#mediation test: intervention_group ->healthy_minus_unhealthy_rt_w2_m_w1 -> CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost
med_model<-lm(healthy_minus_unhealthy_rt_w2_m_w1~intervention_group+age365_z+birthsex_factor,results_by_subj_fulldataset[selected_rows,])
summary(med_model)
ffq_model_medonly <- lm(CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline~healthy_minus_unhealthy_rt_w2_m_w1+age365_z+birthsex_factor,results_by_subj_fulldataset[selected_rows,])
summary(ffq_model_medonly)
ffq_model <- lm(CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline~healthy_minus_unhealthy_rt_w2_m_w1+intervention_group+age365_z+birthsex_factor,results_by_subj_fulldataset[selected_rows,])
summary(ffq_model)

direct_model <- lm(CncrPrevLessCncrProm_2wkAverage_FFQ_NutrientDensity_1wkpost_m_baseline~intervention_group+age365_z+birthsex_factor,results_by_subj_fulldataset[selected_rows,])
summary(direct_model)

set.seed(08311356)
summary(mediation::mediate(model.m=med_model,model.y=ffq_model,mediator="healthy_minus_unhealthy_rt_w2_m_w1",treat="intervention_group",robustSE = TRUE,boot = TRUE,sims = 2000))
```


# Reprise with SSRT

We didn't actually include a hypothesis with respect to SSRT! But, nevertheless, we can take a look.

I think we are best off using `mean_ssrt_0` or `SSRTQuant`, which does not assume 50% breakeven but instead looks at 

```{r}
source("SST_processing.R")


ssrt_scores <- sst_all_data_clean %>% dplyr::select(subid,waveid,runid, trial_n,condition,reaction_time_clean,health_condition_label,SSD_recorded) %>%
  group_by(subid,waveid,runid,health_condition_label) %>%
  get_sst_for_run() %>% 
  group_by(subid,waveid,runid,health_condition_label) %>% 
  summarize(mean_ssrt_0 = mean(ssrt_0[condition=="CorrectStop"],na.rm=TRUE)) %>%
  filter(!is.na(health_condition_label))

ssrt_scores_by_wave <- ssrt_scores %>%
  pivot_wider(id_cols=c("subid","waveid","runid"),names_from="health_condition_label",values_from="mean_ssrt_0",names_prefix = "mean_ssrt_0_") %>%
  rename("session_id"="waveid","SID"="subid") %>% ungroup() %>% dplyr::select(-runid)

ssrt_scores_by_wave$Healthy_m_Unhealthy_mean_ssrt_0<-ssrt_scores_by_wave$mean_ssrt_0_Healthy-ssrt_scores_by_wave$mean_ssrt_0_Unhealthy
ssrt_scores_by_wave$Unhealthy_m_Healthy_mean_ssrt_0<-ssrt_scores_by_wave$Healthy_m_Unhealthy_mean_ssrt_0
all_ppt_data_w_ssrt<-merge(all_ppt_data,ssrt_scores_by_wave,by=c("SID","session_id"),all.x=TRUE,all.y=TRUE)

all_ppt_data_w_ssrt_int <- all_ppt_data_w_ssrt %>% merge(full_dataset_with_groups %>% dplyr::select(SID,intervention_group,age365_z,birthsex_factor),by.x = "SID",by.y="SID",all = FALSE) %>%
  filter(session_id>0)

value_cols <- colnames(all_ppt_data_w_ssrt)[!(colnames(all_ppt_data_w_ssrt) %in% c("SID","session_id")) & sapply(all_ppt_data_w_ssrt,class) %in% c("numeric","logical")]
all_ppt_data_subjs <- all_ppt_data_w_ssrt %>% dplyr::select(SID)

all_ppt_data_w_ssrt_int$PostIntervention <- all_ppt_data_w_ssrt_int$session_id>1
all_ppt_data_w_ssrt_brt<-all_ppt_data_w_ssrt_int %>% filter(intervention_group %in% c("willamette","mckenzie"))
all_ppt_data_w_ssrt_cr <- all_ppt_data_w_ssrt_int %>% filter(intervention_group %in% c("willamette","umpqua"))
```

## Exploratory

When calculating SSRTs for subsets, should we just use the Go trials relating to that subset? Probably, yeah, because if there really is separate SSRT, those are the informative ones we need to know.
```{r}

ssrt_comparison<-all_ppt_data_w_ssrt_int %>%
  filter(session_id>0 & session_id %in% c(1,2))%>%
  select(SID,session_id,intervention_group,
        mean_ssrt_0_Healthy,
         mean_ssrt_0_Unhealthy) %>% 
  pivot_longer(cols = c(-SID,-session_id,-intervention_group), names_to = "measure",values_to = "ssrt_mean") %>%
  mutate(Health = sub("mean_ssrt_0_","",measure),
         session_id = as.factor(session_id)
         )
  
  
  
ggplot(ssrt_comparison,aes(x=intervention_group,group=interaction(Health, session_id,intervention_group),color=session_id,y=ssrt_mean))+geom_boxplot()+facet_wrap(~Health)

```


```{r}
# 
# all_ppt_data_s1 <- all_ppt_data_w_ssrt %>% dplyr::filter(session_id==1) %>% merge(all_ppt_data_subjs,by="SID", sort = TRUE,all.y=TRUE)
# all_ppt_data_s2 <- all_ppt_data_w_ssrt %>% dplyr::filter(session_id==2) %>% merge(all_ppt_data_subjs,by="SID", sort = TRUE,all.y=TRUE)
# all_ppt_data_s3 <- all_ppt_data_w_ssrt %>% dplyr::filter(session_id==3) %>% merge(all_ppt_data_subjs,by="SID", sort = TRUE,all.y=TRUE)
# 
# #w2vw1
# change_scores <- (all_ppt_data_s2[,value_cols] - all_ppt_data_s1[,value_cols])
# colnames(change_scores)<-paste0(colnames(change_scores),"_s2_m_s1")
# all_ppt_data_s2_m_s1 <- cbind(all_ppt_data_subjs,change_scores)
# all_ppt_data_s2_m_s1_int<-all_ppt_data_s2_m_s1 %>% merge(full_dataset_with_groups %>% dplyr::select(SID,intervention_group,age365_z,birthsex_factor),by.x = "SID",by.y="SID",all = FALSE)
# 
# 
# #w3vw1
# change_scores <- (all_ppt_data_s3[,value_cols] - all_ppt_data_s1[,value_cols])
# colnames(change_scores)<-paste0(colnames(change_scores),"_s3_m_s1")
# all_ppt_data_s3_m_s1 <- cbind(all_ppt_data_subjs,change_scores)
# all_ppt_data_s3_m_s1_int<-all_ppt_data_s3_m_s1 %>% merge(full_dataset_with_groups %>% dplyr::select(SID,intervention_group,age365_z,birthsex_factor),by.x = "SID",by.y="SID",all = FALSE)
# 




```

## BRT theory:

  - BR Training (vs. Control) decreases RT to healthy foods and increases RT to unhealthy foods; because the SST task calibrates difficulty level, these two will be difficult to pull apart but we should see an increase Unhealthy RT - Healthy RT
    

### Contrasting health groups
```{r}
library(lme4)

base_model = lmer(Unhealthy_m_Healthy_mean_ssrt_0~session_id+intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_brt %>% filter(session_id %in% c(1,2)))
summary(base_model)
intervention_model = lmer(Unhealthy_m_Healthy_mean_ssrt_0~session_id*intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_brt %>% filter(session_id %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)

```


```{r}
library(lme4)

base_model = lmer(Unhealthy_m_Healthy_mean_ssrt_0~session_id+PostIntervention+intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_brt)
summary(base_model)
intervention_model = lmer(Unhealthy_m_Healthy_mean_ssrt_0~session_id+PostIntervention*intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_brt)

summary(intervention_model)

anova(base_model,intervention_model)

```


### Separated by health condition label

Just look at unhealthy...


```{r}
library(lme4)

base_model = lmer(mean_ssrt_0_Unhealthy~session_id+intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_brt %>% filter(session_id %in% c(1,2)))
intervention_model = lmer(mean_ssrt_0_Unhealthy~session_id*intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_brt %>% filter(session_id %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)

```

## CRT

### Separated by health condition label

Just look at unhealthy...

What we are mainly hypothesizing is an INCREASE in RT in response to Unhealthy food particularly. Do we see that?

```{r}
library(lme4)

base_model = lmer(mean_ssrt_0_Unhealthy~session_id+intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_cr %>% filter(session_id %in% c(1,2)))
intervention_model = lmer(mean_ssrt_0_Unhealthy~session_id*intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_cr %>% filter(session_id %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)

```




### Contrasting health groups
```{r}
library(lme4)

base_model = lmer(Unhealthy_m_Healthy_mean_ssrt_0~session_id+intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_cr %>% filter(session_id %in% c(1,2)))
summary(base_model)
intervention_model = lmer(Unhealthy_m_Healthy_mean_ssrt_0~session_id*intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_cr %>% filter(session_id %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)

```


```{r}
library(lme4)

base_model = lmer(Unhealthy_m_Healthy_mean_ssrt_0~session_id+PostIntervention+intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_brt)
summary(base_model)
intervention_model = lmer(Unhealthy_m_Healthy_mean_ssrt_0~session_id+PostIntervention*intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_brt)

summary(intervention_model)

anova(base_model,intervention_model)

```


## Post-hoc, health groups combined

```{r}


ssrt_scores <- sst_all_data_clean %>% dplyr::select(subid,waveid,runid, trial_n,condition,reaction_time_clean,SSD_recorded) %>%
  group_by(subid,waveid,runid) %>%
  get_sst_for_run() %>% 
  group_by(subid,waveid,runid) %>% 
  summarize(mean_ssrt_0 = mean(ssrt_0[condition=="CorrectStop"],na.rm=TRUE))

ssrt_scores_by_wave <- ssrt_scores %>%
#   pivot_wider(id_cols=c("subid","waveid","runid"),names_from="health_condition_label",values_from="mean_ssrt_0",names_prefix = "mean_ssrt_0_") %>%
   rename("session_id"="waveid","SID"="subid") %>% ungroup() %>% 
  dplyr::select(-runid)



all_ppt_data_w_ssrt<-merge(all_ppt_data,ssrt_scores_by_wave,by=c("SID","session_id"),all.x=TRUE,all.y=TRUE)

#general SSRT, without dividing by health condition, 
all_ppt_data_w_ssrt_int_gen <- all_ppt_data_w_ssrt %>% merge(full_dataset_with_groups %>% dplyr::select(SID,intervention_group,age365_z,birthsex_factor),by.x = "SID",by.y="SID",all = FALSE) %>%
  filter(session_id>0)

all_ppt_data_w_ssrt_int_gen$PostIntervention <- all_ppt_data_w_ssrt_int_gen$session_id>1

```

```{r}

base_model = lmer(mean_ssrt_0~session_id+intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_int_gen %>% filter(session_id %in% c(1,2)))
intervention_model = lmer(mean_ssrt_0~session_id*intervention_group + age365_z + (1 |SID),
     all_ppt_data_w_ssrt_int_gen %>% filter(session_id %in% c(1,2)))

summary(intervention_model)

anova(base_model,intervention_model)
```

# Proportions of healthy vs unhealthy stop trials that are correct

```{r}
sst_all_data_by_run <- sst_all_data_clean %>% mutate(HasResponse=reaction_time>0) %>% group_by(subid, waveid, runid,go_no_go_condition_label,health_condition_label) %>%
  summarise(PropTrialsWithResponse = sum(HasResponse)/n()) %>% pivot_wider(names_from="health_condition_label",values_from = "PropTrialsWithResponse",names_prefix = "P_Responses") %>%
  mutate(P_Responses_healthy_minus_unhealthy = P_ResponsesHealthy - P_ResponsesUnhealthy)

value_cols <- colnames(sst_all_data_by_run)[!(colnames(sst_all_data_by_run) %in% c("subid", "waveid","condition", "runid","go_no_go_condition_label"))]

sst_all_data_by_subj <- sst_all_data_by_run %>% ungroup() %>% select(-runid) %>% pivot_wider(names_from = waveid,id_cols = c(subid,go_no_go_condition_label),values_from = value_cols) %>%
  mutate(P_Responses_healthy_minus_unhealthy_2_m_1 = (P_Responses_healthy_minus_unhealthy_2-P_Responses_healthy_minus_unhealthy_1)
  ) %>% merge(full_dataset_with_groups %>% dplyr::select(SID,intervention_group),by.x = "subid",by.y="SID",all = FALSE)
sst_all_data_by_subj$intervention_group<-factor(sst_all_data_by_subj$intervention_group,levels=c("willamette","umpqua","mckenzie"))

summary(lm(P_Responses_healthy_minus_unhealthy_2_m_1~1,sst_all_data_by_subj %>% filter(go_no_go_condition_label=="Stop")))
summary(lm(P_Responses_healthy_minus_unhealthy_2_m_1~intervention_group,sst_all_data_by_subj %>% filter(go_no_go_condition_label=="Stop")))
```


# Summary

 - Reaction time for healthy (vs unhealthy) correct go trials is significantly shorter in post-intervention compared to pre-intervention
 
```{r}
#should do a line graph with CIs of healthy, unhealthy, post, pre

healthy_unhealthy_rt_graph <- results_by_subj_int %>% 
  pivot_longer(cols=tidyr::contains("ealthy"),names_to = "measure",values_to = "value") %>%
  dplyr::filter(condition=="CorrectGo" & grepl("waveid",measure)) %>%
  separate(col=measure,into=c("measure","wave"),sep="_waveid_") %>%
  mutate(wave=as.numeric(wave))

healthy_unhealthy_rt_graph_means<-healthy_unhealthy_rt_graph %>%
  group_by(condition,intervention_group,measure,wave) %>%
  summarize(subj_mean=mean(value,na.rm=TRUE))


  
library(ggplot2)
# ggplot(healthy_unhealthy_rt_graph_means %>% filter(wave<3 & measure=="healthy_minus_unhealthy_rt" & condition=="CorrectGo"),
#        aes(linetype=measure,x=wave,y=subj_mean,color=intervention_group,group=intervention_group))+geom_line()+geom_point()+labs(
#          y="healthy_minus_unhealthy_rt",
#          title="SST healthy minus unhealthy stimulus Go Reaction Times by intervention group")
# 

ggplot(healthy_unhealthy_rt_graph %>% filter(wave<3 & measure=="healthy_minus_unhealthy_rt" & condition=="CorrectGo"),
       aes(linetype=measure,x=wave,y=value,group=interaction(wave)))+geom_boxplot()+labs(
         y="healthy_minus_unhealthy_rt",
         title="SST healthy minus unhealthy stimulus Go Reaction Times by intervention group")


# ggplot(healthy_unhealthy_rt_graph_means %>% filter(wave<3 & measure %in% c("mean_rt_Healthy","mean_rt_Unhealthy") & condition=="CorrectGo"),
#        aes(linetype=measure,x=wave,y=subj_mean,fill=intervention_group,group=interaction(intervention_group,measure)))+geom_bar(stat = "identity",position = "dodge") + facet_grid(~measure)

```
 
 
 - There is greater change in the BRT group vs. the control group
 
```{r}
library(forcats)

relabel_intervention_groups <- function(igcol){
  trans_col<-fct_recode(igcol,  "control"="willamette", "BRT"="mckenzie", "CR"="umpqua")
  return(trans_col)
}

healthy_unhealthy_rt_graph$intervention_group_labelled <-relabel_intervention_groups(healthy_unhealthy_rt_graph$intervention_group)

ggplot(healthy_unhealthy_rt_graph %>% filter(wave<3 & measure=="healthy_minus_unhealthy_rt" & condition=="CorrectGo"),
       aes(linetype=measure,x=wave,y=value,fill=intervention_group_labelled,group=interaction(wave,intervention_group_labelled)))+geom_boxplot()+labs(
         y="Health > Unhealthy RT",
         title="SST healthy minus unhealthy stimulus Go Reaction Times by intervention group")

```
 
```{r}

results_by_subj_int$intervention_group_labelled <- relabel_intervention_groups(results_by_subj_int$intervention_group)
ggplot(results_by_subj_int %>% filter(condition=="CorrectGo"),
       aes(x=intervention_group_labelled,y=healthy_minus_unhealthy_rt_w2_m_w1,
           fill=intervention_group_labelled,group=interaction(intervention_group_labelled))
       )+geom_boxplot()+labs(
         y="Healthy > Unhealthy RT W2>W1",
         title="SST healthy minus unhealthy stimulus Go Reaction Times by intervention group")


```
 
 
 - This suggests BRT is narrowly effective within its domain, and the efficacy of cognitive reappraisal is also confirmed even when participants are not deliberately trying to do reappraisal
 - Age and birth sex don't change these effects
 - The effect of the intervention on FFQ is mediated by healthy-unhealthy Go response time (p<0.05), but this may not be a convincing result because the effect of the BRT intervention on FFQ, considered in isolation, is not quite significant