Ross AB Data

1. Ross Attentional Bias (AB) Data

1. Load Packages Used in Data Analysis

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.1
## ✔ readr   2.1.2     ✔ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(readxl)
library(itrak)
library(rstatix)

## 
## Attaching package: 'rstatix'

## The following object is masked from 'package:itrak':
## 
##     is_outlier

## The following object is masked from 'package:stats':
## 
##     filter

library(afex)

## Loading required package: lme4

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## ************
## Welcome to afex. For support visit: http://afex.singmann.science/

## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - NEWS: emmeans() for ANOVA models now uses model = 'multivariate' as default.
## - Get and set global package options with: afex_options()
## - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************

## 
## Attaching package: 'afex'

## The following object is masked from 'package:lme4':
## 
##     lmer

library(grid)
library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

library(extrafont)

## Registering fonts with R

loadfonts(device = "postscript")

2. Data Read-In

knitr::opts_chunk$set(echo = TRUE)
#Note: You'll need to change this to match the file path for 
#whatever computer you're on. Also, if you're on pc, the file 
#path "/'s" will need to be changed
df <- readxl::read_excel(path = "/Users/noahwolkowicz/Desktop/CT/West Haven/Postdoc/Postdoc Research/Ross Attentional Bias Study/Data/Ross Data/AB_primary_clean_final (1).xlsx") %>% 
                        mutate(ExperimentName = factor(ExperimentName), 
                               #change data types from character strings to factors
                               Cue = factor(Cue),
                               Probe = factor(Probe),
                               Probe.tf = ifelse(Probe == "congruent", TRUE, FALSE))
glimpse(df) #view data frame

## Rows: 13,120
## Columns: 8
## $ ExperimentName <fct> baseline, baseline, baseline, baseline, baseline, basel…
## $ Subject        <dbl> 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, …
## $ Trial          <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ Accuracy       <dbl> 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1…
## $ RT             <dbl> 1422, 1204, 1277, 2387, 1255, 1369, 1105, 1607, 3166, 7…
## $ Cue            <fct> opioid, pain, opioid, opioid, opioid, pain, pain, opioi…
## $ Probe          <fct> congruent, congruent, congruent, incongruent, congruent…
## $ Probe.tf       <lgl> TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE,…

3. Data Pre-processing

Probe accuracy, incorrect trials, improbable reaction times

-The code below when run illustrates that there are no administrations #With <50% accuracy. Nothing excluded here. No cues at administrations with < 50% accuracy for any timepoint

df %>% 
  group_by(Subject, ExperimentName, Cue) %>% 
  summarise(Percent_Accurate = mean(Accuracy)) %>% 
  filter(Percent_Accurate < 0.50) 

## `summarise()` has grouped output by 'Subject', 'ExperimentName'. You can
## override using the `.groups` argument.

-The code below indicates the overall # of trials and the counts/percents for 1) incorrect trials to be excluded, and 2) improbable reaction time trials (i.e., <200ms OR >1500ms) to be excluded. The code immediately following then excludes this data.

df %>% 
  group_by(ExperimentName) %>% 
  summarise(Overall_Trials = n(), #Overall trial counts
            Count_Inc_Trial = sum(ifelse(Accuracy == 0,1,0)), #Incorrect trial count
            Per_Inc_Trial = 1 - mean(Accuracy), #Percent of trials that are incorrect
            Count_Inc_RT = sum(ifelse(RT < 200 | RT > 1500, 1, 0)), 
            #Count of trials whose reaction times are < 200 or > 1500
            Per_Inc_RT = 1 - mean(ifelse(RT < 200 | RT > 1500, 0, 1)))

-Filtering out improbable reaction times

df <- df %>% 
  filter(Accuracy == 1) %>%  #Exclude incorrect trials
  filter(RT > 200 & RT < 1500) #Exclude RT's <200ms or >1500ms

-Maximum RT after exclusions

max(df$RT)

## [1] 1495

-Minimum RT after exclusions

min(df$RT)

## [1] 273

-Code for Median Absolute Deviation and SD calculations for Reaction Times hidden

