The Beck Depression Inventory

Scale Description

Depressive symptomatolgoy was evaluated using the validated 21-item Beck Depression Inventory (BDI), which is designed to measure characteristic attitudes and symptoms of depression. The score for each item is on a 4-point scale that ranges between 0-3, where 0 denotes less depression (e.g. “I do not feel sad”) and 3 denotes more depression (e.g. “I am so sad and unhappy that I can’t stand it”). The Sleeping Pattern (Item 16) and Changes in Appetite (Item 18) items contain seven options rated, in order, 0, 1a, 1b, 2a, 2b, 3a, 3b, to differentiate between increases and decreases in behavior or motivation.Total score will be between 0-63, where 0 indicates no depression and 63 indicates extreme depression. The score is a total of 21 questions and the following is a typical subscale range distribution:

1-10: These ups and downs are considered normal 11-16: Mild mood disturbance 17-20: Borderline clinical depression 21-30: Moderate depression 31-40: Severe depression over 40: Extreme depression

Data analysis: Beck Depression Inventory

We will perform a linear regression to regress follow-up BDI composite scores on condition, a multiple regression analysis to regress follow-up BDI composite scores on condition controlling for baseline BDI scores, and a mixed effect linear model analysis to assess the effect of study condition on change in BDI scores from baseline to follow-up. In order to analyze the Beck Depression Inventory scores across condition, participants with unknown conditions or who did not complete the inventory (N = 14) were excluded from the final analysis. A composite score was created for BDI scores reported at both baseline and follow-up.

#Analysis of the PRE/POST BDI Scale for the Equanimity Study (categorical: condition 1 and condition 2)

df_eq_bdi <- read.csv(file = "/Users/mlipsett/Desktop/Equanimity/Database/Equanimity_PrePost Qualtrics_FINAL.csv", stringsAsFactors = FALSE, header = TRUE) #Reading in the csv file with BDI (and other) data for the equanimity study

df_eq_bdi <- df_eq_bdi [-c(13, 38),] #Removing PIDS 210 and 502 whose intervention condition is unknown / participant not found in masterlist

df_eq_bdi_missingpost <- df_eq_bdi [-c(1, 19, 28, 44, 55, 56, 58, 64, 74, 76, 82, 88),]# Ddata frame removing participants whose condition is unknown or who did not complete post-qualtrics measures. Note: this is only for reference - this data is not being used. Rather missing data at post is being imputed with the seties mean (see below). 

df_eq_bdi_control <- read.csv(file = "/Users/mlipsett/Desktop/Equanimity/Database/Missing Data/Equanimity_Participant Data Issues.csv", stringsAsFactors = FALSE, header = TRUE) #Reading in the csv file with missing data information 

library(dplyr)

#Subsetting the data to include only Pre/Post BDI data 
df_eq_bdi <- dplyr::select(df_eq_bdi, PID,  Condition_Group,    Post_Qual_Completed, BDI_II_sadness,    BDI_II_pessimism,   BDI_II_failure, BDI_II_pleasure,    BDI_II_guilt,   BDI_II_punishment,  BDI_II_self.dislike,    BDI_II_self.criticalness,   BDI_II_suicide, BDI_II_crying,  BDI_II_agitation,   BDI_II_interest,    BDI_II_indecisvness,    BDI_II_worthlessness,   BDI_II_energy,  BDI_II_sleep,   BDI_II_irritability,    BDI_II_appetite,    BDI_II_concentration,   BDI_II_fatigue, BDI_II_sex,Post_BDI_II_sadness, Post_BDI_II_pessimism,  Post_BDI_II_failure,    Post_BDI_II_pleasure,   Post_BDI_II_guilt,  Post_BDI_II_punishment, Post_BDI_II_self.dislike,   Post_BDI_II_self.criticalness,  Post_BDI_II_suicide,    Post_BDI_II_crying, Post_BDI_II_agitation,  Post_BDI_II_interest,   Post_BDI_II_indecisvness,   Post_BDI_II_worthlessness,  Post_BDI_II_energy, Post_BDI_II_sleep,  Post_BDI_II_irritability,   Post_BDI_II_appetite,   Post_BDI_II_concentration,  Post_BDI_II_fatigue,    Post_BDI_II_sex,    Post_CESD_felt_bothered)

results='hide'

In order to analyze the Beck Depression Inventory scores across condition, participants with unknown experimental conditions / who were not able to be identified (N = 2) were excluded from the final analysis. The total sample size for final analysis: (N = 113). A composite score was created for BDI scores reported at both baseline and follow-up.

Creating Composite Scores for BDI

#Composite Score - Pre BDI
df_eq_bdi$bdi_pre <- (df_eq_bdi$BDI_II_sadness + df_eq_bdi$BDI_II_pessimism + df_eq_bdi$BDI_II_failure + df_eq_bdi$BDI_II_pleasure + df_eq_bdi$BDI_II_guilt + df_eq_bdi$BDI_II_punishment + df_eq_bdi$BDI_II_self.dislike + df_eq_bdi$BDI_II_self.criticalness + df_eq_bdi$BDI_II_suicide + df_eq_bdi$BDI_II_crying + df_eq_bdi$BDI_II_agitation + df_eq_bdi$BDI_II_interest + df_eq_bdi$BDI_II_indecisvness + df_eq_bdi$BDI_II_worthlessness + df_eq_bdi$BDI_II_energy + df_eq_bdi$BDI_II_sleep + df_eq_bdi$BDI_II_irritability + df_eq_bdi$BDI_II_appetite + df_eq_bdi$BDI_II_concentration + df_eq_bdi$BDI_II_fatigue + df_eq_bdi$BDI_II_sex)

#Composite Score - Post BDI
df_eq_bdi$bdi_post <- (df_eq_bdi$Post_BDI_II_sadness + df_eq_bdi$Post_BDI_II_pessimism + df_eq_bdi$Post_BDI_II_failure + df_eq_bdi$Post_BDI_II_pleasure + df_eq_bdi$Post_BDI_II_guilt + df_eq_bdi$Post_BDI_II_punishment + df_eq_bdi$Post_BDI_II_self.dislike + df_eq_bdi$Post_BDI_II_self.criticalness + df_eq_bdi$Post_BDI_II_suicide + df_eq_bdi$Post_BDI_II_crying + df_eq_bdi$Post_BDI_II_agitation + df_eq_bdi$Post_BDI_II_interest + df_eq_bdi$Post_BDI_II_indecisvness + df_eq_bdi$Post_BDI_II_worthlessness + df_eq_bdi$Post_BDI_II_energy + df_eq_bdi$Post_BDI_II_sleep + df_eq_bdi$Post_BDI_II_irritability + df_eq_bdi$Post_BDI_II_appetite + df_eq_bdi$Post_BDI_II_concentration + df_eq_bdi$Post_BDI_II_fatigue + df_eq_bdi$Post_BDI_II_sex)

results='hide'

Subsetting for composite scores

df_eq_bdi <- dplyr::select(df_eq_bdi, PID,  Condition_Group, bdi_pre, bdi_post) 

#Recoding Condition_group as a factor
df_eq_bdi <- df_eq_bdi %>% 
  mutate(group_fct = factor(Condition_Group,
                            labels = c("Group 1", "Group 2")))

results='hide'

Importing Individual Differences Variables at Baseline

df_eq_qual <- read.csv(file = "/Users/mlipsett/Desktop/Equanimity/Database/Equanimity_PrePost Qualtrics_FINAL.csv", stringsAsFactors = FALSE, header = TRUE) #Reading in the csv file with qualtrics data for the equanimity study

df_eq_qual <- df_eq_qual [-c(13, 38),]#Removing participants whose condition is unknown

library(dplyr)

#Subsetting the data to include only baseline measures for individual differences 
df_eq_baseID <- dplyr::select(df_eq_qual, PID,  Condition_Group,    Post_Qual_Completed,
                      
SES_1_self_placement,   SES_2,  SES_3,  SES_4,  SES_5,  SES_6,  SES_7,  SES_8,  SES_9,  SES_10, #SES Ladder

UCLA_in_tune,   UCLA_lack_companionship,    UCLA_no_one_to_turn_to, UCLA_feel_alone,    UCLA_part_of_group, UCLA_common,    UCLA_not_close, UCLA_unshared_interest, UCLA_outgoing,  UCLA_close_to_people,   UCLA_feel_left_out, UCLA_not_meaningful_relationships,  UCLA_no_one_knows_you,  UCLA_isolated,  UCLA_can_find_companionshp, UCLA_people_understand, UCLA_feel_shy,  UCLA_no_one_with_you,   UCLA_talk_to_ppl,   UCLA_turn_to_ppl, #UCLA loneliness inventory 

DEMO_gender,    DEMO_age,   DEMO_dom_hand,  DEMO_born_U.S., DEMO_U.S._citizenship,  DEMO_no_lived_U.S., DEMO_lang_at_home,  DEMO_lang_at_home,  DEMO_ethnic_all,    DEMO_gen_status,    DEMO_relationship,  DEMO_college_member,    DEMO_college_member,    DEMO_college_member,    DEMO_college_member,    DEMO_college_member,    DEMO_college_member,    DEMO_college_member,    DEMO_employ_status, DEMO_housing_all,   DEMO_no_roommates,  DEMO_educ_mother,   DEMO_educ_father) #Demographics

Creating Composite Scores for Subjective Social Standing SSS (aka SES Ladder), UCLA loneliness inventory,

