Bayes Factor Analyses

Load packages

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(here) 

## here() starts at /Users/sarahdaniels/Desktop/Total_Analyses_5

library(janitor) 

## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(haven) 
library(naniar) 
library(ggpubr) 
library(report) 
library(ggplot2)
library(reshape2)

## 
## Attaching package: 'reshape2'
## 
## The following object is masked from 'package:tidyr':
## 
##     smiths

library(lme4)

## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

library(sjPlot)
library(parameters)
library(mediation)

## Loading required package: MASS
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Loading required package: mvtnorm
## Loading required package: sandwich
## mediation: Causal Mediation Analysis
## Version: 4.5.0

library(lavaan)

## This is lavaan 0.6-15
## lavaan is FREE software! Please report any bugs.

library(lmerTest)

## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step

library(modEvA)
library(rsconnect)
library(effectsize)

## 
## Attaching package: 'effectsize'
## 
## The following object is masked from 'package:mvtnorm':
## 
##     standardize

BF_data <- read.csv("Full_data_all.csv")

IUS

Post

IUS_BP <- BF_data %>% 
  dplyr::select("ID", "Group", "A_PRE_IUS_total", "B_POST_IUS_total")
IUS_BP_long <- IUS_BP %>%
  pivot_longer(cols = c(A_PRE_IUS_total, B_POST_IUS_total),
               names_to = "Time",
               values_to = "IUS_Score")

full_lmer_IUSbp <- lmer(IUS_Score ~ Group * Time + (1|ID), data = IUS_BP_long, REML = TRUE)
null_lmer_IUSbp <- update(full_lmer_IUSbp, formula = ~ . -Time:Group) # null means no interaction effect
BF_BIC_IUSbp <- exp((BIC(null_lmer_IUSbp) - BIC(full_lmer_IUSbp))/2)
BF_BIC_IUSbp # for the interaction

## [1] 2612.411

M2_lmer_IUSbp <- lmer(IUS_Score ~ Time + Group + (1|ID), data = IUS_BP_long, REML = TRUE)
null_lmer_IUSbp <- update(M2_lmer_IUSbp, formula = ~ . -Group) # null means no group effect
BF_BIC_IUSbp <- exp((BIC(null_lmer_IUSbp) - BIC(M2_lmer_IUSbp))/2)
BF_BIC_IUSbp # for group

## [1] 0.02741644

M3_lmer_IUSbp <- lmer(IUS_Score ~ Time + Group + (1|ID), data = IUS_BP_long, REML = TRUE)
null_lmer_IUSbp <- update(M3_lmer_IUSbp, formula = ~ . -Time) # null means no time effect
BF_BIC_IUSbp <- exp((BIC(null_lmer_IUSbp) - BIC(M3_lmer_IUSbp))/2)
BF_BIC_IUSbp # for the time

## [1] 2.401982e+14

IUS_BP_long_I <- IUS_BP %>%
  filter(Group == "C_Intervention") %>% 
  pivot_longer(cols = c(A_PRE_IUS_total, B_POST_IUS_total),
               names_to = "Time",
               values_to = "IUS_Score")
IUS_BP_long_C <- IUS_BP %>%
  filter(Group == "B_Controls") %>% 
  pivot_longer(cols = c(A_PRE_IUS_total, B_POST_IUS_total),
               names_to = "Time",
               values_to = "IUS_Score")
IUS_BP_long_EC <- IUS_BP %>%
  filter(Group == "A_ECs") %>% 
  pivot_longer(cols = c(A_PRE_IUS_total, B_POST_IUS_total),
               names_to = "Time",
               values_to = "IUS_Score")

full_lmer_IUSbp_I <- lmer(IUS_Score ~ Time + (1|ID), data = IUS_BP_long_I, REML = TRUE)
null_lmer_IUSbp_I <- update(full_lmer_IUSbp_I, formula = ~ . -Time)
BF_BIC_IUSbp_I <- exp((BIC(null_lmer_IUSbp_I) - BIC(full_lmer_IUSbp_I))/2)
BF_BIC_IUSbp_I

## [1] 2825592191

M2_lmer_IUSbp_C <- lmer(IUS_Score ~ Time + (1|ID), data = IUS_BP_long_C, REML = TRUE)
null_lmer_IUSbp_C <- update(M2_lmer_IUSbp_C, formula = ~ . -Time)
BF_BIC_IUSbp_C <- exp((BIC(null_lmer_IUSbp_C) - BIC(M2_lmer_IUSbp_C))/2)
BF_BIC_IUSbp_C

## [1] 3029998

M3_lmer_IUSbp_EC <- lmer(IUS_Score ~ Time + (1|ID), data = IUS_BP_long_EC, REML = TRUE)
null_lmer_IUSbp_EC <- update(M3_lmer_IUSbp_EC, formula = ~ . -Time)
BF_BIC_IUSbp_EC <- exp((BIC(null_lmer_IUSbp_EC) - BIC(M3_lmer_IUSbp_EC))/2)
BF_BIC_IUSbp_EC

## [1] 0.2040359

1 Week

IUS_B1W <- BF_data %>% 
  dplyr::select("ID", "Group", "A_PRE_IUS_total", "C_W1_IUS_total")
IUS_B1W_long <- IUS_B1W %>%
  pivot_longer(cols = c(A_PRE_IUS_total, C_W1_IUS_total),
               names_to = "Time",
               values_to = "IUS_Score")

full_lmer_IUSb1w <- lmer(IUS_Score ~ Group * Time + (1|ID), data = IUS_B1W_long, REML = TRUE)
null_lmer_IUSb1w <- update(full_lmer_IUSb1w, formula = ~ . -Time:Group)
BF_BIC_IUSb1w <- exp((BIC(null_lmer_IUSb1w) - BIC(full_lmer_IUSb1w))/2)
BF_BIC_IUSb1w

## [1] 241.0079

