Code
cat('\014') # clean terminalCode
rm(list = ls()) # clean workspace
library(tidyverse)
library(afex)
library(lmerTest)
library(emmeans)
library(easystats)
library(optimx)Ratings Pain Empathy 2023, 2 devices, 2 versions
Behavioural Response LMMs
Álvaro Rivera-Rei 
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'))
}
}
}Code
answers_df_biosemi <- read_csv('painEmpathy_2023_answers_clean_2_versions_pilot_biosemi.csv', show_col_types = FALSE) |>
mutate(System = 'BioSemi')
answers_df_emotiv <- read_csv('painEmpathy_2023_answers_clean_2_versions_pilot_emotiv.csv', show_col_types = FALSE) |>
mutate(System = 'EMOTIV')
answers_df <- rbind(answers_df_biosemi, answers_df_emotiv) |>
mutate(Session = factor(Session)) |>
mutate(num_id = factor(num_id)) |>
mutate_if(is.character, as.factor)
write.csv(answers_df, 'painEmpathy_2023_answers_biosemi_emotiv.csv', row.names = FALSE)Answers Summary
Code
addmargins(xtabs(~ Sex + System, data = unique(answers_df[c('System', 'Sex', 'num_id')]))) System
Sex BioSemi EMOTIV Sum
female 19 8 27
male 11 22 33
Sum 30 30 60Code
summary(answers_df) ID Version Sex Block trialN
pilot07 : 121 v1:2842 female:2716 Min. :1.000 Min. : 1.0
piloto62: 113 v2:3062 male :3188 1st Qu.:1.000 1st Qu.:16.0
piloto64: 112 Median :1.000 Median :32.0
Piloto75: 112 Mean :1.499 Mean :31.8
pilot52 : 110 3rd Qu.:2.000 3rd Qu.:47.0
pilot57 : 108 Max. :2.000 Max. :64.0
(Other) :5228
Image Valence Answer Question rT
hp-27 : 33 EASE:2918 Min. : 1.00 displeasing:2969 Min. : 220
2.2.jpg : 31 PAIN:2986 1st Qu.: 1.00 painful :2935 1st Qu.: 3356
51.2.jpg: 30 Median :48.00 Median : 4224
fp-21 : 30 Mean :41.31 Mean : 4834
fp-8 : 30 3rd Qu.:72.00 3rd Qu.: 5527
hnp-28 : 30 Max. :97.00 Max. :94218
(Other) :5720
Session num_id Stimulus Subject
1:5804 7 : 121 no_pain:2918 Piloto75_v2_2024-07-12_19-17-27: 66
2: 41 62 : 113 pain :2986 piloto62_v2_2024-05-17_20-26-13: 65
3: 59 64 : 112 piloto64_v2_2024-05-24_17-09-35: 64
75 : 112 pilot07_v1_2023-11-02_14-27-38 : 62
52 : 110 Piloto74_v2_2024-07-12_17-44-31: 62
57 : 108 pilot57_v1_2024-01-12_17-00-54 : 61
(Other):5228 (Other) :5524
System
BioSemi:2920
EMOTIV :2984
Pain Rating:
Code
painrat_data <- subset(answers_df, Question == 'painful')
pain_rating_lmer <- lmer(Answer ~ System*Stimulus*Version + (Stimulus*Version|num_id) + (1|Image), painrat_data,
control = lmerControl(optimizer = 'optimx', optCtrl = list(method = 'nlminb')))
afex_plot(
pain_rating_lmer,
x = ~Stimulus,
trace = ~Version,
panel = ~System,
id = 'num_id',
error_arg = list(width = .2),
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)
)
Pain Linear Mixed Model
Code
options(width = 120)
summary(pain_rating_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Answer ~ System * Stimulus * Version + (Stimulus * Version | num_id) + (1 | Image)
Data: painrat_data
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))
REML criterion at convergence: 25193.3
Scaled residuals:
Min 1Q Median 3Q Max
-5.1132 -0.4324 -0.0632 0.4190 4.4521
Random effects:
Groups Name Variance Std.Dev. Corr
Image (Intercept) 48.51 6.965
num_id (Intercept) 135.80 11.653
Stimuluspain 253.92 15.935 -0.75
Versionv2 36.08 6.007 -0.61 0.45
Stimuluspain:Versionv2 33.82 5.816 0.08 -0.39 -0.48
Residual 255.88 15.996
Number of obs: 2935, groups: Image, 312; num_id, 60
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 12.6201 2.4616 77.8500 5.127 2.1e-06 ***
SystemEMOTIV 2.5275 3.2643 60.8450 0.774 0.4417
Stimuluspain 63.4092 3.3908 76.5528 18.700 < 2e-16 ***
Versionv2 -2.4630 2.0068 103.5477 -1.227 0.2225
SystemEMOTIV:Stimuluspain -10.0487 4.4847 59.2824 -2.241 0.0288 *
SystemEMOTIV:Versionv2 -0.4553 2.3520 56.3436 -0.194 0.8472
Stimuluspain:Versionv2 2.3590 2.5985 118.6599 0.908 0.3658
SystemEMOTIV:Stimuluspain:Versionv2 -0.1471 2.9046 58.0856 -0.051 0.9598
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) SyEMOTIV Stmlsp Vrsnv2 SyEMOTIV:S SEMOTIV:V Stm:V2
SystmEMOTIV -0.666
Stimuluspan -0.738 0.493
Versionv2 -0.593 0.347 0.430
SysEMOTIV:S 0.494 -0.742 -0.664 -0.253
SyEMOTIV:V2 0.393 -0.579 -0.287 -0.598 0.423
Stmlspn:Vr2 0.262 -0.121 -0.472 -0.645 0.246 0.357
SEMOTIV:S:V -0.145 0.207 0.292 0.374 -0.428 -0.622 -0.580Pain LMM Random Effects
The marginal R2 considers only the variance of the fixed effects (without the random effects), while the conditional R2 takes both the fixed and random effects into account (i.e., the total model).
The adjusted ICC only relates to the random effects, the unadjusted ICC also takes the fixed effects variances into account, more precisely, the fixed effects variance is added to the denominator of the formula to calculate the ICC.
Code
r2(pain_rating_lmer)# R2 for Mixed Models
Conditional R2: 0.803
Marginal R2: 0.686Code
icc(pain_rating_lmer)# Intraclass Correlation Coefficient
Adjusted ICC: 0.375
Unadjusted ICC: 0.118Code
ranova(pain_rating_lmer, reduce.terms = FALSE)ANOVA-like table for random-effects: Single term deletions
Model:
Answer ~ System + Stimulus + Version + (Stimulus * Version | num_id) + (1 | Image) + System:Stimulus + System:Version + Stimulus:Version + System:Stimulus:Version
npar logLik AIC LRT Df Pr(>Chisq)
<none> 20 -12597 25233
(Stimulus * Version | num_id) 10 -12928 25877 663.57 10 < 2.2e-16 ***
(1 | Image) 19 -12693 25425 193.47 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pain LMM Fixed Effects
ANOVA and post-hoc tests.
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)
System 642 642 1 58.615 2.5074 0.1187
Stimulus 186400 186400 1 78.912 728.4791 <2e-16 ***
Version 419 419 1 125.791 1.6392 0.2028
System:Stimulus 1575 1575 1 58.669 6.1541 0.0160 *
System:Version 21 21 1 55.783 0.0825 0.7750
Stimulus:Version 298 298 1 149.853 1.1641 0.2823
System:Stimulus:Version 1 1 1 58.086 0.0026 0.9598
---
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
--------------------------------------------------------------------------
System | 0.02 | [0.00, 0.15] | small
Stimulus | 0.90 | [0.86, 0.93] | large
Version | 4.98e-03 | [0.00, 0.06] | very small
System:Stimulus | 0.08 | [0.00, 0.23] | medium
System:Version | 0.00 | [0.00, 0.00] | very small
Stimulus:Version | 1.08e-03 | [0.00, 0.03] | very small
System: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 12.5 1.50 80.1 9.55 15.5
pain 72.1 1.43 82.6 69.21 74.9
Results are averaged over the levels of: System, Version
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain -59.5 2.21 78.1 -26.982 <.0001
Results are averaged over the levels of: System, Version
Degrees-of-freedom method: kenward-roger
____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
System Stimulus emmean SE df lower.CL upper.CL
BioSemi no_pain 11.4 2.03 68.7 7.33 15.4
EMOTIV no_pain 13.7 2.04 69.3 9.62 17.8
BioSemi pain 76.0 1.94 70.5 72.11 79.8
EMOTIV pain 68.2 1.94 70.1 64.29 72.0
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
BioSemi no_pain - EMOTIV no_pain -2.30 2.76 58.0 -0.835 0.8378
BioSemi no_pain - BioSemi pain -64.59 3.00 68.0 -21.497 <.0001
BioSemi no_pain - EMOTIV pain -56.77 2.81 134.3 -20.203 <.0001
EMOTIV no_pain - BioSemi pain -62.29 2.82 135.4 -22.126 <.0001
EMOTIV no_pain - EMOTIV pain -54.47 3.01 68.1 -18.120 <.0001
BioSemi pain - EMOTIV pain 7.82 2.61 58.1 2.996 0.0204
Results are averaged over the levels of: Version
Degrees-of-freedom method: kenward-roger
P value adjustment: tukey method for comparing a family of 4 estimates Unpleasantness Rating:
Code
unplsnt_data <- subset(answers_df, Question == 'displeasing')
unpleasantness_rating_lmer <- lmer(Answer ~ System*Stimulus*Version + (Stimulus*Version|num_id) + (1|Image), unplsnt_data)
afex_plot(
unpleasantness_rating_lmer,
x = ~Stimulus,
trace = ~Version,
panel = ~System,
id = 'num_id',
error_arg = list(width = .2),
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)
)
Unpleasantness Linear Mixed Model
Code
options(width = 120)
summary(unpleasantness_rating_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Answer ~ System * Stimulus * Version + (Stimulus * Version | num_id) + (1 | Image)
Data: unplsnt_data
REML criterion at convergence: 25824.8
Scaled residuals:
Min 1Q Median 3Q Max
-4.5582 -0.5186 -0.0625 0.4714 3.6452
Random effects:
Groups Name Variance Std.Dev. Corr
Image (Intercept) 58.40 7.642
num_id (Intercept) 190.74 13.811
Stimuluspain 273.93 16.551 -0.66
Versionv2 54.39 7.375 -0.77 0.44
Stimuluspain:Versionv2 63.35 7.959 0.48 -0.35 -0.62
Residual 281.74 16.785
Number of obs: 2969, groups: Image, 316; num_id, 60
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 17.576 2.865 71.324 6.135 4.28e-08 ***
SystemEMOTIV 2.410 3.817 56.624 0.632 0.53025
Stimuluspain 49.024 3.561 76.134 13.767 < 2e-16 ***
Versionv2 -5.172 2.266 110.062 -2.282 0.02439 *
SystemEMOTIV:Stimuluspain -9.164 4.671 56.913 -1.962 0.05466 .
SystemEMOTIV:Versionv2 -1.214 2.652 58.100 -0.458 0.64875
Stimuluspain:Versionv2 7.747 2.930 124.968 2.644 0.00923 **
SystemEMOTIV:Stimuluspain:Versionv2 -3.300 3.314 61.001 -0.996 0.32328
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) SyEMOTIV Stmlsp Vrsnv2 SyEMOTIV:S SEMOTIV:V Stm:V2
SystmEMOTIV -0.672
Stimuluspan -0.671 0.441
Versionv2 -0.681 0.418 0.449
SysEMOTIV:S 0.448 -0.665 -0.660 -0.266
SyEMOTIV:V2 0.475 -0.697 -0.299 -0.604 0.441
Stmlspn:Vr2 0.424 -0.246 -0.480 -0.682 0.250 0.390
SEMOTIV:S:V -0.290 0.421 0.290 0.404 -0.426 -0.665 -0.580Unpleasantness LMM Random Effects
The marginal R2 considers only the variance of the fixed effects (without the random effects), while the conditional R2 takes both the fixed and random effects into account (i.e., the total model).
The adjusted ICC only relates to the random effects, the unadjusted ICC also takes the fixed effects variances into account, more precisely, the fixed effects variance is added to the denominator of the formula to calculate the ICC.
Code
r2(unpleasantness_rating_lmer)# R2 for Mixed Models
Conditional R2: 0.739
Marginal R2: 0.543Code
icc(unpleasantness_rating_lmer)# Intraclass Correlation Coefficient
Adjusted ICC: 0.429
Unadjusted ICC: 0.196Code
ranova(unpleasantness_rating_lmer, reduce.terms = FALSE)ANOVA-like table for random-effects: Single term deletions
Model:
Answer ~ System + Stimulus + Version + (Stimulus * Version | num_id) + (1 | Image) + System:Stimulus + System:Version + Stimulus:Version + System:Stimulus:Version
npar logLik AIC LRT Df Pr(>Chisq)
<none> 20 -12912 25865
(Stimulus * Version | num_id) 10 -13333 26686 840.65 10 < 2.2e-16 ***
(1 | Image) 19 -13025 26088 225.11 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Unpleasantness LMM Fixed Effects
ANOVA and post-hoc tests.
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)
System 626 626 1 58.441 2.2222 0.14141
Stimulus 119007 119007 1 80.132 422.4035 < 2e-16 ***
Version 1210 1210 1 134.302 4.2938 0.04016 *
System:Stimulus 1834 1834 1 57.615 6.5097 0.01341 *
System:Version 587 587 1 56.800 2.0836 0.15438
Stimulus:Version 1837 1837 1 165.046 6.5217 0.01156 *
System:Stimulus:Version 279 279 1 61.001 0.9916 0.32328
---
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
--------------------------------------------------------------------------
System | 0.02 | [0.00, 0.14] | small
Stimulus | 0.84 | [0.77, 0.88] | large
Version | 0.02 | [0.00, 0.10] | small
System:Stimulus | 0.08 | [0.00, 0.24] | medium
System:Version | 0.02 | [0.00, 0.14] | small
Stimulus:Version | 0.03 | [0.00, 0.10] | small
System:Stimulus:Version | 0.00 | [0.00, 0.00] | very 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 15.9 1.66 80.0 12.6 19.2
pain 63.4 1.81 75.7 59.8 67.0
Results are averaged over the levels of: System, Version
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain - pain -47.5 2.31 80.6 -20.545 <.0001
Results are averaged over the levels of: System, 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 41.0 1.57 83.1 37.9 44.1
v2 38.3 1.32 90.3 35.6 40.9
Results are averaged over the levels of: System, Stimulus
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
v1 - v2 2.73 1.32 136 2.066 0.0407
Results are averaged over the levels of: System, Stimulus
Degrees-of-freedom method: kenward-roger
____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
System Stimulus emmean SE df lower.CL upper.CL
BioSemi no_pain 15.0 2.25 69.0 10.5 19.5
EMOTIV no_pain 16.8 2.25 68.7 12.3 21.3
BioSemi pain 67.9 2.47 66.4 63.0 72.8
EMOTIV pain 58.9 2.48 67.2 53.9 63.8
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
BioSemi no_pain - EMOTIV no_pain -1.80 3.05 58.0 -0.592 0.9341
BioSemi no_pain - BioSemi pain -52.90 3.13 69.0 -16.880 <.0001
BioSemi no_pain - EMOTIV pain -43.89 3.35 132.0 -13.106 <.0001
EMOTIV no_pain - BioSemi pain -51.09 3.34 130.8 -15.292 <.0001
EMOTIV no_pain - EMOTIV pain -42.08 3.14 69.4 -13.409 <.0001
BioSemi pain - EMOTIV pain 9.01 3.38 58.1 2.669 0.0472
Results are averaged over the levels of: Version
Degrees-of-freedom method: kenward-roger
P value adjustment: tukey method for comparing a family of 4 estimates
____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
Stimulus Version emmean SE df lower.CL upper.CL
no_pain v1 18.8 2.12 86.1 14.56 23.0
pain v1 63.2 2.00 90.0 59.24 67.2
no_pain v2 13.0 1.62 105.0 9.79 16.2
pain v2 63.5 2.02 84.3 59.53 67.6
Results are averaged over the levels of: System
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
no_pain v1 - pain v1 -44.442 2.68 95.2 -16.602 <.0001
no_pain v1 - no_pain v2 5.779 1.81 143.2 3.188 0.0094
no_pain v1 - pain v2 -44.760 2.84 87.6 -15.736 <.0001
pain v1 - no_pain v2 50.221 2.47 100.6 20.353 <.0001
pain v1 - pain v2 -0.318 1.76 144.3 -0.181 0.9979
no_pain v2 - pain v2 -50.539 2.53 93.8 -19.987 <.0001
Results are averaged over the levels of: System
Degrees-of-freedom method: kenward-roger
P value adjustment: tukey method for comparing a family of 4 estimates