Code
cat('\014') # clean terminalCode
rm(list = ls()) # clean workspace
library(tidyverse)
library(readxl)
library(afex)
library(emmeans)
library(effectsize)
library(GGally)HEP Intero Pernía 2026
Frontal Channles, Average Reference
Á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'))
}
}
}
a_posteriori_glmer <- function(glmer_obj, sig_level = .05) {
anova_glmer <- car::Anova(glmer_obj)
factors <- as.list(row.names(anova_glmer))
for (j in 1:length(factors)) {
if (grepl(':', factors[[j]])) {
factors[[j]] <- unlist(strsplit(factors[[j]], ':'))
}
}
p_values <- anova_glmer$`Pr(>Chisq)`
print(anova_glmer)
for (i in 1:length(p_values)) {
if (p_values[i] <= sig_level) {
cat(rep('_', 60), '\n', sep = '')
print(emmeans(glmer_obj, factors[[i]], contr = 'pairwise'))
}
}
}Code
struct_df <- read_excel('../heart_accuracy/IAcc_rho/ESTRATEGIAS FENO DE INTERO CARDIACA edit.xlsx') |>
rename(Structure = Estructura) |>
mutate(num_id = parse_number(ID),
num_id = as.factor(num_id))
sex_df <- read_csv('../heart_accuracy/IAcc_rho/rho_wide.csv', show_col_types = FALSE) |>
rename(Sex = Sexo) |>
mutate(num_id = as.factor(Participante),
Sex = if_else(str_detect(Sex, 'F'), 'female', 'male')) |>
mutate_if(is.character, as.factor)
hep_name <- 'avRef_hep_160to310_intero_pernia_2026.txt'
hep_data <- read.table(hep_name, header = TRUE, strip.white = TRUE, sep = "\t") |>
rename(Component = mlabel,
Amplitude = value,
Block = binlabel,
Electrode = chlabel) |>
mutate(Electrode = recode(Electrode, `E089-F1a` = 'F1', `E085-Fz` = 'Fz', `E076-F2a` = 'F2',
`E088-FC1a` = 'FC1', `E087-FCz` = 'FCz', `E075-FC2a` = 'FC2'),
Electrode = factor(Electrode, levels = c('F1', 'Fz', 'F2', 'FC1', 'FCz', 'FC2')),
Block = factor(Block, levels = c('regular_beat', 'irregular_beat', 'intero_pre', 'neck_pulse', 'intero_post')),
id = str_remove(ERPset, 'avRef_'),
id = str_remove(id , '_txal_hep'),
num_id = parse_number(id),
num_id = as.factor(num_id),
Group = if_else(str_detect(id, 'N'), 'non_obese', 'obese')) |>
left_join(struct_df[c('Structure', 'num_id')], by = 'num_id') |>
left_join(sex_df[c('Sex', 'num_id')], by = 'num_id') |>
mutate_if(is.character, as.factor)
write.csv(hep_data, 'avRef_hep_160to310_intero_pernia_2026_clean.csv', row.names = FALSE)HEP data general description
Code
summary(hep_data) worklat Component Amplitude chindex Electrode
[160.0 310.0]:2160 HEP:2160 Min. :-7.6537 Min. :75.00 F1 :360
1st Qu.:-1.2104 1st Qu.:76.00 Fz :360
Median :-0.7064 Median :86.00 F2 :360
Mean :-0.7913 Mean :83.33 FC1:360
3rd Qu.:-0.2857 3rd Qu.:88.00 FCz:360
Max. : 4.6160 Max. :89.00 FC2:360
bini Block ERPset id
Min. :1 regular_beat :432 avRef_N101_txal_hep: 30 N101 : 30
1st Qu.:2 irregular_beat:432 avRef_N102_txal_hep: 30 N102 : 30
Median :3 intero_pre :432 avRef_N103_txal_hep: 30 N103 : 30
Mean :3 neck_pulse :432 avRef_N104_txal_hep: 30 N104 : 30
3rd Qu.:4 intero_post :432 avRef_N105_txal_hep: 30 N105 : 30
Max. :5 avRef_N106_txal_hep: 30 N106 : 30
(Other) :1980 (Other):1980
num_id Group Structure Sex
101 : 30 non_obese:1020 extension :720 female:1470
102 : 30 obese :1140 imaginativa:900 male : 690
103 : 30 pecho :510
104 : 30 NAs : 30
105 : 30
106 : 30
(Other):1980 
HEP by Group
Code
options(width = 120)
rosner_hep <- EnvStats::rosnerTest(hep_data$Amplitude, 20)
hep_outliers <- rosner_hep$all.stats$Obs.Num[rosner_hep$all.stats$Outlier == TRUE]
EnvStats::print(rosner_hep)
Results of Outlier Test
-------------------------
Test Method: Rosner's Test for Outliers
Hypothesized Distribution: Normal
Data: hep_data$Amplitude
Sample Size: 2160
Test Statistics: R.1 = 7.952102
R.2 = 6.354710
R.3 = 6.065021
R.4 = 6.103766
R.5 = 6.091766
R.6 = 5.789777
R.7 = 5.737571
R.8 = 5.660630
R.9 = 5.359307
R.10 = 5.089175
R.11 = 4.896714
R.12 = 4.759464
R.13 = 4.468194
R.14 = 4.479133
R.15 = 4.245291
R.16 = 4.257305
R.17 = 4.223327
R.18 = 4.091434
R.19 = 3.987319
R.20 = 3.904166
Test Statistic Parameter: k = 20
Alternative Hypothesis: Up to 20 observations are not
from the same Distribution.
Type I Error: 5%
Number of Outliers Detected: 17
i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
1 0 -0.7913172 0.8629591 -7.653656 1854 7.952102 4.223858 TRUE
2 1 -0.7881387 0.8504184 4.616024 638 6.354710 4.223751 TRUE
3 2 -0.7906429 0.8426155 -5.901124 843 6.065021 4.223643 TRUE
4 3 -0.7882737 0.8355902 4.311973 639 6.103766 4.223536 TRUE
5 4 -0.7906393 0.8285280 -5.837838 647 6.091766 4.223428 TRUE
6 5 -0.7882972 0.8215506 -5.544892 2138 5.789777 4.223321 TRUE
7 6 -0.7860889 0.8153191 3.891862 642 5.737571 4.223213 TRUE
8 7 -0.7882617 0.8092469 -5.369109 646 5.660630 4.223106 TRUE
9 8 -0.7861330 0.8033835 3.519446 641 5.359307 4.222998 TRUE
10 9 -0.7881347 0.7981847 3.273967 640 5.089175 4.222890 TRUE
11 10 -0.7900241 0.7935449 3.095738 637 4.896714 4.222782 TRUE
12 11 -0.7918322 0.7892870 -4.548415 648 4.759464 4.222674 TRUE
13 12 -0.7900834 0.7852950 -4.298934 2132 4.468194 4.222567 TRUE
14 13 -0.7884491 0.7818156 -4.290305 2144 4.479133 4.222459 TRUE
15 14 -0.7868172 0.7783322 -4.091064 2059 4.245291 4.222351 TRUE
16 15 -0.7852768 0.7752347 -4.085687 2131 4.257305 4.222242 TRUE
17 16 -0.7837374 0.7721295 -4.044693 645 4.223327 4.222134 TRUE
18 17 -0.7822158 0.7690874 -3.928886 644 4.091434 4.222026 FALSE
19 18 -0.7807467 0.7662538 -3.836045 845 3.987319 4.221918 FALSE
20 19 -0.7793197 0.7635805 2.201825 2140 3.904166 4.221810 FALSECode
if (length(hep_outliers) != 0) {
hep_data_trim <- hep_data[-hep_outliers, ]
} else {
hep_data_trim <- hep_data
}
hepg_lmer <- lmer(Amplitude ~ Group * Sex * Block + (Block | num_id) + (1 | Electrode), hep_data_trim)
afex_plot(
hepg_lmer,
x = 'Block',
trace = 'Group',
panel = 'Sex',
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)
)
hepg_lmer@calllmer(formula = Amplitude ~ Group * Sex * Block + (Block | num_id) +
(1 | Electrode), data = hep_data_trim)
Code
anova(hepg_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Group 0.05177 0.05177 1 68.009 0.1730 0.67879
Sex 0.00029 0.00029 1 68.009 0.0010 0.97518
Block 1.36807 0.34202 4 63.566 1.1428 0.34459
Group:Sex 0.56375 0.56375 1 68.009 1.8837 0.17443
Group:Block 2.81699 0.70425 4 63.566 2.3531 0.06331 .
Sex:Block 1.40811 0.35203 4 63.566 1.1762 0.32977
Group:Sex:Block 0.55601 0.13900 4 63.566 0.4645 0.76154
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(hepg_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
------------------------------------------------------------------
Group | 0.00 | [0.00, 0.00] | very small
Sex | 0.00 | [0.00, 0.00] | very small
Block | 8.26e-03 | [0.00, 0.01] | very small
Group:Sex | 0.01 | [0.00, 0.11] | small
Group:Block | 0.07 | [0.00, 0.18] | medium
Sex:Block | 0.01 | [0.00, 0.02] | small
Group:Sex:Block | 0.00 | [0.00, 0.00] | very small
- Interpretation rule: field2013Code
a_posteriori_lmer(hepg_lmer)HEP by Structure, Intero Pre
Code
heps_lmer <- lmer(Amplitude ~ Group * Sex * Structure + (1 | num_id) + (1 | Electrode), subset(hep_data_trim, Block == 'intero_pre'))
afex_plot(
heps_lmer,
x = 'Structure',
trace = 'Group',
panel = 'Sex',
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)
)
heps_lmer@calllmer(formula = Amplitude ~ Group * Sex * Structure + (1 | num_id) +
(1 | Electrode), data = subset(hep_data_trim, Block == "intero_pre"))
Code
anova(heps_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Group 0.16919 0.16919 1 57.136 0.5906 0.4454
Sex 0.16553 0.16553 1 57.136 0.5778 0.4503
Structure 0.05514 0.02757 2 57.174 0.0962 0.9084
Group:Sex 0.36579 0.36579 1 57.135 1.2769 0.2632
Group:Structure 0.29252 0.14626 2 57.174 0.5106 0.6029
Sex:Structure 0.60101 0.30050 2 57.174 1.0490 0.3569
Group:Sex:Structure 0.32310 0.16155 2 57.174 0.5639 0.5721Code
interpret(omega_squared(heps_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
----------------------------------------------------------------------
Group | 0.00 | [0.00, 0.00] | very small
Sex | 0.00 | [0.00, 0.00] | very small
Structure | 0.00 | [0.00, 0.00] | very small
Group:Sex | 4.66e-03 | [0.00, 0.10] | very small
Group:Structure | 0.00 | [0.00, 0.00] | very small
Sex:Structure | 1.63e-03 | [0.00, 0.02] | very small
Group:Sex:Structure | 0.00 | [0.00, 0.00] | very small
- Interpretation rule: field2013Code
a_posteriori_lmer(heps_lmer)