#Composite Score - Pre SSS
df_eq_baseID$SSS_pre <- (df_eq_baseID$SES_1_self_placement +    df_eq_baseID$SES_2 +    df_eq_baseID$SES_3 +    df_eq_baseID$SES_4 +    df_eq_baseID$SES_5 +    df_eq_baseID$SES_6 +    df_eq_baseID$SES_7 +    df_eq_baseID$SES_8 +    df_eq_baseID$SES_9 +    df_eq_baseID$SES_10)

#Reverse Scoring UCLA items 1, 5, 6, 9, 10, 15, 16, 19, 20
df_eq_baseID$UCLA_in_tune_r <- 5-df_eq_baseID$UCLA_in_tune
df_eq_baseID$UCLA_part_of_group_r <- 5-df_eq_baseID$UCLA_part_of_group
df_eq_baseID$UCLA_common_r <- 5-df_eq_baseID$UCLA_common
df_eq_baseID$UCLA_outgoing_r <- 5-df_eq_baseID$UCLA_outgoing
df_eq_baseID$UCLA_close_to_people_r <- 5-df_eq_baseID$UCLA_close_to_people
df_eq_baseID$UCLA_can_find_companionshp_r <- 5-df_eq_baseID$UCLA_can_find_companionshp
df_eq_baseID$UCLA_people_understand_r <- 5-df_eq_baseID$UCLA_people_understand
df_eq_baseID$UCLA_talk_to_ppl_r <- 5-df_eq_baseID$UCLA_talk_to_ppl
df_eq_baseID$UCLA_turn_to_ppl_r <- 5-df_eq_baseID$UCLA_turn_to_ppl

#Composite Score - Pre UCLA
df_eq_baseID$UCLA_pre <- (df_eq_baseID$UCLA_in_tune_r +  df_eq_baseID$UCLA_part_of_group_r + df_eq_baseID$UCLA_common_r + df_eq_baseID$UCLA_outgoing_r + df_eq_baseID$UCLA_close_to_people_r + df_eq_baseID$UCLA_can_find_companionshp_r + df_eq_baseID$UCLA_people_understand_r + df_eq_baseID$UCLA_talk_to_ppl_r + df_eq_baseID$UCLA_turn_to_ppl_r + df_eq_baseID$UCLA_lack_companionship + df_eq_baseID$UCLA_no_one_to_turn_to + df_eq_baseID$UCLA_feel_alone + df_eq_baseID$UCLA_not_close + df_eq_baseID$UCLA_unshared_interest + df_eq_baseID$UCLA_feel_left_out + df_eq_baseID$UCLA_not_meaningful_relationships + df_eq_baseID$UCLA_no_one_knows_you + df_eq_baseID$UCLA_isolated + df_eq_baseID$UCLA_feel_shy + df_eq_baseID$UCLA_no_one_with_you) 

#Demographics variables - NEED TO FIGURE OUT HOW / IF TO AGGREGATE
#DEMO_gender, DEMO_age, DEMO_dom_hand, DEMO_born_U.S.,  DEMO_U.S._citizenship,  DEMO_no_lived_U.S., DEMO_lang_at_home, DEMO_lang_at_home, DEMO_ethnic_all, DEMO_gen_status, DEMO_relationship, DEMO_employ_status, DEMO_housing_all, DEMO_no_roommates, DEMO_educ_mother, DEMO_educ_father

#Subsetting the data to include only PID and Baseline composite data 
df_eq_baseID <- dplyr::select(df_eq_baseID, PID, SSS_pre, UCLA_pre, DEMO_gender, DEMO_age, DEMO_dom_hand, DEMO_born_U.S.,   DEMO_U.S._citizenship,  DEMO_no_lived_U.S., DEMO_lang_at_home, DEMO_lang_at_home, DEMO_ethnic_all, DEMO_gen_status, DEMO_relationship, DEMO_employ_status, DEMO_housing_all, DEMO_no_roommates, DEMO_educ_mother, DEMO_educ_father)

write.csv(df_eq_baseID, "df_eq_baseID", row.names = F)

#Merging pre-IDs dataset into GPA dataset

df_eq_bdi <- df_eq_bdi %>% 
  mutate(PID = as.numeric(PID))
         
df_eq_baseID <- df_eq_baseID %>% 
  mutate(PID = as.numeric(PID))
         
df_eq_bdi <- left_join(df_eq_bdi, df_eq_baseID, by = "PID")

Tests to assess group differences at baseline.

Independent samples t-tests to assess group differences at baseline - bdi_pre

library(MASS)       # load the MASS package 

model_ID_bdi <- t.test(bdi_pre ~ group_fct, var.equal=TRUE, data=df_eq_bdi) # t.test to show BDI pre randomized sucessfully 

model_ID_bdi

## 
##  Two Sample t-test
## 
## data:  bdi_pre by group_fct
## t = -0.23711, df = 111, p-value = 0.813
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.596866  2.828069
## sample estimates:
## mean in group Group 1 mean in group Group 2 
##              13.08929              13.47368

dstats <- df_eq_bdi %>% group_by(group_fct) %>% summarise(mean = mean(bdi_pre), sd = sd(bdi_pre), n = n()) # code to generate the mean, standard deviation and sample size for each of your two groups

dstats

## # A tibble: 2 x 4
##   group_fct  mean    sd     n
##   <fct>     <dbl> <dbl> <int>
## 1 Group 1    13.1  8.10    56
## 2 Group 2    13.5  9.10    57

mean_diff <- dstats[1,2] - dstats[2,2] #code to generate the mean difference between bdi_pre for two groups

mean_diff

##         mean
## 1 -0.3843985

Write-Up: Group differences at baseline - bdi_pre

We performed an independent samples t-test to determine if there is a mean difference between Group 1 and Group 2 in terms of the scores on pre-BDI composite scores (to assess group differences at baseline). There was no significant difference between Group 1 (M = 13.09, SD =8.09, n=56) and Group 2 Group 1 (M = 13.47, SD =9.10, n=57), t(111) = -0.24, p = 0.81, d = -0.38, 95% CI [-3.60, 2.82]. Randomization was successful in terms of raw BDI scores at pre-intervention measurement.

Independent samples t-tests to assess group differences at baseline - SSS_pre

library(MASS)       # load the MASS package 

model_ID_SSS <- t.test(SSS_pre ~ group_fct, var.equal=TRUE, data=df_eq_bdi) # t.test to show SSS randomized sucessfully 

model_ID_SSS

## 
##  Two Sample t-test
## 
## data:  SSS_pre by group_fct
## t = -0.23392, df = 111, p-value = 0.8155
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.9673074  0.7630467
## sample estimates:
## mean in group Group 1 mean in group Group 2 
##              4.178571              4.280702

dstats <- df_eq_bdi %>% group_by(group_fct) %>% summarise(mean = mean(SSS_pre), sd = sd(SSS_pre), n = n()) # code to generate the mean, standard deviation and sample size for each of your two groups

dstats

## # A tibble: 2 x 4
##   group_fct  mean    sd     n
##   <fct>     <dbl> <dbl> <int>
## 1 Group 1    4.18  2.18    56
## 2 Group 2    4.28  2.45    57

mean_diff <- dstats[1,2] - dstats[2,2] #code to generate the mean difference between bdi_pre for two groups

mean_diff

##         mean
## 1 -0.1021303

Write-Up: Group differences at baseline - SSS_pre

We performed an independent samples t-test to determine if there is a mean difference between Group 1 and Group 2 in terms of the scores on subjective social standing (to assess group differences at baseline). There was no significant difference between Group 1 (M = 4.18 SD = 2.18, n=56) and Group 2 Group 1 (M = 4.28, SD =2.45, n=57), t(111) = -0.24, p = 0.82, d = -0.10, 95% CI [-0.97, 0.76]. Randomization was successful in terms of subjective social standing at pre-intervention measurement.

1-sample t-test / chi-squared test to assess group differences at baseline - DEMO_gender

library(MASS)       # load the MASS package 

model_ID_gender <- table(df_eq_bdi$DEMO_gender, df_eq_bdi$group_fct) 

#tbl                 # the contingency table 

chisq.test(model_ID_gender)

## 
##  Pearson's Chi-squared test
## 
## data:  model_ID_gender
## X-squared = 2.6004, df = 4, p-value = 0.6268

#randomized sucessfully

Write-Up: Group differences at baseline - gender

We performed a 1-sample t-test / chi-squared test on gender to assess group differences at baseline. Randomization at baseline was successful, X2 (4, N = 113) = 2.60, p = 0.63.

Independent samples t-tests to assess group differences at baseline - DEMO_age

library(MASS)       # load the MASS package 

model_ID_age <- t.test(DEMO_age ~ group_fct, var.equal=TRUE, data=df_eq_bdi) # t.test to show age randomized sucessfully 

model_ID_age

## 
##  Two Sample t-test
## 
## data:  DEMO_age by group_fct
## t = -0.027178, df = 99, p-value = 0.9784
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1451128  0.1411913
## sample estimates:
## mean in group Group 1 mean in group Group 2 
##              18.09804              18.10000

dstats <- df_eq_bdi %>% group_by(group_fct) %>% summarise(mean = mean(DEMO_age), sd = sd(DEMO_age), n = n()) # code to generate the mean, standard deviation and sample size for each of your two groups

dstats

## # A tibble: 2 x 4
##   group_fct  mean    sd     n
##   <fct>     <dbl> <dbl> <int>
## 1 Group 1      NA    NA    56
## 2 Group 2      NA    NA    57

mean_diff <- dstats[1,2] - dstats[2,2] #code to generate the mean difference between bdi_pre for two groups

mean_diff

##   mean
## 1   NA

Write-Up: Group differences at baseline - DEMO_age

