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(1, 13, 19, 28, 38, 44, 55, 56, 58, 64, 74, 76, 82, 88),]#Removing participants whose condition is unknown or who did not complete post-qualtrics measures

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 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. A chi-square test was performed for baseline composite BDI score to assess group differences at baseline (chi-squared = 26.968, df = 31, p-value = 0.6738), confirming that randomization was successful.

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'

Performing a 1-sample t-test / chi-squared test to find out if their are group differences at baseline

library(MASS)       # load the MASS package 

tbl = table(df_eq_bdi$bdi_pre, df_eq_bdi$group_fct) 

#tbl                 # the contingency table 

chisq.test(tbl)

## 
##  Pearson's Chi-squared test
## 
## data:  tbl
## X-squared = 26.968, df = 31, p-value = 0.6738

#randomized sucessfully

Write-Up: Group differences at baseline.

Performed a 1-sample t-test / chi-squared test on pre-BDI composite scores to assess group differences at baseline. Randmoization at baseline was successful, X2 (31, N = 101) = 26.97, p = 0.67.

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, fill = group_fct)) +
    geom_density(lwd    = 1, alpha = .3)

hist_bdipost

boxplot(bdi_post~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 101 13.43 8.67     12   12.51 7.41   1  40    39 0.94     0.54 0.86

#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 52 12.84615 8.134614     12 12.21429 7.4130   1  36    35
## X12    2 Group 2    1 49 14.04082 9.253286     13 12.92683 8.8956   2  40    38
##          skew    kurtosis       se
## X11 0.6737865 -0.09896884 1.128068
## X12 1.0709432  0.58996439 1.321898

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

##    vars   n  mean  sd median trimmed  mad min max range skew kurtosis   se
## X1    1 101 11.99 8.7     10   11.01 7.41   0  46    46 1.24     1.86 0.87

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

##     item  group1 vars  n     mean       sd median  trimmed    mad min max range
## X11    1 Group 1    1 52 11.19231 7.697819    8.5 10.45238 6.6717   0  36    36
## X12    2 Group 2    1 49 12.83673 9.655195   10.0 11.78049 7.4130   1  46    45
##          skew  kurtosis       se
## X11 0.9179824 0.4436953 1.067495
## X12 1.2881765 1.7590460 1.379314

#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)

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+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 101 1.07 0.29   1.11    1.09 0.28 0.3 1.61  1.31 -0.53    -0.01 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 52 1.051975 0.3074462 1.1139 1.076460 0.2773945 0.3010
## X12    2 Group 2    1 49 1.098147 0.2705688 1.1461 1.098898 0.2845109 0.4771
##        max  range       skew    kurtosis         se
## X11 1.5682 1.2672 -0.7406712  0.02835809 0.04263511
## X12 1.6128 1.1357 -0.1166097 -0.70305200 0.03865268

#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 101 1.01 0.31   1.04    1.03 0.32   0 1.67  1.67 -0.43     0.03 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 52 0.9928519 0.3070288 0.9771 1.009700 0.2948891 0.000
## X12    2 Group 2    1 49 1.0365327 0.3177893 1.0414 1.046202 0.3171281 0.301
##        max  range       skew   kurtosis         se
## X11 1.5682 1.5682 -0.5928620  0.4490832 0.04257724
## X12 1.6721 1.3711 -0.2813772 -0.5550168 0.04539847

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), bdi_post, as.numeric(NA))) %>%
  ggplot(., aes(x = factor(group_fct), y = bdi_post)) +
    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 < 36) 

df_eq_bdipost_out2 <- bdi_scale_out %>% 
  filter(group_fct == 'Group 2', bdi_post < 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)

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+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, fill = group_fct)) +
    geom_density(lwd    = 1, alpha = .3)

hist_bdipost_out

