---
title: "SSI Megastudy Analyses"
author: "Benji Kaveladze"
date: "2025-07-26"
output:
  pdf_document: default
  html_document:
    df_print: default
  word_document: default
---

These are the analyses for the megastudy of depression SSIs. For the pre-registered analysis plan, see <https://cdn.clinicaltrials.gov/large-docs/68/NCT06856668/Prot_SAP_ICF_001.pdf>. For the pre-registered analysis scripts, see <https://osf.io/agvh6/files/osfstorage#>.

Reach out to me at benjamin.kaveladze\@northwestern.edu with any questions.

```{r setup, include=FALSE}
#CHANGE BELOW LINE TO YOUR WORKING DIRECTORY!
knitr::opts_knit$set(
  digits=3,
  root.dir = normalizePath('~/megastudy stuff/megastudy_real_data/')
  )
#if you want a pdf of the output, run tinytex::install_tinytex() on your console before you run this code
```

```{r,warning=F,message=F}
#### STEP 1: LOAD DATA ####
#loading all the packages. 
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  lme4,sjPlot,tidyverse,effectsize,ggdist,gridExtra,ggcorrplot,
  broom.mixed,gtsummary,flextable,gt,emmeans,readxl,performance,
  rstatix,lmerTest,randomForest, car,psych,viridis,patchwork,
  kableExtra,plotrix) 

theme_set(theme_ggdist())

#read in data 
load("utw_selectvars.Rdata") #or utw_allvars if using the full data
mega_alltimes_complete <- mega_selectvars #or <- mega_allvars if using the full data

# Keep all rows from non-duplicate participants; For duplicates (based on matching
# IP/gender/race/age at baseline), keep only the first entry per timepoint.
# Then exclude any participants who never completed baseline.

mega_alltimes_preparation1 <- mega_alltimes_complete %>%
  group_by(IP_address_id, time) %>%
  arrange(StartDate, .by_group = TRUE) %>%
  filter( ip_gender_race_age_same_at_bl==0 |
      (ip_gender_race_age_same_at_bl==1 & row_number() == 1 )) %>%
  ungroup()  
  

# above leaves some rows from people who completed wk4 on a different device,
# so here dropping people who didn't complete the baseline
mega_alltimes <- mega_alltimes_preparation1 %>% filter(time == "Baseline") %>%  
    distinct(pid) %>%
    inner_join(mega_alltimes_preparation1, by = "pid")  

#write_csv(mega_alltimes,'utw_selectvars_analytic_sample_(n=7505).csv')

#function for calculating Cohen's D in a longitudinal mixed model
#based on Feingold, A. (2009). https://doi.org/10.1037/a0014699 thanks to Mike Russell for explaining
mixed_cohens_d <- function(model){
  #outcome sd at baseline: 
  # u0i: person mean of outcome minus grand mean of outcome for person i (between residual; or the random intercept)
  #eit : current outcome for person i at time t minus their person mean (within residual)
  outcome_sd <- sqrt(summary(model)$varcor$pid[1] + 
                       summary(model)$sigma^2)
  (tidy_model <- broom.mixed::tidy(model,conf.int=T) %>% 
      mutate(cohens_d= estimate/outcome_sd,
             cohens_d_2.5= conf.low/outcome_sd,
             cohens_d_97.5= conf.high/outcome_sd) %>% 
      filter(effect=='fixed',term!='(Intercept)') %>% 
      select(-c(group,effect,std.error,conf.low,conf.high)) %>%
      mutate(across(where(is.numeric),~round(.,digits=2))) %>% 
      relocate(cohens_d,.after = 'term'))
}

add_cohens_d <- function(model, caption = NULL) {
  tidy(model) %>% 
    mutate(p.value = ifelse(p.value < 0.001,
                            formatC(p.value, format = "e", digits = 2),
                            sprintf("%.3f", p.value))) %>%
    left_join(
      mixed_cohens_d(model) %>% 
        select(term, starts_with('cohens_d')),
      by = 'term'
    ) %>% 
    relocate(c(effect, group), .after = last_col()) %>% 
    select(-std.error, -effect, -group,-df,-statistic) %>% 
    kable(digits = 3, caption = caption)
}

```

------------------------------------------------------------------------

------------------------------------------------------------------------


## Preliminary Analyses for the Primary Analysis

##### Depressive symptoms score across conditions at baseline

```{r}
mega_alltimes %>% filter(time=='Baseline') %>% anova_test(phq9_score~group)
```

##### Depressive symptoms score baseline to wk4 in passive control

```{r}
just_time_phq9 <- lmer(
  phq9_score ~ time + (1|pid), mega_alltimes %>% 
    filter(group=='Passive Control'))
add_cohens_d(just_time_phq9,'Depressive symptoms over time, Control Only')
just_time_phq9 %>% tidy %>% select(term,statistic,df) %>% kable(digits = 3)
icc(just_time_phq9)
```

##### Depressive symptoms score correlation across timepoints

```{r}
mega_alltimes %>% filter(group=='Passive Control') %>% 
  select(pid,phq9_score,time) %>% 
  pivot_wider(names_from='time',values_from = 'phq9_score') %>% 
  filter(!is.na(`Week 4`)) %>% 
  cor_test(Baseline,`Week 4`)
#df = n (1504) -2
```

##### Differential attrition during the SSI or by follow-up
```{r}
#attrited before post-test measures
finished_post_test<- mega_alltimes$pid[!is.na(mega_alltimes$`Q559_Page Submit`)]
mega_alltimes <- mega_alltimes %>% 
  mutate(attrited_b4_post_test= ifelse(pid %in% finished_post_test,0,1))

#attrited before follow-up measures
finished_followup<- mega_alltimes$pid[!is.na(mega_alltimes$`Q423_Page Submit`)]
mega_alltimes <- mega_alltimes %>% 
  mutate(attrited_b4_followup= ifelse(pid %in% finished_followup,0,1))

mega_alltimes %>% filter(time=='Baseline') %>% 
  count(attrited_b4_post_test) %>% mutate(prop=n/sum(n))
mega_alltimes %>% filter(time=='Baseline') %>% 
  count(attrited_b4_followup) %>% mutate(prop=n/sum(n))

paste('We dropped data from ',
      1-(round(mega_alltimes %>% filter(time=='Baseline') %>%nrow()/
  mega_alltimes_complete %>% filter(time=='Baseline') %>% nrow(),2)*100),
  '% of the sample due to responses coming from the same IP address, gender, age, and race at baseline.')

paste('The sample used for analyses included',
      mega_alltimes %>% filter(time=='Baseline') %>% nrow(),
      'participants at baseline.')

paste(mega_alltimes %>% filter(time=='Baseline',!is.na(`Q1525_Page Submit`)) %>% nrow(),
      'Participants finished all post-test measures.')

paste(mega_alltimes %>% filter(time=='Week 4',Progress==100) %>% nrow(),
      'Participants finished the Week Four follow-up.',
      (round(100*(mega_alltimes %>% filter(time=='Week 4',Progress==100) %>% nrow())/
       (mega_alltimes %>% filter(time=='Baseline') %>% nrow()),2)),"% of those who were randomized")

# Test 1 - Any Differential Attrition Rate?

# after completing post-test measures?
mega_alltimes %>%
  filter(time=='Baseline') %>% 
  group_by(group) %>%
  summarise(
    Attrition_Rate_before_posttest  = mean(attrited_b4_post_test) * 100) %>%
  arrange(desc(Attrition_Rate_before_posttest))

chisq.test(with(mega_alltimes %>% filter(time=='Baseline'),
                 table(attrited_b4_post_test,group)))

# before follow-up
mega_alltimes %>%
  filter(time=='Baseline') %>% 
  group_by(group) %>%
  summarise(
    Attrition_Rate_before_Followup  = mean(attrited_b4_followup) * 100) %>%
  arrange(desc(Attrition_Rate_before_Followup))

chisq.test(with(mega_alltimes %>% filter(time=='Baseline'),
                 table(attrited_b4_followup,group)))
```

##### Duration of survey and attention checks
```{r}
mega_alltimes %>% filter(time=='Baseline')%>% 
  count(att_check==4) %>% mutate(prop=n/sum(n))

mega_alltimes %>% filter(time=='Baseline') %>% 
  count(`Duration (in seconds)`<300)%>% mutate(prop=n/sum(n))
```

--- 

# Primary Analysis (Pre-registered)

-   A mixed-effects linear regression model to test if depressive symptoms change differently from baseline to four-week follow-up in each of the 12 active conditions compared to the passive control condition. Timepoint is a two-level factor (baseline and week 4 follow up). The model includes a participant identifier as a random intercept.
-   Note that in this analysis each condition is compared to the **control condition**. That also means you can't use the error bars in the plot to compare SSIs to each other.


------------------------------------------------------------------------


```{r,warning=F,message=F}
m1_phq9 <- lmer(phq9_score ~ time * group + (1|pid),
                data=mega_alltimes)

phq9_cohens_d <- mixed_cohens_d(m1_phq9)
add_cohens_d(m1_phq9,caption='Predicting: change in depressive symptoms over four weeks')
m1_phq9 %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')

#tab_model(m1_phq9, p.val = "satterthwaite", show.df = TRUE, dv.labels ='Depression Mixed Effects Model')

#plot showing Change in Depression over Four Weeks Compared to a passive control
wk4_int_terms <- m1_phq9 %>% tidy() %>% filter(str_detect(term,'Week 4:')) %>% pull(term)

plot_data <- plot_model(m1_phq9, sort.est = T, terms = wk4_int_terms)$data
plot_data$color_override <- ifelse((plot_data$conf.low <= 0 & plot_data$conf.high >= 0), "darkgrey", NA)

(phq9_plot <- plot_model(
  m1_phq9, sort.est = T, terms = wk4_int_terms, auto_label = F,
  show.values = T, value.offset = .4, value.size = 3.2) +
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    scale_color_sjplot("circus", reverse = T) +
    # Override colors for overlapping confidence intervals
    geom_point(data = subset(plot_data, !is.na(color_override)), 
               aes(x = term, y = estimate), 
               color = "darkgrey", size = 2.5) +
    geom_errorbar(data = subset(plot_data, !is.na(color_override)), 
                  aes(x = term, y = estimate, 
                      ymin = conf.low, ymax = conf.high), 
                  color = "darkgrey", width = 0) +
    ggtitle('Change in Depressive Symptoms',
            'Baseline to Week Four, PHQ-9 (range 0-27)'))
```

------------------------------------------------------------------------

------------------------------------------------------------------------

# Secondary Analyses (Pre-registered)

### *H1*: The experimental SSIs will differ from the active comparison SSI

-   Redo the main analysis, except using the active comparison condition as the comparator.

```{r,message=F,warning=F}
#make active comparison SSI (abc project) the comparator
mega_alltimes_active_first <- mega_alltimes %>% 
  mutate(group=fct_relevel(group,'Behavioral Activation',after=0))

m1_phq9_activefirst <- lmer(phq9_score ~ time * group + (1|pid),
                            data=mega_alltimes_active_first)

add_cohens_d(m1_phq9_activefirst,'Predicting: change in depressive symptoms')

wk4_int_terms_activefirst <- m1_phq9_activefirst %>% tidy() %>% filter(str_detect(term,'Week 4:')) %>% pull(term)
plot_model(m1_phq9_activefirst,sort.est = T, terms=wk4_int_terms_activefirst,
           auto_label=F,show.values = T, value.offset =.4,value.size = 3.2)+
      scale_color_sjplot("circus",reverse=T) +
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x))+
  ggtitle('Change in depression by condition', 
          'Baseline to Week 4, Vs. the behavioral activation SSI')
```

------------------------------------------------------------------------

### *H2*: Secondary well-being outcomes will differ between each SSI and the passive control

-   Replicate the primary analysis for each secondary outcome. For outcomes measured at three timepoints, compare baseline to each timepoint.

#### Agency

```{r,warning=F,message=F}
m1_agency <- lmer(agency_score ~ time * group + (1|pid),data=mega_alltimes)

agency_cohens_d <- mixed_cohens_d(m1_agency)
add_cohens_d(m1_agency, 'Predicting: change in agency')
m1_agency %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')


plot_data <- plot_model(m1_agency, sort.est = T, terms = wk4_int_terms)$data
plot_data$color_override <- ifelse((plot_data$conf.low <= 0 & plot_data$conf.high >= 0), "darkgrey", NA)

(agency_wk4_plot <- plot_model(
  m1_agency,sort.est = T,terms=wk4_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_color_sjplot("circus",reverse=T) +
    # Override colors for overlapping confidence intervals
    geom_point(data = subset(plot_data, !is.na(color_override)), 
               aes(x = term, y = estimate), 
               color = "darkgrey", size = 2.5) +
    geom_errorbar(data = subset(plot_data, !is.na(color_override)), 
                  aes(x = term, y = estimate, 
                      ymin = conf.low, ymax = conf.high), 
                  color = "darkgrey", width = 0) +
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Change in Agency',
            'Baseline to Week 4, Hope Pathways (range 3-24)'))

posttest_int_terms <- m1_agency %>% tidy() %>% filter(str_detect(term,'Post-test:')) %>% pull(term)

plot_data <- plot_model(m1_agency, sort.est = T, terms = posttest_int_terms)$data
plot_data$color_override <- ifelse((plot_data$conf.low <= 0 & plot_data$conf.high >= 0), "darkgrey", NA)

(agency_posttest_plot <- plot_model(
  m1_agency,sort.est = T,terms=posttest_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_color_sjplot("circus",reverse=T) +
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
        # Override colors for overlapping confidence intervals
    geom_point(data = subset(plot_data, !is.na(color_override)), 
               aes(x = term, y = estimate), 
               color = "darkgrey", size = 2.5) +
    geom_errorbar(data = subset(plot_data, !is.na(color_override)), 
                  aes(x = term, y = estimate, 
                      ymin = conf.low, ymax = conf.high), 
                  color = "darkgrey", width = 0) +
    ggtitle('Change in Agency',
            'Baseline to Post-test, Hope Pathways (range 3-24)'))

```

#### Hopelessness