We performed an independent samples t-test to determine if there is a mean difference between Group 1 and Group 2 in terms of age (to assess group differences at baseline). There was no significant difference between Group 1 (M = 18.10 SD = , n=56) and Group 2 Group 1 (M = 18.10, SD =, n=57), t(99) = -0.02, p = 0.98, d = 0.00, 95% CI [-0.15, 0.14]. Randomization was successful in terms of age at pre-intervention measurement.

Imputing missing scores for missing BDI post scores

df_eq_bdi$bdi_post_imp <- ifelse(is.na(df_eq_bdi$bdi_post), mean(df_eq_bdi$bdi_post, na.rm=TRUE), df_eq_bdi$bdi_post)

Adding BDI change score

df_eq_bdi$change_bdi <- df_eq_bdi$bdi_pre - df_eq_bdi$bdi_post_imp

Visualizing the data - outliers included

Histograms and boxplots of Pre-BDI - outliers included

library(ggplot2)

hist_bdipre <- df_eq_bdi %>%
  ggplot(aes(x = bdi_pre, fill = group_fct)) +
    geom_density(lwd    = 1, alpha = .3)

hist_bdipre

#boxplot(bdi_pre~Condition_Group,data=df_eq_bdi, main="Pre BDI Scores by Condition", xlab="Condition", ylab="BDI")

boxplot(df_eq_bdi$bdi_pre, main="Distribution of Pre BDI Scores",
   ylab="BDI")

Histograms and boxplots of Post-BDI - outliers included

hist_bdipost <- df_eq_bdi %>%
  ggplot(aes(x = bdi_post_imp, fill = group_fct)) +
    geom_density(lwd    = 1, alpha = .3)

hist_bdipost

boxplot(bdi_post_imp~Condition_Group,data=df_eq_bdi, main="Post BDI Scores by Condition",
   xlab="Condition", ylab="BDI")

#We will run the models with and without these outliers

Outliers detected at pre and post. We will run analysis both with outliers and with outliers excluded.

Descriptive Statistics - Raw scores, outliers included

#Descriptives - Pre BDI
psych::describe(df_eq_bdi$bdi_pre)

##    vars   n  mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 113 13.28 8.58     12    12.4 7.41   1  40    39 0.95     0.52 0.81

#Descriptives - Pre BDI by condition
psych::describeBy(df_eq_bdi$bdi_pre, df_eq_bdi$group_fct, mat = TRUE)

##     item  group1 vars  n     mean       sd median  trimmed   mad min max range
## X11    1 Group 1    1 56 13.08929 8.095513     12 12.52174 7.413   1  36    35
## X12    2 Group 2    1 57 13.47368 9.098872     12 12.36170 7.413   2  40    38
##          skew   kurtosis       se
## X11 0.6421622 -0.1672216 1.081808
## X12 1.1211774  0.7192014 1.205175

#Descriptives - Post BDI 
psych::describe(df_eq_bdi$bdi_post_imp)

##    vars   n  mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 113 11.99 8.22     10    11.1 7.41   0  46    46 1.32     2.45 0.77

#Descriptives - Post BDI by condition
psych::describeBy(df_eq_bdi$bdi_post_imp, df_eq_bdi$group_fct, mat = TRUE)

##     item  group1 vars  n     mean       sd  median trimmed      mad min max
## X11    1 Group 1    1 56 11.24929 7.415513 10.0000 10.5861 7.413000   0  36
## X12    2 Group 2    1 57 12.71791 8.943896 11.9901 11.7430 7.427679   1  46
##     range      skew  kurtosis        se
## X11    36 0.9309647 0.6842241 0.9909396
## X12    45 1.4321965 2.6312880 1.1846480

#Descriptives - Post log-transformed BDI by condition
#psych::describeBy(df_eq_bdi$bdi_post_log, df_eq_bdi$group_fct, mat = TRUE)

Log transforming BDI scores at pre and post

#summary(df_eq_bdi$bdi_pre)

#summary(df_eq_bdi$bdi_post_imp)

df_eq_bdi$bdi_pre_log <- round(log10(df_eq_bdi$bdi_pre+1), digits=4)

hist_bdiprelog <- df_eq_bdi %>%
  ggplot(aes(x = bdi_pre_log, fill = Condition_Group)) +
    geom_density(lwd    = 1, alpha = .3)

hist_bdiprelog

df_eq_bdi$bdi_post_log <- round(log10(df_eq_bdi$bdi_post_imp+1), digits=4)
                                
hist_bdipostlog <- df_eq_bdi %>%
  ggplot(aes(x = bdi_post_log, fill = Condition_Group)) +
    geom_density(lwd    = 1, alpha = .3)

hist_bdipostlog

#### Descriptive Statistics - Log-transformed scores, outliers included

# Descriptives - Log-transformed Pre BDI

psych::describe(df_eq_bdi$bdi_pre_log)

##    vars   n mean   sd median trimmed  mad min  max range skew kurtosis   se
## X1    1 113 1.07 0.29   1.11    1.09 0.28 0.3 1.61  1.31 -0.5    -0.06 0.03

# Descriptives - Log-transformed Pre BDI by condition

psych::describeBy(df_eq_bdi$bdi_pre_log, df_eq_bdi$group_fct, mat = TRUE)

##     item  group1 vars  n     mean        sd median  trimmed       mad    min
## X11    1 Group 1    1 56 1.062841 0.3020580 1.1139 1.086859 0.2570828 0.3010
## X12    2 Group 2    1 57 1.078902 0.2735009 1.1139 1.080651 0.3088256 0.4771
##        max  range        skew   kurtosis         se
## X11 1.5682 1.2672 -0.78576092  0.1796165 0.04036420
## X12 1.6128 1.1357 -0.07356917 -0.7252075 0.03622608

#Descriptives - Log-transformed Post BDI

psych::describe(df_eq_bdi$bdi_post_log)

##    vars   n mean  sd median trimmed  mad min  max range  skew kurtosis   se
## X1    1 113 1.02 0.3   1.04    1.04 0.29   0 1.67  1.67 -0.55     0.39 0.03

#Descriptives - Log-transformed Post BDI by condition

psych::describeBy(df_eq_bdi$bdi_post_log, df_eq_bdi$group_fct, mat = TRUE)

##     item  group1 vars  n     mean        sd median  trimmed       mad   min
## X11    1 Group 1    1 56 1.001477 0.2973139 1.0414 1.018735 0.3040813 0.000
## X12    2 Group 2    1 57 1.047349 0.2954527 1.1136 1.061936 0.2778392 0.301
##        max  range       skew   kurtosis         se
## X11 1.5682 1.5682 -0.6869758  0.7189888 0.03973025
## X12 1.6721 1.3711 -0.4064591 -0.1345192 0.03913366

Identifying and obtaining outliers

Identifying value of outliers

library(ggplot2)

hist_bdipre <- df_eq_bdi %>%
  ggplot(aes(x = bdi_pre, fill = group_fct)) +
    geom_density(lwd    = 1, alpha = .3)

hist_bdipre

box_pre <- boxplot(bdi_pre~Condition_Group,data=df_eq_bdi, main="Pre BDI Scores by Condition",
   xlab="Condition", ylab="BDI")

#Obtaining values of outliers

library(dplyr)
library(ggplot2)

#At PRE

is_outlier <- function(x) {
  return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}

df_eq_bdi %>%
  group_by(group_fct) %>%
  mutate(outlier = ifelse(is_outlier(bdi_pre), bdi_pre, as.numeric(NA))) %>%
  ggplot(., aes(x = factor(group_fct), y = bdi_pre)) +
    geom_boxplot() +
    geom_text(aes(label = outlier), na.rm = TRUE, hjust = -0.3)

#At POST
is_outlier <- function(x) {
  return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}

df_eq_bdi %>%
  group_by(group_fct) %>%
  mutate(outlier = ifelse(is_outlier(bdi_post_imp), bdi_post_imp, as.numeric(NA))) %>%
  ggplot(., aes(x = factor(group_fct), y = bdi_post_imp)) +
    geom_boxplot() +
    geom_text(aes(label = outlier), na.rm = TRUE, hjust = -0.3)

Creating dataframe with outliers removed

#Outliers at PRE
df_eq_bdi_out1 <- df_eq_bdi %>% 
  filter(group_fct == 'Group 1', bdi_pre < 35) 

df_eq_bdi_out2 <- df_eq_bdi %>% 
  filter(group_fct == 'Group 2', bdi_pre < 37) 

bdi_scale_out <- rbind(df_eq_bdi_out1, df_eq_bdi_out2)

#Outliers at POST
df_eq_bdipost_out1 <- bdi_scale_out %>% 
  filter(group_fct == 'Group 1', bdi_post_imp < 36) 

df_eq_bdipost_out2 <- bdi_scale_out %>% 
  filter(group_fct == 'Group 2', bdi_post_imp < 40) 

bdi_scale_out <- rbind(df_eq_bdipost_out1, df_eq_bdipost_out2)

#Check to see if pre and post are same subjects.

Log transforming BDI scores at pre and post - outliers removed

#summary(bdi_scale_out$bdi_pre)

#summary(bdi_scale_out$bdi_post_imp)

bdi_scale_out$bdi_pre_log_out <- round(log10(bdi_scale_out$bdi_pre+1), digits=4)

hist_bdiprelog_out <- bdi_scale_out %>%
  ggplot(aes(x = bdi_pre_log_out, fill = Condition_Group)) +
    geom_density(lwd    = 1, alpha = .3)

hist_bdiprelog_out

bdi_scale_out$bdi_post_log_out <- round(log10(bdi_scale_out$bdi_post_imp+1), digits=4)
                                
