Targets and Standards, Auditory Oddball MIST

Infinity Reference or Reference Electrode Standardization Technique (REST)

cat("\014")     # clean terminal



rm(list = ls()) # clean workspace
try(dev.off(), silent = TRUE) # close all plots
library(afex)
library(emmeans)
library(ggplot2)
library(ggridges)
library(ggdist)
library(dplyr)
library(reshape2)
library(GGally)
library(forcats)

theme_set(
  theme_minimal()
)
a_posteriori <- function(afex_aov, sig_level = .05) {
  factors  <- as.list(rownames(afex_aov$anova_table))
  for (j in 1:length(factors)) {
    if (grepl(":", factors[[j]])) {
      factors[[j]] <- unlist(strsplit(factors[[j]], ":"))
    }
  }
  p_values <- afex_aov$anova_table$`Pr(>F)`
  for (i in 1:length(p_values)) {
    if (p_values[i] <= sig_level) {
      print(emmeans(afex_aov, factors[[i]], contr = "pairwise"))
    }
  }
}

group_id   <- read.table("~/Insync/Drive/00EEG/Proyectos/Huepe/fdcyt_2017/Registro-Evaluaciones-FDCYT-DH-2017 - General ANONIMO.csv",
                         sep = ",", header = TRUE, col.names = c("full.id", "Subject", "Sex", "Group", "Stress"))
stress_response_name  <- paste('~/Insync/Drive/00EEG/Proyectos/Huepe/fdcyt_2017/gabormeta_huepe/behavioural/', 'stress_response.csv', sep='/')
stress_response_data  <- read.csv(stress_response_name, header = TRUE)
master_dir <- '~/Insync/Drive/00EEG/Proyectos/Huepe/fdcyt_2017/mist_oddball'
data_dir   <- paste(master_dir, 'results',  sep = '/')
target_and_standard_name <- paste(data_dir, 'rest_350to550_mist_oddball.txt', sep='/')
target_and_standard_data <- read.table(target_and_standard_name, header = TRUE, strip.white = TRUE, sep = "\t")
target_and_standard_data$Subject  <- readr::parse_number(target_and_standard_data$ERPset)
target_and_standard_data <- merge(target_and_standard_data, group_id, by = "Subject")
target_and_standard_data <- merge(target_and_standard_data, stress_response_data, all.x = TRUE)
names(target_and_standard_data)[names(target_and_standard_data) == "value"]    <- "uvolts"
names(target_and_standard_data)[names(target_and_standard_data) == "binlabel"] <- "stimulus"
names(target_and_standard_data)[names(target_and_standard_data) == "chlabel"]  <- "channel"
target_and_standard_data$session  <- readr::parse_number(target_and_standard_data$stimulus)
target_and_standard_data$stimulus <- gsub("\\_.*", "", target_and_standard_data$stimulus)
target_and_standard_data$Sex      <- factor(target_and_standard_data$Sex)
target_and_standard_data$Group    <- factor(target_and_standard_data$Group)
target_and_standard_data$Stress   <- factor(target_and_standard_data$Stress)
target_and_standard_data$session  <- factor(target_and_standard_data$session)
target_and_standard_data$channel  <- factor(target_and_standard_data$channel)
levels(target_and_standard_data$Sex)     <- list(female  = "F", male  = "M")
levels(target_and_standard_data$Group)   <- list(non_vulnerable  = "CN", vulnerable  = "EX")
levels(target_and_standard_data$Stress)  <- list(not_stressed  = "NS", stressed  = "SS")
levels(target_and_standard_data$session) <- list(first  = "1", second  = "2")
levels(target_and_standard_data$channel) <- list(P1  = "E005-P1a", Pz  = "E019-Pz", P2 = "E032-P2a", PO3 = "E017-PO3a", POz = "E021-POz", PO4 = "E030-PO4a")
target_and_standard_data$anteroposterior[target_and_standard_data$channel %in% c('P1' , 'Pz' , 'P2')]  <- 'Parietal'
target_and_standard_data$anteroposterior[target_and_standard_data$channel %in% c('PO3', 'POz', 'PO4')] <- 'Parieto-occipital'
target_and_standard_data$mediolateral[target_and_standard_data$channel %in% c('P1', 'PO3')]            <- 'Left'
target_and_standard_data$mediolateral[target_and_standard_data$channel %in% c('Pz', 'POz')]            <- 'Midline'
target_and_standard_data$mediolateral[target_and_standard_data$channel %in% c('P2', 'PO4')]            <- 'Right'
target_and_standard_data$Subject         <- factor(target_and_standard_data$Subject)
target_and_standard_data$Group           <- factor(target_and_standard_data$Group)
target_and_standard_data$stimulus        <- factor(target_and_standard_data$stimulus, levels = c('target', 'standard'))
target_and_standard_data$anteroposterior <- factor(target_and_standard_data$anteroposterior)
target_and_standard_data$mediolateral    <- factor(target_and_standard_data$mediolateral, levels = c('Left', 'Midline', 'Right'))
target_and_standard_data$stress_response <- factor(target_and_standard_data$stress_response)
target_and_standard_data$eliminate       <- factor(target_and_standard_data$eliminate)
write.csv(target_and_standard_data,  paste(data_dir, '/targets_and_standards_data_clean.csv', sep = ''),  row.names = FALSE)

