OSCR

Changes in Aperiodic (1/f Slope) Activity During a Picture-Word Interference Task: Effects of Congruency and Sequence Manipulations.

Code ↔ Paper

2 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 2 matches
  1. [1] § Method › Cluster‐Based Permutation Analyses (FWER‐Corrected) ↔ code_afterR1.Rmd, lines 1193–1254 · score 0.57 · sign flip, cluster statistic, FWER, permutation
  2. [2] § Method › Exploratory Analyses (No Correction for Multiple Comparisons) ↔ code_afterR1.Rmd, lines 1193–1254 · score 0.56 · global stimulus induced, Baseline correction, threshold, post, interactions

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R Markdown · 1,689 lines · 61 KB · no license · 2 matches

  1. ---
  2. title: "Code to reproduce the spectral slope analyses from the project (section 3.2.3)."
  3. output: html_document
  4. date: "2024-04-04"
  5. ---
  7. Load libraries
  9. ```{r setup, include=FALSE}
  10. knitr::opts_chunk$set(echo = TRUE)
  12. setwd("C:/Users/patry/Documents/Paranauka/Projekty/Illinois/virginia/r1/")
  13. rm(list = ls())
  14. options(scipen = 999)
  16. library(plyr)
  17. library(tidyverse)
  18. library(psych)
  19. library(stringr)
  20. library(data.table)
  21. library(readxl)
  22. library(eegUtils)
  24. # Time intervals
  25. time_intervals <- c(-160, 0, 160, 320, 480, 640, 800, 960, 1120, 1280)
  26. # Define the correct order of the levels
  27. time_chunk_levels <- c("[-160,0]", "(0,160]", "(160,320]", "(320,480]", "(480,640]", "(640,800]", "(800,960]", "(960,1120]", "(1120,1280]")
  28. # Shading for plots
  29. shading_data <- data.frame(
  30. xmin = time_intervals[-length(time_intervals)], # All but the last element
  31. xmax = time_intervals[-1], # All but the first element
  32. ymin = -Inf, # Extend shading to the bottom
  33. ymax = Inf, # Extend shading to the top
  34. brightness = c(0.0, 0.6,0.1, 0.6,0.1, 0.6,0.1, 0.6,0.1) # Different brightness values
  35. )
  37. ```
  40. Raw data vizualization
  42. ```{r}
  44. load("data_raw.RData")
  46. # Collapse across non-focal dimensions (mean within cells)
  47. reg_pcongr2 = data.table::dcast(setDT(reg), ID + channel + time + PrevCongr + Congr ~ ., value.var = c("value"), mean, na.rm = TRUE)
  48. setnames(reg_pcongr2, old = ".", new = "value")
  50. # Harmonize factors & keep only levels of interest
  51. reg_pcongr2$Congr <- factor(reg_pcongr2$Congr, levels = c("Con", "Inc"))
  52. reg_pcongr2$PrevCongr <- factor(reg_pcongr2$PrevCongr, levels = c("pCon", "pInc"))
  54. reg_pcongr3 = reg_pcongr2[reg_pcongr2$channel %in% c("Fz", "Cz", "Pz", "Oz")]
  55. reg_pcongr3$channel <- factor(reg_pcongr3$channel, levels = c("Fz", "Cz", "Pz", "Oz"))
  57. # Define ROIs
  58. reg_pcongr2 <- reg_pcongr2 %>%
  59. mutate(
  60. ROI = case_when(
  61. channel %in% c("F1","F2","Fz") ~ "Frontal",
  62. channel %in% c("FC1","FC2","FCz") ~ "Fronto-Central",
  63. channel %in% c("C1","C2","Cz") ~ "Central",
  64. channel %in% c("CP1","CP2","CPz") ~ "Centro-Parietal",
  65. channel %in% c("P1","P2","Pz") ~ "Parietal",
  66. channel %in% c("O1","O2","Oz") ~ "Occipital",
  67. TRUE ~ NA_character_ # for any channels not listed
  68. )
  69. )
  70. reg_pcongr4 <- reg_pcongr2 %>%
  71. group_by(ID, time, PrevCongr, Congr, ROI) %>%
  72. summarise(value = mean(value, na.rm = TRUE)) %>%
  73. ungroup()
  75. reg_pcongr4$ROI <- factor(reg_pcongr4$ROI, levels = c("Frontal", "Fronto-Central", "Central", "Centro-Parietal", "Parietal", "Occipital"))
  76. reg_pcongr4 = na.omit(reg_pcongr4)
  78. ### Plot (midline electrodes)
  79. ggplot(reg_pcongr3, aes(x = time, y = value, color = Congr, linetype = PrevCongr)) +
  80. facet_wrap(channel ~ ., ncol = 1) +
  82. # Mean lines only
  83. stat_summary(fun = "mean", geom = "line", size = .2) +
  85. # Manual color: Con = blue, Inc = red
  86. scale_color_manual(values = c("Con" = "blue", "Inc" = "red")) +
  88. # Theme and layout
  89. theme_bw(base_size = 14) +
  90. theme(
  91. panel.grid = element_blank(),
  92. strip.background = element_rect(fill = "white"),
  93. strip.text = element_text(size = 12, face = "bold"),
  94. axis.title = element_text(size = 12),
  95. axis.text = element_text(size = 10),
  96. plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
  97. legend.position = "top"
  98. ) +
  100. # Reference lines
  101. #geom_hline(yintercept = 0, linetype = "dashed") +
  102. geom_vline(xintercept = 0, linetype = "dashed") +
  103. #geom_vline(xintercept = 500, linetype = "dashed", color = "blue") +
  105. # Axes and labels
  106. coord_cartesian(ylim = c(-1.5,-.5)) +
  107. labs(
  108. x = "Time (ms)",
  109. y = expression("Spectral Slope Values [a.u.]"),
  110. color = "Congruency", # Current trial
  111. linetype = "Previous Congr." # Previous trial
  112. ) +
  113. ggtitle("")
  115. # ggsave(file.path(folder_name, "raw_overview.jpg"), width = 8, height = 15, dpi = 600, scale = .5)
  116. # ggsave(file.path(folder_name, "raw_overview.pdf"), width = 8, height = 15, dpi = 600, scale = .5)
  118. ### Plot (electrode averages)
  119. ggplot(reg_pcongr4[reg_pcongr4$ROI %in% c("Frontal", "Central", "Occipital"),], aes(x = time, y = value, color = Congr, linetype = PrevCongr)) +
  120. facet_wrap(ROI ~ ., ncol = 1) +
  122. # Mean lines only
  123. stat_summary(fun = "mean", geom = "line", size = .2) +
  125. # Manual color: Con = blue, Inc = red
  126. scale_color_manual(values = c("Con" = "blue", "Inc" = "red")) +
  128. # Theme and layout
  129. theme_bw(base_size = 14) +
  130. theme(
  131. panel.grid = element_blank(),
  132. strip.background = element_rect(fill = "white"),
  133. strip.text = element_text(size = 12, face = "bold"),
  134. axis.title = element_text(size = 12),
  135. axis.text = element_text(size = 10),
  136. plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
  137. legend.position = "top"
  138. ) +
  140. # Reference lines
  141. #geom_hline(yintercept = 0, linetype = "dashed") +
  142. geom_vline(xintercept = 0, linetype = "dashed") +
  143. #geom_vline(xintercept = 500, linetype = "dashed", color = "blue") +
  145. # Axes and labels
  146. coord_cartesian(ylim = c(-1.5,-.5)) +
  147. labs(
  148. x = "Time (ms)",
  149. y = expression("Spectral Slope Values [a.u.]"),
  150. color = "Congruency", # Current trial
  151. linetype = "Previous Congr." # Previous trial
  152. ) +
  153. ggtitle("")
  155. # ggsave(file.path(folder_name, "raw_overview_av.jpg"), width = 8, height = 15, dpi = 600, scale = .5)
  156. # ggsave(file.path(folder_name, "raw_overview_av.pdf"), width = 8, height = 15, dpi = 600, scale = .5)
  159. rm(reg, reg_pcongr2, reg_pcongr3, reg_pcongr4)
  161. ```
  164. Create data sets and rename conditions
  166. ```{r}
  168. load("data_Pca.RData")
  169. reg <- scores
  171. #Rename and collapse across conditions
  172. # Create a named vector for bin labels
  173. bin_labels <- c(
  174. "Novel-RSI1000-Changed-pCon-Con",
  175. "Novel-RSI1000-Changed-pCon-Inc",
  176. "Novel-RSI1000-Changed-pInc-Con",
  177. "Novel-RSI1000-Changed-pInc-Inc",
  178. "Novel-RSI1000-Repeated-pCon-Con",
  179. "Novel-RSI1000-Repeated-pCon-Inc",
  180. "Novel-RSI1000-Repeated-pInc-Con",
  181. "Novel-RSI1000-Repeated-pInc-Inc",
  182. "Novel-RSI2000-Changed-pCon-Con",
  183. "Novel-RSI2000-Changed-pCon-Inc",
  184. "Novel-RSI2000-Changed-pInc-Con",
  185. "Novel-RSI2000-Changed-pInc-Inc",
  186. "Novel-RSI2000-Repeated-pCon-Con",
  187. "Novel-RSI2000-Repeated-pCon-Inc",
  188. "Novel-RSI2000-Repeated-pInc-Con",
  189. "Novel-RSI2000-Repeated-pInc-Inc",
  190. "Novel-RSI3000-Changed-pCon-Con",
  191. "Novel-RSI3000-Changed-pCon-Inc",
  192. "Novel-RSI3000-Changed-pInc-Con",
  193. "Novel-RSI3000-Changed-pInc-Inc",
  194. "Novel-RSI3000-Repeated-pCon-Con",
  195. "Novel-RSI3000-Repeated-pCon-Inc",
  196. "Novel-RSI3000-Repeated-pInc-Con",
  197. "Novel-RSI3000-Repeated-pInc-Inc",
  198. "Novel-RSI5000-Changed-pCon-Con",
  199. "Novel-RSI5000-Changed-pCon-Inc",
  200. "Novel-RSI5000-Changed-pInc-Con",
  201. "Novel-RSI5000-Changed-pInc-Inc",
  202. "Novel-RSI5000-Repeated-pCon-Con",
  203. "Novel-RSI5000-Repeated-pCon-Inc",
  204. "Novel-RSI5000-Repeated-pInc-Con",
  205. "Novel-RSI5000-Repeated-pInc-Inc",
  206. "Frequent-RSI1000-Changed-pCon-Con",
  207. "Frequent-RSI1000-Changed-pCon-Inc",
  208. "Frequent-RSI1000-Changed-pInc-Con",
  209. "Frequent-RSI1000-Changed-pInc-Inc",
  210. "Frequent-RSI1000-Repeated-pCon-Con",
  211. "Frequent-RSI1000-Repeated-pCon-Inc",
  212. "Frequent-RSI1000-Repeated-pInc-Con",
  213. "Frequent-RSI1000-Repeated-pInc-Inc",
  214. "Frequent-RSI2000-Changed-pCon-Con",
  215. "Frequent-RSI2000-Changed-pCon-Inc",
  216. "Frequent-RSI2000-Changed-pInc-Con",
  217. "Frequent-RSI2000-Changed-pInc-Inc",
  218. "Frequent-RSI2000-Repeated-pCon-Con",
  219. "Frequent-RSI2000-Repeated-pCon-Inc",
  220. "Frequent-RSI2000-Repeated-pInc-Con",
  221. "Frequent-RSI2000-Repeated-pInc-Inc",
  222. "Frequent-RSI3000-Changed-pCon-Con",
  223. "Frequent-RSI3000-Changed-pCon-Inc",
  224. "Frequent-RSI3000-Changed-pInc-Con",
  225. "Frequent-RSI3000-Changed-pInc-Inc",
  226. "Frequent-RSI3000-Repeated-pCon-Con",
  227. "Frequent-RSI3000-Repeated-pCon-Inc",
  228. "Frequent-RSI3000-Repeated-pInc-Con",
  229. "Frequent-RSI3000-Repeated-pInc-Inc",
  230. "Frequent-RSI5000-Changed-pCon-Con",
  231. "Frequent-RSI5000-Changed-pCon-Inc",
  232. "Frequent-RSI5000-Changed-pInc-Con",
  233. "Frequent-RSI5000-Changed-pInc-Inc",
  234. "Frequent-RSI5000-Repeated-pCon-Con",
  235. "Frequent-RSI5000-Repeated-pCon-Inc",
  236. "Frequent-RSI5000-Repeated-pInc-Con",
  237. "Frequent-RSI5000-Repeated-pInc-Inc"
  238. )
  239. names(bin_labels) <- as.character(1:64)
  240. # First, ensure the Bin column is numeric for matching
  241. reg$Bin <- as.numeric(as.character(reg$Bin))
  242. # Now, replace the numeric Bins with their corresponding labels
  243. reg$BinLabel <- bin_labels[as.character(reg$Bin)]
  244. # Separate the 'BinLabel' column into new columns
  245. reg <- separate(reg, col = BinLabel, into = c("StimulusRep", "RSI", "CategoryRep", "PrevCongr", "Congr"), sep = "-")
  247. #Reorder columns
  248. reg <- reg[,c("ID", "time", "Bin", "StimulusRep", "CategoryRep", "RSI", "PrevCongr", "Congr", "MR1", "MR2", "MR3", "MR4", "MR5")]
  250. #Narrow down the time window for analyses
  251. reg <- reg[reg$time %in% c(-160:1280),]
  253. #To the long format
  254. reg <- melt(reg, id.vars = c("ID", "time", "Bin", "StimulusRep", "CategoryRep", "RSI", "PrevCongr", "Congr"))
  256. #Factor order
  257. rotFit$variance_explained_percent
  258. unique(reg$variable)
  260. # Get the factor names ordered by variance explained
  261. ordered_factors <- names(sort(rotFit$variance_explained_percent, decreasing = TRUE))
  262. # Reorder reg$variable based on the ordered factors
  263. reg$variable <- factor(reg$variable, levels = ordered_factors)
  264. # Verify the order
  265. levels(reg$variable)
  267. # Rename levels of reg$variable
  268. levels(reg$variable) <- c("OCCIPITAL (MR3)", "FRONTAL (MR2)", "CENTRAL (MR1)", "TEMPORAL L (MR5)", "TEMPORAL R (MR4)")
  271. #Exclude factors with too small variance
  272. reg <- reg[!reg$variable %in% c("TEMPORAL L (MR5)", "TEMPORAL R (MR4)"),]
  273. reg <- droplevels(reg)
  276. # Reorder factor levels
  277. reg$variable <- factor(reg$variable, levels = c("FRONTAL (MR2)", "CENTRAL (MR1)", "OCCIPITAL (MR3)"))
  278. # Verify the new order
  279. levels(reg$variable)
  281. # Drop the last two elements
  282. variance_explained <- rotFit$variance_explained[1:3]
  283. # Rename the variables
  284. names(variance_explained) <- c("CENTRAL", "FRONTAL", "OCCIPITAL")
  285. # Reorder them: frontal, central, occipital
  286. variance_explained <- variance_explained[c("FRONTAL", "CENTRAL", "OCCIPITAL")]
  287. # Format the reordered variance explained as a text string for the caption
  288. caption_text <- paste(
  289. "Variance Explained:",
  290. paste(names(variance_explained), sprintf("%.2f%%", variance_explained), collapse = ", ")
  291. )
  293. ```
  295. Info for stats & plots
  297. ```{r}
  299. reg$value_to_stats = reg$value
  301. ```
  303. Global Changes in Spectral Slope Across Time Windows and Scalp Regions
  305. ```{r}
  307. ### Average across bins
  308. reg_av = data.table::dcast(setDT(reg), ID + time + variable ~ ., value.var = c("value_to_stats"), mean, na.rm = TRUE)
  309. setnames(reg_av, old = ".", new = "value")
  311. reg_av$time_chunk <- cut(reg_av$time, breaks = time_intervals, include.lowest = TRUE, ordered_result = TRUE, dig.lab = 5)
  312. reg_av$time_chunk <- factor(reg_av$time_chunk, levels = time_chunk_levels, ordered = FALSE)
  314. ### Statistics
  315. dat = reg_av
  316. dat = dat[complete.cases(dat$time_chunk)]
  317. dat <- setDT(dcast(dat, ID + variable + time_chunk ~ ., mean, value.var = c("value"), na.rm = TRUE))
  318. dat = data.table(droplevels(dat))
  319. dat = dat[complete.cases(dat$time_chunk)]
  321. # Create an empty data.table to store the t-test results
  322. pvals <- data.table(
  323. variable = character(), time_chunk = character(),
  324. t.value = numeric(), p.value = numeric(), cohen_d = numeric()
  325. )
  327. # Get unique levels
  328. unique_time <- unique(dat$time_chunk)
  329. unique_var <- unique(dat$variable)
  331. for (tpnt in unique_time) {
  332. for (var in unique_var) {
  333. subset_data <- dat[variable == var & time_chunk %in% c(tpnt, "[-160,0]")]
  335. # Get paired data
  336. baseline <- subset_data[time_chunk == "[-160,0]"][order(ID)]$.
  337. tp_data <- subset_data[time_chunk == tpnt][order(ID)]$.
  339. # Ensure same IDs are used for pairing
  340. common_ids <- intersect(
  341. subset_data[time_chunk == "[-160,0]"]$ID,
  342. subset_data[time_chunk == tpnt]$ID
  343. )
  344. baseline <- subset_data[time_chunk == "[-160,0]" & ID %in% common_ids][order(ID)]$.
  345. tp_data <- subset_data[time_chunk == tpnt & ID %in% common_ids][order(ID)]$.
  347. # Paired t-test
  348. t_test_result <- t.test(tp_data, baseline,
  349. alternative = "two.sided", paired = TRUE, conf.level = 0.95)
  351. # Cohen's d for paired samples
  352. diff_scores <- tp_data - baseline
  353. cohen_d <- mean(diff_scores) / sd(diff_scores)
  355. # Append results
  356. pvals <- rbind(
  357. pvals,
  358. data.table(
  359. variable = var, time_chunk = tpnt,
  360. t.value = t_test_result$statistic,
  361. p.value = t_test_result$p.value,
  362. cohen_d = cohen_d
  363. )
  364. )
  365. }
  366. }
  368. # Significance flags
  369. pvals$crit1 <- ifelse(pvals$p.value <= .05, 1, NA)
  370. pvals$crit2 <- ifelse(pvals$p.value <= .01, 1, NA)
  372. # Convert time_chunk to a factor with the specified levels
  373. pvals$time_chunk <- factor(pvals$time_chunk, levels = time_chunk_levels)
  374. pvals$variable <- factor(pvals$variable, levels = c("FRONTAL (MR2)", "CENTRAL (MR1)", "OCCIPITAL (MR3)"))
  375. writexl::write_xlsx(pvals, file.path("ttest_results.xlsx"))
  377. ### Plot with stats
  379. # Ensure shading_data$brightness is a factor
  380. shading_data$brightness <- as.factor(shading_data$brightness)
  382. ggplot(reg_av, aes(x = time, y = value)) +
  384. # Facet by variable
  385. facet_wrap(. ~ variable, ncol = 1, drop = FALSE) +
  387. # Add shaded regions
  388. geom_rect(data = shading_data,
  389. aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
  390. fill = factor(brightness)),
  391. alpha = 0.5, inherit.aes = FALSE, show.legend = FALSE) +
  392. scale_fill_manual(values = scales::alpha(
  393. c("lightgrey", "lightgrey", "grey", "lightgrey", "grey",
  394. "lightgrey", "grey", "lightgrey", "grey", "lightgrey", "grey"), 0.5)) +
  397. # Mean ± SE ribbon and mean line
  398. stat_summary(fun.data = mean_se, geom = "ribbon", alpha = .15, fill = "grey50") +
  399. stat_summary(fun = "mean", geom = "line", size = 1, color = "black") +
  401. # APA-style theme
  402. theme_bw(base_size = 12) +
  403. theme(
  404. panel.grid = element_blank(),
  405. strip.background = element_rect(fill = "white"),
  406. strip.text = element_text(size = 12, face = "bold"),
  407. axis.title = element_text(size = 12),
  408. axis.text = element_text(size = 10),
  409. plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
  410. legend.position = "none"
  411. ) +
  413. # Axes and labels
  414. geom_hline(yintercept = 0, linetype = "dashed", color = "black", size = 0.3) +
  415. geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.3) +
  416. labs(
  417. x = "Time (ms)",
  418. y = expression("Spectral Slope Values [a.u.]"),
  419. title = ""
  420. ) +
  422. # Horizontal lines for significant time points (black for crit1)
  423. geom_segment(
  424. data = reg_av %>%
  425. inner_join(pvals[pvals$crit1 == 1], by = c("time_chunk", "variable")),
  426. aes(x = time - 20, xend = time + 20,
  427. y = min(value, na.rm = TRUE) - 0.1,
  428. yend = min(value, na.rm = TRUE) - 0.1),
  429. inherit.aes = FALSE, color = "black", size = 1.5
  430. ) +
  432. # Horizontal lines for significant time points (dark green for crit2)
  433. geom_segment(
  434. data = reg_av %>%
  435. inner_join(pvals[pvals$crit2 == 1], by = c("time_chunk", "variable")),
  436. aes(x = time - 20, xend = time + 20,
  437. y = min(value, na.rm = TRUE) - 0.1,
  438. yend = min(value, na.rm = TRUE) - 0.1),
  439. inherit.aes = FALSE, color = "darkgreen", size = 1.5
  440. ) +
  442. # Manual y-axis range
  443. coord_cartesian(ylim = c(min(reg_av$value, na.rm = TRUE) - 0.3, NA))
  446. # ggsave(file.path(folder_name, "global_ttest.jpg"), width = 8, height = 15, dpi = 600, scale = .5)
  447. # ggsave(file.path(folder_name, "global_ttest.pdf"), width = 8, height = 15, dpi = 600, scale = .5)
  449. ```
  451. Repeated-Measures ANOVA Assessing the Effects of Experimental Manipulation on Spectral Slope across Time Windows and Scalp Regions
  453. ```{r}
  455. ### Average across bins
  456. reg_pcongr2 = data.table::dcast(setDT(reg), ID + variable + time + PrevCongr + Congr ~ ., value.var = c("value_to_stats"), mean, na.rm = TRUE)
  457. setnames(reg_pcongr2, old = ".", new = "value")
  459. reg_pcongr2$Congr <- factor(reg_pcongr2$Congr, levels = c("Con", "Inc"))
  460. reg_pcongr2$PrevCongr <- factor(reg_pcongr2$PrevCongr, levels = c("pCon", "pInc"))
  461. reg_pcongr2$time_chunk <- cut(reg_pcongr2$time, breaks = time_intervals, include.lowest = TRUE, ordered_result = TRUE, dig.lab = 5)
  462. reg_pcongr2$time_chunk <- factor(reg_pcongr2$time_chunk, levels = time_chunk_levels, ordered = FALSE)
  464. ### Statistics
  465. dat = reg_pcongr2
  466. dat = data.table(droplevels(dat))
  467. dat = setDT(dat) %>% dcast(ID + variable + time_chunk + PrevCongr + Congr ~ ., mean, value.var = c("value"), na.rm = TRUE)
  468. dat = dat[complete.cases(dat$time_chunk)]
  469. dat$ID = as.factor(dat$ID)
  470. colnames(dat)[6] <- "value"
  472. # Create an empty data.table to store the ANOVA results
  473. pvals <- data.table(
  474. time_chunk = numeric(), variable = character(),
  475. Congr_F = numeric(), Congr_p = numeric(), Congr_eta2p = numeric(),
  476. PrevCongr_F = numeric(), PrevCongr_p = numeric(), PrevCongr_eta2p = numeric(),
  477. PrevCongrCongr_F = numeric(), PrevCongrCongr_p = numeric(), PrevCongrCongr_eta2p = numeric()
  478. )
  480. # Get unique levels
  481. unique_variables <- unique(dat$variable)
  482. unique_time <- droplevels(unique(dat$time_chunk)[-1])
  484. for (tpnt in unique_time) {
  485. for (var in unique_variables) {
  486. subset_data <- dat %>% filter(time_chunk == tpnt, variable == var)
  487. result <- summary(stats::aov(data = subset_data, value ~ PrevCongr * Congr + Error(ID/(PrevCongr * Congr))))
  489. # Congr
  490. congr_ss_effect <- result$`Error: ID:Congr`[[1]]$`Sum Sq`[1]
  491. congr_ss_error <- result$`Error: ID:Congr`[[1]]$`Sum Sq`[2]
  492. congr_f <- result$`Error: ID:Congr`[[1]]$`F value`[1]
  493. congr_p <- result$`Error: ID:Congr`[[1]]$`Pr(>F)`[1]
  494. congr_eta2p <- congr_ss_effect / (congr_ss_effect + congr_ss_error)
  496. # PrevCongr
  497. prevcongr_ss_effect <- result$`Error: ID:PrevCongr`[[1]]$`Sum Sq`[1]
  498. prevcongr_ss_error <- result$`Error: ID:PrevCongr`[[1]]$`Sum Sq`[2]
  499. prevcongr_f <- result$`Error: ID:PrevCongr`[[1]]$`F value`[1]
  500. prevcongr_p <- result$`Error: ID:PrevCongr`[[1]]$`Pr(>F)`[1]
  501. prevcongr_eta2p <- prevcongr_ss_effect / (prevcongr_ss_effect + prevcongr_ss_error)
  503. # Interaction
  504. inter_ss_effect <- result$`Error: ID:PrevCongr:Congr`[[1]]$`Sum Sq`[1]
  505. inter_ss_error <- result$`Error: ID:PrevCongr:Congr`[[1]]$`Sum Sq`[2]
  506. prevcongrcongr_f <- result$`Error: ID:PrevCongr:Congr`[[1]]$`F value`[1]
  507. prevcongrcongr_p <- result$`Error: ID:PrevCongr:Congr`[[1]]$`Pr(>F)`[1]
  508. prevcongrcongr_eta2p <- inter_ss_effect / (inter_ss_effect + inter_ss_error)
  510. # Append to results table
  511. pvals <- rbind(
  512. pvals,
  513. data.table(
  514. time_chunk = tpnt, variable = var,
  515. Congr_F = congr_f, Congr_p = congr_p, Congr_eta2p = congr_eta2p,
  516. PrevCongr_F = prevcongr_f, PrevCongr_p = prevcongr_p, PrevCongr_eta2p = prevcongr_eta2p,
  517. PrevCongrCongr_F = prevcongrcongr_f, PrevCongrCongr_p = prevcongrcongr_p, PrevCongrCongr_eta2p = prevcongrcongr_eta2p
  518. )
  519. )
  520. }
  521. }
  523. # Add significance flags
  524. pvals$Congr_crit1 <- ifelse(pvals$Congr_p <= .05, 1, NA)
  525. pvals$Congr_crit2 <- ifelse(pvals$Congr_p <= .01, 1, NA)
  526. pvals$PrevCongr_crit1 <- ifelse(pvals$PrevCongr_p <= .05, 1, NA)
  527. pvals$PrevCongr_crit2 <- ifelse(pvals$PrevCongr_p <= .01, 1, NA)
  528. pvals$PrevCongrCongr_crit1 <- ifelse(pvals$PrevCongrCongr_p <= .05, 1, NA)
  529. pvals$PrevCongrCongr_crit2 <- ifelse(pvals$PrevCongrCongr_p <= .01, 1, NA)
  531. # Convert time_chunk to a factor with the specified levels
  532. pvals$time_chunk <- factor(pvals$time_chunk, levels = time_chunk_levels, ordered = FALSE)
  533. pvals$variable <- factor(pvals$variable, levels = c("FRONTAL (MR2)", "CENTRAL (MR1)", "OCCIPITAL (MR3)"))
  534. writexl::write_xlsx(pvals, file.path("anova_results.xlsx"))
  537. ###DIFFERENCE WAVES TO PLOT
  539. ###Congr
  541. # Calculate difference between Incongruent and Congruent
  542. aov_congr <- reg_pcongr2 %>%
  543. pivot_wider(names_from = Congr, values_from = value) %>%
  544. mutate(value = Inc - Con) %>%
  545. select(time, time_chunk, variable, value)
  547. ggplot(aov_congr, aes(x = time, y = value)) +
  548. facet_wrap(. ~ variable, ncol = 1, drop = FALSE) +
  550. # Add shaded regions
  551. geom_rect(data = shading_data,
  552. aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
  553. fill = factor(brightness)),
  554. alpha = 0.5, inherit.aes = FALSE, show.legend = FALSE) +
  555. scale_fill_manual(values = scales::alpha(
  556. c("lightgrey", "lightgrey", "grey", "lightgrey", "grey",
  557. "lightgrey", "grey", "lightgrey", "grey", "lightgrey", "grey"), 0.5)) +
  559. # Mean ± SE ribbon and mean line
  560. stat_summary(fun.data = mean_se, geom = "ribbon", alpha = .1, fill = "grey50") +
  561. stat_summary(fun = "mean", geom = "line", size = 1.2, color = "black") +
  563. # Theme and layout
  564. theme_bw(base_size = 14) +
  565. theme(
  566. panel.grid = element_blank(),
  567. strip.background = element_rect(fill = "white"),
  568. strip.text = element_text(size = 12, face = "bold"),
  569. axis.title = element_text(size = 12),
  570. axis.text = element_text(size = 10),
  571. plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
  572. legend.position = "none"
  573. ) +
  575. # Axes and labels
  576. geom_hline(yintercept = 0, linetype = "dashed", color = "black", size = 0.3) +
  577. geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.3) +
  578. coord_cartesian(ylim = c(-.12, .1)) +
  579. labs(
  580. x = "Time (ms)",
  581. y = expression("Spectral Slope Difference [a.u.]")
  582. ) +
  583. #ggtitle("Wavelet: Difference in 1/f slope (Incongruent – Congruent)") +
  585. # Horizontal lines for significant time points (black for crit1)
  586. geom_segment(
  587. data = aov_congr %>%
  588. inner_join(pvals[pvals$Congr_crit1 == 1], by = c("variable", "time_chunk")),
  589. aes(x = time - 20, xend = time + 20,
  590. y = -.1,
  591. yend = -.1),
  592. inherit.aes = FALSE, color = "black", size = 1.5
  593. ) +
  595. # Horizontal lines for significant time points (dark green for crit2)
  596. geom_segment(
  597. data = aov_congr %>%
  598. inner_join(pvals[pvals$Congr_crit2 == 1], by = c("variable", "time_chunk")),
  599. aes(x = time - 20, xend = time + 20,
  600. y = -.12,
  601. yend = - .12),
  602. inherit.aes = FALSE, color = "darkgreen", size = 1.5
  603. )
  605. # ggsave(file.path(folder_name, "anova_congr.jpg"), width = 8, height = 15, dpi = 600, scale = .5)
  606. # ggsave(file.path(folder_name, "anova_congr.pdf"), width = 8, height = 15, dpi = 600, scale = .5)
  608. ###Previous Congr
  610. # Calculate difference between Incongruent and Congruent
  611. aov_pcongr <- reg_pcongr2 %>%
  612. pivot_wider(names_from = PrevCongr, values_from = value) %>%
  613. mutate(value = pInc - pCon) %>%
  614. select(time, time_chunk, variable, value)
  616. ggplot(aov_pcongr, aes(x = time, y = value)) +
  617. facet_wrap(. ~ variable, ncol = 1, drop = FALSE) +
  619. # Add shaded regions
  620. geom_rect(data = shading_data,
  621. aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
  622. fill = factor(brightness)),
  623. alpha = 0.5, inherit.aes = FALSE, show.legend = FALSE) +
  624. scale_fill_manual(values = scales::alpha(
  625. c("lightgrey", "lightgrey", "grey", "lightgrey", "grey",
  626. "lightgrey", "grey", "lightgrey", "grey", "lightgrey", "grey"), 0.5)) +
  628. # Mean ± SE ribbon and mean line
  629. stat_summary(fun.data = mean_se, geom = "ribbon", alpha = .1, fill = "grey50") +
  630. stat_summary(fun = "mean", geom = "line", size = 1.2, color = "black") +
  632. # Theme and layout
  633. theme_bw(base_size = 14) +
  634. theme(
  635. panel.grid = element_blank(),
  636. strip.background = element_rect(fill = "white"),
  637. strip.text = element_text(size = 12, face = "bold"),
  638. axis.title = element_text(size = 12),
  639. axis.text = element_text(size = 10),
  640. plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
  641. legend.position = "none"
  642. ) +
  644. # Axes and labels
  645. geom_hline(yintercept = 0, linetype = "dashed", color = "black", size = 0.3) +
  646. geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.3) +
  647. coord_cartesian(ylim = c(-.12, .1)) +
  648. labs(
  649. x = "Time (ms)",
  650. y = expression("Spectral Slope Difference [a.u.]")
  651. ) +
  652. #ggtitle("Wavelet: Difference in 1/f slope (Incongruent – Congruent)") +
  654. # Horizontal lines for significant time points (black for crit1)
  655. geom_segment(
  656. data = aov_pcongr %>%
  657. inner_join(pvals[pvals$PrevCongr_crit1 == 1], by = c("variable", "time_chunk")),
  658. aes(x = time - 20, xend = time + 20,
  659. y = -.1,
  660. yend = -.1),
  661. inherit.aes = FALSE, color = "black", size = 1.5
  662. ) +
  664. # Horizontal lines for significant time points (dark green for crit2)
  665. geom_segment(
  666. data = aov_pcongr %>%
  667. inner_join(pvals[pvals$PrevCongr_crit2 == 1], by = c("variable", "time_chunk")),
  668. aes(x = time - 20, xend = time + 20,
  669. y = -.12,
  670. yend = - .12),
  671. inherit.aes = FALSE, color = "darkgreen", size = 1.5
  672. )
  674. # ggsave(file.path(folder_name, "anova_pcongr.jpg"), width = 8, height = 15, dpi = 600, scale = .5)
  675. # ggsave(file.path(folder_name, "anova_pcongr.pdf"), width = 8, height = 15, dpi = 600, scale = .5)
  677. ###Interaction
  679. aov_pcongr_by_congr <- reg_pcongr2 %>%
  680. pivot_wider(names_from = PrevCongr, values_from = value) %>%
  681. mutate(value = pInc - pCon) %>%
  682. select(time, time_chunk, variable, Congr, value)
  684. ggplot(aov_pcongr_by_congr, aes(x = time, y = value, color = Congr)) +
  685. facet_wrap(. ~ variable, ncol = 1, drop = FALSE) +
  687. # Add shaded regions
  688. geom_rect(data = shading_data,
  689. aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
  690. fill = factor(brightness)),
  691. alpha = 0.5, inherit.aes = FALSE, show.legend = FALSE) +
  692. scale_fill_manual(values = scales::alpha(
  693. c("lightgrey", "lightgrey", "grey", "lightgrey", "grey",
  694. "lightgrey", "grey", "lightgrey", "grey", "lightgrey", "grey"), 0.5)) +
  696. # Mean ± SE and mean line per Congr level
  697. stat_summary(fun.data = mean_se, geom = "ribbon", aes(fill = Congr), alpha = .1, color = NA) +
  698. stat_summary(fun = mean, geom = "line", size = 1.2) +
  700. # Color mapping
  701. scale_color_manual(values = c("Con" = "blue", "Inc" = "red")) +
  703. # Theme
  704. theme_bw(base_size = 14) +
  705. theme(
  706. panel.grid = element_blank(),
  707. strip.background = element_rect(fill = "white"),
  708. strip.text = element_text(size = 12, face = "bold"),
  709. axis.title = element_text(size = 12),
  710. axis.text = element_text(size = 10),
  711. plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
  712. legend.position = "none"
  713. ) +
  715. # Axes and labels
  716. geom_hline(yintercept = 0, linetype = "dashed", color = "black", size = 0.3) +
  717. geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.3) +
  718. coord_cartesian(ylim = c(-.12, .1)) +
  719. labs(
  720. x = "Time (ms)",
  721. y = expression("Spectral Slope Difference [a.u.]"),
  722. color = "Congruency"
  723. ) +
  725. # Significance for crit1
  726. geom_segment(
  727. data = aov_pcongr_by_congr %>%
  728. inner_join(pvals[pvals$PrevCongrCongr_crit1 == 1], by = c("variable", "time_chunk")),
  729. aes(x = time - 20, xend = time + 20, y = -.1, yend = -.1),
  730. inherit.aes = FALSE, color = "black", size = 1.5
  731. ) +
  733. # Significance for crit2
  734. geom_segment(
  735. data = aov_pcongr_by_congr %>%
  736. inner_join(pvals[pvals$PrevCongrCongr_crit2 == 1], by = c("variable", "time_chunk")),
  737. aes(x = time - 20, xend = time + 20, y = -.12, yend = -.12),
  738. inherit.aes = FALSE, color = "darkgreen", size = 1.5
  739. )
  741. # ggsave(file.path(folder_name, "anova_pcongrcongr.jpg"), width = 8, height = 15, dpi = 600, scale = .5)
  742. # ggsave(file.path(folder_name, "anova_pcongrcongr.pdf"), width = 8, height = 15, dpi = 600, scale = .5)
  745. aov_pcongr_by_congr1 <- reg_pcongr2 %>%
  746. pivot_wider(names_from = Congr, values_from = value) %>%
  747. mutate(value = Inc - Con) %>%
  748. select(time, time_chunk, variable, PrevCongr, value)
  751. ggplot(aov_pcongr_by_congr1, aes(x = time, y = value, color = PrevCongr)) +
  752. facet_wrap(. ~ variable, ncol = 1, drop = FALSE) +
  754. # Add shaded regions
  755. geom_rect(data = shading_data,
  756. aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
  757. fill = factor(brightness)),
  758. alpha = 0.5, inherit.aes = FALSE, show.legend = FALSE) +
  759. scale_fill_manual(values = scales::alpha(
  760. c("lightgrey", "lightgrey", "grey", "lightgrey", "grey",
  761. "lightgrey", "grey", "lightgrey", "grey", "lightgrey", "grey"), 0.5)) +
  763. # Mean ± SE and mean line per Congr level
  764. stat_summary(fun.data = mean_se, geom = "ribbon", aes(fill = PrevCongr), alpha = .1, color = NA) +
  765. stat_summary(fun = mean, geom = "line", size = 1.2) +
  767. # Color mapping
  768. scale_color_manual(values = c("pCon" = "blue", "pInc" = "red")) +
  770. # Theme
  771. theme_bw(base_size = 14) +
  772. theme(
  773. panel.grid = element_blank(),
  774. strip.background = element_rect(fill = "white"),
  775. strip.text = element_text(size = 12, face = "bold"),
  776. axis.title = element_text(size = 12),
  777. axis.text = element_text(size = 10),
  778. plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
  779. legend.position = "none"
  780. ) +
  782. # Axes and labels
  783. geom_hline(yintercept = 0, linetype = "dashed", color = "black", size = 0.3) +
  784. geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.3) +
  785. coord_cartesian(ylim = c(-.12, .1)) +
  786. labs(
  787. x = "Time (ms)",
  788. y = expression("Spectral Slope Difference [a.u.]"),
  789. color = "Congruency"
  790. ) +
  792. # Significance for crit1
  793. geom_segment(
  794. data = aov_pcongr_by_congr %>%
  795. inner_join(pvals[pvals$PrevCongrCongr_crit1 == 1], by = c("variable", "time_chunk")),
  796. aes(x = time - 20, xend = time + 20, y = -.1, yend = -.1),
  797. inherit.aes = FALSE, color = "black", size = 1.5
  798. ) +
  800. # Significance for crit2
  801. geom_segment(
  802. data = aov_pcongr_by_congr %>%
  803. inner_join(pvals[pvals$PrevCongrCongr_crit2 == 1], by = c("variable", "time_chunk")),
  804. aes(x = time - 20, xend = time + 20, y = -.12, yend = -.12),
  805. inherit.aes = FALSE, color = "darkgreen", size = 1.5
  806. )
  808. # ggsave(file.path(folder_name, "anova_pcongrcongr1.jpg"), width = 8, height = 15, dpi = 600, scale = .5)
  809. # ggsave(file.path(folder_name, "anova_pcongrcongr1.pdf"), width = 8, height = 15, dpi = 600, scale = .5)
  811. ###Interaction combined
  812. head(aov_pcongr_by_congr)
  813. head(aov_pcongr_by_congr1)
  815. colnames(aov_pcongr_by_congr)[4] <- "Cond"
  816. colnames(aov_pcongr_by_congr1)[4] <- "Cond"
  818. aov_pcongr_by_congr2 <- rbind(aov_pcongr_by_congr, aov_pcongr_by_congr1)
  820. ggplot(aov_pcongr_by_congr2, aes(x = time, y = value, color = Cond)) +
  821. facet_wrap(. ~ variable, ncol = 1, drop = FALSE) +
  823. # Add shaded regions
  824. geom_rect(data = shading_data,
  825. aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
  826. fill = factor(brightness)),
  827. alpha = 0.5, inherit.aes = FALSE, show.legend = FALSE) +
  828. scale_fill_manual(values = scales::alpha(
  829. c("lightgrey", "lightgrey", "grey", "lightgrey", "grey",
  830. "lightgrey", "grey", "lightgrey", "grey", "lightgrey", "grey"), 0.5)) +
  832. # Mean ± SE and mean line per Congr level
  833. stat_summary(fun.data = mean_se, geom = "ribbon", aes(fill = Cond), alpha = .1, color = NA) +
  834. stat_summary(fun = mean, geom = "line", size = 1.2) +
  836. # Color mapping
  837. scale_color_manual(values = c("pCon" = "#56B4E9", "pInc" = "#0072B2", "Con" = "#E69F00", "Inc" = "#D55E00")) +
  839. # Theme
  840. theme_bw(base_size = 14) +
  841. theme(
  842. panel.grid = element_blank(),
  843. strip.background = element_rect(fill = "white"),
  844. strip.text = element_text(size = 12, face = "bold"),
  845. axis.title = element_text(size = 12),
  846. axis.text = element_text(size = 10),
  847. plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
  848. legend.position = "none"
  849. ) +
  851. # Axes and labels
  852. geom_hline(yintercept = 0, linetype = "dashed", color = "black", size = 0.3) +
  853. geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.3) +
  854. coord_cartesian(ylim = c(-.12, .1)) +
  855. labs(
  856. x = "Time (ms)",
  857. y = expression("Spectral Slope Difference [a.u.]"),
  858. color = "Congruency"
  859. ) +
  861. # Significance for crit1
  862. geom_segment(
  863. data = aov_pcongr_by_congr %>%
  864. inner_join(pvals[pvals$PrevCongrCongr_crit1 == 1], by = c("variable", "time_chunk")),
  865. aes(x = time - 20, xend = time + 20, y = -.1, yend = -.1),
  866. inherit.aes = FALSE, color = "black", size = 1.5
  867. ) +
  869. # Significance for crit2
  870. geom_segment(
  871. data = aov_pcongr_by_congr %>%
  872. inner_join(pvals[pvals$PrevCongrCongr_crit2 == 1], by = c("variable", "time_chunk")),
  873. aes(x = time - 20, xend = time + 20, y = -.12, yend = -.12),
  874. inherit.aes = FALSE, color = "darkgreen", size = 1.5
  875. )
  877. # ggsave(file.path(folder_name, "anova_pcongrcongr2.jpg"), width = 8, height = 15, dpi = 600, scale = .5)
  878. # ggsave(file.path(folder_name, "anova_pcongrcongr2.pdf"), width = 8, height = 15, dpi = 600, scale = .5)
  881. ```
  883. Follow-ups for rmANOVA
  885. ```{r}
  886. #pCon vs pInc within Congr level, for selected ROI × time bin
  887. unique(time_chunk_levels) # sanity-check available bins
  889. #Boxplot to show distribution
  890. plot_box_all_conditions <- function(dat_roi,
  891. ylab = "Spectral Slope [a.u.]",
  892. ylim = c()) {
  894. dat_roi$Cond <- interaction(dat_roi$Congr, dat_roi$PrevCongr,
  895. sep = "-",
  896. lex.order = TRUE)
  898. dat_roi$Cond <- factor(dat_roi$Cond,
  899. levels = c("Con-pCon",
  900. "Con-pInc",
  901. "Inc-pCon",
  902. "Inc-pInc"))
  904. ggplot(dat_roi, aes(x = Cond, y = value, fill = Congr)) +
  906. facet_grid(. ~ variable + time_chunk) +
  908. geom_boxplot(
  909. width = .6,
  910. alpha = .6,
  911. outlier.shape = NA
  912. ) +
  914. geom_jitter(
  915. width = .08,
  916. alpha = .35,
  917. size = 1.5
  918. ) +
  920. stat_summary(
  921. fun = mean,
  922. geom = "point",
  923. size = 3.5,
  924. shape = 21,
  925. fill = "white",
  926. color = "black"
  927. ) +
  929. scale_fill_manual(values = c("Con" = "red", "Inc" = "blue")) +
  931. theme_bw(base_size = 14) +
  932. theme(
  933. panel.grid = element_blank(),
  934. strip.background = element_rect(fill = "white"),
  935. strip.text = element_text(size = 12, face = "bold"),
  936. axis.title = element_text(size = 12),
  937. axis.text = element_text(size = 10)
  938. ) +
  940. labs(
  941. x = "Condition (Congruency × Previous Congruency)",
  942. y = ylab
  943. ) +
  945. coord_cartesian(ylim = ylim)
  946. }
  948. # ==============================================================================
  949. # 1) FRONTAL (MR2), time_chunk = (160,320]
  950. # ==============================================================================
  952. # -- Select the ROI × time bin -------------------------------------------------
  953. dat1 <- subset(dat, time_chunk == "(160,320]" & variable == "FRONTAL (MR2)")
  954. stopifnot(nrow(dat1) > 0)
  956. # -- Quick 2×2 cell means (Con/Inc × pCon/pInc) --------------------------------
  957. means_2x2 <- data.table::dcast(
  958. data.table::as.data.table(dat1),
  959. Congr + PrevCongr ~ .,
  960. value.var = "value",
  961. fun.aggregate = mean, na.rm = TRUE
  962. )
  963. data.table::setnames(means_2x2, ".", "mean")
  964. print(means_2x2)
  966. # -- Paired t-tests: pCon vs pInc within Congr = "Con" -------------------------
  967. w_con <- dat1[dat1$Congr == "Con", c("ID","PrevCongr","value")]
  968. w_con <- tidyr::pivot_wider(w_con, names_from = PrevCongr, values_from = value)
  969. w_con <- subset(w_con, !is.na(pCon) & !is.na(pInc)) # ensure the same IDs contribute to both levels
  970. w_con <- w_con[order(w_con$ID), ]
  971. t.test(w_con$pInc, w_con$pCon, paired = TRUE, alternative = "two.sided", conf.level = 0.95)
  973. # -- Paired t-tests: pCon vs pInc within Congr = "Inc" -------------------------
  974. w_inc <- dat1[dat1$Congr == "Inc", c("ID","PrevCongr","value")]
  975. w_inc <- tidyr::pivot_wider(w_inc, names_from = PrevCongr, values_from = value)
  976. w_inc <- subset(w_inc, !is.na(pCon) & !is.na(pInc))
  977. w_inc <- w_inc[order(w_inc$ID), ]
  978. t.test(w_inc$pInc, w_inc$pCon, paired = TRUE, alternative = "two.sided", conf.level = 0.95)
  980. # -- Plot (within-subject SE via Rmisc::summarySEwithin) ----------------------
  981. library(Rmisc) # summarySEwithin
  982. dat1$PrevCongr <- factor(dat1$PrevCongr, levels = c("pCon","pInc"))
  983. dat1$Congr <- factor(dat1$Congr, levels = c("Con","Inc"))
  985. summary_dat <- summarySEwithin(
  986. dat1,
  987. measurevar = "value",
  988. withinvars = c("variable","time_chunk","PrevCongr","Congr"),
  989. idvar = "ID"
  990. )
  992. ggplot(summary_dat, aes(x = PrevCongr, y = value, color = Congr, group = Congr)) +
  993. facet_grid(. ~ variable + time_chunk) +
  994. geom_line(position = position_dodge(0.2), size = 0.8) +
  995. geom_point(position = position_dodge(0.2), size = 3) +
  996. geom_errorbar(aes(ymin = value - se, ymax = value + se),
  997. width = 0.1, position = position_dodge(0.2)) +
  998. scale_color_manual(values = c("Con" = "red", "Inc" = "blue")) +
  999. theme_bw(base_size = 14) +
  1000. theme(
  1001. panel.grid = element_blank(),
  1002. strip.background = element_rect(fill = "white"),
  1003. strip.text = element_text(size = 12, face = "bold"),
  1004. axis.title = element_text(size = 12),
  1005. axis.text = element_text(size = 10),
  1006. plot.title = element_text(size = 14, face = "bold", hjust = 0.5)
  1007. ) +
  1008. labs(x = "Previous Congruency", y = "Spectral Slope [a.u.]", title = NULL) +
  1009. coord_cartesian(ylim = c(-1.35, -0.68))
  1010. # ggsave(file.path(folder_name, "ttest_frontal_160.jpg"), width = 7, height = 6, dpi = 600, scale = .5)
  1011. # ggsave(file.path(folder_name, "ttest_frontal_160.pdf"), width = 7, height = 6, dpi = 600, scale = .5)
  1013. plot_box_all_conditions(dat1)
  1015. # ==============================================================================
  1016. # 2) CENTRAL (MR1), time_chunk = (160,320]
  1017. # ==============================================================================
  1019. dat1 <- subset(dat, time_chunk == "(160,320]" & variable == "CENTRAL (MR1)")
  1020. stopifnot(nrow(dat1) > 0)
  1022. means_2x2 <- data.table::dcast(
  1023. data.table::as.data.table(dat1),
  1024. Congr + PrevCongr ~ .,
  1025. value.var = "value",
  1026. fun.aggregate = mean, na.rm = TRUE
  1027. )
  1028. data.table::setnames(means_2x2, ".", "mean")
  1029. print(means_2x2)
  1031. w_con <- dat1[dat1$Congr == "Con", c("ID","PrevCongr","value")]
  1032. w_con <- tidyr::pivot_wider(w_con, names_from = PrevCongr, values_from = value)
  1033. w_con <- subset(w_con, !is.na(pCon) & !is.na(pInc))
  1034. w_con <- w_con[order(w_con$ID), ]
  1035. t.test(w_con$pInc, w_con$pCon, paired = TRUE, alternative = "two.sided", conf.level = 0.95)
  1037. w_inc <- dat1[dat1$Congr == "Inc", c("ID","PrevCongr","value")]
  1038. w_inc <- tidyr::pivot_wider(w_inc, names_from = PrevCongr, values_from = value)
  1039. w_inc <- subset(w_inc, !is.na(pCon) & !is.na(pInc))
  1040. w_inc <- w_inc[order(w_inc$ID), ]
  1041. t.test(w_inc$pInc, w_inc$pCon, paired = TRUE, alternative = "two.sided", conf.level = 0.95)
  1043. library(Rmisc)
  1044. dat1$PrevCongr <- factor(dat1$PrevCongr, levels = c("pCon","pInc"))
  1045. dat1$Congr <- factor(dat1$Congr, levels = c("Con","Inc"))
  1047. summary_dat <- summarySEwithin(dat1, "value",
  1048. withinvars = c("variable","time_chunk","PrevCongr","Congr"),
  1049. idvar = "ID")
  1051. ggplot(summary_dat, aes(x = PrevCongr, y = value, color = Congr, group = Congr)) +
  1052. facet_grid(. ~ variable + time_chunk) +
  1053. geom_line(position = position_dodge(0.2), size = 0.8) +
  1054. geom_point(position = position_dodge(0.2), size = 3) +
  1055. geom_errorbar(aes(ymin = value - se, ymax = value + se),
  1056. width = 0.1, position = position_dodge(0.2)) +
  1057. scale_color_manual(values = c("Con" = "red", "Inc" = "blue")) +
  1058. theme_bw(base_size = 14) +
  1059. theme(
  1060. panel.grid = element_blank(),
  1061. strip.background = element_rect(fill = "white"),
  1062. strip.text = element_text(size = 12, face = "bold"),
  1063. axis.title = element_text(size = 12),
  1064. axis.text = element_text(size = 10),
  1065. plot.title = element_text(size = 14, face = "bold", hjust = 0.5)
  1066. ) +
  1067. labs(x = "Previous Congruency", y = "Spectral Slope [a.u.]", title = NULL) +
  1068. coord_cartesian(ylim = c(-1.35, -0.68))
  1069. # ggsave(file.path(folder_name, "ttest_central_160.jpg"), width = 7, height = 6, dpi = 600, scale = .5)
  1070. # ggsave(file.path(folder_name, "ttest_central_160.pdf"), width = 7, height = 6, dpi = 600, scale = .5)
  1072. plot_box_all_conditions(dat1)
  1074. # ==============================================================================
  1075. # 3) OCCIPITAL (MR3), time_chunk = (480,640]
  1076. # ==============================================================================
  1078. dat1 <- subset(dat, time_chunk == "(480,640]" & variable == "OCCIPITAL (MR3)")
  1079. stopifnot(nrow(dat1) > 0)
  1081. means_2x2 <- data.table::dcast(
  1082. data.table::as.data.table(dat1),
  1083. Congr + PrevCongr ~ .,
  1084. value.var = "value",
  1085. fun.aggregate = mean, na.rm = TRUE
  1086. )
  1087. data.table::setnames(means_2x2, ".", "mean")
  1088. print(means_2x2)
  1090. w_con <- dat1[dat1$Congr == "Con", c("ID","PrevCongr","value")]
  1091. w_con <- tidyr::pivot_wider(w_con, names_from = PrevCongr, values_from = value)
  1092. w_con <- subset(w_con, !is.na(pCon) & !is.na(pInc))
  1093. w_con <- w_con[order(w_con$ID), ]
  1094. t.test(w_con$pInc, w_con$pCon, paired = TRUE, alternative = "two.sided", conf.level = 0.95)
  1096. w_inc <- dat1[dat1$Congr == "Inc", c("ID","PrevCongr","value")]
  1097. w_inc <- tidyr::pivot_wider(w_inc, names_from = PrevCongr, values_from = value)
  1098. w_inc <- subset(w_inc, !is.na(pCon) & !is.na(pInc))
  1099. w_inc <- w_inc[order(w_inc$ID), ]
  1100. t.test(w_inc$pInc, w_inc$pCon, paired = TRUE, alternative = "two.sided", conf.level = 0.95)
  1102. library(Rmisc)
  1103. dat1$PrevCongr <- factor(dat1$PrevCongr, levels = c("pCon","pInc"))
  1104. dat1$Congr <- factor(dat1$Congr, levels = c("Con","Inc"))
  1106. summary_dat <- summarySEwithin(dat1, "value",
  1107. withinvars = c("variable","time_chunk","PrevCongr","Congr"),
  1108. idvar = "ID")
  1110. ggplot(summary_dat, aes(x = PrevCongr, y = value, color = Congr, group = Congr)) +
  1111. facet_grid(. ~ variable + time_chunk) +
  1112. geom_line(position = position_dodge(0.2), size = 0.8) +
  1113. geom_point(position = position_dodge(0.2), size = 3) +
  1114. geom_errorbar(aes(ymin = value - se, ymax = value + se),
  1115. width = 0.1, position = position_dodge(0.2)) +
  1116. scale_color_manual(values = c("Con" = "red", "Inc" = "blue")) +
  1117. theme_bw(base_size = 14) +
  1118. theme(
  1119. panel.grid = element_blank(),
  1120. strip.background = element_rect(fill = "white"),
  1121. strip.text = element_text(size = 12, face = "bold"),
  1122. axis.title = element_text(size = 12),
  1123. axis.text = element_text(size = 10),
  1124. plot.title = element_text(size = 14, face = "bold", hjust = 0.5)
  1125. ) +
  1126. labs(x = "Previous Congruency", y = "Spectral Slope [a.u.]", title = NULL) +
  1127. coord_cartesian(ylim = c(-1.35, -0.68))
  1128. # ggsave(file.path(folder_name, "ttest_occ_480.jpg"), width = 7, height = 6, dpi = 600, scale = .5)
  1129. # ggsave(file.path(folder_name, "ttest_occ_480.pdf"), width = 7, height = 6, dpi = 600, scale = .5)
  1131. plot_box_all_conditions(dat1)
  1133. # ==============================================================================
  1134. # 4) FRONTAL (MR2), time_chunk = (960,1120]
  1135. # ==============================================================================
  1137. dat1 <- subset(dat, time_chunk == "(960,1120]" & variable == "FRONTAL (MR2)")
  1138. stopifnot(nrow(dat1) > 0)
  1140. means_2x2 <- data.table::dcast(
  1141. data.table::as.data.table(dat1),
  1142. Congr + PrevCongr ~ .,
  1143. value.var = "value",
  1144. fun.aggregate = mean, na.rm = TRUE
  1145. )
  1146. data.table::setnames(means_2x2, ".", "mean")
  1147. print(means_2x2)
  1149. w_con <- dat1[dat1$Congr == "Con", c("ID","PrevCongr","value")]
  1150. w_con <- tidyr::pivot_wider(w_con, names_from = PrevCongr, values_from = value)
  1151. w_con <- subset(w_con, !is.na(pCon) & !is.na(pInc))
  1152. w_con <- w_con[order(w_con$ID), ]
  1153. t.test(w_con$pInc, w_con$pCon, paired = TRUE, alternative = "two.sided", conf.level = 0.95)
  1155. w_inc <- dat1[dat1$Congr == "Inc", c("ID","PrevCongr","value")]
  1156. w_inc <- tidyr::pivot_wider(w_inc, names_from = PrevCongr, values_from = value)
  1157. w_inc <- subset(w_inc, !is.na(pCon) & !is.na(pInc))
  1158. w_inc <- w_inc[order(w_inc$ID), ]
  1159. t.test(w_inc$pInc, w_inc$pCon, paired = TRUE, alternative = "two.sided", conf.level = 0.95)
  1161. library(Rmisc)
  1162. dat1$PrevCongr <- factor(dat1$PrevCongr, levels = c("pCon","pInc"))
  1163. dat1$Congr <- factor(dat1$Congr, levels = c("Con","Inc"))
  1165. summary_dat <- summarySEwithin(dat1, "value",
  1166. withinvars = c("variable","time_chunk","PrevCongr","Congr"),
  1167. idvar = "ID")
  1169. ggplot(summary_dat, aes(x = PrevCongr, y = value, color = Congr, group = Congr)) +
  1170. facet_grid(. ~ variable + time_chunk) +
  1171. geom_line(position = position_dodge(0.2), size = 0.8) +
  1172. geom_point(position = position_dodge(0.2), size = 3) +
  1173. geom_errorbar(aes(ymin = value - se, ymax = value + se),
  1174. width = 0.1, position = position_dodge(0.2)) +
  1175. scale_color_manual(values = c("Con" = "red", "Inc" = "blue")) +
  1176. theme_bw(base_size = 14) +
  1177. theme(
  1178. panel.grid = element_blank(),
  1179. strip.background = element_rect(fill = "white"),
  1180. strip.text = element_text(size = 12, face = "bold"),
  1181. axis.title = element_text(size = 12),
  1182. axis.text = element_text(size = 10),
  1183. plot.title = element_text(size = 14, face = "bold", hjust = 0.5)
  1184. ) +
  1185. labs(x = "Previous Congruency", y = "Spectral Slope [a.u.]", title = NULL) +
  1186. coord_cartesian(ylim = c(-1.35, -0.68))
  1187. # ggsave(file.path(folder_name, "ttest_frontal_960.jpg"), width = 7, height = 6, dpi = 600, scale = .5)
  1188. # ggsave(file.path(folder_name, "ttest_frontal_960.pdf"), width = 7, height = 6, dpi = 600, scale = .5)
  1190. plot_box_all_conditions(dat1)
  1192. ```
  1193. Cluster-based permutations for global & experimental effects with figures
  1195. ```{r}
  1197. ## =========================================================
  1198. ## Unified Z-based cluster permutation pipeline + plotting
  1199. ## For:
  1200. ## A) Global stimulus-induced change (baseline-corrected; time>=0 vs time<0)
  1201. ## B) 2×2 effects (Congr, PrevCongr, Interaction) via within-subject contrasts
  1202. ##
  1203. ## Cluster permutation:
  1204. ## - Observed statistic: Z(t) = (obs_mean(t) - perm_mu(t)) / perm_sd(t)
  1205. ## - Permutations: sign-flip across subjects (within-subject scheme)
  1206. ## - Cluster-forming threshold: |Z| >= 1.6449
  1207. ## - Cluster statistic: max |Z| within cluster (max intensity)
  1208. ## - FWER control: null = max cluster-stat across time per permutation
  1209. ##
  1210. ## Outputs:
  1211. ## - Cluster table (t_start/t_end, max|Z|, mean Z, p_cluster, sig)
  1212. ## - Plots: Z(t) time course + shaded significant clusters (pFWER <= .05)
  1213. ## =========================================================
  1215. library(data.table)
  1216. library(ggplot2)
  1217. library(patchwork)
  1219. ## =========================================================
  1220. ## 0) GLOBAL SETTINGS (edit here)
  1221. ## =========================================================
  1222. Z_THR <- 1.6449
  1223. N_PERM <- 1000
  1224. SEED <- 123
  1225. POST_RANGE <- 0:1277
  1226. TIME_RANGE <- 0:1277
  1227. COMP_ORDER <- c("FRONTAL (MR2)", "CENTRAL (MR1)", "OCCIPITAL (MR3)")
  1229. ## =========================================================
  1230. ## 1) CORE HELPERS: clustering + permutation Z
  1231. ## =========================================================
  1233. clusters_maxint_1d <- function(stat, thr = Z_THR) {
  1234. above <- abs(stat) >= thr
  1235. if (!any(above)) return(list(idx = list(), cl_stat = numeric(0)))
  1237. r <- rle(above)
  1238. ends <- cumsum(r$lengths)
  1239. starts <- ends - r$lengths + 1
  1240. keep <- which(r$values)
  1242. idx <- lapply(keep, function(k) starts[k]:ends[k])
  1243. cl_stat <- vapply(idx, function(ix) max(abs(stat[ix]), na.rm = TRUE), numeric(1))
  1245. list(idx = idx, cl_stat = cl_stat)
  1246. }
  1248. perm_cluster_maxint_1d_Z <- function(X, n_perm = N_PERM, seed = SEED,
  1249. z_thr = Z_THR, keep_zperm = FALSE) {
  1250. # X: N × T matrix (subjects × time) of values (baseline-corrected or contrasts)
  1252. set.seed(seed)
  1253. N <- nrow(X)
  1254. Tt <- ncol(X)
  1256. # Observed mean map
  1257. obs_mean <- colMeans(X, na.rm = TRUE)
  1259. # Permuted mean maps (n_perm × T)
  1260. perm_means <- matrix(NA_real_, nrow = n_perm, ncol = Tt)
  1261. for (p in seq_len(n_perm)) {
  1262. signs <- sample(c(-1, 1), size = N, replace = TRUE)
  1263. perm_means[p, ] <- colMeans(X * signs, na.rm = TRUE)
  1264. }
  1266. # Pointwise standardization parameters from permutations
  1267. perm_mu <- colMeans(perm_means, na.rm = TRUE)
  1268. perm_sd <- apply(perm_means, 2, sd, na.rm = TRUE)
  1269. perm_sd[perm_sd == 0 | !is.finite(perm_sd)] <- NA_real_
  1271. # Observed Z map
  1272. z_obs <- (obs_mean - perm_mu) / perm_sd
  1274. # Optional permuted Z maps (for diagnostics / pointwise descriptive p(t))
  1275. z_perm <- NULL
  1276. if (isTRUE(keep_zperm)) {
  1277. z_perm <- sweep(perm_means, 2, perm_mu, "-")
  1278. z_perm <- sweep(z_perm, 2, perm_sd, "/")
  1279. }
  1281. # Observed clusters + stats
  1282. cl_obs <- clusters_maxint_1d(z_obs, thr = z_thr)
  1284. # Null distribution: max cluster stat per permutation
  1285. max_null <- numeric(n_perm)
  1286. for (p in seq_len(n_perm)) {
  1287. zp <- if (!is.null(z_perm)) z_perm[p, ] else (perm_means[p, ] - perm_mu) / perm_sd
  1288. cl_p <- clusters_maxint_1d(zp, thr = z_thr)
  1289. max_null[p] <- if (length(cl_p$cl_stat)) max(cl_p$cl_stat) else 0
  1290. }
  1292. # Cluster p-values (FWER-controlled)
  1293. p_cluster <- if (length(cl_obs$cl_stat)) {
  1294. vapply(cl_obs$cl_stat, function(v) (sum(max_null >= v) + 1) / (n_perm + 1), numeric(1))
  1295. } else numeric(0)
  1297. list(
  1298. z_obs = z_obs,
  1299. z_perm = z_perm,
  1300. clusters = cl_obs$idx,
  1301. cluster_stat = cl_obs$cl_stat,
  1302. cluster_p = p_cluster,
  1303. max_null = max_null,
  1304. z_thr = z_thr,
  1305. n_perm = n_perm,
  1306. seed = seed
  1307. )
  1308. }
  1310. summarize_clusters <- function(time_num, res, variable, effect = NA_character_) {
  1311. if (length(res$clusters) == 0) {
  1312. return(data.table(
  1313. variable = variable,
  1314. effect = effect,
  1315. cluster = integer(),
  1316. t_start = integer(),
  1317. t_end = integer(),
  1318. cl_maxabs_z = numeric(),
  1319. mean_z = numeric(),
  1320. mean_abs_z = numeric(),
  1321. p_cluster = numeric(),
  1322. sig = logical()
  1323. ))
  1324. }
  1326. out <- data.table(
  1327. variable = variable,
  1328. effect = effect,
  1329. cluster = seq_along(res$clusters),
  1330. t_start = time_num[vapply(res$clusters, min, integer(1))],
  1331. t_end = time_num[vapply(res$clusters, max, integer(1))],
  1332. cl_maxabs_z = res$cluster_stat,
  1333. mean_z = vapply(res$clusters, function(ix) mean(res$z_obs[ix], na.rm = TRUE), numeric(1)),
  1334. mean_abs_z = vapply(res$clusters, function(ix) mean(abs(res$z_obs[ix]), na.rm = TRUE), numeric(1)),
  1335. p_cluster = res$cluster_p
  1336. )
  1337. out[, sig := (p_cluster <= 0.05)]
  1338. out[]
  1339. }
  1341. plot_z_timecourse_with_clusters <- function(time_num, z_obs, clusters, cluster_p,
  1342. z_thr = Z_THR,
  1343. title = NULL,
  1344. subtitle = NULL) {
  1345. df <- data.table(time = time_num, z = as.numeric(z_obs))
  1347. shade <- NULL
  1348. if (length(clusters) > 0) {
  1349. sig_idx <- which(cluster_p <= 0.05)
  1350. if (length(sig_idx) > 0) {
  1351. shade <- rbindlist(lapply(sig_idx, function(k) {
  1352. ix <- clusters[[k]]
  1353. data.table(xmin = time_num[min(ix)], xmax = time_num[max(ix)])
  1354. }))
  1355. }
  1356. }
  1358. ggplot(df, aes(x = time, y = z)) +
  1359. { if (!is.null(shade)) geom_rect(
  1360. data = shade,
  1361. aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf),
  1362. inherit.aes = FALSE,
  1363. alpha = .15
  1364. )
  1365. } +
  1366. geom_hline(yintercept = c(-z_thr, z_thr), linetype = "dotted") +
  1367. geom_hline(yintercept = 0, linetype = "solid") +
  1368. geom_line(linewidth = 1) +
  1369. theme_bw(base_size = 14) +
  1370. theme(panel.grid = element_blank()) +
  1371. labs(x = "Time (ms)", y = "Observed Z(t)", title = title, subtitle = subtitle)
  1372. }
  1374. ## =========================================================
  1375. ## 2) A) GLOBAL: baseline-corrected post-stimulus change
  1376. ## Input: reg_av must have columns ID, time, variable, value
  1377. ## =========================================================
  1379. run_global_baselinecorr_Z <- function(reg_av,
  1380. post_time_range = POST_RANGE,
  1381. n_perm = N_PERM,
  1382. seed = SEED,
  1383. z_thr = Z_THR) {
  1385. DT <- as.data.table(reg_av)
  1386. DT[, time := as.integer(time)]
  1387. setorder(DT, ID, variable, time)
  1389. # Baseline per subject × variable (time < 0)
  1390. BL <- DT[time < 0, .(baseline = mean(value, na.rm = TRUE)), by = .(ID, variable)]
  1392. # Post-stimulus only
  1393. POST <- DT[time >= 0]
  1394. POST <- merge(POST, BL, by = c("ID", "variable"), all.x = TRUE)
  1396. # Baseline-corrected values
  1397. POST[, value_bc := value - baseline]
  1399. # Restrict post-stim time range
  1400. if (!is.null(post_time_range)) POST <- POST[time %in% post_time_range]
  1402. vars <- unique(POST$variable)
  1404. # Cluster table
  1405. cluster_table <- rbindlist(lapply(vars, function(v) {
  1407. tmp <- POST[variable == v, .(ID, time, val = value_bc)]
  1408. wide <- dcast(tmp, ID ~ time, value.var = "val")
  1410. time_cols <- setdiff(names(wide), "ID")
  1411. time_num <- as.integer(time_cols)
  1412. o <- order(time_num)
  1413. time_cols <- time_cols[o]
  1414. time_num <- time_num[o]
  1416. X <- as.matrix(wide[, ..time_cols])
  1418. res <- perm_cluster_maxint_1d_Z(X, n_perm = n_perm, seed = seed, z_thr = z_thr, keep_zperm = FALSE)
  1420. summarize_clusters(time_num, res, variable = v, effect = "global_bc")
  1421. }), use.names = TRUE, fill = TRUE)
  1423. # Z curves for plotting
  1424. z_curves <- lapply(vars, function(v) {
  1426. tmp <- POST[variable == v, .(ID, time, val = value_bc)]
  1427. wide <- dcast(tmp, ID ~ time, value.var = "val")
  1429. time_cols <- setdiff(names(wide), "ID")
  1430. time_num <- as.integer(time_cols)
  1431. o <- order(time_num)
  1432. time_cols <- time_cols[o]
  1433. time_num <- time_num[o]
  1435. X <- as.matrix(wide[, ..time_cols])
  1437. res <- perm_cluster_maxint_1d_Z(X, n_perm = n_perm, seed = seed, z_thr = z_thr, keep_zperm = FALSE)
  1439. list(variable = v, time = time_num, res = res)
  1440. })
  1441. names(z_curves) <- vars
  1443. list(
  1444. POST = POST,
  1445. cluster_table = cluster_table,
  1446. z_curves = z_curves,
  1447. settings = list(n_perm = n_perm, seed = seed, z_thr = z_thr)
  1448. )
  1449. }
  1451. ## =========================================================
  1452. ## 3) B) 2×2 EFFECTS: Congr / PrevCongr / Interaction
  1453. ## Input: reg_pcongr2 must have columns:
  1454. ## ID, variable, time, PrevCongr, Congr, value
  1455. ## =========================================================
  1457. run_2x2_effects_Z <- function(reg_pcongr2,
  1458. time_range = TIME_RANGE,
  1459. n_perm = N_PERM,
  1460. seed = SEED,
  1461. z_thr = Z_THR) {
  1463. DT <- as.data.table(reg_pcongr2)
  1464. DT[, Congr := factor(Congr, levels = c("Con","Inc"))]
  1465. DT[, PrevCongr := factor(PrevCongr, levels = c("pCon","pInc"))]
  1466. DT[, time := as.integer(time)]
  1467. setorder(DT, ID, variable, time, PrevCongr, Congr)
  1469. if (!is.null(time_range)) DT <- DT[time %in% time_range]
  1471. # Wide 2×2 cells
  1472. W <- dcast(
  1473. DT,
  1474. ID + variable + time ~ PrevCongr + Congr,
  1475. value.var = "value"
  1476. )
  1478. required_cols <- c("pCon_Con","pCon_Inc","pInc_Con","pInc_Inc")
  1479. miss <- setdiff(required_cols, names(W))
  1480. if (length(miss) > 0) stop("Missing condition columns after dcast: ", paste(miss, collapse = ", "))
  1482. # Within-subject contrasts
  1483. EFF <- W[, .(
  1484. eff_Congr = 0.5 * ((pCon_Inc - pCon_Con) + (pInc_Inc - pInc_Con)),
  1485. eff_Prev = 0.5 * ((pInc_Con - pCon_Con) + (pInc_Inc - pCon_Inc)),
  1486. eff_Int = (pCon_Inc - pCon_Con) - (pInc_Inc - pInc_Con)
  1487. ), by = .(ID, variable, time)]
  1489. effects <- c("eff_Congr","eff_Prev","eff_Int")
  1490. vars <- unique(EFF$variable)
  1492. # Cluster table
  1493. cluster_table <- rbindlist(lapply(vars, function(v) {
  1494. rbindlist(lapply(effects, function(eff) {
  1496. tmp <- EFF[variable == v, .(ID, time, val = get(eff))]
  1497. wide <- dcast(tmp, ID ~ time, value.var = "val")
  1499. time_cols <- setdiff(names(wide), "ID")
  1500. time_num <- as.integer(time_cols)
  1501. o <- order(time_num)
  1502. time_cols <- time_cols[o]
  1503. time_num <- time_num[o]
  1505. X <- as.matrix(wide[, ..time_cols])
  1507. res <- perm_cluster_maxint_1d_Z(X, n_perm = n_perm, seed = seed, z_thr = z_thr, keep_zperm = FALSE)
  1509. summarize_clusters(time_num, res, variable = v, effect = eff)
  1510. }))
  1511. }), use.names = TRUE, fill = TRUE)
  1513. # Z curves for plotting
  1514. z_curves <- list()
  1515. for (v in vars) {
  1516. for (eff in effects) {
  1518. tmp <- EFF[variable == v, .(ID, time, val = get(eff))]
  1519. wide <- dcast(tmp, ID ~ time, value.var = "val")
  1521. time_cols <- setdiff(names(wide), "ID")
  1522. time_num <- as.integer(time_cols)
  1523. o <- order(time_num)
  1524. time_cols <- time_cols[o]
  1525. time_num <- time_num[o]
  1527. X <- as.matrix(wide[, ..time_cols])
  1529. res <- perm_cluster_maxint_1d_Z(X, n_perm = n_perm, seed = seed, z_thr = z_thr, keep_zperm = FALSE)
  1531. key <- paste(v, eff, sep = " | ")
  1532. z_curves[[key]] <- list(variable = v, effect = eff, time = time_num, res = res)
  1533. }
  1534. }
  1536. list(
  1537. EFF = EFF,
  1538. cluster_table = cluster_table,
  1539. z_curves = z_curves,
  1540. settings = list(n_perm = n_perm, seed = seed, z_thr = z_thr)
  1541. )
  1542. }
  1544. ## =========================================================
  1545. ## 4) FIGURE BUILDERS
  1546. ## =========================================================
  1548. make_effect_3panel_Z <- function(effects_out,
  1549. effect = c("eff_Congr","eff_Prev","eff_Int"),
  1550. components = COMP_ORDER,
  1551. z_thr = NULL,
  1552. title = NULL) {
  1554. effect <- match.arg(effect)
  1555. if (is.null(z_thr)) z_thr <- effects_out$settings$z_thr
  1556. if (is.null(title)) title <- paste0("Observed Z(t): ", effect, " (shaded clusters pFWER ≤ .05)")
  1558. panels <- lapply(components, function(v) {
  1559. key <- paste(v, effect, sep = " | ")
  1560. if (!key %in% names(effects_out$z_curves)) stop("Missing z_curve for: ", key)
  1562. item <- effects_out$z_curves[[key]]
  1564. plot_z_timecourse_with_clusters(
  1565. time_num = item$time,
  1566. z_obs = item$res$z_obs,
  1567. clusters = item$res$clusters,
  1568. cluster_p = item$res$cluster_p,
  1569. z_thr = z_thr,
  1570. title = v,
  1571. subtitle = NULL
  1572. ) +
  1573. theme(
  1574. plot.title = element_text(size = 12, face = "bold"),
  1575. axis.title.x = element_text(size = 11),
  1576. axis.title.y = element_text(size = 11)
  1577. )
  1578. })
  1580. wrap_plots(panels, ncol = 1) +
  1581. plot_annotation(title = title) &
  1582. theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
  1583. }
  1585. make_global_3panel_Z <- function(global_out,
  1586. components = COMP_ORDER,
  1587. z_thr = NULL,
  1588. title = NULL) {
  1590. if (is.null(z_thr)) z_thr <- global_out$settings$z_thr
  1591. if (is.null(title)) title <- "Global stimulus-induced changes vs baseline (Observed Z(t); shaded clusters pFWER ≤ .05)"
  1593. panels <- lapply(components, function(v) {
  1594. if (!v %in% names(global_out$z_curves)) stop("Missing global z_curve for: ", v)
  1596. item <- global_out$z_curves[[v]]
  1598. plot_z_timecourse_with_clusters(
  1599. time_num = item$time,
  1600. z_obs = item$res$z_obs,
  1601. clusters = item$res$clusters,
  1602. cluster_p = item$res$cluster_p,
  1603. z_thr = z_thr,
  1604. title = v,
  1605. subtitle = NULL
  1606. ) +
  1607. theme(
  1608. plot.title = element_text(size = 12, face = "bold"),
  1609. axis.title.x = element_text(size = 11),
  1610. axis.title.y = element_text(size = 11)
  1611. )
  1612. })
  1614. wrap_plots(panels, ncol = 1) +
  1615. plot_annotation(title = title) &
  1616. theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
  1617. }
  1619. ## =========================================================
  1620. ## 5) RUN PIPELINES
  1621. ## =========================================================
  1623. # ---- A) Global (baseline-corrected) ----
  1624. global_out <- run_global_baselinecorr_Z(
  1625. reg_av,
  1626. post_time_range = POST_RANGE,
  1627. n_perm = N_PERM,
  1628. seed = SEED,
  1629. z_thr = Z_THR
  1630. )
  1632. global_out$cluster_table[]
  1633. global_out$cluster_table[sig == TRUE]
  1635. # ---- B) 2×2 effects ----
  1636. effects_out <- run_2x2_effects_Z(
  1637. reg_pcongr2,
  1638. time_range = TIME_RANGE,
  1639. n_perm = N_PERM,
  1640. seed = SEED,
  1641. z_thr = Z_THR
  1642. )
  1644. effects_out$cluster_table[]
  1645. effects_out$cluster_table[sig == TRUE]
  1647. ## =========================================================
  1648. ## 6) MAKE THE 3-PANEL FIGURES
  1649. ## =========================================================
  1651. fig_global_3panel <- make_global_3panel_Z(
  1652. global_out,
  1653. components = COMP_ORDER,
  1654. title = "Global post-stimulus deviations from baseline (Observed Z(t); shaded clusters pFWER ≤ .05)"
  1655. )
  1657. fig_congr_3panel <- make_effect_3panel_Z(
  1658. effects_out,
  1659. effect = "eff_Congr",
  1660. components = COMP_ORDER,
  1661. title = "Main effect of Congruency (Observed Z(t); shaded clusters pFWER ≤ .05)"
  1662. )
  1664. fig_prev_3panel <- make_effect_3panel_Z(
  1665. effects_out,
  1666. effect = "eff_Prev",
  1667. components = COMP_ORDER,
  1668. title = "Main effect of Previous Congruency (Observed Z(t); shaded clusters pFWER ≤ .05)"
  1669. )
  1671. fig_int_3panel <- make_effect_3panel_Z(
  1672. effects_out,
  1673. effect = "eff_Int",
  1674. components = COMP_ORDER,
  1675. title = "Congruency × Previous Congruency (Observed Z(t); shaded clusters pFWER ≤ .05)"
  1676. )
  1678. fig_global_3panel
  1679. ggsave("fit_global.pdf", width = 6, height = 15, dpi = 600, scale = .5)
  1680. fig_congr_3panel
  1681. ggsave("fit_congr.pdf", width = 6, height = 15, dpi = 600, scale = .5)
  1682. fig_prev_3panel
  1683. ggsave("fit_prev.pdf", width = 6, height = 15, dpi = 600, scale = .5)
  1684. fig_int_3panel
  1685. ggsave("fit_int.pdf", width = 6, height = 15, dpi = 600, scale = .5)
  1688. ```