#The code here calculates Median Absolute Deviation (MADS) and standard #deviations for outlier screening.
MAD_SD_df <- df %>% 
  group_by(Subject) %>% 
  summarise(med = median(RT), #1) Gets median RT for each subject, 
            MAD = mad(RT), #2) gets MAD for each subject, 
            MAD3 = 3*MAD, #3) calculates 3*MAD for each subject,
            MAD_Exclude_Score = med + MAD3, #4) adds the 3*MAD to each subject's median (i.e., the RT score at which a trial would be excluded);
            Mean = mean(RT), #5) it calculates mean reaction time, 
            SD = sd(RT), #6) gets standard deviation of reaction times
            SD3 = 3*SD, #7) gets 3x the standard deviation of reaction times
            SD_Exclude_Score = SD3 + Mean) #8) adds the 3*SD to each subject's mean reaction time (i.e., the RT score at which a trial would be excluded);
glimpse(MAD_SD_df)

## Rows: 29
## Columns: 9
## $ Subject           <dbl> 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 11…
## $ med               <dbl> 547.0, 546.5, 687.0, 667.0, 531.0, 697.5, 520.0, 611…
## $ MAD               <dbl> 109.7124, 82.2843, 232.7682, 149.7426, 93.4038, 115.…
## $ MAD3              <dbl> 329.1372, 246.8529, 698.3046, 449.2278, 280.2114, 34…
## $ MAD_Exclude_Score <dbl> 876.1372, 793.3529, 1385.3046, 1116.2278, 811.2114, …
## $ Mean              <dbl> 573.2624, 570.3297, 747.9965, 727.1134, 566.8553, 72…
## $ SD                <dbl> 138.37599, 133.59148, 245.55486, 224.86856, 156.1705…
## $ SD3               <dbl> 415.1280, 400.7744, 736.6646, 674.6057, 468.5116, 45…
## $ SD_Exclude_Score  <dbl> 988.3903, 971.1042, 1484.6611, 1401.7191, 1035.3669,…

df <- left_join(df, MAD_SD_df, by = "Subject") #Binds the variables calculated above to the data

#The code below shows overall trials (after excluding for accuracy, RT as above),
#and then shows the #/% of trials that will be excluded based on MAD. The code 
#immediately below then excludes these folks.

-MAD & SD Calculation Results

df %>% 
  group_by(ExperimentName) %>% 
  summarise(Ovr_Trl_Aftr_Exc = n(), #Overall trial counts after excluding for above-noted criteria
            Count_MAD3_Over = sum(ifelse(RT > MAD_Exclude_Score, 1,0)), #Count of trials that had RT > 3MAD above individual median
            Per_MAD3_Over = mean(ifelse(RT > MAD_Exclude_Score, 1, 0)),
            Count_SD3_Over = sum(ifelse(RT > SD_Exclude_Score, 1,0)), #Count of trials that had RT > 3MAD above individual median
            Per_SD3_Over = mean(ifelse(RT > SD_Exclude_Score, 1, 0)))

-Excluding Folks based on SD for Reaction Times

df <- df %>% 
  filter(RT < SD_Exclude_Score) #Exclude folks whose RT are >3SD above their individual means

4. Attentional Bias (AB) Calculations

#The code below first orders the data in numerical order by 1) subject, then 2) trial #.
#Then, it groups the data per 1) subject, then 2) administration, and finally 3) cue type.
#Lastly, it runs the summarize_bias() function to calculate mean and trial-level bias scores,
#Which are first stored in a separate dataframe and then bound to the original data
names(df)

##  [1] "ExperimentName"    "Subject"           "Trial"            
##  [4] "Accuracy"          "RT"                "Cue"              
##  [7] "Probe"             "Probe.tf"          "med"              
## [10] "MAD"               "MAD3"              "MAD_Exclude_Score"
## [13] "Mean"              "SD"                "SD3"              
## [16] "SD_Exclude_Score"

AB_df <- df %>% 
  arrange(Subject, Trial) %>% 
  group_by(Subject, ExperimentName, Cue) %>% 
  itrak::summarize_bias(data = ., RT = RT, congruent = Probe.tf, method = "nearest")

df <- left_join(df, AB_df, by = c("ExperimentName", "Cue", "Subject"))
Grouped_Output <- df %>% 
  group_by(Subject, ExperimentName, Cue) %>% 
  summarise(Mean_Bias = mean(mean_bias),
            Mean_Toward = mean(mean_toward),
            Mean_Away = mean(mean_away),
            Peak_Toward = mean(peak_toward),
            Peak_Away = mean(peak_away),
            Variability = mean(variability)) %>% 
  arrange(Subject, ExperimentName)

## `summarise()` has grouped output by 'Subject', 'ExperimentName'. You can
## override using the `.groups` argument.