1 General Description

options(width = 100)
mytable1 <- xtabs(~ Stress + Group, data = target_and_standard_data) / length(unique(target_and_standard_data$channel)) / length(unique(target_and_standard_data$stimulus))  / length(unique(target_and_standard_data$session))
ftable(addmargins(mytable1))

             Group non_vulnerable vulnerable Sum
Stress                                          
not_stressed                   12         15  27
stressed                       13         15  28
Sum                            25         30  55

summary(target_and_standard_data[c('uvolts', 'Sex', 'Group', 'Stress', 'stress_response', 'eliminate', 'stimulus', 'channel', 'anteroposterior', 'mediolateral', 'Subject')])

     uvolts             Sex                 Group              Stress    stress_response
 Min.   :-5.44060   female:672   non_vulnerable:600   not_stressed:648   no  :576       
 1st Qu.: 0.03013   male  :648   vulnerable    :720   stressed    :672   yes :552       
 Median : 1.39330                                                        NA's:192       
 Mean   : 2.27034                                                                       
 3rd Qu.: 3.92770                                                                       
 Max.   :14.10599                                                                       
                                                                                        
 eliminate       stimulus   channel            anteroposterior  mediolateral    Subject    
 no  :1032   target  :660   P1 :220   Parietal         :660    Left   :440   1      :  24  
 yes :  96   standard:660   Pz :220   Parieto-occipital:660    Midline:440   2      :  24  
 NA's: 192                  P2 :220                            Right  :440   4      :  24  
                            PO3:220                                          5      :  24  
                            POz:220                                          6      :  24  
                            PO4:220                                          7      :  24  
                                                                             (Other):1176  

2 ERP plots

2.1 Targets and Standards, Central-Occipital ROI

There were 2 Oddball Task sessions, before and after the Montreal Imaging Stress Task.

2.2 Full topographic map:

3 Targets and Standards by topography

Measurement window: [350.0 550.0] ms

options(width = 100)
target_and_standard_rain_central <- ggplot(target_and_standard_data[target_and_standard_data$anteroposterior == "Parietal", ], aes(y = uvolts, x = mediolateral, color = stimulus, fill = stimulus)) +
  ggtitle("Targets and Standards Parietal") +
  ylab("uvolts") +
  stat_halfeye(
    trim   = FALSE, 
    adjust = .75, 
    .width = 0, 
    justification = -.15, 
    alpha  = .4,
    point_colour = NA) + 
  # theme(legend.position='none')
  geom_boxplot(width = .15, alpha = .2, outlier.shape = NA) +
  geom_point(size = 2, alpha = .4, position = position_jitter(width = .05, height = 0)) 
