Medicating conflict monitoring? Treatment with serotonergic antidepressants predicts stop-N200 but not P300 amplitudes in major depressive disorder

1 Preparation

1.1 Packages

library(tidyverse) # to data wrangle

## Warning: package 'tidyverse' was built under R version 4.5.2

## Warning: package 'readr' was built under R version 4.5.2

library(table1) # to create Table 1

## Warning: package 'table1' was built under R version 4.5.2

library(psych) # to assist in calculating and exporting Spearman's correlations
library(corrplot) # to create figure of correlation matrix
library(rstatix) # to simplify outlier detection, normality assumptions, and ANOVA analyses
library(ggplot2) # to create boxplots, QQ plots, etc.
library(patchwork) # to combine plots

## Warning: package 'patchwork' was built under R version 4.5.2

library(emmeans) # to calculate post-hoc tests
library(MatchIt) # for propensity matching calculating
library(optmatch) # for propensity matching cont. 
library(cobalt) # for Love plots for the propensity matching
library(lmtest) # for propensity matching ATE

## Warning: package 'lmtest' was built under R version 4.5.2

library(sandwich) # for propensity matching ATE

## Warning: package 'sandwich' was built under R version 4.5.2

library(broom) # for propensity matching ATE (tidy)

1.2 Working directory, residualizing incorrect event-related potentials (ERPs), and reading dataframes to analyze

setwd("C:\\Users\\WVaught\\OneDrive - Laureate Institute for Brain Research\\Documents\\01. EEG CORE\\05. Papers\\SSRI EEG\\01. Submission\\data")

dfT<-read.csv("T1000_redcap_wide-S.T0-2022-01-24.csv") # T1000 demographics and self-report
dfRT <- read.csv("df2visualize_all_residN2P3.csv") # T500 RTs
dfRT <- dfRT %>% select(c(id, ss_meanrt_go, ss_meanrt_goerror))
dfEEG<-read.csv("SST_P3_N2_ERN_091923.csv") # T500 SST ERPs 
dfMerged <- merge(dfT,dfEEG,by = 'id') # merging SST ERPs and T1000 demographics
# filter a dataframe for a total number of correct stop trials of > 20; https://doi.org/10.1016/S0167-8760(96)00743-X
dfMerged <- filter(dfMerged, AllCorrectSeg > 20)
dfMerged <- filter(dfMerged, Dx !="NA") # filtering based on available EEG data and exclusion criteria - a current substance dependence besides AUD
dfMerged <- merge(dfMerged, dfRT, by = 'id') # merging SST ERPs/demographics with SST RTs
dfMerged <- select(dfMerged, c(id, Age, Dx, Rx, RxMofA, Education_4level, Ethnicity, Sex, Income, PROMIS_AnxietyTscore, PROMIS_AlcoUseTscore, PROMIS_DepressTscore, RACE_AfricanAmerican, RACE_AlaskanNative, RACE_Asian, RACE_NativeAmerican, RACE_Other, RACE_PacificIslander, RACE_White, AllCorrectSeg, AllIncorrectSeg, GoCorrectSeg,GoIncorrectSeg, GoTargetRespSeg, GoTargetIncorrRespSeg, P3_FCz_AllCorrectAvg, P3_FCz_AllIncorrectAvg, P3_Cz_AllCorrectAvg, P3_Cz_AllIncorrectAvg, P3_Pz_AllCorrectAvg, P3_Pz_AllIncorrectAvg, N2_FCz_AllCorrectAvg, N2_FCz_AllIncorrectAvg, ERN_FCz_GoTargetRespAvg, ERN_FCz_GoTargetIncorrectAvg, ERN_FCz_AllIncorrectRespAvg, ss_meanrt_go, ss_meanrt_goerror))
dfMerged$Ethnicity <- str_replace_all(dfMerged$Ethnicity, c("Not Hispanic or Latino" = "Not Hispanic or Latine", "Hispanic or Latino" = "Hispanic or Latine"))
# filter a dataframe for only MDD/COM participants and selecting for SERT- mechanistic medications
dfMerged <- dfMerged %>%
  filter(RxMofA == "SSRI" | RxMofA =="SNRI" | RxMofA == "None")
dfMergedMDD <- dfMerged %>%
  filter(Dx !="HC") %>%
  filter(RxMofA == "SSRI" | RxMofA =="SNRI" | RxMofA =="None")
dfMerged$N2_FCz_AllCorrectAvg <- as.numeric(dfMerged$N2_FCz_AllCorrectAvg)
dfMerged$N2_FCz_AllIncorrectAvg <- as.numeric(dfMerged$N2_FCz_AllIncorrectAvg)
dfMergedMDD$N2_FCz_AllCorrectAvg <- as.numeric(dfMergedMDD$N2_FCz_AllCorrectAvg)
dfMergedMDD$N2_FCz_AllIncorrectAvg <- as.numeric(dfMergedMDD$N2_FCz_AllIncorrectAvg)

# regressing out the correct signal from the incorrect signal according to procedure introduced by Meyer et al., 2017 - Considering ERP difference scores as individual difference measures: Issues with subtraction and alternative approaches, https://doi.org/10.1111/psyp.12664
dfMerged$N2_FCz_AllCorrectAvg <- as.numeric(dfMerged$N2_FCz_AllCorrectAvg)
dfMerged$N2_FCz_AllIncorrectAvg <- as.numeric(dfMerged$N2_FCz_AllIncorrectAvg)
# residualized incorrect signal in N200 @ FCz
residual_incorrectFCz_N2 <- lm(N2_FCz_AllIncorrectAvg ~ N2_FCz_AllCorrectAvg, data = dfMerged)
dfMerged$N2resid <- resid(residual_incorrectFCz_N2)
# residualized incorrect signal in P300 @ Cz
residual_incorrectCz_P3 <- lm(P3_Cz_AllIncorrectAvg ~ P3_Cz_AllCorrectAvg, data = dfMerged)
dfMerged$P3resid_Cz <- resid(residual_incorrectCz_P3)
# residualized incorrect signal in P300 @ Pz
residual_incorrectPz_P3 <- lm(P3_Pz_AllIncorrectAvg ~ P3_Pz_AllCorrectAvg, data = dfMerged)
dfMerged$P3resid_Pz <- resid(residual_incorrectPz_P3)
# complete case dfMerged to only get participants who have ERN data avaliable
dfMergedERN <- dfMerged[complete.cases(dfMerged),]
# filter df to only include participants with > 6 error trials as advised by Olvet and Hajcak, 2009; https://doi.org/10.1111/j.1469-8986.2009.00848.x
dfMergedERN <- filter(dfMergedERN, GoTargetIncorrRespSeg > 5)
# residualized incorrect signal in N200 @ FCz - repeat for ps who have ERN
residual_incorrectFCz_N2 <- lm(N2_FCz_AllIncorrectAvg ~ N2_FCz_AllCorrectAvg, data = dfMergedERN)
dfMergedERN$N2resid <- resid(residual_incorrectFCz_N2)
# residualized incorrect signal in P300 @ Cz - repeat for ps who have ERN
residual_incorrectCz_P3 <- lm(P3_Cz_AllIncorrectAvg ~ P3_Cz_AllCorrectAvg, data = dfMergedERN)
dfMergedERN$P3resid_Cz <- resid(residual_incorrectCz_P3)
# residualized incorrect signal in P300 @ Pz - repeat for ps who have ERN
residual_incorrectPz_P3 <- lm(P3_Pz_AllIncorrectAvg ~ P3_Pz_AllCorrectAvg, data = dfMergedERN)
dfMergedERN$P3resid_Pz <- resid(residual_incorrectPz_P3)
# residualized incorrect signal in ERN @ FCz 
residual_incorrectFCz_ERN <- lm(ERN_FCz_GoTargetIncorrectAvg ~ ERN_FCz_GoTargetRespAvg, data = dfMergedERN)
dfMergedERN$ERNresid <- resid(residual_incorrectFCz_ERN)

2 Tests of Normality

2.1 Shapiro-Wilk tests, describe function, histograms, and QQ plots in all participants

shapiro.test(dfMerged$Age)

## 
##  Shapiro-Wilk normality test
## 
## data:  dfMerged$Age
## W = 0.92957, p-value = 4.141e-09

describe(dfMerged$Age)

##    vars   n  mean    sd median trimmed   mad   min   max range skew kurtosis
## X1    1 233 34.96 11.37  34.05    34.5 15.31 18.33 55.61 37.28 0.27    -1.26
##      se
## X1 0.74

hist(dfMerged$Age)

qqnorm(dfMerged$Age)
qqline(dfMerged$Age)

shapiro.test(dfMerged$Income)

## 
##  Shapiro-Wilk normality test
## 
## data:  dfMerged$Income
## W = 0.50145, p-value < 2.2e-16

describe(dfMerged$Income)

##    vars   n     mean       sd median  trimmed   mad min    max  range skew
## X1    1 196 56399.19 81349.01  37500 44456.96 37065   0 750000 750000 5.88
##    kurtosis      se
## X1    44.48 5810.64

hist(dfMerged$Income)

qqnorm(dfMerged$Income)
qqline(dfMerged$Income)

shapiro.test(dfMerged$N2_FCz_AllCorrectAvg)

## 
##  Shapiro-Wilk normality test
## 
## data:  dfMerged$N2_FCz_AllCorrectAvg
## W = 0.95033, p-value = 3.65e-07

describe(dfMerged$N2_FCz_AllCorrectAvg)