M2_lmer_IUSb1w <- lmer(IUS_Score ~ Time + Group + (1|ID), data = IUS_B1W_long, REML = TRUE)
null_lmer_IUSb1w <- update(M2_lmer_IUSb1w, formula = ~ . -Group)
BF_BIC_IUSb1w <- exp((BIC(null_lmer_IUSb1w) - BIC(M2_lmer_IUSb1w))/2)
BF_BIC_IUSb1w

## [1] 0.02231583

M3_lmer_IUSb1w <- lmer(IUS_Score ~ Time + Group + (1|ID), data = IUS_B1W_long, REML = TRUE)
null_lmer_IUSb1w <- update(M3_lmer_IUSb1w, formula = ~ . -Time)
BF_BIC_IUSb1w <- exp((BIC(null_lmer_IUSb1w) - BIC(M3_lmer_IUSb1w))/2)
BF_BIC_IUSb1w

## [1] 19872.83

By Group

IUS_B1w_long_I <- IUS_B1W %>%
  filter(Group == "C_Intervention") %>% 
  pivot_longer(cols = c(A_PRE_IUS_total, C_W1_IUS_total),
               names_to = "Time",
               values_to = "IUS_Score")
IUS_B1w_long_C <- IUS_B1W %>%
  filter(Group == "B_Controls") %>% 
  pivot_longer(cols = c(A_PRE_IUS_total, C_W1_IUS_total),
               names_to = "Time",
               values_to = "IUS_Score")
IUS_B1w_long_EC <- IUS_B1W %>%
  filter(Group == "A_ECs") %>% 
  pivot_longer(cols = c(A_PRE_IUS_total, C_W1_IUS_total),
               names_to = "Time",
               values_to = "IUS_Score")

full_lmer_IUSb1w_I <- lmer(IUS_Score ~ Time + (1|ID), data = IUS_B1w_long_I, REML = TRUE)
null_lmer_IUSb1w_I <- update(full_lmer_IUSb1w_I, formula = ~ . -Time)
BF_BIC_IUSb1w_I <- exp((BIC(null_lmer_IUSb1w_I) - BIC(full_lmer_IUSb1w_I))/2)
BF_BIC_IUSb1w_I

## [1] 151951.7

M2_lmer_IUSb1w_C <- lmer(IUS_Score ~ Time + (1|ID), data = IUS_B1w_long_C, REML = TRUE)
null_lmer_IUSb1w_C <- update(M2_lmer_IUSb1w_C, formula = ~ . -Time)
BF_BIC_IUSb1w_C <- exp((BIC(null_lmer_IUSb1w_C) - BIC(M2_lmer_IUSb1w_C))/2)
BF_BIC_IUSb1w_C

## [1] 4.847665

M3_lmer_IUSb1w_EC <- lmer(IUS_Score ~ Time + (1|ID), data = IUS_B1w_long_EC, REML = TRUE)
null_lmer_IUSb1w_EC <- update(M3_lmer_IUSb1w_EC, formula = ~ . -Time)
BF_BIC_IUSb1w_EC <- exp((BIC(null_lmer_IUSb1w_EC) - BIC(M3_lmer_IUSb1w_EC))/2)
BF_BIC_IUSb1w_EC

## [1] 0.395221

1 Month

IUS_B1M <- BF_data %>% 
  dplyr::select("ID", "Group", "A_PRE_IUS_total", "D_M1_IUS_total")
IUS_B1M_long <- IUS_B1M %>%
  pivot_longer(cols = c(A_PRE_IUS_total, D_M1_IUS_total),
               names_to = "Time",
               values_to = "IUS_Score")

full_lmer_IUSb1m <- lmer(IUS_Score ~ Group * Time + (1|ID), data = IUS_B1M_long, REML = TRUE)
null_lmer_IUSb1m <- update(full_lmer_IUSb1m, formula = ~ . -Time:Group)
BF_BIC_IUSb1m <- exp((BIC(null_lmer_IUSb1m) - BIC(full_lmer_IUSb1m))/2)
BF_BIC_IUSb1m

## [1] 2671.184

M2_lmer_IUSb1m <- lmer(IUS_Score ~ Time + Group + (1|ID), data = IUS_B1M_long, REML = TRUE)
null_lmer_IUSb1m <- update(M2_lmer_IUSb1m, formula = ~ . -Group)
BF_BIC_IUSb1m <- exp((BIC(null_lmer_IUSb1m) - BIC(M2_lmer_IUSb1m))/2)
BF_BIC_IUSb1m

## [1] 0.02837185

M3_lmer_IUSb1m <- lmer(IUS_Score ~ Time + Group + (1|ID), data = IUS_B1M_long, REML = TRUE)
null_lmer_IUSb1m <- update(M3_lmer_IUSb1m, formula = ~ . -Time)
BF_BIC_IUSb1m <- exp((BIC(null_lmer_IUSb1m) - BIC(M3_lmer_IUSb1m))/2)
BF_BIC_IUSb1m

## [1] 197.3966

By Group

IUS_B1m_long_I <- IUS_B1M %>%
  filter(Group == "C_Intervention") %>% 
  pivot_longer(cols = c(A_PRE_IUS_total, D_M1_IUS_total),
               names_to = "Time",
               values_to = "IUS_Score")
IUS_B1m_long_C <- IUS_B1M %>%
  filter(Group == "B_Controls") %>% 
  pivot_longer(cols = c(A_PRE_IUS_total, D_M1_IUS_total),
               names_to = "Time",
               values_to = "IUS_Score")
IUS_B1m_long_EC <- IUS_B1M %>%
  filter(Group == "A_ECs") %>% 
  pivot_longer(cols = c(A_PRE_IUS_total, D_M1_IUS_total),
               names_to = "Time",
               values_to = "IUS_Score")

full_lmer_IUSb1m_I <- lmer(IUS_Score ~ Time + (1|ID), data = IUS_B1m_long_I, REML = TRUE)
null_lmer_IUSb1m_I <- update(full_lmer_IUSb1m_I, formula = ~ . -Time)
BF_BIC_IUSb1m_I <- exp((BIC(null_lmer_IUSb1m_I) - BIC(full_lmer_IUSb1m_I))/2)
BF_BIC_IUSb1m_I

## [1] 94753.09