suppressWarnings(print(target_and_standard_rain_central))

target_and_standard_rain_parietal <- ggplot(target_and_standard_data[target_and_standard_data$anteroposterior == "Parieto-occipital", ], aes(y = uvolts, x = mediolateral, color = stimulus, fill = stimulus)) +
  ggtitle("Targets and Standards Parieto-occipital") +
  ylab("uvolts") +
  stat_halfeye(
    trim   = FALSE, 
    adjust = .75, 
    .width = 0, 
    justification = -.15, 
    alpha  = .4,
    point_colour = NA) + 
  # theme(legend.position='none')
  geom_boxplot(width = .15, alpha = .2, outlier.shape = NA) +
  geom_point(size = 2, alpha = .4, position = position_jitter(width = .05, height = 0)) 
suppressWarnings(print(target_and_standard_rain_parietal))

target_and_standard_rep_anova = aov_ez("Subject", "uvolts", target_and_standard_data, within = c("stimulus", "anteroposterior", "mediolateral", "session"))
target_and_standard_afex_plot <-
  afex_plot(
    target_and_standard_rep_anova,
    x = "anteroposterior",
    trace = "stimulus",
    panel = "session",
    error = "within",
    error_arg = list(width = .25),
    dodge = -.5,
    mapping = c("color"),
    point_arg = list(size = 4)
  )
suppressWarnings(print(target_and_standard_afex_plot))

nice(target_and_standard_rep_anova)

Anova Table (Type 3 tests)

Response: uvolts
                                          Effect           df   MSE          F   ges p.value
1                                       stimulus        1, 54 32.38 189.36 ***  .473   <.001
2                                anteroposterior        1, 54  3.92  13.69 ***  .008   <.001
3                                   mediolateral 1.87, 101.19  1.22       2.01 <.001    .143
4                                        session        1, 54 12.34  21.70 ***  .038   <.001
5                       stimulus:anteroposterior        1, 54  4.18     4.08 *  .002    .048
6                          stimulus:mediolateral  1.68, 90.86  0.99       2.43 <.001    .102
7                   anteroposterior:mediolateral  1.69, 91.19  0.51     4.44 * <.001    .019
8                               stimulus:session        1, 54 11.75     5.98 *  .010    .018
9                        anteroposterior:session        1, 54  0.67       0.04 <.001    .835
10                          mediolateral:session 1.95, 105.42  0.24       0.15 <.001    .855
11         stimulus:anteroposterior:mediolateral  1.55, 83.74  0.71     4.62 * <.001    .020
12              stimulus:anteroposterior:session        1, 54  0.57     5.56 * <.001    .022
13                 stimulus:mediolateral:session  1.38, 74.48  0.42       1.46 <.001    .237
14          anteroposterior:mediolateral:session  1.48, 80.06  0.24       1.93 <.001    .162
15 stimulus:anteroposterior:mediolateral:session  1.34, 72.32  0.27       0.02 <.001    .945
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

Sphericity correction method: GG 

emmeans(target_and_standard_rep_anova, "stimulus", contr = "pairwise")

$emmeans
 stimulus emmean    SE df lower.CL upper.CL
 target    4.426 0.334 54    3.757    5.094
 standard  0.115 0.135 54   -0.156    0.387

Results are averaged over the levels of: session, mediolateral, anteroposterior 
Confidence level used: 0.95 

$contrasts
 contrast          estimate    SE df t.ratio p.value
 target - standard     4.31 0.313 54  13.761  <.0001

Results are averaged over the levels of: session, mediolateral, anteroposterior 

emmeans(target_and_standard_rep_anova, "anteroposterior", contr = "pairwise")

$emmeans
 anteroposterior   emmean    SE df lower.CL upper.CL
 Parietal            2.47 0.203 54     2.06     2.88
 Parieto.occipital   2.07 0.213 54     1.64     2.50