##    vars   n  mean   sd median trimmed  mad    min  max range  skew kurtosis
## X1    1 233 -4.47 4.96  -3.62   -4.03 4.34 -20.18 5.24 25.42 -0.88     0.78
##      se
## X1 0.33

hist(dfMerged$N2_FCz_AllCorrectAvg)

qqnorm(dfMerged$N2_FCz_AllCorrectAvg)
qqline(dfMerged$N2_FCz_AllCorrectAvg)

shapiro.test(dfMerged$P3_Pz_AllCorrectAvg)

## 
##  Shapiro-Wilk normality test
## 
## data:  dfMerged$P3_Pz_AllCorrectAvg
## W = 0.99036, p-value = 0.1245

describe(dfMerged$P3_Pz_AllCorrectAvg)

##    vars   n mean   sd median trimmed  mad   min   max range skew kurtosis  se
## X1    1 233 5.19 4.54   5.24    5.07 4.64 -5.22 19.62 24.84 0.28    -0.02 0.3

hist(dfMerged$P3_Pz_AllCorrectAvg)

qqnorm(dfMerged$P3_Pz_AllCorrectAvg)
qqline(dfMerged$P3_Pz_AllCorrectAvg)

shapiro.test(dfMerged$ERN_FCz_GoTargetRespAvg)

## 
##  Shapiro-Wilk normality test
## 
## data:  dfMerged$ERN_FCz_GoTargetRespAvg
## W = 0.85287, p-value = 3.987e-14

describe(dfMerged$ERN_FCz_GoTargetRespAvg)

##    vars   n  mean   sd median trimmed mad   min   max range skew kurtosis   se
## X1    1 233 -1.45 3.49  -1.27   -1.38 2.3 -15.9 21.74 37.64 0.85    10.85 0.23

hist(dfMerged$ERN_FCz_GoTargetRespAvg)

qqnorm(dfMerged$ERN_FCz_GoTargetRespAvg)
qqline(dfMerged$ERN_FCz_GoTargetRespAvg)

shapiro.test(dfMerged$PROMIS_DepressTscore)

## 
##  Shapiro-Wilk normality test
## 
## data:  dfMerged$PROMIS_DepressTscore
## W = 0.96895, p-value = 5.575e-05

describe(dfMerged$PROMIS_DepressTscore)

##    vars   n  mean    sd median trimmed   mad  min  max range  skew kurtosis
## X1    1 233 57.91 10.17   58.9   58.48 10.23 34.2 79.8  45.6 -0.48    -0.25
##      se
## X1 0.67

hist(dfMerged$PROMIS_DepressTscore)

qqnorm(dfMerged$PROMIS_DepressTscore)
qqline(dfMerged$PROMIS_DepressTscore)

shapiro.test(dfMerged$PROMIS_AnxietyTscore)

## 
##  Shapiro-Wilk normality test
## 
## data:  dfMerged$PROMIS_AnxietyTscore
## W = 0.95847, p-value = 2.835e-06

describe(dfMerged$PROMIS_AnxietyTscore)

##    vars   n  mean   sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 233 59.24 9.71   60.8   60.04 9.49 32.9 80.5  47.6 -0.7      0.1 0.64

hist(dfMerged$PROMIS_AnxietyTscore)

qqnorm(dfMerged$PROMIS_AnxietyTscore)
qqline(dfMerged$PROMIS_AnxietyTscore)

shapiro.test(dfMerged$PROMIS_AlcoUseTscore)

## 
##  Shapiro-Wilk normality test
## 
## data:  dfMerged$PROMIS_AlcoUseTscore
## W = 0.91755, p-value = 5.046e-10

describe(dfMerged$PROMIS_AlcoUseTscore)

##    vars   n mean   sd median trimmed mad  min  max range  skew kurtosis   se
## X1    1 231 47.4 6.95     50   47.78 4.6 33.4 62.8  29.4 -0.51    -0.19 0.46

hist(dfMerged$PROMIS_AlcoUseTscore)

qqnorm(dfMerged$PROMIS_AlcoUseTscore)
qqline(dfMerged$PROMIS_AlcoUseTscore)

3 Bivariate Correlations

3.1 Tibbles for Correlation Matrices

# create tibble for all ps
vars4matrixAll <- tibble(
  Age = dfMerged$Age,
  "Corr N200 FCz" = dfMerged$N2_FCz_AllCorrectAvg,
  "Corr P300 Pz" = dfMerged$P3_Pz_AllCorrectAvg,
  "CRN FCz" = dfMerged$ERN_FCz_GoTargetRespAvg,
  "Log Go RTs" = log(dfMerged$ss_meanrt_go),
  "PROMIS Alcohol Use" = dfMerged$PROMIS_AlcoUseTscore,
  "PROMIS Anxiety" = dfMerged$PROMIS_AnxietyTscore,
  "PROMIS Depression" = dfMerged$PROMIS_DepressTscore
)

vars4matrixMDD <- tibble(
  Age = dfMergedMDD$Age,
  "Corr N200 FCz" = dfMergedMDD$N2_FCz_AllCorrectAvg,
  "Corr P300 Cz" = dfMergedMDD$P3_Cz_AllCorrectAvg,
  "Corr P300 Pz" = dfMergedMDD$P3_Pz_AllCorrectAvg,
  "CRN FCz" = dfMergedMDD$ERN_FCz_GoTargetRespAvg,
  "Log Go RTs" = log(dfMergedMDD$ss_meanrt_go),
  "PROMIS Alcohol Use" = dfMergedMDD$PROMIS_AlcoUseTscore,
  "PROMIS Anxiety" = dfMergedMDD$PROMIS_AnxietyTscore,
  "PROMIS Depression" = dfMergedMDD$PROMIS_DepressTscore
)

3.2 Correlation Matrix for All Participants & N200, P300, CRN ERPs (n = 241)

# correlate the measures in the tibble for all participants
rAll <- corr.test(vars4matrixAll, method = "spearman", adjust = "none")$r
pAll <- corr.test(vars4matrixAll, method = "spearman", adjust = "none")$p
# create correlation matrix of the correlations for all participants
corrplot(rAll, p.mat = pAll, method = 'ellipse',diag = FALSE, col = COL2('PuOr'), tl.col = "black", tl.cex = 1.0, cl.cex = 1.00, type = "lower", sig.level = c(0.001,0.01,0.05), pch.cex = 1.00, insig = "label_sig", pch.col = "black")

3.3 Scatterplots for Significant Associations for All Ps & Some ERPs

# linear model of PROMIS_Anx ~ corrP3
model <- lm(PROMIS_AnxietyTscore ~ P3_Pz_AllCorrectAvg, data = dfMerged)
eqn <- sprintf(
  "italic(y) == %.3g %.3g * italic(x) * ',' ~~ italic(r)^2 ~ '=' ~ %.2g",
  coef(model)[1],
  coef(model)[2],
  summary(model)$r.squared
)
parse(text=eqn)

## expression(italic(y) == 61.1 - 0.366 * italic(x) * "," ~ ~italic(r)^2 ~ 
##     "=" ~ 0.029)

# scatterplot with lm fit
scatterAnx_corrP3_Pz <-ggplot(dfMerged, aes(x = P3_Pz_AllCorrectAvg, y = PROMIS_AnxietyTscore)) + 
  geom_point() + 
  stat_smooth(method = lm, se = TRUE,colour = "black") + 
  annotate(
    "text",
    x = Inf, y = -Inf,
    label = eqn, parse = TRUE,
    hjust = 1.1, vjust = -1.5
  ) + theme_classic()
scatterAnx_corrP3_Pz <- scatterAnx_corrP3_Pz + ggtitle("PROMIS Anxiety and Correct-Stop P300 Mean Amplitudes @ Pz") + xlab("Correct-Stop P300 Mean Amplitudes @ Pz (μV)") + ylab("PROMIS Anxiety (T score)")

3.4 Correlation Matrix for MDD/COM Participants & N200, P300, CRN ERPs (n = 185)

rMDD <- corr.test(vars4matrixMDD, method = "spearman",adjust = "none")$r
pMDD <- corr.test(vars4matrixMDD, method = "spearman",adjust = "none")$p
corrplot(rMDD, p.mat = pMDD, method = 'ellipse', diag = FALSE, type = 'lower', sig.level = c(0.001,0.01,0.05), pch.cex = 0.75, insig = 'label_sig', pch.col = "black")

4 Table 1

4.1 Tibbles for Table 1

