Code
cat('\014') # clean terminalCode
rm(list = ls()) # clean workspace
library(tidyverse)
library(readxl)
library(afex)
library(emmeans)
library(easystats)ERPs Pain Empathy 2015 with SASS, stationary device, version1
Recordings with BioSemi ActiveTwo, 64 channels
Code
my_dodge <- .3
my_jitter <- .2
theme_set(
theme_minimal()
)
a_posteriori_aov_ez <- function(aov_ez_object, sig_level = .05) {
factors <- as.list(rownames(aov_ez_object$anova_table))
for (j in 1:length(factors)) {
if (grepl(':', factors[[j]])) {
factors[[j]] <- unlist(strsplit(factors[[j]], ':'))
}
}
p_values <- aov_ez_object$anova_table$`Pr(>F)`
for (i in 1:length(p_values)) {
if (p_values[i] <= sig_level) {
cat(rep('_', 60), '\n', sep = '')
print(emmeans(aov_ez_object, factors[[i]], contr = 'pairwise'))
}
}
}
a_posteriori_lmer <- function(lmer_obj, sig_level = .05) {
anova_lmer <- anova(lmer_obj)
factors <- as.list(row.names(anova_lmer))
for (j in 1:length(factors)) {
if (grepl(':', factors[[j]])) {
factors[[j]] <- unlist(strsplit(factors[[j]], ':'))
}
}
p_values <- anova_lmer$`Pr(>F)`
for (i in 1:length(p_values)) {
if (p_values[i] <= sig_level) {
cat(rep('_', 60), '\n', sep = '')
print(emmeans(lmer_obj, factors[[i]], contr = 'pairwise'))
}
}
}
xclude <- c(8, 23, 59)Code
SASS_data <- read_xlsx(file.path('..', 'bad channels empathy huepe.xlsx')) |>
rename(SASS = SASS_TOTAL) |>
mutate(num_id = parse_number(ID_gsheet)) |>
filter(to_kill ==0)
erp_averages <- list.files(path = './data/ave_volt', full.names = TRUE)
painEmpathy_data <- vroom::vroom(erp_averages, show_col_types = FALSE) |>
separate(ERPset, c('ID', 'Sex'), sep = '_', remove = FALSE) |>
rename(Component = mlabel,
Amplitude = value,
Stimulus = binlabel,
Electrode = chlabel) |>
mutate(Stimulus = case_match(Stimulus, 'NoPain' ~ 'no_pain', 'Pain' ~ 'pain'),
worklat = gsub('.0', '', worklat, fixed = TRUE),
num_id = parse_number(ID),
chindex = factor(chindex),
bini = factor(bini)) |>
mutate_if(is.character, as.factor) |>
mutate(Component = factor(Component, levels = c('N2', 'N3', 'P2', 'P3', 'LPP')),
Electrode = factor(Electrode, levels = c('FC1', 'Fz', 'FC2',
'CP1', 'Cz', 'CP2',
'O1' , 'Oz', 'O2')),
) |>
left_join(SASS_data[c('num_id', 'SASS')], by = 'num_id') |>
mutate(c_SASS = SASS - mean(SASS, na.rm = TRUE)) |>
filter(!(num_id %in% xclude))Warning: Expected 2 pieces. Additional pieces discarded in 1770 rows [1, 2, 3, 4, 5, 6,
7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].Code
write.csv(painEmpathy_data, file.path('data/painEmpathy_2015_erp_clean_by_subject.csv'), row.names = FALSE)ERP plots
BDF files were pre-processed in Matlab using the EEGLAB 2023.1 toolbox (Delorme and Makeig 2004), the ERPLAB 10.0 toolbox (Lopez-Calderon and Luck 2014), and automated with in-house scripts. EEG data was high-pass filtered at 0.1 Hz with a 12 dB/oct roll-off. A working time window was selected from 2 seconds before the first Stimulus to 2 seconds after the last one. Defective channels were identified by eye inspection and omitted from the preprocessing steps. Artifacts originating from eye blinks or movements, heart beats, muscle contraction, channel noise, or electrical interference were identified and removed by means of independent Component analysis (ICA) and using the ICLabel 1.4 classifier (Pion-Tonachini, Kreutz-Delgado, and Makeig 2019), Components with an score between 0.8 and 1 in the aforementioned artifactual categories were removed. At this point defective channels were spherically interpolated. Data was re-referenced to infinity with the REST 1.2 toolbox (Dong et al. 2017), low-pass filtered at 35 Hz with a 12 dB/oct roll-off, cut in [-200ms 1000ms] epochs around stimuli, and baseline-corrected relative to the mean voltage in the pre-Stimulus time.
Topographic layout:
ERPs General description
Code
ggplot(
painEmpathy_data, aes(x = Amplitude, fill = Component, color = Component)) +
geom_histogram(alpha = .4) +
facet_wrap(~Component, ncol = 1) +
theme(strip.text.x = element_blank())
addmargins(xtabs(~ Sex, data = unique(painEmpathy_data[c('Sex', 'num_id')])))Sex
female male Sum
34 25 59 Code
summary(painEmpathy_data) worklat Component Amplitude chindex Electrode
[180 220]:354 N2 :354 Min. :-17.296 11 :236 FC1 :236
[190 230]:354 N3 :354 1st Qu.: -2.119 27 :236 Fz :236
[300 380]:354 P2 :354 Median : 1.027 29 :236 FC2 :236
[330 450]:354 P3 :354 Mean : 1.171 38 :236 O1 :236
[600 800]:354 LPP:354 3rd Qu.: 4.332 46 :236 Oz :236
Max. : 19.694 64 :236 O2 :236
(Other):354 (Other):354
bini Stimulus ERPset ID
1:885 no_pain:885 S01_female_painEmpathy: 30 S01 : 30
2:885 pain :885 S03_female_painEmpathy: 30 S03 : 30
S04_male_painEmpathy : 30 S04 : 30
S06_male_painEmpathy : 30 S06 : 30
S07_male_painEmpathy : 30 S07 : 30
S11_female_painEmpathy: 30 S11 : 30
(Other) :1590 (Other):1590
Sex num_id SASS c_SASS
female:1020 Min. : 1.00 Min. :27.00 Min. :-15.5172
male : 750 1st Qu.:27.00 1st Qu.:36.00 1st Qu.: -6.5172
Median :43.00 Median :42.00 Median : -0.5172
Mean :42.22 Mean :42.52 Mean : 0.0000
3rd Qu.:58.00 3rd Qu.:49.00 3rd Qu.: 6.4828
Max. :99.00 Max. :55.00 Max. : 12.4828
NA's :30 NA's :30 Code
n2_data <- subset(painEmpathy_data, Component == 'N2')
n3_data <- subset(painEmpathy_data, Component == 'N3')
p2_data <- subset(painEmpathy_data, Component == 'P2')
p3_data <- subset(painEmpathy_data, Component == 'P3')
lpp_data <- subset(painEmpathy_data, Component == 'LPP')
In the front: N2 & N3:
N2, average from [190 230] ms interval
Electrodes FC1, Fz, FC2
Code
painEmpathy_n2_lmer <- lmer(Amplitude ~ c_SASS + Stimulus * Sex + c_SASS : Sex + (Stimulus | num_id) + (1 | num_id : Electrode), n2_data)
afex_plot(
painEmpathy_n2_lmer,
x = 'Stimulus',
trace = 'Sex',
# panel = 'Sex',
id = 'num_id',
error_arg = list(width = .1),
dodge = my_dodge,
data_arg = list(
position =
position_jitterdodge(
jitter.width = my_jitter,
jitter.height = 0,
dodge.width = my_dodge ## needs to be same as dodge
)),
mapping = c('color'),
point_arg = list(size = 3)
)
Code
anova(painEmpathy_n2_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
c_SASS 0.64725 0.64725 1 53.999 1.8305 0.1817
Stimulus 0.54772 0.54772 1 55.999 1.5490 0.2185
Sex 0.14635 0.14635 1 54.001 0.4139 0.5227
Stimulus:Sex 0.00622 0.00622 1 55.999 0.0176 0.8950
c_SASS:Sex 0.37604 0.37604 1 53.999 1.0635 0.3070Code
omega_squared(painEmpathy_n2_lmer)# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI
----------------------------------------------
c_SASS | 0.01 | [0.00, 1.00]
Stimulus | 9.38e-03 | [0.00, 1.00]
Sex | 0.00 | [0.00, 1.00]
Stimulus:Sex | 0.00 | [0.00, 1.00]
c_SASS:Sex | 1.13e-03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Code
a_posteriori_lmer(painEmpathy_n2_lmer)N3, average from [300 380] ms interval
Electrodes FC1, Fz, FC2
Code
painEmpathy_n3_lmer <- lmer(Amplitude ~ c_SASS + Stimulus * Sex + c_SASS : Sex + (Stimulus | num_id) + (1 | num_id : Electrode), n3_data)
afex_plot(
painEmpathy_n3_lmer,
x = 'Stimulus',
trace = 'Sex',
# panel = 'Sex',
id = 'num_id',
error_arg = list(width = .1),
dodge = my_dodge,
data_arg = list(
position =
position_jitterdodge(
jitter.width = my_jitter,
jitter.height = 0,
dodge.width = my_dodge ## needs to be same as dodge
)),
mapping = c('color'),
point_arg = list(size = 3)
)
Code
anova(painEmpathy_n3_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
c_SASS 0.29293 0.29293 1 53.999 0.8573 0.3586
Stimulus 0.49114 0.49114 1 55.999 1.4373 0.2356
Sex 0.05516 0.05516 1 54.015 0.1614 0.6894
Stimulus:Sex 0.06549 0.06549 1 55.999 0.1917 0.6632
c_SASS:Sex 0.11460 0.11460 1 53.999 0.3354 0.5649Code
omega_squared(painEmpathy_n3_lmer)# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI
----------------------------------------------
c_SASS | 0.00 | [0.00, 1.00]
Stimulus | 7.48e-03 | [0.00, 1.00]
Sex | 0.00 | [0.00, 1.00]
Stimulus:Sex | 0.00 | [0.00, 1.00]
c_SASS:Sex | 0.00 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Code
a_posteriori_lmer(painEmpathy_n3_lmer)In the back: P2 & P3:
P2, average [180 220] ms interval
Electrodes O1, Oz, O2
Code
painEmpathy_p2_lmer <- lmer(Amplitude ~ c_SASS + Stimulus * Sex + c_SASS : Sex + (Stimulus | num_id) + (1 | num_id : Electrode), p2_data)
afex_plot(
painEmpathy_p2_lmer,
x = 'Stimulus',
trace = 'Sex',
# panel = 'Sex',
id = 'num_id',
error_arg = list(width = .1),
dodge = my_dodge,
data_arg = list(
position =
position_jitterdodge(
jitter.width = my_jitter,
jitter.height = 0,
dodge.width = my_dodge ## needs to be same as dodge
)),
mapping = c('color'),
point_arg = list(size = 3)
)
Code
anova(painEmpathy_p2_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
c_SASS 0.48729 0.48729 1 54.002 0.9372 0.3373
Stimulus 0.00153 0.00153 1 56.004 0.0029 0.9570
Sex 0.61429 0.61429 1 54.001 1.1814 0.2819
Stimulus:Sex 0.01606 0.01606 1 56.004 0.0309 0.8611
c_SASS:Sex 1.25955 1.25955 1 54.002 2.4224 0.1255Code
omega_squared(painEmpathy_p2_lmer)# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI
----------------------------------------------
c_SASS | 0.00 | [0.00, 1.00]
Stimulus | 0.00 | [0.00, 1.00]
Sex | 3.23e-03 | [0.00, 1.00]
Stimulus:Sex | 0.00 | [0.00, 1.00]
c_SASS:Sex | 0.02 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Code
a_posteriori_lmer(painEmpathy_p2_lmer)P3, average from [330 450] ms interval
Electrodes O1, Oz, O2
Code
painEmpathy_p3_lmer <- lmer(Amplitude ~ c_SASS + Stimulus * Sex + c_SASS : Sex + (Stimulus | num_id) + (1 | num_id : Electrode), p3_data)
afex_plot(
painEmpathy_p3_lmer,
x = 'Stimulus',
trace = 'Sex',
# panel = 'Sex',
id = 'num_id',
error_arg = list(width = .1),
dodge = my_dodge,
data_arg = list(
position =
position_jitterdodge(
jitter.width = my_jitter,
jitter.height = 0,
dodge.width = my_dodge ## needs to be same as dodge
)),
mapping = c('color'),
point_arg = list(size = 3)
)
Code
anova(painEmpathy_p3_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
c_SASS 0.6133 0.6133 1 54.000 0.6249 0.43270
Stimulus 5.4503 5.4503 1 56.000 5.5531 0.02197 *
Sex 0.8801 0.8801 1 53.998 0.8967 0.34789
Stimulus:Sex 0.0135 0.0135 1 56.000 0.0137 0.90708
c_SASS:Sex 0.4592 0.4592 1 54.000 0.4678 0.49692
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
omega_squared(painEmpathy_p3_lmer)# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI
----------------------------------------------
c_SASS | 0.00 | [0.00, 1.00]
Stimulus | 0.07 | [0.00, 1.00]
Sex | 0.00 | [0.00, 1.00]
Stimulus:Sex | 0.00 | [0.00, 1.00]
c_SASS:Sex | 0.00 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Code
a_posteriori_lmer(painEmpathy_p3_lmer)____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
Stimulus emmean SE df lower.CL upper.CL
no_pain 6.56 0.630 54.1 5.30 7.82
pain 6.01 0.635 54.1 4.74 7.28
Results are averaged over the levels of: Sex
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain 0.548 0.232 56 2.357 0.0220
Results are averaged over the levels of: Sex
Degrees-of-freedom method: kenward-roger Late Centro-parietal: LPP
LPP, average from [600 800] ms interval
Electrodes CP1, Cz, CP2
Code
painEmpathy_lpp_lmer <- lmer(Amplitude ~ c_SASS + Stimulus * Sex + c_SASS : Sex + (Stimulus | num_id) + (1 | num_id : Electrode), lpp_data)
afex_plot(
painEmpathy_lpp_lmer,
x = 'Stimulus',
trace = 'Sex',
# panel = 'Sex',
id = 'num_id',
error_arg = list(width = .1),
dodge = my_dodge,
data_arg = list(
position =
position_jitterdodge(
jitter.width = my_jitter,
jitter.height = 0,
dodge.width = my_dodge ## needs to be same as dodge
)),
mapping = c('color'),
point_arg = list(size = 3)
)
Code
options(width = 120)
anova(painEmpathy_lpp_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
c_SASS 0.2409 0.2409 1 54.006 0.4601 0.5005
Stimulus 23.6045 23.6045 1 56.001 45.0925 1.022e-08 ***
Sex 1.0532 1.0532 1 54.086 2.0121 0.1618
Stimulus:Sex 0.0149 0.0149 1 56.001 0.0285 0.8666
c_SASS:Sex 0.1510 0.1510 1 54.006 0.2884 0.5934
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
omega_squared(painEmpathy_lpp_lmer)# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI
----------------------------------------------
c_SASS | 0.00 | [0.00, 1.00]
Stimulus | 0.43 | [0.27, 1.00]
Sex | 0.02 | [0.00, 1.00]
Stimulus:Sex | 0.00 | [0.00, 1.00]
c_SASS:Sex | 0.00 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Code
a_posteriori_lmer(painEmpathy_lpp_lmer)____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
Stimulus emmean SE df lower.CL upper.CL
no_pain 1.68 0.276 54.1 1.13 2.23
pain 2.89 0.305 54.5 2.28 3.50
Results are averaged over the levels of: Sex
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain -1.21 0.18 56 -6.715 <.0001
Results are averaged over the levels of: Sex
Degrees-of-freedom method: kenward-roger Mass Univariate analysis (cluster-mass), full time window
Groppe, Urbach, and Kutas (2011)
max_dist value of 50 corresponds to an approximate distance of 5.24 cm (assuming a 56 cm great circle circumference head and that your electrode coordinates are based on an idealized spherical head with radius of 85.000000). Min/Max distances between all pairs of channels (in chanlocs units): 25.654041/169.962716 Median (semi-IQR) distance between all pairs of channels (in chanlocs units): 113.111431 (33.320704) Mean (SD) # of neighbors per channel: 6.5 (1.5) Median (semi-IQR) # of neighbors per channel: 7.0 (1.5) Min/max # of neighbors per channel: 3 to 8
Found 1 significant positive cluster (out of 11):
p < .0001, mass = 7996.
Found 1 significant negative cluster (out of 12):
p = .0042, mass = -2396.
Multivariate pattern analysis (decoding)
No processing besides a 0.1 Hz high-pass filter (12 dB/oct roll-off) and reference to infinity (Dong et al. 2017).
Parietal and occipital electrodes.
Done in the ADAM 1.14 toolbox (Fahrenfort 2020) with FieldTrip 20211209 (Oostenveld et al. 2011).
Exemplar ERPs on Pz
Decoding over time
Temporal generalization plot
Activation pattern
Temporal generalization over time
Individual results
References
Delorme, Arnaud, and Scott Makeig. 2004. “EEGLAB: An Open Source Toolbox for Analysis of Single-Trial EEG Dynamics Including Independent Component Analysis.” Journal of Neuroscience Methods 134 (March): 9–21. https://doi.org/10.1016/J.JNEUMETH.2003.10.009.
Dong, Li, Fali Li, Qiang Liu, Xin Wen, Yongxiu Lai, Peng Xu, and Dezhong Yao. 2017. “MATLAB Toolboxes for Reference Electrode Standardization Technique (REST) of Scalp EEG.” Frontiers in Neuroscience 11 (October): 601. https://doi.org/10.3389/fnins.2017.00601.
Fahrenfort, Johannes Jacobus. 2020. “Multivariate Methods to Track the Spatiotemporal Profile of Feature-Based Attentional Selection Using EEG.” Neuromethods 151: 129–56. https://doi.org/10.1007/7657_2019_26/COVER.
Groppe, David M., Thomas P. Urbach, and Marta Kutas. 2011. “Mass Univariate Analysis of Event-Related Brain Potentials/Fields i: A Critical Tutorial Review.” Psychophysiology 48 (December): 1711–25. https://doi.org/10.1111/J.1469-8986.2011.01273.X.
Lopez-Calderon, Javier, and Steven J. Luck. 2014. “ERPLAB: An Open-Source Toolbox for the Analysis of Event-Related Potentials.” Frontiers in Human Neuroscience 8 (April): 75729. https://doi.org/10.3389/FNHUM.2014.00213/BIBTEX.
Oostenveld, Robert, Pascal Fries, Eric Maris, and Jan Mathijs Schoffelen. 2011. “FieldTrip: Open Source Software for Advanced Analysis of MEG, EEG, and Invasive Electrophysiological Data.” Computational Intelligence and Neuroscience 2011. https://doi.org/10.1155/2011/156869.
Pion-Tonachini, Luca, Ken Kreutz-Delgado, and Scott Makeig. 2019. “ICLabel: An Automated Electroencephalographic Independent Component Classifier, Dataset, and Website.” NeuroImage 198: 181–97. https://doi.org/10.1016/j.neuroimage.2019.05.026.