Results are averaged over the levels of: session, mediolateral, stimulus 
Confidence level used: 0.95 

$contrasts
 contrast                     estimate    SE df t.ratio p.value
 Parietal - Parieto.occipital    0.403 0.109 54   3.700  0.0005

Results are averaged over the levels of: session, mediolateral, stimulus 

emmeans(target_and_standard_rep_anova, "session", contr = "pairwise")

$emmeans
 session emmean    SE df lower.CL upper.CL
 first     2.72 0.223 54     2.27     3.17
 second    1.82 0.223 54     1.37     2.27

Results are averaged over the levels of: mediolateral, anteroposterior, stimulus 
Confidence level used: 0.95 

$contrasts
 contrast       estimate    SE df t.ratio p.value
 first - second    0.901 0.193 54   4.658  <.0001

Results are averaged over the levels of: mediolateral, anteroposterior, stimulus 

p300_data_name    <- paste(data_dir, 'rest_350to550_mist_oddball_difference.txt', sep='/')
p300_data         <- read.table(p300_data_name, header = TRUE, strip.white = TRUE, sep = "\t")
p300_data$Subject <- readr::parse_number(p300_data$ERPset)
p300_data         <- merge(p300_data, group_id, by = "Subject")
p300_data         <- merge(p300_data, stress_response_data, all.x = TRUE)
names(p300_data)[names(p300_data) == "value"]    <- "uvolts"
names(p300_data)[names(p300_data) == "binlabel"] <- "session"
names(p300_data)[names(p300_data) == "chlabel"]  <- "channel"
p300_data$session  <- readr::parse_number(p300_data$session)
p300_data$num_id   <- readr::parse_number(p300_data$ERPset)
p300_data$Sex      <- factor(p300_data$Sex)
p300_data$Group    <- factor(p300_data$Group)
p300_data$Stress   <- factor(p300_data$Stress)
p300_data$session  <- factor(p300_data$session)
p300_data$channel  <- factor(p300_data$channel)
levels(p300_data$Sex)     <- list(female  = "F", male  = "M")
levels(p300_data$Group)   <- list(non_vulnerable  = "CN", vulnerable  = "EX")
levels(p300_data$Stress)  <- list(not_stressed  = "NS", stressed  = "SS")
levels(p300_data$session) <- list(first  = "1", second  = "2")
levels(p300_data$channel) <- list(P1  = "E005-P1a", Pz  = "E019-Pz", P2 = "E032-P2a", PO3 = "E017-PO3a", POz = "E021-POz", PO4 = "E030-PO4a")
p300_data$anteroposterior[p300_data$channel %in% c('P1' , 'Pz' , 'P2')]  <- 'Parietal'
p300_data$anteroposterior[p300_data$channel %in% c('PO3', 'POz', 'PO4')] <- 'Parieto-occipital'
p300_data$mediolateral[p300_data$channel %in% c('P1', 'PO3')]            <- 'Left'
p300_data$mediolateral[p300_data$channel %in% c('Pz', 'POz')]            <- 'Midline'
p300_data$mediolateral[p300_data$channel %in% c('P2', 'PO4')]            <- 'Right'
p300_data$Subject         <- factor(p300_data$Subject)
p300_data$Group           <- factor(p300_data$Group)
p300_data$anteroposterior <- factor(p300_data$anteroposterior)
p300_data$mediolateral    <- factor(p300_data$mediolateral, levels = c('Left', 'Midline', 'Right'))
p300_data$stress_response <- factor(p300_data$stress_response)
p300_data$eliminate       <- factor(p300_data$eliminate)
write.csv(p300_data,  paste(data_dir, '/auditory_oddball_data_clean.csv', sep = ''),  row.names = FALSE)

4 General P3 Description

options(width = 100)
mytable1 <- xtabs(~ Stress + Group, data = p300_data) / length(unique(p300_data$channel))  / length(unique(p300_data$session))
ftable(addmargins(mytable1))

             Group non_vulnerable vulnerable Sum