code_afterR1.Rmd, no license · at the source

Overview

Authors: Virginia Tronelli1,2, Patrycja Kałamała3, Gabriele Gratton4,5, Monica Fabiani4,5, Mate Gyurkovics6, Kathy A Low5, Maurizio Codispoti2, Andrea De Cesarei2
ORCID iDs: Virginia Tronelli, Patrycja Kałamała, Gabriele Gratton, Monica Fabiani, Mate Gyurkovics, Kathy A Low, Maurizio Codispoti, Andrea De Cesarei
  1. Center for Mind/Brain Sciences (CIMeC), Universitá Degli Studi di Trento‐Rovereto, Trento, Italy
  2. Department of Psychology, Universitá Degli Studi di Bologna, Bologna, Italy
  3. Centre for Cognitive Science, Jagiellonian University, Kraków, Poland
  4. Department of Psychology, University of Illinois Urbana‐Champaign, Champaign, Illinois, USA
  5. Beckman Institute for Advanced Science and Technology, University of Illinois Urbana‐Champaign, Champaign, Illinois, USA
  6. School of Psychology, University of East Anglia, Norwich, UK
Institutions: University of Trento (Italy); University of Bologna (Italy); Jagiellonian University (Poland); University of Illinois Urbana-Champaign (United States); University of East Anglia (United Kingdom)
Journal: Psychophysiology, volume 63, issue 8, article e70376
Dates: received 17 November 2025; accepted 29 July 2026; published online 8 August 2026; in print August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1111/psyp.70376 · PMID 42568317 · PMCID PMC13451924 · OpenAlex W7201979157
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: EEG (modality), human (organism), cognitive (subfield)
Methods: Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions, Preprocessing, Evoked potentials, Complexity, fMRI & imaging
Keywords: aperiodic 1/f EEG, cognitive control, congruency effect, congruency sequence effect, neural noise, spatial distribution
MeSH: Attention*, Cerebral Cortex*, Evoked Potentials*, Executive Function*, Inhibition, Psychological*, Pattern Recognition, Visual*, Psychomotor Performance*, Adult, Electroencephalography, Female, Humans, Male, Young Adult (* major topic)
Topic: Neural and Behavioral Psychology Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: National Institute on Aging (RF1AG062666); NIA NIH HHS (RF1AG062666, RF1 AG062666)
Citations: cited by 1 paper (Europe PMC); 75 references in the paper