# mutate to create groups for Table 1
dfMerged <- dfMerged %>%
  mutate(Group = case_when(
    Dx == "HC" & Sex == "Female" ~ "HC F",
    Dx == "HC" & Sex == "Male" ~ "HC M",
    Dx == "MDD" & Sex == "Female" & Rx == "NonSero" ~ "Unmed MDD F",
    Dx == "MDD" & Sex == "Male" & Rx == "NonSero" ~ "Unmed MDD M",
    Dx == "COM" & Sex == "Female" & Rx == "NonSero" ~ "Unmed COM F",
    Dx == "COM" & Sex == "Male" & Rx == "NonSero" ~ "Unmed COM M",
    Dx == "MDD" & Sex == "Female" & Rx == "Sero" ~ "Med MDD F",
    Dx == "MDD" & Sex == "Male" & Rx == "Sero" ~ "Med MDD M",
    Dx == "COM" & Sex == "Female" & Rx == "Sero" ~ "Med COM F",
    Dx == "COM" & Sex == "Male" & Rx == "Sero" ~ "Med COM M"
))
# tibble time for Table 1
vars4table1 <- tibble(
  Age = dfMerged$Age,
  Education = dfMerged$Education_4level,
  Income = dfMerged$Income,
  "Alaskan Native" = dfMerged$RACE_AlaskanNative,
  "Asian" = dfMerged$RACE_Asian,
  "Black" = dfMerged$RACE_AfricanAmerican,
  "Hispanic or Latine" = dfMerged$Ethnicity,
  "Native American" = dfMerged$RACE_NativeAmerican,
  "Other" = dfMerged$RACE_Other,
  "Pacific Islander" = dfMerged$RACE_PacificIslander,
  "White" = dfMerged$RACE_White,
  "Corr N200 FCz" = dfMerged$N2_FCz_AllCorrectAvg,
  "Corr P300 Pz" = dfMerged$P3_Pz_AllCorrectAvg,
  "CRN FCz" = dfMerged$ERN_FCz_GoTargetRespAvg,
  "Go RTs" = dfMerged$ss_meanrt_go,
  "PROMIS Alcohol Use" = dfMerged$PROMIS_AlcoUseTscore,
  "PROMIS Anxiety" = dfMerged$PROMIS_AnxietyTscore,
  "PROMIS Depression" = dfMerged$PROMIS_DepressTscore,
  Group = dfMerged$Group
)

4.2 Table 1

# Age
describeAge <- describeBy(vars4table1$Age, group = vars4table1$Group)
aovAge <- aov(Age ~ Group, data = vars4table1)
summary(aovAge)

##              Df Sum Sq Mean Sq F value Pr(>F)
## Group         9   1187   131.9   1.021  0.424
## Residuals   223  28800   129.2

# Income
describeIncome <- describeBy(vars4table1$Income, group = vars4table1$Group)
kwIncome <- kruskal.test(Income ~ Group, data = vars4table1)
kwIncome

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Income by Group
## Kruskal-Wallis chi-squared = 8.5917, df = 9, p-value = 0.4758

# PROMIS Anxiety
describePROMISanx <- describeBy(vars4table1$`PROMIS Anxiety`, group = vars4table1$Group)
aovPROMISanx <- aov(`PROMIS Anxiety` ~ Group, data = vars4table1)
summary(aovPROMISanx)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## Group         9  11646  1294.0   28.21 <2e-16 ***
## Residuals   223  10231    45.9                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# PROMIS Depression
describePROMISdep <- describeBy(vars4table1$`PROMIS Depression`, group = vars4table1$Group)
aovPROMISdep <- aov(`PROMIS Depression` ~ Group, data = vars4table1)
summary(aovPROMISdep)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## Group         9  12381  1375.7   26.42 <2e-16 ***
## Residuals   223  11611    52.1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# PROMIS Alcohol Use
describePROMISalcUse <- describeBy(vars4table1$`PROMIS Alcohol Use`, group = vars4table1$Group)
aovPROMISalcUse <- aov(`PROMIS Alcohol Use` ~ Group, data = vars4table1)
summary(aovPROMISalcUse)

##              Df Sum Sq Mean Sq F value  Pr(>F)   
## Group         9   1144  127.10    2.82 0.00371 **
## Residuals   221   9961   45.07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2 observations deleted due to missingness

# N200
describeN200 <- describeBy(vars4table1$`Corr N200 FCz`, group = vars4table1$Group)
aovN200 <- aov(`Corr N200 FCz` ~ Group, data = vars4table1)
summary(aovN200)

##              Df Sum Sq Mean Sq F value Pr(>F)
## Group         9    326   36.23   1.501  0.149
## Residuals   223   5384   24.14

# P300 
describeP300 <- describeBy(vars4table1$`Corr P300 Pz`, group = vars4table1$Group)
aovP300 <- aov(`Corr P300 Pz` ~ Group, data = vars4table1)
summary(aovP300)

##              Df Sum Sq Mean Sq F value   Pr(>F)    
## Group         9    572   63.54   3.362 0.000696 ***
## Residuals   223   4214   18.90                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# CRN
describeCRN <- describeBy(vars4table1$`CRN FCz`, group = vars4table1$Group)
aovCRN <- aov(`CRN FCz` ~ Group, data = vars4table1)
summary(aovCRN)

##              Df Sum Sq Mean Sq F value Pr(>F)
## Group         9  152.1   16.90   1.407  0.186
## Residuals   223 2679.1   12.01

# Go RTs
describeRTs <- describeBy(vars4table1$`Go RTs`, group = vars4table1$Group)
aovRTs <- aov(`Go RTs` ~ Group, data = vars4table1)
summary(aovRTs)

##              Df Sum Sq Mean Sq F value Pr(>F)  
## Group         9  0.238 0.02644    1.67 0.0975 .
## Residuals   223  3.531 0.01584                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Alaskan Native
tableAN <- table(vars4table1$`Alaskan Native`, vars4table1$Group) # 0 participants
# Asian
tableAs <- table(vars4table1$Asian, vars4table1$Group)
chisq.test(tableAs)

## Warning in chisq.test(tableAs): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tableAs
## X-squared = 8.332, df = 9, p-value = 0.5011

# Black
tableBl <- table(vars4table1$Black, vars4table1$Group)
chisq.test(tableBl)

## Warning in chisq.test(tableBl): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tableBl
## X-squared = 7.5109, df = 9, p-value = 0.5841

# Hispanic
tableHL <- table(vars4table1$`Hispanic or Latine`, vars4table1$Group)
chisq.test(tableHL)

## Warning in chisq.test(tableHL): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tableHL
## X-squared = 2.9378, df = 9, p-value = 0.9667

# Native American
tableNA <- table(vars4table1$`Native American`, vars4table1$Group)
chisq.test(tableNA)

## Warning in chisq.test(tableNA): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tableNA
## X-squared = 10.541, df = 9, p-value = 0.3085

# Other 
tableOt <- table(vars4table1$Other, vars4table1$Group)
chisq.test(tableOt)

## Warning in chisq.test(tableOt): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tableOt
## X-squared = 7.9402, df = 9, p-value = 0.5402

# Pacific Islander
tablePI <- table(vars4table1$`Pacific Islander`, vars4table1$Group)
chisq.test(tablePI)

## Warning in chisq.test(tablePI): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tablePI
## X-squared = 4.314, df = 9, p-value = 0.8896

# White
tableWh <- table(vars4table1$White, vars4table1$Group)
chisq.test(tableWh)

## Warning in chisq.test(tableWh): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tableWh
## X-squared = 9.1886, df = 9, p-value = 0.4201

# Education 
tableEdu <- table(vars4table1$Education, vars4table1$Group)
chisq.test(tableEdu)

## Warning in chisq.test(tableEdu): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tableEdu
## X-squared = 23.571, df = 27, p-value = 0.654

# Missing data
perColumnNAs <- colSums(is.na(vars4table1))
perColumnNAs

##                Age          Education             Income     Alaskan Native 
##                  0                  2                 37                  0 
##              Asian              Black Hispanic or Latine    Native American 
##                  0                  0                  1                  0 
##              Other   Pacific Islander              White      Corr N200 FCz 
##                  0                  0                  0                  0 
##       Corr P300 Pz            CRN FCz             Go RTs PROMIS Alcohol Use 
##                  0                  0                  0                  2 
##     PROMIS Anxiety  PROMIS Depression              Group 
##                  0                  0                  0

5 Outliers Removal

5.1 9 Participants Removed from Correct Stop-N200 Amplitudes for All Participants

outliers <- dfMerged %>%
  group_by(Sex, Dx) %>%
  identify_outliers(N2_FCz_AllCorrectAvg)
corrN2_fcz_outliers_HC <- dfMerged %>%
  left_join(outliers, by = colnames(dfMerged))
corrN2_fcz_outliers_HC <- corrN2_fcz_outliers_HC %>%
  mutate(
    is.outlier = replace_na(is.outlier, FALSE),
    is.extreme = replace_na(is.extreme, FALSE)
  )
corrN2_FCz_outliersRemoved_HC <- corrN2_fcz_outliers_HC %>% filter(!is.outlier)

5.2 5 Participants Removed from Correct Stop-P300 Amplitudes for All Participants

outliers <- dfMerged %>%
  group_by(Sex, Dx) %>%
  identify_outliers(P3_Pz_AllCorrectAvg)
corrP3_pz_outliers_HC <- dfMerged %>%
  left_join(outliers, by = colnames(dfMerged))
corrP3_pz_outliers_HC <- corrP3_pz_outliers_HC %>%
  mutate(
    is.outlier = replace_na(is.outlier, FALSE),
    is.extreme = replace_na(is.extreme, FALSE)
  )
corrP3_Pz_outliersRemoved_HC <- corrP3_pz_outliers_HC %>% filter(!is.outlier)

5.3 13 Participants Removed from CRN Amplitudes for All Participants

outliers <- dfMerged %>%
  group_by(Sex, Dx) %>% 
  identify_outliers(ERN_FCz_GoTargetRespAvg)
CRN_FCz_outliers_HC <- dfMerged %>%
  left_join(outliers, by = colnames(dfMerged))
CRN_FCz_outliers_HC <- CRN_FCz_outliers_HC %>%
  mutate(
    is.outlier = replace_na(is.outlier, FALSE),
    is.extreme = replace_na(is.extreme, FALSE)
  )
CRN_FCz_outliersRemoved_HC <- CRN_FCz_outliers_HC %>% filter(!is.outlier)