Stress                                          
not_stressed                   12         15  27
stressed                       13         15  28
Sum                            25         30  55

without_cortisol  <- unique(p300_data$full.id[is.na(p300_data$stress_response)])
noquote("Subjects without cortisol info:")

[1] Subjects without cortisol info:

cat(paste(without_cortisol[order(sub('.*_', '', without_cortisol))], collapse = '\n'))

EX_NS_FE_001
EX_SS_MA_002
EX_SS_FE_004
EX_NS_MA_005
EX_NS_FE_006
EX_SS_FE_007
EX_NS_FE_008
EX_SS_FE_009

noquote("Cortisol outliers:")

[1] Cortisol outliers:

cat(paste(unique(na.omit(p300_data$full.id[p300_data$eliminate == 'yes'])), collapse = '\n'))

CN_NS_FE_037
EX_NS_FE_015
EX_NS_MA_028
EX_SS_FE_024

mytable2 <- xtabs(~ stress_response + Group, data = p300_data[p300_data$eliminate == 'no', ]) / length(unique(p300_data$channel))  / length(unique(p300_data$session))
ftable(addmargins(mytable2))

                Group non_vulnerable vulnerable Sum
stress_response                                    
no                                14          9  23
yes                               10         10  20
Sum                               24         19  43

summary(p300_data[c('uvolts', 'Sex', 'Group', 'Stress', 'stress_response', 'eliminate', 'channel', 'anteroposterior', 'mediolateral', 'Subject')])

     uvolts           Sex                 Group              Stress    stress_response eliminate 
 Min.   :-3.312   female:336   non_vulnerable:300   not_stressed:324   no  :288        no  :516  
 1st Qu.: 2.234   male  :324   vulnerable    :360   stressed    :336   yes :276        yes : 48  
 Median : 3.922                                                        NA's: 96        NA's: 96  
 Mean   : 4.310                                                                                  
 3rd Qu.: 6.035                                                                                  
 Max.   :15.993                                                                                  
                                                                                                 
 channel            anteroposterior  mediolateral    Subject   
 P1 :110   Parietal         :330    Left   :220   1      : 12  
 Pz :110   Parieto-occipital:330    Midline:220   2      : 12  
 P2 :110                            Right  :220   4      : 12  
 PO3:110                                          5      : 12  
 POz:110                                          6      : 12  
 PO4:110                                          7      : 12  
                                                  (Other):588  

5 ERP plots, P3 = Target - Standard

5.1 P3 Parietal-Occipital ROI

There were 2 Oddball Task sessions, before and after the Montreal Imaging Stress Task.

5.2 Full topographic P3 map:

5.3 Scalp P3 Map:

6 P3 by topography (no outliers)

Measurement window: [350.0 550.0] ms

options(width = 100)
p300_data <- p300_data[!is.na(p300_data$stress_response), ]
p300_data <- p300_data[p300_data$eliminate == "no", ]
p300_rep_anova = aov_ez("Subject", "uvolts", p300_data, within = c("anteroposterior", "mediolateral", "session"))
p300_rain <- ggplot(p300_rep_anova$data$long, aes(y = uvolts, x = session, color = anteroposterior, fill = anteroposterior)) +
  ggtitle("P300") +
  ylab("uvolts") +
  stat_halfeye(
    trim   = FALSE, 
    adjust = .75, 
    .width = 0, 
    justification = -.15, 
    alpha  = .4,
    point_colour = NA) + 
  # theme(legend.position='none')
  # geom_boxplot(width = .15, alpha = .2, outlier.shape = NA) +
  geom_point(size = 2, alpha = .4, position = position_jitter(width = .05, height = 0)) 
suppressWarnings(print(p300_rain))

p300_afex_plot <-
  afex_plot(
    p300_rep_anova,
    x = "session",
    trace = "anteroposterior",
    panel = "mediolateral",
    error = "within",
    error_arg = list(width = .1),
    dodge = -.5,
    mapping = c("color"),
    point_arg = list(size = 4)
  )
