Code
cat('\014') # clean terminalCode
rm(list = ls()) # clean workspace
library(tidyverse)
library(afex)
library(emmeans)
library(easystats)ERPs Pain Empathy 2024, MoBI device, all stimuli
Recordings with EMOTIV EPOC Flex Gel, 32 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(66, 71 ,72, 73)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')) |>
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
8 22 30 Code
summary(answers_df) ID Version Sex Block trialN
piloto62: 113 v1:1417 female: 821 1:1504 Min. : 1.00
piloto64: 112 v2:1567 male :2163 2:1480 1st Qu.:16.00
Piloto75: 112 Median :32.00
pilot52 : 110 Mean :31.74
pilot57 : 108 3rd Qu.:47.00
Piloto70: 108 Max. :64.00
(Other) :2321
Image Valence Answer Question rT
67.1.jpg: 18 EASE:1484 Min. : 1.0 displeasing:1490 Min. : 220
19.2.jpg: 17 PAIN:1500 1st Qu.: 1.0 painful :1494 1st Qu.: 3132
27.1.jpg: 17 Median :41.0 Median : 3864
33.1.jpg: 17 Mean :39.5 Mean : 4678
51.2.jpg: 17 3rd Qu.:68.0 3rd Qu.: 5245
hnp-4 : 17 Max. :97.0 Max. :37608
(Other) :2881
Session num_id Stimulus
Min. :1 Min. :51.00 no_pain:1484
1st Qu.:1 1st Qu.:58.00 pain :1500
Median :1 Median :65.00
Mean :1 Mean :67.37
3rd Qu.:1 3rd Qu.:77.00
Max. :1 Max. :88.00
Subject
Piloto75_v2_2024-07-12_19-17-27: 66
piloto62_v2_2024-05-17_20-26-13: 65
piloto64_v2_2024-05-24_17-09-35: 64
Piloto74_v2_2024-07-12_17-44-31: 62
pilot57_v1_2024-01-12_17-00-54 : 61
Piloto68_v2_2024-07-10_17-08-26: 61
(Other) :2605 Pain rating:
Code
painrat_data <- subset(answers_df, Question == 'painful')
pain_rating_lmer <- lmer(Answer ~ Stimulus + (Stimulus|num_id), painrat_data)
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 130040 130040 1 29.41 409.11 < 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.88, 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 13.9 1.76 29 10.3 17.5
pain 68.3 1.77 29 64.7 71.9
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain -54.4 2.69 29 -20.225 <.0001
Degrees-of-freedom method: kenward-roger Unpleasantness rating:
Code
unplsnt_data <- subset(answers_df, Question == 'displeasing')
unpleasantness_rating_lmer <- lmer(Answer ~ Stimulus + (Stimulus|num_id), unplsnt_data)
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 58547 58547 1 28.557 182.31 6.337e-14 ***
---
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.86 | [0.77, 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.10 29 12.3 20.9
pain 58.8 2.61 29 53.4 64.1
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain -42.2 3.12 29 -13.502 <.0001
Degrees-of-freedom method: kenward-roger ERP plots
EDF 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
8 22 30 Code
summary(painEmpathy_data) worklat Component Amplitude chindex Electrode
[130 160]:180 N1 :180 Min. :-7.5999 2 :120 FC1 :120
[145 165]:180 N2 :180 1st Qu.:-1.7716 6 :120 Fz :120
[240 310]:180 P1 :180 Median : 0.5714 16 :120 FC2 :120
[260 360]:180 P3 :180 Mean : 0.8813 18 :120 O1 :120
[600 800]:180 LPP:180 3rd Qu.: 3.1963 19 :120 Oz :120
Max. :16.1020 29 :120 O2 :120
(Other):180 (Other):180
bini Stimulus ERPset num_id Sex
5:450 no_pain:450 pilot51: 30 Min. :51.0 female:240
6:450 pain :450 pilot52: 30 1st Qu.:58.0 male :660
pilot53: 30 Median :66.0
pilot54: 30 Mean :67.5
pilot55: 30 3rd Qu.:77.0
pilot56: 30 Max. :88.0
(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 [145 165] 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.15121 0.15121 1 29 0.72 0.4031Code
interpret(omega_squared(painEmpathy_n1_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_n1_lmer)N2, average from [240 310] 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.02837 0.02837 1 28.998 0.1387 0.7123Code
interpret(omega_squared(painEmpathy_n2_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_n2_lmer)In the back: P1 & P3:
P1, average [130 160] 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.048256 0.048256 1 29 0.2638 0.6114Code
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 [260 360] 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.61753 0.61753 1 29 3.8292 0.06006 .
---
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.08 | [0.00, 1.00] | medium
- One-sided CIs: upper bound fixed at [1.00].
- Interpretation rule: field2013Code
a_posteriori_lmer(painEmpathy_p3_lmer)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 6.531 6.531 1 29 25.396 2.273e-05 ***
---
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.44 | [0.21, 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 0.71 0.107 29 0.49 0.929
pain 1.49 0.158 29 1.17 1.810
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain -0.778 0.154 29 -5.039 <.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.