hist_bdipostlog_out <- bdi_scale_out %>%
  ggplot(aes(x = bdi_post_log_out, fill = Condition_Group)) +
    geom_density(lwd    = 1, alpha = .3)

hist_bdipostlog_out

Visualizing the data - outliers removed

Histograms and boxplots of Pre-BDI - outliers removed

library(ggplot2)

hist_bdipre_out <- bdi_scale_out %>%
  ggplot(aes(x = bdi_pre, fill = group_fct)) +
    geom_density(lwd    = 1, alpha = .3)

hist_bdipre_out

#boxplot(bdi_pre~Condition_Group,data=bdi_scale_out, main="Pre BDI Scores by Condition - Outliers removed", xlab="Condition", ylab="BDI")

boxplot(bdi_scale_out$bdi_pre, main="Distribution of Pre BDI Scores - Outliers removed",
   ylab="BDI")

Histograms and boxplots of Post-BDI - outliers removed

hist_bdipost_out <- bdi_scale_out %>%
  ggplot(aes(x = bdi_post_imp, fill = group_fct)) +
    geom_density(lwd    = 1, alpha = .3)

hist_bdipost_out

boxplot(bdi_post_imp~Condition_Group,data=bdi_scale_out, main="Post BDI Scores by Condition - Outliers Removed",
   xlab="Condition", ylab="BDI")

Descriptive Statistics - outliers removed

#Descriptives - Pre BDI - outliers removed
psych::describe(bdi_scale_out$bdi_pre)

##    vars   n  mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 108 12.41 7.56     12   11.74 7.41   1  36    35 0.77     0.08 0.73

#Descriptives - Pre BDI by condition - outliers removed
psych::describeBy(bdi_scale_out$bdi_pre, bdi_scale_out$group_fct, mat = TRUE)

##     item  group1 vars  n     mean       sd median  trimmed    mad min max range
## X11    1 Group 1    1 54 12.46296 7.447205   12.0 11.97727 7.4130   1  30    29
## X12    2 Group 2    1 54 12.35185 7.736601   10.5 11.47727 6.6717   2  36    34
##          skew   kurtosis       se
## X11 0.5062473 -0.5133451 1.013436
## X12 0.9841672  0.4850413 1.052818

#Descriptives - Post BDI - outliers removed
psych::describe(bdi_scale_out$bdi_post_imp)

##    vars   n  mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 108 11.16 6.64     10   10.71 5.93   0  28    28 0.57    -0.49 0.64

#Descriptives - Post BDI by condition - outliers removed
psych::describeBy(bdi_scale_out$bdi_post_imp, bdi_scale_out$group_fct, mat = TRUE)

##     item  group1 vars  n     mean       sd   median  trimmed    mad min max
## X11    1 Group 1    1 54 10.53630 6.435911  9.50000 10.13546 6.6717   0  28
## X12    2 Group 2    1 54 11.77631 6.847805 11.49505 11.38456 6.6717   1  27
##     range      skew   kurtosis        se
## X11    28 0.5704695 -0.5222228 0.8758166
## X12    26 0.5315128 -0.5957325 0.9318682

Descriptive Statistics - log-transformed, outliers excluded

#Descriptives - log-transformed Pre BDI - outliers removed
psych::describe(bdi_scale_out$bdi_pre_log_out)

##    vars   n mean   sd median trimmed  mad min  max range  skew kurtosis   se
## X1    1 108 1.05 0.28   1.11    1.07 0.28 0.3 1.57  1.27 -0.58        0 0.03

#Descriptives - log-transformed Pre BDI by condition - outliers removed
psych::describeBy(bdi_scale_out$bdi_pre_log_out, bdi_scale_out$group_fct, mat = TRUE)

##     item  group1 vars  n     mean        sd median  trimmed       mad    min
## X11    1 Group 1    1 54 1.047278 0.2958013 1.1139 1.072332 0.2367712 0.3010
## X12    2 Group 2    1 54 1.053950 0.2576072 1.0603 1.057541 0.2330647 0.4771
##        max  range       skew   kurtosis         se
## X11 1.4914 1.1904 -0.8318641  0.2057096 0.04025346
## X12 1.5682 1.0911 -0.1468762 -0.7271017 0.03505590

#Descriptives - log-transformed Post BDI - outliers removed
psych::describe(bdi_scale_out$bdi_post_log_out)

##    vars   n mean   sd median trimmed  mad min  max range  skew kurtosis   se
## X1    1 108 1.01 0.28   1.04    1.03 0.29   0 1.46  1.46 -0.77     0.57 0.03

#Descriptives - log-transformed Post BDI by condition - outliers removed
psych::describeBy(bdi_scale_out$bdi_post_log_out, bdi_scale_out$group_fct, mat = TRUE)

##     item  group1 vars  n      mean        sd median  trimmed       mad   min
## X11    1 Group 1    1 54 0.9833241 0.2865470 1.0207 1.004482 0.3109012 0.000
## X12    2 Group 2    1 54 1.0335537 0.2736819 1.0964 1.057116 0.2355851 0.301
##        max  range       skew    kurtosis         se
## X11 1.4624 1.4624 -0.8079165  0.87785646 0.03899410
## X12 1.4472 1.1462 -0.6951117 -0.01701206 0.03724339

Re-shaping the data to long-format

The interaction models will be conducted with the long-format version of the data.

# Re-shaping data to long format - outliers included 

##bdi <- dplyr::select(df_eq_bdi, PID, Condition_Group, group_fct, bdi_pre, bdi_post_imp, SSS_pre, UCLA_pre, DEMO_gender, DEMO_age, DEMO_dom_hand, DEMO_born_U.S.,  DEMO_U.S._citizenship,  DEMO_no_lived_U.S., DEMO_lang_at_home, DEMO_lang_at_home, DEMO_ethnic_all, DEMO_gen_status, DEMO_relationship, DEMO_employ_status, DEMO_housing_all, DEMO_no_roommates, DEMO_educ_mother, DEMO_educ_father)

bdi <- dplyr::select(df_eq_bdi, PID, Condition_Group, group_fct, bdi_pre, bdi_post_imp, change_bdi, SSS_pre, UCLA_pre, DEMO_gender, DEMO_age, DEMO_dom_hand, DEMO_born_U.S.,    DEMO_U.S._citizenship,  DEMO_no_lived_U.S., DEMO_lang_at_home, DEMO_lang_at_home, DEMO_ethnic_all, DEMO_gen_status, DEMO_relationship, DEMO_employ_status, DEMO_housing_all, DEMO_no_roommates, DEMO_educ_mother, DEMO_educ_father)
  
bdi_long <- bdi %>%
  tidyr::pivot_longer(cols = starts_with("bdi"), names_to = "time", names_prefix = NULL,
  names_sep = NULL, names_pattern = NULL, names_ptypes = list(),
  names_repair = "check_unique", values_to = "BDI_score",
  values_drop_na = FALSE, values_ptypes = list()) %>%
  mutate(time = as.factor(time))

bdi_long$time <- recode(bdi_long$time, bdi_pre="time 1_pre")
bdi_long$time <-recode(bdi_long$time, bdi_post_imp="time 2_post")

# Re-shaping data to long format - log-transformed + outliers included 

bdi_log <- dplyr::select(df_eq_bdi, PID, Condition_Group, group_fct, bdi_pre_log, bdi_post_log, SSS_pre, UCLA_pre, DEMO_gender, DEMO_age, DEMO_dom_hand, DEMO_born_U.S.,    DEMO_U.S._citizenship,  DEMO_no_lived_U.S., DEMO_lang_at_home, DEMO_lang_at_home, DEMO_ethnic_all, DEMO_gen_status, DEMO_relationship, DEMO_employ_status, DEMO_housing_all, DEMO_no_roommates, DEMO_educ_mother, DEMO_educ_father)
  
bdi_log_long <- bdi_log %>%
  tidyr::pivot_longer(cols = starts_with("bdi"), names_to = "time", names_prefix = NULL,
  names_sep = NULL, names_pattern = NULL, names_ptypes = list(),
  names_repair = "check_unique", values_to = "BDI_score",
  values_drop_na = FALSE, values_ptypes = list()) %>%
  mutate(time = as.factor(time))

bdi_long$time <- recode(bdi_long$time, bdi_pre="time 1_pre")
bdi_long$time <-recode(bdi_long$time, bdi_post_imp="time 2_post")

#bdi_out_long

bdi_log_out <- dplyr::select(bdi_scale_out, PID, Condition_Group, group_fct, bdi_pre_log, bdi_post_log, SSS_pre, UCLA_pre, DEMO_gender, DEMO_age, DEMO_dom_hand, DEMO_born_U.S.,    DEMO_U.S._citizenship,  DEMO_no_lived_U.S., DEMO_lang_at_home, DEMO_lang_at_home, DEMO_ethnic_all, DEMO_gen_status, DEMO_relationship, DEMO_employ_status, DEMO_housing_all, DEMO_no_roommates, DEMO_educ_mother, DEMO_educ_father)

# Re-shaping data to long format - outliers removed

##bdi_out <- dplyr::select(bdi_scale_out, PID, Condition_Group, group_fct, bdi_pre, bdi_post_imp, SSS_pre, UCLA_pre, DEMO_gender, DEMO_age, DEMO_dom_hand, DEMO_born_U.S.,  DEMO_U.S._citizenship,  DEMO_no_lived_U.S., DEMO_lang_at_home, DEMO_lang_at_home, DEMO_ethnic_all, DEMO_gen_status, DEMO_relationship, DEMO_employ_status, DEMO_housing_all, DEMO_no_roommates, DEMO_educ_mother, DEMO_educ_father)

