Auditory Oddball P3, PeMyCreP

Reference at Infinity, REST (Reference Electrode Standardization Technique).

cat("\014")     # clean terminal



rm(list = ls()) # clean workspace
try(dev.off(), silent = TRUE) # close all plots
library(tidyverse)
library(afex)
library(emmeans)
library(ggdist)
library(GGally)
library(readxl)
library(data.table)
library(performance)

exclude_bad_eeg <- TRUE
my_dodge <- .3
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"))
      cat(rep("_", 100), '\n', sep = "")
    }
  }
}

eeg_check <- read_excel(file.path('..', 'bad channels oddball pemycrep 2022.xlsx'))
eeg_check <- eeg_check %>%
  mutate(badchan_num = ifelse(badchan == '0', 0, sapply(strsplit(badchan, " "), length)))
bad_eeg   <- eeg_check$name[eeg_check$commentary != 'ok']
data_dir  <- file.path('..', 'results')
# target_and_standard_name <- file.path(data_dir, 'average_voltage_275_to_425_auditory_oddball_standard_and_target.txt')
target_and_standard_name <- file.path(data_dir, 'average_voltage_300_to_400_auditory_oddball_standard_and_target.txt')
target_and_standard_data <- read.table(target_and_standard_name, header = TRUE, strip.white = TRUE, sep = "\t")
names(target_and_standard_data)[names(target_and_standard_data) == "value"] <- "uvolts"
names(target_and_standard_data)[names(target_and_standard_data) == "binlabel"] <- "stimulus"
target_and_standard_data$num_id <- readr::parse_number(target_and_standard_data$ERPset)
target_and_standard_data$vulnerability[ grepl("nVul", target_and_standard_data$ERPset)] <- "Invulnerable"
target_and_standard_data$vulnerability[!grepl("nVul", target_and_standard_data$ERPset)] <- "Vulnerable"
target_and_standard_data$belief[ grepl("nCr", target_and_standard_data$ERPset)]         <- "Unbeliever"
target_and_standard_data$belief[!grepl("nCr", target_and_standard_data$ERPset)]         <- "Believer"
target_and_standard_data$sex[ grepl("F", target_and_standard_data$ERPset)]              <- "Female"
target_and_standard_data$sex[!grepl("F", target_and_standard_data$ERPset)]              <- "Male"
target_and_standard_data$anteroposterior[target_and_standard_data$chlabel %in% c('C1', 'Cz', 'C2')]                <- 'Central'
target_and_standard_data$anteroposterior[target_and_standard_data$chlabel %in% c('E112-CP1a', 'CPz', 'E034-CP2a')] <- 'Centro-parietal'
target_and_standard_data$mediolateral[target_and_standard_data$chlabel %in% c('C1', 'E112-CP1a')]                  <- 'Left'
target_and_standard_data$mediolateral[target_and_standard_data$chlabel %in% c('Cz', 'CPz')]                        <- 'Midline'
target_and_standard_data$mediolateral[target_and_standard_data$chlabel %in% c('C2', 'E034-CP2a')]                  <- 'Right'
target_and_standard_data$num_id          <- factor(target_and_standard_data$num_id)
target_and_standard_data$vulnerability   <- factor(target_and_standard_data$vulnerability)
target_and_standard_data$sex             <- factor(target_and_standard_data$sex)
target_and_standard_data$belief          <- factor(target_and_standard_data$belief)
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'))
just_sex <- subset(target_and_standard_data[c('ERPset', 'sex', 'chindex', 'stimulus')], chindex == 1 & stimulus == 'Target')
connectivity_file <- 'erplab/median_connectivity_auditory_oddball.csv'
connectivity_data <- read_csv(connectivity_file, show_col_types = FALSE) |> 
  left_join(y = just_sex[c('ERPset', 'sex')], by = c('ERPset')) |>
  mutate(log_median_connectivity = log(median_connectivity)) |> 
  mutate(stimulus = factor(stimulus, levels = c('target', 'standard'))) |> 
  mutate(num_id = parse_number(ERPset)) |> 
  mutate_if(is.character, as.factor) |> 
  mutate(region_b = factor(region_b, levels = c('cingulate_ipsi', 'cingulate_contra')))