```{r,warning=F,message=F}
m1_hopeless <- lmer(hopeless_score ~ time * group + (1|pid),data=mega_alltimes)

hopeless_cohens_d <- mixed_cohens_d(m1_hopeless)
add_cohens_d(m1_hopeless,'Predicting: change in hopelessness')
m1_hopeless %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')

plot_data <- plot_model(m1_hopeless, sort.est = T, terms = wk4_int_terms)$data
plot_data$color_override <- ifelse((plot_data$conf.low <= 0 & plot_data$conf.high >= 0), "darkgrey", NA)

(hopelessness_wk4_plot <- plot_model(
  m1_hopeless,sort.est = T,terms=wk4_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    scale_color_sjplot("circus",reverse=T) +
    # Override colors for overlapping confidence intervals
    geom_point(data = subset(plot_data, !is.na(color_override)), 
               aes(x = term, y = estimate), 
               color = "darkgrey", size = 2.5) +
    geom_errorbar(data = subset(plot_data, !is.na(color_override)), 
                  aes(x = term, y = estimate, 
                      ymin = conf.low, ymax = conf.high), 
                  color = "darkgrey", width = 0) +
    ggtitle('Change in Hopelessness',
            'Baseline to Week 4, BHS-4 (range 4-12)'))

plot_data <- plot_model(m1_hopeless, sort.est = T, terms = posttest_int_terms)$data
plot_data$color_override <- ifelse((plot_data$conf.low <= 0 & plot_data$conf.high >= 0), "darkgrey", NA)

(hopelessness_posttest_plot <- plot_model(
  m1_hopeless,sort.est = T,terms=posttest_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    scale_color_sjplot("circus",reverse=T) +
    # Override colors for overlapping confidence intervals
    geom_point(data = subset(plot_data, !is.na(color_override)), 
               aes(x = term, y = estimate), 
               color = "darkgrey", size = 2.5) +
    geom_errorbar(data = subset(plot_data, !is.na(color_override)), 
                  aes(x = term, y = estimate, 
                      ymin = conf.low, ymax = conf.high), 
                  color = "darkgrey", width = 0) +
    ggtitle('Change in Hopelessness',
            'Baseline to Post-test, BHS-4 (range 4-12)'))

```

#### Frequency of Actions and Thoughts Scale (positive actions and thoughts)

```{r,warning=F,message=F}
m1_fats <- lmer(fats_score ~ time * group + (1|pid),data=mega_alltimes)

fats_cohens_d <- mixed_cohens_d(m1_fats)
add_cohens_d(m1_fats,'Predicting: change in positive actions and thoughts')
m1_fats %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')

plot_data <- plot_model(m1_fats, sort.est = T, terms = wk4_int_terms)$data
plot_data$color_override <- ifelse((plot_data$conf.low <= 0 & plot_data$conf.high >= 0), "darkgrey", NA)

(fats_wk4_plot <- plot_model(
  m1_fats,sort.est = T,terms=wk4_int_terms,auto_label=F,show.values = T, 
  value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    scale_color_sjplot("circus") +
    # Override colors for overlapping confidence intervals
    geom_point(data = subset(plot_data, !is.na(color_override)), 
               aes(x = term, y = estimate), 
               color = "darkgrey", size = 2.5) +
    geom_errorbar(data = subset(plot_data, !is.na(color_override)), 
                  aes(x = term, y = estimate, 
                      ymin = conf.low, ymax = conf.high), 
                  color = "darkgrey", width = 0) +
    ggtitle('Change in Positive Actions and Thinking','Baseline to Week 4, FATS (range 0-48)'))
```

#### Depression Change Expectancies

```{r,warning=F,message=F}
m1_dep_expect <- lmer(dep_expectancies_score ~ time * group + (1|pid),data=mega_alltimes)

dep_expect_cohens_d <- mixed_cohens_d(m1_dep_expect)
add_cohens_d(m1_dep_expect,'Predicting: change in depression change expectancies')
m1_dep_expect %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')


plot_data <- plot_model(m1_dep_expect, sort.est = T, terms = wk4_int_terms)$data
plot_data$color_override <- ifelse((plot_data$conf.low <= 0 & plot_data$conf.high >= 0), "darkgrey", NA)

(dep_expect_wk4_plot <- plot_model(
  m1_dep_expect,sort.est = T,terms=wk4_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    scale_color_sjplot("circus") +
    # Override colors for overlapping confidence intervals
    geom_point(data = subset(plot_data, !is.na(color_override)), 
               aes(x = term, y = estimate), 
               color = "darkgrey", size = 2.5) +
    geom_errorbar(data = subset(plot_data, !is.na(color_override)), 
                  aes(x = term, y = estimate, 
                      ymin = conf.low, ymax = conf.high), 
                  color = "darkgrey", width = 0) +
    ggtitle('Change in Depression Change Expectancies','Baseline to Week 4, DCES (range 3-15)'))

plot_data <- plot_model(m1_dep_expect, sort.est = T, terms = posttest_int_terms)$data
plot_data$color_override <- ifelse((plot_data$conf.low <= 0 & plot_data$conf.high >= 0), "darkgrey", NA)

(dep_expect_posttest_plot <- plot_model(
  m1_dep_expect,sort.est = T,terms=posttest_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    scale_color_sjplot("circus",reverse = T) +
    # Override colors for overlapping confidence intervals
    geom_point(data = subset(plot_data, !is.na(color_override)), 
               aes(x = term, y = estimate), 
               color = "darkgrey", size = 2.5) +
    geom_errorbar(data = subset(plot_data, !is.na(color_override)), 
                  aes(x = term, y = estimate, 
                      ymin = conf.low, ymax = conf.high), 
                  color = "darkgrey", width = 0) +
    ggtitle('Change in Depression Change Expectancies','Baseline to Post-test, DCES (range 3-15)'))

```

#### Depression Readiness for change

```{r,warning=F,message=F}
m1_readiness <- lmer(readiness_score ~ time * group + (1|group/pid),data=mega_alltimes)

readiness_cohens_d <- mixed_cohens_d(m1_readiness)
add_cohens_d(m1_readiness,'Predicting: change in readiness for change')
m1_readiness %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')


plot_data <- plot_model(m1_readiness, sort.est = T, terms = wk4_int_terms)$data
plot_data$color_override <- ifelse((plot_data$conf.low <= 0 & plot_data$conf.high >= 0), "darkgrey", NA)

(readiness_wk4_plot <- plot_model(
  m1_readiness,sort.est = T,terms=wk4_int_terms,auto_label=F,show.values = T, 
  value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    scale_color_sjplot("circus") +
    # Override colors for overlapping confidence intervals
    geom_point(data = subset(plot_data, !is.na(color_override)), 
               aes(x = term, y = estimate), 
               color = "darkgrey", size = 2.5) +
    geom_errorbar(data = subset(plot_data, !is.na(color_override)), 
                  aes(x = term, y = estimate, 
                      ymin = conf.low, ymax = conf.high), 
                  color = "darkgrey", width = 0) +
    ggtitle('Readiness for Change','Baseline to Week 4, range 2-8'))

plot_data <- plot_model(m1_readiness, sort.est = T, terms = posttest_int_terms)$data
plot_data$color_override <- ifelse((plot_data$conf.low <= 0 & plot_data$conf.high >= 0), "darkgrey", NA)

(readiness_posttest_plot <- plot_model(
  m1_readiness,sort.est = T,terms=posttest_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    scale_color_sjplot("circus",reverse=T) +
    # Override colors for overlapping confidence intervals
    geom_point(data = subset(plot_data, !is.na(color_override)), 
               aes(x = term, y = estimate), 
               color = "darkgrey", size = 2.5) +
    geom_errorbar(data = subset(plot_data, !is.na(color_override)), 
                  aes(x = term, y = estimate, 
                      ymin = conf.low, ymax = conf.high), 
                  color = "darkgrey", width = 0) +
    ggtitle('Readiness for Change','Baseline to Post-test, range 2-8'))
```


------------------------------------------------------------------------

### *H3*: Change in secondary well-being outcomes will differ between each experimental SSI condition and the active comparison SSI

-   Redo the above analysis (*H2*), but with the active comparison SSI as the comparison condition.

#### Agency
```{r,warning=F,message=F}
m1_agency_activefirst <- lmer(agency_score ~ time * group + (1|pid),data=mega_alltimes_active_first)
(agency_cohens_d_activefirst <- add_cohens_d(m1_agency_activefirst,'Predicting: change in agency'))
m1_agency_activefirst %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')


posttest_int_terms_activefirst <- m1_agency_activefirst%>% tidy() %>% filter(str_detect(term,'Post-test:')) %>% pull(term)

(agency_wk4_activefirst_plot <- plot_model(
  m1_agency_activefirst,sort.est = T,terms=wk4_int_terms_activefirst,auto_label=F,
  show.values = T,value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Change in Agency','Baseline to Week 4, vs. Behavioral Activation'))

(agency_posttest_activefirst_plot <- plot_model(
  m1_agency_activefirst,sort.est = T,
  terms=posttest_int_terms_activefirst,auto_label=F,
  show.values = T,value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Change in Agency','Baseline to Post-test, vs. Behavioral Activation'))
```

##### Hopelessness
```{r}

m1_hopeless_activefirst <- lmer(hopeless_score ~ time * group + (1|pid),data=mega_alltimes_active_first)
(hopeless_cohens_d_activefirst <- add_cohens_d(m1_hopeless_activefirst,'Predicting: change in hopelessness'))
m1_hopeless_activefirst %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')


(hopeless_wk4_activefirst_plot <- plot_model(
  m1_hopeless_activefirst,sort.est = T,terms=wk4_int_terms_activefirst,auto_label=F,
  show.values = T,value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Change in Hopelessness','Baseline to Week 4, vs. Behavioral Activation'))

(hopeless_posttest_activefirst_plot <- plot_model(
  m1_hopeless_activefirst,sort.est = T,
  terms=posttest_int_terms_activefirst,auto_label=F,
  show.values = T,value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Change in Hopelessness','Baseline to Post-test, vs. Behavioral Activation'))
```

#### Positive Actions and Thoughts (FATS)
```{r}
m1_fats_activefirst <- lmer(fats_score ~ time * group + (1|pid),data=mega_alltimes_active_first)
add_cohens_d(m1_fats_activefirst,'Predicting: change in thoughts and actions')
m1_fats_activefirst %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')


(fats_wk4_activefirst_plot <- plot_model(
  m1_fats_activefirst,sort.est = T,terms=wk4_int_terms_activefirst,auto_label=F,
  show.values = T,value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Change in Positive Actions and Thoughts','Baseline to Week 4, vs. Behavioral Activation'))
```

#### Depression Change Expectancies
```{r}
m1_dep_expect_activefirst <- lmer(dep_expectancies_score ~ time * group + (1|pid),data=mega_alltimes_active_first)
(dep_expect_cohens_d_activefirst <- add_cohens_d(m1_dep_expect_activefirst,'Predicting: change in depression change expectancies'))
m1_dep_expect_activefirst %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')


(dep_exp_wk4_activefirst_plot <- plot_model(
  m1_dep_expect_activefirst,sort.est = T,terms=wk4_int_terms_activefirst,auto_label=F,
  show.values = T,value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Change in Depression Change Expectancies','Baseline to Week 4, vs. Behavioral Activation'))

(dep_exp_posttest_activefirst_plot <- plot_model(
  m1_dep_expect_activefirst,sort.est = T,
  terms=posttest_int_terms_activefirst,auto_label=F,
  show.values = T,value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Change in Depression Change Expectancies','Baseline to Post-test, vs. Behavioral Activation'))

```

#### Depression Readiness for change
```{r}
m1_readiness_activefirst <- lmer(readiness_score ~ time * group + (1|group/pid),data=mega_alltimes_active_first)
(readiness_cohens_d_activefirst <- add_cohens_d(m1_readiness_activefirst,'Predicting: change in readiness for change'))
m1_readiness_activefirst %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')


(readiness_wk4_activefirst_plot <- plot_model(
  m1_readiness_activefirst,sort.est = T,terms=wk4_int_terms_activefirst,auto_label=F,
  show.values = T,value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Change in Readiness for Change','Baseline to Week 4, vs. Behavioral Activation'))

(readiness_posttest_activefirst_plot <- plot_model(
  m1_readiness_activefirst,sort.est = T,
  terms=posttest_int_terms_activefirst,auto_label=F,
  show.values = T,value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Change in Readiness for Change','Baseline to Post-test, vs. Behavioral Activation'))
```

------------------------------------------------------------------------

### *H4*: The experimental SSIs will differ from the active comparison SSI in acceptability

-   A linear regression model testing between-group differences in intervention acceptability for each acceptability outcome at each timepoint the outcome was collected.

    -   For these analyses, drop data from the passive control because it does not make sense to compare satisfaction between a mental health intervention and an educational program about trout.

#### Credibility / Expectancy Questionnaire

```{r,warning=F,message=F}

#Credibility / Expectancy Questionnaire scores
df_for_acceptability <- mega_alltimes_active_first %>% 
  filter(group!='Passive Control') %>% 
  mutate(group=droplevels(group))

m1_ceq <- lm(CEQ_score~group,data=df_for_acceptability)
(ceq_cohens_d <- rstatix::cohens_d(CEQ_score~group,data=df_for_acceptability) %>% select(-magnitude,-.y.)) %>%kable()

acceptability_terms <- m1_ceq %>% tidy() %>% 
  filter(term!='(Intercept)') %>% pull(term)
(ceq_plot <- plot_model(
  m1_ceq,sort.est = T,terms=acceptability_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Credibility and Expected Benefit','At Post-test, range 4-56'))

ba_stats <- mega_alltimes %>% 
  filter(time=='Post-test'&group == "Behavioral Activation") %>%
    summarize(
      mean_CEQ_score = mean(CEQ_score, na.rm = TRUE),
      se = sd(CEQ_score, na.rm = TRUE) / sqrt(n()),
      ci_lower = mean_CEQ_score - 1.96 * se,
      ci_upper = mean_CEQ_score + 1.96 * se)

(ceq_pointplot <- mega_alltimes %>% filter(time=='Post-test'&group!='Passive Control') %>% 
    group_by(group) %>% 
    summarize(
      mean_CEQ_score = mean(CEQ_score, na.rm = TRUE),
      se = sd(CEQ_score, na.rm = TRUE) / sqrt(n()),
      ci_lower = mean_CEQ_score - 1.96 * se,
      ci_upper = mean_CEQ_score + 1.96 * se) %>% 
    ggplot(aes(x = reorder(group, mean_CEQ_score), y = mean_CEQ_score)) +
    geom_point(size=1.5) +
    geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper),
                  color='green3',width = 0.2,show.legend = F) + 
  # Highlight the CI range of Behavioral Activation
  annotate("rect",
           xmin = -Inf, xmax = Inf,
           ymin = ba_stats$ci_lower, ymax = ba_stats$ci_upper,
           fill = "gray90", alpha = 0.5)+
    coord_flip() +labs(x = "", y = "Mean CEQ Score")+
    ggtitle("Credibility and Expected Benefit","At Post-test, range 4-56"))

```