bdi_out <- dplyr::select(bdi_scale_out, PID, Condition_Group, group_fct, bdi_pre, change_bdi, bdi_post_imp, SSS_pre, UCLA_pre, DEMO_gender, DEMO_age, DEMO_dom_hand, DEMO_born_U.S.,    DEMO_U.S._citizenship,  DEMO_no_lived_U.S., DEMO_lang_at_home, DEMO_lang_at_home, DEMO_ethnic_all, DEMO_gen_status, DEMO_relationship, DEMO_employ_status, DEMO_housing_all, DEMO_no_roommates, DEMO_educ_mother, DEMO_educ_father)
  
bdi_out_long <- bdi_out %>%
  tidyr::pivot_longer(cols = starts_with("bdi"), names_to = "time", names_prefix = NULL,
  names_sep = NULL, names_pattern = NULL, names_ptypes = list(),
  names_repair = "check_unique", values_to = "BDI_score",
  values_drop_na = FALSE, values_ptypes = list()) %>%
  mutate(time = as.factor(time))

bdi_out_long$time <- recode(bdi_out_long$time, bdi_pre="time 1_pre")
bdi_out_long$time <-recode(bdi_out_long$time, bdi_post_imp="time 2_post")

#bdi_out_long

##bdi_log_out <- dplyr::select(bdi_scale_out, PID, Condition_Group, group_fct, bdi_pre_log, bdi_post_log, SSS_pre, UCLA_pre, DEMO_gender, DEMO_age, DEMO_dom_hand, DEMO_born_U.S.,  DEMO_U.S._citizenship,  DEMO_no_lived_U.S., DEMO_lang_at_home, DEMO_lang_at_home, DEMO_ethnic_all, DEMO_gen_status, DEMO_relationship, DEMO_employ_status, DEMO_housing_all, DEMO_no_roommates, DEMO_educ_mother, DEMO_educ_father) #log transformed scores 

bdi_log_out <- dplyr::select(bdi_scale_out, PID, Condition_Group, group_fct, bdi_pre_log, bdi_post_log, change_bdi, SSS_pre, UCLA_pre, DEMO_gender, DEMO_age, DEMO_dom_hand, DEMO_born_U.S.,    DEMO_U.S._citizenship,  DEMO_no_lived_U.S., DEMO_lang_at_home, DEMO_lang_at_home, DEMO_ethnic_all, DEMO_gen_status, DEMO_relationship, DEMO_employ_status, DEMO_housing_all, DEMO_no_roommates, DEMO_educ_mother, DEMO_educ_father) #log transformed scores 
  
bdi_log_out_long <- bdi_log_out %>%
  tidyr::pivot_longer(cols = starts_with("bdi"), names_to = "time", names_prefix = NULL,
  names_sep = NULL, names_pattern = NULL, names_ptypes = list(),
  names_repair = "check_unique", values_to = "BDI_score",
  values_drop_na = FALSE, values_ptypes = list())

bdi_log_out_long$time <- recode(bdi_log_out_long$time, bdi_pre="time 1_pre")
bdi_log_out_long$time <-recode(bdi_log_out_long$time, bdi_post_imp="time 2_post")

INTERACTION MODEL - Raw data

We performed mixed effects models to model time and condition as fixed effects, in order to look at the condition main effect, time main effect, and time by condition interaction. BDI = b0 + b1Time + b2Condition + b3(Time∗Condition)

Running Mixed Model - raw data

#Load libraries
library(lme4)
library(lmerTest)

# Re-leveling time factor for graph - Q: but why does this change the p-value for "group_fctGroup 2"?
bdi_long <- within(data=bdi_long, time <- relevel(time, ref = 2))

model_rawa = lmer(BDI_score ~ group_fct*time + (1|PID), data=bdi_long)
#summary(model_rawa)
#coef(summary(model_rawa))

require(devtools)
#install_version("ggplot2", version = "3.3.0", repos = "http://cran.us.r-project.org")
#install.packages("jtools")
library(ggplot2)
library(jtools) # for summ()

summ(model_rawa)

Observations	226
Dependent variable	BDI_score
Type	Mixed effects linear regression

AIC	1572.00
BIC	1592.52
Pseudo-R² (fixed effects)	0.01
Pseudo-R² (total)	0.51

Fixed Effects
	Est.	S.E.	t val.	d.f.	p
(Intercept)	13.09	1.13	11.63	177.32	0.00
group_fctGroup 2	0.38	1.58	0.24	177.32	0.81
timetime 2_post	-1.84	1.12	-1.64	111.00	0.10
group_fctGroup 2:timetime 2_post	1.08	1.58	0.69	111.00	0.49
p values calculated using Satterthwaite d.f.

Random Effects
Group	Parameter	Std. Dev.
PID	(Intercept)	5.97
Residual		5.94

Grouping Variables
Group	# groups	ICC
PID	113	0.50

Graphing of Mixed Model - raw data

ggplot(fortify(model_rawa), aes(time, BDI_score, color=group_fct)) +
  stat_summary(fun.data=mean_se, geom="pointrange") +
  theme_classic() +
  labs(title="Beck Depression Inventory Scores (Raw)",x="Timepoint", y = "BDI Score")+
  stat_summary(fun.y=mean, geom="line", aes(group = group_fct)) +
  scale_color_manual(values=c('darkorange1','lightcyan4'))

Write-up of Mixed Model - raw data

Using a linear mixed effects model, we compared Group 1 and Group 2 on raw BDI scores (N=113). There was no significant main effect of condition (t(177) = .24, p =.81) time (t(111) = -1.63, p =.10) for raw BDI scores. With raw BDI, there was not a significant difference between change in BDI from time 1 (pre) to time 2 (post) (t(111) = .69, p =.49), specifically the Group 1 change score from pre to post was not different (m = -1.84, SE = ?) from the Group 2 (m = -0.75, SE = ?). Group 1 raw BDI scores at pre (m = 13.09, SE = 1.08) Group 2 raw BDI scores at pre (m = 13.47, SE = 1.21) Group 1 raw BDI scores at post (m = 11.25, SE = .99) Group 2 raw BDI scores at post (m = 12.72, SE = 1.19)

INTERACTION MODEL - Log-transformed

Running Mixed Model - Log-transformed

# Re-leveling time factor for graph 
bdi_log_long <- within(data=bdi_log_long, time <- relevel(time, ref = 2))

model_loga = lmer(BDI_score ~ group_fct*time + (1|PID), data=bdi_log_long)
#summary(model_loga)
#coef(summary(model_loga))
summ(model_loga)

Observations	226
Dependent variable	BDI_score
Type	Mixed effects linear regression

AIC	88.64
BIC	109.16
Pseudo-R² (fixed effects)	0.01
Pseudo-R² (total)	0.44

Fixed Effects
	Est.	S.E.	t val.	d.f.	p
(Intercept)	1.06	0.04	27.22	186.69	0.00
group_fctGroup 2	0.02	0.05	0.29	186.69	0.77
timebdi_post_log	-0.06	0.04	-1.48	111.00	0.14
group_fctGroup 2:timebdi_post_log	0.03	0.06	0.51	111.00	0.61
p values calculated using Satterthwaite d.f.

Random Effects
Group	Parameter	Std. Dev.
PID	(Intercept)	0.19
Residual		0.22

Grouping Variables
Group	# groups	ICC
PID	113	0.43

Graphing of Mixed Model - log-transformed

ggplot(fortify(model_loga), aes(time, BDI_score, color=group_fct)) +
  stat_summary(fun.data=mean_se, geom="pointrange") +
  theme_classic() +
  labs(title="Beck Depression Inventory Scores (Log-transformed)",x="Timepoint", y = "BDI Score")+
  stat_summary(fun.y=mean, geom="line", aes(group = group_fct)) +
  scale_color_manual(values=c('darkorange1','lightcyan4'))

Write-Up for - Log-transformed

BDI scores were positively skewed so we also created a log transformed BDI variable. Using a linear mixed effects model, we compared Group 1 and Group 2 on log-transformed BDI scores (N=113). There was no significant main effect of condition (t(187) =.83, p =.41) time (t(111) = 1.48, p = 0.14) for log-transformed BDI scores. With log-transformed BDI, there was not a significant difference between change in BDI from time 1 (pre) to time 2 (post) (t(111) = -0.51, p =0.61), specifically the Group 1 change score from pre to post was not different (m = , SE = ?) from the Group 2 (m = , SE = ?). Group 1 raw BDI scores at pre (m = , SE = ) Group 2 raw BDI scores at pre (m = , SE = ) Group 1 raw BDI scores at post (m = , SE = ) Group 2 raw BDI scores at post (m = , SE = )

INTERACTION MODEL - Outliers removed

Running Mixed Model - outliers removed

# Re-leveling time factor for graph 
bdi_out_long <- within(data=bdi_out_long, time <- relevel(time, ref = 2))

lmeModel_out = lmer(BDI_score ~ group_fct*time + (1|PID), data=bdi_out_long)
#summary(lmeModel_out)
#coef(summary(lmeModel_out))
summ(lmeModel_out)

Observations	216
Dependent variable	BDI_score
Type	Mixed effects linear regression

AIC	1437.53
BIC	1457.78
Pseudo-R² (fixed effects)	0.01
Pseudo-R² (total)	0.47

Fixed Effects
	Est.	S.E.	t val.	d.f.	p
(Intercept)	12.46	0.97	12.84	174.96	0.00
group_fctGroup 2	-0.11	1.37	-0.08	174.96	0.94
timetime 2_post	-1.93	1.01	-1.91	106.00	0.06
group_fctGroup 2:timetime 2_post	1.35	1.43	0.95	106.00	0.35
p values calculated using Satterthwaite d.f.