suppressWarnings(print(p300_afex_plot))

nice(p300_rep_anova)

Anova Table (Type 3 tests)

Response: uvolts
                                Effect          df   MSE      F   ges p.value
1                      anteroposterior       1, 42  7.85 3.98 +  .007    .053
2                         mediolateral 1.75, 73.43  1.85   2.34  .002    .111
3                              session       1, 42 24.01 7.03 *  .036    .011
4         anteroposterior:mediolateral 1.48, 62.33  1.32 2.73 +  .001    .088
5              anteroposterior:session       1, 42  1.35 4.18 *  .001    .047
6                 mediolateral:session 1.43, 59.91  0.80   1.09 <.001    .325
7 anteroposterior:mediolateral:session 1.35, 56.59  0.58   0.02 <.001    .939
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

Sphericity correction method: GG 

# emmeans(p300_rep_anova, "anteroposterior", contr = "pairwise")
emmeans(p300_rep_anova, "session", contr = "pairwise")

$emmeans
 session emmean    SE df lower.CL upper.CL
 first     4.82 0.489 42     3.83     5.80
 second    3.67 0.330 42     3.01     4.34

Results are averaged over the levels of: mediolateral, anteroposterior 
Confidence level used: 0.95 

$contrasts
 contrast       estimate    SE df t.ratio p.value
 first - second     1.14 0.431 42   2.650  0.0113

Results are averaged over the levels of: mediolateral, anteroposterior 

7 P3 by group, average ROI (no outliers)

options(width = 100)
p300_anova_groups = aov_ez("Subject", "uvolts", p300_data, between = c("stress_response", "Group"), within = c("session"), fun_aggregate = mean)

Contrasts set to contr.sum for the following variables: stress_response, Group

p300_rain_groups <- ggplot(p300_anova_groups$data$long, aes(y = uvolts, x = session, color = Group, fill = Group)) +
  ggtitle("P300") +
  ylab("uvolts") +
  stat_halfeye(
    trim   = FALSE, 
    adjust = .75, 
    .width = 0, 
    justification = -.15, 
    alpha  = .4,
    point_colour = NA) + 
  # theme(legend.position='none')
  # geom_boxplot(width = .15, alpha = .2, outlier.shape = NA) +
  geom_point(size = 2, alpha = .4, position = position_jitter(width = .05, height = 0)) 
suppressWarnings(print(p300_rain_groups))

p300_afex_plot_goups <-
  afex_plot(
    p300_anova_groups,
    x     = "Group",
    trace = "stress_response",
    panel = "session",
    error = "between",
    error_arg = list(width = .1),
    dodge = -.5,
    mapping = c("color"),
    point_arg = list(size = 4)
  )
suppressWarnings(print(p300_afex_plot_goups))

nice(p300_anova_groups)

Anova Table (Type 3 tests)

Response: uvolts
                         Effect    df   MSE      F  ges p.value
1               stress_response 1, 39 11.27   0.71 .013    .405
2                         Group 1, 39 11.27   0.43 .008    .518
3         stress_response:Group 1, 39 11.27   0.36 .007    .552
4                       session 1, 39  3.90 7.28 * .046    .010
5       stress_response:session 1, 39  3.90   1.30 .009    .261
6                 Group:session 1, 39  3.90   1.79 .012    .189
7 stress_response:Group:session 1, 39  3.90   1.18 .008    .284
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

emmeans(p300_anova_groups, "session", contr = "pairwise")

$emmeans
 session emmean    SE df lower.CL upper.CL
 first     4.75 0.507 39     3.73     5.78
 second    3.59 0.325 39     2.93     4.24

Results are averaged over the levels of: stress_response, Group 
Confidence level used: 0.95 

$contrasts
 contrast       estimate    SE df t.ratio p.value
 first - second     1.17 0.432 39   2.699  0.0102

Results are averaged over the levels of: stress_response, Group