#### Star rating at post-test

```{r,message=F,warning=F}
m1_star <- lm(star_rating_1~group,data=df_for_acceptability %>% 
                filter(time=='Post-test'))
(star_cohens_d <- rstatix::cohens_d(star_rating_1~group,data=df_for_acceptability %>% filter(time=='Post-test'))%>% select(-magnitude,-.y.)) %>% kable
(star_plot <- plot_model(
  m1_star,sort.est = T,terms=acceptability_terms,
  auto_label=F,show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Star Rating','At Post-test, compared to a Behavioral Activation SSI'))

ba_stats <- mega_alltimes %>% 
  filter(time=='Post-test'&group == "Behavioral Activation") %>%
    summarize(
      mean_star_rating = mean(star_rating_1, na.rm = TRUE),
      se = sd(star_rating_1, na.rm = TRUE) / sqrt(n()),
      ci_lower = mean_star_rating - 1.96 * se,
      ci_upper = mean_star_rating + 1.96 * se)

(star_pointplot <- mega_alltimes %>% filter(time=='Post-test'&group!='Passive Control') %>% 
    group_by(group) %>% 
    summarize(
      mean_star_score = mean(star_rating_1, na.rm = TRUE),
      se = sd(star_rating_1, na.rm = TRUE) / sqrt(n()),
      ci_lower = mean_star_score - 1.96 * se,
      ci_upper = mean_star_score + 1.96 * se) %>% 
    ggplot(aes(x = reorder(group, mean_star_score), y = mean_star_score)) +
    geom_point(size=1.5) +
    geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper),
                  color='green3',width = 0.2,show.legend = F) + 
      # Highlight the CI range of Behavioral Activation
  annotate("rect",
           xmin = -Inf, xmax = Inf,
           ymin = ba_stats$ci_lower, ymax = ba_stats$ci_upper,
           fill = "gray90", alpha = 0.5)+
    coord_flip() +labs(x = "", y = "Mean Star Rating")+
    ggtitle("Star Rating","At Post-test, range 1-5"))


```

#### Star Rating at Week 4

```{r,message=F,warning=F}
m1_star_wk4 <- lm(star_rating_1~group,data=df_for_acceptability %>% 
                    filter(time=='Week 4'))
rstatix::cohens_d(star_rating_1~group,data=df_for_acceptability %>% 
                    filter(time=='Week 4')) %>% select(-magnitude,-.y.) %>%kable()
(star_plot_wk4 <- plot_model(m1_star_wk4,sort.est = T,terms=acceptability_terms,auto_label=F,show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Star Rating','At Week 4, compared to a Behavioral Activation SSI'))
```

#### Frequency of Aha! moments at post-test

```{r,message=F,warning=F}
m1_aha <- lm(aha~group,data=df_for_acceptability %>% filter(time=='Post-test'))
(aha_cohens_d <- rstatix::cohens_d(aha~group,data=df_for_acceptability%>% filter(time=='Post-test')) %>% select(-magnitude,-.y.)) %>%kable()
(aha_plot <- plot_model(m1_aha,sort.est = T,terms=acceptability_terms,auto_label=F,show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Probability of an Aha! Moment during SSI','At Post-test, compared to a Behavioral Activation SSI'))

ba_stats <- mega_alltimes %>% 
  filter(time=='Post-test'&group == "Behavioral Activation") %>%
    summarize(
      mean_aha_score = mean(aha, na.rm = TRUE),
      se = sd(aha, na.rm = TRUE) / sqrt(n()),
      ci_lower = mean_aha_score - 1.96 * se,
      ci_upper = mean_aha_score + 1.96 * se)

(aha_pointplot <- mega_alltimes %>% 
    filter(time=='Post-test'&group!='Passive Control') %>% 
    group_by(group) %>% 
    summarize(
        aha_prop = mean(aha, na.rm = TRUE),
        n = sum(!is.na(aha)),
        se = sqrt(aha_prop * (1 - aha_prop) / n),
        ci_lower = aha_prop - 1.96 * se,
        ci_upper = aha_prop + 1.96 * se
    ) %>% 
    ggplot(aes(x = reorder(group, aha_prop), y = aha_prop)) +
    geom_point(size=1.5) +
    geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper),
                  color='palegreen2',width = 0.2,show.legend = F) + 
    labs(x = "", y = "Proportion who had an Aha! Moment During SSI")+
      # Highlight the CI range of Behavioral Activation
  annotate("rect",
           xmin = -Inf, xmax = Inf,
           ymin = ba_stats$ci_lower, ymax = ba_stats$ci_upper,
           fill = "gray90", alpha = 0.5)+
    scale_y_continuous(labels = scales::percent)+
    coord_flip() +
    ggtitle("Proportion of Aha! Moments","At Post-test"))
```

#### Frequency of Aha! moments at week 4

```{r,message=F,warning=F}
m1_aha_wk4 <- lm(aha_wk4~group,data=df_for_acceptability %>% filter(time=='Week 4')) 
rstatix::cohens_d(aha_wk4~group,data=df_for_acceptability%>% filter(time=='Week 4')) %>% select(-magnitude,-.y.) %>%kable()
(aha_plot_wk4 <- plot_model(
  m1_aha_wk4,sort.est = T,terms=acceptability_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+scale_x_discrete(
    labels = function(x) gsub(".*group", "", x)) +    
    ggtitle('Probability of an Aha! Moment during SSI','At Week 4'))
```

#### Significance of Aha! moments at post-test

```{r,message=F,warning=F}
m1_aha_significance <- lm(aha_significant~group,data=df_for_acceptability %>% 
                                filter(time=='Post-test'))

(aha_significance_cohens_d <- rstatix::cohens_d(
  aha_significant~group,data=df_for_acceptability) %>% select(-magnitude,-.y.)) %>%kable()
(aha_significance_plot <- plot_model(
  m1_aha_significance,sort.est = T,terms=acceptability_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Significance of Aha! Moment','At Post-test, compared to a Behavioral Activation SSI'))

ba_stats <- mega_alltimes %>% 
  filter(time=='Post-test'&group == "Behavioral Activation") %>%
    summarize(
      mean_aha_score = mean(aha_significant, na.rm = TRUE),
      se = sd(aha_significant, na.rm = TRUE) / sqrt(n()),
      ci_lower = mean_aha_score - 1.96 * se,
      ci_upper = mean_aha_score + 1.96 * se)

(aha_signif_pointplot <- mega_alltimes %>% filter(time=='Post-test'&group!='Passive Control') %>% 
  group_by(group) %>% 
    summarize(mean_aha_signif_score = mean(aha_significant, na.rm = TRUE),
      se = sd(aha_significant, na.rm = TRUE) / sqrt(n()),
      ci_lower = mean_aha_signif_score - 1.96 * se,
      ci_upper = mean_aha_signif_score + 1.96 * se) %>%
  ggplot(aes(x = reorder(group, mean_aha_signif_score), y = mean_aha_signif_score)) +
  geom_point(size=1.5) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper),
                color='palegreen3',width = 0.2,show.legend = F) + 
      # Highlight the CI range of Behavioral Activation
  annotate("rect",
           xmin = -Inf, xmax = Inf,
           ymin = ba_stats$ci_lower, ymax = ba_stats$ci_upper,
           fill = "gray90", alpha = 0.5)+
    coord_flip() +labs(x = "", y = "Average Aha! Moment Significance")+
    ggtitle("Significance of Aha! Moment","At Post-test, range 0-5"))
```

#### Significance of Aha! moments at week 4

```{r,message=F,warning=F}
m1_aha_significance_wk4 <- lm(
  aha_significant~group,data=df_for_acceptability %>% 
                                filter(time=='Week 4'))
rstatix::cohens_d(aha_significant~group,data=df_for_acceptability %>% 
                    filter(time=='Week 4')) %>%select(-magnitude,-.y.) %>% kable()
(aha_significance_plot_wk4 <- plot_model(
  m1_aha_significance_wk4,sort.est = T,terms=acceptability_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Significance of Aha! Moment','At Week 4, compared to a Behavioral Activation SSI'))
```

#### How much Aha! changed behavior at 4-week follow-up

```{r,message=F,warning=F}
m1_aha_howchanged <- lm(aha_howchanged~group,data=df_for_acceptability)
rstatix::cohens_d(aha_howchanged~group,data=df_for_acceptability) %>%select(-magnitude,-.y.) %>% kable()
(aha_howchanged_plot <- plot_model(
  m1_aha_howchanged,sort.est = T,terms=acceptability_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('How Much Aha! Moment Changed behavior','Four Weeks Later, compared to a Behavioral Activation SSI'))
```

------------------------------------------------------------------------

### *H5*: Expecting one can improve one's depression more at baseline will be associated with different change

-   A linear mixed-effects regression model with depressive symptoms as the DV, timepoint, baseline depression expectancies, and the interaction between the two as predictors, and a participant identifier nested in experimental condition as a random intercept.

```{r,message=F,warning=F}
m1_expect_on_phq9 <- lmer(
  phq9_score ~ time * dep_expect_bl + (1|group/pid),
  data=mega_alltimes)
summary(m1_expect_on_phq9)

mega_alltimes %>% select(pid,time,dep_expect_bl,phq9_score) %>% 
  pivot_wider(names_from=time,values_from = phq9_score,names_prefix='phq9_score_') %>% 
  mutate(phq9_change=`phq9_score_Week 4`-phq9_score_Baseline) %>%
  ggplot(aes(x=dep_expect_bl,y=phq9_change))+geom_smooth()+
  xlab('Depression change expectancies at Baseline')+
  ylab('Change in PHQ-9 Score from Baseline to Week 4')+
  ggtitle('Change Expectancies Predicting Change in Depression')
```

------------------------------------------------------------------------

### *H6*: The active SSIs will be better than the Passive SSI on average for all mental health outcomes

-   Re-run the main analysis with experimental condition re-coded into a two-level factor in which any of the 12 active SSIs are coded as “active.” The random intercept, participant ID is nested within experimental condition.

#### Depression

```{r,message=F,warning=F}
mega_alltimes <- mega_alltimes %>% mutate(passive_vs_all=ifelse(
  group=='Passive Control','Passive Control','Active Conditions\n(Average)')%>% factor %>% fct_rev())

m1_phq9_all_vs_passive <- lmer(phq9_score ~ time * passive_vs_all + (1|group/pid),data=mega_alltimes)
add_cohens_d(m1_phq9_all_vs_passive,'Predicting: change in depressive symptoms')

(liness <- plot_model(m1_phq9_all_vs_passive,type='pred',terms=c('time','passive_vs_all'),ci.lvl = .95)+
    geom_line()+ggtitle('Predicted Depressive Symptoms Scores')+
    xlab(element_blank())+ 
    ylab('Depression (PHQ-9)')+ 
    guides(color=guide_legend(title="Active vs. Passive")))
#ggsave("phq9_changeplot_lines.png",last_plot(),height=6,width=6)

mega_alltimes %>% filter(time!='Post-test') %>% 
  mutate(time=fct_drop(time)) %>% 
  ggplot(aes(x=time,y=phq9_score,fill=passive_vs_all, col=passive_vs_all))+ 
  geom_boxplot(alpha=0.8,outlier.shape = NA)+ 
  geom_jitter(alpha=0.1,size=.3)+
  ylab('Depression (PHQ-9)')+ xlab(element_blank())+
  ggtitle('Depression Score by Time and Condition',
          'Passive Control vs All Active Conditions Averaged')+
  scale_fill_manual(values = c("grey", 'darkolivegreen1'))

```

##### Distributions of change in depression over time #####

```{r,message=F,warning=F}
mega_alltimes %>% select(pid,phq9_score,time,passive_vs_all) %>% 
    pivot_wider(names_from='time',values_from = 'phq9_score') %>% 
    mutate(phq_change=`Week 4`-Baseline) %>% ggplot(aes(x=phq_change,color=passive_vs_all))+geom_density()

##### drop in phq over 5 points 
mega_alltimes %>% select(pid,time,dep_expect_bl,phq9_score,group) %>% 
  pivot_wider(names_from=time,values_from = phq9_score,names_prefix='phq9_score_') %>% 
  mutate(phq_change=`phq9_score_Week 4`-phq9_score_Baseline) %>% group_by(group)%>% 
  count(!is.na(phq_change)&phq_change<=-5,name='count_downby5') %>% 
  mutate(pct_phq_went_down_5=count_downby5/sum(count_downby5)) %>%
  filter(`!is.na(phq_change) & phq_change <= -5`=='TRUE')

##### rise in phq over 5 points 
mega_alltimes %>% select(pid,time,phq9_score,group) %>% 
  pivot_wider(names_from=time,values_from = phq9_score,names_prefix='phq9_score_') %>% 
  mutate(phq_change=`phq9_score_Week 4`-phq9_score_Baseline) %>% group_by(group)%>% 
  count(!is.na(phq_change)&phq_change>=5,name='count_upby5') %>% 
  mutate(pct_phq_went_up_5=count_upby5/sum(count_upby5)) %>%
  filter(`!is.na(phq_change) & phq_change >= 5`=='TRUE')

```

It kinda looks like most of the difference is around the mean. Fewer people got a little worse, and more people got a little better. But it's slight!

#### Agency

```{r,message=F,warning=F}

m1_agency_all_vs_passive <- lmer(agency_score ~ time * passive_vs_all + (1|group/pid),data=mega_alltimes)
add_cohens_d(m1_agency_all_vs_passive,'Predicting: change in agency')

#Change in average SS effect on Agency from post-test to week 4
mixed_cohens_d(m1_agency_all_vs_passive) %>% 
  filter(str_detect(term,'Week 4:')) %>% pull(cohens_d) / 
  mixed_cohens_d(m1_agency_all_vs_passive) %>% 
  filter(str_detect(term,'Post-test:')) %>% pull(cohens_d) 

(liness_agency <- plot_model(m1_agency_all_vs_passive,type='pred',terms=c('time','passive_vs_all'),ci.lvl = .95)+
    geom_line()+ggtitle('Predicted Agency Scores')+
    xlab(element_blank())+ 
    ylab('Agency')+ 
    guides(color=guide_legend(title="Active vs. Passive")))

mega_alltimes %>% 
  ggplot(aes(x=time,y=agency_score,fill=passive_vs_all, col=passive_vs_all))+ 
  geom_boxplot(alpha=0.8,outlier.shape = NA)+ 
  geom_jitter(alpha=0.1,size=.3)+
  ylab('Agency')+ xlab(element_blank())+
  ggtitle('Agency by Time and Condition',
          'Passive Control vs All Active Conditions Averaged')+
  scale_fill_manual(values = c("grey",'darkolivegreen1'))

```