Random Effects
Group	Parameter	Std. Dev.
PID	(Intercept)	4.84
Residual		5.24

Grouping Variables
Group	# groups	ICC
PID	108	0.46

Graphing of Mixed Model - Outliers Removed

ggplot(fortify(lmeModel_out), aes(time, BDI_score, color=group_fct)) +
  stat_summary(fun.data=mean_se, geom="pointrange") +
  theme_classic() +
  labs(title="Beck Depression Inventory Scores (Outliers Removed)",x="Timepoint", y = "BDI Score")+
  stat_summary(fun.y=mean, geom="line", aes(group = group_fct)) +
  scale_color_manual(values=c('darkorange1','lightcyan4'))

Write-Up for Mixed Model - outliers removed

Using a linear mixed effects model, we compared Group 1 and Group 2 on raw BDI scores with outliers removed (N=108). There was no significant main effect of condition (t(175) = -0.08, p=0.94) time (t(106) = -1/91, p = 0.06) for raw BDI scores. With raw BDI, there was not a significant difference between change in BDI from time 1 (pre) to time 2 (post) (t(106) = 0.95, p = 0.35), specifically the Group 1 change score from pre to post was not different (m = , SE = ?) from the Group 2 (m = , SE = ?) when outliers were removed. Group 1 raw BDI scores at pre (m = 12.46, SE = 1.01) Group 2 raw BDI scores at pre (m = 12.35, SE = 1.05) Group 1 raw BDI scores at post (m = 10.54, SE = 0.88) Group 2 raw BDI scores at post (m = 11.78, SE = 0.93)

INTERACTION MODEL - Log-transformed + outliers removed

Running Mixed Model - Log-transformed + outliers removed

# Re-leveling time factor for graph 
bdi_log_out_long <- within(data=bdi_log_out_long, time <- relevel(as.factor(time), ref = 2)) 

lmeModel2_log = lmer(BDI_score ~ group_fct*time + (1|PID), data=bdi_log_out_long)
#summary(lmeModel2_log)
#coef(summary(lmeModel2_log))
summ(lmeModel2_log)

Observations	216
Dependent variable	BDI_score
Type	Mixed effects linear regression

AIC	67.84
BIC	88.10
Pseudo-R² (fixed effects)	0.01
Pseudo-R² (total)	0.42

Fixed Effects
	Est.	S.E.	t val.	d.f.	p
(Intercept)	1.05	0.04	27.61	180.71	0.00
group_fctGroup 2	0.01	0.05	0.12	180.71	0.90
timebdi_post_log	-0.06	0.04	-1.56	106.00	0.12
group_fctGroup 2:timebdi_post_log	0.04	0.06	0.75	106.00	0.45
p values calculated using Satterthwaite d.f.

Random Effects
Group	Parameter	Std. Dev.
PID	(Intercept)	0.18
Residual		0.21

Grouping Variables
Group	# groups	ICC
PID	108	0.42

Graphing of Mixed Model - Log-transformed + outliers removed

ggplot(fortify(lmeModel2_log), aes(time, BDI_score, color=group_fct)) +
  stat_summary(fun.data=mean_se, geom="pointrange") +
  theme_classic() +
  labs(title="Beck Depression Inventory Scores (Outliers Removed & Log-Transformed)",x="Timepoint", y = "BDI Score")+
  stat_summary(fun.y=mean, geom="line", aes(group = group_fct)) +
  scale_color_manual(values=c('darkorange1','lightcyan4'))

Write-Up for - Log-transformed + outliers removed

BDI scores were positively skewed so we also created a log transformed BDI variable. Using a linear mixed effects model, we compared Group 1 and Group 2 on log-transformed BDI scores with outliers removed (N=108). There was no significant main effect of condition (t(181) = .94, p =.35) time (t(106) = 1.56, p =.12) for raw BDI scores. With log-transformed BDI, there was not a significant difference between change in BDI from time 1 (pre) to time 2 (post) (t(106) = -0.75, p =.45), specifically the Group 1 change score from pre to post was not different (m = , SE = ?) from the Group 2 (m = , SE = ?) when outliers were removed. Group 1 raw BDI scores at pre (m = , SE = ) Group 2 raw BDI scores at pre (m = , SE = ) Group 1 raw BDI scores at post (m = , SE = ) Group 2 raw BDI scores at post (m = , SE = )

MODELS CONTROLLING FOR INDIVIDUAL DIFFERENCES

Saying-is-beliving audio completion

Two-way interaction with saying-is-believing audio completion as an interaction term (group x time x compliance), where saying-is-believing audio completion is a dichotomous variable (yes = complied, no = didn’t comply) - outliers removed

Does the two-way group x time interaction differ depending on if you got the full dose (i.e. completed and returned the saying-is-believing audio) or not?

#Merging df with Missing Data
df_bdi_total <- merge(bdi_out_long,df_eq_bdi_control,by="PID")

Model for Audio completed

# Re-leveling time factor for graph 
df_bdi_total <- within(data=df_bdi_total, time <- relevel(time, ref = 2))

lmeModel_comply <- lmer(BDI_score ~ group_fct*time + 
                          group_fct*Audio_Received + 
                          group_fct*time*Audio_Received + #to see if the effect of the intervention across time as a function of                           group varies as a function of whether or not audio was received 
                          (1|PID), data=df_bdi_total)
#summary(lmeModel_comply)
#coef(summary(lmeModel_comply))
summ(lmeModel_comply)

Observations	216
Dependent variable	BDI_score
Type	Mixed effects linear regression

AIC	1423.74
BIC	1457.49
Pseudo-R² (fixed effects)	0.03
Pseudo-R² (total)	0.48

Fixed Effects
	Est.	S.E.	t val.	d.f.	p
(Intercept)	9.91	2.80	3.54	202.58	0.00
group_fctGroup 2	3.09	4.04	0.76	188.29	0.45
timetime 1_pre	-1.67	3.00	-0.56	104.85	0.58
Audio_Received	0.70	2.94	0.24	205.52	0.81
group_fctGroup 2:timetime 1_pre	4.00	4.24	0.94	104.85	0.35
group_fctGroup 2:Audio_Received	-2.08	4.26	-0.49	190.83	0.63
timetime 1_pre:Audio_Received	3.95	3.18	1.24	104.85	0.22
group_fctGroup 2:timetime 1_pre:Audio_Received	-5.73	4.50	-1.27	104.85	0.21
p values calculated using Satterthwaite d.f.

Random Effects
Group	Parameter	Std. Dev.
PID	(Intercept)	4.87
Residual		5.20

Grouping Variables
Group	# groups	ICC
PID	107	0.47

Subsetting the data into saying-is-believing audio completed and saying-is-believing audio not completed in order to create two different two-by-two graphs via running individual 2x2 interactions

# Graph 1 - Subsetting for saying-is-believing audio completed
df_bdi_total_audio <- subset(df_bdi_total, Audio_Received_dich=="yes",)

# Re-leveling time factor for graph 
df_bdi_total_audio <- within(data=df_bdi_total_audio, time <- relevel(time, ref = 2))

lmeModel_comply_audio <- lmer(BDI_score ~ group_fct*time + (1|PID), data=df_bdi_total_audio)
summ(lmeModel_comply_audio)

Observations	192
Dependent variable	BDI_score
Type	Mixed effects linear regression

AIC	1275.59
BIC	1295.14
Pseudo-R² (fixed effects)	0.01
Pseudo-R² (total)	0.46

Fixed Effects
	Est.	S.E.	t val.	d.f.	p
(Intercept)	12.90	1.01	12.74	155.99	0.00
group_fctGroup 2	-0.73	1.45	-0.50	155.99	0.62
timetime 2_post	-2.29	1.06	-2.16	94.00	0.03
group_fctGroup 2:timetime 2_post	1.73	1.51	1.14	94.00	0.26
p values calculated using Satterthwaite d.f.

Random Effects
Group	Parameter	Std. Dev.
PID	(Intercept)	4.77
Residual		5.24

Grouping Variables
Group	# groups	ICC
PID	96	0.45

Graphing of Mixed Model - Audio completed

ggplot(fortify(lmeModel_comply_audio), aes(time, BDI_score, color=group_fct)) +
  stat_summary(fun.data=mean_se, geom="pointrange") +
  theme_classic() +
  labs(title="Beck Depression Inventory Scores (Outliers Removed, Audio Completed)",x="Timepoint", y = "BDI Score")+
  stat_summary(fun.y=mean, geom="line", aes(group = group_fct)) +
  scale_color_manual(values=c('darkorange1','lightcyan4'))

Model for Audio NOT completed

# Graph 2 - Subsetting for saying-is-believing audio NOT completed
df_bdi_total_NOaudio <- subset(df_bdi_total, Audio_Received_dich=="no",)

df_bdi_total_NOaudio <- within(df_bdi_total_NOaudio, time <- relevel(time, ref = 2))

lmeModel_comply_NOaudio <- lmer(BDI_score ~ group_fct*time + (1|PID), data=df_bdi_total_NOaudio)
summ(lmeModel_comply_NOaudio)

Observations	24
Dependent variable	BDI_score
Type	Mixed effects linear regression

AIC	152.66
BIC	159.73
Pseudo-R² (fixed effects)	0.15
Pseudo-R² (total)	0.63

Fixed Effects
	Est.	S.E.	t val.	d.f.	p