5.4 7 Participants Removed from Correct Stop-N200 Amplitudes for MDD/COM Participants

outliers <- dfMergedMDD %>%
  group_by(Sex, Dx, Rx) %>%
  identify_outliers(N2_FCz_AllCorrectAvg)
corrN2_FCz_outliers_MDD <- dfMergedMDD %>%
  left_join(outliers, by = colnames(dfMergedMDD))
corrN2_FCz_outliers_MDD <- corrN2_FCz_outliers_MDD %>%
  mutate(
    is.outlier = replace_na(is.outlier, FALSE),
    is.extreme = replace_na(is.extreme, FALSE)
  )
corrN2_FCz_outliersRemoved_MDD <- corrN2_FCz_outliers_MDD %>% filter(!is.outlier)

5.5 2 Participants Removed from Correct Stop-P300 Amplitudes for MDD/COM Participants

outliers <- dfMergedMDD %>%
  group_by(Sex, Dx, Rx) %>%
  identify_outliers(P3_Pz_AllCorrectAvg)
corrP3_Pz_outliers_MDD <- dfMergedMDD %>%
  left_join(outliers, by = colnames(dfMergedMDD))
corrP3_Pz_outliers_MDD <- corrP3_Pz_outliers_MDD %>%
  mutate(
    is.outlier = replace_na(is.outlier, FALSE),
    is.extreme = replace_na(is.extreme, FALSE)
  )
corrP3_Pz_outliersRemoved_MDD <- corrP3_Pz_outliers_MDD %>% filter(!is.outlier)

5.6 10 Participants Removed from CRN Amplitudes for MDD/COM Participants

outliers <- dfMergedMDD %>%
  group_by(Sex, Dx, Rx) %>%
  identify_outliers(ERN_FCz_GoTargetRespAvg)
CRN_FCz_outliers_MDD <- dfMergedMDD %>%
  left_join(outliers, by = colnames(dfMergedMDD))
CRN_FCz_outliers_MDD <- CRN_FCz_outliers_MDD %>%
  mutate(
    is.outlier = replace_na(is.outlier, FALSE),
    is.extreme = replace_na(is.extreme, FALSE)
  )
CRN_FCz_outliersRemoved_MDD <- CRN_FCz_outliers_MDD %>% filter(!is.outlier)

6 Statistical analyses between HC and MDD/COM samples

6.1 Correct Stop-N200 Amplitudes for All Participants

qqPlot_N2_FCz_HC <- ggplot(corrN2_FCz_outliersRemoved_HC, aes(sample = N2_FCz_AllCorrectAvg)) + 
  geom_qq() +
  geom_qq_line()
qqPlot_N2_FCz_HC <- qqPlot_N2_FCz_HC + ylab("Sample") + xlab("Theoretical") + theme_classic()

# ANOVA
aovCorrN2FCz_HC <- corrN2_FCz_outliersRemoved_HC %>% 
  anova_test(N2_FCz_AllCorrectAvg ~ Dx * Sex, type = 3) %>%
  adjust_pvalue(p.col = "p", output.col = "p.adj", method = "bonferroni")
aovCorrN2FCz_HC # significant main effect of Sex

## ANOVA Table (type III tests)
## 
##   Effect DFn DFd     F     p p<.05   ges p.adj
## 1     Dx   2 218 0.119 0.888       0.001 1.000
## 2    Sex   1 218 6.119 0.014     * 0.027 0.042
## 3 Dx:Sex   2 218 3.832 0.023     * 0.034 0.069

# pairwise comparison (inspect main effects)
corrN2_FCz_outliersRemoved_HC %>%
  pairwise_t_test(
    N2_FCz_AllCorrectAvg ~ Sex,
    p.adjust.method = "bonferroni"
  )

## # A tibble: 1 × 9
##   .y.            group1 group2    n1    n2       p p.signif   p.adj p.adj.signif
## * <chr>          <chr>  <chr>  <int> <int>   <dbl> <chr>      <dbl> <chr>       
## 1 N2_FCz_AllCor… Female Male     148    76 0.00164 **       0.00164 **

vioPlotCorrN200FCz_HC <- ggplot(corrN2_FCz_outliersRemoved_HC, aes (x = Sex, y = N2_FCz_AllCorrectAvg, fill = Sex)) + geom_violin(trim = FALSE) + labs(x = "Sex", y = "Amplitude (μV)") + scale_fill_manual(values = c("#5f6266","#f8f9fa"), guide = "none")
vioPlotCorrN200FCz_HC <- vioPlotCorrN200FCz_HC + geom_jitter(shape = 16, position = position_jitter(0.2)) + theme_classic()

# creating dataframe with summary statistics 
corrN2_FCz_HC_Sum <- corrN2_FCz_outliersRemoved_HC %>%
  group_by(Sex) %>%
  summarize(mean=mean(N2_FCz_AllCorrectAvg),
            se=sd(N2_FCz_AllCorrectAvg)/sqrt(length(N2_FCz_AllCorrectAvg))
            )
# creating bar graph (with error bars)
corrN2_FCz_HC_barGraph <- ggplot(corrN2_FCz_HC_Sum, aes(x = Sex, y = mean, fill = Sex)) + 
  geom_col(colour = 'black') + 
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.1) + theme_classic() + scale_fill_manual(values = c("#5f6266","#f8f9fa"), guide = "none")

patchwork_corrN2_FCz_All <-
  qqPlot_N2_FCz_HC +
  corrN2_FCz_HC_barGraph +
  vioPlotCorrN200FCz_HC
patchwork_corrN2_FCz_All

6.2 Correct Stop-P300 Amplitudes for All Participants

qqPlot_corrP3_Pz_HC <- ggplot(corrP3_Pz_outliersRemoved_HC, aes(sample = P3_Pz_AllCorrectAvg)) + 
  geom_qq() +
  geom_qq_line()
qqPlot_corrP3_Pz_HC <- qqPlot_corrP3_Pz_HC + ylab("Sample") + xlab("Theoretical") + theme_classic()

# ANOVA
aovCorrP3Pz_HC <- corrP3_Pz_outliersRemoved_HC %>% 
  anova_test(P3_Pz_AllCorrectAvg ~ Dx * Sex, type = 3) %>%
  adjust_pvalue(p.col = "p", output.col = "p.adj", method = "bonferroni")
aovCorrP3Pz_HC # significant main effect of Dx

## ANOVA Table (type III tests)
## 
##   Effect DFn DFd     F        p p<.05   ges    p.adj
## 1     Dx   2 222 8.156 0.000382     * 0.068 0.001146
## 2    Sex   1 222 3.368 0.068000       0.015 0.204000
## 3 Dx:Sex   2 222 1.269 0.283000       0.011 0.849000

# pairwise comparison (inspect main effects)
corrP3_Pz_outliersRemoved_HC %>%
  pairwise_t_test(
    P3_Pz_AllCorrectAvg ~ Dx,
    p.adjust.method = "bonferroni"
  )

## # A tibble: 3 × 9
##   .y.            group1 group2    n1    n2       p p.signif   p.adj p.adj.signif
## * <chr>          <chr>  <chr>  <int> <int>   <dbl> <chr>      <dbl> <chr>       
## 1 P3_Pz_AllCorr… COM    HC       127    43 1.17e-4 ***      3.51e-4 ***         
## 2 P3_Pz_AllCorr… COM    MDD      127    58 7.97e-1 ns       1   e+0 ns          
## 3 P3_Pz_AllCorr… HC     MDD       43    58 1.4 e-3 **       4.2 e-3 **

vioPlotCorrP300Pz_HC <- ggplot(corrP3_Pz_outliersRemoved_HC, aes (x = Dx, y = P3_Pz_AllCorrectAvg, fill = Dx)) + geom_violin(trim = FALSE) + labs(x = "Diagnosis", y = "Amplitude (μV)")
vioPlotCorrP300Pz_HC <- vioPlotCorrP300Pz_HC + geom_jitter(shape = 16, position = position_jitter(0.2)) + theme_classic() + scale_fill_manual(values = c("#5f6266","#a5a9ae","#f8f9fa"), guide = "none")

# creating dataframe with summary statistics 
corrP3_Pz_HC_Sum <- corrP3_Pz_outliersRemoved_HC %>%
  group_by(Dx) %>%
  summarize(mean=mean(P3_Pz_AllCorrectAvg),
            se=sd(P3_Pz_AllCorrectAvg)/sqrt(length(P3_Pz_AllCorrectAvg))
            )
# creating bar graph (with error bars)
corrP3_Pz_HC_barGraph <- ggplot(corrP3_Pz_HC_Sum, aes(x = Dx, y = mean, fill = Dx)) + 
  geom_col(colour = 'black') + 
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.1) + theme_classic() + scale_fill_manual(values = c("#5f6266","#a5a9ae","#f8f9fa"), guide = "none")

patchwork_corrP3_Pz_All <- qqPlot_corrP3_Pz_HC +
  corrP3_Pz_HC_barGraph +
  vioPlotCorrP300Pz_HC
patchwork_corrP3_Pz_All

6.3 CRN Amplitudes for All Participants

qqPlot_CRN_FCz_HC <- ggplot(CRN_FCz_outliersRemoved_HC, aes(sample = ERN_FCz_GoTargetRespAvg)) + 
  geom_qq() +
  geom_qq_line()
qqPlot_CRN_FCz_HC + ylab("Sample") + xlab("Theoretical") + ggtitle("QQ Plot of CRN Mean Amplitudes @ FCz in All Participants") + theme_classic()