-Illustration of AB output

Grouped_Output

5. Mixed ANOVA’s (DV = RT, WIV = Pre/Post cue, BIV = Condition)

ANOVA_df <- read_csv("/Users/noahwolkowicz/Desktop/CT/West Haven/Postdoc/Postdoc Research/Ross Attentional Bias Study/Data/AB_all.csv")

## Rows: 28 Columns: 28
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (28): ptid, Condition, SEX, MOUD, drug_BS, drug_mean_towrd, drug_mean_aw...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ANOVA_df_long <- ANOVA_df %>% 
  pivot_longer(cols = c(drug_BS:post_pain_var), names_to = "Time_Point", values_to = "RT_Values") %>%
  mutate(across(-RT_Values, factor))
glimpse(ANOVA_df)

## Rows: 28
## Columns: 28
## $ ptid                  <dbl> 100, 101, 102, 103, 104, 106, 107, 108, 109, 110…
## $ Condition             <dbl> 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, …
## $ SEX                   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 2, …
## $ MOUD                  <dbl> 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, …
## $ drug_BS               <dbl> -16.0235294, -15.4226974, 4.1907895, 51.8103617,…
## $ drug_mean_towrd       <dbl> 99.31429, 115.40000, 98.97778, 230.68966, 92.500…
## $ drug_mean_away        <dbl> 99.70588, 92.37778, 189.96970, 152.00000, 82.000…
## $ drug_peak_toward      <dbl> 234, 231, 348, 503, 220, 278, 289, 368, 468, 140…
## $ drug_peak_away        <dbl> 391, 361, 504, 539, 246, 281, 232, 229, 601, 213…
## $ drug_var              <dbl> 66.41176, 66.08696, 101.19481, 172.47619, 84.015…
## $ pain_BS               <dbl> -23.3751987, -11.3683398, -15.6986842, 47.468750…
## $ pain_mean_towrd       <dbl> 103.51724, 72.74286, 100.44681, 219.00000, 63.77…
## $ pain_mean_away        <dbl> 111.64286, 96.18919, 164.83871, 146.75758, 123.7…
## $ pain_peak_toward      <dbl> 231, 214, 279, 521, 200, 219, 323, 315, 661, 115…
## $ pain_peak_away        <dbl> 337, 336, 637, 479, 327, 218, 377, 168, 499, 231…
## $ pain_var              <dbl> 81.38571, 77.01408, 105.55844, 165.55556, 76.507…
## $ post_drug_BS          <dbl> -6.6707152, 29.8648649, 63.8740530, 1.8104167, 2…
## $ post_drug_mean_towrd  <dbl> 59.47619, 106.82609, 242.87500, 127.12903, 113.9…
## $ post_drug_mean_away   <dbl> 83.42857, 107.50000, 202.48000, 187.25806, 71.22…
## $ post_drug_peak_toward <dbl> 196, 326, 515, 413, 247, 198, 467, 297, 459, 158…
## $ post_drug_peak_away   <dbl> 181, 312, 569, 486, 240, 128, 190, 266, 560, 152…
## $ post_drug_var         <dbl> 65.85526, 88.81333, 171.40625, 144.27869, 77.900…
## $ post_pain_BS          <dbl> -5.242105, -11.419958, 11.807460, -9.224858, 26.…
## $ post_pain_mean_towrd  <dbl> 83.93750, 84.09091, 261.13793, 156.86486, 103.92…
## $ post_pain_mean_away   <dbl> 78.52174, 81.31707, 202.05882, 183.50000, 82.218…
## $ post_pain_peak_toward <dbl> 227, 173, 728, 490, 352, 209, 183, 373, 469, 112…
## $ post_pain_peak_away   <dbl> 261, 179, 507, 475, 223, 194, 261, 319, 456, 185…
## $ post_pain_var         <dbl> 70.85714, 72.18667, 202.54839, 166.23438, 88.602…

Plot function

Ross_O_Plot <- function(i) {i %>% 
  ggplot(aes(x = Time_Point, y = RT_Values, col = Condition)) +
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.25)) +
  geom_jitter(width = .05, alpha = .30) +
  xlab("Time Point") +
  scale_x_discrete(labels = c("Pre", "Post")) +
  scale_color_manual(labels = c("AB Control", "CBM"), values = c("red", "blue")) +
  ylab("Attentional Bias (ms)") +
  theme_classic() +
  theme(legend.position = "none",
        text = element_text(family = "Times New Roman",size = 12))
}