M2_lmer_IUSb1m_C <- lmer(IUS_Score ~ Time + (1|ID), data = IUS_B1m_long_C, REML = TRUE)
null_lmer_IUSb1m_C <- update(M2_lmer_IUSb1m_C, formula = ~ . -Time)
BF_BIC_IUSb1m_C <- exp((BIC(null_lmer_IUSb1m_C) - BIC(M2_lmer_IUSb1m_C))/2)
BF_BIC_IUSb1m_C

## [1] 1.316734

M3_lmer_IUSb1m_EC <- lmer(IUS_Score ~ Time + (1|ID), data = IUS_B1m_long_EC, REML = TRUE)
null_lmer_IUSb1m_EC <- update(M3_lmer_IUSb1m_EC, formula = ~ . -Time)
BF_BIC_IUSb1m_EC <- exp((BIC(null_lmer_IUSb1m_EC) - BIC(M3_lmer_IUSb1m_EC))/2)
BF_BIC_IUSb1m_EC

## [1] 8.081472

GM

Post

GM_BP <- BF_data %>% 
  dplyr::select("ID", "Group", "A_PRE_GM", "B_POST_GM")
GM_BP_long <- GM_BP %>%
  pivot_longer(cols = c(A_PRE_GM, B_POST_GM),
               names_to = "Time",
               values_to = "GM_Score")

full_lmer_GMbp <- lmer(GM_Score ~ Group * Time + (1|ID), data = GM_BP_long, REML = TRUE)
null_lmer_GMbp <- update(full_lmer_GMbp, formula = ~ . -Time:Group) # null means no interaction effect
BF_BIC_GMbp <- exp((BIC(null_lmer_GMbp) - BIC(full_lmer_GMbp))/2)
BF_BIC_GMbp # for the interaction

## [1] 0.3769363

M2_lmer_GMbp <- lmer(GM_Score ~ Time + Group + (1|ID), data = GM_BP_long, REML = TRUE)
null_lmer_GMbp <- update(M2_lmer_GMbp, formula = ~ . -Group) # null means no group effect
BF_BIC_GMbp <- exp((BIC(null_lmer_GMbp) - BIC(M2_lmer_GMbp))/2)
BF_BIC_GMbp # for group

## [1] 0.002426745

M3_lmer_GMbp <- lmer(GM_Score ~ Time + Group + (1|ID), data = GM_BP_long, REML = TRUE)
null_lmer_GMbp <- update(M3_lmer_GMbp, formula = ~ . -Time) # null means no time effect
BF_BIC_GMbp <- exp((BIC(null_lmer_GMbp) - BIC(M3_lmer_GMbp))/2)
BF_BIC_GMbp # for the time

## [1] 3566215

By Group

GM_BP_long_I <- GM_BP %>%
  filter(Group == "C_Intervention") %>% 
  pivot_longer(cols = c(A_PRE_GM, B_POST_GM),
               names_to = "Time",
               values_to = "GM_Score")
GM_BP_long_C <- GM_BP %>%
  filter(Group == "B_Controls") %>% 
  pivot_longer(cols = c(A_PRE_GM, B_POST_GM),
               names_to = "Time",
               values_to = "GM_Score")
GM_BP_long_EC <- GM_BP %>%
  filter(Group == "A_ECs") %>% 
  pivot_longer(cols = c(A_PRE_GM, B_POST_GM),
               names_to = "Time",
               values_to = "GM_Score")

full_lmer_GMbp_I <- lmer(GM_Score ~ Time + (1|ID), data = GM_BP_long_I, REML = TRUE)
null_lmer_GMbp_I <- update(full_lmer_GMbp_I, formula = ~ . -Time)
BF_BIC_GMbp_I <- exp((BIC(null_lmer_GMbp_I) - BIC(full_lmer_GMbp_I))/2)
BF_BIC_GMbp_I

## [1] 15556.71

M2_lmer_GMbp_C <- lmer(GM_Score ~ Time + (1|ID), data = GM_BP_long_C, REML = TRUE)
null_lmer_GMbp_C <- update(M2_lmer_GMbp_C, formula = ~ . -Time)
BF_BIC_GMbp_C <- exp((BIC(null_lmer_GMbp_C) - BIC(M2_lmer_GMbp_C))/2)
BF_BIC_GMbp_C

## [1] 470.2298

M3_lmer_GMbp_EC <- lmer(GM_Score ~ Time + (1|ID), data = GM_BP_long_EC, REML = TRUE)
null_lmer_GMbp_EC <- update(M3_lmer_GMbp_EC, formula = ~ . -Time)
BF_BIC_GMbp_EC <- exp((BIC(null_lmer_GMbp_EC) - BIC(M3_lmer_GMbp_EC))/2)
BF_BIC_GMbp_EC

## [1] 0.03192482

1 Week

GM_B1W <- BF_data %>% 
  dplyr::select("ID", "Group", "A_PRE_GM", "C_W1_GM")
GM_B1W_long <- GM_B1W %>%
  pivot_longer(cols = c(A_PRE_GM, C_W1_GM),
               names_to = "Time",
               values_to = "GM_Score")

full_lmer_GMb1w <- lmer(GM_Score ~ Group * Time + (1|ID), data = GM_B1W_long, REML = TRUE)
null_lmer_GMb1w <- update(full_lmer_GMb1w, formula = ~ . -Time:Group)
BF_BIC_GMb1w <- exp((BIC(null_lmer_GMb1w) - BIC(full_lmer_GMb1w))/2)
BF_BIC_GMb1w

## [1] 0.005563526

M2_lmer_GMb1w <- lmer(GM_Score ~ Time + Group + (1|ID), data = GM_B1W_long, REML = TRUE)
null_lmer_GMb1w <- update(M2_lmer_GMb1w, formula = ~ . -Group)
BF_BIC_GMb1w <- exp((BIC(null_lmer_GMb1w) - BIC(M2_lmer_GMb1w))/2)
BF_BIC_GMb1w

## [1] 0.004970869