boxplot(bdi_post~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 96 12.45 7.54     12   11.79 7.41   1  36    35 0.75     0.12 0.77

#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 50 12.16000 7.410307     12 11.65000 7.413   1  30    29
## X12    2 Group 2    1 46 12.76087 7.755098     12 11.94737 7.413   2  36    34
##          skew   kurtosis       se
## X11 0.5216878 -0.4850062 1.047976
## X12 0.9368242  0.4603544 1.143427

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

##    vars  n  mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 96 11.05 7.04      9   10.56 5.93   0  28    28 0.58    -0.73 0.72

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

##     item  group1 vars  n     mean       sd median  trimmed    mad min max range
## X11    1 Group 1    1 50 10.42000 6.679515      8  9.95000 6.6717   0  28    28
## X12    2 Group 2    1 46 11.73913 7.430971     10 11.34211 6.6717   1  27    26
##          skew   kurtosis       se
## X11 0.6014985 -0.6536063 0.944626
## X12 0.5029581 -0.9547173 1.095637

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 96 1.05 0.28   1.11    1.07 0.28 0.3 1.57  1.27 -0.62     0.06 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 50 1.034732 0.3003285 1.1139 1.060830 0.2570828 0.3010
## X12    2 Group 2    1 46 1.070111 0.2534375 1.1139 1.073579 0.2773945 0.4771
##        max  range       skew    kurtosis         se
## X11 1.4914 1.1904 -0.7895186  0.04834864 0.04247286
## X12 1.5682 1.0911 -0.1971168 -0.69492372 0.03736733

#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 96    1 0.29      1    1.02 0.33   0 1.46  1.46 -0.62     0.16 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 50 0.972902 0.2954922 0.9542 0.993570 0.2845109 0.000
## X12    2 Group 2    1 46 1.019633 0.2947549 1.0414 1.038929 0.2910344 0.301
##        max  range       skew   kurtosis         se
## X11 1.4624 1.4624 -0.6987991  0.5945920 0.04178891
## X12 1.4472 1.1462 -0.5179689 -0.5091765 0.04345925

MODELS FOR MAIN EFFECT OF CONDITION

Running a simple regression of post BDI score on treatment (Main effect of condition) - outliers included

#OPTION 1: Two-Sample t-test
ttest_BDI_Model = t.test(bdi_post ~ group_fct, data=df_eq_bdi, var.equal=T)
#ttest_BDI_Model
#ttest_BDI_Model$statistic^2

#OPTION 2: Linear regression (same output as above)
#### Running a standard linear regression model 
lmModelPostOnly<- lm(formula = bdi_post ~ group_fct, data = df_eq_bdi)
#summary(lmModelPostOnly) #Note: The intercept is the mean score of BDI_post for Group 1. The difference between the means for the post BDI for Group 1 and Group 2 (i.e. the coefficient) is 1.64.

#coef(lmModelPostOnly)
#str(lmModelPostOnly)

broom::tidy(lmModelPostOnly)

## # A tibble: 2 x 5
##   term             estimate std.error statistic  p.value
##   <chr>               <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)         11.2       1.21     9.27  4.18e-15
## 2 group_fctGroup 2     1.64      1.73     0.949 3.45e- 1

library (sjPlot)
tab_model(lmModelPostOnly)

	bdi post
Predictors	Estimates	CI	p
(Intercept)	11.19	8.80 – 13.59	<0.001
group_fct [Group 2]	1.64	-1.79 – 5.08	0.345
Observations	101
R² / R² adjusted	0.009 / -0.001

# Primary Index model results
#lmModelPostOnly_results <- apa_print(lmModelPostOnly)

#lmModelPostOnly_results$table %>% 
  #apa_table(caption = "Model Regressing Follow-Up BDI Score on Condition",
            #note = "* p < 0.05; ** p < 0.01; *** p < 0.001",
            #escape = TRUE)

#### OPTION 3: ANOVA of post score on treatment - same output as above

#anovaModelPostOnly = aov(bdi_post ~ group_fct, df_eq_bdi)
#summary(anovaModelPostOnly)
#coef(anovaModelPostOnly)

# t2 from t.test = F-value (.90) of ANOVA, with identical p-values (.035). t2 from t.test = F-value (.90) of ANOVA, with identical p-values (.035). The t for the coefficient (and associated p-value) is identical to the t from the t-test, while the overall model F in the regression is the same as in the ANOVA. This equivalence is explicit in R, as the aov function actually calls the lm function, and the class of an aov object is both ‘aov’ and ‘lm’.

#Compared Group 1 and Group 2 on Post BDI composite score using an ANOVA. With raw BDI-post, there was not a significant difference between groups (F() = , p =), specifically Group 1 was not different (N=52, mean =11.19, sd =7.7, se =1.07) from Group 2 (N=49, mean =12.83, sd = 9.66, se = 1.38).

Write-Up: (Main effect of Condition - raw scores / outliers included) Simple regression of post BDI score on treatment

Compared Group 1 and Group 2 on Post BDI composite score using a two sample t-test. With raw BDI-post, there was not a significant difference between groups (t(99) = -0.95, p =.35), specifically Group 1 was not different (N=52, mean =11.19, sd =7.7, se =1.07) from Group 2 (N=49, mean =12.83, sd = 9.66, se = 1.38).Condition did not significantly predict depression outcomes.

Regression for log-transformed Post BDI across conditions

model_BDI_post_log <- lm(bdi_post_log ~ group_fct, data = df_eq_bdi) 
#Testing for main effect. (Group 1 will be the reference group.) 

#summary(model_BDI_post_log)

#str(model_BDI_post_log)

#broom::tidy(model_BDI_post_log)

library (sjPlot)
tab_model(model_BDI_post_log)

	bdi post log
Predictors	Estimates	CI	p
(Intercept)	0.99	0.91 – 1.08	<0.001
group_fct [Group 2]	0.04	-0.08 – 0.17	0.484
Observations	101
R² / R² adjusted	0.005 / -0.005

# Primary Index model results
#model_BDI_post_log_results <- apa_print(model_BDI_post)

#model_BDI_post_log_results$table %>% 
  #apa_table(caption = "Model Regressing Log-transformed Follow-Up BDI Score on Condition",
            #note = "* p < 0.05; ** p < 0.01; *** p < 0.001",
            #escape = TRUE)

Write-Up (Main effect of Condition - log-transformed scores / outliers included): Regression of time on log-transformed BDI composite scores

Compared Group 1 and Group 2 on log-transformed Post BDI composite score using a simple regression. With log-transformed BDI-post, there was not a significant difference between groups (t(99) = -0.95, p =.35). Specifically, Group 1 was not different (N=52, mean =.99 , sd =0.31 , se =0.04) from Group 2 (N=52, mean =.99 , sd =0.31 , se =0.04).Condition did not significantly predict depression outcomes.

Running a simple regression of post BDI score on treatment (Main effect of condition) - outliers removed

model_BDI_out <- lm(bdi_post ~ group_fct, data = bdi_scale_out) 

#summary(model_BDI_out)

library (sjPlot)
tab_model(model_BDI_out)

	bdi post
Predictors	Estimates	CI	p
(Intercept)	10.42	8.44 – 12.40	<0.001
group_fct [Group 2]	1.32	-1.54 – 4.18	0.362
Observations	96
R² / R² adjusted	0.009 / -0.002

Write-Up: (Main effect of Condition - raw scores / outliers removed) Simple regression of post BDI score on treatment

Removed outliers and Compared Group 1 and Group 2 on Post raw BDI composite score using a regression model. With raw BDI-post, there was not a significant difference between groups (t(94) = .92, p = .36). Specifically, Group 1 was not different (N=50, mean =10.42 , sd =6.68 , se =0.94) from Group 2 (N=46, mean =11.73 , sd =7.43 , se =1.10).Condition did not significantly predict depression outcomes.

Running a simple regression of log-transformed post BDI score on treatment (Main effect of condition) - outliers removed

model_BDI_out_log <- lm(bdi_post_log_out ~ group_fct, data = bdi_scale_out)

#summary(model_BDI_out_log)

library (sjPlot)
tab_model(model_BDI_out_log)

	bdi post log out
Predictors	Estimates	CI	p
(Intercept)	0.97	0.89 – 1.06	<0.001
group_fct [Group 2]	0.05	-0.07 – 0.17	0.440
Observations	96
R² / R² adjusted	0.006 / -0.004

Write-Up: (Main effect of Condition - log-transformed / outliers removed) Simple regression of post BDI score on treatment

Removed outliers and Compared Group 1 and Group 2 on Post log-transformed BDI composite score using a regression model. With raw BDI-post, there was not a significant difference between groups (t(94) = .78, p = .44). Specifically, Group 1 was not different (N=50, mean =0.97 , sd =0.3 , se = .04) from Group 2 (N=46, mean = 1.02, sd =0.3 , se =.04).Condition did not significantly predict depression outcomes.

MODELS FOR MAIN EFFECT OF TIME

Running the regression model for regressing baseline BDI scores on post BDI scores (Outliers included)

model_maineffect <- lm(bdi_post ~ bdi_pre, data = df_eq_bdi) 
#Testing for main effect. 

#summary(model_maineffect)

library (sjPlot)
tab_model(model_maineffect)

	bdi post
Predictors	Estimates	CI	p
(Intercept)	4.91	2.19 – 7.62	0.001
bdi_pre	0.53	0.36 – 0.70	<0.001
Observations	101
R² / R² adjusted	0.277 / 0.269

Write-Up: Main Effect of time on raw BDI composite scores (outliers included)

Running a regressing of baseline BDI scores on post BDI scores revealed a significant main effect of time (t(99) = 6.15, p =.00). Specifically, across both intervention and control group, BDI scores at pre (N=101, mean=13.43, sd=8.67, se=.86) were significantly higher than BDI scores at post (N=101, mean=11.99, sd =8.7, se=.87).

Running the regression model for regressing log-transformed baseline BDI scores on log-transformed post BDI scores (Outliers included)

model_maineffect_log <- lm(bdi_post_log ~ bdi_pre_log, data = df_eq_bdi)  #Testing for main effect of time. 

#summary(model_maineffect_log)

library (sjPlot)
tab_model(model_maineffect_log)

	bdi post log
Predictors	Estimates	CI	p
(Intercept)	0.48	0.27 – 0.69	<0.001
bdi_pre_log	0.50	0.31 – 0.69	<0.001
Observations	101
R² / R² adjusted	0.214 / 0.207

Write-Up: Main Effect of time on log-transformed BDI composite scores (outliers included)

Running a regressing of baseline BDI scores on post BDI scores revealed a significant main effect of time (t(99) = 5.2, p =.00). Specifically, across both intervention and control group, BDI scores at pre () were significantly higher than BDI scores at post ().

Running the regression model for regressing baseline BDI scores on post BDI scores (Outliers removed)

model_maineffect_out <- lm(bdi_post ~ bdi_pre, data = bdi_scale_out)  #Testing for main effect of time. 

#summary(model_maineffect_out)

library (sjPlot)
tab_model(model_maineffect_out)

	bdi post
Predictors	Estimates	CI	p
(Intercept)	5.33	2.92 – 7.75	<0.001
bdi_pre	0.46	0.29 – 0.63	<0.001
Observations	96
R² / R² adjusted	0.242 / 0.234

Write-Up: Main Effect of time on log-transformed BDI composite scores (outliers removed)

Removing outliers and running a regressing of baseline BDI scores on post BDI scores revealed a significant main effect of time (t(94) = 5.48, p =.00). Specifically, across both intervention and control groups, BDI scores at pre () were significantly higher than BDI scores at post ().

MIXED MODELS FOR TIME x CONDITION

Mixed Model controlling for baseline BDI (raw scores, outliers included)

model_BDI_pretopost <- lm(bdi_post ~ group_fct + bdi_pre, data = df_eq_bdi) #Regression for Post BDI across conditions - accounting for Pre BDI. group 1 will be the reference group. 

#summary(model_BDI_pretopost)

#str(model_BDI_pretopost)

#broom::tidy(model_BDI_pretopost)

tab_model(model_BDI_pretopost)

	bdi post
Predictors	Estimates	CI	p
(Intercept)	4.47	1.46 – 7.47	0.004
group_fct [Group 2]	1.02	-1.93 – 3.97	0.495
bdi_pre	0.52	0.35 – 0.69	<0.001
Observations	101
R² / R² adjusted	0.280 / 0.265

# Primary Index model results
#model_BDI_pretopost_results <- apa_print(model_BDI_pretopost)

#model_BDI_pretopost_results$table %>% apa_table(caption = "Model Regressing Follow-Up BDI Score on Condition - Controlling for Baseline BDI", note = "* p < 0.05; ** p < 0.01; *** p < 0.001", escape = TRUE)

#Write-Up: Mixed Model controlling for baseline BDI (raw scores, outliers included) Running a linear regression model, regressing follow-up BDI scores onto condition while holding baseline BDI constant did not reveal a significant difference between the two groups for depression outcomes (t(98)= .69, p = .50). Participants in group 1 did not have significantly different follow-up depression scores than group 2. Condition did not significantly predict follow-up depression when basleine depression was held constant.

Mixed Model controlling for baseline BDI with outliers removed

model_BDI_pretopost_out <- lm(bdi_post ~ group_fct + bdi_pre, data = bdi_scale_out) #Regression for Post BDI across conditions - accounting for Pre BDI. group 1 will be the reference group. 

#summary(model_BDI_pretopost_out)

tab_model(model_BDI_pretopost_out)

	bdi post
Predictors	Estimates	CI	p
(Intercept)	4.87	2.20 – 7.54	<0.001
group_fct [Group 2]	1.04	-1.46 – 3.55	0.410
bdi_pre	0.46	0.29 – 0.62	<0.001
Observations	96
R² / R² adjusted	0.247 / 0.231

Write-Up: Mixed Model controlling for baseline BDI with outliers removed

Running a linear regression model with outliers removed, regressing follow-up BDI scores onto condition while holding baseline BDI constant did not reveal a significant difference between the two groups for depression outcomes (t(93)= .83, p = .41). Participants in group 1 did not have significantly different follow-up depression scores than group 2. Condition did not significantly predict follow-up depression when basleine depression was held constant.

INTERACTION MODELS - LONG-FORMAT

Re-shaping the data to long-format

# Re-shaping data to long format

bdi_out <- dplyr::select(bdi_scale_out, PID, Condition_Group, group_fct, bdi_pre, bdi_post)
  
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

bdi_log_out <- dplyr::select(bdi_scale_out, PID, Condition_Group, group_fct, bdi_pre_log, bdi_post_log) #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())

Repeated measures ANOVA w/ outliers removed (long-format approach) - OPTION 1

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

	BDI score
Predictors	Estimates	CI	p
(Intercept)	10.42	8.38 – 12.46	<0.001
group_fct [Group 2]	1.32	-1.63 – 4.27	0.379
time [bdi_pre]	1.74	-1.15 – 4.63	0.236
group_fct [Group 2] * time [bdi_pre]	-0.72	-4.89 – 3.45	0.734
Observations	192
R² / R² adjusted	0.014 / -0.002

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)