aovCRNFCz_HC <- CRN_FCz_outliersRemoved_HC %>%
  anova_test(ERN_FCz_GoTargetRespAvg ~ Dx * Sex, type = 3) %>%
  adjust_pvalue(p.col = "p", output.col = "p.adj", method = "bonferroni")
aovCRNFCz_HC # no significant effects

## ANOVA Table (type III tests)
## 
##   Effect DFn DFd     F     p p<.05   ges p.adj
## 1     Dx   2 214 0.213 0.809       0.002     1
## 2    Sex   1 214 0.409 0.523       0.002     1
## 3 Dx:Sex   2 214 0.162 0.850       0.002     1

7 Statistical analyses within the MDD/COM sample

7.1 Correct Stop-N200 Amplitudes for MDD/COM Participants

qqPlot_N2_FCz_MDD <- ggplot(corrN2_FCz_outliersRemoved_MDD, aes(sample = N2_FCz_AllCorrectAvg)) + 
  geom_qq() +
  geom_qq_line()
qqPlot_N2_FCz_MDD <- qqPlot_N2_FCz_MDD + ylab("Sample") + xlab("Theoretical") + theme_classic()

# using https://www.youtube.com/watch?v=7vwewIsSAls for learning this technique
corrN2_FCz_outliersRemoved_MDD$RxCoded <- str_replace_all(corrN2_FCz_outliersRemoved_MDD$Rx, c("NonSero"="0","Sero"="1"))
corrN2_FCz_outliersRemoved_MDD$RxCoded <- as.numeric(corrN2_FCz_outliersRemoved_MDD$RxCoded)
# remove participants who have incomplete PROMIS measures (or other covariates)
corrN2_FCz_outliersRemoved_MDD <- corrN2_FCz_outliersRemoved_MDD[complete.cases(corrN2_FCz_outliersRemoved_MDD[ , c('PROMIS_AnxietyTscore','PROMIS_DepressTscore')]), ]
# match using glm, optimal full matching, and average treatment effect (ATE)
psa_full_corr_N2_FCz <- matchit(RxCoded ~ Sex + PROMIS_AnxietyTscore + PROMIS_DepressTscore + Dx + Age + ss_meanrt_go, data = corrN2_FCz_outliersRemoved_MDD, distance = "glm", method = "full", estimand = "ATE")
summary(psa_full_corr_N2_FCz)

## 
## Call:
## matchit(formula = RxCoded ~ Sex + PROMIS_AnxietyTscore + PROMIS_DepressTscore + 
##     Dx + Age + ss_meanrt_go, data = corrN2_FCz_outliersRemoved_MDD, 
##     method = "full", distance = "glm", estimand = "ATE")
## 
## Summary of Balance for All Data:
##                      Means Treated Means Control Std. Mean Diff. Var. Ratio
## distance                    0.4953        0.4613          0.3731     0.7220
## SexFemale                   0.7647        0.6452          0.2644          .
## SexMale                     0.2353        0.3548         -0.2644          .
## PROMIS_AnxietyTscore       62.4847       62.8129         -0.0487     0.9679
## PROMIS_DepressTscore       61.1506       61.6376         -0.0663     0.9800
## DxCOM                       0.6941        0.6774          0.0360          .
## DxMDD                       0.3059        0.3226         -0.0360          .
## Age                        37.1681       34.5722          0.2288     0.9285
## ss_meanrt_go                0.8004        0.7920          0.0652     0.9624
##                      eCDF Mean eCDF Max
## distance                0.0988   0.2096
## SexFemale               0.1195   0.1195
## SexMale                 0.1195   0.1195
## PROMIS_AnxietyTscore    0.0256   0.0949
## PROMIS_DepressTscore    0.0295   0.1001
## DxCOM                   0.0167   0.0167
## DxMDD                   0.0167   0.0167
## Age                     0.0692   0.1908
## ss_meanrt_go            0.0362   0.0912
## 
## Summary of Balance for Matched Data:
##                      Means Treated Means Control Std. Mean Diff. Var. Ratio
## distance                    0.4787        0.4775          0.0133     0.9571
## SexFemale                   0.6966        0.7135         -0.0373          .
## SexMale                     0.3034        0.2865          0.0373          .
## PROMIS_AnxietyTscore       63.0637       62.7707          0.0435     1.1017
## PROMIS_DepressTscore       62.1073       61.6264          0.0655     1.1477
## DxCOM                       0.7291        0.6603          0.1482          .
## DxMDD                       0.2709        0.3397         -0.1482          .
## Age                        35.9189       35.5428          0.0332     0.8656
## ss_meanrt_go                0.8026        0.8023          0.0022     0.9848
##                      eCDF Mean eCDF Max Std. Pair Dist.
## distance                0.0104   0.0562          0.0488
## SexFemale               0.0169   0.0169          0.3623
## SexMale                 0.0169   0.0169          0.3623
## PROMIS_AnxietyTscore    0.0404   0.1208          1.1382
## PROMIS_DepressTscore    0.0240   0.0712          1.1864
## DxCOM                   0.0688   0.0688          0.9658
## DxMDD                   0.0688   0.0688          0.9658
## Age                     0.0310   0.0847          0.7747
## ss_meanrt_go            0.0370   0.1327          1.0340
## 
## Sample Sizes:
##               Control Treated
## All             93.     85.  
## Matched (ESS)   76.25   64.11
## Matched         93.     85.  
## Unmatched        0.      0.  
## Discarded        0.      0.

lovePlotCorrN2_FCz <- love.plot(bal.tab(psa_full_corr_N2_FCz),
                                stat = c("m","v"),
                                grid = TRUE,
                                stars = "std",
                                thresholds = c(m=0.10,v=1.10))
lovePlotCorrN2_FCz

psa_full_matchedDataCorrN2_FCz <- match.data(psa_full_corr_N2_FCz)

# weighted regression with weights from matched data
weightedModelCorrN2FCz_MDD <- lm(N2_FCz_AllCorrectAvg ~ RxCoded + Sex + PROMIS_AnxietyTscore + PROMIS_DepressTscore  + Dx + Age + ss_meanrt_go, data = psa_full_matchedDataCorrN2_FCz, weights = weights)
summary(weightedModelCorrN2FCz_MDD) # Rx unqiuely predicts correct-stop N200 amplitudes

## 
## Call:
## lm(formula = N2_FCz_AllCorrectAvg ~ RxCoded + Sex + PROMIS_AnxietyTscore + 
##     PROMIS_DepressTscore + Dx + Age + ss_meanrt_go, data = psa_full_matchedDataCorrN2_FCz, 
##     weights = weights)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.8662  -2.6360   0.0676   2.6258  10.3044 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -2.0923887  4.4479046  -0.470 0.638657    
## RxCoded               2.4650831  0.6724645   3.666 0.000329 ***
## SexMale              -1.3102779  0.7560764  -1.733 0.084910 .  
## PROMIS_AnxietyTscore -0.0179323  0.0707494  -0.253 0.800217    
## PROMIS_DepressTscore -0.0000636  0.0567516  -0.001 0.999107    
## DxMDD                 1.1754514  0.8146077   1.443 0.150870    
## Age                  -0.0310329  0.0298293  -1.040 0.299655    
## ss_meanrt_go         -0.7550471  2.6664714  -0.283 0.777396    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.464 on 170 degrees of freedom
## Multiple R-squared:  0.09856,    Adjusted R-squared:  0.06144 
## F-statistic: 2.655 on 7 and 170 DF,  p-value: 0.01245

# code for information about SE, CI, and p-values
vcov_cr <- vcovCL(
  weightedModelCorrN2FCz_MDD,
  cluster = psa_full_matchedDataCorrN2_FCz$subclass
)
coeftest(weightedModelCorrN2FCz_MDD, vcov = vcov_cr)

## 
## t test of coefficients:
## 
##                         Estimate  Std. Error t value Pr(>|t|)   
## (Intercept)          -2.0924e+00  4.6782e+00 -0.4473 0.655251   
## RxCoded               2.4651e+00  8.2490e-01  2.9883 0.003221 **
## SexMale              -1.3103e+00  8.1805e-01 -1.6017 0.111079   
## PROMIS_AnxietyTscore -1.7932e-02  7.4898e-02 -0.2394 0.811066   
## PROMIS_DepressTscore -6.3597e-05  6.2432e-02 -0.0010 0.999188   
## DxMDD                 1.1755e+00  9.3976e-01  1.2508 0.212724   
## Age                  -3.1033e-02  2.9714e-02 -1.0444 0.297791   
## ss_meanrt_go         -7.5505e-01  2.5510e+00 -0.2960 0.767603   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tidyWeightedModel_CorrN2_FCz_MDD <- tidy(
  weightedModelCorrN2FCz_MDD,
  vcov = vcov_cr,
  conf.int = TRUE
)
tidyWeightedModel_CorrN2_FCz_MDD

## # A tibble: 8 × 7
##   term                   estimate std.error statistic p.value conf.low conf.high
##   <chr>                     <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)          -2.09         4.45    -0.470   6.39e-1 -10.9       6.69  
## 2 RxCoded               2.47         0.672    3.67    3.29e-4   1.14      3.79  
## 3 SexMale              -1.31         0.756   -1.73    8.49e-2  -2.80      0.182 
## 4 PROMIS_AnxietyTscore -0.0179       0.0707  -0.253   8.00e-1  -0.158     0.122 
## 5 PROMIS_DepressTscore -0.0000636    0.0568  -0.00112 9.99e-1  -0.112     0.112 
## 6 DxMDD                 1.18         0.815    1.44    1.51e-1  -0.433     2.78  
## 7 Age                  -0.0310       0.0298  -1.04    3.00e-1  -0.0899    0.0279
## 8 ss_meanrt_go         -0.755        2.67    -0.283   7.77e-1  -6.02      4.51