M3_lmer_GMb1w <- lmer(GM_Score ~ Time + Group + (1|ID), data = GM_B1W_long, REML = TRUE)
null_lmer_GMb1w <- update(M3_lmer_GMb1w, formula = ~ . -Time)
BF_BIC_GMb1w <- exp((BIC(null_lmer_GMb1w) - BIC(M3_lmer_GMb1w))/2)
BF_BIC_GMb1w

## [1] 58.41477

By Group

GM_B1w_long_I <- GM_B1W %>%
  filter(Group == "C_Intervention") %>% 
  pivot_longer(cols = c(A_PRE_GM, C_W1_GM),
               names_to = "Time",
               values_to = "GM_Score")
GM_B1w_long_C <- GM_B1W %>%
  filter(Group == "B_Controls") %>% 
  pivot_longer(cols = c(A_PRE_GM, C_W1_GM),
               names_to = "Time",
               values_to = "GM_Score")
GM_B1w_long_EC <- GM_B1W %>%
  filter(Group == "A_ECs") %>% 
  pivot_longer(cols = c(A_PRE_GM, C_W1_GM),
               names_to = "Time",
               values_to = "GM_Score")

full_lmer_GMb1w_I <- lmer(GM_Score ~ Time + (1|ID), data = GM_B1w_long_I, REML = TRUE)
null_lmer_GMb1w_I <- update(full_lmer_GMb1w_I, formula = ~ . -Time)
BF_BIC_GMb1w_I <- exp((BIC(null_lmer_GMb1w_I) - BIC(full_lmer_GMb1w_I))/2)
BF_BIC_GMb1w_I

## [1] 199.4642

M2_lmer_GMb1w_C <- lmer(GM_Score ~ Time + (1|ID), data = GM_B1w_long_C, REML = TRUE)
null_lmer_GMb1w_C <- update(M2_lmer_GMb1w_C, formula = ~ . -Time)
BF_BIC_GMb1w_C <- exp((BIC(null_lmer_GMb1w_C) - BIC(M2_lmer_GMb1w_C))/2)
BF_BIC_GMb1w_C

## [1] 0.3694616

M3_lmer_GMb1w_EC <- lmer(GM_Score ~ Time + (1|ID), data = GM_B1w_long_EC, REML = TRUE)
null_lmer_GMb1w_EC <- update(M3_lmer_GMb1w_EC, formula = ~ . -Time)
BF_BIC_GMb1w_EC <- exp((BIC(null_lmer_GMb1w_EC) - BIC(M3_lmer_GMb1w_EC))/2)
BF_BIC_GMb1w_EC

## [1] 0.06464413

1 Month

GM_B1M <- BF_data %>% 
  dplyr::select("ID", "Group", "A_PRE_GM", "D_M1_GM")
GM_B1M_long <- GM_B1M %>%
  pivot_longer(cols = c(A_PRE_GM, D_M1_GM),
               names_to = "Time",
               values_to = "GM_Score")

full_lmer_GMb1m <- lmer(GM_Score ~ Group * Time + (1|ID), data = GM_B1M_long, REML = TRUE)
null_lmer_GMb1m <- update(full_lmer_GMb1m, formula = ~ . -Time:Group)
BF_BIC_GMb1m <- exp((BIC(null_lmer_GMb1m) - BIC(full_lmer_GMb1m))/2)
BF_BIC_GMb1m

## [1] 0.01232158

M2_lmer_GMb1m <- lmer(GM_Score ~ Time + Group + (1|ID), data = GM_B1M_long, REML = TRUE)
null_lmer_GMb1m <- update(M2_lmer_GMb1m, formula = ~ . -Group)
BF_BIC_GMb1m <- exp((BIC(null_lmer_GMb1m) - BIC(M2_lmer_GMb1m))/2)
BF_BIC_GMb1m

## [1] 0.004456573

M3_lmer_GMb1m <- lmer(GM_Score ~ Time + Group + (1|ID), data = GM_B1M_long, REML = TRUE)
null_lmer_GMb1m <- update(M3_lmer_GMb1m, formula = ~ . -Time)
BF_BIC_GMb1m <- exp((BIC(null_lmer_GMb1m) - BIC(M3_lmer_GMb1m))/2)
BF_BIC_GMb1m

## [1] 216.1675

By Group

GM_B1m_long_I <- GM_B1M %>%
  filter(Group == "C_Intervention") %>% 
  pivot_longer(cols = c(A_PRE_GM, D_M1_GM),
               names_to = "Time",
               values_to = "GM_Score")
GM_B1m_long_C <- GM_B1M %>%
  filter(Group == "B_Controls") %>% 
  pivot_longer(cols = c(A_PRE_GM, D_M1_GM),
               names_to = "Time",
               values_to = "GM_Score")
GM_B1m_long_EC <- GM_B1M %>%
  filter(Group == "A_ECs") %>% 
  pivot_longer(cols = c(A_PRE_GM, D_M1_GM),
               names_to = "Time",
               values_to = "GM_Score")

full_lmer_GMb1m_I <- lmer(GM_Score ~ Time + (1|ID), data = GM_B1m_long_I, REML = TRUE)
null_lmer_GMb1m_I <- update(full_lmer_GMb1m_I, formula = ~ . -Time)
BF_BIC_GMb1m_I <- exp((BIC(null_lmer_GMb1m_I) - BIC(full_lmer_GMb1m_I))/2)
BF_BIC_GMb1m_I

## [1] 331.0491

M2_lmer_GMb1m_C <- lmer(GM_Score ~ Time + (1|ID), data = GM_B1m_long_C, REML = TRUE)
null_lmer_GMb1m_C <- update(M2_lmer_GMb1m_C, formula = ~ . -Time)
BF_BIC_GMb1m_C <- exp((BIC(null_lmer_GMb1m_C) - BIC(M2_lmer_GMb1m_C))/2)
BF_BIC_GMb1m_C

## [1] 0.6108577