Abstract

Aperiodic neural activity (1/f EEG) has been proposed to reflect the balance between excitatory and inhibitory (E/I) inputs, with steeper spectral slopes reflecting increased inhibition and flatter slopes indicating excitation. This activity is also thought to reflect the temporal coordination of neural firing, offering insights into fundamental brain dynamics. Recent studies have shown that the 1/f slope is sensitive to stimulus onset, characterized by initial inhibitory shifts followed by excitatory rebounds, which may reflect cognitive control mechanisms involved in suppressing distractions and preparing goal‐directed responses. However, previous research has relied on fixed temporal windows and insufficient control of ERP contamination, limiting our understanding of rapid control dynamics. Here we used newly developed time‐resolved analyses to study 1/f spectral slope modulation during a Picture‐Word Interference task, focusing on two canonical cognitive control markers: the Congruency Effect (CE) and Congruency Sequence Effect (CSE). Forty‐nine participants categorized pictures while ignoring congruent or incongruent words. Behaviorally, we replicated robust CE and CSE patterns. Spectral slope analyses showed that incongruent trials elicited steeper slopes—consistent with increased inhibition—particularly in frontal and central regions, reflecting conflict‐related control engagement. Moreover, CSE analyses revealed dynamic slope modulations across frontal, central, and occipital components over time, suggesting control adjustments influenced by previous trial congruency. These results provide the first within‐trial time‐resolved evidence that aperiodic 1/f EEG activity can track both immediate conflict resolution and cognitive adjustments, offering a temporally sensitive neural marker of cognitive control, albeit with effects that are small in magnitude on average.

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 2 matches between paragraphs and lines of code.