# ANOVA
aovCorrN2FCz_MDD <- corrN2_FCz_outliersRemoved_MDD %>% 
  anova_test(N2_FCz_AllCorrectAvg ~ Rx * Sex * Dx, type = 3) %>%
  adjust_pvalue(p.col = "p", output.col = "p.adj", method = "bonferroni")
aovCorrN2FCz_MDD # Rx and Sex predict N200 amplitudes, but Rx survies bonferroni correction

## ANOVA Table (type III tests)
## 
##      Effect DFn DFd     F     p p<.05      ges p.adj
## 1        Rx   1 170 7.880 0.006     * 4.40e-02 0.042
## 2       Sex   1 170 6.652 0.011     * 3.80e-02 0.077
## 3        Dx   1 170 0.645 0.423       4.00e-03 1.000
## 4    Rx:Sex   1 170 0.754 0.386       4.00e-03 1.000
## 5     Rx:Dx   1 170 0.082 0.775       4.84e-04 1.000
## 6    Sex:Dx   1 170 0.009 0.923       5.54e-05 1.000
## 7 Rx:Sex:Dx   1 170 0.143 0.706       8.38e-04 1.000

# adding on HC participants from corrN2_FCz_outliersRemoved_HC into dataframe
corrN2_FCz_outliersRemoved_HConly <- corrN2_FCz_outliersRemoved_HC %>%
  filter(Dx == "HC")
corrN2_FCz_outliersRemoved_MDDvHC <- bind_rows(corrN2_FCz_outliersRemoved_HConly,corrN2_FCz_outliersRemoved_MDD)
corrN2_FCz_outliersRemoved_MDDvHC <- corrN2_FCz_outliersRemoved_MDDvHC %>%
  mutate(Rx = case_when(
    Dx == "HC" ~ "Non",
    Rx == "NonSero" ~ "Unmed",
    Rx == "Sero" ~ "Med"
  ))
# creating violin plot of all three groups
vioPlotCorrN200FCz_MDD <- ggplot(corrN2_FCz_outliersRemoved_MDDvHC, aes (x = Rx, y = N2_FCz_AllCorrectAvg, fill = Rx)) + geom_violin(trim = FALSE) + labs(x = "Medication", y = "Amplitude (μV)")
vioPlotCorrN200FCz_MDD <- vioPlotCorrN200FCz_MDD + geom_jitter(shape = 16, position = position_jitter(0.2)) + theme_classic() + scale_fill_manual(values = c("#5f6266","#a5a9ae","#f8f9fa"), guide = "none")

# creating dataframe with summary statistics 
corrN2_FCz_MDD_Sum <- corrN2_FCz_outliersRemoved_MDDvHC %>%
  group_by(Rx) %>%
  summarize(mean=mean(N2_FCz_AllCorrectAvg),
            se=sd(N2_FCz_AllCorrectAvg)/sqrt(length(N2_FCz_AllCorrectAvg)))
# creating bar graph (with error bars)
corrN2_FCz_MDD_barGraph <- ggplot(corrN2_FCz_MDD_Sum, aes(x = Rx, y = mean, fill = Rx)) + 
  geom_col(colour = 'black') + 
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.2) + theme_classic() + scale_fill_manual(values = c("#5f6266","#a5a9ae","#f8f9fa"), guide = "none")

patchwork_corrN2_FCz_MDD <- qqPlot_N2_FCz_MDD +
  corrN2_FCz_MDD_barGraph + 
  vioPlotCorrN200FCz_MDD
patchwork_corrN2_FCz_MDD

7.2 Correct Stop-P300 Amplitudes for MDD/COM Participants

qqPlot_corrP3_Pz_MDD <- ggplot(corrP3_Pz_outliersRemoved_MDD, aes(sample = P3_Pz_AllCorrectAvg)) + 
  geom_qq() +
  geom_qq_line()
qqPlot_corrP3_Pz_MDD <- qqPlot_corrP3_Pz_MDD + ylab("Sample") + xlab("Theoretical") + theme_classic()

# using https://www.youtube.com/watch?v=7vwewIsSAls for learning this technique
corrP3_Pz_outliersRemoved_MDD$RxCoded <- str_replace_all(corrP3_Pz_outliersRemoved_MDD$RxMofA, c("None"="0","SSRI"="1","SNRI"="1"))
corrP3_Pz_outliersRemoved_MDD$RxCoded <- as.numeric(corrP3_Pz_outliersRemoved_MDD$RxCoded)
# remove participants who have incomplete PROMIS measures (or other covariates)
corrP3_Pz_outliersRemoved_MDD <- corrP3_Pz_outliersRemoved_MDD[complete.cases(corrP3_Pz_outliersRemoved_MDD[ , c('PROMIS_AnxietyTscore','PROMIS_DepressTscore')]), ]
# match using glm, optimal full matching, and average treatment effect (ATE)
psa_full_corr_P3_Pz <- matchit(RxCoded ~ Sex + PROMIS_AnxietyTscore + PROMIS_DepressTscore + Dx + Age + ss_meanrt_go, data = corrP3_Pz_outliersRemoved_MDD, distance = "glm", method = "full", estimand = "ATE")
summary(psa_full_corr_P3_Pz)

## 
## Call:
## matchit(formula = RxCoded ~ Sex + PROMIS_AnxietyTscore + PROMIS_DepressTscore + 
##     Dx + Age + ss_meanrt_go, data = corrP3_Pz_outliersRemoved_MDD, 
##     method = "full", distance = "glm", estimand = "ATE")
## 
## Summary of Balance for All Data:
##                      Means Treated Means Control Std. Mean Diff. Var. Ratio
## distance                    0.5036        0.4804          0.3062     0.7326
## SexFemale                   0.7667        0.6667          0.2233          .
## SexMale                     0.2333        0.3333         -0.2233          .
## PROMIS_AnxietyTscore       62.5044       62.6473         -0.0210     1.0401
## PROMIS_DepressTscore       61.3433       61.6602         -0.0428     1.0035
## DxCOM                       0.7000        0.6774          0.0488          .
## DxMDD                       0.3000        0.3226         -0.0488          .
## Age                        36.5713       34.5411          0.1774     0.9664
## ss_meanrt_go                0.8023        0.7943          0.0621     0.8981
##                      eCDF Mean eCDF Max
## distance                0.0826   0.1871
## SexFemale               0.1000   0.1000
## SexMale                 0.1000   0.1000
## PROMIS_AnxietyTscore    0.0254   0.0878
## PROMIS_DepressTscore    0.0314   0.1004
## DxCOM                   0.0226   0.0226
## DxMDD                   0.0226   0.0226
## Age                     0.0532   0.1620
## ss_meanrt_go            0.0343   0.0789
## 
## Summary of Balance for Matched Data:
##                      Means Treated Means Control Std. Mean Diff. Var. Ratio
## distance                    0.4924        0.4917          0.0085     0.9689
## SexFemale                   0.7149        0.7250         -0.0224          .
## SexMale                     0.2851        0.2750          0.0224          .
## PROMIS_AnxietyTscore       62.5081       62.5782         -0.0103     1.1191
## PROMIS_DepressTscore       61.6658       61.6993         -0.0045     1.0453
## DxCOM                       0.6931        0.6794          0.0295          .
## DxMDD                       0.3069        0.3206         -0.0295          .
## Age                        36.0001       35.1362          0.0755     0.9451
## ss_meanrt_go                0.7917        0.8030         -0.0875     1.0047
##                      eCDF Mean eCDF Max Std. Pair Dist.
## distance                0.0087   0.0383          0.0350
## SexFemale               0.0100   0.0100          0.3850
## SexMale                 0.0100   0.0100          0.3850
## PROMIS_AnxietyTscore    0.0306   0.1066          1.0757
## PROMIS_DepressTscore    0.0247   0.0856          1.1683
## DxCOM                   0.0137   0.0137          0.8381
## DxMDD                   0.0137   0.0137          0.8381
## Age                     0.0357   0.1257          0.8127
## ss_meanrt_go            0.0396   0.1603          1.1273
## 
## Sample Sizes:
##               Control Treated
## All             93.      90. 
## Matched (ESS)   79.11    75.8
## Matched         93.      90. 
## Unmatched        0.       0. 
## Discarded        0.       0.

lovePlotCorrP3_Pz <- love.plot(bal.tab(psa_full_corr_P3_Pz),
                                stat = c("m","v"),
                                grid = TRUE,
                                stars = "std",
                                thresholds = c(m=0.10,v=1.10))
lovePlotCorrP3_Pz

psa_full_matchedDataCorrP3_Pz <- match.data(psa_full_corr_P3_Pz)

# weighted regression with weights from matched data
weightedModelCorrP3Pz_MDD <- lm(P3_Pz_AllCorrectAvg ~ RxCoded + ss_meanrt_go + Sex + PROMIS_AnxietyTscore + PROMIS_DepressTscore  + Dx + Age, data = psa_full_matchedDataCorrP3_Pz, weights = weights)
summary(weightedModelCorrP3Pz_MDD) # Rx predicts correct-stop P300 amplitudes