Mixed Model using the long form of data.

library(lme4)
lmeModel = lmer(BDI_score ~ group_fct*time + (1|PID), data=bdi_out_long)
anova(lmeModel)

## Analysis of Variance Table
##                Df Sum Sq Mean Sq F value
## group_fct       1 15.138  15.138  0.5536
## time            1 93.521  93.521  3.4202
## group_fct:time  1  6.180   6.180  0.2260

INTERACTION MODELS - WIDE FORMAT

Running an interaction model - raw scores, outliers included

model_BDI_pretopostint <- lm(bdi_post ~ group_fct + bdi_pre + group_fct:bdi_pre, data = df_eq_bdi) #Interaction model 

#summary(model_BDI_pretopostint)

#str(model_BDI_pretopostint)

#broom::tidy(model_BDI_pretopostint)

tab_model(model_BDI_pretopostint)

	bdi post
Predictors	Estimates	CI	p
(Intercept)	3.30	-0.57 – 7.16	0.094
group_fct [Group 2]	3.24	-2.22 – 8.71	0.242
bdi_pre	0.61	0.36 – 0.87	<0.001
group_fct [Group 2] * bdi_pre	-0.17	-0.51 – 0.18	0.340
Observations	101
R² / R² adjusted	0.287 / 0.265

# Primary Index model results
#model_BDI_pretopostint_results <- apa_print(model_BDI_pretopostint)

