Code
cat('\014') # clean terminalCode
rm(list = ls())
library(tidyverse)
library(GGally)
library(afex)
library(emmeans)
library(lmerTest)
library(easystats)Robot 2023, Linear Mixed-Effects Model
Recordings with Emotiv EPOC Flex (gel)
Álvaro Rivera-Rei 
Code
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
robot_files <- c('erplab/average_voltage_100_to_150_robot_2023.txt',
'erplab/average_voltage_290_to_350_robot_2023.txt',
'erplab/average_voltage_600_to_800_robot_2023.txt')
robot_data <- read_tsv(robot_files, show_col_types = FALSE)
robot_data <- robot_data |>
rename(component = mlabel,
amplitude = value,
stimulus = binlabel,
electrode = chlabel) |>
separate(stimulus, c('valence', 'face'), sep = '_', remove = FALSE) |>
mutate(valence = factor(valence, levels = c('neutral', 'happy'))) |>
mutate(face = factor(face, levels = c('robot', 'human'))) |>
mutate(component = factor(component, levels = c('P1', 'P3', 'Late'))) |>
separate(worklat, c('start', 'end'), sep = ' ', remove = FALSE) |>
mutate(start = parse_number(start)) |>
mutate(end = parse_number(end)) |>
mutate(num_id = parse_number(ERPset)) |>
mutate_if(is.character, as.factor) |>
mutate(electrode = factor(electrode, levels = c('P3', 'Pz', 'P4', 'O1', 'Oz', 'O2'))) |>
mutate(antero_posterior = ifelse(electrode %in% c('P3', 'Pz', 'P4'), 'parietal', 'occipital')) |>
mutate(antero_posterior = factor(antero_posterior, levels = c('parietal', 'occipital'))) |>
mutate(medio_lateral = ifelse(electrode %in% c('P3', 'O1'), 'left', ifelse(electrode %in% c('Pz', 'Oz'), 'central', 'right'))) |>
mutate(medio_lateral = factor(medio_lateral, levels = c('left', 'central', 'right'))) |>
mutate_at(c('num_id', 'chindex', 'bini'), as.factor)
write_csv(robot_data, 'mrs_robot_2023_data_clean.csv')
n170_file <- 'erplab/average_voltage_150_to_200_robot_2023.txt'
n170_data <- read_tsv(n170_file, show_col_types = FALSE) |>
rename(component = mlabel,
amplitude = value,
stimulus = binlabel,
electrode = chlabel) |>
separate(stimulus, c('valence', 'face'), sep = '_', remove = FALSE) |>
mutate(valence = factor(valence, levels = c('neutral', 'happy'))) |>
mutate(face = factor(face, levels = c('robot', 'human'))) |>
separate(worklat, c('start', 'end'), sep = ' ', remove = FALSE) |>
mutate(start = parse_number(start)) |>
mutate(end = parse_number(end)) |>
mutate(num_id = parse_number(ERPset)) |>
mutate_if(is.character, as.factor) |>
mutate(electrode = factor(electrode, levels = c('P7', 'P8', 'PO9', 'PO10'))) |>
mutate(antero_posterior = ifelse(electrode %in% c('P7', 'P8'), 'parietal', 'parieto_occipital')) |>
mutate(antero_posterior = factor(antero_posterior, levels = c('parietal', 'parieto_occipital'))) |>
mutate(medio_lateral = ifelse(electrode %in% c('P7', 'PO9'), 'left', 'right')) |>
mutate(medio_lateral = factor(medio_lateral, levels = c('left', 'right'))) |>
mutate_at(c('num_id', 'chindex', 'bini'), as.factor)
write_csv(n170_data, 'mrs_robot_2023_n170_data_clean.csv')
connectivity_file <- 'erplab/median_connectivity.csv'
connectivity_data <- read_csv(connectivity_file, show_col_types = FALSE) |>
mutate(log_median_connectivity = log(median_connectivity)) |>
mutate(face = factor(face, levels = c('robot', 'human'))) |>
mutate(num_id = parse_number(ERPset)) |>
mutate_if(is.character, as.factor) |>
mutate(region_b = factor(region_b, levels = c('inf_temp_ipsi', 'inf_temp_contra')))
write_csv(connectivity_data, 'mrs_robot_2023_connectivity_data_clean.csv')
p1_data <- subset(robot_data, component == 'P1')
p3_data <- subset(robot_data, component == 'P3')
late_data <- subset(robot_data, component == 'Late')ERP plots
EDF files were pre-processed in Matlab using the EEGLAB 2024.0 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 30Hz 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.
Face & Valence ERPs:
Topographic layout:
General description
Code
summary(robot_data) worklat start end component amplitude
[100.0 150.0]:552 Min. :100 Min. :150.0 P1 :552 Min. :-5.184
[290.0 350.0]:552 1st Qu.:100 1st Qu.:150.0 P3 :552 1st Qu.: 1.218
[600.0 800.0]:552 Median :290 Median :350.0 Late:552 Median : 3.093
Mean :330 Mean :433.3 Mean : 3.540
3rd Qu.:600 3rd Qu.:800.0 3rd Qu.: 5.561
Max. :600 Max. :800.0 Max. :20.933
chindex electrode bini stimulus valence face
13:368 P3:368 1:414 happy_human :414 neutral:828 robot:828
16:184 Pz:368 2:414 happy_robot :414 happy :828 human:828
17:368 P4:368 3:414 neutral_human:414
18:184 O1:184 4:414 neutral_robot:414
19:184 Oz:184
22:368 O2:184
ERPset num_id antero_posterior medio_lateral
MRS_04_ct_FT9: 36 4 : 36 parietal :1104 left :552
MRS_07_ct_FT9: 36 7 : 36 occipital: 552 central:552
MRS_08_ct_FT9: 36 8 : 36 right :552
MRS_09_ct_FT9: 36 9 : 36
MRS_10_ct_FT9: 36 10 : 36
MRS_11_ct_FT9: 36 11 : 36
(Other) :1440 (Other):1440 Code
ggplot(
robot_data, aes(x = amplitude, fill = component, color = component)) +
geom_histogram(alpha = .4) +
facet_wrap(~component)
robot_data_wide <- robot_data[c('ERPset', 'face', 'valence', 'medio_lateral', 'component', 'amplitude')] |>
pivot_wider(names_from = component, values_from = amplitude)
components <- c('P1', 'P3', 'Late')
components_pairs <- ggpairs(robot_data_wide,
aes(colour = face, alpha = .5, linewidth = NA),
columns = components,
lower = list(continuous = wrap('points', alpha = .4)),
progress = FALSE
)
suppressWarnings(print(components_pairs))
P1 by Face & Valence, average 100 to 150 ms
Electrodes O1, Oz, O2
P1 Scalp voltage distribution:
P1 ANOVA:
Code
robot_p1_lmer <- lmer(amplitude ~ face * valence + (face * valence | num_id ) + (1 | num_id : electrode), p1_data)
afex_plot(robot_p1_lmer,
x = 'valence',
trace = 'face',
id = 'num_id',
error_arg = list(width = .15, lwd = .75),
dodge = .5,
data_arg = list(
position =
position_jitterdodge(
jitter.width = .2,
dodge.width = .5 ## needs to be same as dodge
)),
mapping = c('color'),, data_alpha = .3,
point_arg = list(size = 3)
)
Code
# options(width = 90)
summary(robot_p1_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: amplitude ~ face * valence + (face * valence | num_id) + (1 |
num_id:electrode)
Data: p1_data
REML criterion at convergence: 2080.5
Scaled residuals:
Min 1Q Median 3Q Max
-2.5119 -0.4296 -0.0151 0.4282 3.2060
Random effects:
Groups Name Variance Std.Dev. Corr
num_id:electrode (Intercept) 1.0279 1.0139
num_id (Intercept) 16.1091 4.0136
facehuman 6.7231 2.5929 -0.65
valencehappy 3.4432 1.8556 -0.32 0.60
facehuman:valencehappy 8.2721 2.8761 0.23 -0.68 -0.89
Residual 0.7334 0.8564
Number of obs: 552, groups: num_id:electrode, 138; num_id, 46
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.7307 0.6025 45.0194 9.512 2.42e-12 ***
facehuman -1.2218 0.3960 45.0013 -3.086 0.00347 **
valencehappy -0.2727 0.2924 44.9794 -0.933 0.35596
facehuman:valencehappy 0.2869 0.4484 44.9803 0.640 0.52549
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) fachmn vlnchp
facehuman -0.642
valencehppy -0.324 0.586
fchmn:vlnch 0.233 -0.681 -0.872Code
ranova(robot_p1_lmer, reduce.terms = FALSE)ANOVA-like table for random-effects: Single term deletions
Model:
amplitude ~ face + valence + (face * valence | num_id) + (1 | num_id:electrode) + face:valence
npar logLik AIC LRT Df Pr(>Chisq)
<none> 16 -1040.3 2112.5
(face * valence | num_id) 6 -1266.0 2544.0 451.50 10 < 2.2e-16 ***
(1 | num_id:electrode) 15 -1114.6 2259.2 148.71 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
anova(robot_p1_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
face 9.9032 9.9032 1 45.017 13.5026 0.0006313 ***
valence 0.5702 0.5702 1 44.991 0.7775 0.3825905
face:valence 0.3003 0.3003 1 44.980 0.4095 0.5254933
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(robot_p1_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
---------------------------------------------------------------
face | 0.21 | [0.04, 0.40] | large
valence | 0.00 | [0.00, 0.00] | very small
face:valence | 0.00 | [0.00, 0.00] | very small
- Interpretation rule: field2013Code
a_posteriori_lmer(robot_p1_lmer)____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
face emmean SE df lower.CL upper.CL
robot 5.59 0.572 45 4.44 6.75
human 4.52 0.425 45 3.66 5.37
Results are averaged over the levels of: valence
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
robot - human 1.08 0.293 45 3.675 0.0006
Results are averaged over the levels of: valence
Degrees-of-freedom method: kenward-roger P1 Assumptions:
Code
check_model(robot_p1_lmer)
P3 by Face & Valence, average 290 to 350 ms
Electrodes P3, Pz, P4
P3 Scalp voltage distribution:
P3 ANOVA:
Code
robot_p3_lmer <- lmer(amplitude ~ face * valence + (face * valence | num_id ) + (1 | num_id : electrode), p3_data)
afex_plot(robot_p3_lmer,
x = 'valence',
trace = 'face',
id = 'num_id',
error_arg = list(width = .15, lwd = .75),
dodge = .5,
data_arg = list(
position =
position_jitterdodge(
jitter.width = .2,
dodge.width = .5 ## needs to be same as dodge
)),
mapping = c('color'),, data_alpha = .3,
point_arg = list(size = 3)
)
Code
summary(robot_p3_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: amplitude ~ face * valence + (face * valence | num_id) + (1 |
num_id:electrode)
Data: p3_data
REML criterion at convergence: 2084.4
Scaled residuals:
Min 1Q Median 3Q Max
-2.68286 -0.47421 0.02464 0.46602 2.11002
Random effects:
Groups Name Variance Std.Dev. Corr
num_id:electrode (Intercept) 1.3136 1.1461
num_id (Intercept) 5.1775 2.2754
facehuman 2.9953 1.7307 -0.54
valencehappy 2.9473 1.7168 -0.46 0.45
facehuman:valencehappy 5.4521 2.3350 0.29 -0.36 -0.82
Residual 0.8202 0.9056
Number of obs: 552, groups: num_id:electrode, 138; num_id, 46
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 4.78787 0.35779 45.00022 13.382 < 2e-16 ***
facehuman -0.86703 0.27749 45.00046 -3.125 0.00311 **
valencehappy -0.48001 0.27561 44.99896 -1.742 0.08840 .
facehuman:valencehappy 0.02077 0.37722 44.99911 0.055 0.95633
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) fachmn vlnchp
facehuman -0.526
valencehppy -0.456 0.460
fchmn:vlnch 0.295 -0.414 -0.801Code
ranova(robot_p3_lmer, reduce.terms = FALSE)ANOVA-like table for random-effects: Single term deletions
Model:
amplitude ~ face + valence + (face * valence | num_id) + (1 | num_id:electrode) + face:valence
npar logLik AIC LRT Df Pr(>Chisq)
<none> 16 -1042.2 2116.4
(face * valence | num_id) 6 -1183.5 2378.9 282.49 10 < 2.2e-16 ***
(1 | num_id:electrode) 15 -1126.0 2282.1 167.64 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
anova(robot_p3_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
face 8.6896 8.6896 1 45.000 10.595 0.002157 **
valence 6.3949 6.3949 1 44.999 7.797 0.007656 **
face:valence 0.0025 0.0025 1 44.999 0.003 0.956328
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(robot_p3_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
---------------------------------------------------------------
face | 0.17 | [0.02, 0.36] | large
valence | 0.13 | [0.01, 0.32] | medium
face:valence | 0.00 | [0.00, 0.00] | very small
- Interpretation rule: field2013Code
a_posteriori_lmer(robot_p3_lmer)____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
face emmean SE df lower.CL upper.CL
robot 4.55 0.319 45 3.90 5.19
human 3.69 0.316 45 3.05 4.33
Results are averaged over the levels of: valence
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
robot - human 0.857 0.263 45 3.255 0.0022
Results are averaged over the levels of: valence
Degrees-of-freedom method: kenward-roger
____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
valence emmean SE df lower.CL upper.CL
neutral 4.35 0.308 45 3.73 4.98
happy 3.88 0.294 45 3.29 4.48
Results are averaged over the levels of: face
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
neutral - happy 0.47 0.168 45 2.792 0.0077
Results are averaged over the levels of: face
Degrees-of-freedom method: kenward-roger P3 Assumptions:
Code
check_model(robot_p3_lmer)
Late by Face & Valence, average 600 to 800 ms
Electrodes P3, Pz, P4
Late Scalp voltage distribution:
Late ANOVA:
Code
robot_late_lmer <- lmer(amplitude ~ face * valence + (face * valence | num_id ) + (1 | num_id : electrode), late_data)
afex_plot(robot_late_lmer,
x = 'valence',
trace = 'face',
id = 'num_id',
error_arg = list(width = .15, lwd = .75),
dodge = .5,
data_arg = list(
position =
position_jitterdodge(
jitter.width = .2,
dodge.width = .5 ## needs to be same as dodge
)),
mapping = c('color'),, data_alpha = .3,
point_arg = list(size = 3)
)
Code
summary(robot_late_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: amplitude ~ face * valence + (face * valence | num_id) + (1 |
num_id:electrode)
Data: late_data
REML criterion at convergence: 1708.1
Scaled residuals:
Min 1Q Median 3Q Max
-3.3074 -0.5147 -0.0266 0.4850 3.4208
Random effects:
Groups Name Variance Std.Dev. Corr
num_id:electrode (Intercept) 0.1418 0.3765
num_id (Intercept) 2.2325 1.4942
facehuman 2.8872 1.6992 -0.44
valencehappy 1.3154 1.1469 -0.45 0.50
facehuman:valencehappy 4.2956 2.0726 0.29 -0.77 -0.70
Residual 0.5347 0.7312
Number of obs: 552, groups: num_id:electrode, 138; num_id, 46
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.90908 0.23116 44.98732 3.933 0.000287 ***
facehuman 1.32028 0.26555 44.99196 4.972 1.01e-05 ***
valencehappy -0.22630 0.19064 45.00333 -1.187 0.241435
facehuman:valencehappy -0.04292 0.32997 45.00151 -0.130 0.897092
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) fachmn vlnchp
facehuman -0.461
valencehppy -0.465 0.495
fchmn:vlnch 0.309 -0.761 -0.700Code
ranova(robot_late_lmer, reduce.terms = FALSE)ANOVA-like table for random-effects: Single term deletions
Model:
amplitude ~ face + valence + (face * valence | num_id) + (1 | num_id:electrode) + face:valence
npar logLik AIC LRT Df Pr(>Chisq)
<none> 16 -854.07 1740.1
(face * valence | num_id) 6 -1009.17 2030.3 310.199 10 < 2.2e-16 ***
(1 | num_id:electrode) 15 -864.08 1758.2 20.028 1 7.631e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
anova(robot_late_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
face 29.0597 29.0597 1 44.992 54.3507 2.87e-09 ***
valence 1.6806 1.6806 1 44.989 3.1433 0.08301 .
face:valence 0.0090 0.0090 1 45.002 0.0169 0.89709
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
interpret(omega_squared(robot_late_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
---------------------------------------------------------------
face | 0.53 | [0.33, 0.67] | large
valence | 0.04 | [0.00, 0.21] | small
face:valence | 0.00 | [0.00, 0.00] | very small
- Interpretation rule: field2013Code
a_posteriori_lmer(robot_late_lmer)____________________________________________________________NOTE: Results may be misleading due to involvement in interactions$emmeans
face emmean SE df lower.CL upper.CL
robot 0.796 0.205 45 0.383 1.21
human 2.095 0.208 45 1.677 2.51
Results are averaged over the levels of: valence
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
robot - human -1.3 0.176 45 -7.372 <.0001
Results are averaged over the levels of: valence
Degrees-of-freedom method: kenward-roger Late Assumptions:
Code
check_model(robot_late_lmer)
Factorial Mass Univariate analysis (FMUT, cluster-mass)
Groppe, Urbach, and Kutas (2011), 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
Interaction
Face
robot_2023: 4 significant Stimulus cluster(s) out of 28 cluster 1 F-masss: 2836 cluster 1 p-value: 0.0036 cluster 2 F-masss: 1305 cluster 2 p-value: 0.0207 cluster 3 F-masss: 2185 cluster 3 p-value: 0.0065 cluster 4 F-masss: 14412 cluster 4 p-value: 0.0001
Valence
robot_2023: 1 significant Stimulus cluster(s) out of 40 cluster 1 F-masss: 2384 cluster 1 p-value: 0.0011
Multivariate pattern analysis (decoding)
Done in the ADAM 1.14 toolbox (Fahrenfort 2020) with FieldTrip 20211209 (Oostenveld et al. 2011).
Exemplar ERPs on O2
Decoding over time
Individual results
Activation pattern
Temporal generalization plot
Temporal generalization over time
Coda: N170 by Face & Valence, average 150 to 200 ms
Electrodes P7, P8, PO9, PO10
Code
summary(n170_data) worklat start end component
[150.0 200.0]:736 Min. :150 Min. :200 N170:736
1st Qu.:150 1st Qu.:200
Median :150 Median :200
Mean :150 Mean :200
3rd Qu.:150 3rd Qu.:200
Max. :150 Max. :200
amplitude chindex electrode bini stimulus
Min. :-14.282380 14:184 P7 :184 1:184 happy_human :184
1st Qu.: -2.164432 15:184 P8 :184 2:184 happy_robot :184
Median : -0.003614 20:184 PO9 :184 3:184 neutral_human:184
Mean : -0.224474 21:184 PO10:184 4:184 neutral_robot:184
3rd Qu.: 1.938663
Max. : 10.430920
valence face ERPset num_id
neutral:368 robot:368 MRS_04_ct_FT9: 16 4 : 16
happy :368 human:368 MRS_07_ct_FT9: 16 7 : 16
MRS_08_ct_FT9: 16 8 : 16
MRS_09_ct_FT9: 16 9 : 16
MRS_10_ct_FT9: 16 10 : 16
MRS_11_ct_FT9: 16 11 : 16
(Other) :640 (Other):640
antero_posterior medio_lateral
parietal :368 left :368
parieto_occipital:368 right:368
Code
ggplot(
n170_data, aes(x = amplitude, fill = face, color = face)) +
geom_histogram(alpha = .4) +
facet_wrap(~face)
Face & Valence ERPs:
N170 Scalp voltage distribution:
N170 ANOVA:
Code
robot_n170_lmer <- lmer(amplitude ~ face * valence + (face * valence | num_id ) + (1 | num_id : electrode), n170_data)
afex_plot(robot_n170_lmer,
x = 'valence',
trace = 'face',
id = 'num_id',
error_arg = list(width = .15, lwd = .75),
dodge = .5,
data_arg = list(
position =
position_jitterdodge(
jitter.width = .2,
dodge.width = .5 ## needs to be same as dodge
)),
mapping = c('color'),, data_alpha = .3,
point_arg = list(size = 3)
)
Code
summary(robot_n170_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: amplitude ~ face * valence + (face * valence | num_id) + (1 |
num_id:electrode)
Data: n170_data
REML criterion at convergence: 3157.2
Scaled residuals:
Min 1Q Median 3Q Max
-4.9813 -0.4425 0.0024 0.4804 5.0418
Random effects:
Groups Name Variance Std.Dev. Corr
num_id:electrode (Intercept) 2.249 1.500
num_id (Intercept) 4.957 2.226
facehuman 1.175 1.084 -0.31
valencehappy 1.815 1.347 0.04 0.57
facehuman:valencehappy 3.155 1.776 0.15 -0.62 -0.85
Residual 2.134 1.461
Number of obs: 736, groups: num_id:electrode, 184; num_id, 46
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -0.04447 0.36274 45.00052 -0.123 0.903
facehuman -0.21485 0.22074 45.00283 -0.973 0.336
valencehappy -0.12470 0.25028 45.00152 -0.498 0.621
facehuman:valencehappy -0.04091 0.33908 45.00310 -0.121 0.905
Correlation of Fixed Effects:
(Intr) fachmn vlnchp
facehuman -0.347
valencehppy -0.102 0.536
fchmn:vlnch 0.202 -0.659 -0.793Code
ranova(robot_n170_lmer, reduce.terms = FALSE)boundary (singular) fit: see help('isSingular')ANOVA-like table for random-effects: Single term deletions
Model:
amplitude ~ face + valence + (face * valence | num_id) + (1 | num_id:electrode) + face:valence
npar logLik AIC LRT Df Pr(>Chisq)
<none> 16 -1578.6 3189.2
(face * valence | num_id) 6 -1650.3 3312.7 143.46 10 < 2.2e-16 ***
(1 | num_id:electrode) 15 -1663.9 3357.7 170.50 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
anova(robot_n170_lmer)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
face 4.1942 4.1942 1 44.998 1.9657 0.1678
valence 1.8681 1.8681 1 44.998 0.8755 0.3544
face:valence 0.0311 0.0311 1 45.003 0.0146 0.9045Code
interpret(omega_squared(robot_n170_lmer, alternative = 'two.sided'), rules = 'field2013')# Effect Size for ANOVA (Type III)
Parameter | Omega2 (partial) | 95% CI | Interpretation
---------------------------------------------------------------
face | 0.02 | [0.00, 0.16] | small
valence | 0.00 | [0.00, 0.00] | very small
face:valence | 0.00 | [0.00, 0.00] | very small
- Interpretation rule: field2013Code
a_posteriori_lmer(robot_n170_lmer)N170 Assumptions:
Code
check_model(robot_n170_lmer)
Encore: Occipital ↔︎ Infero-Temporal beta (13-30 Hz) Functional Connectivity (ROIconnect)
Multivariate Interaction Measure, MIM (Pellegrini et al. 2023).
1 sec post-stimulus.
Code
summary(connectivity_data) ERPset face region_a region_b
MRS_04_ct_FT9: 8 robot:184 left_occipital :184 inf_temp_ipsi :184
MRS_07_ct_FT9: 8 human:184 right_occipital:184 inf_temp_contra:184
MRS_08_ct_FT9: 8
MRS_09_ct_FT9: 8
MRS_10_ct_FT9: 8
MRS_11_ct_FT9: 8
(Other) :320
median_connectivity log_median_connectivity num_id
Min. :0.05162 Min. :-2.964 Min. : 4.00
1st Qu.:0.07457 1st Qu.:-2.596 1st Qu.:17.00
Median :0.08328 Median :-2.485 Median :30.50
Mean :0.08568 Mean :-2.475 Mean :29.65
3rd Qu.:0.09481 3rd Qu.:-2.356 3rd Qu.:42.00
Max. :0.15935 Max. :-1.837 Max. :53.00
Code
connectivity_data_wide_occipital <- connectivity_data[c('ERPset', 'region_a', 'region_b', 'face', 'log_median_connectivity')] |>
pivot_wider(names_from = region_a, values_from = log_median_connectivity)
connectivity_data_wide_inf_temp <- connectivity_data[c('ERPset', 'region_a', 'region_b', 'face', 'log_median_connectivity')] |>
pivot_wider(names_from = region_b, values_from = log_median_connectivity)
connectivity_data_wide_face <- connectivity_data[c('ERPset', 'region_a', 'region_b', 'face', 'log_median_connectivity')] |>
pivot_wider(names_from = face, values_from = log_median_connectivity)
face <- c('robot' , 'human')
parietal <- c('inf_temp_ipsi' , 'inf_temp_contra')
occipital <- c('left_occipital', 'right_occipital')
occipital_pairs <- ggpairs(connectivity_data_wide_occipital,
aes(colour = face, alpha = .5, linewidth = NA),
columns = occipital,
lower = list(continuous = wrap('points', alpha = .4)),
progress = FALSE
)
inf_temp_pairs <- ggpairs(connectivity_data_wide_inf_temp,
aes(colour = face, alpha = .5, linewidth = NA),
columns = parietal,
lower = list(continuous = wrap('points', alpha = .4)),
progress = FALSE
)
face_pairs <- ggpairs(connectivity_data_wide_face ,
aes(colour = region_b, alpha = .5, linewidth = NA),
columns = face,
lower = list(continuous = wrap('points', alpha = .4)),
progress = FALSE
)
suppressWarnings(print(occipital_pairs))
suppressWarnings(print(inf_temp_pairs))
suppressWarnings(print(face_pairs))
Connectivity ANOVA (glmer):
Code
connectivity_glmer <- glmer(median_connectivity ~ face*region_a*region_b + (face|num_id), family = Gamma(log), connectivity_data)
afex_plot(connectivity_glmer,
x = 'region_b',
trace = 'face',
panel = 'region_a',
error = 'within',
error_arg = list(width = .3, lwd = .75),
dodge = .5,
data_arg = list(
position =
position_jitterdodge(
jitter.width = .2,
dodge.width = .5 ## needs to be same as dodge
)),
mapping = c('color'),, data_alpha = .3,
point_arg = list(size = 3)
)
Code
options(width = 110)
summary(connectivity_glmer)Warning in vcov(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RXWarning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RXGeneralized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Gamma ( log )
Formula: median_connectivity ~ face * region_a * region_b + (face | num_id)
Data: connectivity_data
AIC BIC logLik deviance df.resid
-2292.4 -2245.5 1158.2 -2316.4 356
Scaled residuals:
Min 1Q Median 3Q Max
-2.13733 -0.53481 -0.05061 0.50604 3.13368
Random effects:
Groups Name Variance Std.Dev. Corr
num_id (Intercept) 0.01262 0.1123
facehuman 0.01518 0.1232 -0.59
Residual 0.01503 0.1226
Number of obs: 368, groups: num_id, 46
Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept) -2.510976 0.024515 -102.425 <2e-16 ***
facehuman 0.053852 0.031358 1.717 0.0859 .
region_aright_occipital -0.016505 0.025562 -0.646 0.5185
region_binf_temp_contra 0.034745 0.025562 1.359 0.1741
facehuman:region_aright_occipital -0.016083 0.036150 -0.445 0.6564
facehuman:region_binf_temp_contra 0.003375 0.036150 0.093 0.9256
region_aright_occipital:region_binf_temp_contra 0.023414 0.036150 0.648 0.5172
facehuman:region_aright_occipital:region_binf_temp_contra -0.014902 0.051123 -0.291 0.7707
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) fachmn rgn_r_ rgn___ fch:__ fc:___ r__:__
facehuman -0.655
rgn_rght_cc -0.521 0.408
rgn_bnf_tm_ -0.521 0.408 0.500
fchmn:rgn__ 0.369 -0.576 -0.707 -0.354
fchmn:rg___ 0.369 -0.576 -0.354 -0.707 0.500
rgn_rg_:___ 0.369 -0.288 -0.707 -0.707 0.500 0.500
fchm:__:___ -0.261 0.408 0.500 0.500 -0.707 -0.707 -0.707
optimizer (Nelder_Mead) convergence code: 0 (OK)
unable to evaluate scaled gradient
Model failed to converge: degenerate Hessian with 1 negative eigenvaluesCode
anova(connectivity_glmer)Analysis of Variance Table
npar Sum Sq Mean Sq F value
face 1 0.058375 0.058375 3.8844
region_a 1 0.025244 0.025244 1.6798
region_b 1 0.181478 0.181478 12.0759
face:region_a 1 0.012738 0.012738 0.8476
face:region_b 1 0.000382 0.000382 0.0254
region_a:region_b 1 0.005861 0.005861 0.3900
face:region_a:region_b 1 0.001277 0.001277 0.0850Code
a_posteriori_glmer(connectivity_glmer)Warning in vcov.merMod(mod, complete = FALSE): variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RXAnalysis of Deviance Table (Type II Wald chisquare tests)
Response: median_connectivity
Chisq Df Pr(>Chisq)
face 3.8844 1 0.0487365 *
region_a 1.6798 1 0.1949545
region_b 12.0759 1 0.0005108 ***
face:region_a 0.8476 1 0.3572308
face:region_b 0.0254 1 0.8733144
region_a:region_b 0.3900 1 0.5322995
face:region_a:region_b 0.0850 1 0.7706804
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
____________________________________________________________Warning in vcov.merMod(object, complete = FALSE, ...): variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RXNOTE: Results may be misleading due to involvement in interactions$emmeans
face emmean SE df asymp.LCL asymp.UCL
robot -2.50 0.0189 Inf -2.53 -2.46
human -2.45 0.0182 Inf -2.49 -2.42
Results are averaged over the levels of: region_a, region_b
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
robot - human -0.0438 0.0222 Inf -1.971 0.0487
Results are averaged over the levels of: region_a, region_b
Results are given on the log (not the response) scale.
____________________________________________________________Warning in vcov.merMod(object, complete = FALSE, ...): variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RXNOTE: Results may be misleading due to involvement in interactions$emmeans
region_b emmean SE df asymp.LCL asymp.UCL
inf_temp_ipsi -2.50 0.0162 Inf -2.53 -2.46
inf_temp_contra -2.45 0.0162 Inf -2.48 -2.42
Results are averaged over the levels of: face, region_a
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
inf_temp_ipsi - inf_temp_contra -0.0444 0.0128 Inf -3.475 0.0005
Results are averaged over the levels of: face, region_a
Results are given on the log (not the response) scale. Connectivity Assumptions:
Code
check_model(connectivity_glmer)
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.
Pellegrini, Franziska, Arnaud Delorme, Vadim Nikulin, and Stefan Haufe. 2023. “Identifying Good Practices for Detecting Inter-Regional Linear Functional Connectivity from EEG.” NeuroImage 277 (August): 120218. https://doi.org/10.1016/J.NEUROIMAGE.2023.120218.
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.