if (exclude_bad_eeg) {
  target_and_standard_data <- target_and_standard_data[!(target_and_standard_data$ERPset %in% bad_eeg), ]
  connectivity_data        <- connectivity_data[!(connectivity_data$ERPset %in% bad_eeg), ]
  }
write.csv(target_and_standard_data,  file.path(data_dir, 'targets_and_standards_data_clean.csv'),  row.names = FALSE)

1 Participants

options(width = 100)
mytable <- xtabs(~ sex + belief, data = target_and_standard_data) / length(unique(target_and_standard_data$chindex)) / length(unique(target_and_standard_data$stimulus))
ftable(addmargins(mytable))

       belief Believer Unbeliever Sum
sex                                  
Female              24         24  48
Male                16         14  30
Sum                 40         38  78

2 Bad EEG channels

ggplot(eeg_check, aes(badchan_num)) + geom_histogram(bins = 12, color = 'darkblue', fill = 'lightblue') + scale_x_continuous(breaks=0:11)

xtabs(~ badchan_num, data = eeg_check)

badchan_num
 0  1  2  3  4  5  6  7 11 
 7 14 23 15  7  8  3  2  1 

3 ERP plots

3.1 Target and Standard:

3.2 Topographic layout:

4 General description

options(width = 100)
summary(target_and_standard_data[c('uvolts', 'sex', 'vulnerability', 'belief', 'stimulus', 'anteroposterior', 'mediolateral', 'num_id')])

     uvolts            sex           vulnerability        belief        stimulus  
 Min.   :-3.0668   Female:576   Invulnerable:780   Believer  :480   Target  :468  
 1st Qu.:-0.3816   Male  :360   Vulnerable  :156   Unbeliever:456   Standard:468  
 Median : 0.6861                                                                  
 Mean   : 1.7205                                                                  
 3rd Qu.: 3.4862                                                                  
 Max.   :10.9153                                                                  
                                                                                  
        anteroposterior  mediolateral     num_id   
 Central        :468    Left   :312   1      : 12  
 Centro-parietal:468    Midline:312   3      : 12  
                        Right  :312   6      : 12  
                                      15     : 12  
                                      16     : 12  
                                      18     : 12  
                                      (Other):864  

5 Targets and Standards by topography

options(width = 100)
target_and_standard_rep_anova = aov_ez("num_id", "uvolts", target_and_standard_data, within = c("stimulus", "anteroposterior", "mediolateral"))