Ross_P_Plot <- function(i) {i %>% 
  ggplot(aes(x = Time_Point, y = RT_Values, col = Condition)) +
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.25)) +
  geom_jitter(width = .05, alpha = .30) +
  xlab("Time Point") +
  scale_x_discrete(labels = c("Pre", "Post")) +
  scale_color_manual(labels = c("AB Control", "CBM"), values = c("red", "blue")) +
  theme_classic() +
  theme(axis.title.y = element_blank(), 
        text = element_text(family = "Times New Roman",size = 12))
  }

5a. Drug Mean Bias

#This package/function is a wrapper for car::Anova, allowing specification of SS type
#What I've done is create separate df's containing only the within variables of interest (+ID and condition),
#and run separate ANOVAs on these
Drug_M_AOV <- ANOVA_df_long %>% filter(Time_Point %in% c("drug_BS", "post_drug_BS"))
afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Drug_M_AOV, type = 3)

## Contrasts set to contr.sum for the following variables: Condition

## Anova Table (Type 3 tests)
## 
## Response: RT_Values
##                 Effect    df     MSE      F  ges p.value
## 1            Condition 1, 26  918.40   1.38 .024    .251
## 2           Time_Point 1, 26 1040.86 3.44 + .066    .075
## 3 Condition:Time_Point 1, 26 1040.86   0.46 .009    .504
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

glimpse(ANOVA_df_long)

## Rows: 672
## Columns: 6
## $ ptid       <fct> 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100,…
## $ Condition  <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ SEX        <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ MOUD       <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ Time_Point <fct> drug_BS, drug_mean_towrd, drug_mean_away, drug_peak_toward,…
## $ RT_Values  <dbl> -16.023529, 99.314286, 99.705882, 234.000000, 391.000000, 6…

emmeans::emmeans(afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Drug_M_AOV, type = 3, return = afex_options("return_aov")), specs = "Time_Point")

## Contrasts set to contr.sum for the following variables: Condition

##  Time_Point   emmean   SE df lower.CL upper.CL
##  drug_BS       -6.87 5.38 26   -17.92     4.19
##  post_drug_BS   9.17 6.43 26    -4.06    22.39
## 
## Results are averaged over the levels of: Condition 
## Confidence level used: 0.95

Drug_M_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Opioid_MB_Plot <- Ross_O_Plot(Drug_M_AOV)

5b. Pain Mean Bias

Pain_M_AOV <- ANOVA_df_long %>% filter(Time_Point %in% c("pain_BS", "post_pain_BS"))
afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Pain_M_AOV, type = 3)

## Contrasts set to contr.sum for the following variables: Condition

## Anova Table (Type 3 tests)
## 
## Response: RT_Values
##                 Effect    df     MSE    F   ges p.value
## 1            Condition 1, 26  593.35 0.34  .005    .564
## 2           Time_Point 1, 26 1009.01 0.01 <.001    .906
## 3 Condition:Time_Point 1, 26 1009.01 0.20  .005    .659
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Pain_M_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Pain_MB_Plot <- Ross_P_Plot(Pain_M_AOV)

5c. Drug Mean Toward

Drug_M_Twrd_AOV <- ANOVA_df_long %>% filter(Time_Point %in% c("drug_mean_towrd", "post_drug_mean_towrd"))
afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Drug_M_Twrd_AOV, type = 3)

## Contrasts set to contr.sum for the following variables: Condition

## Anova Table (Type 3 tests)
## 
## Response: RT_Values
##                 Effect    df     MSE       F  ges p.value
## 1            Condition 1, 26 4123.24 8.56 ** .191    .007
## 2           Time_Point 1, 26 1609.41    0.09 .001    .762
## 3 Condition:Time_Point 1, 26 1609.41    1.18 .013    .286
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Drug_M_Twrd_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Drug_M_Twrd_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Drug_M_Twrd_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

Opioid_MT_Plot <- Ross_O_Plot(Drug_M_Twrd_AOV)

5d. Pain Mean Toward

Pain_M_twrd_AOV <- ANOVA_df_long %>% filter(Time_Point %in% c("pain_mean_towrd", "post_pain_mean_towrd"))
afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Pain_M_twrd_AOV, type = 3)

## Contrasts set to contr.sum for the following variables: Condition

