OSCR

Locus coeruleus modulation of neurophysiological sensory selectivity differs in autism and other mental health conditions.

Code ↔ Paper

11 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 11 matches
  1. [1] § Methods › Pupillometry ↔ data_preprocessing_visual_oddball.R, lines 84 to 143 · score 0.88 · dilation speed outliers, median absolute deviation, corrected pupil, gaps, imputed, blink · Report a mismatch
  2. [2] § Methods › Pupillometry ↔ data_preprocessing_auditory_oddball.R, lines 330 to 443 · score 0.87 · dilation speed outliers, median absolute deviation, corrected pupil diameter, blink, preprocessed, eye · Report a mismatch
  3. [3] § Results › Stimulus effects ↔ data_analysis_auditory_oddball.R, lines 60 to 122 · score 0.83 · 100–150 ms, 150–250 ms, 500–1500 ms, Baseline, pupil, P3a · Report a mismatch
  4. [4] § Methods › Electroencephalography ↔ data_analysis_auditory_oddball.R, lines 1150 to 1206 · score 0.71 · 100–150 ms, 150–250 ms, cut, ERPs, P3a, 100 ms · Report a mismatch
  5. [5] § Methods › Measures ↔ SEGA_Preprocessing.Rmd, lines 1783 to 1880 · score 0.66 · matrix reasoning, block design, subtests, vocabulary, WAIS, CBCL · Report a mismatch
  6. [6] § Methods › Measures ↔ data_preprocessing_auditory_oddball.R, lines 1749 to 1793 · score 0.66 · matrix reasoning, block design, subtests, vocabulary, WAIS, CBCL · Report a mismatch
  7. [7] § Methods › Auditory oddball task ↔ data_preprocessing_auditory_oddball.R, lines 614 to 660 · score 0.56 · hand dynamometer, Grip strength, counterbalanced, oddballs, Auditory, blocks · Report a mismatch
  8. [8] § Results › Associations of the pupillometric measures ↔ SEGA_supplements.Rmd, lines 796 to 849 · score 0.56 · exploratory model fit, BIC, linear, predictor, SEPR, ERPs · Report a mismatch
  9. [9] § Results › Covariate effects ↔ SEGA_Preprocessing.Rmd, lines 1783 to 1880 · score 0.55 · verbal IQ, Grip strength, dependent variables, MMN latency, gender, age · Report a mismatch
  10. [10] § Results › Covariate effects ↔ data_analysis_auditory_oddball.R, lines 1670 to 1717 · score 0.51 · verbal IQ, Grip strength, MMN latency, covariates, P3a amplitude, models · Report a mismatch
  11. [11] § Results › Associations of the pupillometric measures ↔ data_analysis_auditory_oddball.R, lines 1446 to 1496 · score 0.51 · exploratory model fit, predictor, pupillometric, SEPR, ERPs, BPS · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 2,188 lines · 81 KB · no license · 4 matches

  1. ### --- LOADS preprocessed auditory oddball data and analysis ERPs and pupillometry variables
  3. # SETUP ####
  4. # REQUIRED PACKAGES
  5. require(performance) # marginal + conditional R2
  6. require(lme4) # linear-mixed-effects models
  7. require(lmerTest, warn.conflicts = FALSE) # linear-mixed-effects models
  8. require(emmeans, warn.conflicts = FALSE) # estimated marginal means (EMMs)
  9. library(ggpubr, warn.conflicts = FALSE) # ggscatter()-function
  10. library(knitr) # dynamic report generation
  11. library(kableExtra) # table formatting
  12. library(tidyverse)
  13. library(cowplot, warn.conflicts = FALSE) # get only legend from plot
  14. library(simr) # power analysis
  15. library(gridExtra) # arrange ggplots
  17. # PATHS
  18. if (Sys.info()["sysname"] == "Linux") {
  19. home_path <- "~"
  20. project_path <- "/PowerFolders/project_sega"
  21. data_path <- "/PowerFolders/project_sega/data/AuditoryOddball"
  22. data_path_eeg <- "/PowerFolders/project_sega/data/AuditoryOddball_EEG"
  23. datapath <- paste0(home_path, data_path) # .csv + .hdf5 input files
  24. datapath_eeg <- paste0(home_path, data_path_eeg) # .txt input files (eeg)
  25. # List all .hdf and .csv files
  26. data_files <- list.files(path = datapath, full.names = TRUE)
  27. }
  29. if (Sys.info()["sysname"] == "Windows") {
  30. home_path <- "C:/Users/Nico"
  31. project_path <- "/PowerFolders/project_sega"
  32. data_path <- "/PowerFolders/project_sega/data/et auditory oddball"
  33. data_path_task<-"/PowerFolders/project_sega/data/task auditory oddball"
  34. data_path_eeg <- "/PowerFolders/project_sega/data/tests/eeg preprocessed auditory oddball"
  35. datapath <- paste0(home_path, data_path) # hdf5 input files
  36. datapath_task <- paste0(home_path, data_path_task) # hdf5 input files
  37. datapath_eeg <- paste0(home_path, data_path_eeg) # .txt input files (eeg)
  38. # List all .hdf and .csv files
  39. data_files <- c(list.files(path = datapath, full.names = TRUE),
  40. list.files(path = datapath_task, full.names = TRUE))
  41. }
  43. if (Sys.info()["sysname"] == "Darwin") {
  44. home_path <- "/Volumes/common/"
  45. project_path <- "/code"
  46. data_path <- "/code/input/AuditoryOddball/eyetracking"
  47. data_path_task <- "/code/input/AuditoryOddball/taskdata"
  48. data_path_eeg <- "/code/input/AuditoryOddball/eeg"
  49. datapath <- paste0(home_path, data_path) # .hdf5 input files
  50. datapath_task <- paste0(home_path, data_path_task) # .csv input files
  51. datapath_eeg <- paste0(home_path, data_path_eeg) # .txt input files (eeg)
  52. data_path_single_trial_eeg <- "/code/input/AuditoryOddball/eeg_single_trial"
  53. datapath_single_eeg <- paste0(home_path, data_path_single_trial_eeg)
  54. # List all .hdf and .csv files
  55. data_files <- c(
  56. list.files(path = datapath, full.names = TRUE),
  57. list.files(path = datapath_task, full.names = TRUE))
  58. }
  60. # Can be used to skip preprocessing and directly read proprocessed data from .rds file:
  61. df_trial <- readRDS(paste0(home_path,project_path,'/data/preprocessed_auditory_ETdata.rds')) # trial pupil data including ALL trials
  62. df <- readRDS(paste0(home_path, project_path,'/data/preprocessed_auditory_df.rds')) # pupil data (all data points for visualizing)
  63. et_erp_subject <- readRDS(paste0(home_path,project_path,'/data/preprocessed_auditory_ET_ERP_subject.rds')) # eeg + pupil (subject level, only oddball(rev) blocks)
  64. et_erp_trial <- readRDS(paste0(home_path,project_path,'/data/preprocessed_auditory_ET_ERP_trial.rds')) # eeg + pupil(single-trial, only oddball(rev) blocks)
  65. MMN_diff <- readRDS(paste0(home_path,project_path,'/data/preprocessed_auditory_MMN_diff.rds'))
  67. # Renaming the variable "trial" into "stimulus" to match terminology in the paper
  68. names(et_erp_subject)[names(et_erp_subject) == "trial"] <- "stimulus"
  69. names(et_erp_trial)[names(et_erp_trial) == "trial"] <- "stimulus"
  70. names(df_trial)[names(df_trial) == "trial"] <- "stimulus"
  71. names(df)[names(df) == "trial"] <- "stimulus"
  73. # DATA ANALYSIS ON SUBJECT LEVEL ####
  74. # Distributions of dependent variables
  75. par(mfrow = c(3,2), mar = c(4, 4, 2, 1))
  76. hist(et_erp_subject$z_rpd,
  77. main = "Distribution of rpd (500-1500 ms)",
  78. xlab = "rpd",
  79. xlim = c(-4, 6),
  80. breaks = 200)
  81. hist(et_erp_subject$z_rpd_low,
  82. main = "Distribution of rpd_low (0-250 ms)",
  83. xlab = "rpd_low",
  84. xlim = c(-4, 6),
  85. breaks = 200)
  86. hist(et_erp_subject$z_rpd_block,
  87. main = "Distribution of rpd_block",
  88. xlab = "rpd_block",
  89. xlim = c(-4, 4),
  90. breaks = 200)
  91. hist(et_erp_subject$z_MMN_amplitude,
  92. main = "Distribution of MMN amplitude (100-150 ms)",
  93. xlab = "MMN amplitude (FC1, FC2, FCz, Fz)",
  94. breaks = 200)
  95. hist(et_erp_subject$z_P3a_amplitude,
  96. main = "Distribution of P3a amplitude (150-250 ms)",
  97. xlab = "P3a amplitude (Cz, FCz)",
  98. ylim = c(0, 30),
  99. xlim = c(-4, 6),
  100. breaks = 200)
  101. hist(et_erp_subject$z_P3b_amplitude,
  102. main = "Distribution of P3b amplitude (250-500 ms)",
  103. xlab = "P3b amplitude (Pz)",
  104. xlim = c(-6, 8),
  105. breaks = 200)
  106. par(mfrow = c(1, 1))
  108. # New data frame: 1 row per subject containing sample description data (consistent across conditions)
  109. sample_description_df <- et_erp_subject %>%
  110. group_by(SEGA_ID) %>%
  111. summarise(across(everything(), ~ {
  112. vals <- na.omit(.x) # removes NAs so we only check actual values
  113. # keeps the column if all non-NA values are identical + returns that consistent value (or NA if all were NA)
  114. if (length(unique(vals)) <= 1) unique(vals)[1] else NA
  115. }), .groups = "drop") %>%
  116. # Keep columns that have at least some non-NA values + removes columns that are completely NA across all subjects
  117. select(where(~ any(!is.na(.))))
  118. sample_description_df$block_baseline_mean <- NULL
  120. # Sample size
  121. nrow(sample_description_df)
  122. table(sample_description_df$group)
  124. # Function returning group mean + sd
  125. fun_return_descriptives <- function(group){
  126. group_df <- sample_description_df[
  127. sample_description_df$group == group, ]
  128. # gender
  129. n <- length(unique(group_df$SEGA_ID))
  130. male <- (length(which(group_df$gender == "männlich")))
  131. female <- (length(which(group_df$gender == "weiblich")))
  132. gender_f_m <- paste(female,"/",male)
  133. # age
  134. age_mean <- round(mean(group_df$age, na.rm = TRUE), digits = 1)
  135. age_sd <- round(sd(group_df$age, na.rm = TRUE), digits = 1)
  136. age <- paste(age_mean, "(",age_sd,")" )
  137. # SRS
  138. SRS_mean <- round(mean(group_df$SRS, na.rm = TRUE), digits = 1)
  139. SRS_sd <- round(sd(group_df$SRS, na.rm = TRUE), digits = 1)
  140. SRS <- paste(SRS_mean, "(",SRS_sd,")")
  141. # CBCL
  142. CBCL_mean <- round(mean(group_df$CBCL, na.rm = TRUE), digits = 1)
  143. CBCL_sd <- round(sd(group_df$CBCL, na.rm = TRUE), digits = 1)
  144. CBCL <- paste(CBCL_mean, "(",CBCL_sd,")")
  145. # SCQ
  146. SCQ_mean <- round(mean(group_df$SCQ, na.rm = TRUE), digits = 1)
  147. SCQ_sd <- round(sd(group_df$SCQ, na.rm = TRUE), digits = 1)
  148. SCQ <- paste(SCQ_mean, "(", SCQ_sd, ")")
  149. # YSR
  150. YSR_mean <- round(mean(group_df$YSR, na.rm = TRUE), digits = 1)
  151. YSR_sd <- round(sd(group_df$YSR, na.rm = TRUE), digits = 1)
  152. YSR <- paste(YSR_mean, "(", YSR_sd, ")")
  153. # SP2
  154. SP2_mean <- round(mean(group_df$SP2, na.rm = TRUE), digits = 1)
  155. SP2_sd <- round(sd(group_df$SP2, na.rm = TRUE), digits = 1)
  156. SP2 <- paste(SP2_mean, "(", SP2_sd, ")")
  157. # grip strength
  158. grip_strength_mean <- round(mean(group_df$z_grip_strength, na.rm = TRUE), digits = 1)
  159. grip_strength_sd <- round(sd(group_df$z_grip_strength, na.rm = TRUE), digits = 1)
  160. grip_strength <- paste(grip_strength_mean, "(", grip_strength_sd, ")")
  161. # verbal IQ
  162. verbal_IQ_mean <- round(mean(group_df$verbal_IQ, na.rm = TRUE), digits = 1)
  163. verbal_IQ_sd <- round(sd(group_df$verbal_IQ, na.rm = TRUE), digits = 1)
  164. verbal_IQ <- paste(verbal_IQ_mean, "(", verbal_IQ_sd, ")")
  165. # non-verbal IQ
  166. non_verbal_IQ_mean <- round(mean(group_df$non_verbal_IQ, na.rm = TRUE), digits = 1)
  167. non_verbal_IQ_sd <- round(sd(group_df$non_verbal_IQ, na.rm = TRUE), digits = 1)
  168. non_verbal_IQ <- paste(non_verbal_IQ_mean, "(", non_verbal_IQ_sd, ")")
  169. # included trials eeg
  170. included_trials_eeg_mean <- round((mean(group_df$included_trials_eeg))*100, digits = 1)
  171. included_trials_eeg_sd <- round((sd(group_df$included_trials_eeg))*100, digits = 1 )
  172. included_trials_eeg <- paste(included_trials_eeg_mean, "(", included_trials_eeg_sd, ")")
  173. # included trials et
  174. included_trials_et_mean <- round((mean(group_df$included_trials_et))*100, digits = 1)
  175. included_trials_et_sd <- round((sd(group_df$included_trials_et))*100, digits = 1)
  176. included_trials_et <- paste(included_trials_et_mean, "(", included_trials_et_sd, ")")
  178. group_description <- data.frame(
  179. n,
  180. gender_f_m,
  181. age,
  182. verbal_IQ,
  183. non_verbal_IQ,
  184. CBCL,
  185. YSR,
  186. SRS,
  187. SCQ,
  188. SP2,
  189. grip_strength,
  190. included_trials_et,
  191. included_trials_eeg)
  192. t(group_description)
  193. }
  195. asd_description <- fun_return_descriptives(group = "ASD")
  196. colnames(asd_description) <- "ASD"
  197. con_description <- fun_return_descriptives(group = "CON")
  198. colnames(con_description) <- "CON"
  199. mhc_description <- fun_return_descriptives(group = "MHC")
  200. colnames(mhc_description) <- "MHC"
  202. ## p-values for descriptive statistics
  203. SCQ_anova <- aov(SCQ ~ group, data = sample_description_df)
  204. SCQ_anova_p <- summary(SCQ_anova)[[1]][["Pr(>F)"]][[1]]
  205. age_anova <- aov(age ~ group, data = sample_description_df)
  206. age_anova_p <- summary(age_anova)[[1]][["Pr(>F)"]][[1]]
  207. CBCL_anova <- aov(CBCL ~ group, data = sample_description_df)
  208. CBCL_anova_p <- summary(CBCL_anova)[[1]][["Pr(>F)"]][[1]]
  209. SRS_anova <- aov(SRS ~ group, data = sample_description_df)
  210. SRS_anova_p <- summary(SRS_anova)[[1]][["Pr(>F)"]][[1]]
  211. YSR_anova <- aov(YSR ~ group, data = sample_description_df)
  212. YSR_anova_p <- summary(YSR_anova)[[1]][["Pr(>F)"]][[1]]
  213. SP2_anova <- aov(SP2 ~ group, data = sample_description_df)
  214. SP2_anova_p <- summary(SP2_anova)[[1]][["Pr(>F)"]][[1]]
  215. grip_strength_anova <- aov(z_grip_strength ~ group, data = sample_description_df)
  216. grip_strength_anova_p <- summary(grip_strength_anova)[[1]][["Pr(>F)"]][[1]]
  217. verbal_IQ_anova <- aov(verbal_IQ ~ group, data = sample_description_df)
  218. verbal_IQ_anova_p <- summary(verbal_IQ_anova)[[1]][["Pr(>F)"]][[1]]
  219. non_verbal_IQ_anova <- aov(non_verbal_IQ ~ group, data = sample_description_df)
  220. non_verbal_IQ_anova_p <- summary(non_verbal_IQ_anova)[[1]][["Pr(>F)"]][[1]]
  221. gender_chi2 <- chisq.test(sample_description_df$gender, sample_description_df$group)
  222. gender_chi2_p <- gender_chi2$p.value
  223. included_trials_et_anova <- aov(included_trials_et ~ group, data = sample_description_df)
  224. included_trials_et_p <- summary(included_trials_et_anova)[[1]][["Pr(>F)"]][[1]]
  225. included_trials_eeg_anova <- aov(included_trials_eeg ~ group, data = sample_description_df)
  226. included_trials_eeg_p <- summary(included_trials_eeg_anova)[[1]][["Pr(>F)"]][[1]]
  228. p_values <- c(
  229. NA,
  230. gender_chi2_p,
  231. age_anova_p,
  232. SRS_anova_p,
  233. CBCL_anova_p,
  234. SCQ_anova_p,
  235. SP2_anova_p,
  236. YSR_anova_p,
  237. grip_strength_anova_p,
  238. verbal_IQ_anova_p,
  239. non_verbal_IQ_anova_p,
  240. included_trials_et_p,
  241. included_trials_eeg_p)
  243. p_value <- sapply(p_values, function(p) {
  244. if (is.na(p)) {
  245. NA
  246. } else if (p < 0.001) {
  247. "< 0.001"
  248. } else {
  249. format(p, scientific = F)
  250. round(p, 3)
  251. }
  252. })
  254. sample_table <- cbind(asd_description, con_description, mhc_description, p_value)
  255. # Redo table row and column names
  256. rownames(sample_table)[rownames(sample_table) == "grip_strength"] <- "grip strength [z]"
  257. rownames(sample_table)[rownames(sample_table) == "gender_f_m"] <- "gender (f/m)"
  258. rownames(sample_table)[rownames(sample_table) == "verbal_IQ"] <- "verbal IQ"
  259. rownames(sample_table)[rownames(sample_table) == "non_verbal_IQ"] <- "non verbal IQ"
  260. rownames(sample_table)[rownames(sample_table) == "included_trials_eeg"] <- "included eeg trials [%]"
  261. rownames(sample_table)[rownames(sample_table) == "included_trials_et"] <- "included et trials [%]"
  262. colnames(sample_table)[colnames(sample_table) == "p_value"] <- "p value"
  264. sample_description_table <- kable(sample_table, caption = "Sample description", digits = 4) %>%
  265. kable_classic(full_width = F, html_font = "Cambria") %>% kable_styling
  266. sample_description_table
  268. # Supplements: No of trials per condition
  269. fun_trials_per_cond <- function(group, variable){
  270. group_df <- droplevels(et_erp_trial[et_erp_trial$group == group, ])
  271. # subjects starting with 500 Hz oddball
  272. ## block 1
  273. df_before_forward_O500_S750 <- group_df %>% filter(
  274. (manipulation == "before" & phase == "oddball_block" & stimulus == "oddball" & pitch == "500") |
  275. (manipulation == "before" & phase == "oddball_block" & stimulus == "standard" & pitch == "750"))
  276. sample_size <- length(unique(df_before_forward_O500_S750$SEGA_ID))
  277. max_trials <- sample_size *100
  278. before_forward_O500_S750 <- paste(round(((sum(!is.na(df_before_forward_O500_S750[variable]), na.rm = T))/max_trials)*100, 2), "%")
  279. ## block 2
  280. df_before_reverse_O750_S500 <- group_df %>% filter(
  281. (manipulation == "before" & phase == "oddball_block_rev" & stimulus == "oddball" & pitch == "750") |
  282. (manipulation == "before" & phase == "oddball_block_rev" & stimulus == "standard" & pitch == "500"))
  283. before_reverse_O750_S500 <- paste(round(((sum(!is.na(df_before_reverse_O750_S500[variable]), na.rm = T))/max_trials)*100, 2), "%")
  284. ## block 3
  285. df_after_forward_O500S750 <- group_df %>% filter(
  286. (manipulation == "after" & phase == "oddball_block" & stimulus == "oddball" & pitch == "500") |
  287. (manipulation == "after" & phase == "oddball_block" & stimulus == "standard" & pitch == "750"))
  288. after_forward_O500S750 <- paste(round(((sum(!is.na(df_after_forward_O500S750[variable]), na.rm = T))/max_trials)*100, 2), "%")
  289. # block 4
  290. df_after_reverse_O750S500 <- group_df %>% filter(
  291. (manipulation == "after" & phase == "oddball_block_rev" & stimulus == "oddball" & pitch == "750") |
  292. (manipulation == "after" & phase == "oddball_block_rev" & stimulus == "standard" & pitch == "500"))
  293. after_reverse_O750S500 <- paste(round(((sum(!is.na(df_after_reverse_O750S500[variable]), na.rm = T))/max_trials)*100, 2), "%")
  294. # subjects starting with 750 Hz oddball
  295. ## block 1
  296. df_before_forward_O750_S500 <- group_df %>% filter(
  297. (manipulation == "before" & phase == "oddball_block" & stimulus == "oddball" & pitch == "750") |
  298. (manipulation == "before" & phase == "oddball_block" & stimulus == "standard" & pitch == "500"))
  299. sample_size <- length(unique(df_before_forward_O750_S500$SEGA_ID))
  300. max_trials <- sample_size *100
  301. before_forward_O750_S500 <- paste(round(((sum(!is.na(df_before_forward_O750_S500[variable]), na.rm = T))/max_trials)*100, 2), "%")
  302. ## block 2
  303. df_before_reverse_O500S750 <- group_df %>% filter(
  304. (manipulation == "before" & phase == "oddball_block_rev" & stimulus == "oddball" & pitch == "500") |
  305. (manipulation == "before" & phase == "oddball_block_rev" & stimulus == "standard" & pitch == "750"))
  306. before_reverse_O500S750 <- paste(round(((sum(!is.na(df_before_reverse_O500S750[variable]), na.rm = T))/max_trials)*100, 2), "%")
  307. ## block 3
  308. df_after_forward_O750_S500 <- group_df %>% filter(
  309. (manipulation == "after" & phase == "oddball_block" & stimulus == "oddball" & pitch == "750") |
  310. (manipulation == "after" & phase == "oddball_block" & stimulus == "standard" & pitch == "500"))
  311. after_forward_O750_S500 <- paste(round(((sum(!is.na(df_after_forward_O750_S500[variable]), na.rm = T))/max_trials)*100, 2), "%")
  312. ## block 4
  313. df_after_reverse_O500S750 <- group_df %>% filter(
  314. (manipulation == "after" & phase == "oddball_block_rev" & stimulus == "oddball" & pitch == "500") |
  315. (manipulation == "after" & phase == "oddball_block_rev" & stimulus == "standard" & pitch == "750"))
  316. after_reverse_O500S750 <- paste(round(((sum(!is.na(df_after_reverse_O500S750[variable]), na.rm = T))/max_trials)*100, 2), "%")
  318. trial_per_condition <- data.frame(
  319. before_forward_O500_S750,
  320. before_reverse_O750_S500,
  321. after_forward_O500S750,
  322. after_reverse_O750S500,
  323. before_forward_O750_S500,
  324. before_reverse_O500S750,
  325. after_forward_O750_S500,
  326. after_reverse_O500S750)
  327. t(trial_per_condition)
  328. }
  330. # Call function for each group for EEG
  331. asd_trials_per_cond_eeg <- fun_trials_per_cond(group = "ASD", variable = "z_MMN_amplitude")
  332. colnames(asd_trials_per_cond_eeg) <- "ASD"
  333. con_trials_per_cond_eeg <- fun_trials_per_cond(group = "CON", variable = "z_MMN_amplitude")
  334. colnames(con_trials_per_cond_eeg) <- "CON"
  335. mhc_trials_per_cond_eeg <- fun_trials_per_cond(group = "MHC", variable = "z_MMN_amplitude")
  336. colnames(mhc_trials_per_cond_eeg) <- "MHC"
  337. # Call function for each group for ET
  338. asd_trials_per_cond_et <- fun_trials_per_cond(group = "ASD", variable = "z_rpd")
  339. colnames(asd_trials_per_cond_et) <- "ASD"
  340. con_trials_per_cond_et <- fun_trials_per_cond(group = "CON", variable = "z_rpd")
  341. colnames(con_trials_per_cond_et) <- "CON"
  342. mhc_trials_per_cond_et <- fun_trials_per_cond(group = "MHC", variable = "z_rpd")
  343. colnames(mhc_trials_per_cond_et) <- "MHC"
  345. condition_table <- cbind(
  346. asd_trials_per_cond_eeg, con_trials_per_cond_eeg, mhc_trials_per_cond_eeg,
  347. asd_trials_per_cond_et, con_trials_per_cond_et, mhc_trials_per_cond_et)
  349. condition_trials_table <- kable(
  350. condition_table, align = "c", caption = "Number of included trials per condition") %>%
  351. kable_classic(full_width = F, html_font = "Cambria") %>%
  352. kable_styling() %>%
  353. add_header_above(c(" " = 1, "Eye Tracking" = 3, "EEG" = 3))
  354. condition_trials_table
  356. ## Plot SCQ
  357. ggplot(sample_description_df, aes(x = group, y = SCQ), col = "group") +
  358. geom_boxplot(fill = "grey") +
  359. geom_signif(comparisons = list(c("ASD", "CON"), c("ASD", "MHC"), c("MHC", "CON")),
  360. map_signif_level=TRUE,
  361. y_position = c(24, 25.5, 21)) +
  362. theme_bw() +
  363. theme_classic() +
  364. ggtitle("SCQ") +
  365. theme(plot.title = element_text(face = "bold"))
  367. ## Plot SRS
  368. ggplot(sample_description_df, aes(x = group, y = SRS), col = "group") +
  369. geom_boxplot(fill = "grey") +
  370. geom_signif(comparisons = list(c("ASD", "CON"), c("ASD", "MHC"), c("MHC", "CON")),
  371. map_signif_level=TRUE,
  372. y_position = c(42, 45, 40)) +
  373. theme_bw() +
  374. theme_classic() +
  375. ggtitle("SRS") +
  376. theme(plot.title = element_text(face = "bold"))
  378. ## Plot CBCL
  379. ggplot(sample_description_df, aes(x = group, y = CBCL), col = "group") +
  380. geom_boxplot(fill = "grey") +
  381. geom_signif(comparisons = list(c("ASD", "CON"), c("ASD", "MHC"), c("MHC", "CON")),
  382. map_signif_level=TRUE,
  383. y_position = c(87, 90, 84)) +
  384. theme_bw() +
  385. theme_classic() +
  386. ggtitle("CBCL") +
  387. theme(plot.title = element_text(face = "bold"))
  389. ## Plot: YSR
  390. ggplot(sample_description_df, aes(x = group, y = YSR), col = "group") +
  391. geom_boxplot(fill = "grey") +
  392. geom_signif(comparisons = list(c("ASD", "CON"), c("ASD", "MHC"), c("MHC", "CON")),
  393. map_signif_level=TRUE,
  394. y_position = c(87, 90, 84)) +
  395. theme_bw() +
  396. theme_classic() +
  397. ggtitle("YSR") +
  398. theme(plot.title = element_text(face="bold"))
  400. ## Plot: SP2-Auditiv
  401. ggplot(sample_description_df, aes(x = group, y = SP2), col = "group") +
  402. geom_boxplot(fill = "grey") +
  403. geom_signif(comparisons = list(c("ASD", "CON"), c("ASD", "MHC"), c("MHC", "CON")),
  404. map_signif_level=TRUE,
  405. y_position = c(43, 46, 40)) +
  406. theme_bw() +
  407. theme_classic() +
  408. ggtitle("SP2-Auditiv") +
  409. theme(plot.title = element_text(face="bold"))
  411. ## Plot: Age
  412. ggplot(sample_description_df, aes(x = group, y = age), col = group) +
  413. geom_boxplot(fill = "grey") +
  414. geom_jitter() +
  415. geom_signif(comparisons = list(c("ASD", "CON"), c("ASD", "MHC"), c("MHC", "CON")),
  416. map_signif_level=TRUE,
  417. y_position = c(19.5, 20, 19)) +
  418. theme_bw() +
  419. theme_classic() +
  420. ggtitle("Age") +
  421. theme(plot.title = element_text(face="bold"))
  423. ## Plot: IQ-verbal
  424. plot_verbal_iq <- ggplot(sample_description_df, aes(x = group, y = verbal_IQ), col = group) +
  425. geom_boxplot(fill = "grey") +
  426. geom_signif(comparisons = list(c("ASD", "CON"), c("ASD", "MHC"), c("MHC", "CON")),
  427. map_signif_level=TRUE,
  428. y_position = c(135, 140, 130)) +
  429. theme_bw() +
  430. theme_classic() +
  431. ggtitle("IQ-verbal") +
  432. theme(plot.title = element_text(face="bold")) +
  433. theme(plot.margin = unit(c(1, 0.5, 0.5, 0.5), "cm"))
  435. ## Plot: IQ-non-verbal
  436. plot_non_verbal_iq <- ggplot(sample_description_df, aes(x = group, y = non_verbal_IQ), col = "group") +
  437. geom_boxplot(fill = "grey") +
  438. geom_signif(comparisons = list(c("ASD", "CON"), c("ASD", "MHC"), c("MHC", "CON")),
  439. map_signif_level=TRUE,
  440. y_position = c(135, 140, 130)) +
  441. theme_bw() +
  442. theme_classic() +
  443. ggtitle("IQ-non-verbal") +
  444. theme(plot.title = element_text(face="bold")) +
  445. theme(plot.margin = unit(c(1, 0.5, 0.5, 0.5), "cm"))
  447. grid.arrange(plot_verbal_iq, plot_non_verbal_iq, ncol = 2)
  449. # Trial duration by group (REVIEWER 2025-10)
  450. ## on TRIAL level
  451. et_erp_trial$stimulus_isi_duration <- et_erp_trial$stimulus_duration + et_erp_trial$ISI_duration
  452. ggplot(et_erp_trial, aes(x = stimulus_isi_duration)) +
  453. geom_histogram(binwidth = 0.020, fill = "grey", color = "white") +
  454. labs(x = "Duration (s)", y = "Count") +
  455. coord_cartesian(xlim = c(1.85, 3.4)) +
  456. facet_wrap(~ group) + # separate histogram per group
  457. theme_minimal() +
  458. theme(
  459. strip.text = element_text(face = "bold"))
  461. et_erp_trial %>%
  462. group_by(group) %>%
  463. summarise(
  464. mean_stimulus_isi = round(mean(stimulus_isi_duration, na.rm = TRUE), 1),
  465. sd_stimulus_isi = round(sd(stimulus_isi_duration, na.rm = TRUE), 1))
  467. stimulus_isi_duration_anova <- aov(stimulus_isi_duration ~ group, data = et_erp_trial)
  468. stimulus_isi_duration_p <- summary(stimulus_isi_duration_anova)[[1]][["Pr(>F)"]][[1]]
  470. ## on PARTICIPANT level
  471. et_erp_trial <- et_erp_trial %>% # same aggregated value repeated in each trial within participant
  472. group_by(SEGA_ID) %>%
  473. mutate(stimulus_isi_duration_aggregated = mean(stimulus_isi_duration, na.rm = TRUE)) %>%
  474. ungroup()
  476. df_stimulus_isi_duration_aggregated <- et_erp_trial %>% # new data frame with 1 participant per row for anova
  477. group_by(SEGA_ID, group) %>%
  478. summarise(stimulus_isi_duration_aggregated = first(stimulus_isi_duration_aggregated),
  479. .groups = "drop")
  481. ggplot(df_stimulus_isi_duration_aggregated, aes(x = stimulus_isi_duration_aggregated)) +
  482. geom_histogram(binwidth = 0.2, fill = "grey", color = "white") +
  483. labs(x = "Duration (s)", y = "Count") +
  484. coord_cartesian(xlim = c(1.85, 3.4)) +
  485. facet_wrap(~ group) + # separate histogram per group
  486. theme_minimal() +
  487. theme(
  488. strip.text = element_text(face = "bold"))
  490. stats_stimulus_isi_duration_aggregated <- df_stimulus_isi_duration_aggregated %>%
  491. group_by(group) %>%
  492. summarise(
  493. mean_stimulus_isi_duration = mean(stimulus_isi_duration_aggregated, na.rm = TRUE),
  494. sd_stimulus_isi_duration = sd(stimulus_isi_duration_aggregated, na.rm = TRUE))
  495. stats_stimulus_isi_duration_aggregated
  497. anova_stimulus_isi_duration_aggregated <- aov(
  498. stimulus_isi_duration_aggregated ~ group, data = df_stimulus_isi_duration_aggregated)
  499. summary(anova_stimulus_isi_duration_aggregated)
  501. ggplot(df_stimulus_isi_duration_aggregated, aes(x = group, y = stimulus_isi_duration_aggregated, fill = group)) +
  502. geom_boxplot(outlier.shape = 16, alpha = 0.6) +
  503. geom_jitter(width = 0.2, alpha = 0.7, color = "black") +
  504. labs(
  505. x = NULL, # remioves x-axs label "group"
  506. y = "Stimulus_ISI_Duration (aggregated on particiopant level)") +
  507. theme_minimal() +
  508. theme(
  509. legend.position = "none",
  510. axis.text.x = element_text(face = "bold", size = "14"))
  512. # Power Analysis
  513. n_size<- 150 # sample size
  514. k_size<-4*100 # number of trials --> auditory oddball:
  515. subj <- factor(1:n_size)
  516. trial_id <- 1:k_size # remains integer variable
  517. condition <- c("before", "after")
  518. group <- c("ASD", "MHC", "CON")
  520. subj_full <- rep(subj, k_size)
  521. trial_full <- rep(trial_id, each=n_size)
  522. condition_full <- rep(condition, each=k_size*n_size/2)
  523. group_full <- rep(group, k_size*n_size/3)
  525. covars <- data.frame(id=subj_full, trial=trial_full, condition=condition_full, group=group_full)
  526. fixed <- c(0.1, 0.2, 0.2, 0.2, 0.1, 0.1) # effect sizes
  527. rand <- list(0.05) # random intercept
  528. res <- 1.2 # residual standard deviation
  530. model <- makeLmer(
  531. y ~ group*condition + (1|id), fixef=fixed, VarCorr=rand, sigma=res, data=covars) # create model
  532. summary(model) # model
  534. power_interaction <- powerSim(model, nsim=1000, test = fcompare(y ~ group + condition))
  535. print(power_interaction)
  536. power_curve_interaction <- powerCurve(model, test = fcompare(y ~ group + condition), along = "id")
  537. plot(power_curve_interaction)
  539. # EXPLORATORY (Reviewer's comment): Correlation SP2-P3a amplitude
  540. ## new df with mean P3a amplitude across conditions per participant
  541. df_SP2_P3a <- et_erp_subject %>%
  542. group_by(SEGA_ID) %>%
  543. summarise(
  544. mean_z_P3a_amplitude = mean(z_P3a_amplitude),
  545. SP2 = first(SP2)) %>% # just take the first value, is the same
  546. ungroup()
  549. corr_Sp2_P3a <- df_SP2_P3a %>%
  550. summarise(
  551. cor = cor(mean_z_P3a_amplitude, SP2, use = "complete.obs"),
  552. p_value = cor.test(mean_z_P3a_amplitude, SP2)$p.value,
  553. conf_low = cor.test(mean_z_P3a_amplitude, SP2)$conf.int[1],
  554. conf_high = cor.test(mean_z_P3a_amplitude, SP2)$conf.int[2]
  555. )
  556. print(corr_Sp2_P3a)
  558. ggplot(df_SP2_P3a, aes(x = SP2, y = mean_z_P3a_amplitude)) +
  559. geom_point(size = 3, color = "steelblue") +
  560. geom_smooth(method = "lm", color = "darkred", se = TRUE) +
  561. labs(
  562. x = "SP2",
  563. y = "Mean z_P3a amplitude",
  564. title = "Correlation between SP2 and Mean z_P3a amplitude"
  565. ) +
  566. annotate(
  567. "text",
  568. x = min(df_SP2_P3a$SP2),
  569. y = max(df_SP2_P3a$mean_z_P3a_amplitude),
  570. label = paste0("r = ", round(corr_test$estimate, 2),
  571. "\n95% CI [", round(corr_test$conf.int[1],2),
  572. ", ", round(corr_test$conf.int[2],2), "]\n",
  573. "p = ", round(corr_test$p.value, 3)),
  574. hjust = 0, vjust = 1, size = 4
  575. ) +
  576. theme_minimal(base_size = 14)
  579. # CORRELATION: BLOCK-BASELINE vs. TRIALBASELINE
  580. corr_block_trial_bl <- et_erp_subject %>%
  581. summarise(
  582. cor = cor(block_baseline_mean, rpd_low, use = "complete.obs"),
  583. p_value = cor.test(block_baseline_mean, rpd_low)$p.value,
  584. conf_low = cor.test(block_baseline_mean, rpd_low)$conf.int[1],
  585. conf_high = cor.test(block_baseline_mean, rpd_low)$conf.int[2]
  586. )
  587. print(corr_block_trial_bl)
  589. # Scatter plot with regression line and annotated correlation
  590. ggplot(et_erp_subject, aes(x = block_baseline_mean, y = rpd_low)) +
  591. geom_point(alpha = 0.6, color = "#0072B2", size = 3) + # points
  592. geom_smooth(method = "lm", se = TRUE, color = "#D55E00") + # regression line + CI
  593. annotate("text", x = Inf, y = -Inf,
  594. hjust = 1.1, vjust = -0.5,
  595. label = sprintf("r = %.2f, p < .001", corr_block_trial_bl$cor), # add correlation
  596. size = 5) +
  597. labs(
  598. x = "separate baseline phase",
  599. y = "BPS",
  600. title = "Correlation between separate baseline phase and BPS"
  601. ) +
  602. theme_minimal(base_size = 14)
  604. # RESULT 1: SEPR ON SUBJECT LEVEL
  605. lmm <- lmer(
  606. z_rpd ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  607. data = et_erp_subject)
  608. anova(lmm)
  609. r2_nakagawa(lmm)
  611. ## Post-hoc: stimulus
  612. contrast(emmeans(lmm, ~ stimulus), "pairwise")
  613. confint(contrast(emmeans(lmm, ~ stimulus), "pairwise"))
  615. # RESULT 2: BPS ON SUBJECT LEVEL
  616. lmm <- lmer(
  617. z_rpd_low ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  618. data = et_erp_subject)
  619. anova(lmm)
  620. r2_nakagawa(lmm)
  622. ## post-hoc: manipulation * group
  623. contrast(emmeans(lmm, ~ manipulation|group), method = "revpairwise")
  624. confint(contrast(emmeans(lmm, ~ manipulation|group), method = "revpairwise"))
  625. emmip(lmm, ~ manipulation | group, linearg = list(linetype = "blank"), CIs = T)
  627. ## Plot WTAS: manipulation * group
  628. emm_data <- emmip(lmm, ~ manipulation | group, CIs = TRUE, plotit = FALSE)
  629. ggplot(emm_data, aes(x = manipulation, y = yvar, color = group)) +
  630. geom_point(size = 3) + # Increased point size
  631. geom_errorbar(aes(ymin = LCL, ymax = UCL), width = 0.2, linewidth = 1) +
  632. facet_wrap(~ group) +
  633. labs(x = "Manipulation", y = "Estimated Marginal Mean") +
  634. theme_bw() +
  635. scale_color_manual(values = c("ASD" = "red", "CON" = "#045D5D", "MHC" = "#A05000")) +
  636. theme(plot.title = element_text(face = "bold"), legend.position = "none") +
  637. theme(axis.text.x = element_text(face = "bold")) +
  638. theme(axis.text.y = element_text(face = "bold")) +
  639. theme(text = element_text(size = 14))
  641. ### exploratory: custom contrast: group comparisons of BPS before manipulation
  642. emm <- emmeans(lmm, ~ group * manipulation)
  643. contrast(emm, method = list(
  644. "ASD.before vs CON.before" = c(-1, 1, 0, 0, 0, 0)), infer = T)
  645. contrast(emm, method = list(
  646. "ASD.before vs MHC.before" = c(-1, 0, 1, 0, 0, 0)), infer = T)
  647. contrast(emm, method = list(
  648. "CON.before vs MHC.before" = c(0, 0, 1, 0, 0, 0)), infer = T)
  650. ## Post-hoc: manipulation * block
  651. contrast(emmeans(lmm, ~ manipulation|block), "revpairwise")
  652. confint(contrast(emmeans(lmm, ~ manipulation|block), "revpairwise"))
  653. emmip(lmm, ~ manipulation | block, linearg = list(linetype = "blank"), CIs = T)
  655. ## Plot: manipulation x block interaction
  656. plot_data <- as.data.frame(emmeans(lmm, ~ block * manipulation))
  657. plot_data <- na.omit(plot_data)
  659. interaction_plot_BPS <- ggplot(plot_data) +
  660. geom_crossbar(aes(
  661. x = manipulation, y = emmean, color = block,
  662. ymin = emmean - 1 * SE, ymax = emmean + 1 * SE),
  663. alpha = 0.8,
  664. position = position_dodge(width = 0.7),
  665. width = 0.3) +
  666. geom_errorbar(aes(
  667. x = manipulation, y = emmean, color = block, fill = block,
  668. ymin = lower.CL, ymax = upper.CL),
  669. position = position_dodge(width = 0.7),
  670. width = 0.2) +
  671. scale_color_manual(values = c("blue", "orange")) +
  672. scale_fill_manual(values = c("blue", "orange")) +
  673. theme_bw() +
  674. labs(title = "Manipulation x block interaction on BPS",
  675. x = "Manipulation",
  676. y = "BPS [z]") +
  677. theme(plot.title = element_text(face = "bold"), legend.position = "none")
  678. print(interaction_plot_BPS)
  680. ## Plot BPS for WTAS 2025: manipulation x block interaction
  681. plot_data <- as.data.frame(emmeans(lmm, ~ block * manipulation))
  682. plot_data <- na.omit(plot_data)
  684. interaction_plot_BPS <- ggplot(plot_data) +
  685. geom_crossbar(aes(
  686. x = manipulation, y = emmean, color = block,
  687. ymin = emmean - 1 * SE, ymax = emmean + 1 * SE),
  688. alpha = 0.8,
  689. position = position_dodge(width = 0.9),
  690. width = 0.3,
  691. linewidth = 0.8) +
  692. geom_errorbar(aes(
  693. x = manipulation, y = emmean, color = block,
  694. ymin = lower.CL, ymax = upper.CL),
  695. position = position_dodge(width = 0.9),
  696. width = 0.3,
  697. linewidth = 0.8) +
  698. scale_color_manual(values = c("blue", "orange")) +
  699. scale_fill_manual(values = c("blue", "orange")) +
  700. theme_bw() +
  701. labs(title = "Manipulation x block interaction on BPS",
  702. x = "Manipulation",
  703. y = "BPS [z]") +
  704. theme(plot.title = element_text(face = "bold"), legend.position = "none") +
  705. theme(axis.text.x = element_text(face = "bold")) +
  706. theme(axis.text.y = element_text(face = "bold")) +
  707. theme(text = element_text(size = 14))
  708. print(interaction_plot_BPS)
  710. ## Plot BPS: manipulation x block interaction (SEGA_Paper)
  711. plot_data <- as.data.frame(emmeans(lmm, ~ block * manipulation))
  712. plot_data <- na.omit(plot_data)
  714. interaction_plot_BPS <- ggplot(plot_data) +
  715. geom_crossbar(aes(
  716. x = manipulation, y = emmean, color = block,
  717. ymin = emmean - 1 * SE, ymax = emmean + 1 * SE),
  718. alpha = 0.8,
  719. position = position_dodge(width = 0.9),
  720. width = 0.3,
  721. linewidth = 0.8) +
  722. geom_errorbar(aes(
  723. x = manipulation, y = emmean, color = block,
  724. ymin = lower.CL, ymax = upper.CL),
  725. position = position_dodge(width = 0.9),
  726. width = 0.3,
  727. linewidth = 0.8) +
  728. scale_color_manual(values = c("blue", "orange")) +
  729. scale_fill_manual(values = c("blue", "orange")) +
  730. theme_bw() +
  731. labs(x = "Manipulation", y = "BPS [z]") +
  732. theme(legend.position = "none") +
  733. theme(axis.text.x = element_text(face = "bold")) +
  734. theme(axis.text.y = element_text(face = "bold")) +
  735. theme(text = element_text(size = 18))
  736. print(interaction_plot_BPS)
  738. ## custom contrast: Manipulation effect in block 3? (before.reverse vs. after.forward)
  739. emm <- emmeans(lmm, ~ block * manipulation)
  740. contrast(emm, method = list(
  741. "reverse.before vs forward.after" = c(0, -1, 1, 0)), infer = T)
  742. emmip(emm, ~ block * manipulation, linearg = list(linetype = "blank"), CI = TRUE) +
  743. theme_bw() + labs(x = "task block", y = "Estimated Marginal Means")
  745. ### custom contrast: Manipulation effect in block 4? (before.reverse vs. after.reverse).
  746. emm <- emmeans(lmm, ~ block * manipulation)
  747. contrast(emm, method = list(
  748. "reverse.before vs reverse.after" = c(0, -1, 0, 1)), infer = T)
  749. emmip(emm, ~ block * manipulation, linearg = list(linetype = "blank"), CI = TRUE) +
  750. theme_bw() + labs(x = "task block", y = "Estimated Marginal Means")
  752. ### custom contrast: Habituation from block 1 to block 2?
  753. emm <- emmeans(lmm, ~ block * manipulation)
  754. contrast(emm, method = list(
  755. "reverse.before vs forward.before" = c(-1, 1, 0, 0)), infer = T)
  757. ### custom contrast: Habituation from block 3 to 4
  758. emm <- emmeans(lmm, ~ block * manipulation)
  759. contrast(emm, method = list(
  760. "reverse.after vs forward.after" = c(0, 0, -1, 1)), infer = T)
  762. # RESULT 3: MMN AMPLITUDE ON SUBJECT LEVEL
  763. lmm <- lmer(
  764. z_MMN_amplitude ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  765. data = et_erp_subject)
  766. anova(lmm)
  767. r2_nakagawa(lmm)
  769. ## Post-hoc: stimulus
  770. contrast(emmeans(lmm, ~ stimulus), "pairwise")
  771. confint(contrast(emmeans(lmm, ~ stimulus), "pairwise"))
  773. ## Post-hoc: manipulation * block
  774. contrast(emmeans(lmm, ~ manipulation|block), "revpairwise")
  775. confint(contrast(emmeans(lmm, ~ manipulation|block), "revpairwise"))
  776. emm <- emmeans(lmm, ~ block * manipulation)
  777. emmip(emm, ~ block * manipulation, linearg = list(linetype = "blank"), CI = TRUE) +
  778. theme_bw() + labs(x = "task block", y = "Estimated Marginal Means")
  780. ## Plot 1: manipulation x block interaction
  781. plot_data <- as.data.frame(emmeans(lmm, ~ block * manipulation))
  782. plot_data <- na.omit(plot_data)
  784. interaction_plot_MMN <- ggplot(plot_data) +
  785. geom_crossbar(aes(
  786. x = manipulation, y = emmean, color = block,
  787. ymin = emmean - 1 * SE, ymax = emmean + 1 * SE),
  788. alpha = 0.8,
  789. position = position_dodge(width = 0.7),
  790. width = 0.3) +
  791. geom_errorbar(aes(
  792. x = manipulation, y = emmean, color = block, fill = block,
  793. ymin = lower.CL, ymax = upper.CL),
  794. position = position_dodge(width = 0.7),
  795. width = 0.2) +
  796. scale_color_manual(values = c("blue", "orange")) +
  797. scale_fill_manual(values = c("blue", "orange")) +
  798. theme_bw() +
  799. labs(title = "Manipulation x block interaction on MMN amplitude",
  800. x = "Manipulation",
  801. y = "MMN_amplitude [z]") +
  802. theme(plot.title = element_text(face = "bold"), legend.position = "none")
  803. print(interaction_plot_MMN)
  805. legend <- cowplot::get_legend(interaction_plot_MMN + theme(legend.position="right"))
  806. legend_plot <- ggdraw() + draw_plot(legend)
  807. print(legend_plot)
  809. ## Slightly modified plot for WTAS 2025
  810. plot_data <- as.data.frame(emmeans(lmm, ~ block * manipulation))
  811. plot_data <- na.omit(plot_data)
  813. interaction_plot_MMN <- ggplot(plot_data) +
  814. geom_crossbar(aes(
  815. x = manipulation, y = emmean, color = block,
  816. ymin = emmean - 1 * SE, ymax = emmean + 1 * SE),
  817. alpha = 0.8,
  818. position = position_dodge(width = 0.9),
  819. width = 0.3,
  820. linewidth = 0.8) +
  821. geom_errorbar(aes(
  822. x = manipulation, y = emmean, color = block,
  823. ymin = lower.CL, ymax = upper.CL),
  824. position = position_dodge(width = 0.9),
  825. width = 0.3,
  826. linewidth = 0.8) +
  827. scale_color_manual(values = c("blue", "orange")) +
  828. scale_fill_manual(values = c("blue", "orange")) +
  829. theme_bw() +
  830. labs(title = "Manipulation x block interaction on MMN amplitude",
  831. x = "Manipulation",
  832. y = "MMN_amplitude [z]") +
  833. theme(plot.title = element_text(face = "bold"), legend.position = "none") +
  834. theme(axis.text.x = element_text(face = "bold")) +
  835. theme(axis.text.y = element_text(face = "bold")) +
  836. theme(text = element_text(size = 14))
  837. print(interaction_plot_MMN)
  839. ## Plot MMN amplitude (SEGA-Paper)
  840. plot_data <- as.data.frame(emmeans(lmm, ~ block * manipulation))
  841. plot_data <- na.omit(plot_data)
  843. interaction_plot_MMN <- ggplot(plot_data) +
  844. geom_crossbar(aes(
  845. x = manipulation, y = emmean, color = block,
  846. ymin = emmean - 1 * SE, ymax = emmean + 1 * SE),
  847. alpha = 0.8,
  848. position = position_dodge(width = 0.9),
  849. width = 0.3,
  850. linewidth = 0.8) +
  851. geom_errorbar(aes(
  852. x = manipulation, y = emmean, color = block,
  853. ymin = lower.CL, ymax = upper.CL),
  854. position = position_dodge(width = 0.9),
  855. width = 0.3,
  856. linewidth = 0.8) +
  857. scale_color_manual(values = c("blue", "orange")) +
  858. scale_fill_manual(values = c("blue", "orange")) +
  859. theme_bw() +
  860. labs(x = "Manipulation", y = "MMN_amplitude [z]") +
  861. theme(legend.position = "none") +
  862. theme(axis.text.x = element_text(face = "bold")) +
  863. theme(axis.text.y = element_text(face = "bold")) +
  864. theme(text = element_text(size = 18))
  865. print(interaction_plot_MMN)
  867. ## Custom contrast: Manipulation effect in block 3? (block 2 vs. 3).
  868. emm <- emmeans(lmm, ~ block * manipulation)
  869. contrast(emm, method = list(
  870. "reverse.before vs forward.after" = c(0, -1, 1, 0)), infer = T)
  872. ## Custom contrast: Manipulation effect in block 4? (block 2 vs 4).
  873. emm <- emmeans(lmm, ~ block * manipulation)
  874. contrast(emm, method = list(
  875. "reverse.before vs reverse.after" = c(0, -1, 0, 1)), infer = T)
  877. ## Correlation with age
  878. corr_by_group <- et_erp_subject %>%
  879. group_by(group) %>%
  880. summarise(cor = cor(age, z_MMN_latency, use = "complete.obs"),
  881. p_value = cor.test(age, z_MMN_latency)$p.value)
  882. print(corr_by_group)
  884. # RESULT 4: MMN LATENCY ON SUBJECT LEVEL
  885. lmm <- lmer(
  886. z_MMN_latency ~ stimulus * manipulation * group * block + (1|SEGA_ID) + gender + age,
  887. data = et_erp_subject)
  888. anova(lmm)
  889. r2_nakagawa(lmm)
  891. ## Correlation with age
  892. corr_by_group <- et_erp_subject %>%
  893. group_by(group) %>%
  894. summarise(cor = cor(age, z_MMN_amplitude, use = "complete.obs"),
  895. p_value = cor.test(age, z_MMN_latency)$p.value)
  896. print(corr_by_group)
  898. # RESULT 5: MMN DIFF AMPLITUDE ON SUBJECT LEVEL
  899. lmm <- lmer(
  900. z_MMN_diff_amplitude ~ manipulation * group * block + (1|SEGA_ID) + gender + age,
  901. data = MMN_diff)
  902. anova(lmm)
  903. r2_nakagawa(lmm)
  905. # RESULT 6: MMN DIFF LATENCY ON SUBJECT LEVEL
  906. lmm <- lmer(
  907. z_MMN_diff_latency ~ manipulation * group * block + (1|SEGA_ID) + gender + age,
  908. data = MMN_diff)
  909. anova(lmm)
  910. r2_nakagawa(lmm)
  912. # RESULT 7: P3A AMPLITUDE ON SUBJECT LEVEL
  913. lmm <- lmer(
  914. z_P3a_amplitude ~ stimulus * manipulation * group * block + (1|SEGA_ID) + gender + age,
  915. data = et_erp_subject)
  916. anova(lmm)
  917. r2_nakagawa(lmm)
  919. ## Post-hoc: stimulus * group * block
  920. emm <- emmeans(lmm, ~ group * stimulus * block)
  921. stimulus_diff <- contrast(emm, method = "pairwise", by = c("group", "block"))
  922. stimulus_diff_emm <- as.emm_list(stimulus_diff)
  923. group_diff <- contrast(stimulus_diff_emm, method = "pairwise", by = "block")
  924. summary(group_diff)
  925. confint(group_diff)
  927. ## WTAS: Plot with stimulus * group
  928. emm_data <- emmip(lmm, ~ stimulus | group, CIs = TRUE, plotit = FALSE)
  929. ggplot(emm_data, aes(x = xvar, y = yvar, color = stimulus)) +
  930. geom_point(size = 3) +
  931. geom_errorbar(aes(ymin = LCL, ymax = UCL), width = 0.2, size = 1) +
  932. facet_wrap(~ group) +
  933. labs(x = "Stimulus", y = "Estimated Marginal Mean") +
  934. theme_bw() +
  935. scale_color_manual(values = c("Oddball" = "red", "Standard" = "black")) +
  936. theme(plot.title = element_text(face = "bold"), legend.position = "none") +
  937. theme(axis.text.x = element_text(face = "bold")) +
  938. theme(axis.text.y = element_text(face = "bold")) +
  939. theme(text = element_text(size = 14))
  941. ## Post-hoc: manipulation * block
  942. emmeans(lmm, ~ manipulation|block)
  943. contrast(emmeans(lmm, ~ manipulation|block), "revpairwise")
  944. confint(contrast(emmeans(lmm, ~ manipulation|block), "revpairwise"))
  945. emm <- emmeans(lmm, ~ block * manipulation)
  946. emmip(emm, ~ block * manipulation, linearg = list(linetype = "blank"), CI = TRUE) +
  947. theme_bw() + labs(x = "task block", y = "Estimated Marginal Means")
  949. ## Plot: manipulation * block
  950. plot_data <- as.data.frame(emmeans(lmm, ~ block * manipulation))
  951. plot_data <- na.omit(plot_data)
  953. interaction_plot_P3a <- ggplot(plot_data, aes(x = manipulation, y = emmean, group = block, color = block)) +
  954. geom_crossbar(aes(ymin = emmean - 1 * SE, ymax = emmean + 1 * SE),
  955. alpha = 0.8,
  956. position = position_dodge(width = 0.9),
  957. width = 0.3) +
  958. geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
  959. position = position_dodge(width = 0.9),
  960. width = 0.2) +
  961. scale_color_manual(values = c("blue", "orange")) +
  962. scale_fill_manual(values = c("blue", "orange")) +
  963. scale_x_discrete(expand = expansion(mult = c(0.5, 0.5))) +
  964. theme_bw() +
  965. labs(title = "Manipulation x block interaction on P3a amplitude",
  966. x = "Manipulation",
  967. y = "P3a_amplitude [z]") +
  968. theme(plot.title = element_text(face = "bold"), legend.position = "none")
  969. print(interaction_plot_P3a)
  971. ## Slightly modified for WTAS 2025
  972. plot_data <- as.data.frame(emmeans(lmm, ~ block * manipulation))
  973. plot_data <- na.omit(plot_data)
  974. interaction_plot_P3a <- ggplot(plot_data, aes(x = manipulation, y = emmean, group = block, color = block)) +
  975. geom_crossbar(aes(ymin = emmean - 1 * SE, ymax = emmean + 1 * SE),
  976. alpha = 0.8,
  977. position = position_dodge(width = 0.9),
  978. width = 0.3,
  979. size = 0.8) +
  980. geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
  981. position = position_dodge(width = 0.9),
  982. width = 0.2,
  983. size = 0.8) +
  984. scale_color_manual(values = c("blue", "orange")) +
  985. scale_fill_manual(values = c("blue", "orange")) +
  986. scale_x_discrete(expand = expansion(mult = c(0.5, 0.5))) +
  987. theme_bw() +
  988. labs(title = "Manipulation x block interaction on P3a amplitude",
  989. x = "Manipulation",
  990. y = "P3a_amplitude [z]") +
  991. theme(plot.title = element_text(face = "bold"), legend.position = "none") +
  992. theme(axis.text.x = element_text(face = "bold")) +
  993. theme(axis.text.y = element_text(face = "bold")) +
  994. theme(text = element_text(size = 14))
  995. print(interaction_plot_P3a)
  997. legend <- cowplot::get_legend(interaction_plot_MMN + theme(legend.position="right"))
  998. legend_plot <- ggdraw() + draw_plot(legend)
  999. print(legend_plot)
  1001. ## Plot P3a: Manipulation effect (SEGA Paper)
  1002. plot_data <- as.data.frame(emmeans(lmm, ~ block * manipulation))
  1003. plot_data <- na.omit(plot_data)
  1004. interaction_plot_P3a <- ggplot(plot_data, aes(x = manipulation, y = emmean, group = block, color = block)) +
  1005. geom_crossbar(aes(ymin = emmean - 1 * SE, ymax = emmean + 1 * SE),
  1006. alpha = 0.8,
  1007. position = position_dodge(width = 0.9),
  1008. width = 0.3,
  1009. size = 0.8) +
  1010. geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
  1011. position = position_dodge(width = 0.9),
  1012. width = 0.2,
  1013. size = 0.8) +
  1014. scale_color_manual(values = c("blue", "orange")) +
  1015. scale_fill_manual(values = c("blue", "orange")) +
  1016. scale_x_discrete(expand = expansion(mult = c(0.5, 0.5))) +
  1017. theme_bw() +
  1018. labs(x = "Manipulation", y = "P3a_amplitude [z]") +
  1019. theme(legend.position = "none") +
  1020. theme(axis.text.x = element_text(face = "bold")) +
  1021. theme(axis.text.y = element_text(face = "bold")) +
  1022. theme(text = element_text(size = 18))
  1023. print(interaction_plot_P3a)
  1025. ## custom contrast: Manipulation effect in block 3? (before.reverse vs. after.forward)
  1026. emm <- emmeans(lmm, ~ block * manipulation)
  1027. contrast(emm, method = list(
  1028. "reverse.before vs forward.after" = c(0, -1, 1, 0)), infer = T)
  1030. ## custom contrast: Manipulation effect in block 4? (before.reverse vs. after.reverse)
  1031. emm <- emmeans(lmm, ~ block * manipulation)
  1032. contrast(emm, method = list(
  1033. "reverse.before vs reverse.after" = c(0, 1, 0, -1)), infer = T)
  1035. # RESULT 8: P3A LATENCY ON SUBJECT LEVEL
  1036. lmm <- lmer(
  1037. z_P3a_latency ~ stimulus * manipulation * group * block + (1|SEGA_ID) + gender + age,
  1038. data = et_erp_subject)
  1039. anova(lmm)
  1040. r2_nakagawa(lmm)
  1042. ## Post-hoc: stimulus
  1043. contrast(emmeans(lmm, ~ stimulus), "pairwise")
  1044. confint(contrast(emmeans(lmm, ~ stimulus), "pairwise"))
  1046. # RESULT 9: P3B AMPLITUDE ON SUBJECT LEVEL
  1047. lmm <- lmer(
  1048. z_P3b_amplitude ~ stimulus * manipulation * group * block + (1|SEGA_ID) + gender + age,
  1049. data = et_erp_subject)
  1050. anova(lmm)
  1051. r2_nakagawa(lmm)
  1053. ## Post-hoc: stimulus
  1054. contrast(emmeans(lmm, ~ stimulus), "pairwise")
  1055. confint(contrast(emmeans(lmm, ~ stimulus), "pairwise"))
  1057. ## Post-hoc: manipulation
  1058. contrast(emmeans(lmm, ~ manipulation), "pairwise")
  1059. confint(contrast(emmeans(lmm, ~ manipulation), "pairwise"))
  1060. emmip(lmm, ~ manipulation, linearg = list(linetype = "blank"), CI = T)
  1063. # RESULT 10: P3B LATENCY ON SUBJECT LEVEL
  1064. lmm <- lmer(
  1065. z_P3b_latency ~ stimulus * manipulation * group * block + (1|SEGA_ID) + gender + age,
  1066. data = et_erp_subject)
  1067. anova(lmm)
  1068. r2_nakagawa(lmm)
  1070. # Post-hoc: stimulus * manipulation * group
  1071. contrast(emmeans(lmm, ~ group|manipulation*stimulus), "pairwise")
  1072. confint(contrast(emmeans(lmm, ~ group|manipulation*stimulus), "pairwise"))
  1073. emmip(lmm, ~ group|manipulation*stimulus, CI = T)
  1074. contrast(emmeans(lmm, ~ manipulation|group*stimulus), "pairwise")
  1076. # RESULT 11: ASSOCIATIONS OF PUPILLOMETRIC MEASURES
  1077. lmm <- lmer(
  1078. z_rpd ~ z_rpd_low * stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  1079. data = et_erp_subject)
  1080. anova(lmm)
  1081. r2_nakagawa(lmm)
  1083. ## Post-hoc: BPS
  1084. emtrends(lmm, ~ z_rpd_low, var = "z_rpd_low")
  1085. summary(emtrends(lmm, ~ z_rpd_low, var = "z_rpd_low"), infer = T)
  1087. # RESULT 12: ASSOCIATIONS: PUPIL DATA- ERPs ON SUBJECT LEVEL
  1088. ## MMN amplitude
  1089. lmm <- lmer(
  1090. z_MMN_amplitude ~ z_rpd * z_rpd_low * stimulus * manipulation * group + (1|SEGA_ID) + age + gender,
  1091. data = et_erp_subject)
  1092. anova(lmm)
  1093. r2_nakagawa(lmm)
  1095. ## Post-hoc: SEPR * BPS * stimulus
  1096. emt <- emtrends(lmm, ~ z_rpd_low * z_rpd | stimulus, var = 'z_rpd', at = list(z_rpd_low = c(-2,-1,0,1,2)))
  1097. summary(emt, infer = T)
  1099. ## P3a amplitude
  1100. lmm <- lmer(
  1101. z_P3a_amplitude ~ z_rpd * z_rpd_low * stimulus * manipulation * group + (1|SEGA_ID) + age + gender,
  1102. data = et_erp_subject)
  1103. anova(lmm)
  1104. r2_nakagawa(lmm)
  1106. ### post-hoc: SEPR x manipulation x group x stimulus
  1107. emt <- emtrends(lmm, ~ z_rpd * manipulation | group * stimulus, var = 'z_rpd')
  1108. summary(emt, infer = T)
  1110. ### post-hoc: BPS x stimulus x group
  1111. emt <- emtrends(lmm, ~ z_rpd_low * group | stimulus, var = 'z_rpd_low', at = list(z_rpd_low = c(-2,-1,0,1,2)))
  1112. summary(emt, infer = T)
  1114. lmm <- lmer(
  1115. z_P3b_amplitude ~ z_rpd * z_rpd_low * stimulus * manipulation * group + (1|SEGA_ID) + age + gender,
  1116. data = et_erp_subject)
  1117. anova(lmm)
  1118. r2_nakagawa(lmm)
  1120. crl <- cor(et_erp_subject[et_erp_subject$stimulus == "Oddball", c(
  1121. "z_rpd_low",
  1122. "z_rpd",
  1123. "z_MMN_amplitude",
  1124. "z_P3a_amplitude",
  1125. "z_P3b_amplitude"
  1126. )],
  1127. use="complete.obs")
  1128. crl
  1129. corrplot::corrplot(
  1130. crl,
  1131. method = "circle",
  1132. title = "Association between pupil data and ERPs in oddball trials",
  1133. mar=c(0,0,1,0))
  1135. crl <- cor(et_erp_subject[et_erp_subject$stimulus == "Standard", c(
  1136. "z_rpd_low",
  1137. "z_rpd",
  1138. "z_MMN_amplitude",
  1139. "z_P3a_amplitude",
  1140. "z_P3b_amplitude"
  1141. )],
  1142. use="complete.obs")
  1143. crl
  1144. corrplot::corrplot(
  1145. crl,
  1146. method = "circle",
  1147. title = "Association between pupil data and ERPs in standard trials",
  1148. mar=c(0,0,1,0))
  1150. # DATA ANALYSIS ON TRIAL LEVEL ####
  1151. # Distributions of dependent variables
  1152. par(mfrow = c(2,2), mar = c(4, 4, 2, 1))
  1153. hist(et_erp_trial$z_rpd,
  1154. main = "Distribution of rpd (500-1500 ms)",
  1155. xlab = "z_rpd",
  1156. xlim = c(-6, 6),
  1157. breaks = 200)
  1158. hist(et_erp_trial$z_rpd_low,
  1159. main = "Distribution of rpd_low (0-250 ms)",
  1160. xlab = "rpd_low",
  1161. xlim = c(-4, 6),
  1162. breaks = 200)
  1163. hist(et_erp_trial$z_MMN_amplitude,
  1164. main = "Distribution of MMN amplitude (100-150 ms)",
  1165. xlab = "MMN amplitude (FC1, FC2, FCz, Fz)",
  1166. xlim = c(-6, 6),
  1167. breaks = 200)
  1168. hist(et_erp_trial$z_P3a_amplitude,
  1169. main = "Distribution of P3a amplitude (150-250 ms)",
  1170. xlab = "P3a amplitude (Cz, FCz)",
  1171. xlim = c(-4, 6),
  1172. breaks = 200)
  1173. par(mfrow = c(1, 1))
  1175. # RESULT 12: SEPR ON TRIAL LEVEL
  1176. lmm <- lmer(
  1177. z_rpd ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  1178. data = et_erp_trial)
  1179. anova(lmm)
  1180. r2_nakagawa(lmm)
  1183. ## Post-hoc: stimulus
  1184. contrast(emmeans(lmm, ~ stimulus), "pairwise")
  1185. confint(contrast(emmeans(lmm, ~ stimulus), "pairwise"))
  1187. # EXPLORATORY: SEPR OVER TASK
  1188. lmm <- lmer(
  1189. z_rpd ~ stimulus * manipulation * group * block * oddball_trial_counter + (1|SEGA_ID),
  1190. data = et_erp_trial)
  1191. anova(lmm)
  1192. r2_nakagawa(lmm)
  1194. emt <- emtrends(lmm, pairwise ~ manipulation, var = "oddball_trial_counter")
  1195. summary(emt, infer = T)
  1197. emt <- emtrends(lmm,~ stimulus * group * block, var = 'oddball_trial_counter')
  1198. summary(emt, infer = T)
  1200. et_erp_trial$binned_trial_number <- cut(
  1201. et_erp_trial$oddball_trial_counter, c(
  1202. 1, 103, 206, 309, 412))
  1203. ggplot(et_erp_trial[et_erp_trial$stimulus == "oddball", ]) + geom_boxplot(aes(binned_trial_number, z_rpd)) +
  1204. facet_grid(cols = vars(group)) +
  1205. ggtitle("SEPR in oddball trials during the task, \nsplit by group") +
  1206. theme(plot.title = element_text(face = "bold", size = 14))
  1208. # EXPLORATORY: SEPR OVER BLOCK
  1209. lmm <- lmer(
  1210. z_rpd ~ stimulus * manipulation * group * block * trial_number_in_block + (1|SEGA_ID),
  1211. data = et_erp_trial)
  1212. anova(lmm)
  1213. r2_nakagawa(lmm)
  1215. emt <- emtrends(lmm, pairwise ~ manipulation, var = "trial_number_in_block")
  1216. summary(emt, infer = T)
  1218. emt <- emtrends(lmm,~ stimulus * group * block, var = 'trial_number_in_block')
  1219. summary(emt, infer = T)
  1221. bin_size <- 20
  1222. df_trial <- df_trial %>%
  1223. mutate(trial_number_in_block_binned = factor(trial_number_in_block%/%bin_size*20))
  1224. ggplot(
  1225. df_trial[df_trial$phase %in% c("oddball_block", "oddball_block_rev") & df_trial$stimulus == "oddball", ],
  1226. aes(x = trial_number_in_block_binned,
  1227. y = scale(rpd))) +
  1228. geom_boxplot() +
  1229. facet_grid(
  1230. rows = vars(block),
  1231. cols = vars(manipulation)) +
  1232. theme_bw() +
  1233. ggtitle("SEPR in oddball trials during blocks, \nsplit by forward + reverse and before + after manipulation") +
  1234. theme(plot.title = element_text(face = "bold", size = 14))
  1236. et_erp_trial$binned_trial_number <- cut(
  1237. et_erp_trial$trial_number_in_block, c(
  1238. 0, 20, 40, 60, 80, 100))
  1239. ggplot(et_erp_trial[et_erp_trial$stimulus == "oddball", ]) + geom_boxplot(aes(binned_trial_number, z_rpd)) +
  1240. facet_grid(cols = vars(group)) +
  1241. ggtitle("SEPR in oddball trials during blocks, split by group") +
  1242. theme(plot.title = element_text(face = "bold", size = 14))
  1244. # RESULT 13: BPS ON TRIAL LEVEL
  1245. lmm <- lmer(
  1246. z_rpd_low ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  1247. data = et_erp_trial)
  1248. anova(lmm)
  1249. r2_nakagawa(lmm)
  1251. ## Post-hoc: manipulation * group
  1252. contrast(emmeans(lmm, ~ manipulation|group), "revpairwise")
  1253. confint(contrast(emmeans(lmm, ~ manipulation|group), "revpairwise"))
  1254. emm <- emmeans(lmm, ~ manipulation * group)
  1255. emmip(emm, ~ manipulation|group, linearg = list(linetype = "blank"), CI = TRUE) +
  1256. theme_bw() + labs(x = "task block", y = "Estimated Marginal Means")
  1258. ## Post-hoc: manipulation * block
  1259. contrast(emmeans(lmm, ~ manipulation|block), "revpairwise")
  1260. confint(contrast(emmeans(lmm, ~ manipulation|block), "revpairwise"))
  1261. emm <- emmeans(lmm, ~ block * manipulation)
  1262. emmip(emm, ~ block * manipulation, linearg = list(linetype = "blank"), CI = TRUE) +
  1263. theme_bw() + labs(x = "task block", y = "Estimated Marginal Means")
  1265. ## custom contrast: manipulation effect in block 3? (block 2 vs. 3).
  1266. emm <- emmeans(lmm, ~ block * manipulation)
  1267. contrast(emm, method = list(
  1268. "reverse.before vs forward.after" = c(0, -1, 1, 0)), infer = T)
  1270. ## custom contrast: manipulation effect in block 4? (block 2 vs 4).
  1271. emm <- emmeans(lmm, ~ block * manipulation)
  1272. contrast(emm, method = list(
  1273. "reverse.before vs reverse.after" = c(0, -1, 0, 1)), infer = T)
  1275. ## Post-hoc: group * block
  1276. contrast(emmeans(lmm, ~ block|group), "pairwise")
  1277. confint(contrast(emmeans(lmm, ~ block|group), "pairwise"))
  1278. emm <- emmeans(lmm, ~ block | group)
  1279. emmip(emm, ~ block | group, linearg = list(linetype = "blank"), CI = TRUE) +
  1280. theme_bw() + labs(x = "task block", y = "Estimated Marginal Means")
  1282. # EXPLORATORY: BPS OVER TASK
  1283. lmm <- lmer(
  1284. z_rpd_low ~ stimulus * manipulation * group * block * oddball_trial_counter + (1|SEGA_ID),
  1285. data = et_erp_trial)
  1286. anova(lmm)
  1287. r2_nakagawa(lmm)
  1289. emt <- emtrends(lmm,~ manipulation * group * block, var = 'oddball_trial_counter')
  1290. summary(emt, infer = T)
  1292. et_erp_trial$binned_trial_number <- cut(
  1293. et_erp_trial$oddball_trial_counter, c(
  1294. 1, 103, 206, 309, 412))
  1295. ggplot(et_erp_trial[et_erp_trial$stimulus == "standard", ]) + geom_boxplot(aes(binned_trial_number, z_rpd_low)) +
  1296. facet_grid(cols = vars(group)) +
  1297. ggtitle("BPS in standard trials during the task, split by group") +
  1298. theme(plot.title = element_text(face = "bold", size = 14))
  1300. # EXPLORATORY: BPS OVER BLOCK
  1301. lmm <- lmer(
  1302. z_rpd_low ~ stimulus * manipulation * group * block * trial_number_in_block + (1|SEGA_ID),
  1303. data = et_erp_trial)
  1304. anova(lmm)
  1305. r2_nakagawa(lmm)
  1307. emt <- emtrends(lmm,~ manipulation * group * block, var = 'trial_number_in_block')
  1308. summary(emt, infer = T)
  1310. bin_size <- 20
  1311. df_trial <- df_trial %>%
  1312. mutate(trial_number_in_block_binned = factor(trial_number_in_block%/%bin_size*20))
  1313. ggplot(
  1314. df_trial[df_trial$phase %in% c("oddball_block", "oddball_block_rev") & df_trial$stimulus == "standard", ],
  1315. aes(x = trial_number_in_block_binned,
  1316. y = scale(rpd_low))) +
  1317. geom_boxplot() +
  1318. facet_grid(
  1319. rows = vars(block),
  1320. cols = vars(manipulation)) +
  1321. theme_bw() +
  1322. ggtitle("BPS in standard trials during blocks, \nsplit by forward + reverse and before + after manipulation") +
  1323. theme(plot.title = element_text(face = "bold", size = 14))
  1325. # RESULT 8: MMN AMPLITUDE ON TRIAL LEVEL
  1326. lmm <- lmer(
  1327. z_MMN_amplitude ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  1328. data = et_erp_trial)
  1329. anova(lmm)
  1330. r2_nakagawa(lmm)
  1332. ## Post-hoc: stimulus
  1333. contrast(emmeans(lmm, ~ stimulus), "pairwise")
  1334. confint(contrast(emmeans(lmm, ~ stimulus), "pairwise"))
  1336. # RESULT 9: MMN LATENCY ON TRIAL LEVEL
  1337. lmm <- lmer(
  1338. z_MMN_latency ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  1339. data = et_erp_trial)
  1340. anova(lmm)
  1341. r2_nakagawa(lmm)
  1343. ## Post-hoc: stimulus
  1344. contrast(emmeans(lmm, ~ stimulus), "pairwise")
  1345. confint(contrast(emmeans(lmm, ~ stimulus), "pairwise"))
  1347. ## Post-hoc: manipulation * block
  1348. contrast(emmeans(lmm, ~ manipulation|block), "revpairwise")
  1349. confint(contrast(emmeans(lmm, ~ manipulation|block), "revpairwise"))
  1350. emm <- emmeans(lmm, ~ block * manipulation)
  1351. emmip(emm, ~ block * manipulation, linearg = list(linetype = "blank"), CI = TRUE) +
  1352. theme_bw() + labs(x = "task block", y = "Estimated Marginal Means")
  1354. ## Plot Manipulation effect on MMN Latency
  1355. plot_data <- as.data.frame(emmeans(lmm, ~ block * manipulation))
  1356. plot_data <- na.omit(plot_data)
  1358. interaction_plot_MNN_latency <- ggplot(plot_data, aes(x = manipulation, y = emmean, group = block, color = block)) +
  1359. geom_crossbar(aes(ymin = emmean - 1 * SE, ymax = emmean + 1 * SE),
  1360. alpha = 0.8,
  1361. position = position_dodge(width = 0.7),
  1362. width = 0.3) +
  1363. geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL),
  1364. position = position_dodge(width = 0.7),
  1365. width = 0.2) +
  1366. scale_color_manual(values = c("blue", "orange")) +
  1367. scale_fill_manual(values = c("blue", "orange")) +
  1368. theme_bw() +
  1369. labs(title = "Manipulation x block interaction on MMN latency",
  1370. x = "Manipulation",
  1371. y = "MMN latency [z]") +
  1372. theme(plot.title = element_text(face = "bold"), legend.position = "none")
  1373. print(interaction_plot_MNN_latency)
  1375. ## Slightly modified plot for WTAS 2025
  1376. plot_data <- as.data.frame(emmeans(lmm, ~ block * manipulation))
  1377. plot_data <- na.omit(plot_data)
  1379. interaction_plot_MNN_latency <- ggplot(plot_data, aes(x = manipulation, y = emmean, group = block, color = block)) +
  1380. geom_crossbar(aes(ymin = emmean - 1 * SE, ymax = emmean + 1 * SE),
  1381. alpha = 0.8,
  1382. position = position_dodge(width = 0.9),
  1383. width = 0.3,
  1384. linewidth = 0.8) +
  1385. geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL),
  1386. position = position_dodge(width = 0.9),
  1387. width = 0.2,
  1388. linewidth = 0.8) +
  1389. scale_color_manual(values = c("blue", "orange")) +
  1390. scale_fill_manual(values = c("blue", "orange")) +
  1391. theme_bw() +
  1392. labs(title = "Manipulation x block interaction on MMN latency",
  1393. x = "Manipulation",
  1394. y = "MMN latency [z]") +
  1395. theme(plot.title = element_text(face = "bold"), legend.position = "none") +
  1396. theme(axis.text = element_text(face = "bold")) +
  1397. theme(text = element_text(size = 14))
  1398. print(interaction_plot_MNN_latency)
  1400. grid.arrange(interaction_plot_BPS,
  1401. interaction_plot_MMN,
  1402. interaction_plot_P3a)
  1404. # RESULT 10: P3A AMPLITUDE ON TRIAL LEVEL
  1405. lmm <- lmer(
  1406. z_P3a_amplitude ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  1407. data = et_erp_trial)
  1408. anova(lmm)
  1409. r2_nakagawa(lmm)
  1411. ## Post-hoc: stimulus * group * block
  1412. emm <- emmeans(lmm, ~ group * stimulus * block)
  1413. stimulus_diff <- contrast(emm, method = "pairwise", by = c("group", "block"))
  1414. stimulus_diff_emm <- as.emm_list(stimulus_diff)
  1415. group_diff <- contrast(stimulus_diff_emm, method = "pairwise", by = "block")
  1416. summary(group_diff)
  1417. confint(group_diff)
  1419. # Post-hoc: manipulation * block
  1420. contrast(emmeans(lmm, ~ manipulation|block), "pairwise")
  1421. confint(contrast(emmeans(lmm, ~ manipulation|block), "pairwise"))
  1422. emm <- emmeans(lmm, ~ block * manipulation)
  1423. emmip(emm, ~ block * manipulation, linearg = list(linetype = "blank"), CI = TRUE) +
  1424. theme_bw() + labs(x = "task block", y = "Estimated Marginal Means")
  1426. ## custom contrast: Manipulation effect in block 3? (block 2 vs. 3).
  1427. emm <- emmeans(lmm, ~ block * manipulation)
  1428. contrast(emm, method = list(
  1429. "reverse.before vs forward.after" = c(0, -1, 1, 0)), infer = T)
  1431. ## custom contrast: Manipulation effect in block 4? (block 2 vs 4).
  1432. contrast(emm, method = list(
  1433. "reverse.before vs reverse.after" = c(0, 1, 0, -1)), infer = T)
  1435. # RESULT 11: P3A LATENCY ON TRIAL LEVEL
  1436. lmm <- lmer(
  1437. z_P3a_latency ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  1438. data = et_erp_trial)
  1439. anova(lmm)
  1440. r2_nakagawa(lmm)
  1442. # Post-hoc: stimulus
  1443. contrast(emmeans(lmm, ~ stimulus), "pairwise")
  1444. confint(contrast(emmeans(lmm, ~ stimulus), "pairwise"))
  1446. # RESULT 14: ASSOCIATION BETWEEN PUPILLOMETRIC MEASURES
  1447. ggscatter(et_erp_trial[et_erp_trial$stimulus == "oddball", ],
  1448. x = "z_rpd_low",
  1449. y = "z_rpd",
  1450. add = "reg.line",
  1451. conf.int = TRUE,
  1452. cor.coef = TRUE,
  1453. cor.method = "pearson",
  1454. xlab = "scaled BPS (rpd_low)",
  1455. ylab = "scaled SEPR (rpd)",
  1456. title = "Association between SEPR + BPS in oddball trials")
  1458. ggscatter(et_erp_trial[et_erp_trial$stimulus == "standard", ],
  1459. x = "z_rpd_low",
  1460. y = "z_rpd",
  1461. add = "reg.line",
  1462. conf.int = TRUE,
  1463. cor.coef = TRUE,
  1464. cor.method = "pearson",
  1465. xlab = "scaled BPS (rpd_low)",
  1466. ylab = "scaled SEPR (rpd)",
  1467. title = "Association between SEPR + BPS in standard trials")
  1469. # Associations of the pupillometric measures
  1470. lmm <- lmer(
  1471. z_rpd ~ z_rpd_low * stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  1472. data = et_erp_trial)
  1473. anova(lmm)
  1474. r2_nakagawa(lmm)
  1476. ### post-hoc: BPS
  1477. summary(emtrends(lmm, ~ z_rpd_low, var = "z_rpd_low"), infer = T)
  1479. ### post-hoc: BPS x manipulation x block
  1480. emt <- emtrends(lmm, ~ manipulation|block, var = "z_rpd_low")
  1481. contrast(emt, "revpairwise")
  1482. confint(contrast(emt, "revpairwise"))
  1484. ### EXPLORATORY (Reviewer's comment, 2025-10-27)
  1485. # lmer()-, but not poly()-function can handle NAs in the predictors, so
  1486. # et_erp_trial_for_model_fit excludes NAs for these exploratory model fit comparisons.
  1487. et_erp_trial_for_model_fit <- na.omit(et_erp_trial[ , c(
  1488. "z_rpd",
  1489. "z_rpd_low",
  1490. "stimulus",
  1491. "manipulation",
  1492. "group",
  1493. "block",
  1494. "SEGA_ID",
  1495. "age",
  1496. "gender")])
  1498. # linear
  1499. lmm_linear_fit <- lmer(
  1500. z_rpd ~ z_rpd_low * stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  1501. data = et_erp_trial_for_model_fit, REML = F)
  1502. anova(lmm_linear_fit)
  1503. r2_nakagawa(lmm_linear_fit)
  1505. # quadratic
  1506. lmm_quadratic_fit <- lmer(
  1507. z_rpd ~ poly(z_rpd_low, 2) * stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  1508. data = et_erp_trial_for_model_fit, REML = F)
  1509. anova(lmm_quadratic_fit)
  1510. r2_nakagawa(lmm_quadratic_fit)
  1511. # cubic
  1512. lmm_cubic_fit <- lmer(
  1513. z_rpd ~ poly(z_rpd_low, 3) * stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender,
  1514. data = et_erp_trial_for_model_fit, REML = F)
  1515. anova(lmm_cubic_fit)
  1516. r2_nakagawa(lmm_cubic_fit)
  1518. table_model_compare <- anova(lmm_linear_fit, lmm_quadratic_fit, lmm_cubic_fit)
  1520. table_model_compare <- cbind(
  1521. c("linear_fit", "quadratic_fit", "cubic_fit"),
  1522. table_model_compare)
  1524. table_formatted_BPS_SEPR <- table_model_compare %>%
  1525. kbl(caption = "Model comparision of BPS-SEPR association",
  1526. col.names = c("","number of parameters",'AIC','BIC','log likelihood','deviance','Chi-squared','df','p-value'),
  1527. row.names = F) %>%
  1528. kable_classic(full_width = F, html_font = "Cambria")
  1530. table_formatted_BPS_SEPR
  1532. # RESULT 15: Correlation: Pupil-ERP
  1533. crl <- cor(et_erp_trial[et_erp_trial$stimulus == "oddball", c(
  1534. "z_rpd_low",
  1535. "z_rpd",
  1536. "z_MMN_amplitude",
  1537. "z_P3a_amplitude"
  1538. )],
  1539. use="complete.obs")
  1540. crl
  1541. corrplot::corrplot(
  1542. crl,
  1543. method = "circle",
  1544. title = "Association between pupil data and ERPs in oddball trials",
  1545. mar=c(0,0,1,0))
  1547. crl <- cor(et_erp_trial[et_erp_trial$stimulus == "standard", c(
  1548. "z_rpd_low",
  1549. "z_rpd",
  1550. "z_MMN_amplitude",
  1551. "z_P3a_amplitude"
  1552. )],
  1553. use="complete.obs")
  1554. crl
  1555. corrplot::corrplot(
  1556. crl,
  1557. method = "circle",
  1558. title = "Association between pupil data and ERPs in standard trials",
  1559. mar=c(0,0,1,0))
  1561. # RESULT 16: Association between pupillometric indices and ERPs
  1562. ## MMN amplitude
  1563. lmm <- lmer(
  1564. z_MMN_amplitude ~ z_rpd * z_rpd_low * stimulus * manipulation * group + (1|SEGA_ID) + age + gender,
  1565. data = et_erp_trial)
  1566. anova(lmm)
  1567. r2_nakagawa(lmm)
  1569. ### post-hoc: SEPR
  1570. summary(emtrends(lmm, ~ z_rpd, var = "z_rpd"), infer = T)
  1572. ### post-hoc: BPS
  1573. summary(emtrends(lmm, ~ z_rpd_low, var = "z_rpd_low"), infer = T)
  1575. ### post-hoc: SEPR * BPS *group
  1576. summary(emtrends(lmm, ~ z_rpd * z_rpd_low|group, var = 'z_rpd', at = list(z_rpd_low = c(-2,0, 2))), infer = T)
  1578. emtr <- emtrends(lmm, ~ group | z_rpd_low, var = "z_rpd", at = list(z_rpd_low = c(-2, 0, 2)))
  1579. contrast_summary <- summary(pairs(emtr), infer = TRUE)
  1580. print(contrast_summary)
  1582. ## P3a amplitude
  1583. lmm <- lmer(
  1584. z_P3a_amplitude ~ z_rpd * z_rpd_low * stimulus * manipulation * group + (1|SEGA_ID) + age + gender,
  1585. data = et_erp_trial)
  1586. anova(lmm)
  1587. r2_nakagawa(lmm)
  1589. ### Post-hoc: SEPR
  1590. summary(emtrends(lmm, ~ z_rpd, var = "z_rpd"), infer = T)
  1592. ## BPS * group
  1593. summary(emtrends(lmm, ~ z_rpd_low|group, var = 'z_rpd_low'), infer = T)
  1595. # RESULT 17: EXPLORATORY ANALYSIS OF MODEL FIT FOR "trial_number_in_block" ON BPS
  1596. linear_fit <- lmer(
  1597. z_rpd_low ~ stimulus * manipulation * group * trial_number_in_block * block + (1|SEGA_ID) + age + gender,
  1598. data = et_erp_trial, REML = F)
  1599. anova(linear_fit)
  1601. quadratic_fit <- lmer(
  1602. z_rpd_low ~ stimulus * manipulation * group * poly(trial_number_in_block, 2) * block + (1|SEGA_ID) + age + gender,
  1603. data = et_erp_trial, REML = F)
  1604. anova(quadratic_fit)
  1605. r2_nakagawa(quadratic_fit)
  1607. cubic_fit <- lmer(
  1608. z_rpd_low ~ stimulus * manipulation * group * poly(trial_number_in_block, 3) * block + (1|SEGA_ID) + age + gender,
  1609. data = et_erp_trial, REML = F)
  1610. anova(cubic_fit)
  1611. r2_nakagawa(cubic_fit)
  1613. table_model_compare <- anova(linear_fit, quadratic_fit, cubic_fit)
  1615. table_model_compare <- cbind(
  1616. c("linear_fit", "quadratic_fit", "cubic_fit"),
  1617. table_model_compare
  1618. )
  1620. table_formatted_BPS <- table_model_compare %>%
  1621. kbl(caption = "Model comparision of trial_number_in_block on BPS",
  1622. col.names = c("","number of parameters",'AIC','BIC','log likelihood','deviance','Chi-squared','df','p-value'),
  1623. row.names = F) %>%
  1624. kable_classic(full_width = F, html_font = "Cambria")
  1626. table_formatted_BPS
  1628. # RESULT 18: EXPLORATIVE ANALYSIS OF MODEL FIT FOR "trial_number_in_block" on SEPR**
  1629. linear_fit <- lmer(
  1630. z_rpd ~ stimulus * manipulation * group * trial_number_in_block * block+ (1|SEGA_ID) + pitch + age + gender,
  1631. data = et_erp_trial, REML = F)
  1632. anova(linear_fit)
  1633. r2_nakagawa(linear_fit)
  1635. quadratic_fit <- lmer(
  1636. z_rpd ~ stimulus * manipulation * group * poly(trial_number_in_block, 2) * block + (1|SEGA_ID) + pitch + age + gender,
  1637. data = et_erp_trial, REML = F)
  1638. anova(quadratic_fit)
  1639. r2_nakagawa(quadratic_fit)
  1641. cubic_fit <- lmer(
  1642. z_rpd ~ stimulus * manipulation * group * poly(trial_number_in_block, 3) * block + (1|SEGA_ID) + pitch + age + gender,
  1643. data = et_erp_trial, REML = F)
  1644. anova(cubic_fit)
  1645. r2_nakagawa(cubic_fit)
  1647. table_model_compare <- anova(linear_fit, quadratic_fit, cubic_fit)
  1649. table_model_compare <- cbind(
  1650. c("linear_fit", "quadratic_fit", "cubic_fit"),
  1651. table_model_compare
  1652. )
  1654. table_formatted_SEPR <- table_model_compare %>%
  1655. kbl(caption = "Model comparision of trial_number_in_block on SEPR",
  1656. col.names = c("","number of parameters",'AIC','BIC','log likelihood','deviance','Chi-squared','df','p-value'),
  1657. row.names = F) %>%
  1658. kable_classic(full_width = F, html_font = "Cambria")
  1660. table_formatted_SEPR
  1662. # RESULT 19: ANALYSIS OF GRIP STRENGTH EFFECT
  1663. ## MNN amplitude
  1664. lmm <- lmer(
  1665. z_MMN_amplitude ~ stimulus * group * block * z_grip_strength + (1|SEGA_ID) + age + gender,
  1666. data = et_erp_subject[et_erp_subject$manipulation == "after", ])
  1667. anova(lmm)
  1668. r2_nakagawa(lmm)
  1670. ## MMN latency
  1671. lmm <- lmer(
  1672. z_MMN_latency ~ stimulus * group * block * z_grip_strength+ (1|SEGA_ID) + age + gender,
  1673. data = et_erp_subject[et_erp_subject$manipulation == "after", ])
  1674. anova(lmm)
  1675. r2_nakagawa(lmm)
  1677. ## P3a amplitude
  1678. lmm <- lmer(
  1679. z_P3a_amplitude ~ stimulus * group * block * z_grip_strength + (1|SEGA_ID) + age + gender,
  1680. data = et_erp_subject[et_erp_subject$manipulation == "after", ])
  1681. anova(lmm)
  1682. r2_nakagawa(lmm)
  1684. ## P3a latency
  1685. lmm <- lmer(
  1686. z_P3a_latency ~ stimulus * group * block * z_grip_strength+ (1|SEGA_ID) + age + gender,
  1687. data = et_erp_subject[et_erp_subject$manipulation == "after", ])
  1688. anova(lmm)
  1689. r2_nakagawa(lmm)
  1691. ## SEPR
  1692. lmm <- lmer(
  1693. z_rpd ~ stimulus * group * block * z_grip_strength + (1|SEGA_ID) + age + gender,
  1694. data = et_erp_subject[et_erp_subject$manipulation == "after", ])
  1695. anova(lmm)
  1696. r2_nakagawa(lmm)
  1698. ## BPS
  1699. lmm <- lmer(
  1700. z_rpd_low ~ stimulus * group * block * z_grip_strength + (1|SEGA_ID) + age + gender,
  1701. data = et_erp_subject[et_erp_subject$manipulation == "after", ])
  1702. anova(lmm)
  1703. r2_nakagawa(lmm)
  1705. # RESULT 20: Covariates Models (on subject + trial level)
  1706. ## MMN amplitude
  1707. lmm <- lmer(
  1708. z_MMN_amplitude ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender + verbal_IQ + non_verbal_IQ,
  1709. data = et_erp_subject)
  1710. anova(lmm)
  1711. r2_nakagawa(lmm)
  1713. lmm <- lmer(
  1714. z_MMN_amplitude ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender + verbal_IQ + non_verbal_IQ,
  1715. data = et_erp_trial)
  1716. anova(lmm)
  1717. r2_nakagawa(lmm)
  1719. ## MMN latency
  1720. lmm <- lmer(
  1721. z_MMN_latency ~ stimulus * manipulation * group * block + (1|SEGA_ID) + gender + age + verbal_IQ + non_verbal_IQ,
  1722. data = et_erp_subject)
  1723. anova(lmm)
  1724. r2_nakagawa(lmm)
  1726. ## Post-hoc: age
  1727. emtrends(lmm, ~ age, var = "age")
  1729. lmm <- lmer(
  1730. z_MMN_latency ~ stimulus * manipulation * group * block + (1|SEGA_ID) + gender + age + verbal_IQ + non_verbal_IQ,
  1731. data = et_erp_trial)
  1732. anova(lmm)
  1733. r2_nakagawa(lmm)
  1735. ## MMN diff amplitude
  1736. lmm <- lmer(
  1737. z_MMN_diff_amplitude ~ manipulation * group * block + (1|SEGA_ID) + gender + age + verbal_IQ + non_verbal_IQ,
  1738. data = MMN_diff)
  1739. anova(lmm)
  1740. r2_nakagawa(lmm)
  1742. ## MMN diff latency
  1743. lmm <- lmer(
  1744. z_MMN_diff_latency ~ manipulation * group * block + (1|SEGA_ID) + gender + age + verbal_IQ + non_verbal_IQ,
  1745. data = MMN_diff)
  1746. anova(lmm)
  1747. r2_nakagawa(lmm)
  1749. ## P3a amplitude
  1750. lmm <- lmer(
  1751. z_P3a_amplitude ~ stimulus * manipulation * group * block + (1|SEGA_ID) + gender + age + verbal_IQ + non_verbal_IQ,
  1752. data = et_erp_subject)
  1753. anova(lmm)
  1754. r2_nakagawa(lmm)
  1756. ## Post-hoc: age
  1757. emtrends(lmm, ~ age, var = "age")
  1759. lmm <- lmer(
  1760. z_P3a_amplitude ~ stimulus * manipulation * group * block + (1|SEGA_ID) + gender + age + verbal_IQ + non_verbal_IQ,
  1761. data = et_erp_trial)
  1762. anova(lmm)
  1763. r2_nakagawa(lmm)
  1766. ## P3a latency
  1767. lmm <- lmer(
  1768. z_P3a_latency ~ stimulus * manipulation * group * block + (1|SEGA_ID) + gender + age + verbal_IQ + non_verbal_IQ,
  1769. data = et_erp_subject)
  1770. anova(lmm)
  1771. r2_nakagawa(lmm)
  1773. lmm <- lmer(
  1774. z_P3a_latency ~ stimulus * manipulation * group * block + (1|SEGA_ID) + gender + age + verbal_IQ + non_verbal_IQ,
  1775. data = et_erp_trial)
  1776. anova(lmm)
  1777. r2_nakagawa(lmm)
  1779. ## SEPR
  1780. lmm <- lmer(
  1781. z_rpd ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender + verbal_IQ + non_verbal_IQ,
  1782. data = et_erp_subject)
  1783. anova(lmm)
  1784. r2_nakagawa(lmm)
  1786. lmm <- lmer(
  1787. z_rpd ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender + verbal_IQ + non_verbal_IQ,
  1788. data = et_erp_trial)
  1789. anova(lmm)
  1790. r2_nakagawa(lmm)
  1792. ## BPS
  1793. lmm <- lmer(
  1794. z_rpd_low ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender + verbal_IQ + non_verbal_IQ,
  1795. data = et_erp_subject)
  1796. anova(lmm)
  1797. r2_nakagawa(lmm)
  1799. ## Post-hoc: age
  1800. emtrends(lmm, ~ age, var = "age")
  1802. lmm <- lmer(
  1803. z_rpd ~ stimulus * manipulation * group * block + (1|SEGA_ID) + age + gender + verbal_IQ + non_verbal_IQ,
  1804. data = et_erp_trial)
  1805. anova(lmm)
  1806. r2_nakagawa(lmm)
  1808. # ANALYSIS OF df (ONLY PUPIL DATA-ALL DATA POINTS) ####
  1809. # Baseline phase
  1810. df_baseline <- df[df$phase %in% c("baseline", "baseline_calibration"), ]
  1811. df_baseline <- df_baseline[is.finite(df_baseline$pd), ]
  1813. # Number of pupil data per baseline trial (white, black and grey slide)
  1814. with(df_baseline[df_baseline$phase == "baseline_calibration", ], table(stimulus))
  1815. # Check: Block counter should be in line with aforementioned table
  1816. with(df_baseline, table(stimulus, block_counter))
  1817. with(df_baseline, by(
  1818. timestamp_exp, interaction(baseline_trial_counter, phase, id),
  1819. mean, na.rm = TRUE))
  1821. # mean_pd (subject baseline) is defined as mean of baseline trial during calibration phase
  1822. df_meanpd <- with(
  1823. df_baseline[df_baseline$phase == "baseline_calibration" &
  1824. df_baseline$stimulus == "baseline" , ], by(as.numeric(pd), id, mean))
  1825. df_meanpd <- data.frame(names(df_meanpd), as.numeric(df_meanpd))
  1826. names(df_meanpd) <- c("id", "mean_pd")
  1827. df_baseline <- merge(df_baseline, df_meanpd, by = "id")
  1829. # Visualization: LAPR group comparison
  1830. # New df with baseline trial only from calibration phase
  1831. df_baseline_calibration <- df_baseline[
  1832. df_baseline$phase == "baseline_calibration", ]
  1834. # Plot 1: LAPR_Subplots with trial_baseline corrected pd
  1835. ggplot(
  1836. # Baseline trial duration is theoretically 5 seconds.
  1837. df_baseline_calibration[df_baseline_calibration$ts_trial < 5.0, ],
  1838. aes(x = ts_trial,
  1839. y = pd,
  1840. group = group,
  1841. color = group)) +
  1842. geom_smooth() +
  1843. theme_bw() +
  1844. xlab("trial duration [s]") +
  1845. ylab("pupil dilation [mm]") +
  1846. labs(title = "LAPR for group and trial") +
  1847. facet_wrap(~ factor(stimulus,
  1848. levels = c("baseline", "baseline_whiteslide", "baseline_blackslide")))
  1850. ggsave(
  1851. "output/baseline_calibration_pd.tiff",
  1852. device = "tiff",
  1853. width = 6,
  1854. height = 4,
  1855. units = "in",
  1856. compression = "lzw",
  1857. dpi = 800
  1858. )
  1860. # Oddball phase
  1861. df_oddball <- df[df$phase %in% c("oddball_block", "oddball_block_rev"), ]
  1862. table(df_oddball$stimulus)
  1864. # merge with baseline PD
  1865. df_oddball <- merge(df_oddball, df_meanpd, by = "id")
  1867. # subject-baseline corrected pd
  1868. df_oddball$subj_corr_rpd <- df_oddball$pd - df_oddball$mean_pd
  1869. hist(df_oddball$subj_corr_rpd, 50)
  1871. # correct for non-finite values
  1872. df_oddball <- df_oddball[is.finite(df_oddball$subj_corr_rpd), ]
  1873. df_oddball <- df_oddball[is.finite(df_oddball$trial_corr_rpd), ]
  1875. # trial type (standard versus oddball)
  1876. df_oddball$trial_type <- substr(df_oddball$stimulus, 1, 7)
  1878. table(df_oddball$block_counter, df_oddball$stimulus)
  1879. df_oddball$manipulation <- factor(ifelse(df_oddball$block_counter < 8, "before",
  1880. ifelse(df_oddball$block_counter > 8, "after", "manipulation")),
  1881. levels = c("before", "after"))
  1883. df_oddball$order <- ifelse(grepl("rev", df_oddball$stimulus), "reverse", "normal")
  1884. table(df_oddball$order, df_oddball$block_counter)
  1886. # define trials in oddball phase
  1887. df_oddball$trial_index_oddballphase <- with(
  1888. # name to indentify individual trials
  1889. df_oddball, interaction(block_counter, .thisRepN, .thisTrialN))
  1890. # all trials including empty trials
  1891. df_oddball$trial_number_oddballphase <- with(
  1892. df_oddball, rep(seq_along(table(trial_index_oddballphase)),
  1893. times = table(trial_index_oddballphase)))
  1894. df_oddball$trial_number_oddballphase <- with(
  1895. df_oddball, rep(seq_along(table(trial_number_oddballphase)),
  1896. # rerun remove empty trials from trial number
  1897. times = table(trial_number_oddballphase)))
  1899. # Number of pupil data per trial number
  1900. table(df_oddball$trial_number)
  1901. hist(df_oddball$trial_number_oddballphase)
  1903. # Plot: Pupillary response across experiment
  1904. ggplot(
  1905. df_oddball,
  1906. aes(x = trial_number_oddballphase, y = rpd)) +
  1907. geom_smooth(method = "lm") + facet_wrap(~stimulus + manipulation)
  1908. dev.off()
  1910. # Plot: Density of manipulation per condition
  1911. ggplot(df_oddball,
  1912. aes(rpd, fill = interaction(manipulation))) +
  1913. geom_density(alpha = 0.2) + facet_wrap(~stimulus)
  1914. dev.off()
  1916. # Plot: Pupil response over a trial (only forward)
  1917. ggplot(
  1918. df_oddball[df_oddball$ts_trial < 1.8 & df_oddball$order == "normal", ],
  1919. aes(x = ts_trial,
  1920. y = scale(rpd),
  1921. group = interaction(group, trial_type),
  1922. color = group,
  1923. linetype = trial_type)) +
  1924. geom_smooth() +
  1925. theme_bw() +
  1926. labs(x = "trial duration (s)",
  1927. y = "standardized pupil response (z)",
  1928. title = "effect of manipulation of pupil response")
  1929. dev.off()
  1931. # Plot: Pupil response over a trial
  1932. ggplot(
  1933. df_oddball[df_oddball$ts_trial < 1.8, ],
  1934. aes(x = ts_trial,
  1935. y = scale(rpd),
  1936. group = interaction(order, manipulation),
  1937. color = order,
  1938. linetype = trial_type)) +
  1939. geom_smooth() +
  1940. theme_bw() +
  1941. labs(x = "trial duration (s)",
  1942. y = "standardized pupil response (z)",
  1943. title = "effect of manipulation of pupil response by block")
  1944. dev.off()
  1947. # Plot: Pupil response over a trial across groups (SEGA Paper Fig 2)
  1948. plot_dgkjp <- ggplot(
  1949. df_oddball[df_oddball$ts_trial < 1.9, ],
  1950. aes(x = ts_trial,
  1951. y = pd,
  1952. group = trial_type,
  1953. color = trial_type)) +
  1954. geom_smooth(formula = y ~ x + poly(x, 5)) +
  1955. theme_bw() +
  1956. theme(axis.text.x = element_text(face = "bold", size = 14),
  1957. axis.text.y = element_text(face = "bold", size = 14),
  1958. axis.title.x = element_text(face = "bold", size = 16),
  1959. axis.title.y = element_text(face = "bold", size = 16)) +
  1960. labs(x = "trial duration [s]",
  1961. y = "standardized pupil response [z]",
  1962. title = "Oddball Effect on pupillary response") +
  1963. scale_color_manual(values = c("red", "black"))
  1964. plot_dgkjp
  1965. dev.off()
  1968. # get legend only and save it
  1969. legend_plot_dgkjp <- get_legend(plot_dgkjp)
  1970. as_ggplot(legend_plot_dgkjp)
  1971. ggsave(
  1972. "output/plot_legend.tiff",
  1973. device = "tiff",
  1974. width = 3,
  1975. height = 2,
  1976. units = "in",
  1977. compression = "lzw",
  1978. dpi = 800
  1979. )
  1981. # Plot: Effect of manipulation on pupillary response
  1982. ggplot(
  1983. df_oddball[df_oddball$ts_trial < 1.8 & df_oddball$order == "normal", ],
  1984. aes(x = ts_trial,
  1985. y = scale(rpd),
  1986. group = interaction(group, trial_type),
  1987. color = group,
  1988. linetype = trial_type)) +
  1989. geom_smooth() +
  1990. theme_bw() +
  1991. labs(x = "trial duration (s)",
  1992. y = "standardized pupil response (z)",
  1993. title = "effect of manipulation of pupil response") +
  1994. facet_wrap(~manipulation)
  1995. dev.off()
  1997. lmm_interaction <- lmerTest::lmer(
  1998. rpd ~ trial_type * manipulation * group + (1|id),
  1999. data = df_oddball)
  2000. anova(lmm_interaction)
  2001. confint(contrast(emmeans::emmeans(
  2002. lmm_interaction, ~ group + trial_type + manipulation), method = 'pairwise'))
  2003. confint(contrast(emmeans::emmeans(
  2004. lmm_interaction, ~ manipulation + trial_type | group), method = 'pairwise'))
  2006. # Pupil response to manipulation and conditions within block
  2007. ggplot(
  2008. df_oddball[df_oddball$ts_trial < 2 & is.finite(df_oddball$rpd), ],
  2009. aes(x = trial_number_in_block,
  2010. y = scale(rpd),
  2011. group = interaction(manipulation, order, trial_type),
  2012. color = interaction(manipulation, order),
  2013. linetype = trial_type)) +
  2014. geom_smooth() + theme_bw() +
  2015. xlab("trial within block") +
  2016. ylab("relative pupil size (z)")
  2017. dev.off()
  2019. # DATA ANALYSIS: gaze data - raw pupil data ####
  2020. hist(df_oddball$left_pupil_measure1)
  2021. hist(df_oddball$right_pupil_measure1)
  2022. # raw gaze data
  2023. hist(
  2024. df_oddball$left_gaze_x[df_oddball$left_gaze_x < 200 &
  2025. df_oddball$left_gaze_x > -200])
  2026. hist(
  2027. df_oddball$right_gaze_x[df_oddball$right_gaze_x < 200 &
  2028. df_oddball$right_gaze_x > -200])
  2029. hist(
  2030. df_oddball$left_gaze_y[df_oddball$left_gaze_y < 200 &
  2031. df_oddball$left_gaze_y > -200])
  2032. hist(
  2033. df_oddball$right_gaze_y[df_oddball$right_gaze_y < 200 &
  2034. df_oddball$right_gaze_y > -200])
  2036. # refined heatmap
  2037. ggplot(
  2038. df_oddball[df_oddball$ts_trial < 2 & is.finite(df_oddball$rpd) &
  2039. is.finite(df_oddball$left_gaze_x) & is.finite(df_oddball$left_gaze_y), ],
  2040. aes(x = left_gaze_x, y = left_gaze_y)) +
  2041. stat_density_2d(aes(fill = after_stat(density)),
  2042. n = 50,
  2043. geom = "raster",
  2044. contour = FALSE) +
  2045. xlim(-200, 200) +
  2046. ylim(-200, 200) +
  2047. scale_fill_gradientn(colours = rev(rainbow(3)))
  2049. # Manipulation phase
  2050. df_manip <- df[df$phase %in% c("manipulation_block"), ]
  2052. table(df_manip$stimulus)
  2053. table(df_manip$manipulation_trial_counter)
  2054. with(df_manip, table(stimulus, manipulation_trial_counter))
  2056. # merge with baseline PD
  2057. df_manip <- merge(df_manip, df_meanpd, by = "id")
  2059. # baselinecorrected pd
  2060. df_manip$rpd <- df_manip$pd - df_manip$mean_pd
  2061. hist(df_manip$rpd)
  2063. # correct for non-finite values
  2064. df_manip <- df_manip[is.finite(df_manip$rpd), ]
  2066. # define different baseline scenarios
  2067. df_manip$trial_type <- with(
  2068. df_manip, ifelse(
  2069. manipulation_trial_counter %in% c(3, 7, 11, 15, 19),
  2070. "baseline_after_squeeze", ifelse(
  2071. manipulation_trial_counter %in% c(5, 9, 13, 17), "baseline_after_relax",
  2072. ifelse(manipulation_trial_counter %in% c(2, 6, 10, 14, 18), "squeeze",
  2073. ifelse(manipulation_trial_counter %in% c(4, 8, 12, 16, 20),
  2074. "relax", "baseline_start")))))
  2076. # 18 is duration of squeeze phase
  2077. ggplot(
  2078. df_manip[df_manip$ts_trial < 18 & df_manip$stimulus != "baseline", ],
  2079. aes(x = ts_trial,
  2080. y = rpd,
  2081. group = stimulus,
  2082. color = stimulus)) +
  2083. geom_smooth() +
  2084. theme_bw() +
  2085. labs(x = "trial duration (s)",
  2086. y = "standardized pupillary response (z)",
  2087. title = "pupil response between manipulation phases")
  2088. ggplot(
  2089. df_manip[df_manip$trial_type != "baseline_start", ],
  2090. aes(x = trial_type,
  2091. y = rpd,
  2092. fill = trial_type)) +
  2093. geom_boxplot() +
  2094. theme_bw() +
  2095. labs(x = "manipulation phase (s)",
  2096. y = "standardized pupillary response (z)",
  2097. title = "pupil response between manipulation phases")
  2099. ggplot(
  2100. df_manip,
  2101. aes(x = manipulation_trial_counter,
  2102. y = rpd,
  2103. group = manipulation_trial_counter,
  2104. fill = trial_type)) +
  2105. geom_boxplot() +
  2106. theme_bw()
  2108. lmm <- lmer(
  2109. rpd ~ trial_type + (1 | manipulation_trial_counter),
  2110. data = df_manip[df_manip$trial_type != "baseline_start", ])
  2111. anova(lmm)
  2112. plot(contrast(emmeans(lmm, ~ trial_type), "pairwise"))
  2114. # Grip strength ####
  2115. hist(df_trial$mean_grip_strength)
  2116. df_trial$mean_grip_strength_z<-scale(df_trial$mean_grip_strength)
  2117. #pupillary response - SEPR
  2118. lmm<-lmer(scale(rpd)~mean_grip_strength_z*manipulation*stimulus*group*block+(1|id),data=df_trial[df_trial$phase %in% c('oddball_block','oddball_block_rev'),])
  2119. anova(lmm) #--> no effect of hand grip strength
  2121. #pupillary response - BPS
  2122. lmm<-lmer(scale(rpd_low)~mean_grip_strength_z*manipulation*stimulus*group*block+(1|id),data=df_trial[df_trial$phase %in% c('oddball_block','oddball_block_rev'),])
  2123. anova(lmm)
  2124. fixef(lmm) #higher hand grip is assoicated with lower bps
  2125. emtrends(lmm,~manipulation|group+block,var = 'mean_grip_strength_z') ##--> higher pupil size after manipulation
  2126. #the hand grip strength has an effect on rpd only in the TD group
  2128. # Visualization
  2129. g1 <- ggplot(
  2130. df_baseline[df_baseline$phase == "baseline_calibration" &
  2131. df_baseline$ts_trial < 6, ],
  2132. aes(x = ts_trial,
  2133. y = pd,
  2134. group = stimulus,
  2135. color = stimulus)) +
  2136. geom_smooth() +
  2137. theme_bw() +
  2138. xlab("time (s)") +
  2139. ylab("pupil size (mm)") +
  2140. labs(title = "change of PD during initial baseline between IDs")
  2141. g2 <- ggplot(
  2142. df_manip[df_manip$ts_trial < 18 & df_manip$stimulus != "baseline", ],
  2143. aes(x = ts_trial,
  2144. y = rpd,
  2145. group = stimulus,
  2146. color = stimulus)) +
  2147. geom_smooth() +
  2148. theme_bw() +
  2149. labs(x = "trial duration (s)",
  2150. y = "standardized pupillary response (z)",
  2151. title = "pupil response during manipulation")
  2152. g3 <- ggplot(
  2153. df_manip[df_manip$trial_type != "baseline_start", ],
  2154. aes(x = trial_type,
  2155. y = rpd,
  2156. fill = trial_type)) +
  2157. geom_boxplot() +
  2158. theme_bw() +
  2159. labs(x = "manipulation phase (s)",
  2160. y = "standardized pupillary response (z)",
  2161. title = "pupil response in manipulation phases")
  2162. g4 <- ggplot(
  2163. df_oddball[df_oddball$ts_trial < 1.8 & df_oddball$order == "normal", ],
  2164. aes(x = ts_trial,
  2165. y = scale(rpd),
  2166. group = trial_type,
  2167. color = trial_type)) +
  2168. geom_smooth() +
  2169. theme_bw() +
  2170. labs(x = "trial duration (s)",
  2171. y = "standardized pupil response (z)",
  2172. title = "effect of manipulation on stimulus-evoked pupillary response
  2173. (before vs. after manipulation)") +
  2174. facet_wrap(~manipulation) +
  2175. scale_color_brewer(palette = "Dark2")
  2176. tiff(
  2177. file = paste0(home_path, project_path,
  2178. "/output/figure_audio_results_testdata.tiff"),
  2179. width = 12,
  2180. height = 16,
  2181. units = "in",
  2182. res = 300,
  2183. compression = "lzw")
  2184. grid.arrange(g1, g2, g3, g4,
  2185. layout_matrix = rbind(c(1, 1), c(2, 3), c(4, 4)))
  2186. dev.off()
  2188. print("The End")

data_analysis_auditory_oddball.R at commit 43374b6, no license · at the source

Overview

Authors: Anna K Müller1, Christina Luckhardt1, Christine M Freitag1, Nico Bast1
ORCID iDs: Anna K Müller, Christina Luckhardt, Christine M Freitag, Nico Bast
  1. Department of Child and Adolescent Psychiatry, Psychosomatics and Psychotherapy, University Medicine Frankfurt, Goethe-University, Deutschordenstraße 50, 60528 Frankfurt, Germany
Institutions: Goethe University Frankfurt (Germany)
Journal: Translational psychiatry, volume 16, issue 1, article 200
Dates: received 7 July 2025; accepted 6 March 2026; published online 24 March 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41398-026-03948-0 · PMID 41872152 · PMCID PMC13039879 · OpenAlex W7140176029
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: EEG (modality), human (organism), autism (population)
Methods: Spectral & time-frequency, Smoothing, state filtering, decompositions, Preprocessing, Evoked potentials, Statistics, Single-unit activity, calcium imaging, Physiology & signal measures
Keywords: Neuroscience, Diagnostic markers, Autism spectrum disorders, Psychology
MeSH: Auditory Perception*, Autism Spectrum Disorder*, Autistic Disorder*, Locus Coeruleus*, Adolescent, Attention, Case-Control Studies, Child, Electroencephalography, Evoked Potentials, Female, Humans, Male, Norepinephrine, Pupil (* major topic)
Topic: Autism Spectrum Disorder Research (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Deutsche Forschungsgemeinschaft (492582254); Daimler und Benz Stiftung (32-01/22)
Citations: not cited yet (Europe PMC); 71 references in the paper

Abstract

Sensory symptoms are common in autism and may result from differences in sensory processing. The locus coeruleus-norepinephrine (LC-NE) system modulates sensory processing by increasing selectivity to salience. To specify this mechanism, we manipulated LC-NE activity and explored sensory selectivity in autistic and non-autistic adolescents. We assessed pupillometry and electroencephalography (EEG) in 52 autistic adolescents (ASD), 55 healthy controls (CON), and 43 adolescents with other mental health conditions (MHC) during a passive auditory oddball task. A handgrip exercise manipulated LC-NE activity. Baseline pupil size (BPS) and stimulus-evoked pupillary response (SEPR) measured LC-NE tonic and phasic activity, respectively. Sensory selectivity was estimated with mismatch negativity (MMN) as change detection and P3a as bottom-up attention. Oddballs versus standards elicited increased SEPR, MMN amplitude, and longer MMN and P3a latencies. Increased P3a amplitude was more consistent in ASD. Across oddballs and standards, the manipulation transiently increased MMN and P3a amplitudes. A manipulation-induced BPS increase was specific to ASD and MHC. Pupillometric measures (BPS, SEPR) modulated sensory selectivity measures (MMN, P3a) differently between clinical groups (ASD, MHC) and CON. Findings indicate an altered orienting response to sensory stimuli in autism. The manipulation did not increase sensory selectivity but temporarily enhanced sensory reactivity to all stimuli with an arousal upregulation in both clinical groups. LC-NE activity was differentially related to sensory selectivity in autism and other mental health conditions. The tonic upregulation is discussed as increased stress susceptibility across autism and other mental health conditions.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 11 matches between paragraphs and lines of code.

sega-studie/code

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 43374b6482e1ff5c52be6dc02f807be66200e4aa, 2 November 2025
Languages: R (10)
Size: 11 files, 10 scripts
Software Heritage: not archived
Found in: the text, “Pupillometry”
Holds: README, 4 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (8 files), emmeans (6 files), lme4 (6 files), lmerTest (6 files), easystats (5 files), ggplot2 (5 files), ggpubr (5 files), cowplot (4 files), data.table (4 files), reshape2 (2 files), psych (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
11 files
At the source: github.com/SEGA-Studie/code/blob/main/SEGA_Preprocessing.Rmd

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 10 scripts, each with its path and the digest of its content;
  • 11 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

Scripts for pupil data preprocessing and full analysis can be found at github, accessible via the link provided in the methods. All data supporting this study are described in the text and in the supplementary material. Due to ethical considerations and data privacy, the raw data that is tied to person identifiable information cannot be publicly shared.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 4 keywords, 15 MeSH terms, 2 funders, 62 references.

Cite

This paper

Müller, A. K., Luckhardt, C., Freitag, C. M., & Bast, N. (2026). Locus coeruleus modulation of neurophysiological sensory selectivity differs in autism and other mental health conditions. Translational psychiatry, 16(1), 200. https://doi.org/10.1038/s41398-026-03948-0

BibTeX

@article{muller2026locus,
author = {Müller, Anna K and Luckhardt, Christina and Freitag, Christine M and Bast, Nico},
title = {{Locus coeruleus modulation of neurophysiological sensory selectivity differs in autism and other mental health conditions}},
journal = {Translational psychiatry},
year = {2026},
month = mar,
volume = {16},
number = {1},
pages = {200},
publisher = {Nature Publishing Group},
issn = {2158-3188},
doi = {10.1038/s41398-026-03948-0},
url = {https://doi.org/10.1038/s41398-026-03948-0},
pmid = {41872152},
pmcid = {PMC13039879}
}

RIS

TY - JOUR
AU - Müller, Anna K
AU - Luckhardt, Christina
AU - Freitag, Christine M
AU - Bast, Nico
TI - Locus coeruleus modulation of neurophysiological sensory selectivity differs in autism and other mental health conditions
T2 - Translational psychiatry
J2 - Transl Psychiatry
PY - 2026
DA - 2026/03/24
VL - 16
IS - 1
SP - 200
SN - 2158-3188
PB - Nature Publishing Group
DO - 10.1038/s41398-026-03948-0
UR - https://doi.org/10.1038/s41398-026-03948-0
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41398-026-03948-0",
"type": "article-journal",
"title": "Locus coeruleus modulation of neurophysiological sensory selectivity differs in autism and other mental health conditions",
"container-title": "Translational psychiatry",
"author": [
{
"family": "Müller",
"given": "Anna K"
},
{
"family": "Luckhardt",
"given": "Christina"
},
{
"family": "Freitag",
"given": "Christine M"
},
{
"family": "Bast",
"given": "Nico"
}
],
"container-title-short": "Transl Psychiatry",
"volume": "16",
"issue": "1",
"page": "200",
"DOI": "10.1038/s41398-026-03948-0",
"PMID": "41872152",
"PMCID": "PMC13039879",
"ISSN": "2158-3188",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41398-026-03948-0",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
24
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1186/s13229-025-00678-w [code]
Locus coeruleus tonic upregulation increases selectivity to inconspicuous auditory information in autistic compared to non-autistic individuals: a combined pupillometry and electroencephalography study.
Journal: Molecular autism
In common: easystats, psych, emmeans, 6 other tools, autism, EEG, 18 references, 4 authors
[2] doi:10.1186/s13229-023-00537-6 [code]
Sensory salience processing moderates attenuated gazes on faces in autism spectrum disorder: a case-control study.
Journal: Molecular autism
In common: psych, emmeans, lmerTest, 4 other tools, autism, 13 references
[3] doi:10.1371/journal.pone.0355165 [code]
Pupillary dynamics during hands-off L2 driving and transitions of control under high cognitive load.
Journal: PloS one
In common: easystats, emmeans, lmerTest, 5 other tools, 7 references
[4] doi:10.1038/s41598-023-28781-z [code]
Changes in pupil dilation and P300 amplitude indicate the possible involvement of the locus coeruleus-norepinephrine (LC-NE) system in psychological flow.
Journal: Scientific reports
In common: easystats, lmerTest, lme4, 4 other tools, EEG, 7 references
[5] doi:10.7554/elife.109678 [code]
Strong mnemonic prediction errors increase cognitive control, attention, and arousal.
Journal: eLife
In common: easystats, psych, emmeans, 5 other tools, EEG, 4 references
[6] doi:10.3758/s13423-022-02124-x [code]
Explicitly predicting outcomes enhances learning of expectancy-violating information.
Journal: Psychonomic bulletin & review
In common: easystats, psych, emmeans, 6 other tools, 3 references
[7] doi:10.3758/s13415-023-01096-2 [code]
Testing locus coeruleus-norepinephrine accounts of working memory, attention control, and fluid intelligence.
Journal: Cognitive, affective & behavioral neuroscience
In common: psych, lmerTest, cowplot, 3 other tools, 6 references
[8] doi:10.1093/scan/nsag009 [code]
Covert neural and autonomic signatures of shared perception.
Journal: Social cognitive and affective neuroscience
In common: easystats, emmeans, lmerTest, 5 other tools, EEG, 2 references
[9] doi:10.7554/elife.84822 [code]
Pupil size reflects activation of subcortical ascending arousal system nuclei during rest.
Journal: eLife
In common: psych, lme4, ggpubr, 3 other tools, 4 references
[10] doi:10.1038/s41598-023-41357-1
Acceleration of inferred neural responses to oddball targets in an individual with bilateral amygdala lesion compared to healthy controls.
Journal: Scientific reports
In common: 8 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1038/s41398-026-03948-0.svg)](https://openscicode.org/badge/doi/10.1038/s41398-026-03948-0)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1038/s41398-026-03948-0.svg
   :target: https://openscicode.org/badge/doi/10.1038/s41398-026-03948-0
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1038/s41398-026-03948-0"><img src="https://openscicode.org/badge/doi/10.1038/s41398-026-03948-0.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.