Code
cat('\014') # clean terminalCode
rm(list = ls()) # clean workspace
library(tidyverse)
library(afex)
library(emmeans)
library(easystats)ERPs Pain Empathy 2024, stationary device, all stimuli
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(666)Code
csv_answers <- list.files(path = '../csv', full.names = TRUE)
answers_df <- vroom::vroom(csv_answers, col_types = c(Sex = 'c'), show_col_types = FALSE) |>
filter(Block > 0) |>
mutate(Block = factor(Block)) |>
filter(!is.na(Question)) |>
mutate(Sex = if_else(Sex == 'f', 'female', 'male')) |>
mutate(Answer = if_else(Answer > 97, 97, Answer)) |>
mutate(Answer = if_else(Answer < 1, 1, Answer)) |>
separate(Subject, c('ID', 'Version'), sep = '_', remove = FALSE, extra = 'drop') |>
mutate(Version = tolower(Version),
num_id = parse_number(ID),
Question = recode(Question, 'unpleasantness' = 'displeasing', 'pain' = 'painful'),
Stimulus = case_match(Valence, 'EASE' ~ 'no_pain', 'PAIN' ~ 'pain')
) |>
mutate_if(is.character, as.factor) |>
relocate(Subject, .after = last_col()) |>
filter(!(num_id %in% xclude))
write.csv(answers_df, 'data/painEmpathy_2024_answers_clean_by_subject.csv', row.names = FALSE)
erp_averages <- list.files(path = './data/ave_volt_by_subject', full.names = TRUE)
painEmpathy_data <- vroom::vroom(erp_averages, show_col_types = FALSE) |>
rename(Component = mlabel,
Amplitude = value,
Stimulus = binlabel,
Electrode = chlabel) |>
mutate(worklat = gsub('.0', '', worklat, fixed = TRUE),
num_id = parse_number(ERPset),
chindex = factor(chindex),
bini = factor(bini)) |>
left_join(unique(answers_df[c('num_id', 'Sex')]), by = 'num_id') |> mutate_if(is.character, as.factor) |>
mutate_if(is.character, as.factor) |>
mutate(Component = factor(Component, levels = c('N1', 'N2', 'P1', 'P3', 'LPP')),
Electrode = factor(Electrode, levels = c('FC1', 'Fz', 'FC2',
'CP1', 'Cz', 'CP2',
'O1' , 'Oz', 'O2')),
) |>
filter(!(num_id %in% xclude))
write.csv(painEmpathy_data, file.path('data/painEmpathy_2024_erp_clean_by_subject.csv'), row.names = FALSE)Answers
Code
addmargins(xtabs(~ Sex, data = unique(answers_df[c('Sex', 'num_id')])))Sex
female male Sum
19 11 30 Code
summary(answers_df) ID Version Sex Block trialN
pilot07: 121 v1:1425 female:1895 1:1453 Min. : 1.00
pilot04: 106 v2:1495 male :1025 2:1467 1st Qu.:16.00
pilot09: 106 Median :32.00
pilot30: 105 Mean :31.85
pilot14: 104 3rd Qu.:47.00
pilot22: 104 Max. :64.00
(Other):2274
Image Valence Answer Question rT
hp-5 : 19 EASE:1434 Min. : 1.00 displeasing:1479 Min. : 267
10.1.jpg: 18 PAIN:1486 1st Qu.: 1.00 painful :1441 1st Qu.: 3711
2.2.jpg : 18 Median :49.00 Median : 4543
fp-28 : 18 Mean :46.73 Mean : 4993
hp-27 : 18 3rd Qu.:87.00 3rd Qu.: 5691
fnp-19 : 17 Max. :97.00 Max. :94218
(Other) :2812
Session num_id Stimulus
Min. :1.000 Min. : 1.00 no_pain:1434
1st Qu.:1.000 1st Qu.: 8.00 pain :1486
Median :1.000 Median :15.00
Mean :1.054 Mean :15.55
3rd Qu.:1.000 3rd Qu.:24.00
Max. :3.000 Max. :31.00
Subject
pilot07_v1_2023-11-02_14-27-38: 62
pilot04_v2 : 60
pilot01_v2 : 59
pilot07_v2_2023-11-02_14-12-05: 59
pilot18_v2_2023-11-14_21-04-52: 57
pilot22_v2_2025-03-17_17-39-34: 56
(Other) :2567 Pain rating:
Code
pain_rating_lmer <- lmer(Answer ~ Stimulus + (Stimulus|num_id), subset(answers_df, Question == 'painful'))
afex_plot(
pain_rating_lmer,
x = 'Stimulus',
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 = 4)
)
Code
anova(pain_rating_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Stimulus 145475 145475 1 29.485 423.99 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(pain_rating_lmer), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
------------------------------------------------------------
Stimulus | 0.93 | [0.89, 1.00] | large
- One-sided CIs: upper bound fixed at [1.00].
- Interpretation rule: field2013Code
a_posteriori_lmer(pain_rating_lmer)____________________________________________________________
$emmeans
Stimulus emmean SE df lower.CL upper.CL
no_pain 12.3 2.25 29 7.73 16.9
pain 81.0 2.08 29 76.76 85.3
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain -68.7 3.34 29 -20.590 <.0001
Degrees-of-freedom method: kenward-roger Unpleasantness rating:
Code
unpleasantness_rating_lmer <- lmer(Answer ~ Stimulus + (Stimulus|num_id), subset(answers_df, Question == 'displeasing'))
afex_plot(
unpleasantness_rating_lmer,
x = 'Stimulus',
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 = 4)
)
Code
anova(unpleasantness_rating_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Stimulus 164574 164574 1 29.193 356.74 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(unpleasantness_rating_lmer), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
------------------------------------------------------------
Stimulus | 0.92 | [0.87, 1.00] | large
- One-sided CIs: upper bound fixed at [1.00].
- Interpretation rule: field2013Code
a_posteriori_lmer(unpleasantness_rating_lmer)____________________________________________________________
$emmeans
Stimulus emmean SE df lower.CL upper.CL
no_pain 16.6 2.29 29 11.9 21.2
pain 74.7 2.44 29 69.7 79.7
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain -58.2 3.08 29 -18.886 <.0001
Degrees-of-freedom method: kenward-roger 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
19 11 30 Code
summary(painEmpathy_data) worklat Component Amplitude chindex Electrode
[140 170]:180 N1 :180 Min. :-11.954 11 :120 FC1 :120
[150 180]:180 N2 :180 1st Qu.: -2.442 27 :120 Fz :120
[270 330]:180 P1 :180 Median : 1.085 29 :120 FC2 :120
[280 380]:180 P3 :180 Mean : 0.700 38 :120 O1 :120
[600 800]:180 LPP:180 3rd Qu.: 3.332 46 :120 Oz :120
Max. : 15.005 64 :120 O2 :120
(Other):180 (Other):180
bini Stimulus ERPset num_id Sex
5:450 no_pain:450 pilot01: 30 Min. : 1.00 female:570
6:450 pain :450 pilot02: 30 1st Qu.: 8.00 male :330
pilot03: 30 Median :15.50
pilot04: 30 Mean :15.77
pilot05: 30 3rd Qu.:24.00
pilot06: 30 Max. :31.00
(Other):720 Code
n1_data <- subset(painEmpathy_data, Component == 'N1',)
n2_data <- subset(painEmpathy_data, Component == 'N2',)
p1_data <- subset(painEmpathy_data, Component == 'P1',)
p3_data <- subset(painEmpathy_data, Component == 'P3',)
lpp_data <- subset(painEmpathy_data, Component == 'LPP',)
In the front: N1 & N2:
N1, average from [150 180] ms interval
Electrodes FC1, Fz, FC2
Code
painEmpathy_n1_lmer <- lmer(Amplitude ~ Stimulus + (Stimulus | num_id) + (1 | num_id : Electrode), n1_data)
afex_plot(
painEmpathy_n1_lmer,
x = 'Stimulus',
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_n1_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Stimulus 0.11264 0.11264 1 29 2.9066 0.09891 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
# cat(rep('_', 60), '\n', sep = '')
interpret(omega_squared(painEmpathy_n1_lmer), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
------------------------------------------------------------
Stimulus | 0.06 | [0.00, 1.00] | small
- One-sided CIs: upper bound fixed at [1.00].
- Interpretation rule: field2013Code
a_posteriori_lmer(painEmpathy_n1_lmer)
# cat(rep('_', 60), '\n', sep = '')
# emmeans(painEmpathy_n1_lmer, pairwise ~ Group | Sex)N2, average from [270 330] ms interval
Electrodes FC1, Fz, FC2
Code
painEmpathy_n2_lmer <- lmer(Amplitude ~ Stimulus + (Stimulus | num_id) + (1 | num_id : Electrode), n2_data)
afex_plot(
painEmpathy_n2_lmer,
x = 'Stimulus',
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)
Stimulus 0.034502 0.034502 1 28.996 1.0197 0.3209Code
interpret(omega_squared(painEmpathy_n2_lmer), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
------------------------------------------------------------
Stimulus | 6.35e-04 | [0.00, 1.00] | very small
- One-sided CIs: upper bound fixed at [1.00].
- Interpretation rule: field2013Code
a_posteriori_lmer(painEmpathy_n2_lmer)In the back: P1 & P3:
P1, average [140 170] ms interval
Electrodes O1, Oz, O2
Code
painEmpathy_p1_lmer <- lmer(Amplitude ~ Stimulus + (Stimulus | num_id) + (1 | num_id : Electrode), p1_data)
afex_plot(
painEmpathy_p1_lmer,
x = 'Stimulus',
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_p1_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Stimulus 0.00035457 0.00035457 1 29.002 0.0098 0.922Code
interpret(omega_squared(painEmpathy_p1_lmer), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
------------------------------------------------------------
Stimulus | 0.00 | [0.00, 1.00] | very small
- One-sided CIs: upper bound fixed at [1.00].
- Interpretation rule: field2013Code
a_posteriori_lmer(painEmpathy_p1_lmer)P3, average from [280 380] ms interval
Electrodes O1, Oz, O2
Code
painEmpathy_p3_lmer <- lmer(Amplitude ~ Stimulus + (Stimulus | num_id) + (1 | num_id : Electrode), p3_data)
afex_plot(
painEmpathy_p3_lmer,
x = 'Stimulus',
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)
Stimulus 0.90451 0.90451 1 29.002 23.407 3.977e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(painEmpathy_p3_lmer), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
------------------------------------------------------------
Stimulus | 0.42 | [0.19, 1.00] | large
- One-sided CIs: upper bound fixed at [1.00].
- Interpretation rule: field2013Code
a_posteriori_lmer(painEmpathy_p3_lmer)____________________________________________________________
$emmeans
Stimulus emmean SE df lower.CL upper.CL
no_pain 5.18 0.586 29 3.99 6.38
pain 4.60 0.566 29 3.45 5.76
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain 0.582 0.12 29 4.838 <.0001
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 ~ Stimulus + (Stimulus | num_id) + (1 | num_id : Electrode), lpp_data)
afex_plot(
painEmpathy_lpp_lmer,
x = 'Stimulus',
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)
Stimulus 1.3775 1.3775 1 29 60.339 1.447e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(painEmpathy_lpp_lmer), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
------------------------------------------------------------
Stimulus | 0.66 | [0.47, 1.00] | large
- One-sided CIs: upper bound fixed at [1.00].
- Interpretation rule: field2013Code
a_posteriori_lmer(painEmpathy_lpp_lmer)____________________________________________________________
$emmeans
Stimulus emmean SE df lower.CL upper.CL
no_pain 1.30 0.146 29 0.997 1.59
pain 2.01 0.179 29 1.645 2.38
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain -0.716 0.0922 29 -7.768 <.0001
Degrees-of-freedom method: kenward-roger 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.
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.
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.