#apa_print(model_BDI_pretopostint_results)$fullresults$group_fct:bdi_pre
  
#model_BDI_pretopostint_results$table %>% 
  #apa_table(caption = "Model Regressing Follow-Up BDI Score on Condition - Controlling for interaction between Baseline BDI and Condition",
           # note = "* p < 0.05; ** p < 0.01; *** p < 0.001",
           # escape = TRUE)

#Write-up: A mixed effect linear model analysis A mixed effect linear model analysis was performed to assess the effect of study condition on change in BDI scores from baseline to follow-up, which showed no significant difference between the two groups for depression outcomes (t(97)=-.96, p=.34). Participants in group 1 did not have significantly different follow-up depression scores than group 2. Condition did not significantly predict follow-up depression when baseline depression was held constant.

Interaction model w/ outliers removed (wide format version)

#create descriptives with outliers removed 

model_BDI_pretopostint_out <- lm(bdi_post ~ group_fct + bdi_pre + group_fct:bdi_pre, data = bdi_scale_out) #Repeated measures ANOVA. #This asks: Does the change in depression over time depend on what group you are in?

#coefficients:
#Int > expected value for BDI post in group 1, when pre BDI. 
#group_fct [Group 2] > change in BDI score as you move from group 1 to group 2.
#bdi_pre > Expected change in bdi post when all other predictors are a zero (would need to center variables for this to meaningful - look at the interactions lab)