#### Hopelessness

```{r,message=F,warning=F}
m1_hopeless_all_vs_passive <- lmer(hopeless_score ~ time * passive_vs_all + (1|group/pid),data=mega_alltimes)
add_cohens_d(m1_hopeless_all_vs_passive,'Predicting: change in hopelessness')

(liness_hopelessness <- plot_model(m1_hopeless_all_vs_passive,type='pred',terms=c('time','passive_vs_all'),ci.lvl = .95)+
    geom_line()+ggtitle('Predicted Hopelessness Scores')+
    xlab(element_blank())+ 
    ylab('Hopelessness')+ 
    guides(color=guide_legend(title="Active vs. Passive")))

mega_alltimes %>% 
  ggplot(aes(x=time,y=hopeless_score,fill=passive_vs_all, col=passive_vs_all))+ 
  geom_boxplot(alpha=0.8,outlier.shape = NA)+ 
  geom_jitter(alpha=0.1,size=.3)+
  ylab('Hopelessness')+ xlab(element_blank())+
  ggtitle('Hopelessness by Time and Condition',
          'Passive Control vs All Active Conditions Averaged')+
  scale_fill_manual(values = c("grey",'darkolivegreen1'))
```

#### Frequency of Actions and Thoughts Scale

```{r,message=F,warning=F}
m1_fats_all_vs_passive <- lmer(fats_score ~ time * passive_vs_all + (1|group/pid),data=mega_alltimes)
add_cohens_d(m1_fats_all_vs_passive,'Predicting: change in actions and thoughts')

(liness_fats <- plot_model(m1_fats_all_vs_passive,type='pred',terms=c('time','passive_vs_all'),ci.lvl = .95)+
    geom_line()+ggtitle('Predicted positive actions and thoughts Scores')+
    xlab(element_blank())+ 
    ylab('FATS score')+ 
    guides(color=guide_legend(title="Active vs. Passive")))

mega_alltimes %>% filter(time!='Post-test') %>% 
  mutate(time=fct_drop(time)) %>% 
  ggplot(aes(x=time,y=fats_score,fill=passive_vs_all, col=passive_vs_all))+ 
  geom_boxplot(alpha=0.8,outlier.shape = NA)+ 
  geom_jitter(alpha=0.1,size=.3)+
  ylab('Agency')+ xlab(element_blank())+
  ggtitle('Positive Actions and Thinking by Time and Condition',
          'Passive Control vs All Active Conditions Averaged')+
  scale_fill_manual(values = c("grey",'darkolivegreen1'))
```

#### Depression Change Expectancies

```{r,message=F,warning=F}
m1_dep_expect_all_vs_passive <- lmer(dep_expectancies_score ~ time * passive_vs_all + (1|group/pid),data=mega_alltimes)
add_cohens_d(m1_dep_expect_all_vs_passive,'Predicting: change in depression change expectancies')

(liness_change_expectancies <- plot_model(m1_dep_expect_all_vs_passive,type='pred',terms=c('time','passive_vs_all'),ci.lvl = .95)+
    geom_line()+ggtitle('Predicted Depression Change Expectancies')+
    xlab(element_blank())+ 
    ylab('Change expectancy')+ 
    guides(color=guide_legend(title="Active vs. Passive")))

mega_alltimes %>%  
  ggplot(aes(x=time,y=dep_expectancies_score,fill=passive_vs_all, col=passive_vs_all))+ 
  geom_boxplot(alpha=0.8,outlier.shape = NA)+ 
  geom_jitter(alpha=0.1,size=.3)+
  ylab('Agency')+ xlab(element_blank())+
  ggtitle('Depression Change Expectancies by Time and Condition',
          'Passive Control vs All Active Conditions Averaged')+
  scale_fill_manual(values = c("grey",'darkolivegreen1'))
```

#### Readiness for change

```{r,message=F,warning=F}
m1_readiness_all_vs_passive <- lmer(readiness_score ~ time * passive_vs_all + (1|group/pid),data=mega_alltimes)
add_cohens_d(m1_readiness_all_vs_passive,'Predicting: change in readiness for change')

(liness_readiness <- plot_model(m1_readiness_all_vs_passive,type='pred',terms=c('time','passive_vs_all'),ci.lvl = .95)+
    geom_line()+ggtitle('Predicted Readiness for Change')+
    xlab(element_blank())+ 
    ylab('Readiness for Change')+ 
    guides(color=guide_legend(title="Active vs. Passive")))

mega_alltimes %>%  
  ggplot(aes(x=time,y=readiness_score,fill=passive_vs_all, col=passive_vs_all))+ 
  geom_boxplot(alpha=0.8,outlier.shape = NA)+ 
  geom_jitter(alpha=0.1,size=.3)+
  ylab('Agency')+ xlab(element_blank())+
  ggtitle('Readiness for Change by Time and Condition',
          'Passive Control vs All Active Conditions Averaged')+
  scale_fill_manual(values = c("grey",'darkolivegreen1'))
```

------------------------------------------------------------------------

------------------------------------------------------------------------

# Sensitivity Analyses (Pre-registered)

------------------------------------------------------------------------

### Including demographic covariates

-   Re-run the primary analysis with demographic covariates: age, gender, disability, race, education, relationship status, political party, employment status, household income, and perceived socioeconomic standing (Social Ladder).

```{r,warning=F,message=F,rows.print = 100}
m1_phq9_w_dem_covars <- lmer(phq9_score ~ time * group + Age + Gender + Race +
                               Education + 
                               `Relationship/Marital Status` + 
                               `Political Party` + `Employment Status` + 
                               `Household Income` + `Social ladder`+
                               disability +(1|pid),data=mega_alltimes)

add_cohens_d(m1_phq9_w_dem_covars,caption='Predicting: change in depressive symptoms')
m1_phq9_w_dem_covars %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')


(phq9_with_demog_covars_plot <- plot_model(m1_phq9_w_dem_covars,sort.est = T,terms=wk4_int_terms,auto_label=F,
                                           show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Change in Depressive Symptoms','With demographic covariates, vs. the control'))
```

------------------------------------------------------------------------

### Including both demographic and baseline psychological covariates

-   Re-run the primary analysis with the demographics above and baseline mental health covariates: depression change expectancies, positive actions and thoughts, readiness for change, agency, and hopelessness.

```{r,warning=F,message=F,rows.print = 100}
m1_phq9_w_all_covars <- lmer(phq9_score ~ `time` * `group` + Age + Gender + Race +
                               Education + `Relationship/Marital Status` + 
                               `Political Party` + `Employment Status` + 
                               `Household Income` + `Social ladder`+
                               disability + dep_expect_bl+ readiness_bl+ 
                               agency_bl + hopeless_bl+ fats_bl + 
                               (1|pid),data=mega_alltimes)

add_cohens_d(m1_phq9_w_all_covars,caption='Predicting: change in depressive symptoms')
m1_phq9_w_all_covars %>% tidy %>% select(term,statistic,df) %>% kable(digits=3,caption='Model T score and df')


(phq9_with_all_covars_plot <- plot_model(m1_phq9_w_all_covars,sort.est = T,terms=wk4_int_terms,auto_label=F,
                                         show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Change in Depressive Symptoms','With demographic and baseline well-being covariates, vs. passive control'))
```

Adding the demographic and baseline covariates really doesn't change the model much.

------------------------------------------------------------------------

------------------------------------------------------------------------

# Robustness Check: Correcting for multiple comparisons (Pre-registered)

-   Using the Benjamini-Hochberg procedure to re-rerun the models correcting for multiple comparisons to produce adjusted p-values.

```{r,message=F,warning=F,rows.print = 40}
#Depression
m1_phq9 %>% tidy() %>% filter(str_detect(term,':')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)%>% print(n=100)

plot_model(m1_phq9,sort.est = T,terms=wk4_int_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
  ggtitle('Depressive Symptoms baseline to week 4','With multiple comparison correction')

#Agency
m1_agency %>% tidy() %>% filter(str_detect(term,':group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj) %>% print(n=100)

plot_model(m1_agency,sort.est = T,terms=wk4_int_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
  ggtitle('Agency baseline to week 4','With multiple comparison correction')

plot_model(m1_agency,sort.est = T,terms=posttest_int_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
  ggtitle('Agency baseline to post-test','With multiple comparison correction')

#Hopelessness
m1_hopeless %>% tidy() %>% filter(str_detect(term,':group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)%>% print(n=100)

plot_model(m1_hopeless,sort.est = T,terms=wk4_int_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
  ggtitle('Hopelessness baseline to week 4','With multiple comparison correction')

plot_model(m1_hopeless,sort.est = T,terms=posttest_int_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
  ggtitle('Hopelessness baseline to post-test','With multiple comparison correction')

#Frequency of actions and Thoughts
m1_fats %>% tidy() %>% filter(str_detect(term,':group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)%>% print(n=100)

plot_model(m1_fats,sort.est = T,terms=wk4_int_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
  ggtitle('FATS baseline to week 4','With multiple comparison correction')

#Depression Change Expectancies
m1_dep_expect %>% tidy() %>% filter(str_detect(term,':group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)%>% print(n=100)

plot_model(m1_dep_expect,sort.est = T,terms=wk4_int_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
  ggtitle('Change Expectancies baseline to week 4',
          'With multiple comparison correction')

plot_model(m1_dep_expect,sort.est = T,terms=posttest_int_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
  ggtitle('Change Expectancies baseline to post-test',
          'With multiple comparison correction')

#Readiness to change
m1_readiness %>% tidy() %>% filter(str_detect(term,':group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)%>% print(n=100)

plot_model(m1_readiness,sort.est = T,terms=wk4_int_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
  ggtitle('Readiness to change baseline to week 4',
          'With multiple comparison correction')

plot_model(m1_readiness,sort.est = T,terms=posttest_int_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
  ggtitle('Readiness to change baseline to post-test',
          'With multiple comparison correction')

#Credibility / Expectancy Questionnaire
m1_ceq %>% tidy() %>% filter(str_detect(term,'group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)%>% print(n=100)

(ceq_plot <- plot_model(m1_ceq,sort.est = T,terms=acceptability_terms,auto_label=F,
                        show.values = T, value.offset = .4,value.size = 3.2,
                        p.adjust = 'hochberg')+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Credibility and Expected Benefit','With multiple comparison correction'))

#Star ratings
m1_star %>% tidy() %>% filter(str_detect(term,'group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)%>% print(n=100)

plot_model(m1_star,sort.est = T,terms=acceptability_terms,auto_label=F,
                        show.values = T, value.offset = .4,value.size = 3.2,
                        p.adjust = 'hochberg')+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Star rating at post-test','With multiple comparison correction')

m1_star_wk4 %>% tidy() %>% filter(str_detect(term,'group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)

plot_model(m1_star_wk4,sort.est = T,terms=acceptability_terms,auto_label=F,
                         show.values = T, value.offset = .4,value.size = 3.2,
                         p.adjust = 'hochberg')+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Star rating at Week 4','With multiple comparison correction')


#Aha moments
m1_aha %>% tidy() %>% filter(str_detect(term,'group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)%>% print(n=100)

plot_model(m1_aha,sort.est = T,terms=acceptability_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
  ggtitle('Aha moment at post-test','With multiple comparison correction')

m1_aha_wk4 %>% tidy() %>% filter(str_detect(term,'group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)

plot_model(m1_aha_wk4,sort.est = T,terms=acceptability_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
  ggtitle('Aha moment at week 4','With multiple comparison correction')

#Aha moment significance
m1_aha_significance %>% tidy() %>% filter(str_detect(term,'group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)%>% print(n=100)

plot_model(m1_aha_significance,sort.est = T,terms=acceptability_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
  ggtitle('Aha moment significance at post-test','With multiple comparison correction')

m1_aha_significance_wk4 %>% tidy() %>% filter(str_detect(term,'group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)%>% print(n=100)

plot_model(m1_aha_significance_wk4,sort.est = T,terms=acceptability_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
  ggtitle('Aha moment significance at week 4','With multiple comparison correction')

#How much aha moment changed behavior four weeks later
m1_aha_howchanged %>% tidy() %>% filter(str_detect(term,'group')) %>% 
  adjust_pvalue(method='BH') %>% select(term,p.value,p.value.adj)

plot_model(m1_aha_howchanged,sort.est = T,terms=acceptability_terms,auto_label=F,
           show.values = T, value.offset = .4,value.size = 3.2,
           p.adjust = 'hochberg')+
  scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
  ggtitle("How much Aha! moment changed one's behavior at week 4",'With multiple comparison correction')

```

------------------------------------------------------------------------

------------------------------------------------------------------------

# Additional Analyses (Not pre-registered)

### *1*. Intervention features predicting outcomes