(Intercept)	7.50	3.07	2.45	15.12	0.03
group_fctGroup 2	7.83	4.34	1.81	15.12	0.09
timetime 2_post	1.67	2.85	0.59	10.00	0.57
group_fctGroup 2:timetime 2_post	-4.00	4.03	-0.99	10.00	0.34
p values calculated using Satterthwaite d.f.

Random Effects
Group	Parameter	Std. Dev.
PID	(Intercept)	5.66
Residual		4.93

Grouping Variables
Group	# groups	ICC
PID	12	0.57

Graphing of Mixed Model - Audio NOT completed

ggplot(fortify(lmeModel_comply_NOaudio), aes(time, BDI_score, color=group_fct)) +
  stat_summary(fun.data=mean_se, geom="pointrange") +
  theme_classic() +
  labs(title="Beck Depression Inventory Scores (Outliers Removed, Audio NOT Completed)",x="Timepoint", y = "BDI Score")+
  stat_summary(fun.y=mean, geom="line", aes(group = group_fct)) +
  scale_color_manual(values=c('darkorange1','lightcyan4'))

Interaction model controlling for gender - outliers removed

# Re-leveling time factor for graph 
#bdi_out_long <- within(data=bdi_out_long, time <- relevel(time, ref = 2))

lmeModel_gender<- lmer(BDI_score ~ group_fct*time + (1|PID) + DEMO_gender, data=bdi_out_long)
#summary(lmeModel_gender)
#coef(summary(lmeModel_gender))

summ(lmeModel_gender)

Observations	216
Dependent variable	BDI_score
Type	Mixed effects linear regression

AIC	1418.00
BIC	1451.75
Pseudo-R² (fixed effects)	0.09
Pseudo-R² (total)	0.48

Fixed Effects
	Est.	S.E.	t val.	d.f.	p
(Intercept)	11.72	2.00	5.85	115.76	0.00
group_fctGroup 2	-0.40	1.35	-0.30	172.04	0.77
timetime 2_post	-1.93	1.01	-1.91	106.00	0.06
DEMO_genderFemale	1.82	2.00	0.91	102.00	0.37
DEMO_genderMale	-1.28	2.10	-0.61	102.00	0.54
DEMO_genderNon-binary	8.22	3.86	2.13	102.00	0.04
DEMO_genderTransgender	8.24	6.17	1.34	102.00	0.18
group_fctGroup 2:timetime 2_post	1.35	1.43	0.95	106.00	0.35
p values calculated using Satterthwaite d.f.

Random Effects
Group	Parameter	Std. Dev.
PID	(Intercept)	4.54
Residual		5.24

Grouping Variables
Group	# groups	ICC
PID	108	0.43

Graphing of Mixed Model - controlling for gender - outliers removed

ggplot(fortify(lmeModel_gender), aes(time, BDI_score, color=group_fct)) +
  stat_summary(fun.data=mean_se, geom="pointrange") +
  theme_classic() +
  labs(title="Beck Depression Inventory Scores (Outliers Removed, Controlling for Gender)",x="Timepoint", y = "BDI Score")+
  stat_summary(fun.y=mean, geom="line", aes(group = group_fct)) +
  scale_color_manual(values=c('darkorange1','lightcyan4'))

Interaction model controlling for SSS - outliers removed

# Re-leveling time factor for graph 
#bdi_out_long <- within(data=bdi_out_long, time <- relevel(time, ref = 2))

lmeModel_SSS <- lmer(BDI_score ~ group_fct*time + SSS_pre + (1|PID), data=bdi_out_long) #(1|PID) is the random effects term
#summary(lmeModel_SSS)
#coef(summary(lmeModel_SSS))
summ(lmeModel_SSS)

Observations	216
Dependent variable	BDI_score
Type	Mixed effects linear regression

AIC	1431.80
BIC	1455.43
Pseudo-R² (fixed effects)	0.07
Pseudo-R² (total)	0.47

Fixed Effects
	Est.	S.E.	t val.	d.f.	p
(Intercept)	15.55	1.40	11.13	135.78	0.00
group_fctGroup 2	-0.03	1.34	-0.02	177.31	0.98
timetime 2_post	-1.93	1.01	-1.91	106.00	0.06
SSS_pre	-0.73	0.24	-3.00	105.00	0.00
group_fctGroup 2:timetime 2_post	1.35	1.43	0.95	106.00	0.35
p values calculated using Satterthwaite d.f.

Random Effects
Group	Parameter	Std. Dev.
PID	(Intercept)	4.56
Residual		5.24

Grouping Variables
Group	# groups	ICC
PID	108	0.43

Graphing of Mixed Model - controlling for SSS - outliers removed

ggplot(fortify(lmeModel_SSS), aes(time, BDI_score, color=group_fct)) +
  stat_summary(fun.data=mean_se, geom="pointrange") +
  theme_classic() +
  labs(title="Beck Depression Inventory Scores (Outliers Removed, Controlling for Subjective Social Status)",x="Timepoint", y = "BDI Score")+
  stat_summary(fun.y=mean, geom="line", aes(group = group_fct)) +
  scale_color_manual(values=c('darkorange1','lightcyan4'))

Write-up…

The DV (BDI score)…Running an interaction model controlling for time (pre and post intervention) and condition; and accounting for subjective social standing (SSS_pre). As subjective social standing increases, BDI score (for both time and condition) significantly decreases (t(105) = -3.00, p = 0.00), such that for every unit increase in social standing, the expected BDI score decreases by 0.73 units.

Comparing models

model_comp_SSS <- anova(lmeModel_out,lmeModel_SSS)

#Note: this warning "refitting model(s) with ML (instead of REML)" means that R was not able to use REML - likely due to missing data. Because the same estimator has to be used for each model, it used ML instead.

model_comp_SSS

## Data: bdi_out_long
## Models:
## lmeModel_out: BDI_score ~ group_fct * time + (1 | PID)
## lmeModel_SSS: BDI_score ~ group_fct * time + SSS_pre + (1 | PID)
##              Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## lmeModel_out  6 1444.1 1464.4 -716.07   1432.1                            
## lmeModel_SSS  7 1437.2 1460.9 -711.62   1423.2 8.8923      1   0.002864 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Including the subjective social status estimator in the model (lmeModel_SSS) outperformed the model without including the subjective social status estimator (x2 (?) = 8.89, p = 0.00). However, this did not significantly impact the parameters of the interaction term, suggesting that subjective social status is not impacting the variance explained by the interaction model.

Three-Way interaction models to account for 1) researcher errors in sending microhits and 2) participant non-compliance with saying-is-believing audio

Three-way interactions reveal potential differences in lower-order two-way interactions. Here, we are testing whether the individual 2x2 interactions (above) are different from each other.

Three way-interaction (group x time x saying-is-believing audio completion), where saying-is-believing audio completion is a dichotomous variable (1 = complied, 0 = didn’t comply). - outliers removed

Tables of means - BDI_score (across pre and post)

# Q: Is the two-way interaction different depending on if you got the full dose or not?

library("apaTables")
#install.packages("here")
library("here")

#df_bdi_total <- df_bdi_total %>% mutate(time  = factor(time ))

# TABLE OF MEANS - BDI_score
table_means <- apa.2way.table(
               iv1 = group_fct, 
               iv2 = Audio_Received_dich, 
               #iv3 = time, # adding to look at the BDI score in the two groups across time (t1, t2)                 separately for the yes/no audio people.
               dv = BDI_score, 
               data = df_bdi_total,
               show.marginal.means = TRUE,
               #table.number = 1,
               #filename = here::here("/Users/mlipsett/Desktop/Equanimity/Database/means.doc")
               ) 

table_means

## 
## 
## Means and standard deviations for BDI_score as a function of a 2(group_fct) X 2(Audio_Received_dich) design 
## 
##            Audio_Received_dich                              
##                             no        yes      Marginal     
##  group_fct                   M   SD     M   SD        M   SD
##    Group 1                8.33 5.33 11.75 7.09    11.38 6.99
##    Group 2               14.17 8.74 11.89 7.10    12.15 7.29
##   Marginal               11.25 7.68 11.82 7.08              
## 
## Note. M and SD represent mean and standard deviation, respectively. 
## Marginal indicates the means and standard deviations pertaining to main effects.

# Merging wide-format df with Missing Data
wide_df_bdi_total <- merge(bdi,df_eq_bdi_control,by="PID")

table_means2 <- apa.2way.table(
               iv1 = group_fct, 
               iv2 = Audio_Received_dich, 
               dv = bdi_post_imp, 
               data = wide_df_bdi_total,
               show.marginal.means = TRUE,
               #table.number = 1,
               #filename = here::here("/Users/mlipsett/Desktop/Equanimity/Database/means.doc")
               ) 

table_means2

## 
## 
## Means and standard deviations for bdi_post_imp as a function of a 2(group_fct) X 2(Audio_Received_dich) design 
## 
##            Audio_Received_dich                              
##                             no        yes      Marginal     
##  group_fct                   M   SD     M   SD        M   SD
##    Group 1                9.17 7.47 11.39 7.41    11.16 7.38
##    Group 2               13.00 6.51 12.70 9.33    12.73 9.02
##   Marginal               11.08 6.97 12.04 8.40              
## 
## Note. M and SD represent mean and standard deviation, respectively. 
## Marginal indicates the means and standard deviations pertaining to main effects.

TABLE OF MEANS - BDI change score (Change score = BDI pre - BDI post so negative scores reflect increases in depression score from pre to post) x saying-is-believing audio completion)

table_means_changescore <- apa.2way.table(
              iv1 = group_fct, 
              iv2 = Audio_Received_dich, 
              dv = change_bdi, 
              data = df_bdi_total,
              show.marginal.means = TRUE,
              table.number = 3,
              filename = "Table3_APA.doc"
              ) 