target_and_standard_rain_central <- ggplot(target_and_standard_data[target_and_standard_data$anteroposterior == "Central", ], aes(y = uvolts, x = mediolateral, color = stimulus, fill = stimulus)) +
  ggtitle("Targets and Standards Central") +
  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 == "Centro-parietal", ], aes(y = uvolts, x = mediolateral, color = stimulus, fill = stimulus)) +
  ggtitle("Targets and Standards Centro-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_parietal))

target_and_standard_afex_plot <-
  afex_plot(
    target_and_standard_rep_anova,
    x = "anteroposterior",
    trace = "stimulus",
    panel = "mediolateral",
    error = "within",
    error_arg = list(width = .15),
    dodge     = my_dodge,
    data_arg  = list(
      position = 
        position_jitterdodge(
          jitter.width  = .1, 
          dodge.width   = my_dodge
        )),
    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, 77 18.45 211.26 ***  .499   <.001
2                       anteroposterior        1, 77  1.12  28.45 ***  .008   <.001
3                          mediolateral 1.67, 128.88  1.05       1.17 <.001    .308
4              stimulus:anteroposterior        1, 77  0.60  17.18 ***  .003   <.001
5                 stimulus:mediolateral 1.65, 126.97  0.80     3.47 *  .001    .043
6          anteroposterior:mediolateral 1.76, 135.18  0.29     2.88 + <.001    .066
7 stimulus:anteroposterior:mediolateral 1.95, 150.33  0.22       0.68 <.001    .507
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

Sphericity correction method: GG 

a_posteriori(target_and_standard_rep_anova)

$emmeans
 stimulus emmean   SE df lower.CL upper.CL
 Target     3.76 0.29 77    3.183    4.340
 Standard  -0.32 0.11 77   -0.539   -0.101

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

$contrasts
 contrast          estimate    SE df t.ratio p.value
 Target - Standard     4.08 0.281 77  14.535  <.0001

Results are averaged over the levels of: mediolateral, anteroposterior 

____________________________________________________________________________________________________
$emmeans
 anteroposterior emmean    SE df lower.CL upper.CL
 Central           1.54 0.171 77     1.20     1.88
 Centro.parietal   1.91 0.174 77     1.56     2.25

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

$contrasts
 contrast                  estimate     SE df t.ratio p.value
 Central - Centro.parietal   -0.369 0.0692 77  -5.334  <.0001

Results are averaged over the levels of: mediolateral, stimulus 

____________________________________________________________________________________________________
$emmeans
 stimulus anteroposterior emmean    SE df lower.CL upper.CL
 Target   Central           3.47 0.297 77    2.881   4.0624
 Standard Central          -0.40 0.115 77   -0.629  -0.1710
 Target   Centro.parietal   4.05 0.295 77    3.464   4.6375
 Standard Centro.parietal  -0.24 0.111 77   -0.461  -0.0199

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

$contrasts
 contrast                                          estimate     SE df t.ratio p.value
 Target Central - Standard Central                    3.872 0.2920 77  13.259  <.0001
 Target Central - Target Centro.parietal             -0.579 0.1103 77  -5.250  <.0001
 Target Central - Standard Centro.parietal            3.712 0.2867 77  12.946  <.0001
 Standard Central - Target Centro.parietal           -4.451 0.2916 77 -15.260  <.0001
 Standard Central - Standard Centro.parietal         -0.159 0.0504 77  -3.162  0.0118
 Target Centro.parietal - Standard Centro.parietal    4.291 0.2785 77  15.409  <.0001

Results are averaged over the levels of: mediolateral 
P value adjustment: tukey method for comparing a family of 4 estimates 

____________________________________________________________________________________________________
$emmeans
 stimulus mediolateral emmean    SE df lower.CL upper.CL
 Target   Left          3.801 0.285 77    3.233   4.3683
 Standard Left         -0.357 0.115 77   -0.587  -0.1274
 Target   Midline       3.878 0.315 77    3.251   4.5046
 Standard Midline      -0.323 0.114 77   -0.550  -0.0964
 Target   Right         3.605 0.300 77    3.007   4.2033
 Standard Right        -0.280 0.111 77   -0.500  -0.0594

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

$contrasts
 contrast                          estimate     SE df t.ratio p.value
 Target Left - Standard Left         4.1580 0.2776 77  14.979  <.0001
 Target Left - Target Midline       -0.0770 0.1346 77  -0.572  0.9926
 Target Left - Standard Midline      4.1241 0.2689 77  15.335  <.0001
 Target Left - Target Right          0.1957 0.1543 77   1.268  0.8012
 Target Left - Standard Right        4.0805 0.2625 77  15.546  <.0001
 Standard Left - Target Midline     -4.2350 0.3150 77 -13.443  <.0001
 Standard Left - Standard Midline   -0.0339 0.0434 77  -0.781  0.9700
 Standard Left - Target Right       -3.9624 0.3102 77 -12.774  <.0001
 Standard Left - Standard Right     -0.0776 0.0590 77  -1.315  0.7757
 Target Midline - Standard Midline   4.2011 0.3040 77  13.820  <.0001
 Target Midline - Target Right       0.2726 0.1027 77   2.655  0.0967
 Target Midline - Standard Right     4.1575 0.2952 77  14.082  <.0001
 Standard Midline - Target Right    -3.9284 0.3024 77 -12.990  <.0001
 Standard Midline - Standard Right  -0.0436 0.0374 77  -1.168  0.8505
 Target Right - Standard Right       3.8848 0.2899 77  13.401  <.0001

Results are averaged over the levels of: anteroposterior 
P value adjustment: tukey method for comparing a family of 6 estimates 

____________________________________________________________________________________________________

# P3_data_name   <- file.path(data_dir, 'average_voltage_275_to_425_auditory_oddball.txt')
P3_data_name   <- file.path(data_dir, 'average_voltage_300_to_400_auditory_oddball.txt')
P3_data        <- read.table(P3_data_name, header = TRUE, strip.white = TRUE, sep = "\t")
names(P3_data)[names(P3_data) == "value"] <- "uvolts"
P3_data$num_id <- readr::parse_number(P3_data$ERPset)
P3_data$vulnerability[ grepl("nVul", P3_data$ERPset)] <- "Invulnerable"
P3_data$vulnerability[!grepl("nVul", P3_data$ERPset)] <- "Vulnerable"
P3_data$belief[ grepl("nCr", P3_data$ERPset)]         <- "Unbeliever"
P3_data$belief[!grepl("nCr", P3_data$ERPset)]         <- "Believer"
P3_data$sex[ grepl("F", P3_data$ERPset)]              <- "Female"
P3_data$sex[!grepl("F", P3_data$ERPset)]              <- "Male"
P3_data$anteroposterior[P3_data$chlabel %in% c('C1', 'Cz', 'C2')]                <- 'Central'
P3_data$anteroposterior[P3_data$chlabel %in% c('E112-CP1a', 'CPz', 'E034-CP2a')] <- 'Centro-parietal'
P3_data$mediolateral[P3_data$chlabel %in% c('C1', 'E112-CP1a')]                  <- 'Left'
P3_data$mediolateral[P3_data$chlabel %in% c('Cz', 'CPz')]                        <- 'Midline'
P3_data$mediolateral[P3_data$chlabel %in% c('C2', 'E034-CP2a')]                  <- 'Right'
P3_data$num_id          <- factor(P3_data$num_id)
P3_data$vulnerability   <- factor(P3_data$vulnerability)
P3_data$sex             <- factor(P3_data$sex)
P3_data$belief          <- factor(P3_data$belief)
P3_data$anteroposterior <- factor(P3_data$anteroposterior)
P3_data$mediolateral    <- factor(P3_data$mediolateral, levels = c('Left', 'Midline', 'Right'))
if (exclude_bad_eeg) {
  P3_data <- P3_data[!(P3_data$ERPset %in% bad_eeg), ]
}
write.csv(P3_data,  file.path(data_dir, 'auditory_oddball_data_clean.csv'),  row.names = FALSE)
P3_data.table <- data.table(P3_data)
P3_by_subject <- unique(P3_data.table[, .(worklat, mlabel, uv_mean = mean(uvolts), bini, binlabel, num_id, vulnerability, belief, sex), by = .(ERPset)])
P3_by_subject$RUT <- sub("_odd.*", "", P3_by_subject$ERPset)
P3_by_subject$RUT <- sub(".*F", "", P3_by_subject$RUT)
P3_by_subject$RUT <- sub(".*M", "", P3_by_subject$RUT)
write.csv(P3_by_subject,  file.path(data_dir, 'P3_by_subject.csv'),  row.names = FALSE)

6 P3 plots

6.1 Wave difference, Target - Standard

6.2 Wave difference scalp, Target - Standard

6.3 Wave difference layout

7 P3 general description

Measurement window: [300.0 400.0] ms

options(width = 100)
summary(P3_data[c('uvolts', 'sex', 'vulnerability', 'belief', 'anteroposterior', 'mediolateral', 'num_id')])

     uvolts           sex           vulnerability        belief           anteroposterior
 Min.   :-3.501   Female:288   Invulnerable:390   Believer  :240   Central        :234   
 1st Qu.: 2.244   Male  :180   Vulnerable  : 78   Unbeliever:228   Centro-parietal:234   
 Median : 3.952                                                                          
 Mean   : 4.081                                                                          
 3rd Qu.: 5.522                                                                          
 Max.   :11.433                                                                          
                                                                                         
  mediolateral     num_id   
 Left   :156   1      :  6  
 Midline:156   3      :  6  
 Right  :156   6      :  6  
               15     :  6  
               16     :  6  
               18     :  6  
               (Other):432  

8 P3 by topography

options(width = 100)
P3_rain <- ggplot(P3_data, aes(y = uvolts, x = mediolateral, color = anteroposterior, fill = anteroposterior)) +
  ggtitle("P3") +
  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(P3_rain))

P3_rep_anova = aov_ez("num_id", "uvolts", P3_data, within = c("anteroposterior", "mediolateral"))
P3_afex_plot <-
  afex_plot(
    P3_rep_anova,
    x = "mediolateral",
    trace = "anteroposterior",
    error = "within",
    error_arg = list(width = .15),
    dodge     = my_dodge,
    data_arg  = list(
      position = 
        position_jitterdodge(
          jitter.width  = .1, 
          dodge.width   = my_dodge
        )),
    mapping   = c('color'),
    point_arg = list(size = 4)
  )
suppressWarnings(print(P3_afex_plot))

nice(P3_rep_anova)

Anova Table (Type 3 tests)

Response: uvolts
                        Effect           df  MSE         F   ges p.value
1              anteroposterior        1, 77 1.20 17.18 ***  .006   <.001
2                 mediolateral 1.65, 126.97 1.61    3.47 *  .003    .043
3 anteroposterior:mediolateral 1.95, 150.33 0.45      0.68 <.001    .507
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

Sphericity correction method: GG 

a_posteriori(P3_rep_anova)

$emmeans
 anteroposterior emmean    SE df lower.CL upper.CL
 Central           3.87 0.292 77     3.29     4.45
 Centro.parietal   4.29 0.278 77     3.74     4.85

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

$contrasts
 contrast                  estimate    SE df t.ratio p.value
 Central - Centro.parietal    -0.42 0.101 77  -4.145  0.0001

Results are averaged over the levels of: mediolateral 

____________________________________________________________________________________________________
$emmeans
 mediolateral emmean    SE df lower.CL upper.CL
 Left           4.16 0.278 77     3.61     4.71
 Midline        4.20 0.304 77     3.60     4.81
 Right          3.88 0.290 77     3.31     4.46

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

$contrasts
 contrast        estimate     SE df t.ratio p.value
 Left - Midline   -0.0431 0.1338 77  -0.322  0.9445
 Left - Right      0.2732 0.1525 77   1.792  0.1792
 Midline - Right   0.3163 0.0989 77   3.198  0.0056

Results are averaged over the levels of: anteroposterior 
P value adjustment: tukey method for comparing a family of 3 estimates 

____________________________________________________________________________________________________

9 P3 by belief and sex

options(width = 100)
P3_anova_sex = aov_ez("num_id", "uvolts", P3_data, between = c("belief", "sex"), fun_aggregate = mean)

Contrasts set to contr.sum for the following variables: belief, sex

mytable <- xtabs(~ sex + belief, data = P3_anova_sex$data$long)
ftable(addmargins(mytable))

       belief Believer Unbeliever Sum
sex                                  
Female              24         24  48
Male                16         14  30
Sum                 40         38  78

P3_afex_plot_sex <-
  afex_plot(
    P3_anova_sex,
    x     = "belief",
    trace = "sex",
    error = "between",
    error_arg = list(width = .15),
    dodge     = my_dodge,
    data_arg  = list(
      position = 
        position_jitterdodge(
          jitter.width  = .1, 
          dodge.width   = my_dodge
        )),
    mapping   = c('color'),
    point_arg = list(size = 4)
  )
suppressWarnings(print(P3_afex_plot_sex))

nice(P3_anova_sex)

Anova Table (Type 3 tests)

Response: uvolts
      Effect    df  MSE    F  ges p.value
1     belief 1, 74 6.03 0.11 .001    .746
2        sex 1, 74 6.03 2.70 .035    .105
3 belief:sex 1, 74 6.03 1.58 .021    .213
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

a_posteriori(P3_anova_sex)

10 Decoding (MVPA)

Fahrenfort, J. J., van Driel, J., van Gaal, S., & Olivers, C. N. L. (2018). From ERPs to MVPA using the Amsterdam Decoding and Modeling toolbox (ADAM). Frontiers in Neuroscience, 12(JUL), 351586. https://doi.org/10.3389/FNINS.2018.00368/BIBTEX

All channels, reference to infinite, 40 Hz low-pass
Dong, L., Li, F., Liu, Q., Wen, X., Lai, Y., Xu, P., & Yao, D. (2017). MATLAB Toolboxes for Reference Electrode Standardization Technique (REST) of Scalp EEG. Frontiers in Neuroscience, 11(OCT), 601. https://doi.org/10.3389/fnins.2017.00601

10.1 Exemplar ERP on Cz

10.2 ERP difference on Cz

10.3 Decoding over time

10.3.1 10-fold cross validation

10.3.2 Train 1&3, test 2&4

10.4 Individual results

10.5 Activation patterns

10.5.1 10-fold cross validation

10.5.2 Train 1&3, test 2&4

10.6 Temporal generalization plot

10.6.1 10-fold cross validation

10.6.2 Train 1&3, test 2&4

10.7 Generalization over time

10.7.1 10-fold cross validation

10.7.2 Train 1&3, test 2&4

10.8 Group decoding differences

10.8.1 10-fold cross validation

10.8.2 Train 1&3, test 2&4

11 Temporo-cingulate theta-alpha (4-12 Hz) Functional Connectivity (ROIconnect)

Multivariate Interaction Measure, MIM [@Pellegrini2023].
800 ms post-stimulus.

options(width = 100)
summary(connectivity_data)

                     ERPset        stimulus             region_a               region_b  
 S001VulCrF20957113_odd :  8   target  :312   left_temporal :312   cingulate_ipsi  :312  
 S003nVulCrM18466555_odd:  8   standard:312   right_temporal:312   cingulate_contra:312  
 S006nVulCrM17923449_odd:  8                                                             
 S015nVulCrM18833005_odd:  8                                                             
 S016VulCrF19423156_odd :  8                                                             
 S018nVulCrF17671904_odd:  8                                                             
 (Other)                :576                                                             
 median_connectivity     sex      log_median_connectivity     num_id     
 Min.   :0.03258     Female:384   Min.   :-3.424          Min.   :  1.0  
 1st Qu.:0.05716     Male  :240   1st Qu.:-2.862          1st Qu.: 47.0  
 Median :0.07289                  Median :-2.619          Median : 76.5  
 Mean   :0.08027                  Mean   :-2.590          Mean   : 72.1  
 3rd Qu.:0.09631                  3rd Qu.:-2.340          3rd Qu.:101.0  
 Max.   :0.21169                  Max.   :-1.553          Max.   :120.0  
                                                                         

connectivity_data_wide_cingulate <- connectivity_data[c('ERPset', 'region_a', 'region_b', 'stimulus', 'log_median_connectivity')] |>
  pivot_wider(names_from = region_b, values_from = log_median_connectivity)
connectivity_data_wide_temporal <- connectivity_data[c('ERPset', 'region_a', 'region_b', 'stimulus', 'log_median_connectivity')] |>
  pivot_wider(names_from = region_a, values_from = log_median_connectivity)
connectivity_data_wide_stimulus   <- connectivity_data[c('ERPset', 'region_a', 'region_b', 'stimulus', 'log_median_connectivity')] |>
  pivot_wider(names_from = stimulus, values_from = log_median_connectivity)
stimulus   <- c('target'         , 'standard')
temporal   <- c('left_temporal'  , 'right_temporal')
cingulate <- c('cingulate_ipsi', 'cingulate_contra')
cingulate_pairs <- ggpairs(connectivity_data_wide_cingulate,
                            aes(colour = stimulus, alpha = .5, linewidth = NA),
                            columns  = cingulate,
                            lower    = list(continuous = wrap('points', alpha = .4)),
                            progress = FALSE
                            )
temporal_pairs <- ggpairs(connectivity_data_wide_temporal,
                            aes(colour = stimulus, alpha = .5, linewidth = NA),
                            columns  = temporal,
                            lower    = list(continuous = wrap('points', alpha = .4)),
                            progress = FALSE
                            )
stimulus_pairs <- ggpairs(connectivity_data_wide_stimulus,
                            aes(colour = region_b, alpha = .5, linewidth = NA),
                            columns  = stimulus,
                            lower    = list(continuous = wrap('points', alpha = .4)),
                            progress = FALSE
                            )
suppressWarnings(print(cingulate_pairs))

suppressWarnings(print(temporal_pairs))

suppressWarnings(print(stimulus_pairs))

11.1 Connectivity ANOVA:

connectivity_rep_anova <- aov_ez('num_id', 'log_median_connectivity', 
                             connectivity_data,
                             within = c('stimulus', 'region_a', 'region_b'))
connectivity_afex_plot <- afex_plot(connectivity_rep_anova,
                                x     = 'region_a',
                                trace = 'stimulus',
                                panel = 'region_b',
                                error = 'within',
                                error_arg = list(width = .3, lwd = .75),
                                dodge = .5,
                                data_arg  = list(
                                  position = 
                                    position_jitterdodge(
                                      jitter.width = .2, 
                                      dodge.width  = .5  ## needs to be same as dodge
                                    )),
                                mapping   = c('color'),, data_alpha = .3,
                                point_arg = list(size = 3)
                                )
suppressWarnings(print(connectivity_afex_plot))

mytable <- xtabs(~ stimulus + region_a + region_b, data = connectivity_rep_anova$data$long)
ftable(mytable)

                        region_b cingulate_ipsi cingulate_contra
stimulus region_a                                               
target   left_temporal                       78               78
         right_temporal                      78               78
standard left_temporal                       78               78
         right_temporal                      78               78

print(connectivity_rep_anova)

Anova Table (Type 3 tests)

Response: log_median_connectivity
                      Effect    df  MSE         F   ges p.value
1                   stimulus 1, 77 0.18 49.48 ***  .109   <.001
2                   region_a 1, 77 0.04    5.32 *  .003    .024
3                   region_b 1, 77 0.00   8.15 ** <.001    .006
4          stimulus:region_a 1, 77 0.03      1.19 <.001    .279
5          stimulus:region_b 1, 77 0.00      0.51 <.001    .476
6          region_a:region_b 1, 77 0.00    3.16 + <.001    .079
7 stimulus:region_a:region_b 1, 77 0.00      0.01 <.001    .914
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

a_posteriori(connectivity_rep_anova)

$emmeans
 stimulus emmean     SE df lower.CL upper.CL
 target    -2.47 0.0381 77    -2.55    -2.40
 standard  -2.71 0.0360 77    -2.78    -2.64

Results are averaged over the levels of: region_b, region_a 
Confidence level used: 0.95 

$contrasts
 contrast          estimate     SE df t.ratio p.value
 target - standard    0.238 0.0339 77   7.034  <.0001

Results are averaged over the levels of: region_b, region_a 

____________________________________________________________________________________________________
$emmeans
 region_a       emmean     SE df lower.CL upper.CL
 left_temporal   -2.61 0.0332 77    -2.68    -2.54
 right_temporal  -2.57 0.0347 77    -2.64    -2.50

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

$contrasts
 contrast                       estimate     SE df t.ratio p.value
 left_temporal - right_temporal  -0.0386 0.0167 77  -2.307  0.0237

Results are averaged over the levels of: region_b, stimulus 

____________________________________________________________________________________________________
$emmeans
 region_b         emmean     SE df lower.CL upper.CL
 cingulate_ipsi    -2.59 0.0334 77    -2.66    -2.53
 cingulate_contra  -2.59 0.0325 77    -2.65    -2.52

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

$contrasts
 contrast                          estimate      SE df t.ratio p.value
 cingulate_ipsi - cingulate_contra -0.00865 0.00303 77  -2.855  0.0055

Results are averaged over the levels of: region_a, stimulus 

____________________________________________________________________________________________________

11.2 Connectivity Assumptions:

norm_connectivity <- check_normality(connectivity_rep_anova)
print(norm_connectivity)

Warning: Non-normality of residuals detected (p < .001).

plot(norm_connectivity, type = 'density')

plot(norm_connectivity, type = 'qq')

hetero_connectivity <- check_heteroscedasticity(connectivity_rep_anova)

Data was changed during ANOVA calculation. Thus, fitted values cannot be added to original data.
fitted(..., append = TRUE) will return data and fitted values.

print(hetero_connectivity)

OK: Error variance appears to be homoscedastic (p = 0.385).

plot(hetero_connectivity)

Data was changed during ANOVA calculation. Thus, residuals cannot be added to original data.
residuals(..., append = TRUE) will return data and residuals.
Data was changed during ANOVA calculation. Thus, fitted values cannot be added to original data.
fitted(..., append = TRUE) will return data and fitted values.

check_sphericity(connectivity_rep_anova)

Warning in min(x) : no non-missing arguments to min; returning Inf

OK: Data seems to be spherical (p > .999).