OSF z3qsg

License: none: the authors keep all their rights
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Languages: R (2)
Size: 6 files, 2 scripts
Software Heritage: not checked
Found in: “Data Availability Statement”
Holds: 2 notebooks
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: data.table (2 files), eegUtils (2 files), psych (2 files), tidyverse (2 files), ggplot2 (1 file), patchwork (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers (HTTP 200)
  • 26 September 2026: the link answers (HTTP 200)
2 files
At the source: osf.io/z3qsg

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;
  • 2 scripts, each with its path and the digest of its content;
  • 2 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 Statement

Data and code necessary to reproduce the statistical analyses (Section 3.2.3) are available at https://osf.io/z3qsg.

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 3, 28 September 2026

  • Publisher: n/a → Wiley

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 6 keywords, 13 MeSH terms, 2 funders, 73 references.

Cite

This paper

Tronelli, V., Kałamała, P., Gratton, G., Fabiani, M., Gyurkovics, M., Low, K. A., Codispoti, M., & De Cesarei, A. (2026). Changes in Aperiodic (1/f Slope) Activity During a Picture-Word Interference Task: Effects of Congruency and Sequence Manipulations. Psychophysiology, 63(8), e70376. https://doi.org/10.1111/psyp.70376

BibTeX

@article{tronelli2026changes,
author = {Tronelli, Virginia and Kałamała, Patrycja and Gratton, Gabriele and Fabiani, Monica and Gyurkovics, Mate and Low, Kathy A and Codispoti, Maurizio and De Cesarei, Andrea},
title = {{Changes in Aperiodic (1/f Slope) Activity During a Picture-Word Interference Task: Effects of Congruency and Sequence Manipulations}},
journal = {Psychophysiology},
year = {2026},
month = aug,
volume = {63},
number = {8},
pages = {e70376},
publisher = {Wiley},
issn = {0048-5772},
doi = {10.1111/psyp.70376},
url = {https://doi.org/10.1111/psyp.70376},
pmid = {42568317},
pmcid = {PMC13451924}
}

RIS

TY - JOUR
AU - Tronelli, Virginia
AU - Kałamała, Patrycja
AU - Gratton, Gabriele
AU - Fabiani, Monica
AU - Gyurkovics, Mate
AU - Low, Kathy A
AU - Codispoti, Maurizio
AU - De Cesarei, Andrea
TI - Changes in Aperiodic (1/f Slope) Activity During a Picture-Word Interference Task: Effects of Congruency and Sequence Manipulations
T2 - Psychophysiology
J2 - Psychophysiology
PY - 2026
DA - 2026/08/01
VL - 63
IS - 8
SP - e70376
SN - 0048-5772
PB - Wiley
DO - 10.1111/psyp.70376
UR - https://doi.org/10.1111/psyp.70376
LA - en
ER -

CSL-JSON

{
"id": "10.1111/psyp.70376",
"type": "article-journal",
"title": "Changes in Aperiodic (1/f Slope) Activity During a Picture-Word Interference Task: Effects of Congruency and Sequence Manipulations",
"container-title": "Psychophysiology",
"author": [
{
"family": "Tronelli",
"given": "Virginia"
},
{
"family": "Kałamała",
"given": "Patrycja"
},
{
"family": "Gratton",
"given": "Gabriele"
},
{
"family": "Fabiani",
"given": "Monica"
},
{
"family": "Gyurkovics",
"given": "Mate"
},
{
"family": "Low",
"given": "Kathy A"
},
{
"family": "Codispoti",
"given": "Maurizio"
},
{
"family": "De Cesarei",
"given": "Andrea"
}
],
"container-title-short": "Psychophysiology",
"volume": "63",
"issue": "8",
"page": "e70376",
"DOI": "10.1111/psyp.70376",
"PMID": "42568317",
"PMCID": "PMC13451924",
"ISSN": "0048-5772",
"publisher": "Wiley",
"URL": "https://doi.org/10.1111/psyp.70376",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
1
]
]
}
}

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.1111/psyp.70272 [code]
How to Improve the Reliability of Aperiodic Parameter Estimates in M/EEG: A Method Comparison.
Journal: Psychophysiology
In common: EEG, 15 references, 4 authors
[2] doi:10.1016/j.isci.2026.116936 [code]
Beyond neural oscillations: Stress-related aperiodic activity and aperiodic-oscillatory spectral covariation.
Journal: iScience
In common: eegUtils, data.table, patchwork, 2 other tools, EEG, 5 references
[3] doi:10.1111/ejn.70255 [code]
A Systematic Review of Aperiodic Neural Activity in Clinical Investigations
Journal: n/a
In common: EEG, 9 references
[4] doi:10.1093/braincomms/fcag351 [code]
Time-resolved aperiodic dynamics in event segmentation in attention-deficit/hyperactivity disorder.
Journal: Brain communications
In common: ggplot2, tidyverse, EEG, cognitive, 6 references
[5] doi:10.7554/elife.108673 [code]
Adaptive behavior is guided by integrated representations of controlled and non-controlled information.
Journal: eLife
In common: EEG, 6 references
[6] doi:10.1111/ejn.70543 [code]
Neural Oscillations Track Subjective and Pupillary Arousal During Naturalistic Movie Viewing.
Journal: The European journal of neuroscience
In common: EEG, cognitive, 7 references
[7] doi:10.1371/journal.pone.0328318
Education shapes the link between EEG aperiodic components and cognitive aging.
Journal: PloS one
In common: EEG, 7 references
[8] doi:10.3390/bioengineering13030323 [code]
Scale-Free Neurodynamics as Functional Fingerprint of Brain Regions.
Journal: Bioengineering (Basel, Switzerland)
In common: 7 references
[9] doi:10.1093/cercor/bhag113 [code]
Long-term reliability and stability of parameterized resting state EEG: evidence from a five-year follow-up.
Journal: Cerebral cortex (New York, N.Y. : 1991)
In common: psych, ggplot2, tidyverse, EEG, 4 references
[10] doi:10.1162/imag.a.1298 [code]
Frontoparietal control-default mode connectivity predicts TMS effects on cognitive control.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: psych, data.table, patchwork, 2 other tools, cognitive, 1 reference

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.

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: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

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