M3_lmer_GMb1m_EC <- lmer(GM_Score ~ Time + (1|ID), data = GM_B1m_long_EC, REML = TRUE)
null_lmer_GMb1m_EC <- update(M3_lmer_GMb1m_EC, formula = ~ . -Time)
BF_BIC_GMb1m_EC <- exp((BIC(null_lmer_GMb1m_EC) - BIC(M3_lmer_GMb1m_EC))/2)
BF_BIC_GMb1m_EC

## [1] 0.05271971

PHQ

1 Week

PHQ_B1W <- BF_data %>% 
  dplyr::select("ID", "Group", "A_PRE_PHQ_total", "C_W1_PHQ_total")
PHQ_B1W_long <- PHQ_B1W %>%
  pivot_longer(cols = c(A_PRE_PHQ_total, C_W1_PHQ_total),
               names_to = "Time",
               values_to = "PHQ_Score")

full_lmer_PHQb1w <- lmer(PHQ_Score ~ Group * Time + (1|ID), data = PHQ_B1W_long, REML = TRUE)
null_lmer_PHQb1w <- update(full_lmer_PHQb1w, formula = ~ . -Time:Group)
BF_BIC_PHQb1w <- exp((BIC(null_lmer_PHQb1w) - BIC(full_lmer_PHQb1w))/2)
BF_BIC_PHQb1w

## [1] 0.02531879

M2_lmer_PHQb1w <- lmer(PHQ_Score ~ Time + Group + (1|ID), data = PHQ_B1W_long, REML = TRUE)
null_lmer_PHQb1w <- update(M2_lmer_PHQb1w, formula = ~ . -Group)
BF_BIC_PHQb1w <- exp((BIC(null_lmer_PHQb1w) - BIC(M2_lmer_PHQb1w))/2)
BF_BIC_PHQb1w

## [1] 0.02038741

M3_lmer_PHQb1w <- lmer(PHQ_Score ~ Time + Group + (1|ID), data = PHQ_B1W_long, REML = TRUE)
null_lmer_PHQb1w <- update(M3_lmer_PHQb1w, formula = ~ . -Time)
BF_BIC_PHQb1w <- exp((BIC(null_lmer_PHQb1w) - BIC(M3_lmer_PHQb1w))/2)
BF_BIC_PHQb1w

## [1] 9.848049

By Group

PHQ_B1w_long_I <- PHQ_B1W %>%
  filter(Group == "C_Intervention") %>% 
  pivot_longer(cols = c(A_PRE_PHQ_total, C_W1_PHQ_total),
               names_to = "Time",
               values_to = "PHQ_Score")
PHQ_B1w_long_C <- PHQ_B1W %>%
  filter(Group == "B_Controls") %>% 
  pivot_longer(cols = c(A_PRE_PHQ_total, C_W1_PHQ_total),
               names_to = "Time",
               values_to = "PHQ_Score")
PHQ_B1w_long_EC <- PHQ_B1W %>%
  filter(Group == "A_ECs") %>% 
  pivot_longer(cols = c(A_PRE_PHQ_total, C_W1_PHQ_total),
               names_to = "Time",
               values_to = "PHQ_Score")

full_lmer_PHQb1w_I <- lmer(PHQ_Score ~ Time + (1|ID), data = PHQ_B1w_long_I, REML = TRUE)
null_lmer_PHQb1w_I <- update(full_lmer_PHQb1w_I, formula = ~ . -Time)
BF_BIC_PHQb1w_I <- exp((BIC(null_lmer_PHQb1w_I) - BIC(full_lmer_PHQb1w_I))/2)
BF_BIC_PHQb1w_I

## [1] 6.213481

M2_lmer_PHQb1w_C <- lmer(PHQ_Score ~ Time + (1|ID), data = PHQ_B1w_long_C, REML = TRUE)
null_lmer_PHQb1w_C <- update(M2_lmer_PHQb1w_C, formula = ~ . -Time)
BF_BIC_PHQb1w_C <- exp((BIC(null_lmer_PHQb1w_C) - BIC(M2_lmer_PHQb1w_C))/2)
BF_BIC_PHQb1w_C

## [1] 0.6586489

M3_lmer_PHQb1w_EC <- lmer(PHQ_Score ~ Time + (1|ID), data = PHQ_B1w_long_EC, REML = TRUE)
null_lmer_PHQb1w_EC <- update(M3_lmer_PHQb1w_EC, formula = ~ . -Time)
BF_BIC_PHQb1w_EC <- exp((BIC(null_lmer_PHQb1w_EC) - BIC(M3_lmer_PHQb1w_EC))/2)
BF_BIC_PHQb1w_EC

## [1] 0.1382616

1 Month

PHQ_B1M <- BF_data %>% 
  dplyr::select("ID", "Group", "A_PRE_PHQ_total", "D_M1_PHQ_total")
PHQ_B1M_long <- PHQ_B1M %>%
  pivot_longer(cols = c(A_PRE_PHQ_total, D_M1_PHQ_total),
               names_to = "Time",
               values_to = "PHQ_Score")

full_lmer_PHQb1m <- lmer(PHQ_Score ~ Group * Time + (1|ID), data = PHQ_B1M_long, REML = TRUE)
null_lmer_PHQb1m <- update(full_lmer_PHQb1m, formula = ~ . -Time:Group)
BF_BIC_PHQb1m <- exp((BIC(null_lmer_PHQb1m) - BIC(full_lmer_PHQb1m))/2)
BF_BIC_PHQb1m

## [1] 0.7440856

M2_lmer_PHQb1m <- lmer(PHQ_Score ~ Time + Group + (1|ID), data = PHQ_B1M_long, REML = TRUE)
null_lmer_PHQb1m <- update(M2_lmer_PHQb1m, formula = ~ . -Group)
BF_BIC_PHQb1m <- exp((BIC(null_lmer_PHQb1m) - BIC(M2_lmer_PHQb1m))/2)
BF_BIC_PHQb1m

## [1] 0.02448161

M3_lmer_PHQb1m <- lmer(PHQ_Score ~ Time + Group + (1|ID), data = PHQ_B1M_long, REML = TRUE)
null_lmer_PHQb1m <- update(M3_lmer_PHQb1m, formula = ~ . -Time)
BF_BIC_PHQb1m <- exp((BIC(null_lmer_PHQb1m) - BIC(M3_lmer_PHQb1m))/2)
BF_BIC_PHQb1m