#summary(model_BDI_pretopostint_out)

tab_model(model_BDI_pretopostint_out)

	bdi post
Predictors	Estimates	CI	p
(Intercept)	3.86	0.50 – 7.22	0.025
group_fct [Group 2]	3.11	-1.75 – 7.97	0.207
bdi_pre	0.54	0.30 – 0.78	<0.001
group_fct [Group 2] * bdi_pre	-0.17	-0.50 – 0.17	0.328
Observations	96
R² / R² adjusted	0.255 / 0.231

#

# ***

#FUTURE TO DO - make a long-format df and run the model like this: 
#lm(bdi_value ~ group_fct * pre_post)

#Graphs don't make sense in wide format
#library(sjPlot)

#plot_model(model = model_BDI_pretopostint_out, type = "int")

#plot this at +- 1 SD around mean . look at interaction lab to change this. 

#once it is in long format, plot relationship between time and depression as a function of group. This will now be two categorical predictors (group and time)
#plot_model(model = model_int, # model fit object
          # type = "int", # interaction 
          # mdrt.values = "meansd") # which values of the moderator variable to use
#or get code from factorial anova lab
# Look at the lab to get the simple slopes. rs partner 
#NEED TO FIX THIS

OTHER MODELS

Running a t-test on the change from pre to post. The target variable is the change score (sometimes called the gain or difference score)

