Code
cat('\014') # clean terminalCode
rm(list = ls()) # clean workspace
library(tidyverse)
library(afex)
library(lmerTest)
library(emmeans)
library(easystats)ERPs Pain Empathy 2024, MoBI device, 2 versions
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_2023_answers_clean_2_versions_pilot_emotiv.csv', row.names = FALSE)
erp_averages <- list.files(path = './data/ave_volt_2_versions', full.names = TRUE)
painEmpathy_data <- vroom::vroom(erp_averages, show_col_types = FALSE) |>
separate(ERPset, c('ID', 'Version'), sep = '_', remove = FALSE) |>
rename(Component = mlabel,
Amplitude = value,
Stimulus = binlabel,
Electrode = chlabel) |>
mutate(worklat = gsub('.0', '', worklat, fixed = TRUE),
num_id = parse_number(ID),
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(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_data_clean_2_versions_pilot_emotiv.csv'), row.names = FALSE)Answers
Code
addmargins(xtabs(~ Version + Sex, data = unique(answers_df[c('Version', 'Sex', 'num_id')])), 2) Sex
Version female male Sum
v1 8 22 30
v2 8 22 30Code
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*Version + (Stimulus*Version|num_id) + (1|Image), painrat_data)
afex_plot(
pain_rating_lmer,
x = 'Stimulus',
trace = 'Version',
id = 'num_id',
error_arg = list(width = .15),
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
options(width = 150)
summary(pain_rating_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Answer ~ Stimulus * Version + (Stimulus * Version | num_id) + (1 | Image)
Data: painrat_data
REML criterion at convergence: 12863
Scaled residuals:
Min 1Q Median 3Q Max
-4.6924 -0.4425 -0.0708 0.4202 4.3810
Random effects:
Groups Name Variance Std.Dev. Corr
Image (Intercept) 51.67 7.188
num_id (Intercept) 83.49 9.137
Stimuluspain 172.03 13.116 -0.78
Versionv2 20.07 4.480 -0.22 0.36
Stimuluspain:Versionv2 15.02 3.876 -0.86 0.35 0.00
Residual 259.85 16.120
Number of obs: 1494, groups: Image, 262; num_id, 30
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 15.179 2.081 42.856 7.294 4.93e-09 ***
Stimuluspain 53.244 2.968 42.425 17.939 < 2e-16 ***
Versionv2 -2.807 1.872 69.562 -1.499 0.138
Stimuluspain:Versionv2 2.258 2.460 167.786 0.918 0.360
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) Stmlsp Vrsnv2
Stimuluspan -0.752
Versionv2 -0.451 0.386
Stmlspn:Vr2 0.083 -0.314 -0.609Code
r2(pain_rating_lmer)# R2 for Mixed Models
Conditional R2: 0.772
Marginal R2: 0.651Code
icc(pain_rating_lmer)# Intraclass Correlation Coefficient
Adjusted ICC: 0.346
Unadjusted ICC: 0.121Code
suppressMessages(ranova(pain_rating_lmer, reduce.terms = FALSE))ANOVA-like table for random-effects: Single term deletions
Model:
Answer ~ Stimulus + Version + (Stimulus * Version | num_id) + (1 | Image) + Stimulus:Version
npar logLik AIC LRT Df Pr(>Chisq)
<none> 16 -6431.5 12895
(Stimulus * Version | num_id) 6 -6561.7 13136 260.45 10 < 2.2e-16 ***
(1 | Image) 15 -6467.6 12965 72.22 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
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 95678 95678 1 35.899 368.201 <2e-16 ***
Version 330 330 1 47.306 1.271 0.2653
Stimulus:Version 219 219 1 167.786 0.842 0.3601
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(pain_rating_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
-------------------------------------------------------------------
Stimulus | 0.91 | [0.84, 0.94] | large
Version | 5.47e-03 | [0.00, 0.11] | very small
Stimulus:Version | 0.00 | [0.00, 0.00] | very small
- Interpretation rule: field2013Code
a_posteriori_lmer(pain_rating_lmer)____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
Stimulus emmean SE df lower.CL upper.CL
no_pain 13.8 1.86 36.8 10.0 17.5
pain 68.1 1.83 36.9 64.4 71.9
Results are averaged over the levels of: Version
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.84 35.5 -19.179 <.0001
Results are averaged over the levels of: Version
Degrees-of-freedom method: kenward-roger Unpleasantness rating:
Code
unplsnt_data <- subset(answers_df, Question == 'displeasing')
unpleasantness_rating_lmer <- lmer(Answer ~ Stimulus*Version + (Stimulus*Version|num_id) + (1|Image), unplsnt_data)
afex_plot(
unpleasantness_rating_lmer,
x = 'Stimulus',
trace = 'Version',
id = 'num_id',
error_arg = list(width = .15),
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
options(width = 150)
summary(unpleasantness_rating_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Answer ~ Stimulus * Version + (Stimulus * Version | num_id) + (1 | Image)
Data: unplsnt_data
REML criterion at convergence: 12890.1
Scaled residuals:
Min 1Q Median 3Q Max
-3.8070 -0.5192 -0.0686 0.4868 3.5526
Random effects:
Groups Name Variance Std.Dev. Corr
Image (Intercept) 36.31 6.026
num_id (Intercept) 184.82 13.595
Stimuluspain 297.11 17.237 -0.63
Versionv2 51.32 7.163 -0.73 0.32
Stimuluspain:Versionv2 78.65 8.868 0.65 -0.32 -0.83
Residual 271.38 16.473
Number of obs: 1490, groups: Image, 264; num_id, 30
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 19.664 2.748 32.401 7.156 3.74e-08 ***
Stimuluspain 40.011 3.556 32.774 11.250 8.57e-13 ***
Versionv2 -5.932 2.061 47.152 -2.878 0.006 **
Stimuluspain:Versionv2 4.179 2.777 49.191 1.505 0.139
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) Stmlsp Vrsnv2
Stimuluspan -0.644
Versionv2 -0.655 0.365
Stmlspn:Vr2 0.518 -0.434 -0.749Code
r2(unpleasantness_rating_lmer)# R2 for Mixed Models
Conditional R2: 0.707
Marginal R2: 0.489Code
icc(unpleasantness_rating_lmer)# Intraclass Correlation Coefficient
Adjusted ICC: 0.427
Unadjusted ICC: 0.219Code
suppressMessages(ranova(unpleasantness_rating_lmer, reduce.terms = FALSE))ANOVA-like table for random-effects: Single term deletions
Model:
Answer ~ Stimulus + Version + (Stimulus * Version | num_id) + (1 | Image) + Stimulus:Version
npar logLik AIC LRT Df Pr(>Chisq)
<none> 16 -6445.0 12922
(Stimulus * Version | num_id) 6 -6683.8 13380 477.49 10 < 2.2e-16 ***
(1 | Image) 15 -6461.8 12954 33.56 1 6.917e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
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 46763 46763 1 31.638 172.3192 2.399e-14 ***
Version 2119 2119 1 45.515 7.8084 0.007585 **
Stimulus:Version 614 614 1 49.191 2.2643 0.138779
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(unpleasantness_rating_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
-------------------------------------------------------------------
Stimulus | 0.84 | [0.72, 0.89] | large
Version | 0.13 | [0.01, 0.31] | medium
Stimulus:Version | 0.02 | [0.00, 0.16] | small
- Interpretation rule: field2013Code
a_posteriori_lmer(unpleasantness_rating_lmer)____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
Stimulus emmean SE df lower.CL upper.CL
no_pain 16.7 2.22 32.5 12.2 21.2
pain 58.8 2.66 31.4 53.4 64.2
Results are averaged over the levels of: Version
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain -42.1 3.21 32.3 -13.122 <.0001
Results are averaged over the levels of: Version
Degrees-of-freedom method: kenward-roger
____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
Version emmean SE df lower.CL upper.CL
v1 39.7 2.10 32.6 35.4 43.9
v2 35.8 1.84 33.6 32.1 39.6
Results are averaged over the levels of: Stimulus
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
v1 - v2 3.84 1.38 46.4 2.784 0.0077
Results are averaged over the levels of: Stimulus
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.5 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
options(width = 90)
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(~ Version + Sex, data = unique(painEmpathy_data[c('Version', 'Sex', 'num_id')])), 2) Sex
Version female male Sum
v1 8 22 30
v2 8 22 30Code
summary(painEmpathy_data) worklat Component Amplitude chindex Electrode bini
[130 160]:360 N1 :360 Min. :-8.5324 2 :240 FC1 :240 1:900
[145 165]:360 N2 :360 1st Qu.:-1.9035 6 :240 Fz :240 2:900
[240 310]:360 P1 :360 Median : 0.6011 16 :240 FC2 :240
[260 360]:360 P3 :360 Mean : 0.8823 18 :240 O1 :240
[600 800]:360 LPP:360 3rd Qu.: 3.0665 19 :240 Oz :240
Max. :16.7597 29 :240 O2 :240
(Other):360 (Other):360
Stimulus ERPset ID Version num_id Sex
no_pain:900 pilot51_v1: 30 pilot51: 60 v1:900 Min. :51.0 female: 480
pain :900 pilot51_v2: 30 pilot52: 60 v2:900 1st Qu.:58.0 male :1320
pilot52_v1: 30 pilot53: 60 Median :66.0
pilot52_v2: 30 pilot54: 60 Mean :67.5
pilot53_v1: 30 pilot55: 60 3rd Qu.:77.0
pilot53_v2: 30 pilot56: 60 Max. :88.0
(Other) :1620 (Other):1440 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 * Version + (Stimulus * Version | num_id) + (1 | num_id : Electrode), n1_data)
afex_plot(
painEmpathy_n1_lmer,
x = 'Stimulus',
trace = 'Version',
id = 'num_id',
error_arg = list(width = .25),
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 = 150)
summary(painEmpathy_n1_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Amplitude ~ Stimulus * Version + (Stimulus * Version | num_id) + (1 | num_id:Electrode)
Data: n1_data
REML criterion at convergence: 1162.2
Scaled residuals:
Min 1Q Median 3Q Max
-3.5338 -0.4072 -0.0469 0.3860 2.8709
Random effects:
Groups Name Variance Std.Dev. Corr
num_id:Electrode (Intercept) 0.7451 0.8632
num_id (Intercept) 2.5004 1.5813
Stimuluspain 0.9471 0.9732 -0.21
Versionv2 1.0268 1.0133 -0.53 0.15
Stimuluspain:Versionv2 1.1906 1.0912 0.28 -0.60 -0.42
Residual 0.5597 0.7481
Number of obs: 360, groups: num_id:Electrode, 90; num_id, 30
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.96397 0.31280 28.99936 -6.279 7.44e-07 ***
Stimuluspain -0.11372 0.20978 28.99977 -0.542 0.592
Versionv2 -0.03102 0.21602 28.99972 -0.144 0.887
Stimuluspain:Versionv2 -0.07574 0.25409 28.99987 -0.298 0.768
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) Stmlsp Vrsnv2
Stimuluspan -0.256
Versionv2 -0.509 0.248
Stmlspn:Vr2 0.284 -0.629 -0.505Code
r2(painEmpathy_n1_lmer)# R2 for Mixed Models
Conditional R2: 0.849
Marginal R2: 0.002Code
icc(painEmpathy_n1_lmer)# Intraclass Correlation Coefficient
Adjusted ICC: 0.849
Unadjusted ICC: 0.847Code
suppressMessages(ranova(painEmpathy_n1_lmer,reduce.terms = FALSE))Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00481871 (tol = 0.002,
component 1)ANOVA-like table for random-effects: Single term deletions
Model:
Amplitude ~ Stimulus + Version + (Stimulus * Version | num_id) + (1 | num_id:Electrode) + Stimulus:Version
npar logLik AIC LRT Df Pr(>Chisq)
<none> 16 -581.12 1194.2
(Stimulus * Version | num_id) 6 -638.81 1289.6 115.378 10 < 2.2e-16 ***
(1 | num_id:Electrode) 15 -627.36 1284.7 92.471 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
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.48346 0.48346 1 29 0.8638 0.3603
Version 0.07575 0.07575 1 29 0.1353 0.7156
Stimulus:Version 0.04973 0.04973 1 29 0.0889 0.7678Code
interpret(omega_squared(painEmpathy_n1_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
-------------------------------------------------------------------
Stimulus | 0.00 | [0.00, 0.00] | very small
Version | 0.00 | [0.00, 0.00] | very small
Stimulus:Version | 0.00 | [0.00, 0.00] | very small
- Interpretation rule: field2013Code
a_posteriori_lmer(painEmpathy_n1_lmer)N1 model check
Code
check_model(painEmpathy_n1_lmer)
N2, average from [240 310] ms interval
Electrodes FC1, Fz, FC2
Code
painEmpathy_n2_lmer <- lmer(Amplitude ~ Stimulus * Version + (Stimulus * Version | num_id) + (1 | num_id : Electrode), n2_data)
afex_plot(
painEmpathy_n2_lmer,
x = 'Stimulus',
trace = 'Version',
id = 'num_id',
error_arg = list(width = .25),
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 = 150)
summary(painEmpathy_n2_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Amplitude ~ Stimulus * Version + (Stimulus * Version | num_id) + (1 | num_id:Electrode)
Data: n2_data
REML criterion at convergence: 1104.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.91477 -0.47243 0.00705 0.43947 2.99247
Random effects:
Groups Name Variance Std.Dev. Corr
num_id:Electrode (Intercept) 1.2308 1.1094
num_id (Intercept) 1.5747 1.2549
Stimuluspain 0.6415 0.8010 0.06
Versionv2 1.1349 1.0653 -0.46 0.43
Stimuluspain:Versionv2 0.8448 0.9191 0.05 -0.81 -0.58
Residual 0.4418 0.6647
Number of obs: 360, groups: num_id:Electrode, 90; num_id, 30
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -3.0909 0.2666 29.0006 -11.594 2.08e-12 ***
Stimuluspain 0.1989 0.1766 29.0009 1.126 0.2694
Versionv2 1.0200 0.2183 28.9999 4.673 6.29e-05 ***
Stimuluspain:Versionv2 -0.4440 0.2186 29.0005 -2.031 0.0515 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) Stmlsp Vrsnv2
Stimuluspan -0.061
Versionv2 -0.435 0.445
Stmlspn:Vr2 0.119 -0.768 -0.599Code
r2(painEmpathy_n2_lmer)# R2 for Mixed Models
Conditional R2: 0.877
Marginal R2: 0.048Code
icc(painEmpathy_n2_lmer)# Intraclass Correlation Coefficient
Adjusted ICC: 0.871
Unadjusted ICC: 0.830Code
suppressMessages(ranova(painEmpathy_n2_lmer, reduce.terms = FALSE))ANOVA-like table for random-effects: Single term deletions
Model:
Amplitude ~ Stimulus + Version + (Stimulus * Version | num_id) + (1 | num_id:Electrode) + Stimulus:Version
npar logLik AIC LRT Df Pr(>Chisq)
<none> 16 -552.15 1136.3
(Stimulus * Version | num_id) 6 -599.16 1210.3 94.023 10 8.502e-16 ***
(1 | num_id:Electrode) 15 -639.79 1309.6 175.284 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
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.0174 0.0174 1 29.001 0.0395 0.84391
Version 9.0749 9.0749 1 29.000 20.5403 9.29e-05 ***
Stimulus:Version 1.8220 1.8220 1 29.000 4.1239 0.05154 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(painEmpathy_n2_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
-------------------------------------------------------------------
Stimulus | 0.00 | [0.00, 0.00] | very small
Version | 0.39 | [0.12, 0.59] | large
Stimulus:Version | 0.09 | [0.00, 0.32] | medium
- Interpretation rule: field2013Code
a_posteriori_lmer(painEmpathy_n2_lmer)____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
Version emmean SE df lower.CL upper.CL
v1 -2.99 0.276 29 -3.56 -2.43
v2 -2.19 0.256 29 -2.72 -1.67
Results are averaged over the levels of: Stimulus
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
v1 - v2 -0.798 0.176 29 -4.532 0.0001
Results are averaged over the levels of: Stimulus
Degrees-of-freedom method: kenward-roger N2 model check
Code
check_model(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 * Version + (Stimulus * Version | num_id) + (1 | num_id : Electrode), p1_data)
afex_plot(
painEmpathy_p1_lmer,
x = 'Stimulus',
trace = 'Version',
id = 'num_id',
error_arg = list(width = .25),
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 = 150)
summary(painEmpathy_p1_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Amplitude ~ Stimulus * Version + (Stimulus * Version | num_id) + (1 | num_id:Electrode)
Data: p1_data
REML criterion at convergence: 1194.9
Scaled residuals:
Min 1Q Median 3Q Max
-2.3851 -0.4599 -0.0062 0.4584 2.5601
Random effects:
Groups Name Variance Std.Dev. Corr
num_id:Electrode (Intercept) 1.6717 1.2929
num_id (Intercept) 14.7963 3.8466
Stimuluspain 1.0553 1.0273 -0.10
Versionv2 1.9940 1.4121 -0.42 -0.08
Stimuluspain:Versionv2 1.0505 1.0250 -0.07 -0.71 -0.20
Residual 0.3945 0.6281
Number of obs: 360, groups: num_id:Electrode, 90; num_id, 30
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.38627 0.71845 29.01290 4.713 5.62e-05 ***
Stimuluspain 0.09051 0.20962 28.98821 0.432 0.669
Versionv2 0.25863 0.27428 29.00602 0.943 0.354
Stimuluspain:Versionv2 0.02277 0.22924 28.98599 0.099 0.922
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) Stmlsp Vrsnv2
Stimuluspan -0.119
Versionv2 -0.408 0.011
Stmlspn:Vr2 -0.027 -0.703 -0.289
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00207968 (tol = 0.002, component 1)Code
r2(painEmpathy_p1_lmer)# R2 for Mixed Models
Conditional R2: 0.974
Marginal R2: 0.001Code
icc(painEmpathy_p1_lmer)# Intraclass Correlation Coefficient
Adjusted ICC: 0.974
Unadjusted ICC: 0.973Code
suppressMessages(ranova(painEmpathy_p1_lmer, reduce.terms = FALSE))ANOVA-like table for random-effects: Single term deletions
Model:
Amplitude ~ Stimulus + Version + (Stimulus * Version | num_id) + (1 | num_id:Electrode) + Stimulus:Version
npar logLik AIC LRT Df Pr(>Chisq)
<none> 16 -597.47 1226.9
(Stimulus * Version | num_id) 6 -729.05 1470.1 263.16 10 < 2.2e-16 ***
(1 | num_id:Electrode) 15 -712.23 1454.5 229.52 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
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.17588 0.17588 1 28.995 0.4459 0.5096
Version 0.40979 0.40979 1 29.004 1.0389 0.3165
Stimulus:Version 0.00389 0.00389 1 28.986 0.0099 0.9216Code
interpret(omega_squared(painEmpathy_p1_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
-------------------------------------------------------------------
Stimulus | 0.00 | [0.00, 0.00] | very small
Version | 1.25e-03 | [0.00, 0.11] | very small
Stimulus:Version | 0.00 | [0.00, 0.00] | very small
- Interpretation rule: field2013Code
a_posteriori_lmer(painEmpathy_p1_lmer)P1 model check
Code
check_model(painEmpathy_p1_lmer)
P3, average from [260 360] ms interval
Electrodes O1, Oz, O2
Code
painEmpathy_p3_lmer <- lmer(Amplitude ~ Stimulus * Version + (Stimulus * Version | num_id) + (1 | num_id : Electrode), p3_data)
afex_plot(
painEmpathy_p3_lmer,
x = 'Stimulus',
trace = 'Version',
id = 'num_id',
error_arg = list(width = .25),
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 = 150)
summary(painEmpathy_p3_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Amplitude ~ Stimulus * Version + (Stimulus * Version | num_id) + (1 | num_id:Electrode)
Data: p3_data
REML criterion at convergence: 1093.3
Scaled residuals:
Min 1Q Median 3Q Max
-3.8948 -0.4237 -0.0040 0.4404 2.3825
Random effects:
Groups Name Variance Std.Dev. Corr
num_id:Electrode (Intercept) 1.0919 1.0449
num_id (Intercept) 4.5871 2.1418
Stimuluspain 1.9172 1.3846 -0.27
Versionv2 1.5780 1.2562 -0.32 0.37
Stimuluspain:Versionv2 1.7037 1.3053 0.32 -0.82 -0.55
Residual 0.2998 0.5475
Number of obs: 360, groups: num_id:Electrode, 90; num_id, 30
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.0809 0.4103 29.0004 12.382 4.21e-13 ***
Stimuluspain -0.3657 0.2656 29.0000 -1.377 0.17919
Versionv2 -1.0159 0.2434 29.0001 -4.173 0.00025 ***
Stimuluspain:Versionv2 0.0423 0.2648 28.9998 0.160 0.87418
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) Stmlsp Vrsnv2
Stimuluspan -0.274
Versionv2 -0.317 0.380
Stmlspn:Vr2 0.307 -0.795 -0.570Code
r2(painEmpathy_p3_lmer)# R2 for Mixed Models
Conditional R2: 0.953
Marginal R2: 0.044Code
icc(painEmpathy_p3_lmer)# Intraclass Correlation Coefficient
Adjusted ICC: 0.951
Unadjusted ICC: 0.909Code
suppressMessages(ranova(painEmpathy_p3_lmer, reduce.terms = FALSE))ANOVA-like table for random-effects: Single term deletions
Model:
Amplitude ~ Stimulus + Version + (Stimulus * Version | num_id) + (1 | num_id:Electrode) + Stimulus:Version
npar logLik AIC LRT Df Pr(>Chisq)
<none> 16 -546.65 1125.3
(Stimulus * Version | num_id) 6 -659.20 1330.4 225.11 10 < 2.2e-16 ***
(1 | num_id:Electrode) 15 -648.99 1328.0 204.69 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
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 1.1047 1.1047 1 29 3.6851 0.06479 .
Version 7.4050 7.4050 1 29 24.7023 2.756e-05 ***
Stimulus:Version 0.0077 0.0077 1 29 0.0255 0.87418
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(painEmpathy_p3_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
-------------------------------------------------------------------
Stimulus | 0.08 | [0.00, 0.31] | medium
Version | 0.43 | [0.16, 0.63] | large
Stimulus:Version | 0.00 | [0.00, 0.00] | very small
- Interpretation rule: field2013Code
a_posteriori_lmer(painEmpathy_p3_lmer)____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
Version emmean SE df lower.CL upper.CL
v1 4.9 0.395 29 4.09 5.71
v2 3.9 0.404 29 3.08 4.73
Results are averaged over the levels of: Stimulus
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
v1 - v2 0.995 0.2 29 4.970 <.0001
Results are averaged over the levels of: Stimulus
Degrees-of-freedom method: kenward-roger P3 model check
Code
check_model(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 * Version + (Stimulus * Version | num_id) + (1 | num_id : Electrode), lpp_data)
afex_plot(
painEmpathy_lpp_lmer,
x = 'Stimulus',
trace = 'Version',
id = 'num_id',
error_arg = list(width = .25),
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 = 150)
summary(painEmpathy_lpp_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Amplitude ~ Stimulus * Version + (Stimulus * Version | num_id) + (1 | num_id:Electrode)
Data: lpp_data
REML criterion at convergence: 951.3
Scaled residuals:
Min 1Q Median 3Q Max
-4.5678 -0.4040 -0.0218 0.4855 2.6483
Random effects:
Groups Name Variance Std.Dev. Corr
num_id:Electrode (Intercept) 0.07779 0.2789
num_id (Intercept) 0.42671 0.6532
Stimuluspain 0.81297 0.9017 -0.48
Versionv2 0.46069 0.6787 -0.68 0.35
Stimuluspain:Versionv2 1.78885 1.3375 0.63 -0.60 -0.78
Residual 0.46976 0.6854
Number of obs: 360, groups: num_id:Electrode, 90; num_id, 30
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.7999 0.1425 28.9990 5.613 4.63e-06 ***
Stimuluspain 0.6460 0.1937 29.0013 3.334 0.00235 **
Versionv2 -0.1790 0.1606 28.9951 -1.114 0.27426
Stimuluspain:Versionv2 0.2343 0.2837 28.9945 0.826 0.41560
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) Stmlsp Vrsnv2
Stimuluspan -0.529
Versionv2 -0.665 0.399
Stmlspn:Vr2 0.580 -0.631 -0.749Code
r2(painEmpathy_lpp_lmer)# R2 for Mixed Models
Conditional R2: 0.636
Marginal R2: 0.117Code
icc(painEmpathy_lpp_lmer)# Intraclass Correlation Coefficient
Adjusted ICC: 0.588
Unadjusted ICC: 0.519Code
suppressMessages(ranova(painEmpathy_lpp_lmer, reduce.terms = FALSE))ANOVA-like table for random-effects: Single term deletions
Model:
Amplitude ~ Stimulus + Version + (Stimulus * Version | num_id) + (1 | num_id:Electrode) + Stimulus:Version
npar logLik AIC LRT Df Pr(>Chisq)
<none> 16 -475.63 983.26
(Stimulus * Version | num_id) 6 -522.52 1057.04 93.785 10 9.482e-16 ***
(1 | num_id:Electrode) 15 -478.77 987.54 6.281 1 0.01221 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
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 11.9154 11.9154 1 29.002 25.3651 2.292e-05 ***
Version 0.1520 0.1520 1 29.000 0.3235 0.5739
Stimulus:Version 0.3204 0.3204 1 28.994 0.6821 0.4156
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(painEmpathy_lpp_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
-------------------------------------------------------------------
Stimulus | 0.44 | [0.17, 0.63] | large
Version | 0.00 | [0.00, 0.00] | very small
Stimulus:Version | 0.00 | [0.00, 0.00] | very small
- Interpretation rule: field2013Code
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 0.71 0.107 29 0.491 0.93
pain 1.47 0.155 29 1.156 1.79
Results are averaged over the levels of: Version
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain -0.763 0.152 29 -5.036 <.0001
Results are averaged over the levels of: Version
Degrees-of-freedom method: kenward-roger LPP model check
Code
check_model(painEmpathy_lpp_lmer)
Factorial Mass Univariate analysis (FMUT, cluster mass)
Fields and Kuperberg (2020)
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): 26.560450/169.891840 Median (semi-IQR) distance between all pairs of channels (in chanlocs units): 119.172654 (30.583399) Mean (SD) # of neighbors per channel: 2.6 (1.2) Median (semi-IQR) # of neighbors per channel: 2.0 (1.0) Min/max # of neighbors per channel: 0 to 4
Effects over time
painEmpathy_fmut_2_versions_pilot_emotiv: 0 significant StimulusXVersion cluster(s) out of 43 1 significant Stimulus cluster(s) out of 35 cluster 1 F-masss: 4245 cluster 1 p-value: 0.0001 2 significant Version cluster(s) out of 38 cluster 1 F-masss: 1581 cluster 1 p-value: 0.0009 cluster 2 F-masss: 1665 cluster 2 p-value: 0.0006
Mass Univariate analysis (cluster-mass), full time window
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): 26.560450/169.891840 Median (semi-IQR) distance between all pairs of channels (in chanlocs units): 119.172654 (30.583399) Mean (SD) # of neighbors per channel: 2.6 (1.2) Median (semi-IQR) # of neighbors per channel: 2.0 (1.0) Min/max # of neighbors per channel: 0 to 4
Groppe, Urbach, and Kutas (2011)
painEmpathy_v1_mass_pilot_emotiv: 1 significant positive cluster(s) out of 11 cluster 1 t-masss: 567 cluster 1 p-value: 0.0000 0 significant negative cluster(s) out of 35 painEmpathy_v2_mass_pilot_emotiv: 1 significant positive cluster(s) out of 16 cluster 1 t-masss: 797 cluster 1 p-value: 0.0002 0 significant negative cluster(s) out of 13
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).
Using 13 central electrodes shared between BioSemi ActiveTwo and Emotiv Epoc Flex devices:
- F3, Fz, F4
- FC1, FC2
- C3, Cz, C4
- CP1 ,CP2
- P3, Pz, P4
Done in the ADAM 1.14 toolbox (Fahrenfort 2020) with FieldTrip 20211209 (Oostenveld et al. 2011).
Exemplar ERPs on Pz
Decoding over time
Decoding difference
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.
Fields, Eric C., and Gina R. Kuperberg. 2020. “Having Your Cake and Eating It Too: Flexibility and Power with Mass Univariate Statistics for ERP Data.” Psychophysiology 57 (February): e13468. https://doi.org/10.1111/PSYP.13468.
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.