## 
## Call:
## lm(formula = P3_Pz_AllCorrectAvg ~ RxCoded + ss_meanrt_go + Sex + 
##     PROMIS_AnxietyTscore + PROMIS_DepressTscore + Dx + Age, data = psa_full_matchedDataCorrP3_Pz, 
##     weights = weights)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.9321 -3.4396  0.2525  2.8984  9.5982 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)           6.37527    4.14286   1.539   0.1256  
## RxCoded               1.61104    0.62939   2.560   0.0113 *
## ss_meanrt_go         -1.66131    2.51244  -0.661   0.5093  
## SexMale              -0.03406    0.71188  -0.048   0.9619  
## PROMIS_AnxietyTscore -0.06295    0.06301  -0.999   0.3192  
## PROMIS_DepressTscore  0.07495    0.05496   1.364   0.1744  
## DxMDD                 0.32862    0.74518   0.441   0.6598  
## Age                  -0.05559    0.02793  -1.991   0.0481 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.248 on 175 degrees of freedom
## Multiple R-squared:  0.06925,    Adjusted R-squared:  0.03202 
## F-statistic:  1.86 on 7 and 175 DF,  p-value: 0.0788

# code for information about SE, CI, and p-values
vcov_cr <- vcovCL(
  weightedModelCorrP3Pz_MDD,
  cluster = psa_full_matchedDataCorrP3_Pz$subclass
)
coeftest(weightedModelCorrP3Pz_MDD, vcov = vcov_cr)

## 
## t test of coefficients:
## 
##                       Estimate Std. Error t value Pr(>|t|)  
## (Intercept)           6.375274   4.593405  1.3879  0.16693  
## RxCoded               1.611038   0.660800  2.4380  0.01577 *
## ss_meanrt_go         -1.661308   2.819925 -0.5891  0.55653  
## SexMale              -0.034060   0.634576 -0.0537  0.95726  
## PROMIS_AnxietyTscore -0.062948   0.075423 -0.8346  0.40508  
## PROMIS_DepressTscore  0.074948   0.063852  1.1738  0.24207  
## DxMDD                 0.328620   0.648247  0.5069  0.61284  
## Age                  -0.055591   0.027099 -2.0514  0.04172 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tidyWeightedModel_CorrP3_Pz_MDD <- tidy(
  weightedModelCorrP3Pz_MDD,
  vcov = vcov_cr,
  conf.int = TRUE
)
tidyWeightedModel_CorrP3_Pz_MDD

## # A tibble: 8 × 7
##   term                 estimate std.error statistic p.value conf.low conf.high
##   <chr>                   <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)            6.38      4.14      1.54    0.126   -1.80   14.6     
## 2 RxCoded                1.61      0.629     2.56    0.0113   0.369   2.85    
## 3 ss_meanrt_go          -1.66      2.51     -0.661   0.509   -6.62    3.30    
## 4 SexMale               -0.0341    0.712    -0.0478  0.962   -1.44    1.37    
## 5 PROMIS_AnxietyTscore  -0.0629    0.0630   -0.999   0.319   -0.187   0.0614  
## 6 PROMIS_DepressTscore   0.0749    0.0550    1.36    0.174   -0.0335  0.183   
## 7 DxMDD                  0.329     0.745     0.441   0.660   -1.14    1.80    
## 8 Age                   -0.0556    0.0279   -1.99    0.0481  -0.111  -0.000474

# ANOVA
aovCorrP3Pz_MDD <- corrP3_Pz_outliersRemoved_MDD %>% 
  anova_test(P3_Pz_AllCorrectAvg ~ Rx * Sex * Dx, type = 3) %>%
  adjust_pvalue(p.col = "p", output.col = "p.adj", method = "bonferroni")
aovCorrP3Pz_MDD # no significant effects

## ANOVA Table (type III tests)
## 
##      Effect DFn DFd     F     p p<.05      ges p.adj
## 1        Rx   1 175 0.102 0.750       0.000584 1.000
## 2       Sex   1 175 0.061 0.806       0.000347 1.000
## 3        Dx   1 175 0.155 0.694       0.000885 1.000
## 4    Rx:Sex   1 175 3.075 0.081       0.017000 0.567
## 5     Rx:Dx   1 175 0.051 0.822       0.000290 1.000
## 6    Sex:Dx   1 175 0.031 0.860       0.000177 1.000
## 7 Rx:Sex:Dx   1 175 0.999 0.319       0.006000 1.000

# adding on HC participants from corrN2_FCz_outliersRemoved_HC into dataframe
corrP3_Pz_outliersRemoved_HConly <- corrP3_Pz_outliersRemoved_HC %>%
  filter(Dx == "HC")
corrP3_Pz_outliersRemoved_MDDvHC <- bind_rows(corrP3_Pz_outliersRemoved_HConly,corrP3_Pz_outliersRemoved_MDD)
corrP3_Pz_outliersRemoved_MDDvHC <- corrP3_Pz_outliersRemoved_MDDvHC %>%
  mutate(Rx = case_when(
    Dx == "HC" ~ "Non",
    Rx == "NonSero" ~ "Unmed",
    Rx == "Sero" ~ "Med"
  ))
# violin plot per Rx (to demonstrate propensity matching results)
vioPlotCorrP300Pz_MDD <- ggplot(corrP3_Pz_outliersRemoved_MDDvHC, aes (x = Rx, y = P3_Pz_AllCorrectAvg, fill = Rx)) + geom_violin(trim = FALSE) + labs(x = "Rx", y = "Amplitude (μV)") +  scale_fill_manual(values = c("#5f6266","#a5a9ae","#f8f9fa"), guide = "none")
vioPlotCorrP300Pz_MDD <- vioPlotCorrP300Pz_MDD + geom_jitter(shape = 16, position = position_jitter(0.2)) + theme_classic()

# creating dataframe with summary statistics 
corrP3_Pz_MDD_Sum <- corrP3_Pz_outliersRemoved_MDDvHC %>%
  group_by(Rx) %>%
  summarize(mean=mean(P3_Pz_AllCorrectAvg),
            se=sd(P3_Pz_AllCorrectAvg)/sqrt(length(P3_Pz_AllCorrectAvg)))
# creating bar graph (with error bars)
corrP3_Pz_MDD_barGraph <- ggplot(corrP3_Pz_MDD_Sum, aes(x = Rx, y = mean, fill = Rx)) + 
  geom_col(colour = 'black') + 
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.2) + theme_classic() + scale_fill_manual(values = c("#5f6266","#a5a9ae","#f8f9fa"), guide = "none")

patchwork_corrP3_Pz_MDD <- qqPlot_corrP3_Pz_MDD + 
corrP3_Pz_MDD_barGraph + vioPlotCorrP300Pz_MDD + plot_layout(ncol = 3)
patchwork_corrP3_Pz_MDD # 900 w x 400 h 

7.3 CRN Amplitudes for MDD/COM Participants

qqPlot_CRN_FCz_MDD <- ggplot(CRN_FCz_outliersRemoved_MDD, aes(sample = ERN_FCz_GoTargetRespAvg)) + 
  geom_qq() +
  geom_qq_line()
qqPlot_CRN_FCz_MDD <- qqPlot_CRN_FCz_MDD + ylab("Sample") + xlab("Theoretical") + ggtitle("QQ Plot of CRN Mean Amplitudes @ FCz in MDD/COM Participants") + theme_classic()

# using https://www.youtube.com/watch?v=7vwewIsSAls for learning this technique
CRN_FCz_outliersRemoved_MDD$RxCoded <- str_replace_all(CRN_FCz_outliersRemoved_MDD$RxMofA, c("None"="0","SSRI"="1","SNRI"="1"))
CRN_FCz_outliersRemoved_MDD$RxCoded <- as.numeric(CRN_FCz_outliersRemoved_MDD$RxCoded)
# remove participants who have incomplete PROMIS measures (or other covariates)
CRN_FCz_outliersRemoved_MDD <- CRN_FCz_outliersRemoved_MDD[complete.cases(CRN_FCz_outliersRemoved_MDD[ , c('PROMIS_AnxietyTscore','PROMIS_DepressTscore')]), ]
# match using glm, optimal full matching, and average treatment effect (ATE)
psa_full_CRN_FCz <- matchit(RxCoded ~ Sex + PROMIS_AnxietyTscore  + PROMIS_DepressTscore + Dx + Age + ss_meanrt_go, data = CRN_FCz_outliersRemoved_MDD, distance = "glm", method = "full", estimand = "ATE")
summary(psa_full_CRN_FCz)