## Anova Table (Type 3 tests)
## 
## Response: RT_Values
##                 Effect    df     MSE      F  ges p.value
## 1            Condition 1, 26 5364.99 5.76 * .147    .024
## 2           Time_Point 1, 26 1535.92   0.74 .006    .398
## 3 Condition:Time_Point 1, 26 1535.92 6.42 * .052    .018
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Pain_M_twrd_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

Pain_M_twrd_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Pain_MT_Plot <- Ross_P_Plot(Pain_M_twrd_AOV)

5e. Drug Mean Away

Drug_M_Away_AOV <- ANOVA_df_long %>% filter(Time_Point %in% c("drug_mean_away", "post_drug_mean_away"))
afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Drug_M_Away_AOV, type = 3)

## Contrasts set to contr.sum for the following variables: Condition

## Anova Table (Type 3 tests)
## 
## Response: RT_Values
##                 Effect    df     MSE      F  ges p.value
## 1            Condition 1, 26 5178.90 4.99 * .143    .034
## 2           Time_Point 1, 26  773.49   0.26 .001    .614
## 3 Condition:Time_Point 1, 26  773.49   2.28 .011    .143
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Drug_M_Away_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Drug_M_Away_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

Opioid_MA_Plot <- Ross_O_Plot(Drug_M_Away_AOV)

5f. Pain Mean Away

Pain_M_Away_AOV <- ANOVA_df_long %>% filter(Time_Point %in% c("pain_mean_away", "post_pain_mean_away"))
afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Pain_M_Away_AOV, type = 3)

## Contrasts set to contr.sum for the following variables: Condition

## Anova Table (Type 3 tests)
## 
## Response: RT_Values
##                 Effect    df     MSE      F  ges p.value
## 1            Condition 1, 26 4457.38 6.11 * .148    .020
## 2           Time_Point 1, 26 1576.59   0.13 .001    .724
## 3 Condition:Time_Point 1, 26 1576.59 3.33 + .032    .079
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Pain_M_Away_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Pain_M_Away_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

Pain_MA_Plot <- Ross_P_Plot(Pain_M_Away_AOV)

5g. Drug Peak Toward

Drug_P_Toward_AOV <- ANOVA_df_long %>% filter(Time_Point %in% c("drug_peak_toward", "post_drug_peak_toward"))
afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Drug_P_Toward_AOV, type = 3)

## Contrasts set to contr.sum for the following variables: Condition

## Anova Table (Type 3 tests)
## 
## Response: RT_Values
##                 Effect    df      MSE      F   ges p.value
## 1            Condition 1, 26 35044.60 5.16 *  .144    .032
## 2           Time_Point 1, 26  6478.72   0.94  .006    .341
## 3 Condition:Time_Point 1, 26  6478.72   0.07 <.001    .793
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Drug_P_Toward_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Drug_P_Toward_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

Opioid_PT_Plot <- Ross_O_Plot(Drug_P_Toward_AOV)

5h. Peak Pain Toward

Pain_P_twrd_AOV <- ANOVA_df_long %>% filter(Time_Point %in% c("pain_peak_toward", "post_pain_peak_toward"))
afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Pain_P_twrd_AOV, type = 3)

## Contrasts set to contr.sum for the following variables: Condition

## Anova Table (Type 3 tests)
## 
## Response: RT_Values
##                 Effect    df      MSE      F  ges p.value
## 1            Condition 1, 26 50428.84 5.32 * .150    .029
## 2           Time_Point 1, 26  8219.99   0.26 .001    .618
## 3 Condition:Time_Point 1, 26  8219.99   1.41 .008    .245
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Pain_P_twrd_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Pain_P_twrd_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

Pain_PT_Plot <- Ross_P_Plot(Pain_P_twrd_AOV)

5i. Drug Peak Away

Drug_P_Away_AOV <- ANOVA_df_long %>% filter(Time_Point %in% c("drug_peak_away", "post_drug_peak_away"))
afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Drug_P_Away_AOV, type = 3)

## Contrasts set to contr.sum for the following variables: Condition

## Anova Table (Type 3 tests)
## 
## Response: RT_Values
##                 Effect    df      MSE      F   ges p.value
## 1            Condition 1, 26 57895.85 5.97 *  .169    .022
## 2           Time_Point 1, 26  7305.64   0.36  .002    .553
## 3 Condition:Time_Point 1, 26  7305.64   0.00 <.001    .953
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Drug_P_Away_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Drug_P_Away_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

Opioid_PA_Plot <- Ross_O_Plot(Drug_P_Away_AOV)