table_means_changescore

## 
## 
## Table 3 
## 
## Means and standard deviations for change_bdi as a function of a 2(group_fct) X 2(Audio_Received_dich) design 
## 
##            Audio_Received_dich                             
##                             no       yes      Marginal     
##  group_fct                   M   SD    M   SD        M   SD
##    Group 1               -1.67 7.15 2.29 6.15     1.86 6.35
##    Group 2                2.33 6.11 0.55 8.46     0.76 8.22
##   Marginal                0.33 6.82 1.44 7.40              
## 
## Note. M and SD represent mean and standard deviation, respectively. 
## Marginal indicates the means and standard deviations pertaining to main effects.

# Merging wide-format df with Missing Data
wide_df_bdi_total <- merge(bdi,df_eq_bdi_control,by="PID")

table_means_changescore2 <- apa.2way.table(
               iv1 = group_fct, 
               iv2 = Audio_Received_dich, 
               dv = change_bdi, 
               data = wide_df_bdi_total,
               show.marginal.means = TRUE,
               table.number = 4,
               filename = "Table4_APA.doc"
               ) 

table_means_changescore2

## 
## 
## Table 4 
## 
## Means and standard deviations for change_bdi as a function of a 2(group_fct) X 2(Audio_Received_dich) design 
## 
##            Audio_Received_dich                              
##                             no       yes       Marginal     
##  group_fct                   M   SD    M    SD        M   SD
##    Group 1               -1.67 7.50 2.18  6.50     1.77 6.65
##    Group 2                2.33 6.41 0.76 10.17     0.93 9.80
##   Marginal                0.33 6.97 1.48  8.50              
## 
## Note. M and SD represent mean and standard deviation, respectively. 
## Marginal indicates the means and standard deviations pertaining to main effects.

Interpreting the table

Change score = BDI pre - BDI post so negative scores reflect increases in depression score from pre to post. For participants who did NOT turn in the saying-is-believing audi - group one had a moderate increase in depression (negative change score) (N = ?, M = -1.67), group 2 had a moderate decrease (positive change score) in depression (N = ?, M = 2.33). For participants who DID turn in their saying is believing audio - group one had a moderate decrease in depression (positive change score) (N = ?, M = 2.18), group 2 had a small decrease (positive change score) in depression (N = ?, M = 0.76).

Graphing the change scores - x saying-is-believing audio completion

# Re-leveling time factor for graph 
df_bdi_total <- within(data=df_bdi_total, time <- relevel(time, ref = 2))

# Change group_fct (is it still categorical?) to factor variables? 
df_bdi_total <- df_bdi_total %>% 
   mutate(group_fct = as.numeric(group_fct))

lmeModel_comply_3way <- lmer(BDI_score ~ group_fct*time*Audio_Received + (1|PID), data=df_bdi_total)

summ(lmeModel_comply_3way)

Observations	216
Dependent variable	BDI_score
Type	Mixed effects linear regression

AIC	1423.74
BIC	1457.49
Pseudo-R² (fixed effects)	0.03
Pseudo-R² (total)	0.48

Fixed Effects
	Est.	S.E.	t val.	d.f.	p
(Intercept)	1.15	6.31	0.18	197.63	0.86
group_fct	7.09	4.04	1.76	188.29	0.08
timetime 2_post	5.67	6.71	0.84	104.85	0.40
Audio_Received	12.46	6.64	1.88	200.88	0.06
group_fct:timetime 2_post	-4.00	4.24	-0.94	104.85	0.35
group_fct:Audio_Received	-7.81	4.26	-1.83	190.83	0.07
timetime 2_post:Audio_Received	-9.69	7.11	-1.36	104.85	0.18
group_fct:timetime 2_post:Audio_Received	5.73	4.50	1.27	104.85	0.21
p values calculated using Satterthwaite d.f.

Random Effects
Group	Parameter	Std. Dev.
PID	(Intercept)	4.87
Residual		5.20

Grouping Variables
Group	# groups	ICC
PID	107	0.47

# Create two different two-by-two graphs (1 for people who got the full dose (or some threshold) and 1 for people who did not)
interactions::probe_interaction(lmeModel_comply_3way, pred = group_fct, modx = time, mod2 = Audio_Received,
                 alpha = .1)

## ██████████████████████████ While Audio_Received (2nd moderator) = 0.00 (0) ██████████████████████████ 
## 
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of group_fct when time = time 1_pre: 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   7.09   4.04     1.76   0.08
## 
## Slope of group_fct when time = time 2_post: 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   3.09   4.04     0.76   0.45
## 
## ██████████████████████████ While Audio_Received (2nd moderator) = 1.00 (1) ██████████████████████████ 
## 
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of group_fct when time = time 1_pre: 
## 
##    Est.   S.E.   t val.      p
## ------- ------ -------- ------
##   -0.72   1.45    -0.50   0.62
## 
## Slope of group_fct when time = time 2_post: 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   1.01   1.45     0.70   0.49

Three-way interaction model controlling for number of microhits received - outliers removed

#Three way-interaction (group x time x number of microhits)
require(devtools)
#install_version("ggplot2", version = "3.3.0", repos = "http://cran.us.r-project.org")
#install.packages("jtools")
#install.packages("interactions")
library(ggplot2)
library(jtools) # for summ()
library(interactions)

# Re-leveling time factor for graph 
df_bdi_total <- within(data=df_bdi_total, time <- relevel(time, ref = 2))

lmeModel_micro <- lmer(BDI_score ~ group_fct*time*micro_prePOST + (1|PID), data=df_bdi_total)
summ(lmeModel_micro)

Observations	194
Dependent variable	BDI_score
Type	Mixed effects linear regression

AIC	1287.44
BIC	1320.12
Pseudo-R² (fixed effects)	0.02
Pseudo-R² (total)	0.50

Fixed Effects
	Est.	S.E.	t val.	d.f.	p
(Intercept)	9.24	8.36	1.10	150.13	0.27
group_fct	0.89	5.40	0.16	150.13	0.87
timetime 1_pre	8.06	8.45	0.95	93.00	0.34
micro_prePOST	0.02	3.62	0.01	150.13	1.00
group_fct:timetime 1_pre	-2.94	5.46	-0.54	93.00	0.59
group_fct:micro_prePOST	0.13	2.33	0.06	150.13	0.96
timetime 1_pre:micro_prePOST	-2.48	3.66	-0.68	93.00	0.50
group_fct:timetime 1_pre:micro_prePOST	0.94	2.36	0.40	93.00	0.69
p values calculated using Satterthwaite d.f.

Random Effects
Group	Parameter	Std. Dev.
PID	(Intercept)	5.08
Residual		5.19

Grouping Variables
Group	# groups	ICC
PID	97	0.49

# Create two different two-by-two graphs (1 for people who got the full dose (or some threshold) and 1 for people who did not)

interactions::probe_interaction(lmeModel_micro, pred = group_fct, modx = time, mod2 = micro_prePOST,
                  alpha = .1)

## ████████████████████████ While micro_prePOST (2nd moderator) = 1.59 (- 1 SD) ████████████████████████ 
## 
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of group_fct when time = time 2_post: 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   1.09   2.09     0.52   0.60
## 
## Slope of group_fct when time = time 1_pre: 
## 
##    Est.   S.E.   t val.      p
## ------- ------ -------- ------
##   -0.36   2.09    -0.17   0.86
## 
## █████████████████████████ While micro_prePOST (2nd moderator) = 2.23 (Mean) █████████████████████████ 
## 
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of group_fct when time = time 2_post: 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   1.17   1.48     0.79   0.43
## 
## Slope of group_fct when time = time 1_pre: 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   0.32   1.48     0.22   0.83
## 
## ████████████████████████ While micro_prePOST (2nd moderator) = 2.86 (+ 1 SD) ████████████████████████ 
## 
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of group_fct when time = time 2_post: 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   1.26   2.09     0.60   0.55
## 
## Slope of group_fct when time = time 1_pre: 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   1.00   2.09     0.48   0.63

INTERACTION MODELS - OPTION #2 - SHOULD I BE USING THIS APPROACH INSTEAD?

Repeated measures ANOVA w/ outliers removed

#model_BDI_pretopostint_out_long <-lm(BDI_score ~ group_fct + time + group_fct:time, data = bdi_out_long)  #outliers removed

#summary(model_BDI_pretopostint_out_long)

#tab_model(model_BDI_pretopostint_out_long) #outliers removed

#library(sjPlot)

#plot_model(model = model_BDI_pretopostint_out_long, type = "pred", terms = c("time", "group_fct"), order.terms = c(2,1)) # HOW DO I SWITCH PRE AND POST IN THE GRAPH?

#str(bdi_out_long)

# Running model with log-transformed data 
#model_BDI_log_pretopostint_out_long <-lm(BDI_score ~ group_fct + time + group_fct:time, data = bdi_log_out_long)  #outliers removed, log-transformed 

#tab_model(model_BDI_log_pretopostint_out_long) #outliers removed, log-transformed

# Write-Up: Repeated measures ANOVA with outliers removed. There is no significant difefrece between Group 1 and Group 2 (t(188)=-.34, p = .73).

Running a Repeated Measures Anova with condition as the between subjects effect and BDI score as the repeated measure (pre-post). - OPTION 2 (from this: https://m-clark.github.io/docs/mixedModels/anovamixed.html#pre-post_design)

#rm_anova_model <- aov(BDI_score ~ group_fct*time + Error(PID), bdi_out_long)
#summary(rm_anova_model)

Equanimity_BDI Depression

Megan Lipsett

2/25/2020