## [1] 22.08378

By Group

PHQ_B1m_long_I <- PHQ_B1M %>%
  filter(Group == "C_Intervention") %>% 
  pivot_longer(cols = c(A_PRE_PHQ_total, D_M1_PHQ_total),
               names_to = "Time",
               values_to = "PHQ_Score")
PHQ_B1m_long_C <- PHQ_B1M %>%
  filter(Group == "B_Controls") %>% 
  pivot_longer(cols = c(A_PRE_PHQ_total, D_M1_PHQ_total),
               names_to = "Time",
               values_to = "PHQ_Score")
PHQ_B1m_long_EC <- PHQ_B1M %>%
  filter(Group == "A_ECs") %>% 
  pivot_longer(cols = c(A_PRE_PHQ_total, D_M1_PHQ_total),
               names_to = "Time",
               values_to = "PHQ_Score")

full_lmer_PHQb1m_I <- lmer(PHQ_Score ~ Time + (1|ID), data = PHQ_B1m_long_I, REML = TRUE)
null_lmer_PHQb1m_I <- update(full_lmer_PHQb1m_I, formula = ~ . -Time)
BF_BIC_PHQb1m_I <- exp((BIC(null_lmer_PHQb1m_I) - BIC(full_lmer_PHQb1m_I))/2)
BF_BIC_PHQb1m_I

## [1] 527.3174

M2_lmer_PHQb1m_C <- lmer(PHQ_Score ~ Time + (1|ID), data = PHQ_B1m_long_C, REML = TRUE)
null_lmer_PHQb1m_C <- update(M2_lmer_PHQb1m_C, formula = ~ . -Time)
BF_BIC_PHQb1m_C <- exp((BIC(null_lmer_PHQb1m_C) - BIC(M2_lmer_PHQb1m_C))/2)
BF_BIC_PHQb1m_C

## [1] 0.6946523

M3_lmer_PHQb1m_EC <- lmer(PHQ_Score ~ Time + (1|ID), data = PHQ_B1m_long_EC, REML = TRUE)
null_lmer_PHQb1m_EC <- update(M3_lmer_PHQb1m_EC, formula = ~ . -Time)
BF_BIC_PHQb1m_EC <- exp((BIC(null_lmer_PHQb1m_EC) - BIC(M3_lmer_PHQb1m_EC))/2)
BF_BIC_PHQb1m_EC

## [1] 0.2919433

GAD

1 Week

GAD_B1W <- BF_data %>% 
  dplyr::select("ID", "Group", "A_PRE_GAD_total", "C_W1_GAD_total")
GAD_B1W_long <- GAD_B1W %>%
  pivot_longer(cols = c(A_PRE_GAD_total, C_W1_GAD_total),
               names_to = "Time",
               values_to = "GAD_Score")

full_lmer_GADb1w <- lmer(GAD_Score ~ Group * Time + (1|ID), data = GAD_B1W_long, REML = TRUE)
null_lmer_GADb1w <- update(full_lmer_GADb1w, formula = ~ . -Time:Group)
BF_BIC_GADb1w <- exp((BIC(null_lmer_GADb1w) - BIC(full_lmer_GADb1w))/2)
BF_BIC_GADb1w

## [1] 0.02352804

M2_lmer_GADb1w <- lmer(GAD_Score ~ Time + Group + (1|ID), data = GAD_B1W_long, REML = TRUE)
null_lmer_GADb1w <- update(M2_lmer_GADb1w, formula = ~ . -Group)
BF_BIC_GADb1w <- exp((BIC(null_lmer_GADb1w) - BIC(M2_lmer_GADb1w))/2)
BF_BIC_GADb1w

## [1] 0.01108104

M3_lmer_GADb1w <- lmer(GAD_Score ~ Time + Group + (1|ID), data = GAD_B1W_long, REML = TRUE)
null_lmer_GADb1w <- update(M3_lmer_GADb1w, formula = ~ . -Time)
BF_BIC_GADb1w <- exp((BIC(null_lmer_GADb1w) - BIC(M3_lmer_GADb1w))/2)
BF_BIC_GADb1w

## [1] 0.3252588

By Group

GAD_B1w_long_I <- GAD_B1W %>%
  filter(Group == "C_Intervention") %>% 
  pivot_longer(cols = c(A_PRE_GAD_total, C_W1_GAD_total),
               names_to = "Time",
               values_to = "GAD_Score")
GAD_B1w_long_C <- GAD_B1W %>%
  filter(Group == "B_Controls") %>% 
  pivot_longer(cols = c(A_PRE_GAD_total, C_W1_GAD_total),
               names_to = "Time",
               values_to = "GAD_Score")
GAD_B1w_long_EC <- GAD_B1W %>%
  filter(Group == "A_ECs") %>% 
  pivot_longer(cols = c(A_PRE_GAD_total, C_W1_GAD_total),
               names_to = "Time",
               values_to = "GAD_Score")

full_lmer_GADb1w_I <- lmer(GAD_Score ~ Time + (1|ID), data = GAD_B1w_long_I, REML = TRUE)
null_lmer_GADb1w_I <- update(full_lmer_GADb1w_I, formula = ~ . -Time)
BF_BIC_GADb1w_I <- exp((BIC(null_lmer_GADb1w_I) - BIC(full_lmer_GADb1w_I))/2)
BF_BIC_GADb1w_I

## [1] 0.76916

M2_lmer_GADb1w_C <- lmer(GAD_Score ~ Time + (1|ID), data = GAD_B1w_long_C, REML = TRUE)
null_lmer_GADb1w_C <- update(M2_lmer_GADb1w_C, formula = ~ . -Time)
BF_BIC_GADb1w_C <- exp((BIC(null_lmer_GADb1w_C) - BIC(M2_lmer_GADb1w_C))/2)
BF_BIC_GADb1w_C

## [1] 0.2132653