-   Before collecting data, we coded the SSIs on several characteristics (<https://osf.io/5kxqz>)

-   This analysis is massively underpowered (n=12), but may be worth a look anyways.

```{r,message=F,warning=F}
#each of these compares an outcome by the presence of a feature

#loading our ratings of each SSI's characteristics
ssi_characteristics <- read_excel('top10 ssi characteristics all three ratings.xlsx')

ceq_avgs <- mega_alltimes %>% filter(time=='Post-test') %>% group_by(group) %>%
  summarize(CEQ_mean=mean(CEQ_score,na.rm=T)) 
#joining study effect sizes and our team's ratings of SSI characteristics
chars_plus_studydata <- ceq_avgs %>% left_join(
    mixed_cohens_d(m1_phq9) %>% 
      filter(str_detect(term,'timeWeek 4:')) %>% 
      select(group=term,PHQ_ES='cohens_d') %>% 
      mutate(group=str_remove(group,'timeWeek 4:group'))) %>% 
  left_join(ssi_characteristics,by=c('group'='SSI')) %>% 
  rename_with(~ str_replace(., " \\(.*", ""))

# Define predictors and outcomes
predictors <- c(
  "CEQ_mean","`Video or audio mostly`", "Engagingness",
  "`Interactive activities`",  "`Personalization / choice`",
  "Effortful", "`Refer to other resources`", "`Meditation practice`",
  "`Cognitive reframing`", "`Behavioral activation`", "`Mindfulness`",
  "`Clinician-created mostly`", "`take-home reminder`")

# Create a dataframe to store results
run_models <- function(outcome) {
  do.call(rbind, lapply(predictors, function(predictor) {
    model <- lm(as.formula(paste(outcome, "~", predictor)), data = chars_plus_studydata)
    tidy(model) %>% 
      filter(term != "(Intercept)") %>% 
      adjust_pvalue(method='BH')# Remove intercept
  }))
}

# Run models for each outcome and store in separate dataframes
(results_PHQ <- run_models("PHQ_ES")) #predicting PHQ_ES
(results_CEQ <- run_models("CEQ_mean")) #predicting CEQ_mean

ggplot(results_PHQ, aes(x = reorder(term, estimate), y = estimate)) +
  geom_bar(stat = "identity", fill = "blue3") +
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error),
                width = 0.2) + coord_flip() +
  labs(x = "Term", y = "Estimate \u00B1 SE",
       title = "Intervention characteristics predicting change in depression")

ggplot(results_CEQ, aes(x = reorder(term, estimate), y = estimate)) +
  geom_bar(stat = "identity", fill = "blue3") +
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error),
                width = 0.2) + coord_flip() +
  labs(x = "Term", y = "Estimate \u00B1 SE",
       title = "SSI characteristics predicting acceptability",'At Post-test')
```

### *2*. Correlation between intervention effects

This analysis is also super underpowered.

```{r,message=F,warning=F}

#Week 4 correlations
ESs_outcomes_wk4 <- phq9_cohens_d %>%  filter(str_detect(term,'Week 4:group')) %>%
  rename(PHQ9_ES=cohens_d) %>% 
  left_join(agency_cohens_d %>% filter(str_detect(term,'Week 4:group')) %>% 
              rename(agency_ES=cohens_d),by='term') %>% 
  left_join(hopeless_cohens_d %>% filter(str_detect(term,'Week 4:group')) %>%  
              rename(hopeless_ES=cohens_d),by='term') %>% 
  left_join(fats_cohens_d %>% filter(str_detect(term,'Week 4:group')) %>%  
              rename(FATS_ES=cohens_d),by='term') %>% 
  left_join(dep_expect_cohens_d %>% filter(str_detect(term,'Week 4:group')) %>%  
              rename(dep_expect_ES=cohens_d),by='term') %>% 
  left_join(readiness_cohens_d %>% filter(str_detect(term,'Week 4:group')) %>%  
              rename(readiness_ES=cohens_d),by='term') %>% 
  mutate(term=str_remove(term,'timeWeek 4:group')) %>% 
  select(ends_with('_ES'))

corr_wk4 <- round(cor(ESs_outcomes_wk4), 2)
ggcorrplot(corr_wk4, hc.order = TRUE, type = "lower",
           outline.col = "white",lab=T,
           title = 'Correlation between intervention effect sizes over 4 weeks')

# Post-test correlations

star_avgs <- mega_alltimes %>% filter(time=='Post-test') %>% group_by(group) %>% 
  summarize(star_mean=mean(star_rating_1,na.rm=T)) 
aha_avgs <- mega_alltimes %>% filter(time=='Post-test') %>% group_by(group) %>% 
  summarize(aha_mean=mean(aha,na.rm=T)) 
aha_significance_avgs <- mega_alltimes %>% filter(time=='Post-test') %>% group_by(group) %>% 
  summarize(aha_significance_mean=mean(aha_significant,na.rm=T)) 

ESs_outcomes_pt <- agency_cohens_d %>% 
  filter(str_detect(term,'Post-test:group')) %>% 
  rename(agency_ES=cohens_d) %>% 
  left_join(hopeless_cohens_d %>% filter(str_detect(term,'Post-test:group')) %>%  
              rename(hopeless_ES=cohens_d),by='term') %>% 
  left_join(dep_expect_cohens_d %>% filter(str_detect(term,'Post-test:group')) %>%
              rename(dep_expect_ES=cohens_d),by='term') %>% 
  left_join(readiness_cohens_d %>% filter(str_detect(term,'Post-test:group')) %>%  
              rename(readiness_ES=cohens_d),by='term') %>% 
  mutate(term=str_remove(term,'timePost-test:group')) %>% 
  rename(group=term) %>% 
  left_join(ceq_avgs) %>% 
  left_join(star_avgs) %>% 
  left_join(aha_avgs) %>% 
  left_join(aha_significance_avgs) %>% 
  select(ends_with('_ES')|ends_with('_mean')) 
#correlating acceptability outcomes 

corr_pt <- round(cor(ESs_outcomes_pt), 2)
ggcorrplot(corr_pt, hc.order = TRUE, type = "lower",
           outline.col = "white",lab=T,
           title = 'Correlation between intervention effect sizes at post-test')
```

```{r,message=F,warning=F}
#Make figures showing effect sizes per group across each outcome
ESs_outcomes_all_with_names <- phq9_cohens_d %>%  filter(str_detect(term,'Week 4:group')) %>%
  rename(PHQ9_Week4_ES=cohens_d) %>% 
  mutate(term=str_remove(term,'.*:group')) %>% 
  left_join(agency_cohens_d %>% filter(str_detect(term,'Week 4:group')) %>% 
              mutate(term=str_remove(term,'.*:group')) %>% 
              rename(Agency_Week4_ES=cohens_d),by='term') %>% 
  left_join(hopeless_cohens_d %>% filter(str_detect(term,'Week 4:group')) %>% 
              mutate(term=str_remove(term,'.*:group')) %>%  
              rename(Hopelessness_Week4_ES=cohens_d),by='term') %>% 
  left_join(fats_cohens_d %>% filter(str_detect(term,'Week 4:group')) %>% 
              mutate(term=str_remove(term,'.*:group')) %>%  
              rename(FATS_Week4_ES=cohens_d),by='term') %>% 
  left_join(dep_expect_cohens_d %>% filter(str_detect(term,'Week 4:group')) %>% 
              mutate(term=str_remove(term,'.*:group')) %>%  
              rename(Change_Expectancy_Week4_ES=cohens_d),by='term') %>% 
  left_join(readiness_cohens_d %>% filter(str_detect(term,'Week 4:group')) %>% 
              mutate(term=str_remove(term,'.*:group')) %>%  
              rename(Readiness_Week4_ES=cohens_d),by='term') %>% 
  left_join(agency_cohens_d %>% filter(str_detect(term,'Post-test:group')) %>% 
              mutate(term=str_remove(term,'.*:group')) %>% 
              rename(Agency_Posttest_ES=cohens_d),by='term') %>% 
  left_join(hopeless_cohens_d %>% filter(str_detect(term,'Post-test:group')) %>% 
              mutate(term=str_remove(term,'.*:group')) %>%  
              rename(Hopelessness_Posttest_ES=cohens_d),by='term') %>% 
  left_join(dep_expect_cohens_d %>% filter(str_detect(term,'Post-test:group')) %>% 
              mutate(term=str_remove(term,'.*:group')) %>%  
              rename(Change_Expectancy_Posttest_ES=cohens_d),by='term') %>% 
  left_join(readiness_cohens_d %>% filter(str_detect(term,'Post-test:group')) %>% 
              mutate(term=str_remove(term,'.*:group')) %>%  
              rename(Readiness_Posttest_ES=cohens_d),by='term') %>% 
  select(term,ends_with('_ES'))%>% rename_with(~str_remove(., '_ES'))


ESs_outcomes_all_with_names_long <- ESs_outcomes_all_with_names %>%
  pivot_longer(
    cols = -term,
    names_to = "Outcome",
    values_to = 'Effect') %>% mutate(
      Effect =ifelse(str_detect(Outcome,'PHQ|Hope'),
                     -1*Effect,Effect)) %>% 
  mutate(Outcome=str_replace(
    Outcome,'PHQ9_Week4','Depression_Week4')) %>% 
  mutate(Outcome=str_replace(
    Outcome,'FATS_Week4','+_Actions/Thinking_Week4')) %>% 
  mutate(Outcome=factor(Outcome,levels=c(
    'Depression_Week4','Agency_Week4','Agency_Posttest',
    'Hopelessness_Week4','Hopelessness_Posttest',
    '+_Actions/Thinking_Week4',
    'Change_Expectancy_Week4','Change_Expectancy_Posttest',
    'Readiness_Week4','Readiness_Posttest')
  ))


#at post-test
(ES_across_outcomes_posttest <- ESs_outcomes_all_with_names_long %>%
    filter(str_detect(Outcome,'Posttest')) %>% 
    mutate(Outcome=fct_relabel(
      Outcome,~str_remove(.x,'_Posttest')))  %>%
    mutate(Outcome=fct_relabel(
      Outcome,~str_replace(.x,'_',' '))) %>%
    ggplot(aes(Outcome,term,fill=Effect)) +
    geom_tile() +
    geom_text(aes(label = Effect),  size = 3) +
    scale_fill_gradient2(
      low = "red",
      high = "blue"
    )  +
    ggtitle("Effect Sizes across Outcomes at Post-test") +
    labs(fill = "Cohen's d") +
    theme(axis.text.x = element_text(angle = 20, hjust = 1))+xlab('')+ylab('')
)

#at Week 4
(ES_across_outcomes_wk4 <-ESs_outcomes_all_with_names_long %>% 
    filter(str_detect(Outcome,'Week4'))%>% 
    mutate(Outcome=fct_relabel(
      Outcome,~str_remove(.x,'_Week4')))%>%
    mutate(Outcome=fct_recode(
      Outcome, 'Positive Actions / Thoughts' = '+_Actions/Thinking'
    ))  %>%
    mutate(Outcome=fct_relabel(
      Outcome,~str_replace(.x,'_',' '))) %>% 
    ggplot(aes(Outcome,term,fill=Effect)) +
    geom_tile() +
    geom_text(aes(label = Effect),  size = 3) +
    scale_fill_gradient2(
      low = "red",
      high = "blue"
    )  + ggtitle("Effect Sizes across Outcomes at Week 4") +
    labs(fill = "Cohen's d") +
    theme(axis.text.x = element_text(angle = 20, hjust = 1))+xlab('')+ylab('')
)

#put all wk4 plots together
combined_plot_wk4_effects <- wrap_plots(
  phq9_plot,
  agency_wk4_plot,
  hopelessness_wk4_plot,
  fats_wk4_plot,
  dep_expect_wk4_plot,
  readiness_wk4_plot,ncol=2) 

#put all posttest plots together
combined_plot_posttest_effects <- wrap_plots(
  agency_posttest_plot,
  hopelessness_posttest_plot,
  dep_expect_posttest_plot,
  readiness_posttest_plot,ncol=2)

#put Effect size across outcomes plots together
combined_effects_x_outcomes <- wrap_plots(ES_across_outcomes_posttest, ES_across_outcomes_wk4)

#put all acceptability plots together
combined_acceptability_plots <- wrap_plots(
  ceq_pointplot,aha_pointplot,aha_signif_pointplot,star_pointplot,ncol=2)+
  plot_annotation(
    theme = theme(plot.title = element_text(hjust = 0.5)))

#ggsave("~/megastudy stuff/mega_figures/combined_plot_posttest_effects.png",combined_plot_posttest_effects,height=12,width=12)
#ggsave("~/megastudy stuff/mega_figures/combined_plot_wk4_effects.png",combined_plot_wk4_effects,height=18,width=12)
#ggsave("~/megastudy stuff/mega_figures/combined_effects_x_outcome.png",combined_effects_x_outcomes,height=6,width=12)
#ggsave("~/megastudy stuff/mega_figures/combined_acceptability_plots.png",combined_acceptability_plots,height=12,width=12)
```

#### Correlation between Post-test and Week Four Effects
```{r}
pt_and_wk4_ESs_for_corr <- ESs_outcomes_all_with_names %>% select(term,
  contains('Agency')|contains('Hope')|contains('Expect')|contains('Readiness'))

(agency_es_cor <- with(pt_and_wk4_ESs_for_corr,cor.test(Agency_Posttest,Agency_Week4)))
(agency_es_cor_plot <- pt_and_wk4_ESs_for_corr %>% ggplot(aes(
  Agency_Posttest,Agency_Week4,color=term,label=term))+
  geom_point() +  ggrepel::geom_text_repel(size=2.8)+
    ggtitle(paste('Agency, r =',agency_es_cor$estimate %>% round(2),', p =',agency_es_cor$p.value %>% round(2)))+
  theme(legend.position="none"))

(hope_es_cor <- with(pt_and_wk4_ESs_for_corr,cor.test(Hopelessness_Posttest,Hopelessness_Week4)))
(hope_es_cor_plot <-pt_and_wk4_ESs_for_corr %>% ggplot(aes(
  Hopelessness_Posttest,Hopelessness_Week4,color=term,label=term))+
  geom_point() +  ggrepel::geom_text_repel(size=2.8)+
    ggtitle(paste('Hopelessness, r =',hope_es_cor$estimate %>% round(2),', p =',hope_es_cor$p.value %>% round(2)))+
  theme(legend.position="none"))

(DCES_es_cor <- with(pt_and_wk4_ESs_for_corr,cor.test(Change_Expectancy_Posttest,Change_Expectancy_Week4)))
(DCES_es_cor_plot <-pt_and_wk4_ESs_for_corr %>% ggplot(aes(
  Change_Expectancy_Posttest,Change_Expectancy_Week4,color=term,label=term))+
  geom_point() + ggrepel::geom_text_repel(size=2.8)+
    ggtitle(paste('Depression Change Expectancy, r =',DCES_es_cor$estimate %>% round(2),', p =',DCES_es_cor$p.value %>% round(2)))+
  theme(legend.position="none"))

(readiness_es_cor <- with(pt_and_wk4_ESs_for_corr,cor.test(Readiness_Posttest,Readiness_Week4)))
(readiness_es_cor_plot <-pt_and_wk4_ESs_for_corr %>% ggplot(aes(Readiness_Posttest,Readiness_Week4,color=term,label=term))+
  geom_point() + ggrepel::geom_text_repel(size=2.8)+
    ggtitle(paste('Readiness for Change, r =',readiness_es_cor$estimate %>% round(2),', p =',readiness_es_cor$p.value %>% round(2)))+
  theme(legend.position="none"))

combined_es_cor <- wrap_plots(
  agency_es_cor_plot,hope_es_cor_plot,DCES_es_cor_plot,readiness_es_cor_plot,ncol=2)+
  plot_annotation(
    title="Correlations between Cohen's D at Post-test and Week 4",
    theme = theme(plot.title = element_text(hjust = 0.5)))

#ggsave("~/megastudy stuff/mega_figures/combined_es_cor_plots.png",combined_es_cor,height=12,width=12)


```

### *3*. Change in letting distressing thoughts go (single-item)

```{r,message=F,warning=F}
##A single item from the five-factor mindfulness questionnaire
m1_let_distressing_thoughts_go <- lmer(ffmq_distressing ~ time * group + (1|group/pid),data=mega_alltimes)
add_cohens_d(m1_let_distressing_thoughts_go,'Predicting: change in "I let distressing thoughts go"')
(let_thoughts_go_cohens_d <- mixed_cohens_d(m1_let_distressing_thoughts_go))
(let_thoughts_go_plot <- plot_model(m1_let_distressing_thoughts_go,sort.est = T,terms=wk4_int_terms,auto_label=F,
                                  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Let distressing thoughts go','Baseline to Week 4'))

```

### *4.* How much do you feel the program you completed four weeks ago, in the first part of this study, has helped you?

```{r}
mega_alltimes %>% filter(!is.na(program_helped)) %>% count(program_helped) %>% 
  mutate(n/sum(n))
mega_alltimes %>% filter(!is.na(program_helped),program_helped!=999) %>% 
  summarize(mean(program_helped),sd(program_helped))

(program_helped_hist <- mega_alltimes %>% filter(!is.na(program_helped)) %>%
  mutate(program_helped=ifelse(
    program_helped==999,"Don't remember SSI",program_helped)) %>%
  mutate(program_helped=as.factor(program_helped)) %>% 
  ggplot(aes(program_helped))+geom_bar(fill='brown3')+
    theme(axis.text.x = element_text(angle = 20, hjust = 1))+
    xlab('The SSI 4 Weeks Ago Helped')+
  ggtitle('How Much the SSI Helped','1 (Not at all) - 5 (A great deal)'))

mega_alltimes  %>% filter(group!='Passive Control') %>% 
  filter(!is.na(program_helped),program_helped!=999) %>% 
  group_by(group) %>% 
  summarize(mean_program_helped=mean(program_helped,na.rm=T),std.error=std.error(program_helped,na.rm=T)) %>% 
  ggplot(aes(x = reorder(group, mean_program_helped), y = mean_program_helped)) +
  geom_bar(stat = "identity", fill = "brown2") +
  geom_errorbar(aes(ymin = mean_program_helped - std.error, ymax = mean_program_helped + std.error),
                width = 0.2) + coord_flip() +labs(x = "Term", y = "Estimate \u00B1 SE")

mega_alltimes %>% filter(group!='Passive Control')%>% group_by(group) %>% count(forgot=program_helped==999) %>% filter(forgot==T) %>% ggplot(aes(x=n,y=reorder(group,n)))+geom_bar(stat='identity',fill='gold')+ggtitle('4 weeks later forgot they completed a program')

```


------------------------------------------------------------------------

------------------------------------------------------------------------

### *5.* Having accessed formal support for mental health as a moderator of SSIs' effects on depression from baseline to week 4

```{r,message=F,warning=F}
m1_received_help_moderating_change_in_phq9 <- lmer(
  phq9_score ~ time * received_help + (1|group/pid),
  data=mega_alltimes %>% filter(group!='Passive Control') %>%
    mutate(received_help=ifelse(
      str_detect(help4mental,'not accessed any')&
        str_detect(help4mental,',',negate=T),1,0)))

summary(m1_received_help_moderating_change_in_phq9)
add_cohens_d(m1_received_help_moderating_change_in_phq9)
```

------------------------------------------------------------------------

------------------------------------------------------------------------

### *6.* Addressing Duplicate IP Addresses

### 1. Complete dataset (n=7,733)

First, we'll re-run the central analyses with the complete dataset (n = 7,733 at baseline)

```{r,message=F,warning=F}

#### Primary analysis: condition predicting change in PHQ #### 

m1_phq9_complete <- lmer(phq9_score ~ time * group + (1|pid),
                data=mega_alltimes_complete)

#plot showing Change in Depression over Four Weeks Compared to a passive control
(phq9_plot_complete <- plot_model(m1_phq9_complete,sort.est = T,terms=wk4_int_terms,auto_label=F,
                         show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Depressive Symptoms, Complete Sample','Baseline to Week 4'))

## agency ##
m1_agency_complete <- lmer(agency_score ~ time * group + (1|pid),data=mega_alltimes_complete)

(agency_wk4_plot_complete <- plot_model(
  m1_agency_complete,sort.est = T,terms=wk4_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Agency, Complete Sample','Baseline to Week 4'))

(agency_posttest_plot_complete <- plot_model(
  m1_agency_complete,sort.est = T,terms=posttest_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Agency, Complete Sample','Baseline to Post-test'))

## hopelessness ##
m1_hopeless_complete <- lmer(hopeless_score ~ time * group + (1|pid),data=mega_alltimes_complete)

(hopelessness_wk4_plot <- plot_model(
  m1_hopeless_complete,sort.est = T,terms=wk4_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Hopelessness, Complete Sample','Baseline to Week 4'))

(hopelessness_posttest_plot_complete <- plot_model(
  m1_hopeless_complete,sort.est = T,terms=posttest_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Hopelessness, Complete Sample','Baseline to Post-test'))

## FATS ##
m1_fats_complete <- lmer(fats_score ~ time * group + (1|pid),data=mega_alltimes_complete)

(fats_wk4_plot_complete <- plot_model(m1_fats_complete,sort.est = T,terms=wk4_int_terms,auto_label=F,
                             show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Positive Actions and Thinking, Complete Sample','Baseline to Week 4'))

## Depression expectancies for Change ##
m1_dep_expect_complete <- lmer(dep_expectancies_score ~ time * group + (1|pid),data=mega_alltimes_complete)

dep_expect_wk4_plot_complete <- plot_model(m1_dep_expect_complete,sort.est = T,terms=wk4_int_terms,auto_label=F,
                                   show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Depression Change Expectancies, complete sample','Baseline to Week 4')

(dep_expect_posttest_plot_complete <- plot_model(
  m1_dep_expect_complete,sort.est = T,terms=posttest_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Depression Change Expectancies, Complete Sample','Baseline to Post-test'))

## Depression Readiness for change ##
m1_readiness_complete <- lmer(readiness_score ~ time * group + (1|group/pid),data=mega_alltimes_complete)

(readiness_wk4_plot <- plot_model(m1_readiness_complete,sort.est = T,terms=wk4_int_terms,auto_label=F,
                                  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Readiness for Change, Complete Sample','Baseline to Week 4'))

(readiness_posttest_plot_complete <- plot_model(
  m1_readiness_complete,sort.est = T,terms=posttest_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Readiness for Change, Complete Sample','Baseline to Post-test'))

#Credibility / Expectancy Questionnaire scores
df_for_acceptability_complete <- mega_alltimes_active_first %>% 
  filter(group!='Passive Control') %>% 
  filter(duplicate_ip=='No') %>% 
  mutate(group=droplevels(group))

m1_ceq_complete <- lm(CEQ_score~group,data=df_for_acceptability_complete)
summary(m1_ceq_complete)
rstatix::cohens_d(CEQ_score~group,data=df_for_acceptability_complete) %>% select(-magnitude,-.y.)

acceptability_terms <- m1_ceq_complete %>% tidy() %>% 
  filter(term!='(Intercept)') %>% pull(term)
(ceq_plot <- plot_model(m1_ceq_complete,sort.est = T,terms=acceptability_terms,auto_label=F,
                        show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Credibility and Expected Benefit, Complete Sample','At Post-test'))

(ceq_pointplot_complete <- mega_alltimes_complete %>% filter(time=='Post-test'&group!='Passive Control') %>% 
    group_by(group) %>% 
    summarize(
      mean_CEQ_score = mean(CEQ_score, na.rm = TRUE),
      se = sd(CEQ_score, na.rm = TRUE) / sqrt(n()),
      ci_lower = mean_CEQ_score - 1.96 * se,
      ci_upper = mean_CEQ_score + 1.96 * se) %>% 
    ggplot(aes(x = reorder(group, mean_CEQ_score), y = mean_CEQ_score)) +
    geom_point(size=1.5) +
    geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper),
                  color='green3',width = 0.2,show.legend = F) + 
    coord_flip() +labs(x = "", y = "Mean CEQ Score")+
    ggtitle("Credibility and Expected Benefit, Complete Sample","At Post-test, range 4-56"))


#Frequency of Aha! moments at post-test
m1_aha_complete <- lm(aha~group,data=df_for_acceptability_complete %>% filter(time=='Post-test'))
summary(m1_aha_complete)
rstatix::cohens_d(aha~group,data=df_for_acceptability_complete%>% filter(time=='Post-test')) %>% select(-magnitude,-.y.)
(aha_plot_complete <- plot_model(m1_aha_complete,sort.est = T,terms=acceptability_terms,auto_label=F,
                        show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Probability of an Aha! Moment during SSI, Complete Sample','At Post-test'))


(aha_pointplot <- mega_alltimes_complete %>% 
    filter(time=='Post-test'&group!='Passive Control') %>% 
    group_by(group) %>% 
    summarize(
      aha_prop = mean(aha, na.rm = TRUE),
      n = sum(!is.na(aha)),
      se = sqrt(aha_prop * (1 - aha_prop) / n)
    ) %>% 
    ggplot(aes(x = reorder(group, aha_prop), y = aha_prop)) +
    geom_point(size=1.5) +
    geom_errorbar(aes(ymin = aha_prop - se, ymax = aha_prop +
                        se),
                  color='palegreen2',width = 0.2,show.legend = F) + 
    labs(x = "", y = "Proportion who had an Aha! Moment During SSI")+
    scale_y_continuous(labels = scales::percent)+
    coord_flip() +
    ggtitle("Proportion of Aha! Moments, Complete Sample","At Post-test"))


#significance of aha moments
m1_aha_significance_complete <- lm(aha_significant~group,
                          data=df_for_acceptability_complete %>% 
                                filter(time=='Post-test'))

rstatix::cohens_d(aha_significant~group,data=df_for_acceptability_complete)  %>% select(-magnitude,-.y.)%>%
  kable()
plot_model(
  m1_aha_significance_complete,sort.est = T,terms=acceptability_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Significance of Aha! Moment, Complete Sample','At Post-test, compared to a Behavioral Activation SSI')

(aha_signif_pointplot <- mega_alltimes_complete %>%
    filter(time=='Post-test'&group!='Passive Control') %>% 
  group_by(group) %>% 
  summarize(mean_aha_signif_score=mean(aha_significant,na.rm=T),std.error=std.error(aha_significant,na.rm=T)) %>% 
  ggplot(aes(x = reorder(group, mean_aha_signif_score), y = mean_aha_signif_score)) +
  geom_point(size=1.5) +
  geom_errorbar(aes(ymin = mean_aha_signif_score - std.error, 
                    ymax = mean_aha_signif_score +
                      std.error),
                color='palegreen3',width = 0.2,show.legend = F) + 
    coord_flip() +labs(x = "", y = "Estimate \u00B1 SE")+
    ggtitle("Significance of Aha! Moment, Complete Sample",
            "At Post-test, range 0-5"))

#star ratings
m1_star_complete <- lm(star_rating_1~group,data=df_for_acceptability_complete %>% 
                filter(time=='Post-test'))
 rstatix::cohens_d(star_rating_1~group,data=df_for_acceptability_complete %>% filter(time=='Post-test')) %>% select(-magnitude,-.y.) %>% kable
(star_plot_complete <- plot_model(
  m1_star_complete,sort.est = T,terms=acceptability_terms,
  auto_label=F,show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Star Rating, Complete Sample','At Post-test, compared to a Behavioral Activation SSI'))

(star_pointplot_complete <- mega_alltimes_complete %>%
    filter(time=='Post-test'&group!='Passive Control') %>% 
  group_by(group) %>% 
  summarize(mean_star_score=mean(star_rating_1,na.rm=T),
            std.error=std.error(star_rating_1,na.rm=T)) %>% 
  ggplot(aes(x = reorder(group, mean_star_score), y = mean_star_score)) +
  geom_point(size=1.5) +
  geom_errorbar(aes(ymin = mean_star_score - std.error, ymax = mean_star_score +
                      std.error),
                color='palegreen1',width = 0.2,show.legend = F) + 
    coord_flip() +labs(x = "", y = "Estimate \u00B1 SE")+
    ggtitle("Star Rating, Complete Sample","At Post-test, range 1-5"))

```


### 2. Remove all Duplicate IPs

Next, we'll remove all responses from anyone who had a duplicate IP. This leaves 91% of the complete data

```{r,message=F,warning=F}
mega_alltimes_nodups <- mega_alltimes_complete %>% filter(duplicate_ip=='No')


#### Primary analysis: condition predicting change in PHQ #### 

m1_phq9_nodups <- lmer(phq9_score ~ time * group + (1|pid),
                data=mega_alltimes_nodups)

#plot showing Change in Depression over Four Weeks Compared to a passive control
(phq9_plot_nodups <- plot_model(m1_phq9_nodups,sort.est = T,terms=wk4_int_terms,auto_label=F,
                         show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Depressive Symptoms, no Duplicate IPs','Baseline to Week 4'))

## agency ##
m1_agency_nodups <- lmer(agency_score ~ time * group + (1|pid),data=mega_alltimes_nodups)

(agency_wk4_plot_nodups <- plot_model(
  m1_agency_nodups,sort.est = T,terms=wk4_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Agency, no Duplicate IPs','Baseline to Week 4'))

(agency_posttest_plot_nodups <- plot_model(
  m1_agency_nodups,sort.est = T,terms=posttest_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Agency, no Duplicate IPs','Baseline to Post-test'))

## hopelessness ##
m1_hopeless_nodups <- lmer(hopeless_score ~ time * group + (1|pid),data=mega_alltimes_nodups)

(hopelessness_wk4_plot <- plot_model(m1_hopeless_nodups,sort.est = T,terms=wk4_int_terms,auto_label=F,
                                     show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Hopelessness, no Duplicate IPs','Baseline to Week 4'))

(hopelessness_posttest_plot_nodups <- plot_model(
  m1_hopeless_nodups,sort.est = T,terms=posttest_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Hopelessness, no Duplicate IPs','Baseline to Post-test'))

## FATS ##
m1_fats_nodups <- lmer(fats_score ~ time * group + (1|pid),data=mega_alltimes_nodups)

(fats_wk4_plot_nodups <- plot_model(m1_fats_nodups,sort.est = T,terms=wk4_int_terms,auto_label=F,
                             show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Positive Actions and Thinking, no Duplicate IPs','Baseline to Week 4'))

## Depression expectancies for Change ##
m1_dep_expect_nodups <- lmer(dep_expectancies_score ~ time * group + (1|pid),data=mega_alltimes_nodups)

dep_expect_wk4_plot_nodups <- plot_model(m1_dep_expect_nodups,sort.est = T,terms=wk4_int_terms,auto_label=F,
                                   show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Depression Change Expectancies, no Duplicate IPs','Baseline to Week 4')

(dep_expect_posttest_plot_nodups <- plot_model(
  m1_dep_expect_nodups,sort.est = T,terms=posttest_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Depression Change Expectancies, no Duplicate IPs',
            'Baseline to Post-test'))

## Depression Readiness for change ##
m1_readiness_nodups <- lmer(readiness_score ~ time * group + (1|group/pid),data=mega_alltimes_nodups)
add_cohens_d(m1_readiness_nodups)
(readiness_wk4_plot <- plot_model(m1_readiness_nodups,sort.est = T,terms=wk4_int_terms,auto_label=F,
                                  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Readiness for Change, no Duplicate IPs','Baseline to Week 4'))

(readiness_posttest_plot_nodups <- plot_model(
  m1_readiness_nodups,sort.est = T,terms=posttest_int_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
    ggtitle('Readiness for Change, no Duplicate IPs','Baseline to Post-test'))

#Credibility / Expectancy Questionnaire scores
df_for_acceptability_nodups <- mega_alltimes_active_first %>% 
  filter(group!='Passive Control') %>% 
  filter(duplicate_ip=='No') %>% 
  mutate(group=droplevels(group))

m1_ceq_nodups <- lm(CEQ_score~group,data=df_for_acceptability_nodups)


acceptability_terms <- m1_ceq_nodups %>% tidy() %>% 
  filter(term!='(Intercept)') %>% pull(term)
(ceq_plot <- plot_model(m1_ceq_nodups,sort.est = T,terms=acceptability_terms,auto_label=F,
                        show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Credibility and Expected Benefit, no Duplicate IPs','At Post-test'))

(ceq_pointplot_nodups <- mega_alltimes_nodups %>% filter(time=='Post-test'&group!='Passive Control') %>% 
    group_by(group) %>% 
    summarize(mean_CEQ_score=mean(CEQ_score,na.rm=T),std.error=std.error(CEQ_score,na.rm=T)) %>% 
    ggplot(aes(x = reorder(group, mean_CEQ_score), y = mean_CEQ_score)) +
    geom_point(size=1.5) +
    geom_errorbar(aes(ymin = mean_CEQ_score - std.error, ymax = mean_CEQ_score +
                        std.error),
                  color='green3',width = 0.2,show.legend = F) + 
    coord_flip() +labs(x = "", y = "Estimate ± SE")+
    ggtitle("Credibility and Expected Benefit, no Duplicate IPs",
            "At Post-test, range 4-56"))


#Frequency of Aha! moments at post-test
m1_aha_nodups <- lm(aha~group,data=df_for_acceptability_nodups %>% filter(time=='Post-test'))

(aha_plot_nodups <- plot_model(m1_aha_nodups,sort.est = T,terms=acceptability_terms,auto_label=F,
                        show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Probability of an Aha! Moment during SSI, no Duplicate IPs','At Post-test'))


(aha_pointplot <- mega_alltimes_nodups %>% 
    filter(time=='Post-test'&group!='Passive Control') %>% 
    group_by(group) %>% 
    summarize(
      aha_prop = mean(aha, na.rm = TRUE),
      n = sum(!is.na(aha)),
      se = sqrt(aha_prop * (1 - aha_prop) / n)
    ) %>% 
    ggplot(aes(x = reorder(group, aha_prop), y = aha_prop)) +
    geom_point(size=1.5) +
    geom_errorbar(aes(ymin = aha_prop - se, ymax = aha_prop +
                        se),
                  color='palegreen2',width = 0.2,show.legend = F) + 
    labs(x = "", y = "Proportion who had an Aha! Moment During SSI")+
    scale_y_continuous(labels = scales::percent)+
    coord_flip() +
    ggtitle("Proportion of Aha! Moments, no duplicate IPs","At Post-test"))


#significance of aha moments
m1_aha_significance_nodups <- lm(aha_significant~group,
                          data=df_for_acceptability_nodups %>% 
                                filter(time=='Post-test'))

rstatix::cohens_d(aha_significant~group,data=df_for_acceptability_nodups) %>%select(-magnitude,-.y.) %>%
  kable()
plot_model(
  m1_aha_significance_nodups,sort.est = T,terms=acceptability_terms,auto_label=F,
  show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Significance of Aha! Moment, complete dataset','At Post-test, compared to a Behavioral Activation SSI')

(aha_signif_pointplot <- mega_alltimes_nodups %>%
    filter(time=='Post-test'&group!='Passive Control') %>% 
  group_by(group) %>% 
  summarize(mean_aha_signif_score=mean(aha_significant,na.rm=T),
            std.error=std.error(aha_significant,na.rm=T)) %>% 
  ggplot(aes(x = reorder(group, mean_aha_signif_score), y = mean_aha_signif_score)) +
  geom_point(size=1.5) +
  geom_errorbar(aes(ymin = mean_aha_signif_score - std.error, 
                    ymax = mean_aha_signif_score +
                      std.error),
                color='palegreen3',width = 0.2,show.legend = F) + 
    coord_flip() +labs(x = "", y = "Estimate \u00B1 SE")+
    ggtitle("Significance of Aha! Moment, no Duplicate IPs","At Post-test, range 0-5"))

#star ratings
m1_star_nodups <- lm(star_rating_1~group,data=df_for_acceptability_nodups %>% 
                filter(time=='Post-test'))
 rstatix::cohens_d(star_rating_1~group,data=df_for_acceptability_nodups  %>%filter(time=='Post-test'))%>% select(-magnitude,-.y.) %>% kable
(star_plot_nodups <- plot_model(
  m1_star_nodups,sort.est = T,terms=acceptability_terms,
  auto_label=F,show.values = T, value.offset = .4,value.size = 3.2)+
    scale_x_discrete(labels = function(x) gsub(".*group", "", x)) +
    ggtitle('Star Rating no Duplicate IPs','At Post-test, compared to a Behavioral Activation SSI'))

(star_pointplot_nodups <- mega_alltimes_nodups %>%
    filter(time=='Post-test'&group!='Passive Control') %>% 
  group_by(group) %>% 
  summarize(mean_star_score=mean(star_rating_1,na.rm=T),
            std.error=std.error(star_rating_1,na.rm=T)) %>% 
  ggplot(aes(x = reorder(group, mean_star_score), y = mean_star_score)) +
  geom_point(size=1.5) +
  geom_errorbar(aes(ymin = mean_star_score - std.error, ymax = mean_star_score +
                      std.error),
                color='palegreen1',width = 0.2,show.legend = F) + 
    coord_flip() +labs(x = "", y = "Estimate \u00B1 SE")+
    ggtitle("Star Rating no Duplicate IPs","At Post-test, range 1-5"))
```


------------------------------------------------------------------------

------------------------------------------------------------------------

### *7.* Accounting for differential attrition with inverse probability weighting (pre-registered)

First, checking for differential attrition at post-test and 4-week follow-up by experimental condition

```{r,include=F}
#STEP 7. ATTRITION ANALYSES - Test if we need to adjust for potential attrition bias
ipwForPostTest <- F
ipwForFollowUp <- F

#attrited before post-test measures
table(mega_alltimes$attrited_b4_post_test)

#attrited before follow-up measures
table(mega_alltimes$attrited_b4_followup)

# Test 1 - Any Differential Attrition Rate?

# after completing post-test measures?
mega_alltimes %>%
  group_by(group) %>%
  summarise(AttritionRate  = mean(attrited_b4_post_test) * 100) %>%
  arrange(desc(AttritionRate))

#differential attrition at post-test
attTest1 <- lm(attrited_b4_post_test ~ 1 + group, data = mega_alltimes %>% 
                 filter(time=='Baseline'))
mainEffects <- grep('group', names(attTest1$coefficients), value = TRUE)
attTest1_R <- car::linearHypothesis(attTest1, mainEffects)
attTest1_R
if (attTest1_R$`Pr(>F)`[2] < .05) {
  ipwForPostTest <- TRUE
}
rm(attTest1, attTest1_R, mainEffects)


# differential attrition at four-week follow-up?
mega_alltimes %>%
  group_by(group) %>%
  summarise(AttritionRate  = mean(attrited_b4_followup,na.rm=T) * 100) %>%
  arrange(desc(AttritionRate))
attTest1_wk4 <- lm(attrited_b4_followup ~ 1 + group, data = mega_alltimes %>% 
                 filter(time=='Baseline'))
mainEffects_wk4 <- grep('group', names(attTest1_wk4$coefficients), value = TRUE)
attTest1_R_wk4 <- car::linearHypothesis(attTest1_wk4, mainEffects_wk4) 
attTest1_R_wk4
if (attTest1_R_wk4$`Pr(>F)`[2] < .05) {
  ipwForFollowUp <- TRUE
}
rm(attTest1_wk4, attTest1_R_wk4, mainEffects_wk4)

#Any Differential Attrition by Covariates?
#couldn't add many covariates because they just end up being aliased.

# post-test
attTest2 <- lm(attrited_b4_post_test ~ 1 + group:(
    dep_expect_bl+readiness_bl+ agency_bl+hopeless_bl), 
  data = mega_alltimes %>% filter(time=='Baseline')) 

#dealing with aliased coefficients
# Get all terms from the original model
all_terms <- attr(terms(attTest2), "term.labels")

ld.vars <- attributes(alias(attTest2)$Complete)$dimnames[[1]]

# Remove linearly dependent terms
valid_terms <- setdiff(all_terms, ld.vars)

# Reconstruct formula
rhs <- paste(valid_terms, collapse = " + ")
formula.new <- as.formula(paste("attrited_b4_post_test ~", rhs))

fit.new <- lm(formula.new, 
              data = mega_alltimes %>% 
                filter(time=='Baseline') )
interactions <- grep(':', names(fit.new$coefficients), value = TRUE)

attTest2_R <- car::linearHypothesis(fit.new, interactions)
attTest2_R
if (attTest2_R$`Pr(>F)`[2] < .05) {
  ipwForPostTest <- TRUE
}
rm(attTest2, attTest2_R, interactions)

# follow-up
attTest2 <- lm(attrited_b4_followup ~ 1 + group:(
    dep_expect_bl+readiness_bl+ agency_bl+hopeless_bl), 
  data = mega_alltimes %>% filter(time=='Baseline')) 

#dealing with aliased coefficients
# Get all terms from the original model
all_terms <- attr(terms(attTest2), "term.labels")

ld.vars <- attributes(alias(attTest2)$Complete)$dimnames[[1]]

# Remove linearly dependent terms
valid_terms <- setdiff(all_terms, ld.vars)

# Reconstruct formula
rhs <- paste(valid_terms, collapse = " + ")
formula.new <- as.formula(paste("attrited_b4_followup ~", rhs))

fit.new <- lm(formula.new, 
              data = mega_alltimes %>% 
                filter(time=='Baseline') )
interactions <- grep(':', names(fit.new$coefficients), value = TRUE)

              
attTest2_R <- car::linearHypothesis(attTest2, interactions, white.adjust = "hc1")

attTest2_R
if (attTest2_R$`Pr(>F)`[2] < .05) {
  ipwForFollowUp <- TRUE
}
rm(attTest2, attTest2_R, interactions)

###### ATTRITION ANALYSES - Get IPW weights ######
mega_alltimes <- mega_alltimes %>%
  mutate(attrited_b4_post_test = as.factor(attrited_b4_post_test),
         attrited_b4_followup = as.factor(attrited_b4_followup))

#i guess random forest only likes the names in this format
mega_renamed <- mega_alltimes %>%  
  mutate(relationship=`Relationship/Marital Status`,
         political_party=`Political Party`,
         employment_status=`Employment Status`,
         household_income=`Household Income`,
         social_ladder=`Social ladder`)

# post-test weights
if (ipwForPostTest == TRUE) {
  rf_post_test <-randomForest(
    attrited_b4_post_test ~ group + Age + Gender + Race + 
      Education + relationship + political_party + employment_status + 
      household_income + social_ladder + disability+
      dep_expect_bl+readiness_bl+ agency_bl+hopeless_bl ,
    data = mega_renamed %>% filter(time=='Baseline'),
    importance = TRUE, ntree = 200, proximity = TRUE,na.action = na.exclude)
  
  pre_post_test <- predict(
    rf_post_test, newdata = mega_renamed %>% filter(time=='Baseline'),
    type='prob') # probability of attrition
  
  cbind(pre_post_test, mega_alltimes$attrited_b4_post_test[mega_alltimes$time=='Baseline']) # double check the model
  
  mega_alltimes$weights_post_test[mega_alltimes$time=='Baseline'] <- 
    1/pre_post_test[,'0'] # take inverse of probability of leaving as weights
  
  mega_alltimes <- mega_alltimes %>% group_by(pid) %>% 
    fill(weights_post_test,.direction='downup') %>% ungroup() %>% 
    mutate(weights_post_test=ifelse(
      attrited_b4_post_test==1,1,weights_post_test)) # if they attrited the weight is 1
  
  rm(rf_post_test, pre_post_test)
}


# follow-up weights
if (ipwForFollowUp == TRUE) {
  rf_followup <-randomForest(
    attrited_b4_followup ~ Age + Gender + Race + 
      Education + relationship + political_party + employment_status + 
      household_income + social_ladder + disability+
      dep_expect_bl+readiness_bl+ agency_bl+hopeless_bl ,
    data = mega_renamed %>% filter(time=='Baseline'),
    importance = TRUE, ntree = 200, proximity = TRUE,na.action = na.exclude)
  
  pre_followup_test <- predict(
    rf_followup, newdata = mega_renamed %>% filter(time=='Baseline'),
    type='prob') # probability of attrition

    cbind(pre_followup_test, mega_alltimes$attrited_b4_followup[mega_alltimes$time=='Baseline']) # double check the model
 
     mega_alltimes$weights_followup[
    mega_alltimes$time=='Baseline'] <- 
    1/pre_followup_test[,'0'] # take inverse of probability of leaving as weights
  
     mega_alltimes <- mega_alltimes %>% group_by(pid) %>% 
    fill(weights_followup,.direction='downup') %>% ungroup() %>% 
    mutate(weights_followup=ifelse(
      attrited_b4_followup==1,1,weights_followup))
 
      rm(rf_followup, pre_followup_test)
}
#not exactly sure which to use for 3-timepoint measures, but going with follow-up

#I think they should be pretty well correlated.
#cor.test(mega_alltimes$weights_followup[mega_alltimes$attrited_b4_followup==0],
#         mega_alltimes$weights_post_test[mega_alltimes$attrited_b4_followup==0])
```

-   Re-run analyses with inverse-probability weighting to account for differential dropout

```{r,message=F,warning=F}

ipw_m1_phq9 <- lmer(phq9_score ~ time * group  +(1|pid),
                    data=mega_alltimes,weights = weights_followup)
add_cohens_d(ipw_m1_phq9,'Predicting: change in depressive symptoms')

(phq9_plot_ipw <- plot_model(ipw_m1_phq9,sort.est = T,terms=wk4_int_terms,auto_label=F,
                            show.values = T, value.offset = .4,value.size = 3.2)+
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
  ggtitle('Change in Depressive Symptoms with IPW'))

phq9_plot +ggtitle('Original depresssion plot, to compare')


ipw_m1_phq9_activefirst <- lmer(
  phq9_score ~ time * group + (1|pid),data=mega_alltimes_active_first %>% 

        left_join(mega_alltimes %>% select(pid,time,weights_followup)),
  weights= weights_followup)
add_cohens_d(ipw_m1_phq9_activefirst)

ipw_m1_phq9_w_demcovars <- lmer(phq9_score ~ time * group + Age + Gender + Race + 
                                  Education + `Relationship/Marital Status` + 
                                  `Political Party` + `Employment Status` + 
                                  `Household Income` + `Social ladder` + 
                                  disability +(1|pid),
                                data=mega_alltimes,weights = weights_followup)
add_cohens_d(ipw_m1_phq9_w_demcovars,'Predicting: change in depressive symptoms')


ipw_m1_phq9_w_all_covars <- lmer(phq9_score ~ time * group + Age + Gender + Race + 
                                   Education + `Relationship/Marital Status` + 
                                   `Political Party` + `Employment Status` + 
                                   `Household Income` + `Social ladder` + 
                                   disability + dep_expect_bl+readiness_bl+
                                   agency_bl+hopeless_bl+ fats_bl + (1|pid),
                                 data=mega_alltimes,weights = weights_followup)
add_cohens_d(ipw_m1_phq9_w_all_covars,'Predicting: change in depressive symptoms')


## agency ##

ipw_m1_agency <- lmer(agency_score ~ time * group + (1|pid),data=mega_alltimes,
                      weights=weights_followup)
add_cohens_d(ipw_m1_agency,'Predicting: change in agency')

(agency_posttest_plot_ipw <- plot_model(ipw_m1_agency,sort.est = T,terms=posttest_int_terms,auto_label=F,
                            show.values = T, value.offset = .4,value.size = 3.2)+
  scale_x_discrete(labels = function(x) gsub(".*:group", "", x)) +
  ggtitle('Change in Agency with IPW'))

agency_posttest_plot+ggtitle('Original agency plot, to compare')


## hopelessness ##

ipw_m1_hopeless <- lmer(hopeless_score ~ time * group + (1|pid),data=mega_alltimes,
                        weights=weights_followup)
add_cohens_d(ipw_m1_hopeless,'Predicting: change in hopelessness')

## FATS ##

ipw_m1_fats <- lmer(fats_score ~ time * group + (1|pid),data=mega_alltimes,
                    weights=weights_followup)
add_cohens_d(ipw_m1_fats,'Predicting: change in thoughts and actions')

## Depression Change Expectancies ##

ipw_m1_dep_expect <- lmer(dep_expectancies_score ~ time * group + (1|pid),data=mega_alltimes,
                          weights=weights_followup)
add_cohens_d(ipw_m1_dep_expect,'Predicting: change in depression change expectancies')

#Acceptability across SSIs

df_for_acceptability_with_weights <- df_for_acceptability%>% 
  left_join(mega_alltimes %>% select(pid,time,weights_followup,weights_post_test))

#Credibility / Expectancy Questionnaire scores
# i computed this by summing the likert scales and the percentages / 10
ipw_m1_ceq <- lm(CEQ_score~group,data= df_for_acceptability_with_weights,
                 weights=df_for_acceptability_with_weights$weights_post_test)
summary(ipw_m1_ceq) 

#Star rating
ipw_m1_star <- lm(star_rating_1~group,data=df_for_acceptability_with_weights %>% 
                    filter(time=='Post-test'),
                  weights=df_for_acceptability_with_weights %>% 
                    filter(time=='Post-test') %>% pull(weights_post_test))
summary(ipw_m1_star)

ipw_m1_star_wk4 <- lm(star_rating_1~group,data=df_for_acceptability_with_weights %>% 
                        filter(time=='Week 4'),
                      weights=df_for_acceptability_with_weights %>% 
                        filter(time=='Week 4') %>% pull(weights_post_test))
summary(ipw_m1_star_wk4)

#Frequency of Aha! moments

ipw_m1_aha <- lm(aha~group,data=df_for_acceptability_with_weights %>% 
                    filter(time=='Post-test'),
                  weights=df_for_acceptability_with_weights %>% 
                    filter(time=='Post-test') %>% pull(weights_post_test))
summary(ipw_m1_aha)

ipw_m1_aha_wk4 <- lm(aha_wk4~group,data=df_for_acceptability_with_weights %>% 
                        filter(time=='Week 4'),
                      weights=df_for_acceptability_with_weights %>% 
                        filter(time=='Week 4') %>% pull(weights_post_test))
summary(ipw_m1_aha_wk4)

#significance of Aha! moments
ipw_m1_aha_significance <- lm(aha_significant~group,data=df_for_acceptability_with_weights %>% 
                                filter(time=='Post-test'),
                              weights=df_for_acceptability_with_weights %>% 
                                filter(time=='Post-test') %>% pull(weights_followup))
summary(ipw_m1_aha_significance)

ipw_m1_aha_significance_wk4 <- lm(aha_significant~group,data=df_for_acceptability_with_weights %>% 
                                    filter(time=='Week 4'),
                                  weights=df_for_acceptability_with_weights %>% 
                                    filter(time=='Week 4') %>% pull(weights_followup))
summary(ipw_m1_aha_significance_wk4)


ipw_m1_aha_significance <- lm(aha_significant~group,data=df_for_acceptability_with_weights %>% 
                   filter(time=='Post-test'),
                 weights=df_for_acceptability_with_weights %>% 
                   filter(time=='Post-test') %>% pull(weights_post_test))
summary(ipw_m1_aha_significance)

ipw_m1_aha_significance_wk4 <- lm(aha_significant~group,data=df_for_acceptability_with_weights %>% 
                       filter(time=='Week 4'),
                     weights=df_for_acceptability_with_weights %>% 
                       filter(time=='Week 4') %>% pull(weights_post_test))
summary(ipw_m1_aha_significance_wk4)

#how much Aha! changed behavior at 4-week follow-up
ipw_m1_aha_howchanged <- lm(aha_howchanged~group,data=df_for_acceptability_with_weights,
                            weights=df_for_acceptability_with_weights %>% 
                              pull(weights_post_test))
summary(ipw_m1_aha_howchanged)

#expectancies at baseline predicting PHQ-9
ipw_m1_expect_on_phq9 <- lmer(
  phq9_score ~ time * dep_expect_bl + (1|group/pid),
  data=mega_alltimes,
  weights=weights_followup)
add_cohens_d(ipw_m1_expect_on_phq9,'Predicting: change in depressive symptoms')

#passive vs all
mega_alltimes <- mega_alltimes %>% mutate(passive_vs_all=ifelse(
  group=='Passive Control','Passive Control','active')%>% factor %>% fct_rev())
ipw_m1_phq9_all_vs_passive <- lmer(phq9_score ~ time * passive_vs_all + (1|group/pid),data=mega_alltimes,
                                   weights=weights_followup)
summary(ipw_m1_phq9_all_vs_passive)
add_cohens_d(ipw_m1_phq9_all_vs_passive,'Predicting: change in depressive symptoms')
```


------------------------------------------------------------------------

------------------------------------------------------------------------

# Participant Characteristics

```{r,message=F,warning=F}

##Summary statistics for outcomes at their first measurement timepoint (baseline or post-test)##
mega_alltimes %>% filter(time=='Baseline') %>% #baseline ones
  select(ends_with('score'),-CEQ_score) %>% 
  get_summary_stats(type = "common") %>% bind_rows(
    mega_alltimes %>% filter(time=='Post-test') %>% #post-test ones
      select(CEQ_score,star_rating_1,aha,aha_significant,aha_howchanged) %>% 
      get_summary_stats(type = "common"))

load('demogs_before_impute.Rdata')
demog_notimputed_filtered <- demog_notimputed %>% 
  left_join(mega_alltimes_complete %>%filter(time=='Baseline') %>%  select(pid,ip_gender_race_age_same_at_bl,IP_address_id)) %>%
  group_by(IP_address_id, time) %>%
  arrange(StartDate, .by_group = TRUE) %>%
  filter( ip_gender_race_age_same_at_bl==0 |
      (ip_gender_race_age_same_at_bl==1 & row_number() == 1 )) %>%
  ungroup()  
  
  
#Baseline demographics of sample
#need to use the data before we imputed missing values
chars <- demog_notimputed_filtered %>%
   mutate(Race=ifelse(str_detect(Race,'not listed'),'Other',Race)) %>% 
   add_count(Race,name='race_count') %>%
   add_count(Ethnicity,name='ethnicity_count') %>%
   add_count(Gender,name='gender_count') %>% 
   add_count(`Political Party`,name='political_count') %>% 
   add_count(`Relationship/Marital Status`,name='relation_count') %>% 
   add_count(`Employment Status`,name='employment_count') %>% 
   add_count(disability,name='disability_count') %>% 
   mutate(Race=fct_relevel(fct_reorder(factor(Race),race_count,.desc=T),
                           "Other", after = Inf),
          Ethnicity=fct_relevel(fct_reorder(factor(Ethnicity),
                                            ethnicity_count,.desc=T),
                           "Other", after = Inf),
          Gender=fct_reorder(factor(Gender),gender_count,.desc=T),
          `Political Party`=fct_reorder(
            factor(`Political Party`),political_count,.desc=T),
          `Relationship/Marital Status`=fct_reorder(
            factor(`Relationship/Marital Status`),relation_count,.desc=T),
          `Employment Status`=fct_reorder(
            factor(`Employment Status`),employment_count,.desc=T),
          disability=fct_reorder(factor(disability),
                                 disability_count,.desc=T)) %>% 
   select(Age,Race,Ethnicity, Gender,
          'Disability or Health Condition'=disability,Education,
          `Political Party`,`Relationship/Marital Status`,`Employment Status`,
          `Household Income`,'Social Ladder (1-10)'=`Social ladder`) %>% 
   tbl_summary(statistic = list(
     'Age' ~ "{mean} ({sd})",'Social Ladder (1-10)' ~ "{mean} ({sd})"),
     digits = list('Age' ~ c(2,2),'Social Ladder (1-10)' ~c(2,2)),missing='no')

times_chars <- chars %>%
  as_gt() %>% 
  gt::tab_options(table.font.names = "Times New Roman")
#gt::gtsave(times_chars,"baseline_characteristics_mega_final.docx")

help_received_past_month_df <- mega_alltimes %>%  filter(!is.na(help4mental)) %>%
  filter(time=='Baseline') %>% 
  mutate(item = str_split(help4mental, ",")) %>%
  unnest(item) %>%
  count(item, sort = TRUE) %>% 
  mutate(`Percent Who Endorsed`=
           (n/nrow(mega_alltimes %>% filter(time=='Baseline'))*100)) 

#received any formal support 
mega_alltimes %>%  
    filter(time=='Baseline') %>% 
  count(str_detect(help4mental,'not accessed')&str_detect(help4mental,',',negate=T)) %>% mutate(n/sum(n)) # this estimate of "no help" is better than the other one

ggplot(help_received_past_month_df, 
       aes(x = reorder(item, n), y = `Percent Who Endorsed`)) +
  geom_col(fill = "navy") +coord_flip() + 
  scale_y_continuous(labels=scales::percent)  +
  labs( title = "Help Used in the Past Month", x = "" )

#creating plots of well-being outcomes at baseline

# Histogram plotting function with auto bins and smart x limits
plot_histogram_with_median <- function(data, variable, label, fill_color) {
    x <- data[[variable]]
    ggplot(data, aes_string(variable)) +
        geom_histogram(fill = fill_color, alpha = 0.7, color = "white", bins = nclass.FD(x), na.rm = TRUE) +
        xlab(label) + 
        ylab(NULL) +
        geom_vline(xintercept = median(x, na.rm = TRUE), 
                   color = 'forestgreen', linewidth = 1, linetype = 'twodash', alpha = 0.7) +
        scale_x_continuous(limits = range(x, na.rm = TRUE) + c(-0.5, 0.5)) +
        theme(axis.title.x = element_text(margin = grid::unit(c(10, 0, 0, 0), "pt")))
}

# Define the outcome variables and their labels
outcomes <- list(
    list(variable = 'phq9_score', label = 'Depressive Symptoms'),
    list(variable = 'agency_score', label = 'Agency'),
    list(variable = 'hopeless_score', label = 'Hopelessness'),
    list(variable = 'fats_score', label = 'Positive Actions and Thinking'),
    list(variable = 'dep_expectancies_score', label = 'Depression Change Expectancies'),
    list(variable = 'readiness_score', label = 'Readiness for Change')
)

# Filter for baseline data
baseline_data <- mega_alltimes%>% filter(time == 'Baseline')  # update if needed

# Use a colorful but distinct palette
colors <- viridis(length(outcomes), option = "plasma")

# Generate the plots
plots <- lapply(seq_along(outcomes), function(i) {
    outcome <- outcomes[[i]]
    plot_histogram_with_median(baseline_data, outcome$variable, outcome$label, colors[i])
})

# Combine and title the plots
combined_plot <- wrap_plots(plots) + 
    plot_annotation(title = "Outcomes at Baseline",
                    theme = theme(plot.title = element_text(hjust = 0.5, size = 16)))

# Show or save
combined_plot
# ggsave("~/megastudy stuff/mega_figures/outcome_dists.png", combined_plot, width = 12, height = 8)

#Correlations between outcomes at baseline
corr <- round(cor(mega_alltimes %>% filter(time=='Baseline') %>% select(
  phq9_score,agency_score,hopeless_score,fats_score,
  dep_expectancies_score,readiness_score)), 2)
ggcorrplot(corr, hc.order = TRUE, type = "lower",
           outline.col = "white",lab=T,title = 'Correlations at Baseline')

```



Thanks for reading :)