## 
## Call:
## matchit(formula = RxCoded ~ Sex + PROMIS_AnxietyTscore + PROMIS_DepressTscore + 
##     Dx + Age + ss_meanrt_go, data = CRN_FCz_outliersRemoved_MDD, 
##     method = "full", distance = "glm", estimand = "ATE")
## 
## Summary of Balance for All Data:
##                      Means Treated Means Control Std. Mean Diff. Var. Ratio
## distance                    0.5064        0.4661          0.4058     0.6623
## SexFemale                   0.8000        0.6778          0.2810          .
## SexMale                     0.2000        0.3222         -0.2810          .
## PROMIS_AnxietyTscore       62.4929       62.9478         -0.0687     0.9314
## PROMIS_DepressTscore       61.3059       61.8344         -0.0727     0.8791
## DxCOM                       0.7176        0.6778          0.0869          .
## DxMDD                       0.2824        0.3222         -0.0869          .
## Age                        36.8536       34.5843          0.1992     0.9677
## ss_meanrt_go                0.7987        0.7899          0.0694     0.9441
##                      eCDF Mean eCDF Max
## distance                0.1020   0.1902
## SexFemale               0.1222   0.1222
## SexMale                 0.1222   0.1222
## PROMIS_AnxietyTscore    0.0327   0.1105
## PROMIS_DepressTscore    0.0396   0.1216
## DxCOM                   0.0399   0.0399
## DxMDD                   0.0399   0.0399
## Age                     0.0582   0.1725
## ss_meanrt_go            0.0337   0.0863
## 
## Summary of Balance for Matched Data:
##                      Means Treated Means Control Std. Mean Diff. Var. Ratio
## distance                    0.4872        0.4853          0.0187     0.9370
## SexFemale                   0.7048        0.7457         -0.0942          .
## SexMale                     0.2952        0.2543          0.0942          .
## PROMIS_AnxietyTscore       63.0101       63.2026         -0.0291     1.0546
## PROMIS_DepressTscore       62.3262       62.0735          0.0348     1.1169
## DxCOM                       0.7196        0.6937          0.0566          .
## DxMDD                       0.2804        0.3063         -0.0566          .
## Age                        36.4052       35.7593          0.0567     0.8771
## ss_meanrt_go                0.8095        0.7981          0.0900     0.8605
##                      eCDF Mean eCDF Max Std. Pair Dist.
## distance                0.0102   0.0514          0.0503
## SexFemale               0.0410   0.0410          0.4481
## SexMale                 0.0410   0.0410          0.4481
## PROMIS_AnxietyTscore    0.0366   0.0857          1.1109
## PROMIS_DepressTscore    0.0338   0.0776          1.0803
## DxCOM                   0.0260   0.0260          0.8312
## DxMDD                   0.0260   0.0260          0.8312
## Age                     0.0347   0.1092          0.8814
## ss_meanrt_go            0.0353   0.0976          1.0215
## 
## Sample Sizes:
##               Control Treated
## All             90.     85.  
## Matched (ESS)   72.93   62.42
## Matched         90.     85.  
## Unmatched        0.      0.  
## Discarded        0.      0.

lovePlotCRN_FCz <- love.plot(bal.tab(psa_full_CRN_FCz),
                                stat = c("m","v"),
                                grid = TRUE,
                                stars = "std",
                                thresholds = c(m=0.10,v=1.10))
lovePlotCRN_FCz

psa_full_matchedDataCRN_FCz <- match.data(psa_full_CRN_FCz)

# weighted regression with weights from matched data
weightedModelCRN_FCz_MDD <- lm(ERN_FCz_GoTargetRespAvg ~ RxCoded + Sex + PROMIS_AnxietyTscore + PROMIS_DepressTscore + Dx + Age + ss_meanrt_go, data = psa_full_matchedDataCRN_FCz, weights = weights)
summary(weightedModelCRN_FCz_MDD) # nothing predicts CRN amplitudes

## 
## Call:
## lm(formula = ERN_FCz_GoTargetRespAvg ~ RxCoded + Sex + PROMIS_AnxietyTscore + 
##     PROMIS_DepressTscore + Dx + Age + ss_meanrt_go, data = psa_full_matchedDataCRN_FCz, 
##     weights = weights)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2205 -1.3646  0.1489  1.5111  6.1243 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)  
## (Intercept)          -1.829044   2.412888  -0.758   0.4495  
## RxCoded              -0.499347   0.365916  -1.365   0.1742  
## SexMale               0.342381   0.417848   0.819   0.4137  
## PROMIS_AnxietyTscore  0.009894   0.038203   0.259   0.7960  
## PROMIS_DepressTscore -0.040569   0.030958  -1.310   0.1918  
## DxMDD                -0.240252   0.461222  -0.521   0.6031  
## Age                   0.004145   0.016204   0.256   0.7984  
## ss_meanrt_go          2.885476   1.466317   1.968   0.0507 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.41 on 167 degrees of freedom
## Multiple R-squared:  0.04654,    Adjusted R-squared:  0.006575 
## F-statistic: 1.165 on 7 and 167 DF,  p-value: 0.3258

# ANOVA
# code for information about SE, CI, and p-values
vcov_cr <- vcovCL(
  weightedModelCRN_FCz_MDD,
  cluster = psa_full_matchedDataCRN_FCz$subclass
)
coeftest(weightedModelCRN_FCz_MDD, vcov = vcov_cr)

## 
## t test of coefficients:
## 
##                        Estimate Std. Error t value Pr(>|t|)  
## (Intercept)          -1.8290440  2.8106585 -0.6508  0.51610  
## RxCoded              -0.4993473  0.4404630 -1.1337  0.25855  
## SexMale               0.3423813  0.4714447  0.7262  0.46871  
## PROMIS_AnxietyTscore  0.0098939  0.0456101  0.2169  0.82853  
## PROMIS_DepressTscore -0.0405690  0.0328077 -1.2366  0.21798  
## DxMDD                -0.2402516  0.3701077 -0.6491  0.51714  
## Age                   0.0041447  0.0177861  0.2330  0.81602  
## ss_meanrt_go          2.8854763  1.4658935  1.9684  0.05068 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tidyWeightedModel_CRN_FCz_MDD <- tidy(
  weightedModelCRN_FCz_MDD,
  vcov = vcov_cr,
  conf.int = TRUE
)
tidyWeightedModel_CRN_FCz_MDD

## # A tibble: 8 × 7
##   term                 estimate std.error statistic p.value conf.low conf.high
##   <chr>                   <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)          -1.83       2.41      -0.758  0.450  -6.59       2.93  
## 2 RxCoded              -0.499      0.366     -1.36   0.174  -1.22       0.223 
## 3 SexMale               0.342      0.418      0.819  0.414  -0.483      1.17  
## 4 PROMIS_AnxietyTscore  0.00989    0.0382     0.259  0.796  -0.0655     0.0853
## 5 PROMIS_DepressTscore -0.0406     0.0310    -1.31   0.192  -0.102      0.0206
## 6 DxMDD                -0.240      0.461     -0.521  0.603  -1.15       0.670 
## 7 Age                   0.00414    0.0162     0.256  0.798  -0.0278     0.0361
## 8 ss_meanrt_go          2.89       1.47       1.97   0.0507 -0.00943    5.78

aovCorrP3Pz_MDD <- CRN_FCz_outliersRemoved_MDD %>% 
  anova_test(ERN_FCz_GoTargetRespAvg ~ Rx * Sex * Dx, type = 3) %>%
  adjust_pvalue(p.col = "p", output.col = "p.adj", method = "bonferroni")
aovCorrP3Pz_MDD # nothing survives bonferroni correction

## ANOVA Table (type III tests)
## 
##      Effect DFn DFd     F     p p<.05      ges p.adj
## 1        Rx   1 167 0.035 0.853       2.08e-04 1.000
## 2       Sex   1 167 1.144 0.286       7.00e-03 1.000
## 3        Dx   1 167 0.115 0.735       6.88e-04 1.000
## 4    Rx:Sex   1 167 7.097 0.008     * 4.10e-02 0.056
## 5     Rx:Dx   1 167 0.002 0.969       9.32e-06 1.000
## 6    Sex:Dx   1 167 0.028 0.867       1.69e-04 1.000
## 7 Rx:Sex:Dx   1 167 1.393 0.240       8.00e-03 1.000

# adding on HC participants from corrN2_FCz_outliersRemoved_HC into dataframe
CRN_FCz_outliersRemoved_HConly <- CRN_FCz_outliersRemoved_HC %>%
  filter(Dx == "HC")
CRN_FCz_outliersRemoved_MDDvHC <- bind_rows(CRN_FCz_outliersRemoved_HConly,CRN_FCz_outliersRemoved_MDD)
CRN_FCz_outliersRemoved_MDDvHC <- CRN_FCz_outliersRemoved_MDDvHC %>%
  mutate(Rx = case_when(
    Dx == "HC" ~ "Non",
    Rx == "NonSero" ~ "Unmed",
    Rx == "Sero" ~ "Med"
  ))
# violin plot 
vioPlotCRNFCz_MDD <- ggplot(CRN_FCz_outliersRemoved_MDDvHC, aes (x = Rx, y = ERN_FCz_GoTargetRespAvg, fill = Rx)) + geom_violin(trim = FALSE) + labs(x = "Medications", y = "Amplitude (μV)")
vioPlotCRNFCz_MDD + geom_jitter(shape = 16, position = position_jitter(0.2)) + theme_classic() + scale_fill_manual(values = c("#5f6266","#a5a9ae","#f8f9fa"), guide = "none")

# creating dataframe with summary statistics 
CRN_FCz_MDD_Sum <- CRN_FCz_outliersRemoved_MDDvHC %>%
  group_by(Rx) %>%
  summarize(mean=mean(ERN_FCz_GoTargetRespAvg),
            se=sd(ERN_FCz_GoTargetRespAvg)/sqrt(length(ERN_FCz_GoTargetRespAvg)))
# creating bar graph (with error bars)
CRN_FCz_MDD_barGraph <- ggplot(CRN_FCz_MDD_Sum, aes(x = Rx, y = mean, fill = Rx)) + 
  geom_col(colour = 'black') + 
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.1) + theme_classic() + scale_fill_manual(values = c("#5f6266","#a5a9ae","#f8f9fa"), guide = "none")