── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ purrr 1.0.4
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.2 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tidyr)library(dplyr)library(lme4)
Loading required package: Matrix
Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
library(emmeans)
Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'
library(rstanarm)
Loading required package: Rcpp
This is rstanarm version 2.32.1
- See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
- Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
- For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores())
library(ggplot2)library(effects)
Loading required package: carData
lattice theme set by effectsTheme()
See ?effectsTheme for details.
library(car)
Attaching package: 'car'
The following object is masked from 'package:rstanarm':
logit
The following object is masked from 'package:dplyr':
recode
The following object is masked from 'package:purrr':
some
library(DHARMa)
This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(robustlmm)packageVersion("ggplot2")
[1] '3.5.2'
library(lmerTest)
Attaching package: 'lmerTest'
The following object is masked from 'package:lme4':
lmer
The following object is masked from 'package:stats':
step
dir_path2 <-"/Users/philippschwarzmann/Desktop/Master HFE/Thesis /Data Analysis /Python_Results/Motion"# Read and combine all CSV files into one data frameMotor_Execution_EEG <-list.files(path = dir_path2, pattern ="\\.csv$", full.names =TRUE) %>%map_dfr(read_csv)
Rows: 9318 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9387 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 1008 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9360 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9360 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9360 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9360 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9360 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9333 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9360 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9360 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9459 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9360 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9339 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 9927 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): segment, channel
dbl (7): participant, block, trial, theta_power, beta_power, theta_rel, beta...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Data prep ERDS
#Clean up and format Motor_Execution_EEG <-Motor_Execution_EEG %>%mutate(participant =as.factor(participant), block =as.factor(block), trial =as.factor(trial), segment =as.factor(segment), channel =as.factor(channel))#Clean Extra blocksMotor_Execution_EEG <- Motor_Execution_EEG %>%# 0. Make sure participant really is numericmutate(participant =as.numeric(as.character(participant)),block =as.numeric(as.character(block))) %>%# 1. Drop the unwanted blocks:# - p 3, 7, 11: drop block 1# - p 12: drop block 2filter(!(participant %in%c(3,7,11) & block ==1),!(participant ==12& block ==2) ) %>%# 2. Shift down the remaining block numbers:mutate(block =case_when( participant %in%c(3,7,11) & block >=2~ block -1, participant ==12& block >=3~ block -1,TRUE~ block )) %>%# 3. Back to a factor with levels 1–5mutate(block =factor(block, levels =1:5), participant =as.factor(participant))Motor_ex_CZ <- Motor_Execution_EEG %>%filter(channel=="Cz") library(dplyr) Motor_ex_CZ <- Motor_ex_CZ %>%mutate(segment =factor(segment, levels =c("baseline", paste0("step", 1:18))))df_erds <- Motor_ex_CZ %>%# 1. extract one baseline-per-trial (mean if multiple)filter(segment =="baseline") %>%group_by(participant, block, trial) %>%summarise(baseline_theta =mean(theta_power, na.rm =TRUE),.groups ="drop") %>%# 2. join back & compute ERDS for every rowright_join(Motor_ex_CZ, by =c("participant", "block", "trial")) %>%mutate(theta_erds =case_when( segment =="baseline"~0, # or NA, if you preferTRUE~100* (theta_power - baseline_theta) / baseline_theta ) ) %>%select(-baseline_theta)q_high <-quantile(df_erds$theta_erds, 0.95, na.rm =TRUE)q_low <-quantile(df_erds$theta_erds, 0.05, na.rm =TRUE)df_erds$theta_erds_wins <-pmax(pmin(df_erds$theta_erds, q_high), q_low)df_erds <- df_erds %>%mutate(step= segment ) %>%filter (segment !="baseline")#transform eeg data for mergedf_erds <- df_erds %>%# (1) Make sure block & trial are numericmutate(block =as.integer(block),trial =as.integer(trial) ) %>%# (2) Within each block, re-rank the trial values to 1–48group_by(block,participant) %>%mutate(trial =dense_rank(trial)) %>%ungroup()df_erds<- df_erds %>%mutate(block=as.factor(block),trial=as.factor(trial))CZ_B1 <- df_erds %>%filter(block==1, step !="baseline", step !="step12",step !="step11",step !="step10",step !="step9",step !="step8",step !="step7") CZ_B2 <- df_erds %>%filter(block==2,step !="baseline") CZ_B3 <- df_erds %>%filter(block==3, step !="baseline")#df_erds <- df_erds %>% rename(step=segment)CZ_B4 <- df_erds %>%filter(block==4, step !="baseline")CZ_B5 <- df_erds %>%filter(block==5, step !="baseline")df_erds_learning_all <-bind_rows(CZ_B1,CZ_B2,CZ_B3)CZ_B3 %>%distinct(participant) %>% nrow
df_learning_acc <- df_learning %>%filter(trial.acc >=0.8)df_test_acc <- df_test %>%filter(trial.acc >=0.8)#Step 1: Exclude universally too-fast or too-slow RTsdf_learning_acc <- df_learning_acc%>%filter(RT >=150, RT <=3000)# Step 2: Exclude outliers based on ±3 SD within each participantdf_learning_acc <- df_learning_acc %>%group_by(id) %>%mutate(meanRT =mean(RT, na.rm =TRUE),sdRT =sd(RT, na.rm =TRUE)) %>%filter(RT >= meanRT -3* sdRT, RT <= meanRT +3* sdRT) %>%ungroup() %>%select(-meanRT, -sdRT) # Clean up if not neededdf_test_acc <- df_test_acc%>%filter(RT >=150, RT <=3000)df_test_acc <- df_test_acc %>%group_by(id) %>%mutate(meanRT =mean(RT, na.rm =TRUE),sdRT =sd(RT, na.rm =TRUE)) %>%filter(RT >= meanRT -3* sdRT, RT <= meanRT +3* sdRT) %>%ungroup() %>%select(-meanRT, -sdRT) # Clean up if not neededdf_learning_acc %>%nrow()
[1] 15773
df_test_acc %>%nrow()
[1] 9556
Models theta learning
# Fit with ML for valid likelihood-ratio test M1_theta_learning <-lmer(theta_erds_wins ~ step * block + (1| id) + (1| trial) + (1| trial)+ (1| block),data = df_learning_acc, REML =TRUE)
fixed-effect model matrix is rank deficient so dropping 18 columns / coefficients
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge: degenerate Hessian with 1 negative eigenvalues
Warning: Model failed to converge with 1 negative eigenvalue: -1.4e+00
# Be careful with nesting; see notes belowM2_theta_learning <-lmer(theta_erds_wins ~ step * block + (1| id) + (1| id:block) , # trials within subjectsdata = df_learning_acc, REML =TRUE)
fixed-effect model matrix is rank deficient so dropping 18 columns / coefficients
# Compare models (likelihood-ratio test + AIC/BIC)anova(M1_theta_learning, M2_theta_learning)
fixed-effect model matrix is rank deficient so dropping 36 columns / coefficients
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge: degenerate Hessian with 1 negative eigenvalues
Warning: Model failed to converge with 1 negative eigenvalue: -9.0e-05
# Be careful with nesting; see notes belowM2_RT_test <-lmer(RT ~ cond*step + step*difficulty + cond*difficulty + cond*difficulty*step + (1| id) + (1| id:block) +# sessions within subjects (1| id:trial),data = df_test_acc, REML =TRUE)
fixed-effect model matrix is rank deficient so dropping 36 columns / coefficients
# Compare models (likelihood-ratio test + AIC/BIC)anova(M1_RT_test, M2_RT_test)