M3_lmer_GADb1w_EC <- lmer(GAD_Score ~ Time + (1|ID), data = GAD_B1w_long_EC, REML = TRUE)
null_lmer_GADb1w_EC <- update(M3_lmer_GADb1w_EC, formula = ~ . -Time)
BF_BIC_GADb1w_EC <- exp((BIC(null_lmer_GADb1w_EC) - BIC(M3_lmer_GADb1w_EC))/2)
BF_BIC_GADb1w_EC

## [1] 0.1446356

1 Month

GAD_B1M <- BF_data %>% 
  dplyr::select("ID", "Group", "A_PRE_GAD_total", "D_M1_GAD_total")
GAD_B1M_long <- GAD_B1M %>%
  pivot_longer(cols = c(A_PRE_GAD_total, D_M1_GAD_total),
               names_to = "Time",
               values_to = "GAD_Score")

full_lmer_GADb1m <- lmer(GAD_Score ~ Group * Time + (1|ID), data = GAD_B1M_long, REML = TRUE)
null_lmer_GADb1m <- update(full_lmer_GADb1m, formula = ~ . -Time:Group)
BF_BIC_GADb1m <- exp((BIC(null_lmer_GADb1m) - BIC(full_lmer_GADb1m))/2)
BF_BIC_GADb1m

## [1] 1.697957

M2_lmer_GADb1m <- lmer(GAD_Score ~ Time + Group + (1|ID), data = GAD_B1M_long, REML = TRUE)
null_lmer_GADb1m <- update(M2_lmer_GADb1m, formula = ~ . -Group)
BF_BIC_GADb1m <- exp((BIC(null_lmer_GADb1m) - BIC(M2_lmer_GADb1m))/2)
BF_BIC_GADb1m

## [1] 0.009551524

M3_lmer_GADb1m <- lmer(GAD_Score ~ Time + Group + (1|ID), data = GAD_B1M_long, REML = TRUE)
null_lmer_GADb1m <- update(M3_lmer_GADb1m, formula = ~ . -Time)
BF_BIC_GADb1m <- exp((BIC(null_lmer_GADb1m) - BIC(M3_lmer_GADb1m))/2)
BF_BIC_GADb1m

## [1] 0.9999603

By Group

GAD_B1m_long_I <- GAD_B1M %>%
  filter(Group == "C_Intervention") %>% 
  pivot_longer(cols = c(A_PRE_GAD_total, D_M1_GAD_total),
               names_to = "Time",
               values_to = "GAD_Score")
GAD_B1m_long_C <- GAD_B1M %>%
  filter(Group == "B_Controls") %>% 
  pivot_longer(cols = c(A_PRE_GAD_total, D_M1_GAD_total),
               names_to = "Time",
               values_to = "GAD_Score")
GAD_B1m_long_EC <- GAD_B1M %>%
  filter(Group == "A_ECs") %>% 
  pivot_longer(cols = c(A_PRE_GAD_total, D_M1_GAD_total),
               names_to = "Time",
               values_to = "GAD_Score")

full_lmer_GADb1m_I <- lmer(GAD_Score ~ Time + (1|ID), data = GAD_B1m_long_I, REML = TRUE)
null_lmer_GADb1m_I <- update(full_lmer_GADb1m_I, formula = ~ . -Time)
BF_BIC_GADb1m_I <- exp((BIC(null_lmer_GADb1m_I) - BIC(full_lmer_GADb1m_I))/2)
BF_BIC_GADb1m_I

## [1] 52.05475

M2_lmer_GADb1m_C <- lmer(GAD_Score ~ Time + (1|ID), data = GAD_B1m_long_C, REML = TRUE)
null_lmer_GADb1m_C <- update(M2_lmer_GADb1m_C, formula = ~ . -Time)
BF_BIC_GADb1m_C <- exp((BIC(null_lmer_GADb1m_C) - BIC(M2_lmer_GADb1m_C))/2)
BF_BIC_GADb1m_C

## [1] 0.4528085

M3_lmer_GADb1m_EC <- lmer(GAD_Score ~ Time + (1|ID), data = GAD_B1m_long_EC, REML = TRUE)
null_lmer_GADb1m_EC <- update(M3_lmer_GADb1m_EC, formula = ~ . -Time)
BF_BIC_GADb1m_EC <- exp((BIC(null_lmer_GADb1m_EC) - BIC(M3_lmer_GADb1m_EC))/2)
BF_BIC_GADb1m_EC

## [1] 0.7801084

FI

1 Week

FI_B1W <- BF_data %>% 
  dplyr::select("ID", "Group", "A_PRE_FI_total", "C_W1_FI_total")
FI_B1W_long <- FI_B1W %>%
  pivot_longer(cols = c(A_PRE_FI_total, C_W1_FI_total),
               names_to = "Time",
               values_to = "FI_Score")

full_lmer_FIb1w <- lmer(FI_Score ~ Group * Time + (1|ID), data = FI_B1W_long, REML = TRUE)
null_lmer_FIb1w <- update(full_lmer_FIb1w, formula = ~ . -Time:Group)
BF_BIC_FIb1w <- exp((BIC(null_lmer_FIb1w) - BIC(full_lmer_FIb1w))/2)
BF_BIC_FIb1w

## [1] 0.02112748

M2_lmer_FIb1w <- lmer(FI_Score ~ Time + Group + (1|ID), data = FI_B1W_long, REML = TRUE)
null_lmer_FIb1w <- update(M2_lmer_FIb1w, formula = ~ . -Group)
BF_BIC_FIb1w <- exp((BIC(null_lmer_FIb1w) - BIC(M2_lmer_FIb1w))/2)
BF_BIC_FIb1w

## [1] 0.003390009

M3_lmer_FIb1w <- lmer(FI_Score ~ Time + Group + (1|ID), data = FI_B1W_long, REML = TRUE)
null_lmer_FIb1w <- update(M3_lmer_FIb1w, formula = ~ . -Time)
BF_BIC_FIb1w <- exp((BIC(null_lmer_FIb1w) - BIC(M3_lmer_FIb1w))/2)
BF_BIC_FIb1w

## [1] 0.1681013

