OSCR

No evidence of neural feature-specific pre-activation during the prediction of an upcoming stimulus.

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] § Methods › Pipeline 3 (Fig. 1e, Theoretical Null) ↔ postprocessing.qmd, lines 651–714 · score 0.55 · theoretical accuracy, confusion matrices, transition matrices, scores, hypothesis, training
  2. [2] § Methods › Pipeline 1 (Fig. 1c & 1e, Replication) ↔ reproduce_and_reorder_omission.py, lines 108–158 · score 0.53 · cross validation, iteration, omission, midminus, midplus, reproduce

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

Quarto · 1,980 lines · 61 KB · no license · 1 match

  1. ---
  2. title: "No evidence of neural feature-specific pre-activation during the prediction of an upcoming stimulus"
  3. subtitle: "A reanalysis of Demarchi et al. (Nature Communications, 2019)"
  4. author:
  5. - name: Abdoun Oussama
  6. orcid: 0000-0001-5123-4956
  7. email: [email hidden]
  8. affiliation:
  9. - name: Centre de Rercherches en Neurosciences de Lyon
  10. city: Providence
  11. - Todorov Dimitrii
  12. - Poublan-couzardot Arnaud
  13. - Tirou Coumarane
  14. - Lutz Antoine
  15. - Vernet Marine
  16. - Quentin Romain
  18. email-obfuscation: javascript
  20. abstract: ""
  22. filters:
  23. - lightbox
  25. lightbox:
  26. match: auto
  27. effect: fade
  28. loop: true
  30. format:
  31. html:
  32. toc: true
  33. toc-depth: 3
  34. toc-location: left
  35. code-fold: true
  36. code-line-numbers: true
  37. code-summary: "Source code"
  38. code-tools:
  39. source: https://github.com/romquentin/predictive_activity/blob/main/postprocessing.qmd
  40. toggle: true
  41. caption: "Code options"
  42. fig-align: center
  43. embed-resources: true
  44. link-external-newwindow: true
  45. link-external-icon: true
  46. other-links:
  47. - text: Github project home page
  48. href: https://github.com/romquentin/predictive_activity
  49. - text: Classification data
  50. href: https://doi.org/10.5281/zenodo.17699784
  52. editor: visual
  53. bibliography: references.bib
  54. ---
  56. # README
  58. **Background information**
  60. See the associated Github project [home page](https://github.com/romquentin/predictive_activity){target="_blank"}.
  62. The present R Quarto document calculates the correlation between decoding accuracy and sequence entropy, performs cluster-based permutation statistical tests and Bayesian statistics on accuracy and correlation time-generalization matrices, and plots all the figures and supplementary figures of the paper.
  64. **Instructions for use**
  66. 1. Obtain the classification data either from [Zenodo](https://doi.org/10.5281/zenodo.17699784){target="_blank"}, or by running the Python scripts on raw MEG data (see the Github project [home page](https://github.com/romquentin/predictive_activity){target="_blank"} for details). Save it in the `./data` folder.
  68. 2. For a first quick run of the script, we recommend that you set the number of cluster permutations to a low value (e.g. 100) in [**Prepare** / Parameters / Global](postprocessing.html#sec-params-global). Set it to 10,000 for stable results matching the ones reported in the publication.
  72. # **Prepare**
  74. ## Libraries
  76. ```{r}
  77. library(expm) # exponentiate matrices
  78. library(patchwork) # combine plots
  79. library(tictoc) # time executions
  80. library(correlation)
  81. library(abind) # combine multi-dimensional arrays
  82. library(BayesFactor)
  84. library(magrittr) # special pipes
  85. library(tidyverse)
  87. library(reticulate)
  88. np <- import("numpy")
  89. os <- import("os")
  90. mne <- import("mne")
  91. ```
  93. ## Parameters
  95. ### Global {#sec-params-global}
  97. ```{r}
  98. path.root <- "./data/results"
  99. list.subj <- list.dirs(path.root, full.names = F, recursive = F)
  100. n.subj <- length(list.subj)
  102. time_offset <- 670
  103. n_permutations <- 100
  104. save <- FALSE
  105. ```
  108. ### Plotting
  110. ```{r}
  112. # Temporal resolution in ms (1000 / sampling frequency)
  113. dt <- 10
  115. # Colors
  116. col.pal.base <- rev(pals::brewer.rdbu(10))
  117. col.pal.diff <- rev(pals::brewer.brbg(10))
  118. col.pal.corr <- pals::ocean.curl(100)
  119. col.lines <- "grey40"
  120. col.clusters <- "grey10"
  122. # Breaks & Labels
  123. # --- function
  124. format_labels_acc <- function(labels) {
  125. labels[seq(2,length(labels),2)] <- ""
  126. labels <- str_replace(labels, "0.", ".")
  127. }
  128. # --- base
  129. acc.breaks.base <- seq(.21,.29, .005)
  130. acc.labels.base <- format_labels_acc(acc.breaks.base)
  131. # --- diff
  132. acc.breaks.diff <- seq(-.02,.02, .005)
  133. acc.labels.diff <- format_labels_acc(acc.breaks.diff)
  134. # --- correlations
  135. corr.breaks <- seq(-1,1,.2)
  136. corr.labels <- str_replace(corr.breaks, "0.", ".")
  138. # GGPLOT theme
  139. theme_set(theme_minimal())
  140. theme_update(text = element_text(family = "Arial"),
  141. plot.tag.position = c(0, 0.95),
  142. plot.tag = element_text(size = 14, hjust = 0, vjust = 1, face = "bold"))
  143. ```
  145. ### Transition matrices
  147. From @demarchi2019
  149. ```{r}
  151. mat.T <- list(RD = matrix(rep(.25, 16), ncol = 4, byrow = T),
  152. MM = matrix(), MP = matrix(), OR = matrix())
  154. mat.T$MM <- matrix(c(c(.25,0,.37,.38),
  155. c(.38,.25,0,.37),
  156. c(.37,.38,.25,0),
  157. c(0,.37,.38,.25)),
  158. ncol = 4, byrow = T)
  160. mat.T$MP <- matrix(c(c(.25,0,.15,.60),
  161. c(.60,.25,0,.15),
  162. c(.15,.60,.25,0),
  163. c(0,.15,.60,.25)),
  164. ncol = 4, byrow = T)
  166. mat.T$OR <- matrix(c(c(.25,0,0,.75),
  167. c(.75,.25,0,0),
  168. c(0,.75,.25,0),
  169. c(0,0,.75,.25)),
  170. ncol = 4, byrow = T)
  172. covCT <- function(mat.C, mat.T) {
  173. covar <- 0
  174. for (i in 1:dim(mat.C)[1]) {
  175. covar <- covar + mean((mat.C[i,]-mean(mat.C[i,]))*(mat.T[i,]-mean(mat.T[i,])))
  176. }
  177. return(covar)
  178. }
  179. ```
  181. ## Functions
  183. ### Matrix to df
  185. ```{r}
  187. matrix_to_df2 <- function(array, row.name, col.name, val.name) {
  188. array %>%
  189. as_tibble() %>%
  190. rownames_to_column(row.name) %>%
  191. pivot_longer(-row.name, names_to = col.name, values_to = val.name) %>%
  192. mutate("{col.name}" := str_remove(!!sym(col.name), "V")) %>%
  193. mutate_all(as.numeric) %>%
  194. return
  195. }
  197. # !!! USE as.data.frame.table() instead
  198. ```
  200. ### Batch data loading
  202. ```{r}
  203. load_loop_subj <- function(root = path.root, folder, condition, diag = F, crop=T) {
  205. # Individual scores
  206. if (folder == "omissions") {
  207. array.scores <- array(dim = c(141,141,n.subj))
  208. } else if (folder == "sounds") {
  209. array.scores <- array(dim = c(33,135,n.subj))
  210. }
  212. i <- 1
  213. for (s in list.subj) {
  214. file <- os$path$join(root, s, folder,
  215. paste0("cv_", condition, "_scores.npy"))
  217. # Proceed to next subject if file doesn't exist
  218. if (!os$path$isfile(file)) {
  219. cat(file, " does not exist !!!\n")
  220. next
  221. }
  223. # Load data otherwise
  224. tmp <- np$load(file)
  226. # Average across CV folds if present
  227. if (length(dim(tmp)) == 3) {
  228. tmp <- tmp %>% apply(c(2,3), mean) # average across CV folds
  229. }
  231. array.scores[,,i] <- tmp
  232. i <- i+1
  233. }
  235. # Average across participants
  236. t_train_offset <- case_when(folder == "omissions" ~ time_offset,
  237. folder == "sounds" ~ 0)
  238. t_test_offset <- time_offset
  240. df.mean <- array.scores %>%
  241. apply(c(1,2), mean) %>%
  242. matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "accuracy") %>%
  243. # --- convert indices to time
  244. mutate(t_train = dt*(t_train-1) - t_train_offset,
  245. t_test = dt*(t_test-1) - t_test_offset)
  247. # --- crop time if requested (default)
  248. if (crop) {
  249. df.mean %<>% filter(t_train>=0, (t_train-333) < dt)
  250. }
  252. return(list(array.scores, df.mean))
  253. }
  254. ```
  256. ### Correlations
  258. ```{r}
  260. ## Correlation function between two vectors
  261. correlations_vecs <- function(vec1, vec2) {
  262. cor(as.vector(vec1), as.vector(vec2))
  263. }
  265. # Bind time-generalization matrices from different entropy conditions into a single 3D array
  266. correlations_extract3D <- function(data, idx.lines, idx.subj) {
  267. arr.tmp <- data[[idx.lines[1],"array"]][[1]][,,idx.subj]
  268. for (i in 2:length(idx.lines)) {
  269. arr.tmp <- abind(arr.tmp, data[[idx.lines[i],"array"]][[1]][,,idx.subj], along = 3)
  270. }
  271. return(arr.tmp)
  272. }
  275. # Replace NAs by 0 and extreme values by +/- .99
  276. correlations_fixval <- function(arr.corr) {
  277. arr.corr[is.na(arr.corr)] <- 0
  278. arr.corr[arr.corr == 1] <- .99
  279. arr.corr[arr.corr == -1] <- -.99
  280. return(arr.corr)
  281. }
  282. ```
  284. ### Cluster permutation stats
  286. ```{r}
  288. get_signif_clust2d <- function(X, n_permutations = 2**12, n_jobs = -1) {
  289. # run permutaiton tests
  290. mne$stats$spatio_temporal_cluster_1samp_test(X,
  291. n_permutations = n_permutations,
  292. n_jobs = as.integer(n_jobs),
  293. out_type = "mask",
  294. verbose = FALSE) -> tmp
  295. names(tmp) <- c("t_obs","clusters","cluster_pv","H0")
  297. # Extract clusters: the final output is a 0/1 matrix, 1 indicating belonging to a significant cluster
  298. # --- get ids of significant clusters
  299. idx.signif <- which(tmp$cluster_pv < .05)
  300. # --- aggregate clusters
  301. signif <- array(0, dim = dim(tmp$t_obs))
  302. for (i in idx.signif) {
  303. signif <- signif + tmp$clusters[[i]]
  304. }
  306. return(1*(signif > 0))
  307. }
  310. # Loop through multiple conditions, stored in a dataframe where each row contains a 3D array in an "array" column
  311. # Other columns are appended to the output as descriptors of the condition
  312. cluster_loop_cond <- function(df, H0 = 0, n = 2**12, n_jobs = -1, crop = FALSE) {
  313. # Initialize significant clusters dataframe
  314. df.clusters <- tibble()
  315. # Loop through all conditions
  316. for (i in 1:nrow(df)) {
  318. # Extract array of subject-level time-generalized accuracy
  319. X <- df[[i,"array"]][[1]]
  320. # --- crop to stimulus N-1 if requested
  321. if (crop) {X <- X[,35:68,]}
  322. # --- reshape to how MNE expects it
  323. X <- X %>% np$moveaxis(as.integer(2), as.integer(0))
  325. df.clusters %<>% bind_rows(
  326. get_signif_clust2d(X-H0, n_jobs = n_jobs, n_permutations = n) %>%
  327. matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "signif") %>%
  328. bind_cols(select(df,-array)[i,])
  329. )
  330. }
  332. return(df.clusters)
  333. }
  334. ```
  336. ### Difference original-reordered
  338. ```{r}
  340. calculate_diff <- function(df, f.map, values_to) {
  341. df %>%
  342. mutate(manip = case_match(manip, "" ~ "original", "_reord" ~ "reordered")) %>%
  343. pivot_wider(names_from = manip, values_from = values_to) %>%
  344. mutate(diff = f.map(original, reordered, .f=\(x,y){x-y})) %>%
  345. pivot_longer(cols = c(original, reordered, diff), names_to = "manip", values_to = values_to)
  346. }
  347. ```
  349. ### Base plot
  351. ```{r}
  352. plot_base <- function(df.plot,
  353. z = "accuracy",
  354. contour = TRUE,
  355. col.contour = col.clusters,
  356. col.pal = rev(pals::brewer.rdbu(10)),
  357. z.breaks = acc.breaks.base,
  358. z.labels = acc.labels.base,
  359. legend.position = "bottom") {
  361. # Onset lines
  362. lines.v.pos <- 330*seq(round(min(df.plot$t_test)/333), round(max(df.plot$t_test)/333), 1)
  363. lines.v.pos <- lines.v.pos[lines.v.pos!=0]
  364. lines.h.pos <- 330*seq(round(min(df.plot$t_train)/333), round(max(df.plot$t_train)/333), 1)
  365. lines.h.pos <- lines.h.pos[lines.h.pos!=0]
  367. # Plot
  368. g <- df.plot %>%
  369. ggplot(aes(x = t_test, y = t_train)) +
  370. # --- heatmap
  371. geom_tile(aes_string(fill = z)) +
  372. # --- sounds onsets
  373. geom_vline(xintercept = lines.v.pos, linetype = 2, color = col.lines) +
  374. geom_hline(yintercept = lines.h.pos, linetype = 2, color = col.lines) +
  375. geom_vline(xintercept = 0, color = col.lines) +
  376. geom_hline(yintercept = 0, color = col.lines) +
  378. scale_x_continuous(breaks = scales::pretty_breaks(12),
  379. expand = c(0,0)) +
  380. scale_y_continuous(breaks = scales::pretty_breaks(4),
  381. expand = c(0,0)) +
  382. scale_fill_stepsn(colors = col.pal,
  383. breaks = z.breaks,
  384. labels = z.labels,
  385. limits = c(min(z.breaks), max(z.breaks)),
  386. oob = scales::squish) +
  388. labs(x = "Test time (ms)", y = "Train time (ms)") +
  390. coord_equal(clip = "off") +
  392. theme(# Text size
  393. plot.title = element_text(size = 12, hjust = 0.5, face = "bold"),
  394. plot.subtitle = element_text(size = 9, hjust = 0.5),
  395. axis.title = element_text(size = 9),
  396. legend.title = element_text(size = 9),
  397. strip.text = element_text(size = 9),
  398. axis.text = element_text(size = 5),
  399. legend.text = element_text(size = 5.5),
  400. # Time ticks
  401. axis.ticks = element_line(color = "black"),
  402. axis.ticks.length = unit(0.1, "cm"),
  403. # Legend
  404. legend.position = legend.position,
  405. legend.box.margin = margin(0,0,0,0),
  406. legend.box.spacing = unit(0, "pt"),
  407. # no grid
  408. panel.grid = element_blank(),
  409. # don't clip facet titles
  410. strip.clip = "off")
  412. if (legend.position %in% c("bottom","top")) {
  413. g <- g +
  414. guides(fill = guide_colorbar(barwidth = 0.5*length(z.breaks-1), barheight = 0.5,
  415. ticks = F, title.vjust = 1))
  416. } else {
  417. g <- g +
  418. guides(fill = guide_colorbar(barheight = 0.5*length(z.breaks-1), barwidth = 0.5,
  419. ticks = F, title.vjust = 1))
  420. }
  422. # --- significant clusters
  423. if (contour) {
  424. g <- g +
  425. geom_contour(aes(z = signif), size = 0.1, color = col.contour)
  426. }
  428. return(g)
  429. }
  430. ```
  432. # **Replication & Reordering (empirical null)**
  434. ## Load data
  436. ### Sounds
  438. ```{r}
  439. #| cache: true
  441. tic()
  443. # Load subject-level data (used for clustering) and derive grand average (used for visualization)
  444. df.stim.mean <- tibble()
  445. df.stim.subj <- tibble()
  446. for (manip in c("","_reord")) {
  448. for (direction in c("rd_to_rd", "rd_to_mm", "rd_to_mp", "rd_to_or")) {
  449. tmp <- load_loop_subj(folder = "sounds",
  450. condition = paste0(direction, manip))
  452. df.stim.mean %<>% bind_rows(tmp[[2]] %>% mutate(manip = manip, direction = direction))
  453. df.stim.subj %<>% bind_rows(tibble(manip = manip, direction = direction, array = list(tmp[[1]])))
  454. }
  455. }
  457. df.stim.mean %<>% mutate(stim = "sounds")
  458. df.stim.subj %<>% mutate(stim = "sounds") %>%
  459. # --- center decoding accuracy on 0
  460. mutate(array = map(array, .f=\(x){x-0.25}))
  462. toc() # ~1s
  463. ```
  465. ### Omissions
  467. ```{r}
  468. #| cache: true
  470. tic()
  472. df.omit.mean <- tibble()
  473. df.omit.subj <- tibble()
  474. for (manip in c("","_reord")) {
  476. for (direction in c("rd_to_rd", "rd_to_mm", "rd_to_mp", "rd_to_or")) {
  478. tmp <- load_loop_subj(folder = "omissions",
  479. condition = paste0(direction, manip))
  481. df.omit.mean %<>% bind_rows(tmp[[2]] %>% mutate(manip = manip, direction = direction))
  482. df.omit.subj %<>% bind_rows(tibble(manip = manip, direction = direction, array = list(tmp[[1]])))
  483. }
  484. }
  486. df.omit.mean %<>% mutate(stim = "omissions") %>% filter(abs(t_test) <= 670)
  487. df.omit.subj %<>% mutate(stim = "omissions") %>%
  488. # --- center decoding accuracy on 0
  489. mutate(array = map(array, .f=\(x){x-0.25})) %>%
  490. # --- crop subject-level data for cluster analysis
  491. mutate(array = map(array, .f=\(x){x[71:104,4:138,]}))
  493. toc() # ~ 2s
  494. ```
  496. ## Correlations with entropy
  498. ```{r}
  499. #| cache: true
  500. #| warning: false
  502. ## Prepare input data
  503. df.tmp <- bind_rows(df.stim.subj, df.omit.subj) %>%
  504. calculate_diff(f.map = map2, values_to = "array") %>%
  505. mutate(regularity = str_remove(direction, "rd_to_"),
  506. regularity = factor(regularity, levels = c("rd","mm","mp","or"))) %>%
  507. select(-direction) %>%
  508. arrange(as.numeric(regularity))
  509. df.tmp
  512. tic()
  513. df.corr.mean <- tibble()
  514. df.corr.subj <- tibble()
  516. ## Loop through condition "manip" and stimulus type "stim
  517. for (mm in c("original","reordered","diff")) {
  519. for (ss in c("sounds", "omissions")) { #
  520. cat("\n\nProcessing ", str_to_upper(mm), " / ", str_to_upper(ss), "...\n", sep="")
  522. # --- Loop through subjects to derive a 3D array of participant-level time-generalized correlations with entropy, from a 3D array of participant-level accuracy matrix
  523. for (idx.subj in 1:n.subj) {
  525. # --- bind time-generalization matrices for the 3 non-random conditions into a single 3D array
  526. # (add accuracy data from the random condition if we are processing sounds)
  527. if (ss == "sounds") {
  528. idx.lines = 1:4
  529. } else if (ss == "omissions") {
  530. idx.lines = 2:4
  531. }
  532. arr.tmp <- correlations_extract3D(data = filter(df.tmp, manip == mm, stim == ss),
  533. idx.lines = idx.lines, idx.subj)
  535. # --- calculate correlation with entropy level coded as 0,1,2(,3), obtaining a 2D array
  536. if (ss == "sounds") {
  537. arr.corr.tmp <- apply(arr.tmp, 1:2, correlations_vecs, vec2 = 0:3)
  538. } else if (ss == "omissions") {
  539. arr.corr.tmp <- apply(arr.tmp, 1:2, correlations_vecs, vec2 = 1:3)
  540. }
  542. # --- append to already processed subjects in a 3D array
  543. if (idx.subj == 1) {
  544. arr.corr <- arr.corr.tmp
  545. } else {
  546. arr.corr %<>% abind(arr.corr.tmp, along = 3)
  547. }
  548. }
  550. # Postprocessing
  551. arr.corr %<>%
  552. correlations_fixval %>% # replace NAs by 0 and extreme values by +/- .99
  553. atanh # apply Fisher's transformation
  555. # --- store participant-level results in a dataframe
  556. df.corr.subj %<>% bind_rows(tibble(manip = mm, stim = ss, array = list(arr.corr)))
  558. # --- calculate groupe average and store in a dataframe
  559. df.corr.mean %<>% bind_rows(
  560. arr.corr %>%
  561. apply(c(1,2), mean) %>% # average across participants
  562. tanh %>% # convert back to correlation coefficients
  563. matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "r") %>%
  564. mutate(manip = mm, stim = ss)
  565. )
  566. }
  567. }
  569. toc() # ~20 s
  570. ```
  572. ## Clusters stats
  574. ### Accuracy
  576. ```{r}
  577. #| cache: true
  578. #| lightbox:
  579. #| group: global
  581. tic()
  583. df.tmp.subj <- bind_rows(df.stim.subj, df.omit.subj)
  584. df.tmp.mean <- bind_rows(df.stim.mean, df.omit.mean)
  586. # Cluster analysis
  587. df.clusters.0 <- df.tmp.subj %>%
  588. calculate_diff(f.map = map2, values_to = "array") %>%
  589. cluster_loop_cond(n = n_permutations) %>%
  590. # --- convert indices to time
  591. mutate(t_test = (dt*(t_test-1)-time_offset)) %>%
  592. mutate(t_train = (dt*(t_train-1)))
  595. # Quick check of clusters with a mask plot
  596. df.clusters.0 %>%
  597. mutate(direction = (direction %>% str_replace("rd_to_", "") %>% str_to_upper)) %>%
  598. ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  599. facet_grid(direction ~ manip) +
  600. geom_tile() +
  601. coord_equal()
  604. # Aggregate all data and clean
  605. df.plot.0 <- df.tmp.mean %>%
  606. calculate_diff(f.map = map2_dbl, values_to = "accuracy") %>%
  607. # --- append cluster significance
  608. left_join(df.clusters.0) %>%
  609. mutate(manip = (recode_factor(manip, original = "ORIGINAL", reordered = "REORDERED",
  610. diff = "ORIGINAL – REORDERED"))) %>%
  611. mutate(direction = toupper(sub("rd_to_", "", direction)),
  612. direction = factor(direction, levels = rev(c("RD","MM","MP","OR"))))
  614. toc() # ~3.5 minutes
  615. ```
  617. ### Correlations
  619. ```{r}
  620. #| cache: true
  621. #| lightbox:
  622. #| group: global
  624. tic()
  626. # Cluster analysis
  627. df.clusters.corr <- df.corr.subj %>%
  628. cluster_loop_cond(n = n_permutations)
  630. # Quick check of clusters with a mask plot
  631. df.clusters.corr %>%
  632. ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  633. facet_grid( ~ manip) +
  634. geom_tile() +
  635. coord_equal()
  637. # --- append cluster significance
  638. df.plot.0.corr <- df.corr.mean %>%
  639. left_join(df.clusters.corr) %>%
  640. mutate(manip = (recode_factor(manip, original = "ORIGINAL", reordered = "REORDERED",
  641. diff = "ORIGINAL – REORDERED"))) %>%
  642. # --- convert indices to time
  643. mutate(t_test = (dt*(t_test-1)-time_offset)) %>%
  644. mutate(t_train = (dt*(t_train-1)))
  646. df.plot.0.corr %>% summary
  648. toc() # ~ 1 min
  649. ```
  651. # **Theoretical null**
  653. ## Load confusion matrices & calculate
  655. ```{r}
  656. #| cache: true
  658. tic()
  660. null_accuracy <- function(mat.confusion, mat.transition) {
  661. # Calculate accuracy under the null hypothesis from the transition and confusion matrices mat.T and mat.C
  662. # The 3/4 factor is used to compensate the unbiased estimator used in cov() by default (see ?cov)
  663. return(0.25 + sum(diag(3/4*cov(t(mat.transition), t(mat.confusion)))))
  664. }
  666. df.maths.subj <- tibble()
  667. df.maths.mean <- tibble()
  668. for (entropy in c("RD","MM","MP","OR")) {
  670. array.scores <- array(0.25, dim = c(33,136,n.subj))
  672. i <- 1
  673. for (s in list.subj) {
  675. # --- load the confusion matrix
  676. tmp.file <- os$path$join(path.root, s, "rd_to_rd_confmats.npz") %>%
  677. np$load()
  679. # --- crop around S0 (decoded stimulus)
  680. mat.C <- tmp.file$f[["arr_0"]][,71:103,71:104,,] %>%
  681. apply(c(2,3,4,5), mean)
  683. # NOTE: because the transition matrix is the same for all subjects, theoretical accuracy can be calculated from the group average of the confusion matrix (faster)
  684. # --- calculate theoretical accuracy
  685. # --- for S-2
  686. array.scores[1:33,1:34,i] <- mat.C %>%
  687. apply(c(1,2), null_accuracy, mat.transition = mat.T[[entropy]]%^%2)
  688. # --- for S-1
  689. array.scores[1:33,35:68,i] <- mat.C %>%
  690. apply(c(1,2), null_accuracy, mat.transition = mat.T[[entropy]]%^%1)
  691. # --- for S
  692. array.scores[1:33,69:102,i] <- mat.C %>%
  693. apply(c(1,2), \(x)mean(diag(x)))
  694. # --- for S+1
  695. array.scores[1:33,103:136,i] <- mat.C %>%
  696. apply(c(1,2), null_accuracy, mat.transition = t(mat.T[[entropy]]))
  697. i <- i+1
  698. }
  700. array.scores <- array.scores[,2:136,]
  702. # Append to other entropy levels in dataframe
  703. # --- subject-level data
  704. df.maths.subj %<>% bind_rows(tibble(array = list(array.scores), entropy = entropy))
  705. # --- group-averaged data
  706. df.maths.mean %<>% bind_rows(array.scores %>%
  707. apply(c(1,2), mean) %>%
  708. matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "accuracy") %>%
  709. mutate(entropy = entropy)
  710. )
  711. }
  713. toc() # ~ 1min
  714. ```
  716. ## Append original data & calculate diff
  718. ```{r}
  720. df.maths.all <-
  721. # --- aggregate original, theoretical and difference between the 2
  722. bind_rows(df.stim.subj %>%
  723. # --- select original analysis
  724. filter(manip == "") %>% mutate(manip = "original") %>%
  725. # --- crop to the same time-window than for theoretical analysis
  726. mutate(entropy = toupper(str_replace(direction,"rd_to_",""))),
  727. df.maths.subj %>%
  728. mutate(manip = "maths") %>%
  729. mutate(array = map(array, .f=\(x){x-0.25}))) %>%
  730. select(manip, entropy, array) %>%
  731. pivot_wider(names_from = manip, values_from = "array") %>%
  732. mutate(diff = map2(original, maths, .f=\(x,y){x-y})) %>%
  733. pivot_longer(cols = c(original, maths, diff), names_to = "manip", values_to = "array")
  735. ```
  737. ## Correlations with entropy
  739. ```{r}
  740. #| cache: true
  741. #| warning: false
  743. tic()
  745. df.corr.maths.mean <- tibble()
  746. df.corr.maths.subj <- tibble()
  748. for (mm in c("original","maths","diff")) {
  749. # --- Loop through subjects to derive a 3D array of participant-level time-generalized correlations with entropy, from a 3D array of participant-level accuracy matrix
  750. for (idx.subj in 1:n.subj) {
  752. # --- bind time-generalization matrices of the 4 conditions into a single 3D array
  753. arr.tmp <- correlations_extract3D(data = df.maths.all %>% filter(manip == mm),
  754. idx.lines = 1:4, idx.subj)
  756. # --- calculate correlation with entropy level coded as 0,1,2,3, obtaining a 2D array
  757. arr.corr.tmp <- apply(arr.tmp, 1:2, correlations_vecs, vec2 = 0:3)
  759. # --- append to already processed subjects in a 3D array
  760. if (idx.subj == 1) {
  761. arr.corr <- arr.corr.tmp
  762. } else {
  763. arr.corr %<>% abind(arr.corr.tmp, along = 3)
  764. }
  765. }
  767. # Postprocessing
  768. arr.corr %<>%
  769. correlations_fixval %>% # replace NAs by 0 and extreme values by +/- .99
  770. atanh # apply Fisher's transformation
  772. # --- store participant-level results in a dataframe of arrays
  773. df.corr.maths.subj %<>% bind_rows(tibble(manip = mm, array = list(arr.corr)))
  775. # --- calculate group average and store in a dataframe
  776. df.corr.maths.mean %<>% bind_rows(
  777. arr.corr %>%
  778. apply(c(1,2), mean) %>% # average across participants
  779. tanh %>% # convert back to correlation coefficients
  780. matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "r") %>%
  781. mutate(manip = mm)
  782. )
  783. }
  785. toc() # ~ 20s
  786. ```
  788. ## Clusters stats
  790. ### Accuracy
  792. ```{r}
  793. #| cache: true
  794. #| lightbox:
  795. #| group: global
  797. tic()
  799. df.clusters.maths <- df.maths.all %>%
  800. cluster_loop_cond(n = n_permutations) %>%
  801. # --- convert indices to time
  802. mutate(t_test = (dt*(t_test-1)-time_offset)) %>%
  803. mutate(t_train = (dt*(t_train-1)))
  805. toc() # ~ 2 min
  807. # Quick check of clusters with a mask plot
  808. df.clusters.maths %>%
  809. ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  810. facet_grid(entropy ~ manip) +
  811. geom_tile() +
  812. coord_equal()
  814. ```
  816. ### Correlations
  818. ```{r}
  819. #| cache: true
  820. #| lightbox:
  821. #| group: global
  823. tic()
  824. df.clusters.corr.maths <- df.corr.maths.subj %>%
  825. cluster_loop_cond(n = n_permutations)
  826. toc() # ~ 20s
  828. # Quick check of clusters with a mask plot
  829. df.clusters.corr.maths %>%
  830. ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  831. facet_grid( ~ manip) +
  832. geom_tile() +
  833. coord_equal()
  835. # --- append cluster significance
  836. df.plot.maths.corr <- df.corr.maths.mean %>%
  837. left_join(df.clusters.corr.maths) %>%
  838. mutate(manip = (recode_factor(manip, original = "ORIGINAL", maths = "THEORETICAL NULL",
  839. diff = "ORIGINAL – THEORETICAL NULL"))) %>%
  840. # --- convert indices to time
  841. mutate(t_test = (dt*(t_test-1)-time_offset)) %>%
  842. mutate(t_train = (dt*(t_train-1)))
  844. df.plot.maths.corr %>% summary
  845. ```
  847. # **Prediction of the most likely ("simple prediction")**
  849. ## Load data
  851. ```{r}
  852. #| cache: true
  854. tic()
  856. # Load subject-level data and derive grand average
  857. df.stim.ml.mean <- tibble()
  858. df.stim.ml.subj <- tibble()
  859. for (manip in c("","_reord")) {
  861. for (direction in c("rd_to_mm", "rd_to_mp", "rd_to_or")) {
  863. tmp <- load_loop_subj(folder = "sounds",
  864. condition = paste0(direction, manip, "_sp"))
  866. df.stim.ml.mean %<>% bind_rows(tmp[[2]] %>% mutate(manip = manip, direction = direction))
  867. df.stim.ml.subj %<>% bind_rows(tibble(manip = manip, direction = direction, array = list(tmp[[1]])))
  868. }
  869. }
  871. df.stim.ml.mean %<>% mutate(stim = "sounds")
  872. df.stim.ml.subj %<>% mutate(stim = "sounds") %>%
  873. # --- center decoding accuracy on 0
  874. mutate(array = map(array, .f=\(x){x-0.25}))
  876. toc() # ~ 1s
  878. df.stim.ml.mean %>% mutate_if(is.character, as.factor) %>% summary
  880. ```
  882. ## Correlations with entropy
  884. ```{r}
  885. #| cache: true
  887. tic()
  889. ## Prepare input data
  890. df.tmp <- df.stim.ml.subj %>%
  891. calculate_diff(f.map = map2, values_to = "array") %>%
  892. mutate(entropy = str_remove(direction, "rd_to_"),
  893. entropy = factor(entropy, levels = c("rd","mm","mp","or"))) %>%
  894. select(-direction) %>%
  895. arrange(as.numeric(entropy))
  898. ## Calculate correlations
  899. df.corr.ml.mean <- tibble()
  900. df.corr.ml.subj <- tibble()
  902. # --- Loop through condition "manip" and stimulus type "stim"
  903. for (mm in c("original","reordered","diff")) {
  905. for (ss in c("sounds")) {
  906. cat("\n\nProcessing ", str_to_upper(mm), " / ", str_to_upper(ss), "...\n", sep="")
  909. # --- Loop through subjects to derive a 3D array of participant-level time-generalized correlations with entropy, from a 3D array of participant-level accuracy matrix
  910. for (idx.subj in 1:n.subj) {
  912. # --- bind time-generalization matrices for the 3 non-random conditions into a single 3D array
  913. arr.tmp <- correlations_extract3D(data = filter(df.tmp, manip == mm, stim == ss),
  914. idx.lines = 1:3, idx.subj)
  916. # --- calculate correlation with entropy level coded as 0,1,2, obtaining a 2D array
  917. arr.corr.tmp <- apply(arr.tmp, 1:2, correlations_vecs, vec2 = 0:2)
  919. # --- append to already processed subjects in a 3D array
  920. if (idx.subj == 1) {
  921. arr.corr <- arr.corr.tmp
  922. } else {
  923. arr.corr %<>% abind(arr.corr.tmp, along = 3)
  924. }
  925. }
  927. # Post-processing
  928. arr.corr <- arr.corr %>%
  929. correlations_fixval %>% # replace NAs by 0 and extreme values by +/- .99
  930. atanh # apply Fisher's transformation
  932. # --- store participant-level results in a dataframe
  933. df.corr.ml.subj %<>% bind_rows(tibble(manip = mm, stim = ss, array = list(arr.corr)))
  935. # --- calculate groupe average and store in a dataframe
  936. df.corr.ml.mean %<>% bind_rows(
  937. arr.corr %>%
  938. apply(c(1,2), mean) %>% # average across participants
  939. tanh %>% # convert back to correlation coefficients
  940. matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "r") %>%
  941. mutate(manip = mm, stim = ss)
  942. )
  943. }
  944. }
  946. toc() # ~ 15s
  947. ```
  949. ## Cluster stats
  951. ### Accuracy
  953. ```{r}
  954. #| cache: true
  955. #| lightbox:
  956. #| group: global
  958. tic()
  960. ## Cluster analysis with subject-level data
  961. # --- calculate original-reordered difference
  962. tmp.diff <- df.stim.ml.subj %>%
  963. calculate_diff(f.map = map2, values_to = "array")
  965. # --- combine & run cluster analysis
  966. df.clusters.ml <- tmp.diff %>% #bind_rows(tmp.diff, tmp.ave) %>%
  967. cluster_loop_cond(n = n_permutations) %>%
  968. # --- convert indices to time
  969. mutate(t_test = (dt*(t_test-1)-time_offset)) %>%
  970. mutate(t_train = (dt*(t_train-1)))
  972. # Quick check of clusters with a mask plot
  973. df.clusters.ml %>%
  974. mutate(direction = (direction %>% str_replace("rd_to_", "") %>% str_to_upper)) %>%
  975. ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  976. facet_grid(direction ~ manip) +
  977. geom_tile() +
  978. coord_equal()
  980. ## Grand average data
  981. # --- calculate original-reordered difference
  982. tmp.diff <- df.stim.ml.mean %>%
  983. calculate_diff(f.map = map2_dbl, values_to = "accuracy")
  984. # --- calculate average of difference across entropy levels
  985. tmp.ave <- tmp.diff %>%
  986. filter(manip == "diff") %>%
  987. group_by(t_train, t_test, manip, stim) %>%
  988. summarise(accuracy = mean(accuracy, na.rm=T)) %>%
  989. mutate(direction = "average of\n(MM,MP,OR)")
  990. # --- combine
  991. df.plot.ml <- tmp.diff #bind_rows(tmp.diff, tmp.ave)
  993. ## Aggregate all data and clean
  994. df.plot.ml %<>%
  995. # --- append cluster significance
  996. left_join(df.clusters.ml) %>%
  997. mutate(manip = (recode_factor(manip, original = "ORIGINAL", reordered = "REORDERED",
  998. diff = "ORIGINAL – REORDERED"))) %>%
  999. mutate(direction = toupper(sub("rd_to_", "", direction))) %>%
  1000. mutate(direction = factor(direction, levels = rev(c("AVERAGE OF\n(MM,MP,OR)","MM","MP","OR"))))
  1002. toc() # ~1.5min
  1003. ```
  1005. ### Correlations
  1007. ```{r}
  1008. #| cache: true
  1009. #| lightbox:
  1010. #| group: global
  1012. tic()
  1014. df.clusters.corr.ml <- df.corr.ml.subj %>%
  1015. cluster_loop_cond(n = n_permutations)
  1017. toc() # ~ 30s
  1020. # Quick check of clusters with a mask plot
  1021. df.clusters.corr.ml %>%
  1022. ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  1023. facet_grid( ~ manip) +
  1024. geom_tile() +
  1025. coord_equal()
  1028. # --- append cluster significance
  1029. df.plot.ml.corr <- df.corr.ml.mean %>%
  1030. filter(between(t_test,35,68)) %>%
  1031. left_join(df.clusters.corr.ml) %>%
  1032. mutate(manip = (recode_factor(manip, original = "ORIGINAL", reordered = "REORDERED",
  1033. diff = "ORIGINAL – REORDERED"))) %>%
  1034. # --- convert indices to time
  1035. mutate(t_test = (dt*(t_test-1)-time_offset)) %>%
  1036. mutate(t_train = (dt*(t_train-1)))
  1039. df.plot.ml.corr %>% summary
  1040. ```
  1042. # **Figures**
  1044. ## Methodological plots
  1046. ### Transition matrices (fig. 1a)
  1048. ```{r}
  1049. #| output: false
  1051. # Convert & combine transition matrices into a dataframe
  1052. df.plot <- tibble()
  1053. for (condition in c("RD","MM","MP","OR")) {
  1054. df.plot %<>% bind_rows(
  1055. mat.T[[condition]] %>%
  1056. matrix_to_df2(row.name = "from", col.name = "to", val.name = "p") %>%
  1057. mutate(condition = condition)
  1058. )
  1059. }
  1060. df.plot %<>%
  1061. mutate(condition = factor(condition, names(mat.T)))
  1063. df.plot.pitch <- bind_rows(tibble(x = 0.1, y = 1:4, pitch = 1:4),
  1064. tibble(y = 0.1, x = 1:4, pitch = 1:4))
  1065. # Plot
  1066. df.plot %>%
  1067. ggplot(aes(x = to, y = from)) +
  1068. facet_wrap(~ condition, nrow = 1, strip.position = "bottom") +
  1069. # --- heatmap
  1070. geom_tile(aes(fill = p), color = "black") +
  1071. # --- display percentages
  1072. geom_text(data = . %>% filter(p<=.25),
  1073. aes(label = paste0(100*p,"%")),
  1074. color = "black", size = 2) +
  1075. geom_text(data = . %>% filter(p>.25),
  1076. aes(label = paste0(100*p,"%")),
  1077. color = "white", size = 2) +
  1079. # --- display color-coded pitch in axes
  1080. geom_point(data = df.plot.pitch,
  1081. aes(x, y, color = pitch), size = 4, show.legend = FALSE) +
  1082. geom_text(data = df.plot.pitch,
  1083. aes(x, y, label = LETTERS[pitch]), color = "white", size = 2, show.legend = FALSE) +
  1085. scale_color_viridis_c(end = 0.8) +
  1086. scale_x_continuous(name = "to\n", position = "top") +
  1087. scale_y_reverse(name = "from\n") +
  1088. scale_fill_gradient(name = "transition probability",
  1089. low = "white", high = "black",
  1090. labels = ~str_replace(., "^0.","."),
  1091. limits = c(0,1)) +
  1092. guides(fill = guide_colorbar(barheight = 0.5,
  1093. ticks = F, title.vjust = 1)) +
  1094. coord_equal(xlim = c(0.5,4.5), ylim = rev(c(0.5,4.5)),
  1095. expand = F, clip = "off") +
  1096. theme(legend.position = "bottom",
  1097. legend.box.margin = margin(t = -20),
  1098. # text size
  1099. axis.title = element_text(size = 9),
  1100. legend.title = element_text(size = 9),
  1101. strip.text = element_text(size = 12),
  1102. # remove axis text
  1103. axis.text = element_blank(),
  1104. # Positions & margins
  1105. axis.title.x = element_text(hjust = 0.08),
  1106. panel.spacing.x = unit(0.8,"cm")) -> g.transitions
  1108. g.transitions
  1109. ```
  1111. ### Reordering (fig. 1b)
  1113. ```{r}
  1114. #| output: false
  1116. # Generate Markov sequences
  1117. n.samples <- 12
  1118. stims <- list(RD = c(), OR = c())
  1119. set.seed(128764) #20476, 99713
  1120. for (entropy in c("RD","OR")) {
  1121. stims[[entropy]] <- vector("numeric", n.samples)
  1122. stims[[entropy]][1] <- sample(1:4, 1)
  1123. for (i in 2:n.samples) {
  1124. stims[[entropy]][i] <- sample(1:4, 1, prob = mat.T[[entropy]][stims[[entropy]][i-1],])
  1125. }
  1126. }
  1128. # Add global sequence indices, class-specific indices and offset according to pitch
  1129. df.plot <- stims %>% as.tibble() %>%
  1130. mutate(idx = as.double(1:n.samples)) %>%
  1131. pivot_longer(-idx, names_to = "entropy", values_to = "class") %>%
  1132. mutate(entropy = factor(entropy, levels = c("OR","RD"))) %>%
  1133. group_by(entropy,class) %>% mutate(rank = rank(idx)) %>%
  1134. mutate(y = as.double(entropy) + class/16-10/64)
  1136. pitch.y <- sort(unique(df.plot$y))
  1138. # Plot
  1139. df.plot %>%
  1140. ggplot(aes(x = idx, y = y, color = class)) +
  1142. # --- pitch lines
  1143. geom_hline(yintercept = pitch.y, color = "grey90") +
  1144. # --- pitch labels
  1145. annotate(geom = "text", label = "pitch",
  1146. x = 12.5, y = c(mean(head(pitch.y,4)), mean(tail(pitch.y,4))),
  1147. hjust = 0.5, vjust = 2.5, size = 2.5, color = "grey50", angle = 90) +
  1149. annotate(geom = "text", label = "low",
  1150. x = 12.5, y = c(min(head(pitch.y,4)), min(tail(pitch.y,4))),
  1151. hjust = 0.74, vjust = 1.2, size = 2, color = "grey50", angle = 90) +
  1153. annotate(geom = "text", label = "high",
  1154. x = 12.5, y = c(max(head(pitch.y,4)), max(tail(pitch.y,4))),
  1155. hjust = 0.25, vjust = 1.2, size = 2, color = "grey50", angle = 90) +
  1157. # --- reordering trajectories
  1158. ggh4x::geom_pointpath(data = df.plot %>% arrange(desc(entropy)),
  1159. aes(x = idx, y = y, group = interaction(class,rank)),
  1160. mult = 0.4, color = "black", linewidth = 0.3,
  1161. arrow = arrow(type = "closed", angle = 25, length=unit(.15, 'cm'))) +
  1162. # --- stimuli
  1163. geom_point(size = 4, show.legend = F) +
  1164. # --- stimuli ranks
  1165. geom_text(aes(label = rank), size = 2, color = "white") +
  1167. scale_x_continuous(name = "trial number", breaks = 1:n.samples) +
  1168. scale_y_continuous(breaks = c(1,2), labels = c("OR","RD"), position = "left") +
  1169. scale_color_viridis_c(end = 0.8) +
  1170. coord_cartesian(clip = "off") +
  1172. theme(axis.title.y = element_blank(),
  1173. axis.text.y = element_text(size = 10, color = "black", margin = margin(r=-2)),
  1174. axis.title.x = element_text(size = 9, color = "black"),
  1175. panel.grid = element_blank(),
  1176. panel.grid.major.x = element_line(color = "grey80", linewidth = 0.2)) -> g.reordering
  1178. g.reordering
  1179. ```
  1181. ### Prediction types (fig. 2a)
  1183. Presented vs. most likely
  1185. ```{r}
  1186. #| output: false
  1188. # Generate Markov sequences
  1189. n.samples <- 6
  1190. stims <- mat.T
  1191. set.seed(12871)
  1192. for (entropy in names(stims)) {
  1193. stims[[entropy]] <- vector("numeric", n.samples)
  1194. stims[[entropy]][1] <- 4#sample(1:4, 1)
  1195. for (i in 2:n.samples) {
  1196. stims[[entropy]][i] <- sample(1:4, 1, prob = mat.T[[entropy]][stims[[entropy]][i-1],])
  1197. }
  1198. }
  1200. # Add sequence index & predictions
  1201. df.plot <- stims %>% as.tibble() %>%
  1202. mutate(idx = as.double(1:n.samples)) %>%
  1203. pivot_longer(-idx, names_to = "entropy", values_to = "class") %>%
  1204. mutate(entropy = factor(entropy, levels = rev(names(stims)))) %>%
  1205. group_by(entropy) %>% mutate(actual = lead(class),
  1206. mostlikely = ifelse(class==1, 4, class-1)) %>%
  1207. pivot_longer(c(actual,mostlikely), names_to = "status", values_to = "class.pred") %>%
  1208. filter(entropy != "RD")
  1210. # Plot
  1211. df.plot %>%
  1212. ggplot(aes(x = idx, y = class)) +
  1213. facet_wrap(~ entropy, ncol = 1, strip.position = "right") +
  1215. geom_segment(data = df.plot %>% filter(status == "mostlikely", idx != n.samples),
  1216. aes(x = idx, xend = idx+1, y = class, yend = class.pred, linetype = as.factor(sign(class.pred-class))),
  1217. #arrow = arrow(type = "closed", length = unit(6,"pt")),
  1218. color = "grey50", show.legend = F) +
  1220. geom_point(aes(fill = class), size = 3, shape = 21, color = "white", show.legend = T) +
  1222. geom_point(data = df.plot %>% filter(status == "mostlikely", idx != n.samples),
  1223. aes(x = idx+1, y = class.pred, color = class.pred),
  1224. size = 2.5, shape = 8) +
  1226. scale_x_continuous(name = "trial number", breaks = 1:n.samples) +
  1227. scale_y_continuous(name = "") +
  1228. scale_color_viridis_c(end = 0.8, guide = "legend", name = "", labels = c("","","","most likely")) +
  1229. scale_fill_viridis_c(end = 0.8, guide = "legend", name = "", labels = c("","","","actually presented")) +
  1231. guides(linetype = F) +
  1232. coord_cartesian(clip = "off") +
  1233. expand_limits(x = c(0.9,6.1), y = c(0.3,4.7)) +
  1235. guides(#color = guide_legend(order = 1, override.aes = list(color = "grey50")),
  1236. color = guide_legend(order = 1),
  1237. fill = guide_legend(order = 2)) +
  1238. theme(plot.caption = element_text(size = 6),
  1239. axis.text.y = element_blank(),
  1240. axis.title.x = element_text(size = 9, color = "black"),
  1241. strip.text = element_text(size = 12),
  1242. panel.spacing = unit(0.2,"cm"),
  1243. panel.grid = element_blank(),
  1244. panel.border = element_rect(fill=NA),
  1245. legend.position = "top",
  1246. legend.box = "vertical",
  1247. legend.box.just = "left",
  1248. legend.margin = margin(b = -10),
  1249. legend.text = element_text(margin = margin(r = -12, l=0, unit = "pt"))
  1250. ) -> g.mostlikely
  1252. g.mostlikely
  1253. ```
  1255. ## Figure 1 & Supplementary Figure 2
  1257. ### Initialize
  1259. ```{r}
  1260. stims <- list(RD = c(), OR = c())
  1261. g.base <- list(sound = c(), omission = c())
  1262. g.diff <- list(sound = c(), omission = c())
  1264. # Create ad hoc dataframes for onset' labels
  1265. df.labels.onset <- list(
  1266. sounds = tibble(label = "0 = sound onset",
  1267. manip = factor(c("REPLICATION\n", "EMPIRICAL NULL\n",
  1268. "REPLICATION – EMPIRICAL NULL\n")),
  1269. direction = factor("OR", levels = names(stims))),
  1270. omissions = tibble(label = "0 = omission onset",
  1271. manip = factor(c("REPLICATION\n", "EMPIRICAL NULL\n",
  1272. "REPLICATION – EMPIRICAL NULL\n")),
  1273. direction = factor("OR", levels = names(stims)))
  1274. )
  1276. # Function to add onset' labels to a plot
  1277. add_onset <- function(g, df.labels) {
  1278. g <- g + geom_text(data = df.labels %>% filter(manip %in% unique(g$data$manip)),
  1279. aes(label = label, x = 0, y = 350),
  1280. hjust = 0, vjust = -0.5, size = 3, fontface = "italic")
  1281. return(g)
  1282. }
  1283. ```
  1285. ### Original & empirical null: accuracy plots
  1287. ```{r}
  1288. #| warning: false
  1289. #| lightbox:
  1290. #| group: global
  1292. for (s in c("sounds","omissions")) { #
  1293. # ORIGINAL & REORDERED
  1294. df.plot.0 %>%
  1295. filter(stim == s, manip != "ORIGINAL – REORDERED") %>%
  1296. mutate(manip = recode_factor(manip,
  1297. ORIGINAL = "REPLICATION\n",
  1298. REORDERED = "EMPIRICAL NULL\n")) %>%
  1299. plot_base(col.pal = col.pal.base,
  1300. z.breaks = acc.breaks.base,
  1301. z.labels = acc.labels.base) +
  1302. facet_grid(direction ~ manip) -> g
  1304. g.base[["acc"]][[s]] <- g %>% add_onset(df.labels.onset[[s]])
  1306. # DIFFERENCE
  1307. g.diff[["acc"]][[s]] <- df.plot.0 %>%
  1308. filter(stim == s, manip == "ORIGINAL – REORDERED") %>%
  1309. mutate(manip = recode_factor(manip,
  1310. `ORIGINAL – REORDERED` = "REPLICATION – EMPIRICAL NULL\n")) %>%
  1311. plot_base(col.pal = col.pal.diff,
  1312. z.breaks = acc.breaks.diff,
  1313. z.labels = acc.labels.diff) +
  1314. facet_grid(direction ~ manip) +
  1315. theme(axis.title.y = element_blank()) -> g
  1317. g.diff[["acc"]][[s]] <- g %>% add_onset(df.labels.onset[[s]])
  1318. }
  1320. # Assemble Figure 1c
  1321. g.fig1.1.acc <- (
  1322. g.base[["acc"]][["sounds"]] +
  1323. (g.diff[["acc"]][["sounds"]] + plot_layout(tag_level = "new")) +
  1324. plot_layout(widths = c(2, 1)) &
  1325. theme(strip.text.x = element_text(face = "bold"))
  1326. )
  1328. g.fig1.1.acc
  1330. # Assemble Supplementary Figure 1a
  1331. g.supfig2.acc <- (
  1332. g.base[["acc"]][["omissions"]] +
  1333. (g.diff[["acc"]][["omissions"]] + plot_layout(tag_level = "new")) +
  1334. plot_layout(widths = c(2, 1)) &
  1335. theme(strip.text.x = element_text(face = "bold"))
  1336. )
  1338. g.supfig2.acc
  1340. ```
  1342. ### Original & empirical null: correlations plots
  1344. ```{r}
  1345. #| warning: false
  1346. #| lightbox:
  1347. #| group: global
  1349. for (s in c("sounds", "omissions")) {
  1350. # ORIGINAL & REORDERED
  1351. df.plot.0.corr %>%
  1352. filter(stim == s, manip != "ORIGINAL – REORDERED") %>%
  1353. mutate(manip = recode_factor(manip,
  1354. ORIGINAL = "REPLICATION\n", # (original dataset)\n
  1355. REORDERED = "EMPIRICAL NULL\n")) %>% # (reordered RD)\n
  1356. mutate(direction = "") %>%
  1357. plot_base(col.pal = col.pal.corr,
  1358. z = "r",
  1359. z.breaks = corr.breaks,
  1360. z.labels = corr.labels) +
  1361. labs(fill = "Correlation") +
  1362. facet_grid(direction ~ manip) -> g
  1364. g.base[["corr"]][[s]] <- g #%>% add_onset(df.labels.onset[[s]] %>% mutate(direction = ""))
  1366. # DIFFERENCE
  1367. df.plot.0.corr %>%
  1368. filter(stim == s, manip == "ORIGINAL – REORDERED") %>%
  1369. mutate(manip = recode_factor(manip,
  1370. `ORIGINAL – REORDERED` = "REPLICATION – EMPIRICAL NULL\n")) %>%
  1371. mutate(direction = "") %>%
  1372. plot_base(col.pal = col.pal.corr,
  1373. z = "r",
  1374. z.breaks = corr.breaks,
  1375. z.labels = corr.labels) +
  1376. facet_grid(. ~ manip) +
  1377. labs(fill = "Correlation") +
  1378. theme(axis.title.y = element_blank()) -> g
  1380. g.diff[["corr"]][[s]] <- g #%>% add_onset(df.labels.onset[[s]] %>% mutate(direction = ""))
  1381. }
  1383. # Assemble Figure 1d
  1384. g.fig1.1.corr <- (
  1385. g.base[["corr"]][["sounds"]] +
  1386. (g.diff[["corr"]][["sounds"]] + plot_layout(tag_level = "new")) +
  1387. plot_layout(widths = c(2, 1), guides = "collect") &
  1388. theme(strip.text = element_blank(),
  1389. legend.position = "bottom",
  1390. legend.margin = margin(t = -20, l = 60))
  1391. )
  1393. g.fig1.1.corr
  1395. # Assemble Supplementary Figure 1a
  1396. g.supfig2.corr <- (
  1397. g.base[["corr"]][["omissions"]] +
  1398. (g.diff[["corr"]][["omissions"]] + plot_layout(tag_level = "new")) +
  1399. plot_layout(widths = c(2, 1)) &
  1400. theme(strip.text.x = element_text(face = "bold"))
  1401. )
  1403. g.supfig2.corr
  1404. ```
  1406. ### Theoretical null: accuracy plots
  1408. ```{r}
  1410. ## Aggregate all data and clean
  1412. df.plot.maths <- bind_rows(
  1413. # --- make datasets compatible
  1414. df.maths.mean %>% mutate(manip = "maths") %>%
  1415. mutate(t_train = dt*(t_train-1) - 0,
  1416. t_test = dt*(t_test-1) - time_offset),
  1417. df.stim.mean %>%
  1418. filter(manip == "") %>%
  1419. mutate(manip = "original") %>%
  1420. mutate(entropy = toupper(sub("rd_to_", "", direction))) %>%
  1421. # mutate(across(starts_with("t_"), ~(dt*(.-1)-700))) %>%
  1422. select(-stim, -direction)
  1423. ) %>%
  1424. # --- calculate difference
  1425. pivot_wider(names_from = manip, values_from = "accuracy") %>%
  1426. mutate(diff = original - maths) %>%
  1427. pivot_longer(cols = c(original, maths, diff),
  1428. names_to = "manip", values_to = "accuracy") %>%
  1429. # --- crop time window
  1430. filter(t_train>=0, (t_train-333) < dt) %>%
  1431. # --- aggregate results of cluster analysis
  1432. left_join(df.clusters.maths) %>%
  1433. # --- clean labels
  1434. mutate(manip = (recode_factor(manip, original = "REPLICATION\n", maths = "THEORETICAL NULL\n",
  1435. diff = "REPLICATION – THEORETICAL NULL\n")),
  1436. direction = factor(entropy, levels = rev(c("RD","MM","MP","OR"))))
  1437. ```
  1439. ```{r}
  1440. #| warning: false
  1441. #| lightbox:
  1442. #| group: global
  1444. ## PLOT
  1445. # Initialize subplots
  1446. g.maths.base <- list()
  1447. g.maths.diff <- list()
  1449. # --- ORIGINAL & THEORETICAL NULL
  1450. df.plot.maths %>%
  1451. filter(manip != "REPLICATION – THEORETICAL NULL\n") %>%
  1452. plot_base(col.pal = col.pal.base,
  1453. z.breaks = acc.breaks.base,
  1454. z.labels = acc.labels.base) +
  1455. facet_grid(direction ~ manip) -> g
  1457. g.maths.base[["acc"]] <- g %>% add_onset(df.labels.onset[["sounds"]] %>%
  1458. mutate(manip = str_replace(manip, "EMPIRICAL", "THEORETICAL")))
  1460. # --- DIFFERENCE
  1461. df.plot.maths %>%
  1462. filter(manip == "REPLICATION – THEORETICAL NULL\n") %>%
  1463. plot_base(col.pal = col.pal.diff,
  1464. z.breaks = acc.breaks.diff,
  1465. z.labels = acc.labels.diff) +
  1466. facet_grid(direction ~ manip) +
  1467. theme(axis.title.y = element_blank()) -> g
  1469. g.maths.diff[["acc"]] <- g %>% add_onset(df.labels.onset[["sounds"]] %>%
  1470. mutate(manip = str_replace(manip, "EMPIRICAL", "THEORETICAL")))
  1472. g.fig1.2.acc <- (
  1473. g.maths.base[["acc"]] +
  1474. (g.maths.diff[["acc"]] + plot_layout(tag_level = "new")) +
  1475. plot_layout(widths = c(2, 1)) &
  1476. theme(strip.text.x = element_text(face = "bold"))
  1477. )
  1479. g.fig1.2.acc
  1480. ```
  1482. ### Theoretical null: correlations plots
  1484. ```{r}
  1485. #| warning: false
  1486. #| lightbox:
  1487. #| group: global
  1489. df.plot <- df.plot.maths.corr %>%
  1490. mutate(manip = str_replace(manip, "ORIGINAL", "REPLICATION"),
  1491. manip = paste0(manip, "\n"))
  1493. # ORIGINAL & THEORETICAL
  1494. df.plot %>%
  1495. filter(manip != "REPLICATION – THEORETICAL NULL\n") %>%
  1496. mutate(direction = "") %>%
  1497. plot_base(col.pal = col.pal.corr,
  1498. z = "r",
  1499. z.breaks = corr.breaks,
  1500. z.labels = corr.labels) +
  1501. labs(fill = "Correlation") +
  1502. facet_grid(direction ~ manip) -> g
  1504. g.maths.base[["corr"]] <- g
  1506. # DIFFERENCE
  1507. df.plot %>%
  1508. filter(manip == "REPLICATION – THEORETICAL NULL\n") %>%
  1509. mutate(direction = "") %>%
  1510. plot_base(col.pal = col.pal.corr,
  1511. z = "r",
  1512. z.breaks = corr.breaks,
  1513. z.labels = corr.labels) +
  1514. facet_grid(. ~ manip) +
  1515. labs(fill = "Correlation") +
  1516. theme(axis.title.y = element_blank()) -> g
  1519. g.maths.diff[["corr"]] <- g
  1521. g.fig1.2.corr <- (
  1522. g.maths.base[["corr"]] +
  1523. (g.maths.diff[["corr"]] + plot_layout(tag_level = "new")) +
  1524. plot_layout(widths = c(2, 1), guides = "collect") &
  1525. theme(strip.text.x = element_blank(),
  1526. legend.position = "bottom",
  1527. legend.margin = margin(t = -20, l = 60))
  1528. )
  1530. g.fig1.2.corr
  1531. ```
  1533. ### Assemble Figure 1
  1535. ```{r}
  1536. #| fig-width: 7
  1537. #| fig-height: 11
  1538. #| warning: false
  1539. #| lightbox:
  1540. #| group: global
  1542. ((g.transitions + theme(panel.spacing.x = unit(5,"mm"), legend.box.margin = margin(t = -30))) +
  1543. g.reordering + theme(axis.text.x = element_text(margin = margin(t = -15))) +
  1544. plot_layout(widths = c(1.7, 1))) /
  1545. ((g.fig1.1.acc / g.fig1.1.corr & theme(axis.text.x = element_text(angle = 45, hjust = 1))) +
  1546. plot_layout(heights = c(4.7, 1))) /
  1547. ((g.fig1.2.acc / g.fig1.2.corr & theme(axis.text.x = element_text(angle = 45, hjust = 1))) +
  1548. plot_layout(heights = c(4.7, 1))) +
  1549. plot_layout(heights = c(1., 4, 4)) +
  1550. plot_annotation(tag_levels = "a") -> g.fig1
  1552. g.fig1
  1554. if (save) {
  1555. g.fig1 %>% ggsave(filename = "./figures/fig1_MattersArising.pdf", width = 7, height = 11, device = cairo_pdf)
  1556. }
  1557. ```
  1559. ### Assemble Supplementary Figure 2 (omissions)
  1561. ```{r}
  1562. #| fig-width: 7
  1563. #| fig-height: 5.5
  1564. #| warning: false
  1566. ((g.supfig2.acc / g.supfig2.corr & theme(axis.text.x = element_text(angle = 45, hjust = 1))) +
  1567. plot_layout(heights = c(4, 1))) +
  1568. plot_annotation(tag_levels = "a") -> g.supfig2
  1570. g.supfig2
  1572. if (save) {
  1573. g.supfig2 %>% ggsave(filename = "./figures/supfig2_MattersArising.pdf", width = 7, height = 5.5, device = cairo_pdf)
  1574. g.supfig2 %>% ggsave(filename = "./figures/supfig2_MattersArising.png", width = 7, height = 5.5, device = ragg::agg_png(), dpi = 600)
  1575. }
  1576. ```
  1578. ## Figure 2
  1580. ### Initialize
  1582. ```{r}
  1583. g.base.ml <- list()
  1584. g.diff.ml <- list()
  1585. ```
  1587. ### Build accuracy plots
  1589. ```{r}
  1590. #| lightbox:
  1591. #| group: global
  1593. for (s in c("sounds")) {
  1594. # ORIGINAL & REORDERED
  1595. df.plot.ml %>%
  1596. filter(t_test >=-333, t_test <= 0) %>%
  1597. filter(stim == s, manip != "ORIGINAL – REORDERED") %>%
  1598. mutate(manip = recode_factor(manip,
  1599. ORIGINAL = "ORIGINAL",
  1600. REORDERED = "EMPIRICAL NULL")) %>%
  1601. plot_base(col.pal = col.pal.base,
  1602. z.breaks = acc.breaks.base,
  1603. z.labels = acc.labels.base,
  1604. legend.position = "right") +
  1605. facet_grid(direction ~ manip) +
  1606. labs(title = paste0("Decoding MOST LIKELY sounds")) +
  1607. theme(plot.title = element_text(hjust = -0.2)) -> g.base.ml[["acc"]][[s]]
  1609. # DIFFERENCE
  1610. df.plot.ml %>%
  1611. filter(t_test >=-333, t_test <= 0) %>%
  1612. filter(!grepl("AVERAGE",direction)) %>%
  1613. filter(stim == s, manip == "ORIGINAL – REORDERED") %>%
  1614. mutate(manip = "DIFF") %>%
  1615. plot_base(col.pal = col.pal.diff,
  1616. z.breaks = acc.breaks.diff,
  1617. z.labels = acc.labels.diff,
  1618. legend.position = "right") +
  1619. facet_grid(direction ~ manip) +
  1620. theme(axis.title.y = element_blank()) -> g.diff.ml[["acc"]][[s]]
  1621. }
  1623. (g.base.ml[["acc"]][["sounds"]] + (g.diff.ml[["acc"]][["sounds"]] + plot_layout(tag_level = "new")) + plot_layout(widths = c(2, 1))) +
  1624. theme(plot.tag.position = c(0, 1)) -> g.fig2.acc
  1626. g.fig2.acc
  1627. ```
  1629. ### Build correlations plots
  1631. ```{r}
  1632. #| lightbox:
  1633. #| group: global
  1635. for (s in c("sounds")) {
  1636. # ORIGINAL & REORDERED
  1637. df.plot.ml.corr %>%
  1638. filter(t_test >=-333, t_test <= 0) %>%
  1639. filter(stim == s, manip != "ORIGINAL – REORDERED") %>%
  1640. mutate(manip = recode_factor(manip,
  1641. ORIGINAL = "ORIGINAL",
  1642. REORDERED = "EMPIRICAL NULL")) %>%
  1643. mutate(direction = "") %>%
  1644. plot_base(col.pal = col.pal.corr,
  1645. z = "r",
  1646. z.breaks = corr.breaks,
  1647. z.labels = corr.labels,
  1648. legend.position = "right") +
  1649. labs(fill = "Correlation") +
  1650. facet_grid(direction ~ manip) -> g.base.ml[["corr"]][[s]]
  1652. # DIFFERENCE
  1653. df.plot.ml.corr %>%
  1654. filter(t_test >=-333, t_test <= 0) %>%
  1655. filter(stim == s, manip == "ORIGINAL – REORDERED") %>%
  1656. mutate(manip = "DIFF") %>%
  1657. mutate(direction = "") %>%
  1658. plot_base(col.pal = col.pal.corr,
  1659. z = "r",
  1660. z.breaks = corr.breaks,
  1661. z.labels = corr.labels,
  1662. legend.position = "right") +
  1663. facet_grid(. ~ manip) +
  1664. labs(fill = "Correlation") +
  1665. theme(axis.title.y = element_blank()) -> g.diff.ml[["corr"]][[s]]
  1666. }
  1668. (g.base.ml[["corr"]][["sounds"]] + g.diff.ml[["corr"]][["sounds"]] + plot_layout(widths = c(2, 1))) +
  1669. # plot_annotation(tag_levels = list("b.")) &
  1670. theme(plot.tag.position = c(0, 1)) -> g.fig2.corr
  1672. g.fig2.corr
  1673. ```
  1675. ### Assemble figure 2
  1677. ```{r}
  1678. #| fig-width: 8
  1679. #| fig-height: 5.5
  1680. #| warning: false
  1681. #| lightbox:
  1682. #| group: global
  1684. (g.mostlikely + coord_fixed(0.5))+
  1685. plot_spacer() +
  1686. ((g.fig2.acc / (g.fig2.corr + plot_layout(tag_level = 'new'))) +
  1687. plot_layout(heights = c(3, 1)) &
  1688. theme(plot.title = element_text(vjust = -8),
  1689. legend.box.margin = margin(l=-5),
  1690. strip.clip = "off",
  1691. strip.text = element_text(size = 8),
  1692. axis.text.x = element_text(angle = 45, hjust = 1))
  1693. ) +
  1694. plot_layout(widths = c(1, 0.1, 2)) +
  1695. plot_annotation(tag_levels = "a") &
  1696. theme(plot.tag.position = c(0, 1),
  1697. plot.tag = element_text(vjust = 3)) -> g.fig2
  1699. g.fig2
  1702. if (save) {
  1703. g.fig2 %>% ggsave(filename = "./figures/fig2_MattersArising.pdf", width = 8, height = 5.5, device = cairo_pdf)
  1704. }
  1705. ```
  1707. ## Supplementary Figure 1
  1709. ### Prepare data
  1711. ```{r}
  1712. # Load subject-level data (used for clustering) and derive grand average (used for visualization) **without cropping**
  1713. tic()
  1714. df.supfig1.mean <- tibble()
  1715. df.supfig1.subj <- tibble()
  1716. for (manip in c("","_reord")) {
  1718. for (direction in c("rd_to_rd", "rd_to_mm", "rd_to_mp", "rd_to_or")) {
  1719. tmp <- load_loop_subj(folder = "sounds",
  1720. condition = paste0(direction, manip),
  1721. crop = FALSE)
  1723. df.supfig1.mean %<>% bind_rows(tmp[[2]] %>% mutate(manip = manip, direction = direction))
  1724. df.supfig1.subj %<>% bind_rows(tibble(manip = manip, direction = direction, array = list(tmp[[1]])))
  1725. }
  1726. }
  1727. toc() # ~1s
  1729. ```
  1731. ### Plot
  1733. ```{r}
  1734. #| fig-width: 4
  1735. #| fig-height: 2.5
  1736. #| lightbox:
  1737. #| group: global
  1739. # Extract the diagonal elements from the temporal generalization matrix
  1740. tmp <- filter(df.supfig1.subj, direction == "rd_to_rd", manip == "")[[1, "array"]][[1]]
  1742. tmp.diag <- matrix(nrow = dim(tmp)[3], ncol = dim(tmp)[1])
  1743. for (i in 1:dim(tmp)[3]) {
  1744. tmp.diag[i,] <- tmp[,,i] %>% diag
  1745. }
  1747. # Plot
  1748. df.supfig1.mean %>%
  1749. filter(t_train == t_test,
  1750. direction == "rd_to_rd",
  1751. manip == "") %>%
  1752. mutate(sd = tmp.diag %>% apply(MARGIN = 2, FUN = sd) %>% `/`(sqrt(33))) %>%
  1753. filter(t_train >= -300, t_train <= 700) %>%
  1755. ggplot(aes(x = t_train, y = accuracy)) +
  1756. geom_hline(yintercept = 0.25, color = "grey50", linetype = 2) +
  1757. geom_ribbon(aes(ymin = accuracy-sd, ymax = accuracy+sd), alpha = .3) +
  1758. geom_line(size = 0.8) +
  1759. scale_x_continuous(breaks = scales::pretty_breaks(10), name = "Time (ms)") +
  1760. scale_y_continuous(breaks = seq(.20, .40, .02), limits = c(.2,.4), name = "Accuracy") +
  1761. coord_cartesian(expand = F)-> g.supfig1
  1763. g.supfig1
  1765. if (save) {
  1766. g.supfig1 %>% ggsave(filename = "./figures/supfig1_MattersArising.pdf", width = 4, height = 2.5, device = cairo_pdf)
  1767. g.supfig1 %>% ggsave(filename = "./figures/supfig1_MattersArising.png", width = 4, height = 2.5, device = ragg::agg_png(), dpi = 600)
  1768. }
  1770. ```
  1772. ## Supplementary Figure 3
  1774. ### Functions & parameters
  1776. ```{r}
  1778. # Function to calculate and extract BF
  1779. getBF01 <- function(y, prior = "medium") {
  1780. ttestBF(y, rscale = prior) %>% extractBF(onlybf = T)
  1781. }
  1783. # Color palette for BF
  1784. col.pal.bf <- rev(pals::brewer.piyg(7)); col.pal.bf[c(4)] <- "grey50";
  1786. ```
  1788. ### Prepare data
  1790. ```{r}
  1792. tic()
  1794. ## BF for EMPIRICAL null, DIFF
  1795. # --- extract array of subject-level correlation between accuracy and entropy
  1796. X <- df.corr.subj[[3,"array"]][[1]]
  1797. # --- calculate BF for each cell of the TG matrix
  1798. bf.corr.diff.empir <- apply(X, c(1,2), getBF01, simplify = TRUE)
  1801. ## BF for THEORETICAL null, DIFF
  1802. # --- extract array of subject-level correlation between accuracy and entropy
  1803. X <- df.corr.maths.subj[[3,"array"]][[1]]
  1804. # --- calculate BF for each cell of the TG matrix
  1805. bf.corr.diff.maths <- apply(X, c(1,2), getBF01, simplify = TRUE)
  1807. ## BF for MOST LIKELY null, DIFF
  1808. # --- extract array of subject-level correlation between accuracy and entropy
  1809. X <- df.corr.ml.subj[[3,"array"]][[1]]
  1810. # --- calculate BF for each cell of the TG matrix
  1811. bf.corr.diff.ml <- apply(X, c(1,2), getBF01, simplify = TRUE)
  1813. toc() # ~1.5min
  1816. ## Post-processing
  1817. df.plot <- bind_rows(
  1818. # --- aggregate data
  1819. left_join(matrix_to_df2(bf.corr.diff.empir, "t_train", "t_test", "BF"), # BF results
  1820. filter(df.clusters.corr, manip == "original", stim == "sounds")) %>% # append correlation clusters
  1821. mutate(approach = "EMPIRICAL\nNULL"),
  1822. left_join(matrix_to_df2(bf.corr.diff.maths, "t_train", "t_test", "BF"), # BF results
  1823. filter(df.clusters.corr, manip == "original", stim == "sounds")) %>% # append correlation clusters
  1824. mutate(approach = "THEORETICAL\nNULL"),
  1825. left_join(matrix_to_df2(bf.corr.diff.ml, "t_train", "t_test", "BF"), # BF results
  1826. filter(df.clusters.corr, manip == "original", stim == "sounds")) %>% # append correlation clusters
  1827. mutate(approach = "MOST LIKELY")) %>%
  1829. # --- order factor levels
  1830. mutate(approach = factor(approach, levels = c("EMPIRICAL\nNULL", "THEORETICAL\nNULL", "MOST LIKELY"))) %>%
  1832. # --- convert indices to time
  1833. mutate(t_test = (dt * (t_test - 1) - time_offset)) %>%
  1834. mutate(t_train = (dt * (t_train - 1))) %>%
  1835. # --- discretized version of BF
  1836. mutate(BFd = cut(BF,
  1837. breaks = c(0, 1/30, 1/10, 1/3, 3, 10, 30, Inf),
  1838. labels = c("very strong in favor of null (BF < 1/30)",
  1839. "strong in favor of null (1/30 < BF < 1/10)",
  1840. "moderate in favor of null (1/10 < BF < 1/3)",
  1841. "inconclusive (1/3 < BF < 3)",
  1842. "moderate in favor of alt. (3 < BF < 10)",
  1843. "strong in favor of alt. (10 < BF < 30)",
  1844. "very strong in favor of alt. (BF > 30)")),
  1845. .after = 3
  1846. ) %>%
  1847. # --- transform continuous BF for easy symetric display on colorbar
  1848. mutate(BF = case_when(BF < 1/3 ~ -1/BF,
  1849. between(BF, 1/3, 3) ~ NA,
  1850. BF > 3 ~ BF))
  1852. df.plot %>% summary
  1853. ```
  1855. ### Plot
  1857. ```{r}
  1858. #| fig-width: 10
  1859. #| fig-height: 3
  1860. #| lightbox:
  1861. #| group: global
  1863. lines.v.pos <- 330*seq(round(min(df.plot$t_test)/333), round(max(df.plot$t_test)/333), 1)
  1864. lines.v.pos <- lines.v.pos[lines.v.pos!=0]
  1865. lines.h.pos <- 330*seq(round(min(df.plot$t_train)/333), round(max(df.plot$t_train)/333), 1)
  1866. lines.h.pos <- lines.h.pos[lines.h.pos!=0]
  1869. df.plot %>%
  1870. ggplot(aes(x = t_test, y = t_train)) +
  1871. facet_grid(approach ~ .) +
  1872. geom_tile(aes(fill = BFd), show.legend = TRUE) +
  1874. geom_vline(xintercept = lines.v.pos, linetype = 2, color = "black") +
  1875. geom_hline(yintercept = lines.h.pos, linetype = 2, color = "black") +
  1876. geom_vline(xintercept = 0, color = "black") +
  1877. geom_hline(yintercept = 0, color = "black") +
  1879. geom_contour(aes(z = signif), size = 0.2, color = col.clusters) +
  1881. scale_x_continuous(breaks = scales::pretty_breaks(12),
  1882. expand = c(0,0)) +
  1883. scale_y_continuous(breaks = scales::pretty_breaks(4),
  1884. expand = c(0,0)) +
  1885. scale_fill_manual(values = rev(col.pal.bf),
  1886. labels = levels(df.plot$BFd),
  1887. drop = FALSE) +
  1889. labs(fill = expression(BF[10])) +
  1890. coord_equal() +
  1892. labs(x = "Test time (ms)", y = "Train time (ms)") -> g.supfig3
  1894. g.supfig3
  1896. if (save) {
  1897. g.supfig3 %>% ggsave(filename = "./figures/supfig3_MattersArising.pdf", width = 10, height = 3, device = cairo_pdf)
  1898. g.supfig3 %>% ggsave(filename = "./figures/supfig3_MattersArising.png", width = 10, height = 3, device = ragg::agg_png(), dpi = 600)
  1899. }
  1901. ```
  1903. ## Supplementary Figure 4
  1905. ```{r}
  1906. #| fig-width: 7
  1907. #| fig-height: 2.5
  1908. #| lightbox:
  1909. #| group: global
  1911. mat.C.all <- array(dim = c(34,34,4,4,n.subj))
  1912. i <- 1
  1914. for (s in list.subj) {
  1916. # --- load the confusion matrix
  1917. tmp.file <- os$path$join(path.root, s, "rd_to_rd_confmats.npz") %>%
  1918. np$load()
  1920. # --- crop around S0 (decoded stimulus)
  1921. mat.C.all[,,,,i] <- tmp.file$f[["arr_0"]][,71:104,71:104,,] %>%
  1922. apply(c(2,3,4,5), mean)
  1924. i <- i+1
  1925. }
  1927. mat.C.mean <- mat.C.all %>% apply(c(1,2,3,4), mean)
  1930. df.confusion <- matrix_to_df2(mat.C.mean[22,22,,], "from", "to", "p") %>%
  1931. mutate(label = paste0("a[",from,to,"]"))
  1932. df.confusion %>%
  1933. ggplot(aes(x = to, y = from, fill = p)) +
  1934. geom_tile(color = "firebrick3") +
  1935. geom_text(aes(label = label), parse = T, color = "firebrick3") +
  1936. scale_x_continuous(name = "predicted class", position = "top") +
  1937. scale_y_reverse(name = "true class") +
  1938. scale_fill_stepsn(name = "p",
  1939. colors = pals::brewer.greys(100),
  1940. # low = "white", high = "black",
  1941. labels = ~str_replace(., "^0.","."),
  1942. breaks = scales::pretty_breaks(6),
  1943. limits = c(.20,.30), oob = scales::squish) +
  1944. coord_equal(xlim = c(0.5,4.5), ylim = rev(c(0.5,4.5)),
  1945. expand = F, clip = "off") -> g.matC
  1947. matrix_to_df2(mat.T$MP, "from", "to", "p") %>%
  1948. mutate(label = paste0("a[",from,to,"]")) %>%
  1949. ggplot(aes(x = to, y = from, fill = p)) +
  1950. geom_tile(color = "firebrick3") +
  1951. geom_text(aes(label = label), parse = T, color = "firebrick3") +
  1952. scale_x_continuous(name = "to", position = "top") +
  1953. scale_y_reverse(name = "from") +
  1954. scale_fill_stepsn(name = "p",
  1955. colors = pals::brewer.greys(100),
  1956. labels = ~str_replace(., "^0.","."),
  1957. breaks = scales::pretty_breaks(6),
  1958. limits = c(.0,.50), oob = scales::squish) +
  1959. coord_equal(xlim = c(0.5,4.5), ylim = rev(c(0.5,4.5)),
  1960. expand = F, clip = "off") -> g.matT
  1963. ((g.matT + labs(caption = "Transition matrix\n(midplus (MP) sequence)") +
  1964. guides(fill = guide_colorbar(barwidth = 0.5, ticks = F))) +
  1965. (g.matC + labs(caption = "Confusion matrix\n(group average, t=210ms)") +
  1966. guides(fill = guide_colorbar(barwidth = 0.5, ticks = F))) &
  1967. theme(plot.caption = element_text(size = 10),
  1968. axis.title = element_text(face = "italic"),
  1969. legend.box.spacing = unit(2,"mm"),
  1970. )
  1971. ) +
  1972. plot_annotation(tag_levels = "a") -> g.supfig4
  1974. g.supfig4
  1976. if (save) {
  1977. g.supfig4 %>% ggsave(filename = "./figures/supfig4_MattersArising.pdf", width = 7, height = 2.5, device = cairo_pdf)
  1978. g.supfig4 %>% ggsave(filename = "./figures/supfig4_MattersArising.png", width = 7, height = 2.5, device = ragg::agg_png(), dpi = 600)
  1979. }
  1980. ```

postprocessing.qmd at commit a454439, no license · at the source

Overview

Authors: Oussama Abdoun1, Dmitrii Todorov1,2, Arnaud Poublan-couzardot1, Coumarane Tirou1, Antoine Lutz1, Marine Vernet3, Romain Quentin1
ORCID iDs: Oussama Abdoun, Dmitrii Todorov, Antoine Lutz, Marine Vernet, Romain Quentin
  1. Université Claude Bernard Lyon 1, INSERM U1028, CNRS UMR5292, Lyon Neuroscience Research Center (CRNL), EDUWELL team,Lyon, France
  2. NCP Team, The Biomedical Imaging Laboratory (LIB), CNRS (UMR 7371), INSERM (U1146), Sorbonne Université,Paris, France
  3. Université Claude Bernard Lyon 1, INSERM U1028, CNRS UMR5292, Lyon Neuroscience Research Center (CRNL), IMPACT team,Lyon, France
Institutions: Centre de Recherche en Neurosciences de Lyon (France); Sorbonne Université (France)
Journal: Nature communications, volume 17, issue 1, article 4639
Dates: received 2 July 2024; accepted 8 May 2026; published online 26 May 2026
Type: Letter · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-73568-1 · PMID 42192112 · PMCID PMC13212915 · OpenAlex W4414358687
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.09.16.676331
Status: code verified
Categories: cognitive (subfield)
Methods: Spectral & time-frequency, Connectivity, Statistics, Machine learning
Keywords: Learning and memory, Computational neuroscience, Auditory system, Perception
Journal subjects: Matters Arising
Topic: Neural and Behavioral Psychology Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Agence Nationale de la Recherche (ANR-11-LABX-0042); Institut National de la Santé et de la Recherche Médicale; Université Claude Bernard Lyon 1 (ED476, ANR-11-LABX-0042)
Citations: cited by 1 paper (Europe PMC); 12 references in the paper

Abstract

No abstract was found for this paper: the paper, at the publisher.

Repositories

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

doi:10.5281/zenodo

License: none: the authors keep all their rights
State: the link is dead, verified on 28 September 2026
Evidence: found in the paper
Software Heritage: not checked
Found in: the text, “Data and preprocessing”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 28 September 2026: the link is dead (HTTP 404)
  • 28 September 2026: the link is dead (HTTP 404)
At the source: doi.org/10.5281/zenodo

MEL-Eduwell-lab/predictive_activity

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: a454439aec8936ff9de2199ab255472c7d140b8c, 25 November 2025
Languages: Python (3), JavaScript (1), Quarto (1)
Size: 101 files, 5 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, environment (environment.yml), documentation, 1 notebook
Not found: license file, CITATION.cff, tests, continuous integration
Tools: MNE-Python (3 files), NumPy (3 files), pandas (3 files), scikit-learn (2 files), BayesFactor (1 file), easystats (1 file), patchwork (1 file), reticulate (1 file), seaborn (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
6 files
At the source: github.com/MEL-Eduwell-lab/predictive_activity

romquentin/predictive_activity

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: a454439aec8936ff9de2199ab255472c7d140b8c, 25 November 2025
Languages: Python (3), JavaScript (1), Quarto (1)
Size: 101 files, 5 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, environment (environment.yml), documentation, 1 notebook
Not found: license file, CITATION.cff, tests, continuous integration
Tools: MNE-Python (3 files), NumPy (3 files), pandas (3 files), scikit-learn (2 files), BayesFactor (1 file), easystats (1 file), patchwork (1 file), reticulate (1 file), seaborn (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
6 files
At the source: github.com/romquentin/predictive_activity

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-73568-1.

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:

  • 3 repositories 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;
  • 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

Datasets cited

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • no repository, dataset or request procedure was recognized in it

Read it in the paper: doi.org/10.1038/s41467-026-73568-1.

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

  • Funding: added Agence Nationale de la Recherche: ANR-11-LABX-0042; Institut National de la Santé et de la Recherche Médicale; Université Claude Bernard Lyon 1: ED476, ANR-11-LABX-0042

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 4 keywords, 9 references.

Cite

This paper

Abdoun, O., Todorov, D., Poublan-couzardot, A., Tirou, C., Lutz, A., Vernet, M., & Quentin, R. (2026). No evidence of neural feature-specific pre-activation during the prediction of an upcoming stimulus. Nature communications, 17(1), 4639. https://doi.org/10.1038/s41467-026-73568-1

BibTeX

@article{abdoun2026no,
author = {Abdoun, Oussama and Todorov, Dmitrii and Poublan-couzardot, Arnaud and Tirou, Coumarane and Lutz, Antoine and Vernet, Marine and Quentin, Romain},
title = {{No evidence of neural feature-specific pre-activation during the prediction of an upcoming stimulus}},
journal = {Nature communications},
year = {2026},
month = may,
volume = {17},
number = {1},
pages = {4639},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-73568-1},
url = {https://doi.org/10.1038/s41467-026-73568-1},
pmid = {42192112},
pmcid = {PMC13212915}
}

RIS

TY - JOUR
AU - Abdoun, Oussama
AU - Todorov, Dmitrii
AU - Poublan-couzardot, Arnaud
AU - Tirou, Coumarane
AU - Lutz, Antoine
AU - Vernet, Marine
AU - Quentin, Romain
TI - No evidence of neural feature-specific pre-activation during the prediction of an upcoming stimulus
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/05/26
VL - 17
IS - 1
SP - 4639
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-73568-1
UR - https://doi.org/10.1038/s41467-026-73568-1
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-73568-1",
"type": "article-journal",
"title": "No evidence of neural feature-specific pre-activation during the prediction of an upcoming stimulus",
"container-title": "Nature communications",
"author": [
{
"family": "Abdoun",
"given": "Oussama"
},
{
"family": "Todorov",
"given": "Dmitrii"
},
{
"family": "Poublan-couzardot",
"given": "Arnaud"
},
{
"family": "Tirou",
"given": "Coumarane"
},
{
"family": "Lutz",
"given": "Antoine"
},
{
"family": "Vernet",
"given": "Marine"
},
{
"family": "Quentin",
"given": "Romain"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "4639",
"DOI": "10.1038/s41467-026-73568-1",
"PMID": "42192112",
"PMCID": "PMC13212915",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-73568-1",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
26
]
]
}
}

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.1038/s41467-026-74824-0 [code]
Learning regularities in noise engages both neural predictive activity and representational changes.
Journal: Nature communications
In common: reticulate, easystats, MNE-Python, 4 other tools, cognitive, 4 references, 3 authors
[2] doi:10.1093/nc/niag029 [code]
A data-driven approach to identifying and evaluating connectivity-based neural correlates of conscious visual perception.
Journal: Neuroscience of consciousness
In common: BayesFactor, easystats, MNE-Python, 6 other tools, cognitive
[3] doi:10.1038/s41597-026-07350-9 [code]
An open multi-center MEG-EEG dataset for studying conscious visual perception.
Journal: Scientific data
In common: BayesFactor, easystats, MNE-Python, 5 other tools, 1 reference
[4] doi:10.1038/s41597-026-07377-y [code]
An open-access multi-site fMRI dataset for investigating conscious visual perception.
Journal: Scientific data
In common: BayesFactor, easystats, MNE-Python, 5 other tools, 1 reference
[5] doi:10.1167/jov.26.8.4 [code]
The neural processes of illusory occlusion in object recognition.
Journal: Journal of vision
In common: BayesFactor, MNE-Python, seaborn, 4 other tools, cognitive, 1 reference
[6] doi:10.1016/j.isci.2026.117285 [code]
Working memory demands modulate memory brain state engagement.
Journal: iScience
In common: BayesFactor, easystats, MNE-Python, 4 other tools, cognitive
[7] doi:10.1038/s41467-026-75662-w [code]
Distinct Roles of Deep and Superficial Cortical Layers in Tone Prediction, Comparison, and Adaptation in Human Auditory Cortices.
Journal: Nature communications
In common: MNE-Python, seaborn, tidyverse, 3 other tools, cognitive, 2 references
[8] doi:10.3390/s26103131 [code]
Neuromagnetism "On the Cheap": Evaluating a Combined Cylindrical Shield and Partial-Coverage OPM-MEG System for Detecting Sensorimotor Responses in Humans.
Journal: Sensors (Basel, Switzerland)
In common: BayesFactor, reticulate, MNE-Python, 4 other tools
[9] doi:10.1002/hbm.70605 [code]
BrainEnrich: Revealing Biological Insights for Imaging-Derived Phenotypes Through Transcriptomic Enrichment.
Journal: Human brain mapping
In common: reticulate, easystats, patchwork, 5 other tools
[10] doi:10.7554/elife.107088 [code]
Development of auditory and spontaneous movement responses to music over the first postnatal year.
Journal: eLife
In common: easystats, MNE-Python, seaborn, 4 other tools, 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.