# Models with the pre/post design taken into account

df_eq_bdi$change_score <- df_eq_bdi$bdi_post - df_eq_bdi$bdi_pre

ttestChange = t.test(change_score ~ group_fct, df_eq_bdi, var.equal=T)
ttestChange

## 
##  Two Sample t-test
## 
## data:  change_score by group_fct
## t = -0.2659, df = 99, p-value = 0.7909
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.806026  2.906497
## sample estimates:
## mean in group Group 1 mean in group Group 2 
##             -1.653846             -1.204082

t_sqrd <- ttestChange$statistic^2
t_sqrd

##          t 
## 0.07070288

Running and ANCOVA

#ancova <- aov(bdi_post~ bdi_pre + group_fct, df_eq_bdi)
#ancovalm <- lm(bdi_post~ bdi_pre + group_fct, df_eq_bdi)
#rbind(coef(ancova), coef(ancovalm)

#FUTURE TO DO: Filtering out participants who do not meet the clinical cut off for mild (14+) to severe depression.
#Descriptives - Pre BDI
#psych::describeBy(df_eq_bdi$bdi_pre, df_eq_bdi$group_fct, mat = TRUE)

#df_eq_bdi_dep <- df_eq_bdi %>% 
  #filter(bdi_pre > 13, )

#Regression for Post BDI (w/ mild-severe pre BDI only) across conditions

#model_BDI_dep <- lm(bdi_post ~ group_fct, data = df_eq_bdi_dep) 
#So, group 1 will be the reference group. 

#summary(model_BDI_dep)

#Removing outliers 
#df_eq_bdi_dep_outrmv <- df_eq_bdi_outrmv %>% 
  #filter(bdi_pre > 13, )

#model_BDI_dep_outrmv <- lm(bdi_post ~ group_fct, data = df_eq_bdi_dep_outrmv) 
#So, group 1 will be the reference group. 

#summary(model_BDI_dep_outrmv)

Equanimity_BDI Depression

Megan Lipsett

2/25/2020