By Group

FI_B1w_long_I <- FI_B1W %>%
  filter(Group == "C_Intervention") %>% 
  pivot_longer(cols = c(A_PRE_FI_total, C_W1_FI_total),
               names_to = "Time",
               values_to = "FI_Score")
FI_B1w_long_C <- FI_B1W %>%
  filter(Group == "B_Controls") %>% 
  pivot_longer(cols = c(A_PRE_FI_total, C_W1_FI_total),
               names_to = "Time",
               values_to = "FI_Score")
FI_B1w_long_EC <- FI_B1W %>%
  filter(Group == "A_ECs") %>% 
  pivot_longer(cols = c(A_PRE_FI_total, C_W1_FI_total),
               names_to = "Time",
               values_to = "FI_Score")

full_lmer_FIb1w_I <- lmer(FI_Score ~ Time + (1|ID), data = FI_B1w_long_I, REML = TRUE)
null_lmer_FIb1w_I <- update(full_lmer_FIb1w_I, formula = ~ . -Time)
BF_BIC_FIb1w_I <- exp((BIC(null_lmer_FIb1w_I) - BIC(full_lmer_FIb1w_I))/2)
BF_BIC_FIb1w_I

## [1] 1.596137

M2_lmer_FIb1w_C <- lmer(FI_Score ~ Time + (1|ID), data = FI_B1w_long_C, REML = TRUE)
null_lmer_FIb1w_C <- update(M2_lmer_FIb1w_C, formula = ~ . -Time)
BF_BIC_FIb1w_C <- exp((BIC(null_lmer_FIb1w_C) - BIC(M2_lmer_FIb1w_C))/2)
BF_BIC_FIb1w_C

## [1] 0.07294422

M3_lmer_FIb1w_EC <- lmer(FI_Score ~ Time + (1|ID), data = FI_B1w_long_EC, REML = TRUE)
null_lmer_FIb1w_EC <- update(M3_lmer_FIb1w_EC, formula = ~ . -Time)
BF_BIC_FIb1w_EC <- exp((BIC(null_lmer_FIb1w_EC) - BIC(M3_lmer_FIb1w_EC))/2)
BF_BIC_FIb1w_EC

## [1] 0.1252424

1 Month

FI_B1M <- BF_data %>% 
  dplyr::select("ID", "Group", "A_PRE_FI_total", "D_M1_FI_total")
FI_B1M_long <- FI_B1M %>%
  pivot_longer(cols = c(A_PRE_FI_total, D_M1_FI_total),
               names_to = "Time",
               values_to = "FI_Score")

full_lmer_FIb1m <- lmer(FI_Score ~ Group * Time + (1|ID), data = FI_B1M_long, REML = TRUE)
null_lmer_FIb1m <- update(full_lmer_FIb1m, formula = ~ . -Time:Group)
BF_BIC_FIb1m <- exp((BIC(null_lmer_FIb1m) - BIC(full_lmer_FIb1m))/2)
BF_BIC_FIb1m

## [1] 0.01319748

M2_lmer_FIb1m <- lmer(FI_Score ~ Time + Group + (1|ID), data = FI_B1M_long, REML = TRUE)
null_lmer_FIb1m <- update(M2_lmer_FIb1m, formula = ~ . -Group)
BF_BIC_FIb1m <- exp((BIC(null_lmer_FIb1m) - BIC(M2_lmer_FIb1m))/2)
BF_BIC_FIb1m

## [1] 0.004394219

M3_lmer_FIb1m <- lmer(FI_Score ~ Time + Group + (1|ID), data = FI_B1M_long, REML = TRUE)
null_lmer_FIb1m <- update(M3_lmer_FIb1m, formula = ~ . -Time)
BF_BIC_FIb1m <- exp((BIC(null_lmer_FIb1m) - BIC(M3_lmer_FIb1m))/2)
BF_BIC_FIb1m

## [1] 7.371853

By Group

FI_B1m_long_I <- FI_B1M %>%
  filter(Group == "C_Intervention") %>% 
  pivot_longer(cols = c(A_PRE_FI_total, D_M1_FI_total),
               names_to = "Time",
               values_to = "FI_Score")
FI_B1m_long_C <- FI_B1M %>%
  filter(Group == "B_Controls") %>% 
  pivot_longer(cols = c(A_PRE_FI_total, D_M1_FI_total),
               names_to = "Time",
               values_to = "FI_Score")
FI_B1m_long_EC <- FI_B1M %>%
  filter(Group == "A_ECs") %>% 
  pivot_longer(cols = c(A_PRE_FI_total, D_M1_FI_total),
               names_to = "Time",
               values_to = "FI_Score")

full_lmer_FIb1m_I <- lmer(FI_Score ~ Time + (1|ID), data = FI_B1m_long_I, REML = TRUE)
null_lmer_FIb1m_I <- update(full_lmer_FIb1m_I, formula = ~ . -Time)
BF_BIC_FIb1m_I <- exp((BIC(null_lmer_FIb1m_I) - BIC(full_lmer_FIb1m_I))/2)
BF_BIC_FIb1m_I

## [1] 2.89944

M2_lmer_FIb1m_C <- lmer(FI_Score ~ Time + (1|ID), data = FI_B1m_long_C, REML = TRUE)
null_lmer_FIb1m_C <- update(M2_lmer_FIb1m_C, formula = ~ . -Time)
BF_BIC_FIb1m_C <- exp((BIC(null_lmer_FIb1m_C) - BIC(M2_lmer_FIb1m_C))/2)
BF_BIC_FIb1m_C

## [1] 0.8673062

M3_lmer_FIb1m_EC <- lmer(FI_Score ~ Time + (1|ID), data = FI_B1m_long_EC, REML = TRUE)
null_lmer_FIb1m_EC <- update(M3_lmer_FIb1m_EC, formula = ~ . -Time)
BF_BIC_FIb1m_EC <- exp((BIC(null_lmer_FIb1m_EC) - BIC(M3_lmer_FIb1m_EC))/2)
BF_BIC_FIb1m_EC

## [1] 0.1167154