5j. Peak Pain Away

Pain_P_Away_AOV <- ANOVA_df_long %>% filter(Time_Point %in% c("pain_peak_away", "post_pain_peak_away"))
afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Pain_P_Away_AOV, type = 3)

## Contrasts set to contr.sum for the following variables: Condition

## Anova Table (Type 3 tests)
## 
## Response: RT_Values
##                 Effect    df      MSE      F  ges p.value
## 1            Condition 1, 26 51254.93 4.24 * .128    .050
## 2           Time_Point 1, 26  5918.46   0.95 .004    .338
## 3 Condition:Time_Point 1, 26  5918.46   0.97 .004    .333
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Pain_P_Away_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Pain_P_Away_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

Pain_PA_Plot <- Ross_P_Plot(Pain_P_Away_AOV)

5k. Drug Variability

Drug_V_AOV <- ANOVA_df_long %>% filter(Time_Point %in% c("drug_var", "post_drug_var"))
afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Drug_V_AOV, type = 3)

## Contrasts set to contr.sum for the following variables: Condition

## Anova Table (Type 3 tests)
## 
## Response: RT_Values
##                 Effect    df     MSE      F  ges p.value
## 1            Condition 1, 26 2823.56 7.62 * .199    .010
## 2           Time_Point 1, 26  510.11   0.66 .004    .425
## 3 Condition:Time_Point 1, 26  510.11 3.06 + .018    .092
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Drug_V_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Drug_V_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

Opioid_V_Plot <- Ross_O_Plot(Drug_V_AOV)

Opioid_V_Plot <- Opioid_V_Plot + scale_y_continuous(limits = c(30, 350), breaks = seq(0, 350, by = 100))

5l. Pain Variability

Pain_V_AOV <- ANOVA_df_long %>% filter(Time_Point %in% c("pain_var", "post_pain_var"))
afex::aov_car(formula = RT_Values ~ Time_Point*Condition + Error(ptid/(Time_Point)), data = Pain_V_AOV, type = 3)

## Contrasts set to contr.sum for the following variables: Condition

## Anova Table (Type 3 tests)
## 
## Response: RT_Values
##                 Effect    df     MSE      F  ges p.value
## 1            Condition 1, 26 3911.63 6.03 * .154    .021
## 2           Time_Point 1, 26 1085.28   0.77 .006    .388
## 3 Condition:Time_Point 1, 26 1085.28 5.24 * .042    .030
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Pain_V_AOV %>% select(Condition, Time_Point, RT_Values) %>% group_by(Condition, Time_Point) %>% summarise(M = mean(RT_Values), Med = median(RT_Values), SD = sd(RT_Values))

## `summarise()` has grouped output by 'Condition'. You can override using the
## `.groups` argument.

Pain_V_Plot <- Ross_P_Plot(Pain_V_AOV)

Big Ol’ Plots

library(cowplot)
setwd("/Users/noahwolkowicz/Desktop/CT/West Haven/Postdoc/Postdoc Research/Ross Attentional Bias Study/Data/Ross Data/Figures/")

cowplot::plot_grid(Opioid_MB_Plot, Pain_MB_Plot, 
                   Opioid_MT_Plot, Pain_MT_Plot,
                   Opioid_MA_Plot, Pain_MA_Plot,
                   nrow = 3, ncol = 2)

x <- cowplot::plot_grid(Opioid_MB_Plot, Pain_MB_Plot, 
                   Opioid_MT_Plot, Pain_MT_Plot,
                   Opioid_MA_Plot, Pain_MA_Plot,
                   nrow = 3, ncol = 2)

MT_Plot <- cowplot::plot_grid(Opioid_MT_Plot, Pain_MT_Plot,
                   nrow = 1, ncol = 2, rel_widths = c(1,1.45)) 
TV_Plot <- cowplot::plot_grid(Opioid_V_Plot, Pain_V_Plot,
                   nrow = 1, ncol = 2, rel_widths = c(1,1.45)) 

#save_plot("/Users/noahwolkowicz/Desktop/CT/West Haven/Postdoc/Postdoc Research/Ross Attentional Bias Study/Data/Ross Data/Figures/MT_Plot.png", MT_Plot)
#save_plot("/Users/noahwolkowicz/Desktop/CT/West Haven/Postdoc/Postdoc Research/Ross Attentional Bias Study/Data/Ross Data/Figures/TV_Plot.png", TV_Plot)