OSCR

Corpus Callosum Dysgenesis impairs metacognition: evidence from multi-modality and multi-cohort replications

Code ↔ Paper

4 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 4 matches
  1. [1] § Methods › Random Dot Kinematogram (RDK) › Experiment 1 ↔ Scripts/Analysis_Script_Commented_Github.r, lines 678–726 · score 0.64 · posterior predictive checks, HMeta, ratio, binning, fitted, models
  2. [2] § Methods › Random Dot Kinematogram (RDK) › Experiment 1 ↔ Scripts/Analysis_Script_Commented_Github.r, lines 228–273 · score 0.57 · 50–100 %, position, coherently, confidence
  3. [3] § Results › Experiment 3 › Behavioral Analysis ↔ Scripts/Analysis_Script_Commented_Github.r, lines 1264–1351 · score 0.52 · monocular trials, lateralized trials, binocular trials, scores, VR, NT
  4. [4] § Methods › Random Dot Kinematogram (RDK) › Experiment 1 ↔ Scripts/Analysis_Script_Commented_Github.r, lines 2–43 · score 0.51 · perceptual accuracy, RDK task, calibration, online, efficiency, NT

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 · 1,351 lines · 66 KB · CC-BY-NC-SA-4.0 · 4 matches

  2. # Barnby et al. (2026): Corpus Callosum Dysgenesis impairs metacognition:
  3. # evidence from multi-modality and multi-cohort replications
  4. #
  5. # Analysis script covering three experiments:
  6. # Experiment 1: Online RDK (CCD vs NT, computer-based)
  7. # Experiment 2: In-lab RDK (CCD vs NT, MRI environment)
  8. # Experiment 3: VR RDK (CCD vs NT, binocular / monocular / lateralized)
  9. #
  10. # Key outcome measures: perceptual accuracy, confidence, metacognitive
  11. # efficiency (hierarchical meta-d' / d'; HMeta-D method, Maniscalco & Lau 2012)
  13. # ── 1. SETUP ------------------------------------------------------------
  15. library(ggpubr); library(tidyverse); library(patchwork)
  16. library(doParallel);library(rstanarm); library(brms)
  17. library(magrittr); library(reshape2); library(rjags)
  18. library(coda); library(lattice); library(lme4)
  19. library(broom); library(ggmcmc); library(foreach)
  20. library(sjPlot); library(tidyplots); library(tidybayes)
  22. source('Scripts/RDKUtilityFuncs.R')
  24. # Colour palettes used throughout all figures
  25. colpal <- c('#931F1D', '#FE938C', '#475B5A', '#81D6E3') # Exp 1 & 2 (4 groups)
  26. colpal_vr <- c('#B80C09', '#0E1C36') # Exp 3 (CCD / NT)
  28. # ── 2. LOAD DATA ------------------------------------------------------------
  30. # Experiments 1 (online) and 2 (MRI/lab): conventional RDK task
  31. Experiment_1_CCD <- read.csv('Data/Experiment1_CCD.csv')
  32. Experiment_1_NT <- read.csv('Data/Experiment1_NT.csv')
  33. Experiment_2_CCD <- read.csv('Data/Experiment2_CCD.csv')
  34. Experiment_2_NT <- read.csv('Data/Experiment2_NT.csv')
  36. # Experiment 3: VR task with binocular / monocular / lateralized presentations
  37. data_vr <- read.csv('Data/Experiment3.csv') %>% select(-X, -X.1)
  39. # ── 3. DATA QUALITY: CALIBRATION & EXCLUSIONS ------------------------------------------------------------
  40. # Participants were excluded if their coherence threshold during calibration
  41. # exceeded 0.5, if confidence ratings showed no variance, or if self-reported
  42. # error counts were implausibly high. Trials on which participants reported a
  43. # button-press mistake are also removed.
  45. ## --- 3a. Visual diagnostics ----
  46. # Inspect coherence trajectories across calibration trials to flag outliers
  47. (Experiment_1_NT %>% check_coherence('NT (Online)', legend = F) +
  48. Experiment_1_CCD %>% check_coherence('CCD (Online)', legend = F)) /
  49. (Experiment_2_NT %>% check_coherence('NT (MRI)', legend = F) +
  50. Experiment_2_CCD %>% check_coherence('CCD (MRI)', legend = F))
  52. # Spot-check specific outlier participants identified during calibration review
  53. Experiment_2_CCD %>%
  54. filter(ID %in% c('1001', '1027')) %>%
  55. check_coherence('CCD (MRI)', legend = F)
  57. # Inspect confidence-rating distributions per participant (ridge plots)
  58. # — applied to each of the four groups; code is identical so only one shown
  59. for (dat in list(Experiment_1_CCD, Experiment_2_CCD, Experiment_1_NT, Experiment_2_NT)) {
  60. P_dat <- dat %>%
  61. filter(conf != 0, type == 'main', mistake %in% c(0, NA, FALSE)) %>%
  62. ggplot(aes(conf, factor(ID), fill = ID, group = ID)) +
  63. ggridges::geom_density_ridges(scale = 1, jittered_points = TRUE,
  64. position = "raincloud", alpha = 0.7, scale = 0.2) +
  65. labs(x = 'Confidence', y = 'ID') +
  66. theme_bw() +
  67. theme(text = element_text(size = 18), legend.position = 'none')
  68. print(P_dat)
  69. }
  71. # Count self-reported button-press errors per participant (used to set
  72. # the ">20 mistakes" exclusion threshold applied below)
  73. for (dat in list(Experiment_1_CCD, Experiment_2_CCD, Experiment_1_NT, Experiment_2_NT)) {
  74. P_dat <- dat %>%
  75. filter(conf != 0, type == 'main') %>%
  76. group_by(ID) %>%
  77. mutate(mistake_sum = ifelse(!is.na(mistake) & mistake == TRUE, 1, 0),
  78. mistake_sum = sum(mistake_sum)) %>%
  79. dplyr::select(mistake_sum, ID) %>%
  80. unique() %>%
  81. ggplot(aes(factor(ID), mistake_sum, fill = ID)) +
  82. geom_col() +
  83. coord_cartesian(ylim = c(1, 25)) +
  84. labs(x = 'ID', y = 'Mistake') +
  85. theme_bw() +
  86. theme(text = element_text(size = 18), legend.position = 'none')
  87. print(P_dat)
  88. }
  90. ## --- 3b. Apply exclusions ----
  91. # Reasons documented per participant:
  92. # 2004 : coherence threshold > 0.5 (Exp 1 CCD)
  93. # 5f20… : coherence threshold > 0.5 (Exp 1 NT)
  94. # 1001 : coherence threshold > 0.5 (Exp 2 CCD)
  95. # 1006 : zero variance in confidence + >20 mistakes (Exp 2 CCD)
  96. # 1027 : >20 mistakes (Exp 2 CCD)
  97. # 3002 : zero variance in confidence (Exp 2 NT)
  98. # 1004 : self-reported mistakes not marked in data (Exp 2 CCD)
  99. # 2008 : uncertain phenotype (Exp 1 CCD); retained in sensitivity set
  101. Experiment_1_CCDcheck <- Experiment_1_CCD %>% filter(ID != '2004')
  102. Experiment_1_CCDcheck_w_exc <- Experiment_1_CCD %>% filter(ID != '2004', ID != '2008')
  103. Experiment_1_NTcheck <- Experiment_1_NT %>% filter(ID != '5f20a850b6d5ed3d186bc3c3')
  104. Experiment_2_CCDcheck <- Experiment_2_CCD %>% filter(!ID %in% c('1001', '1006', '1027'))
  105. Experiment_2_NTcheck <- Experiment_2_NT %>% filter(ID != '3002')
  107. # Remove mistake-flagged trials across all groups
  108. Experiment_1_CCDcheck <- Experiment_1_CCDcheck %>% filter(mistake == FALSE)
  109. Experiment_1_NTcheck <- Experiment_1_NTcheck %>% filter(mistake == FALSE)
  110. Experiment_2_CCDcheck <- Experiment_2_CCDcheck %>% filter(mistake == FALSE)
  111. Experiment_2_NTcheck <- Experiment_2_NTcheck %>% filter(mistake == FALSE)
  113. # Final sample sizes (Experiments 1 & 2)
  114. cat("Exp1 CCD (main):", length(unique(Experiment_1_CCDcheck$ID)), "\n")
  115. cat("Exp1 CCD (w/ 2008 excluded):", length(unique(Experiment_1_CCDcheck_w_exc$ID)), "\n")
  116. cat("Exp1 NT:", length(unique(Experiment_1_NTcheck$ID)), "\n")
  117. cat("Exp2 CCD:", length(unique(Experiment_2_CCDcheck$ID)), "\n")
  118. cat("Exp2 NT:", length(unique(Experiment_2_NTcheck$ID)), "\n")
  120. checkBothRDK <- rbind(
  121. Experiment_1_CCDcheck %>% mutate(correct = as.numeric(correct)),
  122. Experiment_1_NTcheck %>% mutate(correct = as.numeric(correct)),
  123. Experiment_2_CCDcheck %>% mutate(correct = as.numeric(correct)),
  124. Experiment_2_NTcheck %>% mutate(correct = as.numeric(correct))
  125. )
  127. # Sensitivity dataset (Exp 1 CCD with participant 2008 also excluded)
  128. checkBothRDK_w_exc <- rbind(
  129. Experiment_1_CCDcheck_w_exc %>% mutate(correct = as.numeric(correct)),
  130. Experiment_1_NTcheck %>% mutate(correct = as.numeric(correct)),
  131. Experiment_2_CCDcheck %>% mutate(correct = as.numeric(correct)),
  132. Experiment_2_NTcheck %>% mutate(correct = as.numeric(correct))
  133. )
  135. # Binary flag for the high-coherence (kmed × 2) condition
  136. checkBothRDK$highCoh <- ifelse(checkBothRDK$kmed == "kmed x 2", 1, 0)
  137. checkBothRDK_w_exc$highCoh <- ifelse(checkBothRDK_w_exc$kmed == "kmed x 2", 1, 0)
  139. # ── 4. ENVIRONMENT CHECK & DEMOGRAPHICS -------------------------
  140. # We test whether the testing platform (computer / MRI / VR) explains variance
  141. # in accuracy or confidence. If not, experiments can be analysed jointly.
  143. # Accuracy by environment
  144. accuracy_data_vr <- data_vr %>%
  145. rename(ID = PID, environment = Environment) %>%
  146. filter(TrialType == 'Main', Presentation == 'Binocular') %>%
  147. group_by(ID, environment, CCD) %>%
  148. summarise(Accuracy = mean(Correct), .groups = 'drop')
  149. accuracy_data_vr$ID <- as.character(accuracy_data_vr$ID)
  151. accuracy_data_non_vr <- checkBothRDK %>%
  152. filter(type == 'main') %>%
  153. group_by(ID, environment, CCD) %>%
  154. summarise(Accuracy = mean(correct), .groups = 'drop')
  156. accuracy_data <- rbind(accuracy_data_vr, accuracy_data_non_vr)
  158. # Linear model: environment and CCD group as predictors of accuracy
  159. model_acc <- lm(Accuracy ~ environment + CCD, data = accuracy_data)
  160. summary(aov(Accuracy ~ environment + CCD, data = accuracy_data))
  161. plot_model(model_acc, title = "Effect of environment on participant accuracy",
  162. show.values = TRUE, show.p = TRUE)
  163. tab_model(model_acc, show.re.var = TRUE,
  164. dv.labels = "Effect of environment on participant accuracy",
  165. file = "Stats/participant_accuracy_lmer.html")
  167. # Confidence by environment (scale confidence to z-scores before aggregating)
  168. data_vr$Confidence <- scale(data_vr$ReportedVeryConfident)
  169. confidence_data_vr <- data_vr %>%
  170. filter(TrialType == 'Main') %>%
  171. group_by(PID, Environment, CCD) %>%
  172. summarise(Confidence = mean(Confidence), .groups = 'drop') %>%
  173. mutate(PID = as.character(PID)) %>%
  174. rename(ID = PID, environment = Environment)
  176. checkBothRDKConfidenceTrialsOnly <- checkBothRDK %>%
  177. filter(conf != 0) %>%
  178. mutate(Confidence = scale(conf))
  180. confidence_data_non_vr <- checkBothRDKConfidenceTrialsOnly %>%
  181. filter(type == 'main') %>%
  182. group_by(ID, environment, CCD) %>%
  183. summarise(Confidence = mean(Confidence), .groups = 'drop')
  185. confidence_data <- rbind(confidence_data_vr, confidence_data_non_vr)
  186. summary(aov(Confidence ~ environment + CCD, data = confidence_data))
  188. # ── 5. BEHAVIOURAL RESULTS: EXPERIMENTS 1 & 2 -----------------------------------
  189. # Primary outcome: does accuracy / confidence scale with coherence level,
  190. # and does this scaling differ between CCD and NT?
  192. lmer_data <- checkBothRDK %>% filter(type == 'main', kmed != 'kmed', mistake %in% c(0,NA,FALSE), conf != 0)
  193. lmer_data_exc <- checkBothRDK_w_exc %>% filter(type == 'main', kmed != 'kmed', mistake %in% c(0,NA,FALSE), conf != 0)
  194. data_vr_tests <- data_vr %>%
  195. rename(ID=PID, type=TrialType, coherence=ActiveCoherence, group=Group) %>%
  196. filter(type=='Main') %>%
  197. group_by(ID, PresentationType) %>%
  198. mutate(
  199. conf = scale(ReportedVeryConfident),
  200. #conf = 1/(1+exp(-scale(conf))),
  201. kmed = ifelse(HighCoherenceTrial==1, 'x2','x0.5'),
  202. kmed = factor(kmed, levels = c('x0.5', 'x2'), ordered = T)
  203. ) %>%
  204. group_by(ID, group)
  206. # --- 5a. Reaction-time density plots (sanity check) ----
  207. CCDcheck1 <- ggplot(
  208. checkBothRDK %>%
  209. filter(type == 'main', kmed != 'kmed', RT < 5000, mistake %in% c(0,NA,FALSE)) %>%
  210. mutate(correct = ifelse(correct == 1, 'Correct: Yes', 'Correct: No')),
  211. aes(RT, fill = group)) +
  212. geom_density(show.legend = TRUE, alpha = 0.7) +
  213. scale_fill_manual(values = colpal) +
  214. facet_wrap(kmed ~ correct) +
  215. labs(x = 'Reaction Time (ms)') +
  216. theme_minimal() +
  217. theme(axis.title.y = element_blank())
  219. # Non-parametric tests on RT medians (group × environment)
  220. kruskal.test(RT ~ group, data = rbind(Experiment_1_CCDcheck, Experiment_1_NTcheck))
  221. kruskal.test(RT ~ group, data = rbind(Experiment_2_CCDcheck, Experiment_2_NTcheck))
  222. kruskal.test(RT ~ group, data = rbind(Experiment_2_CCDcheck, Experiment_1_CCDcheck))
  223. kruskal.test(RT ~ group, data = rbind(Experiment_1_NTcheck, Experiment_2_NTcheck))
  225. # --- 5b. Per-participant accuracy and confidence by coherence level --------
  226. # Compute individual means across the two coherence conditions
  228. ## EXPERIMENT 1 & 2 -------
  230. av_cor <- checkBothRDK %>%
  231. filter(type == 'main', kmed != 'kmed', mistake %in% c(0,NA,FALSE)) %>%
  232. group_by(group, ID, kmed) %>%
  233. mutate(av_cor = mean(corAdjusted),
  234. kmed = ifelse(kmed == 'kmed x 2', 'x2', 'x0.5')) %>%
  235. group_by(group, kmed) %>%
  236. dplyr::select(ID, av_cor, kmed) %>%
  237. distinct()
  239. av_co <- checkBothRDK %>%
  240. filter(type == 'main', kmed != 'kmed', conf != 0, mistake %in% c(0,NA,FALSE)) %>%
  241. group_by(group, ID, kmed) %>%
  242. mutate(av_co = mean(conf),
  243. kmed = ifelse(kmed == 'kmed x 2', 'x2', 'x0.5')) %>%
  244. group_by(group, kmed) %>%
  245. dplyr::select(ID, av_co, kmed) %>%
  246. distinct()
  248. # Shared aesthetics for accuracy and confidence plots
  249. dot_layer <- list(
  250. geom_jitter(shape = 21, width = 0.1, alpha = 0.1, colour = 'black'),
  251. stat_summary(),
  252. stat_summary(geom = 'line', show.legend = FALSE),
  253. scale_color_manual(values = colpal),
  254. scale_fill_manual(values = colpal),
  255. theme_minimal(),
  256. theme(panel.grid = element_blank(), text = element_text(size = 18))
  257. )
  259. CCDcheck3 <- ggplot(av_cor, aes(kmed, av_cor, color = group, group = group,
  260. fill = group)) +
  261. dot_layer +
  262. labs(x = 'Coherence', y = 'Correct') +
  263. coord_cartesian(ylim = c(0.5, 1)) +
  264. theme(legend.position = 'none')
  266. CCDcheck4 <- ggplot(av_co, aes(kmed, av_co, color = group, group = group,
  267. fill = group)) +
  268. dot_layer +
  269. labs(x = 'Coherence', y = 'Confidence') +
  270. coord_cartesian(ylim = c(50, 100)) +
  271. theme(legend.position = 'none')
  273. (CCDcheck3 | CCDcheck4)
  275. ## EXPERIMENT 3 --------
  277. # Per-participant accuracy and confidence by coherence × presentation type
  278. av_cor_vr <- data_vr_tests %>%
  279. group_by(group, ID, kmed, PresentationType) %>%
  280. summarise(av_cor = mean(Correct), .groups='drop')
  282. av_co_vr <- data_vr_tests %>%
  283. group_by(group, ID, kmed, PresentationType) %>%
  284. summarise(av_co = mean(conf), .groups='drop')
  286. # Accuracy and confidence plots (faceted by presentation type)
  287. vr_cor_plot <- ggplot(av_cor_vr, aes(kmed, av_cor, color = group, group = group, fill = group)) +
  288. geom_jitter(shape = 21, width = 0.1, alpha = 0.1, colour = 'black') +
  289. stat_summary() + stat_summary(geom = 'line', show.legend = FALSE) +
  290. labs(x = 'Coherence', y = 'Correct') +
  291. facet_wrap(~PresentationType) +
  292. scale_color_manual(values = colpal_vr) + scale_fill_manual(values = colpal_vr) +
  293. coord_cartesian(ylim = c(0.5, 1)) + theme_minimal() +
  294. theme(legend.position = 'none', axis.title.x = element_blank())
  296. vr_co_plot <- ggplot(av_co_vr, aes(kmed, av_co, color = group, group = group, fill = group)) +
  297. geom_jitter(shape = 21, width = 0.1, alpha = 0.1, colour = 'black', show.legend = FALSE) +
  298. stat_summary() + stat_summary(geom = 'line', show.legend = FALSE) +
  299. labs(x = 'Coherence', y = 'Confidence') +
  300. facet_wrap(~PresentationType) +
  301. scale_color_manual(values = colpal_vr, name = 'Group') +
  302. scale_fill_manual(values = colpal_vr) +
  303. theme_minimal() +
  304. theme(legend.position = 'none', strip.text.x = element_blank())
  306. (vr_cor_plot / vr_co_plot) & theme(text = element_text(size = 18), panel.grid = element_blank())
  308. # --- 5c. Paired t-tests: coherence effect within each group ----------
  310. for (grp in c('CCD (Online)', 'NT (Online)', 'CCD (MRI)', 'NT (MRI)')) {
  311. cat("\n--- Accuracy:", grp, "---\n")
  312. print(t.test(av_cor ~ kmed, data = av_cor %>% filter(group == grp), paired = TRUE))
  313. cat("\n--- Confidence:", grp, "---\n")
  314. print(t.test(av_co ~ kmed, data = av_co %>% filter(group == grp), paired = TRUE))
  315. }
  317. # Group × coherence ANOVAs (Exp 1 and Exp 2 separately)
  318. aov(av_cor ~ group*kmed, data = av_cor %>%
  319. filter(group %in% c('NT (Online)', 'CCD (Online)'))) %>% summary()
  320. aov(av_cor ~ group*kmed, data = av_cor %>%
  321. filter(group %in% c('NT (MRI)', 'CCD (MRI)'))) %>% summary()
  323. aov(av_co ~ group*kmed, data = av_co %>%
  324. filter(group %in% c('NT (Online)', 'CCD (Online)'))) %>% summary()
  325. aov(av_co ~ group*kmed, data = av_co %>%
  326. filter(group %in% c('NT (MRI)', 'CCD (MRI)'))) %>% summary()
  328. # --- 5d. Linear mixed-effects models (trial-level; primary inference) --------
  329. # kmed × group interaction on accuracy and confidence, with random intercept
  330. # per participant. Sensitivity models include/exclude participant 2008.
  332. ### Experiment 1 ------
  334. run_lmer(correct ~ kmed + (1|ID), lmer_data %>% filter(group == 'CCD (Online)'))
  335. run_lmer(correct ~ kmed + (1|ID), lmer_data %>% filter(group == 'NT (Online)'))
  336. run_lmer(correct ~ kmed * group + (1|ID), lmer_data %>% filter(group %in% c('NT (Online)', 'CCD (Online)')))
  338. run_lmer(conf ~ kmed + (1|ID), lmer_data %>% filter(group == 'CCD (Online)'))
  339. run_lmer(conf ~ kmed + (1|ID), lmer_data %>% filter(group == 'NT (Online)'))
  340. run_lmer(conf ~ kmed * group + (1|ID), lmer_data %>% filter(group %in% c('NT (Online)', 'CCD (Online)')))
  342. ### Experiment 2 --------
  344. run_lmer(correct ~ kmed + (1|ID), lmer_data %>% filter(group == 'CCD (MRI)'))
  345. run_lmer(correct ~ kmed + (1|ID), lmer_data %>% filter(group == 'NT (MRI)'))
  346. run_lmer(correct ~ kmed * group + (1|ID), lmer_data %>% filter(group %in% c('NT (MRI)', 'CCD (MRI)')))
  348. run_lmer(conf ~ kmed + (1|ID), lmer_data %>% filter(group == 'CCD (MRI)'))
  349. run_lmer(conf ~ kmed + (1|ID), lmer_data %>% filter(group == 'NT (MRI)'))
  350. run_lmer(conf ~ kmed * group + (1|ID), lmer_data %>% filter(group %in% c('NT (MRI)', 'CCD (MRI)')))
  352. ### Experiment 3 ---------
  354. run_lmer(Correct ~ kmed * group + (1|ID), data_vr_tests)
  355. run_lmer(conf ~ kmed * group + (1|ID), data_vr_tests)
  357. #### BINO ----
  359. run_lmer(Correct ~ kmed + (1|ID), data_vr_tests %>% filter(group=='CCD', PresentationType=='Binocular'))
  360. run_lmer(Correct ~ kmed + (1|ID), data_vr_tests %>% filter(group=='NT', PresentationType=='Binocular'))
  362. run_lmer(Correct ~ kmed * group + (1|ID), data_vr_tests %>% filter(PresentationType=='Binocular'))
  364. run_lmer(conf ~ kmed + (1|ID), data_vr_tests %>% filter(group=='CCD', PresentationType=='Binocular'))
  365. run_lmer(conf ~ kmed + (1|ID), data_vr_tests %>% filter(group=='NT', PresentationType=='Binocular'))
  367. run_lmer(conf ~ kmed * group + (1|ID), data_vr_tests %>% filter(PresentationType=='Binocular'))
  369. #### LAT ----
  371. run_lmer(Correct ~ kmed + (1|ID), data_vr_tests %>% filter(group=='CCD', PresentationType=='Lateralized'))
  372. run_lmer(Correct ~ kmed + (1|ID), data_vr_tests %>% filter(group=='NT', PresentationType=='Lateralized'))
  374. run_lmer(Correct ~ kmed * group + (1|ID), data_vr_tests %>% filter(PresentationType=='Lateralized'))
  376. run_lmer(conf ~ kmed + (1|ID), data_vr_tests %>% filter(group=='CCD', PresentationType=='Lateralized'))
  377. run_lmer(conf ~ kmed + (1|ID), data_vr_tests %>% filter(group=='NT', PresentationType=='Lateralized'))
  379. run_lmer(conf ~ kmed * group + (1|ID), data_vr_tests %>% filter(PresentationType=='Lateralized'))
  381. #### MONO ----
  383. run_lmer(Correct ~ kmed + (1|ID), data_vr_tests %>% filter(group=='CCD', PresentationType=='Monocular'))
  384. run_lmer(Correct ~ kmed + (1|ID), data_vr_tests %>% filter(group=='NT', PresentationType=='Monocular'))
  386. run_lmer(Correct ~ kmed * group + (1|ID), data_vr_tests %>% filter(PresentationType=='Monocular'))
  388. run_lmer(conf ~ kmed + (1|ID), data_vr_tests %>% filter(group=='CCD', PresentationType=='Monocular'))
  389. run_lmer(conf ~ kmed + (1|ID), data_vr_tests %>% filter(group=='NT', PresentationType=='Monocular'))
  391. run_lmer(conf ~ kmed * group + (1|ID), data_vr_tests %>% filter(PresentationType=='Monocular'))
  393. # --- 5e. Confidence Bin Spread --------
  395. LAB_CCD_metaDat <- HMetaGroupPrep(checkBothRDK, 'CCD (MRI)')
  396. LAB_CCD_metaDat_hc <- HMetaGroupPrep(checkBothRDK %>% filter(kmed == 'kmed x 2'), 'CCD (MRI)')
  397. LAB_CCD_metaDat_lc <- HMetaGroupPrep(checkBothRDK %>% filter(kmed == 'kmed x 0.5'), 'CCD (MRI)')
  399. # Create data set for ONL CCD participants
  400. ONL_CCD_metaDat <- HMetaGroupPrep(checkBothRDK, 'CCD (Online)')
  401. ONL_CCD_metaDat_hc <- HMetaGroupPrep(checkBothRDK %>% filter(kmed == 'kmed x 2'), 'CCD (Online)')
  402. ONL_CCD_metaDat_lc <- HMetaGroupPrep(checkBothRDK %>% filter(kmed == 'kmed x 0.5'), 'CCD (Online)')
  404. # Create data set for ONL NT participants
  405. ONL_NT_metaDat <- HMetaGroupPrep(checkBothRDK, 'NT (Online)')
  406. ONL_NT_metaDat_hc <- HMetaGroupPrep(checkBothRDK %>% filter(kmed == 'kmed x 2'), 'NT (Online)')
  407. ONL_NT_metaDat_lc <- HMetaGroupPrep(checkBothRDK %>% filter(kmed == 'kmed x 0.5'), 'NT (Online)')
  409. # Create data set for LAB CCD participants
  410. LAB_NT_metaDat <- HMetaGroupPrep(checkBothRDK, 'NT (MRI)')
  411. LAB_NT_metaDat_hc <- HMetaGroupPrep(checkBothRDK %>% filter(kmed == 'kmed x 2'), 'NT (MRI)')
  412. LAB_NT_metaDat_lc <- HMetaGroupPrep(checkBothRDK %>% filter(kmed == 'kmed x 0.5'), 'NT (MRI)')
  414. # Confidence spread (normalised response count distributions by coherence)
  415. conf_spread <- data.frame(
  416. conf = c(6:1, -1:-6),
  417. cond = c(rep('lab', 72), rep('online', 72)),
  418. type = rep(c(rep('x2',24), rep('Both',24), rep('x0.5',24)), 2),
  419. s1 = c(rowSums(LAB_CCD_metaDat_hc$S1mat), rowSums(LAB_NT_metaDat_hc$S1mat),
  420. rowSums(LAB_CCD_metaDat$S1mat), rowSums(LAB_NT_metaDat$S1mat),
  421. rowSums(LAB_CCD_metaDat_lc$S1mat), rowSums(LAB_NT_metaDat_lc$S1mat),
  422. rowSums(ONL_CCD_metaDat_hc$S1mat), rowSums(ONL_NT_metaDat_hc$S1mat),
  423. rowSums(ONL_CCD_metaDat$S1mat), rowSums(ONL_NT_metaDat$S1mat),
  424. rowSums(ONL_CCD_metaDat_lc$S1mat), rowSums(ONL_NT_metaDat_lc$S1mat)),
  425. s2 = c(rev(rowSums(LAB_CCD_metaDat_hc$S2mat)), rev(rowSums(LAB_NT_metaDat_hc$S2mat)),
  426. rev(rowSums(LAB_CCD_metaDat$S2mat)), rev(rowSums(LAB_NT_metaDat$S2mat)),
  427. rev(rowSums(LAB_NT_metaDat_lc$S2mat)), rev(rowSums(LAB_NT_metaDat_lc$S2mat)),
  428. rev(rowSums(ONL_CCD_metaDat_hc$S2mat)), rev(rowSums(ONL_NT_metaDat_hc$S2mat)),
  429. rev(rowSums(ONL_CCD_metaDat$S2mat)), rev(rowSums(ONL_NT_metaDat$S2mat)),
  430. rev(rowSums(ONL_NT_metaDat_lc$S2mat)), rev(rowSums(ONL_NT_metaDat_lc$S2mat)))
  431. ) %>%
  432. pivot_longer(s1:s2, names_to = 'Dec', values_to = 'Count') %>%
  433. mutate(group = c(rep('CCD',144), rep('NT',144))) %>%
  434. group_by(type, Dec, cond, group) %>%
  435. mutate(type = factor(type, levels = c('x0.5','Both','x2')),
  436. Count = Count / sum(Count)) %>%
  437. ggplot(aes(conf, Count, fill = group)) +
  438. geom_col(position = 'dodge') +
  439. labs(x = 'Confidence Bin (1=low, 6=high)', y = 'Normalised Count') +
  440. scale_x_continuous(breaks = c(6,1,-1,-6)) +
  441. facet_wrap(~type) +
  442. theme_bw(base_size = 24) +
  443. theme(legend.position = c(0.10,0.8), legend.background = element_blank(),
  444. legend.title = element_blank(), panel.grid = element_blank())
  446. conf_spread
  448. # ── 6. HIERARCHICAL META-D' FITTING ---------------
  449. # Metacognitive efficiency (meta-d'/d') estimated using the hierarchical Bayesian
  450. # method (HMeta-D; Fleming 2017). Response matrices are built with HMetaGroupPrep()
  451. # then passed to metad_group(), which is slow and should only be run once.
  453. # --- 6a. Prepare response count matrices ----
  454. LAB_CCD_metaDat <- HMetaGroupPrep(checkBothRDK, 'CCD (MRI)')
  455. ONL_CCD_metaDat <- HMetaGroupPrep(checkBothRDK, 'CCD (Online)')
  456. ONL_CCD_metaDat_exc <- HMetaGroupPrep(checkBothRDK_w_exc, 'CCD (Online)')
  458. ONL_NT_metaDat <- HMetaGroupPrep(checkBothRDK, 'NT (Online)')
  459. LAB_NT_metaDat <- HMetaGroupPrep(checkBothRDK, 'NT (MRI)')
  461. # VR: three presentation types × with/without sensitivity exclusions
  462. exc_ids <- c('4003', '4004', '4007', '4008')
  463. VR_CCD_metaDat <- HMetaGroupPrep_vr(data_vr %>% filter(BinocularTrial == 1), 'CCD')
  464. VR_CCD_metaDat_L <- HMetaGroupPrep_vr(data_vr %>% filter(LateralizedTrial == 1), 'CCD')
  465. VR_CCD_metaDat_M <- HMetaGroupPrep_vr(data_vr %>% filter(MonocularTrial == 1), 'CCD')
  466. VR_CCD_metaDat_exc <- HMetaGroupPrep_vr(data_vr %>% filter(BinocularTrial == 1, !PID %in% exc_ids), 'CCD')
  467. VR_CCD_metaDat_L_exc <- HMetaGroupPrep_vr(data_vr %>% filter(LateralizedTrial == 1, !PID %in% exc_ids), 'CCD')
  468. VR_CCD_metaDat_M_exc <- HMetaGroupPrep_vr(data_vr %>% filter(MonocularTrial == 1, !PID %in% exc_ids), 'CCD')
  470. VR_NT_metaDat <- HMetaGroupPrep_vr(data_vr %>% filter(BinocularTrial == 1), 'NT')
  471. VR_NT_metaDat_L <- HMetaGroupPrep_vr(data_vr %>% filter(LateralizedTrial == 1), 'NT')
  472. VR_NT_metaDat_M <- HMetaGroupPrep_vr(data_vr %>% filter(MonocularTrial == 1), 'NT')
  474. # --- 6b. Fit HMeta-D (run once; results saved to disk) --------
  475. ## NOTE: outputs are stochastic. Expect minor numeric deviation from published work.
  477. do_mratio <- 0
  478. if(do_mratio == 1){
  480. source('HMeta-d-master/R/Function_metad_group.R')
  482. outputOnline_NT <- metad_group(nR_S1 = list(ONL_NT_metaDat[[1]]), nR_S2 = list(ONL_NT_metaDat[[2]]));
  483. outputLab_NT <- metad_group(nR_S1 = list(LAB_NT_metaDat[[1]]), nR_S2 = list(LAB_NT_metaDat[[2]]));
  485. outputOnline_CCD <- metad_group(nR_S1 = list(ONL_CCD_metaDat[[1]]), nR_S2 = list(ONL_CCD_metaDat[[2]]));
  486. outputLab_CCD <- metad_group(nR_S1 = list(LAB_CCD_metaDat[[1]]), nR_S2 = list(LAB_CCD_metaDat[[2]]));
  488. outputVR_NT <- metad_group(nR_S1 = list(VR_NT_metaDat[[1]]), nR_S2 = list(VR_NT_metaDat[[2]]));
  489. outputVR_NT_L <- metad_group(nR_S1 = list(VR_NT_metaDat_L[[1]]), nR_S2 = list(VR_NT_metaDat_L[[2]]));
  490. outputVR_NT_M <- metad_group(nR_S1 = list(VR_NT_metaDat_M[[1]]), nR_S2 = list(VR_NT_metaDat_M[[2]]));
  492. outputVR_CCD <- metad_group(nR_S1 = list(VR_CCD_metaDat[[1]]), nR_S2 = list(VR_CCD_metaDat[[2]]));
  493. outputVR_CCD_L <- metad_group(nR_S1 = list(VR_CCD_metaDat_L[[1]]),nR_S2 = list(VR_CCD_metaDat_L[[2]]));
  494. outputVR_CCD_M <- metad_group(nR_S1 = list(VR_CCD_metaDat_M[[1]]),nR_S2 = list(VR_CCD_metaDat_M[[2]]));
  496. outputVR_CCD_exc <- metad_group(nR_S1 = list(VR_CCD_metaDat_exc[[1]]), nR_S2 = list(VR_CCD_metaDat_exc[[2]]));
  497. outputVR_CCD_L_exc <- metad_group(nR_S1 = list(VR_CCD_metaDat_L_exc[[1]]), nR_S2 = list(VR_CCD_metaDat_L_exc[[2]]));
  498. outputVR_CCD_M_exc <- metad_group(nR_S1 = list(VR_CCD_metaDat_M_exc[[1]]), nR_S2 = list(VR_CCD_metaDat_M_exc[[2]]));
  500. # Save the output
  501. saveRDS(outputLab_CCD, 'HMetaDFit/Mratio/LAB_CCD.rdata')
  502. saveRDS(outputOnline_CCD, 'HMetaDFit/Mratio/ONL_CCD.rdata')
  503. saveRDS(outputOnline_CCD_exc, 'HMetaDFit/Mratio/ONL_CCD_exc.rdata')
  504. saveRDS(outputVR_CCD, 'HMetaDFit/Mratio/VR_CCD.rdata')
  505. saveRDS(outputVR_CCD_L, 'HMetaDFit/Mratio/VR_CCD_L.rdata')
  506. saveRDS(outputVR_CCD_M, 'HMetaDFit/Mratio/VR_CCD_M.rdata')
  507. saveRDS(outputVR_CCD_exc, 'HMetaDFit/Mratio/VR_CCD_exc.rdata')
  508. saveRDS(outputVR_CCD_L_exc, 'HMetaDFit/Mratio/VR_CCD_L_exc.rdata')
  509. saveRDS(outputVR_CCD_M_exc, 'HMetaDFit/Mratio/VR_CCD_M_exc.rdata')
  510. saveRDS(outputLab_NT, 'HMetaDFit/Mratio/LAB_NT.rdata')
  511. saveRDS(outputOnline_NT, 'HMetaDFit/Mratio/ONL_NT.rdata')
  512. saveRDS(outputVR_NT, 'HMetaDFit/Mratio/VR_NT.rdata')
  513. saveRDS(outputVR_NT_L, 'HMetaDFit/Mratio/VR_NT_L.rdata')
  514. saveRDS(outputVR_NT_M, 'HMetaDFit/Mratio/VR_NT_M.rdata')
  516. }
  518. # --- 6c. Load pre-fitted models ----
  519. outputLAB_CCD <- readRDS('HMetaDFit/Mratio/LAB_CCD.rdata')
  520. outputONL_CCD <- readRDS('HMetaDFit/Mratio/ONL_CCD.rdata')
  521. outputVR_CCD <- readRDS('HMetaDFit/Mratio/VR_CCD.rdata')
  522. outputVR_CCD_L <- readRDS('HMetaDFit/Mratio/VR_CCD_L.rdata')
  523. outputVR_CCD_M <- readRDS('HMetaDFit/Mratio/VR_CCD_M.rdata')
  524. outputLAB_NT <- readRDS('HMetaDFit/Mratio/LAB_NT.rdata')
  525. outputONL_NT <- readRDS('HMetaDFit/Mratio/ONL_NT.rdata')
  526. outputVR_NT <- readRDS('HMetaDFit/Mratio/VR_NT.rdata')
  527. outputVR_NT_L <- readRDS('HMetaDFit/Mratio/VR_NT_L.rdata')
  528. outputVR_NT_M <- readRDS('HMetaDFit/Mratio/VR_NT_M.rdata')
  530. # --- 6d. Extract cleaned MCMC samples and posterior summaries ----
  531. Results_LAB_CCD <- HMeta_post_clean(outputLAB_CCD, 'CCD (LAB)', F)
  532. Results_ONL_CCD <- HMeta_post_clean(outputONL_CCD, 'CCD (ONL)', F)
  533. Results_ONL_CCD_exc <- HMeta_post_clean(outputONL_CCD_exc, 'CCD (ONL)', F)
  534. Results_VR_CCD <- HMeta_post_clean(outputVR_CCD, 'CCD (VR)', F)
  535. Results_VR_CCD_L <- HMeta_post_clean(outputVR_CCD_L, 'CCD (VR L)',F)
  536. Results_VR_CCD_M <- HMeta_post_clean(outputVR_CCD_M, 'CCD (VR M)',F)
  537. Results_ONL_NT <- HMeta_post_clean(outputONL_NT, 'NT (ONL)', F)
  538. Results_LAB_NT <- HMeta_post_clean(outputLAB_NT, 'NT (LAB)', F)
  539. Results_VR_NT <- HMeta_post_clean(outputVR_NT, 'NT (VR)', F)
  540. Results_VR_NT_L <- HMeta_post_clean(outputVR_NT_L, 'NT (VR L)', F)
  541. Results_VR_NT_M <- HMeta_post_clean(outputVR_NT_M, 'NT (VR M)', F)
  543. # Stack MCMC draws for Experiments 1 & 2
  544. mcmcGroup <- rbind(Results_LAB_CCD[[1]], Results_ONL_CCD[[1]],
  545. Results_ONL_NT[[1]], Results_LAB_NT[[1]])
  547. # Stack MCMC draws for Experiment 3 (all three presentation types)
  548. mcmcGroup_vr <- rbind(Results_VR_CCD[[1]], Results_VR_NT[[1]],
  549. Results_VR_CCD_L[[1]], Results_VR_NT_L[[1]],
  550. Results_VR_CCD_M[[1]], Results_VR_NT_M[[1]])
  552. statGroup <- data.frame(
  553. CCD_EXP1 = c(unlist(extract_stat(Results_ONL_CCD, "mu_logMratio")),
  554. unlist(extract_stat(Results_ONL_CCD, "sigma_logMratio"))),
  555. NT_EXP1 = c(unlist(extract_stat(Results_ONL_NT, "mu_logMratio")),
  556. unlist(extract_stat(Results_ONL_NT, "sigma_logMratio"))),
  557. CCD_EXP2 = c(unlist(extract_stat(Results_LAB_CCD, "mu_logMratio")),
  558. unlist(extract_stat(Results_LAB_CCD, "sigma_logMratio"))),
  559. NT_EXP2 = c(unlist(extract_stat(Results_LAB_NT, "mu_logMratio")),
  560. unlist(extract_stat(Results_LAB_NT, "sigma_logMratio"))),
  561. Variable = c('mu', 'lower mu', 'upper mu', 'sigma', 'lower sigma', 'upper sigma')
  562. )
  563. statGroup
  565. # Print VR group means for quick reference
  566. cat("VR CCD mu_Mratio:", exp(Results_VR_CCD[[3]]$mean[Results_VR_CCD[[3]]$name == "mu_logMratio"]), "\n")
  567. cat("VR NT mu_Mratio:", exp(Results_VR_NT[[3]]$mean[Results_VR_NT[[3]]$name == "mu_logMratio"]), "\n")
  568. cat("VR CCD lat mu_Mratio:",exp(Results_VR_CCD_L[[3]]$mean[Results_VR_CCD_L[[3]]$name == "mu_logMratio"]),"\n")
  569. cat("VR NT lat mu_Mratio:",exp(Results_VR_NT_L[[3]]$mean[Results_VR_NT_L[[3]]$name == "mu_logMratio"]),"\n")
  570. cat("VR CCD mon mu_Mratio:",exp(Results_VR_CCD_M[[3]]$mean[Results_VR_CCD_M[[3]]$name == "mu_logMratio"]),"\n")
  571. cat("VR NT mon mu_Mratio:",exp(Results_VR_NT_M[[3]]$mean[Results_VR_NT_M[[3]]$name == "mu_logMratio"]),"\n")
  573. # --- 6e. Posterior distribution plots (mu_logMratio) -----------
  574. # Shared theme for HMeta-D histogram panels
  575. hmeta_theme <- list(
  576. coord_cartesian(xlim = c(0, 1.5)),
  577. scale_y_continuous(expand = c(0, 0)),
  578. ggdist::theme_tidybayes(),
  579. theme(panel.background = element_rect(fill = 'white'),
  580. plot.margin = margin(0.5, 0.5, 0.5, 0.5, 'cm'))
  581. )
  583. # Experiments 1 & 2
  584. p1 <- mcmcGroup %>%
  585. filter(Parameter == "mu_logMratio") %>%
  586. ggplot(aes(exp(value), fill = group)) +
  587. geom_histogram(binwidth = 0.03, alpha = 0.9, position = 'dodge') +
  588. scale_fill_manual(values = colpal) +
  589. labs(y = "Sample Count", x = expression(paste(mu["ratio"]))) +
  590. hmeta_theme +
  591. theme(legend.position = c(0.15, 0.70), legend.background = element_blank(),
  592. legend.title = element_blank())
  594. # Experiment 3: binocular, lateralized, monocular panels
  595. p1b <- make_vr_hist(c('CCD (VR)', 'NT (VR)'))
  596. p1c <- make_vr_hist(c('CCD (VR L)', 'NT (VR L)'))
  597. p1d <- make_vr_hist(c('CCD (VR M)', 'NT (VR M)'))
  599. p1; p1b; p1c; p1d
  601. # ── 7. INDIVIDUAL-LEVEL META-D' ESTIMATES -----------------------------
  602. # Compute individual d', meta-d', and M-ratio values using make_results()
  603. # and visualise spread across participants and presentation types.
  605. # NOTE: make_results() is slow; gate with do == 1 to re-run
  607. do_meta_d_d <- 0
  608. if (do_meta_d_d == 1) {
  610. source('HMeta-d-master/R/Function_metad_group.R')
  612. ccd_onl <- make_results(ONL_CCD_metaDat, 'CCD', 'ONL')
  613. ccd_lab <- make_results(LAB_CCD_metaDat, 'CCD', 'LAB')
  614. ccd_b <- make_results(VR_CCD_metaDat, 'CCD', 'BINO')
  615. ccd_l <- make_results(VR_CCD_metaDat_L, 'CCD', 'LAT')
  616. ccd_m <- make_results(VR_CCD_metaDat_M, 'CCD', 'MONO')
  618. nt_onl <- make_results(ONL_NT_metaDat, 'NT', 'ONL')
  619. nt_lab <- make_results(LAB_NT_metaDat, 'NT', 'LAB')
  620. nt_b <- make_results(VR_NT_metaDat, 'NT', 'BINO')
  621. nt_l <- make_results(VR_NT_metaDat_L, 'NT', 'LAT')
  622. nt_m <- make_results(VR_NT_metaDat_M, 'NT', 'MONO')
  624. meta_d_df <- rbind(ccd_onl, ccd_lab, ccd_b, ccd_l, ccd_m,
  625. nt_onl, nt_lab, nt_b, nt_l, nt_m)
  627. saveRDS(meta_d_df, 'HMetaDFit/Meta_d_d/meta_d_df.rdata')
  629. }
  631. meta_d_df <- readRDS('HMetaDFit/Meta_d_d/meta_d_df.rdata')
  633. # --- 7a. LM tests of group differences ---------
  634. # Online and Lab for M-ratio, meta-d', d'
  635. for (meas in c('Mratio', 'meta_d', 'd1')) {
  636. for (env in c('ONL', 'LAB', 'BINO')) {
  637. cat("\n", meas, env, "\n")
  638. print(summary(lm(mean ~ group,
  639. data = meta_d_df %>%
  640. filter(str_detect(name, meas),
  641. !name %in% c('mu_logMratio', 'sigma_logMratio', 'mu_meta_d', 'mu_d1'),
  642. type == env))))
  643. }
  644. }
  646. # --- 7b. Group-level summary plots -----
  647. d_prime_plot_1 <- ggplot(meta_d_df %>% filter(str_detect(name, "d1"), !name %in% c('mu_d1'), type %in% c('ONL','LAB')),
  648. aes(type, mean, fill = group))
  649. meta_d_plot_1 <- ggplot(meta_d_df %>% filter(str_detect(name, "meta_d"), !name %in% c('mu_meta_d'), type %in% c('ONL','LAB')),
  650. aes(type, mean, fill = group))
  651. m_ratio_plot_1 <- ggplot(meta_d_df %>% filter(str_detect(name, "Mratio"), !name %in% c('mu_logMratio', 'sigma_logMratio'), type %in% c('ONL','LAB')),
  652. aes(type, mean, fill = group))
  654. (d_prime_plot_1 | meta_d_plot_1 | m_ratio_plot_1) &
  655. stat_summary(shape = 21, size = 1.2) &
  656. coord_cartesian(ylim = c(0, 3)) &
  657. scale_y_continuous(expand = c(0,0), breaks = seq(0, 3, 0.5)) &
  658. scale_fill_manual(values = colpal) &
  659. theme_bw(base_size = 18) &
  660. theme(legend.position = 'none', panel.grid = element_blank(), axis.title = element_blank())
  662. d_prime_plot <- ggplot(meta_d_df %>% filter(str_detect(name, "d1"), !str_detect(name, "mu_d1"), type %in% c('BINO','MONO','LAT')),
  663. aes(type, abs(mean), fill = group))
  664. meta_d_plot <- ggplot(meta_d_df %>% filter(str_detect(name, "meta_d"), !str_detect(name, "mu_meta_d"), type %in% c('BINO','MONO','LAT')),
  665. aes(type, abs(mean), fill = group))
  666. m_ratio_plot <- ggplot(meta_d_df %>% filter(str_detect(name, "Mratio"), !str_detect(name, "mu_logMratio"), !str_detect(name, "sigma_logMratio"), type %in% c('BINO','MONO','LAT')),
  667. aes(type, abs(mean), fill = group))
  669. (d_prime_plot | meta_d_plot | m_ratio_plot) &
  670. stat_summary(geom = 'line', aes(group = group)) &
  671. stat_summary(shape = 21, size = 1.2) &
  672. coord_cartesian(ylim = c(0, 2)) &
  673. scale_y_continuous(expand = c(0,0), breaks = seq(0, 2, 0.5)) &
  674. scale_fill_manual(values = colpal_vr[1:2]) &
  675. theme_bw(base_size = 18) &
  676. theme(legend.position = 'none', panel.grid = element_blank(), axis.title = element_blank())
  678. # ── 8. MODEL RECOVERY & POSTERIOR PREDICTIVE CHECK ---------------------
  679. # Validate the HMeta-D model by recovering the fitted response counts from
  680. # the posterior and correlating them with the empirical counts.
  682. ## NOTE: outcomes are stochastic. Expect some deviation from published
  683. ## results if running again.
  685. # --- 8a. Run recovery -----
  687. do_recovery <- 0
  688. if(do_recovery==1){
  690. source('HMeta-d-master/R/Function_metad_group.R')
  692. fit_CCD_lab_rec <- run_recovery(LAB_CCD_metaDat, 'CCD', 'LAB')
  693. fit_NT_lab_rec <- run_recovery(LAB_NT_metaDat, 'NT', 'LAB')
  694. fit_CCD_onl_rec <- run_recovery(ONL_CCD_metaDat, 'CCD', 'ONL')
  695. fit_NT_onl_rec <- run_recovery(ONL_NT_metaDat, 'NT', 'ONL')
  696. fit_CCD_bin_rec <- run_recovery(VR_CCD_metaDat, 'CCD', 'BINO')
  697. fit_NT_bin_rec <- run_recovery(VR_NT_metaDat, 'NT', 'BINO')
  698. fit_CCD_lat_rec <- run_recovery(VR_CCD_metaDat_L,'CCD', 'LAT')
  699. fit_NT_lat_rec <- run_recovery(VR_NT_metaDat_L, 'NT', 'LAT')
  700. fit_CCD_mon_rec <- run_recovery(VR_CCD_metaDat_M,'CCD', 'MONO')
  701. fit_NT_mon_rec <- run_recovery(VR_NT_metaDat_M, 'NT', 'MONO')
  703. recovered_dfs <- rbind(fit_CCD_lab_rec[[1]], fit_NT_lab_rec[[1]],
  704. fit_CCD_onl_rec[[1]], fit_NT_onl_rec[[1]],
  705. fit_CCD_bin_rec[[1]], fit_NT_bin_rec[[1]]) %>%
  706. mutate(type = case_when(type == 'ONL' ~ 'EXP1',
  707. type == 'LAB' ~ 'EXP2',
  708. type == 'BINO' ~ 'EXP3',
  709. .default = type))
  711. saveRDS(recovered_dfs, 'HMetaDFit/Recovery/recovered_dfs.rdata')
  713. }
  715. recovered_dfs <- readRDS('HMetaDFit/Recovery/recovered_dfs.rdata')
  717. # Correlation of real vs recovered M-ratio per participant
  718. ggplot(recovered_dfs %>% filter(str_detect(name,'Mratio'), !str_detect(name,'log')),
  719. aes(mean_real, mean_rec)) +
  720. geom_point(alpha = 0.7) +
  721. geom_smooth(method = 'lm', colour = 'black', se = FALSE) +
  722. ggpubr::stat_cor(method = 'spearman', label.x.npc = 0.25, label.y.npc = 0,
  723. colour = 'firebrick', size = 5) +
  724. facet_wrap(group ~ type, scales = 'free', nrow = 2) +
  725. labs(x = 'Real', y = 'Recovered') +
  726. theme_minimal(base_size = 18) + theme(panel.grid = element_blank())
  728. # --- 8b. Posterior predictive check -----------
  729. # Compare recovered and real response count distributions across confidence bins
  730. post_pred_check <- data.frame(
  731. conf = rep(c(rep(c(6:1,-1:-6), 4), rep(c(2:1,-1:-2), 2)), 2),
  732. cond = rep(c(rep('EXP2',24), rep('EXP1',24), rep('EXP3',8)), 2),
  733. type = c(rep('Rec',56), rep('Real',56)),
  734. group = rep(c(rep(c(rep('CCD',12), rep('NT',12)), 2), c(rep('CCD',4), rep('NT',4))), 2),
  735. s1 = c(rowSums(fit_CCD_lab_rec[[2]][[1]]), rowSums(fit_NT_lab_rec[[2]][[1]]),
  736. rowSums(fit_CCD_onl_rec[[2]][[1]]), rowSums(fit_NT_onl_rec[[2]][[1]]),
  737. rowSums(fit_CCD_bin_rec[[2]][[1]]), rowSums(fit_NT_bin_rec[[2]][[1]]),
  738. rowSums(LAB_CCD_metaDat[[1]]), rowSums(LAB_NT_metaDat[[1]]),
  739. rowSums(ONL_CCD_metaDat[[1]]), rowSums(ONL_NT_metaDat[[1]]),
  740. rowSums(VR_CCD_metaDat[[1]]), rowSums(VR_NT_metaDat[[1]])),
  741. s2 = c(rowSums(fit_CCD_lab_rec[[2]][[2]]), rowSums(fit_NT_lab_rec[[2]][[2]]),
  742. rowSums(fit_CCD_onl_rec[[2]][[2]]), rowSums(fit_NT_onl_rec[[2]][[2]]),
  743. rowSums(fit_CCD_bin_rec[[2]][[2]]), rowSums(fit_NT_bin_rec[[2]][[2]]),
  744. rowSums(LAB_CCD_metaDat[[2]]), rowSums(LAB_NT_metaDat[[2]]),
  745. rowSums(ONL_CCD_metaDat[[2]]), rowSums(ONL_NT_metaDat[[2]]),
  746. rowSums(VR_CCD_metaDat[[2]]), rowSums(VR_NT_metaDat[[2]]))
  747. ) %>%
  748. pivot_longer(s1:s2, names_to = 'Dec', values_to = 'Count') %>%
  749. group_by(type, Dec, cond, group) %>%
  750. mutate(Count = Count / sum(Count))
  752. # Compute per-bin Pearson r between real and recovered counts for each experiment
  753. bin_correlations <- bind_rows(
  754. bind_rows(fit_to_count(fit_CCD_onl_rec[[2]][[1]],"Rec","EXP1","s1"),
  755. fit_to_count(ONL_CCD_metaDat[[1]], "Real","EXP1","s1"),
  756. fit_to_count(fit_NT_onl_rec[[2]][[2]], "Rec","EXP1","s2"),
  757. fit_to_count(ONL_NT_metaDat[[2]], "Real","EXP1","s2")) %>% mutate(exp='EXP1'),
  758. bind_rows(fit_to_count(fit_CCD_lab_rec[[2]][[1]],"Rec","EXP2","s1"),
  759. fit_to_count(LAB_CCD_metaDat[[1]], "Real","EXP2","s1"),
  760. fit_to_count(fit_NT_lab_rec[[2]][[2]], "Rec","EXP2","s2"),
  761. fit_to_count(LAB_NT_metaDat[[2]], "Real","EXP2","s2")) %>% mutate(exp='EXP2'),
  762. bind_rows(fit_to_count(fit_CCD_bin_rec[[2]][[1]],"Rec","EXP3","s1"),
  763. fit_to_count(VR_CCD_metaDat[[1]], "Real","EXP3","s1"),
  764. fit_to_count(fit_NT_bin_rec[[2]][[2]], "Rec","EXP3","s2"),
  765. fit_to_count(VR_NT_metaDat[[2]], "Real","EXP3","s2")) %>% mutate(exp='EXP3')
  766. ) %>%
  767. pivot_wider(names_from = Type, values_from = Count) %>%
  768. group_by(Bin, exp) %>%
  769. summarise(test_results = list(broom::tidy(cor.test(Real, Rec, method="pearson"))),
  770. .groups = 'drop') %>%
  771. unnest(test_results) %>%
  772. dplyr::select(Bin, exp, pearson_r = estimate, p_value = p.value,
  773. conf_low = conf.low, conf_high = conf.high) %>%
  774. mutate(conf_level = case_when(
  775. exp == 'EXP3' & Bin == 1 ~ 2, exp == 'EXP3' & Bin == 2 ~ 1,
  776. exp == 'EXP3' & Bin == 3 ~ -1, exp == 'EXP3' & Bin == 4 ~ -2,
  777. exp %in% c('EXP1','EXP2') & Bin %in% 1:6 ~ 7 - Bin,
  778. exp %in% c('EXP1','EXP2') & Bin %in% 7:12 ~ 6 - Bin,
  779. .default = as.numeric(Bin)),
  780. conf_level = as.numeric(conf_level))
  782. bin_correlations %>%
  783. filter(exp=='EXP3')
  785. post_pred_check %>%
  786. pivot_wider(names_from = type, values_from = Count) %>%
  788. ggplot(aes(x = Real, y = Rec, color = factor(conf))) +
  789. geom_point(alpha = 0.5) +
  790. geom_abline(intercept = 0, slope = 1, linetype = "dashed") + # Identity line
  791. facet_wrap(~cond) +
  792. labs(x = "Real",
  793. y = "Rec.")+
  794. theme_bw(base_size = 24)+
  795. theme(legend.position = 'top',
  796. legend.background = element_blank(),
  797. legend.title = element_blank(),
  798. panel.grid = element_blank(),
  799. strip.background.x = element_blank())
  801. # Plot raw distribution comparison + per-bin r (stacked)
  802. post_pred_check_plot <- ggplot(post_pred_check, aes(conf, Count, fill = type)) +
  803. geom_col(position = 'dodge') +
  804. scale_fill_manual(values = c('#2D3142','#F6AE2D')) +
  805. scale_x_continuous(breaks = c(6,1,-1,-6)) +
  806. scale_y_continuous(breaks = seq(0.1,0.5,0.1), expand = c(0,0)) +
  807. facet_wrap(group ~ cond, scales = 'free') +
  808. labs(y = 'Normalised Count') +
  809. theme_bw(base_size = 20) +
  810. theme(legend.position = 'none', panel.grid = element_blank(),
  811. strip.background.x = element_blank(), axis.title.x = element_blank())
  813. cor_plot_post_pred <- ggplot(bin_correlations, aes(conf_level, pearson_r)) +
  814. geom_col(fill = "black") +
  815. geom_hline(yintercept = 0.8, color = "red") +
  816. labs(x = "Confidence Bin (1=low, 6=high)", y = "Pearson r") +
  817. facet_wrap(~exp, scales = 'free') +
  818. scale_x_continuous(breaks = c(6,1,-1,-6)) +
  819. scale_y_continuous(breaks = c(0, 0.5, 1)) +
  820. theme_bw(base_size = 20) +
  821. theme(panel.grid = element_blank(), strip.background.x = element_blank(),
  822. strip.text.x = element_blank())
  824. (post_pred_check_plot / cor_plot_post_pred) & plot_layout(heights = c(5,1))
  826. # ── 9. PERMUTATION TESTS ON META-D' GROUP DIFFERENCES -------------------
  827. # We use a permutation procedure to test whether the observed NT − CCD difference
  828. # in mu_logMratio (and d', meta-d') exceeds the null distribution
  829. # generated by label-shuffling.
  831. # --- 9a. Extract observed group differences ----
  832. # Helper: pull exp(mu_logMratio) from a make_results() output dataframe
  833. pull_mu <- function(df, param = 'mu_logMratio') df %>% filter(name == param) %>% mutate(mean = exp(mean))
  834. pull_meta_d <- function(df, param = 'mu_meta_d') df %>% filter(name == param) %>% mutate(mean = abs(mean))
  835. pull_d1 <- function(df, param = 'mu_d1') df %>% filter(name == param) %>% mutate(mean = abs(mean))
  837. mu_log_exp1_NT <- pull_mu(Results_ONL_NT$Fit); mu_log_exp1_CCD <- pull_mu(Results_ONL_CCD$Fit)
  838. mu_log_exp2_NT <- pull_mu(Results_LAB_NT$Fit); mu_log_exp2_CCD <- pull_mu(Results_LAB_CCD$Fit)
  839. mu_log_exp3_NT <- pull_mu(Results_VR_NT$Fit); mu_log_exp3_CCD <- pull_mu(Results_VR_CCD$Fit)
  840. mu_log_exp3_NT_L <- pull_mu(Results_VR_NT_L$Fit); mu_log_exp3_CCD_L <- pull_mu(Results_VR_CCD_L$Fit)
  841. mu_log_exp3_NT_M <- pull_mu(Results_VR_NT_M$Fit); mu_log_exp3_CCD_M <- pull_mu(Results_VR_CCD_M$Fit)
  843. meta_d_log_exp1_NT <- pull_meta_d(nt_onl);meta_d_log_exp1_CCD <- pull_meta_d(ccd_onl);
  844. meta_d_log_exp2_NT <- pull_meta_d(nt_lab);meta_d_log_exp2_CCD <- pull_meta_d(ccd_lab);
  845. meta_d_log_exp3_NT <- pull_meta_d(nt_b); meta_d_log_exp3_CCD <- pull_meta_d(ccd_b);
  846. meta_d_log_exp3_NT_L <- pull_meta_d(nt_l); meta_d_log_exp3_CCD_L <- pull_meta_d(ccd_l);
  847. meta_d_log_exp3_NT_M <- pull_meta_d(nt_m); meta_d_log_exp3_CCD_M <- pull_meta_d(ccd_m);
  849. d1_log_exp1_NT <- pull_d1(nt_onl);d1_log_exp1_CCD <- pull_d1(ccd_onl);
  850. d1_log_exp2_NT <- pull_d1(nt_lab);d1_log_exp2_CCD <- pull_d1(ccd_lab);
  851. d1_log_exp3_NT <- pull_d1(nt_b); d1_log_exp3_CCD <- pull_d1(ccd_b);
  852. d1_log_exp3_NT_L <- pull_d1(nt_l); d1_log_exp3_CCD_L <- pull_d1(ccd_l);
  853. d1_log_exp3_NT_M <- pull_d1(nt_m); d1_log_exp3_CCD_M <- pull_d1(ccd_m);
  855. # Compute NT − CCD differences for M-ratio, d', and meta-d' in each context
  856. diffExp1 <- as.numeric(mu_log_exp1_NT[2] - mu_log_exp1_CCD[2])
  857. diffExp2 <- as.numeric(mu_log_exp2_NT[2] - mu_log_exp2_CCD[2])
  858. diffExp3 <- as.numeric(mu_log_exp3_NT[2] - mu_log_exp3_CCD[2])
  859. diffExp3_L <- as.numeric(mu_log_exp3_NT_L[2] - mu_log_exp3_CCD_L[2])
  860. diffExp3_M <- as.numeric(mu_log_exp3_NT_M[2] - mu_log_exp3_CCD_M[2])
  862. diffExp1_d1 <- as.numeric(d1_log_exp1_NT[2] - d1_log_exp1_CCD[2])
  863. diffExp2_d1 <- as.numeric(d1_log_exp2_NT[2] - d1_log_exp2_CCD[2])
  864. diffExp3_d1 <- as.numeric(d1_log_exp3_NT[2] - d1_log_exp3_CCD[2])
  865. diffExp3_L_d1 <- as.numeric(d1_log_exp3_NT_L[2] - d1_log_exp3_CCD_L[2])
  866. diffExp3_M_d1 <- as.numeric(d1_log_exp3_NT_M[2] - d1_log_exp3_CCD_M[2])
  868. diffExp1_metad <- as.numeric(meta_d_log_exp1_NT[2] - meta_d_log_exp1_CCD[2])
  869. diffExp2_metad <- as.numeric(meta_d_log_exp2_NT[2] - meta_d_log_exp2_CCD[2])
  870. diffExp3_metad <- as.numeric(meta_d_log_exp3_NT[2] - meta_d_log_exp3_CCD[2])
  871. diffExp3_L_metad <- as.numeric(meta_d_log_exp3_NT_L[2] - meta_d_log_exp3_CCD_L[2])
  872. diffExp3_M_metad <- as.numeric(meta_d_log_exp3_NT_M[2] - meta_d_log_exp3_CCD_M[2])
  874. cat("mu difference Lab:", diffExp1, "\n")
  875. cat("mu difference Onl:", diffExp2, "\n")
  876. cat("mu difference VR:", diffExp3, "\n")
  877. cat("mu difference VR-L:", diffExp3_L, "\n")
  878. cat("mu difference VR-M:", diffExp3_M, "\n")
  880. # --- 9b. Prepare data for permutation ----
  881. data_for_permute_exp1_CCD <- checkBothRDK %>% filter(group=='CCD (Online)') %>% dplyr::select(ID, Trial, conf, dotDirection, type, kmed, referenceSelection, group)
  882. data_for_permute_exp1_NT <- checkBothRDK %>% filter(group=='NT (Online)') %>% dplyr::select(ID, Trial, conf, dotDirection, type, kmed, referenceSelection, group)
  884. data_for_permute_exp2_CCD <- checkBothRDK %>% filter(group=='CCD (MRI)') %>% dplyr::select(ID, Trial, conf, dotDirection, type, kmed, referenceSelection, group)
  885. data_for_permute_exp2_NT <- checkBothRDK %>% filter(group=='NT (MRI)') %>% dplyr::select(ID, Trial, conf, dotDirection, type, kmed, referenceSelection, group)
  887. data_for_permute_exp3_CCD <- data_vr %>% filter(BinocularTrial == 1, Group=='CCD')
  888. data_for_permute_exp3_NT <- data_vr %>% filter(BinocularTrial == 1, Group=='NT')
  889. data_for_permute_exp3_CCD_L <- data_vr %>% filter(LateralizedTrial == 1, Group=='CCD')
  890. data_for_permute_exp3_NT_L <- data_vr %>% filter(LateralizedTrial == 1, Group=='NT')
  891. data_for_permute_exp3_CCD_M <- data_vr %>% filter(MonocularTrial == 1, Group=='CCD')
  892. data_for_permute_exp3_NT_M <- data_vr %>% filter(MonocularTrial == 1, Group=='NT')
  894. # --- 9c. Run permutation (slow; only re-run when needed) ----
  896. ## NOTE: outcomes are stochastic. Expect some deviation from published
  897. ## results if running again.
  898. do_perm <- 0
  899. if(do_perm==1){
  901. source('HMeta-d-master/R/Function_metad_group.R')
  902. nreps <- 500
  904. onl_CCDvsNTp <- permute_metad_group(data_for_permute_exp1_CCD, data_for_permute_exp1_NT, cores = 5, nreps = nreps);
  905. onl_CCDvsNTp_mrat <- permute_metad_group_mrat(data_for_permute_exp1_CCD, data_for_permute_exp1_NT, cores = 10, nreps = nreps);
  907. lab_CCDvsNTp <- permute_metad_group(data_for_permute_exp2_CCD, data_for_permute_exp2_NT, cores = 10, nreps = nreps);
  908. lab_CCDvsNTp_mrat <- permute_metad_group_mrat(data_for_permute_exp2_CCD, data_for_permute_exp2_NT, cores = 10, nreps = nreps);
  910. vr_CCDvsNTp <- permute_metad_group_vr(data_for_permute_exp3_CCD, data_for_permute_exp3_NT,cores = 10, nreps = nreps);
  911. vr_CCDvsNTp_L<- permute_metad_group_vr(data_for_permute_exp3_CCD_L, data_for_permute_exp3_NT_L,cores = 10, nreps = nreps);
  912. vr_CCDvsNTp_M<- permute_metad_group_vr(data_for_permute_exp3_CCD_M, data_for_permute_exp3_NT_M,cores = 10, nreps = nreps);
  914. saveRDS(onl_CCDvsNTp, 'HMetaDFit/Permutation/CCDvsNT_ONL.rdata')
  915. saveRDS(onl_CCDvsNTp_mrat, 'HMetaDFit/Permutation/CCDvsNT_ONL_mrat.rdata')
  916. saveRDS(lab_CCDvsNTp, 'HMetaDFit/Permutation/CCDvsNT_LAB.rdata')
  917. saveRDS(lab_CCDvsNTp_mrat, 'HMetaDFit/Permutation/CCDvsNT_LAB_mrat.rdata')
  918. saveRDS(vr_CCDvsNTp, 'HMetaDFit/Permutation/CCDvsNT_VR.rdata')
  919. saveRDS(vr_CCDvsNTp_L,'HMetaDFit/Permutation/CCDvsNT_VR_L.rdata')
  920. saveRDS(vr_CCDvsNTp_M,'HMetaDFit/Permutation/CCDvsNT_VR_M.rdata')
  922. }
  924. # --- 9d. Load permutation null distributions ----
  925. onl_CCDvsNTp <- readRDS('HMetaDFit/Permutation/CCDvsNT_ONL.rdata')
  926. onl_CCDvsNTp_mrat <- readRDS('HMetaDFit/Permutation/CCDvsNT_ONL_mrat.rdata')
  927. lab_CCDvsNTp <- readRDS('HMetaDFit/Permutation/CCDvsNT_LAB.rdata')
  928. lab_CCDvsNTp_mrat <- readRDS('HMetaDFit/Permutation/CCDvsNT_LAB_mrat.rdata')
  929. vr_CCDvsNTp <- readRDS('HMetaDFit/Permutation/CCDvsNT_VR.rdata')
  930. vr_CCDvsNTp_L <- readRDS('HMetaDFit/Permutation/CCDvsNT_VR_L.rdata')
  931. vr_CCDvsNTp_M <- readRDS('HMetaDFit/Permutation/CCDvsNT_VR_M.rdata')
  933. # --- 9e. Effect sizes and permutation p-values ------
  935. # perm_effect_stats() returns the p-value and U3 superiority index
  936. cat("\n--- Online (Exp 1) ---\n")
  937. perm_effect_stats(round(diffExp1, 2), onl_CCDvsNTp_mrat[,4])
  938. perm_effect_stats(round(diffExp1_d1, 2), onl_CCDvsNTp[,7])
  939. perm_effect_stats(round(diffExp1_metad, 2), onl_CCDvsNTp[,10])
  941. cat("\n--- Lab (Exp 2) ---\n")
  942. perm_effect_stats(round(diffExp2, 2), lab_CCDvsNTp_mrat[,4])
  943. perm_effect_stats(round(diffExp2_d1, 2), lab_CCDvsNTp[,7])
  944. perm_effect_stats(round(diffExp2_metad, 2), lab_CCDvsNTp[,10])
  946. cat("\n--- VR binocular (Exp 3) ---\n")
  947. perm_effect_stats(round(diffExp3, 2), vr_CCDvsNTp_mrat[,4])
  948. perm_effect_stats(round(diffExp3_d1, 2), vr_CCDvsNTp[,7])
  949. perm_effect_stats(round(diffExp3_meta_d, 2), vr_CCDvsNTp[,10])
  951. cat("\n--- VR lateralized (Exp 3) ---\n")
  952. perm_effect_stats(round(diffExp3_L, 2), vr_CCDvsNTp_mrat_L[,4])
  953. perm_effect_stats(round(diffExp3_L_d1, 2), vr_CCDvsNTp_L[,7])
  954. perm_effect_stats(round(diffExp3_L_metad, 2), vr_CCDvsNTp_L[,10])
  956. cat("\n--- VR monocular (Exp 3) ---\n")
  957. perm_effect_stats(round(diffExp3_M, 2), vr_CCDvsNTp_mrat_M[,4])
  958. perm_effect_stats(round(diffExp3_M_d1, 2), vr_CCDvsNTp_M[,7])
  959. perm_effect_stats(round(diffExp3_M_metad, 2), vr_CCDvsNTp_M[,10])
  961. # --- 9f. Permutation histogram plots -------
  962. # One panel per experiment, x-axis = null M-ratio difference, vertical line = observed
  963. p_p1 <- onl_CCDvsNTp %>% as.data.frame() %>% rename(Diff = 4) %>%
  964. ggplot(aes(Diff)) + geom_histogram(binwidth = 0.01) +
  965. labs(x = expression(paste("Randomised ",mu['ratio'])), y = 'Sample Count') +
  966. geom_label(x = diffExp1 + 0.07, y = 10, label = round(diffExp1,2), fill='#81D6E3', colour='white') +
  967. geom_vline(xintercept = diffExp1, color = '#81D6E3', size = 1.5) +
  968. theme_minimal() + theme(panel.grid = element_blank(), axis.text = element_text(size=24), axis.title = element_text(size=24))
  970. p_p2 <- lab_CCDvsNTp %>% as.data.frame() %>% rename(Diff = 4) %>%
  971. ggplot(aes(Diff)) + geom_histogram(binwidth = 0.01) +
  972. labs(x = expression(paste("Randomised ",mu['ratio'])), y = 'Sample Count') +
  973. geom_label(x = diffExp2 + 0.07, y = 10, label = round(diffExp2,2), fill='#475B5A', colour='white') +
  974. geom_vline(xintercept = diffExp2, color = '#475B5A', size = 1.5) +
  975. theme_minimal() + theme(panel.grid = element_blank(), axis.text = element_text(size=24), axis.title = element_text(size=24))
  977. p_p_b <- make_vr_perm_plot(vr_CCDvsNTp, diffExp3)
  978. p_p_l <- make_vr_perm_plot(vr_CCDvsNTp_L, diffExp3_L)
  979. p_p_m <- make_vr_perm_plot(vr_CCDvsNTp_M, diffExp3_M)
  981. (p_p1 | p_p2)
  982. (p_p_b | p_p_l | p_p_m)
  984. # ── 10. SUPPLEMENT: REACTION TIME ANALYSIS -----------------------------
  985. # Supplementary: RT distributions and LME tests by coherence and correctness.
  986. # Decision trials (conf == 0) and confidence trials (conf != 0) treated separately.
  988. library(ggridges)
  989. Experiment_1_both <- rbind(Experiment_1_CCD, Experiment_1_NT)
  990. Experiment_2_both <- rbind(Experiment_2_CCD, Experiment_2_NT)
  992. ## Non Confidence Trials -----
  994. ### Split by coherence ----
  996. Experiment_1_both %>%
  997. filter(type == 'main', kmed != 'kmed', ID != '5f20a850b6d5ed3d186bc3c3', conf == 0) %>%
  998. group_by(group, kmed, ID) %>%
  999. summarise(Subj_Acc = mean(correct)*100, Subj_RT = mean(RT, na.rm=TRUE), .groups='drop') %>%
  1000. group_by(group, kmed) %>%
  1001. summarise(Mean_Acc=mean(Subj_Acc), SD_Acc=sd(Subj_Acc),
  1002. Mean_RT=mean(Subj_RT), SD_RT=sd(Subj_RT),
  1003. SE_RT=SD_RT/sqrt(n()), n=n(), .groups='drop')
  1005. Experiment_2_both %>%
  1006. filter(type == 'main', kmed != 'kmed', conf == 0) %>%
  1007. group_by(group, kmed, ID) %>%
  1008. summarise(Subj_Acc = mean(correct)*100, Subj_RT = mean(RT, na.rm=TRUE), .groups='drop') %>%
  1009. group_by(group, kmed) %>%
  1010. summarise(Mean_Acc=mean(Subj_Acc), SD_Acc=sd(Subj_Acc),
  1011. Mean_RT=mean(Subj_RT), SD_RT=sd(Subj_RT),
  1012. SE_RT=SD_RT/sqrt(n()), n=n(), .groups='drop')
  1014. run_rt_lmer(Experiment_1_both, 'Computer', conf_filter = 0) # Exp 1: decision RT × coherence
  1015. run_rt_lmer(Experiment_2_both, 'MRI', conf_filter = 0) # Exp 2: decision RT × coherence
  1017. ### Split by correctness ----
  1019. Experiment_1_both %>%
  1020. filter(type == 'main', kmed != 'kmed', ID != '5f20a850b6d5ed3d186bc3c3', conf == 0) %>%
  1021. group_by(group, correct, ID) %>%
  1022. summarise(Subj_Acc = mean(correct)*100, Subj_RT = mean(RT, na.rm=TRUE), .groups='drop') %>%
  1023. group_by(group, correct) %>%
  1024. summarise(Mean_Acc=mean(Subj_Acc), SD_Acc=sd(Subj_Acc),
  1025. Mean_RT=mean(Subj_RT), SD_RT=sd(Subj_RT),
  1026. SE_RT=SD_RT/sqrt(n()), n=n(), .groups='drop')
  1028. Experiment_2_both %>%
  1029. filter(type == 'main', kmed != 'kmed', conf == 0) %>%
  1030. group_by(group, correct, ID) %>%
  1031. summarise(Subj_Acc = mean(correct)*100, Subj_RT = mean(RT, na.rm=TRUE), .groups='drop') %>%
  1032. group_by(group, correct) %>%
  1033. summarise(Mean_Acc=mean(Subj_Acc), SD_Acc=sd(Subj_Acc),
  1034. Mean_RT=mean(Subj_RT), SD_RT=sd(Subj_RT),
  1035. SE_RT=SD_RT/sqrt(n()), n=n(), .groups='drop')
  1037. lmerTest::lmer(log(RT) ~ CCD * correct + (1|ID),
  1038. data = Experiment_1_both %>% filter(type=='main', conf==0, kmed!='kmed', environment=='Computer')) %>% summary()
  1039. lmerTest::lmer(log(RT) ~ CCD * correct + (1|ID),
  1040. data = Experiment_2_both %>% filter(type=='main', conf==0, kmed!='kmed', environment=='MRI')) %>% summary()
  1042. # VR RT models by presentation type
  1043. data_vr %>%
  1044. filter(TrialType == 'Main') %>%
  1045. group_by(Group, PresentationType, HighCoherenceTrialLabel, PID) %>%
  1046. summarise(Subj_Acc = mean(Correct)*100, Subj_RT = mean(ReactionTime, na.rm=TRUE), .groups='drop') %>%
  1047. group_by(Group, PresentationType, HighCoherenceTrialLabel) %>%
  1048. summarise(Mean_Acc=mean(Subj_Acc), SD_Acc=sd(Subj_Acc),
  1049. Mean_RT=mean(Subj_RT), SD_RT=sd(Subj_RT),
  1050. SE_RT=SD_RT/sqrt(n()), n=n(), .groups='drop')
  1052. for (ptype in c('Binocular','Lateralized','Monocular')) {
  1053. cat("\n=== RT VR", ptype, "===\n")
  1054. lmerTest::lmer(log(ReactionTime) ~ CCD * HighCoherenceTrialLabel + (1|PID),
  1055. data = data_vr %>% filter(TrialType=='Main', PresentationType==ptype,
  1056. HighCoherenceTrialLabel %in% c('x0.5','x2'))) %>% summary() %>% print()
  1057. lmerTest::lmer(log(ReactionTime) ~ CCD * Correct + (1|PID),
  1058. data = data_vr %>% filter(TrialType=='Main', PresentationType==ptype,
  1059. HighCoherenceTrialLabel %in% c('x0.5','x2'))) %>% summary() %>% print()
  1060. }
  1062. # ── 11. SUPPLEMENT: CALIBRATION ASSESSMENT ---------------------------------
  1063. # Show that k_med staircase converges similarly in both groups and environments.
  1065. ## Non-VR Based ----
  1067. ID_baseK <- checkBothRDK %>%
  1068. filter(type == 'calibration', Trial > 100) %>%
  1069. group_by(ID, CCD, environment) %>%
  1070. summarise(baseK = median(coherence))
  1072. allNonVRK <- ID_baseK %>%
  1073. group_by(CCD, environment) %>%
  1074. summarise(Coherence = mean(baseK),
  1075. CoherenceSD = sd(baseK))
  1076. allNonVRK
  1077. # Test whether base k varied across computer or MRI versions of the same task
  1078. t.test(x = (ID_baseK %>% filter(CCD == 1, environment == 'Computer'))$baseK,
  1079. y = (ID_baseK %>% filter(CCD == 0, environment == 'Computer'))$baseK)
  1080. t.test(x = (ID_baseK %>% filter(CCD == 1, environment == 'MRI'))$baseK,
  1081. y = (ID_baseK %>% filter(CCD == 0, environment == 'MRI'))$baseK)
  1083. # Base K between conditions ------
  1084. coh_tog <- rbind(
  1085. ID_baseK,
  1086. data_vr %>%
  1087. filter(TrialType == 'Main') %>%
  1088. group_by(PID, Group) %>%
  1089. summarise(Coherence = mean(BaseCoherence)) %>%
  1090. mutate(environment = 'VR',
  1091. Group = ifelse(Group=='CCD', 1, 0),
  1092. PID = as.character(PID)) %>%
  1093. dplyr::select(ID = PID, CCD = Group, environment, baseK = Coherence)
  1094. )
  1096. aov_tog_ccd <- aov(baseK ~ environment, coh_tog %>% filter(CCD==1))
  1097. summary(aov_tog_ccd)
  1098. TukeyHSD(aov_tog_ccd)
  1100. aov_tog_nt <- aov(baseK ~ environment, coh_tog %>% filter(CCD==0))
  1101. summary(aov_tog_nt)
  1102. TukeyHSD(aov_tog_nt)
  1104. ## VR Based -----
  1105. allVRK <- data_vr %>%
  1106. filter(TrialType == 'Main') %>%
  1107. group_by(CCD, PresentationType) %>%
  1108. summarise(Coherence = mean(BaseCoherence),
  1109. CoherenceSD = sd(BaseCoherence))
  1111. binoK <- data_vr %>%
  1112. filter(BinocularTrial == 1, TrialType == 'Main') %>%
  1113. group_by(PID, CCD) %>%
  1114. summarise(Coherence = mean(BaseCoherence))
  1116. monoK <- data_vr %>%
  1117. filter(MonocularTrial == 1, TrialType == 'Main') %>%
  1118. group_by(PID, CCD) %>%
  1119. summarise(Coherence = mean(BaseCoherence))
  1121. latK <- data_vr %>%
  1122. filter(LateralizedTrial == 1, TrialType == 'Main') %>%
  1123. group_by(PID, CCD) %>%
  1124. summarise(Coherence = mean(BaseCoherence))
  1126. # Test whether base k varied across presentations
  1127. t.test((binoK %>% filter(CCD == 1))$Coherence, y = (binoK %>% filter(CCD == 0))$Coherence)
  1128. t.test((monoK %>% filter(CCD == 1))$Coherence, y = (monoK %>% filter(CCD == 0))$Coherence)
  1129. t.test((latK %>% filter(CCD == 1))$Coherence, y = (latK %>% filter(CCD == 0))$Coherence)
  1131. #ANOVA
  1132. coh_aov_ccd <- lmerTest::lmer(Coherence ~ PresentationType + (1|PID),
  1133. data =
  1134. data_vr %>%
  1135. filter(TrialType == 'Main', Group == 'CCD') %>%
  1136. group_by(PID, PresentationType) %>%
  1137. mutate(PresentationType = factor(PresentationType, levels = c( 'Monocular', 'Binocular', 'Lateralized'))) %>%
  1138. summarise(Coherence = mean(BaseCoherence)))
  1139. anova(coh_aov_ccd)
  1140. summary(coh_aov_ccd)
  1142. coh_aov_nt <- lmerTest::lmer(Coherence ~ PresentationType + (1|PID),
  1143. data =
  1144. data_vr %>%
  1145. filter(TrialType == 'Main', Group == 'NT') %>%
  1146. group_by(PID, PresentationType) %>%
  1147. summarise(Coherence = mean(BaseCoherence)))
  1148. anova(coh_aov_nt)
  1149. summary(coh_aov_nt)
  1151. coh_aov_comb <- lmerTest::lmer(Coherence ~ PresentationType * Group + (1|PID),
  1152. data =
  1153. data_vr %>%
  1154. filter(TrialType == 'Main') %>%
  1155. group_by(PID, PresentationType, Group) %>%
  1156. summarise(Coherence = mean(BaseCoherence)))
  1157. anova(coh_aov_comb)
  1158. summary(coh_aov_comb)
  1161. # Combine all experiments into one long format
  1162. check_coh_RDK <- rbind(
  1163. Experiment_1_CCD, Experiment_2_CCD, Experiment_1_NT, Experiment_2_NT
  1164. ) %>%
  1165. mutate(prestype = NA, vr = 0) %>%
  1166. dplyr::select(ID, Trial, coherence, group, prestype, type, vr) %>%
  1167. rbind(
  1168. data_vr %>%
  1169. rename(ID = PID, group = Group, coherence = ActiveCoherence,
  1170. prestype = PresentationType, type = TrialType) %>%
  1171. mutate(vr = 1) %>%
  1172. dplyr::select(ID, Trial, coherence, group, prestype, type, vr)
  1173. )
  1175. # Coherence trajectory plot: Experiments 1 & 2
  1176. CCDcheck2 <-
  1177. ggplot(check_coh_RDK %>% filter(type == 'calibration', vr == 0)) +
  1178. stat_summary(aes(Trial, coherence, fill = group), geom='ribbon', alpha=0.1, colour=NA) +
  1179. stat_summary(aes(Trial, coherence, colour = group), geom='line') +
  1180. scale_color_manual(values = colpal) + scale_fill_manual(values = colpal) +
  1181. labs(x = 'Trial', y = 'Coherence') +
  1182. coord_cartesian(ylim = c(0, 0.5), xlim = c(0, 118)) +
  1183. scale_x_continuous(expand = c(0,0)) +
  1184. theme_minimal(base_size = 20) +
  1185. theme(legend.position = 'top', panel.grid.minor = element_blank(), legend.title = element_blank())
  1187. bino_check <- make_vr_cal_plot('Binocular', 'Bino.')
  1188. lat_check <- make_vr_cal_plot('Lateralized','Lat.')
  1189. mono_check <- make_vr_cal_plot('Monocular', 'Mono.')
  1191. coh_vr_check_full <- ggarrange(bino_check, lat_check, mono_check, nrow = 1)
  1193. CCDcheck2; coh_vr_check_full
  1195. # Correlation: does mean confidence predict k_med? (supplementary validity check)
  1196. check_coh_lm <- checkBothRDK %>%
  1197. group_by(ID) %>%
  1198. filter(type == 'main', conf != 0) %>%
  1199. mutate(kmed_av = ifelse(kmed == 'kmed x 0.5', coherence*2, coherence/2),
  1200. conf = mean(conf)) %>%
  1201. dplyr::select(kmed_av, conf, ID, group) %>%
  1202. distinct()
  1204. ggplot(check_coh_lm, aes(conf, kmed_av, colour = group)) +
  1205. geom_point() + geom_smooth(method = 'lm') +
  1206. coord_cartesian(ylim = c(0, 0.5)) +
  1207. scale_color_manual(values = colpal) +
  1208. stat_cor(show.legend = FALSE, method = 'spearman') +
  1209. stat_cor(show.legend = FALSE, aes(colour = NULL), data = check_coh_lm,
  1210. label.y.npc = 0.50, method = 'spearman') +
  1211. theme_bw(base_size = 20) +
  1212. theme(legend.position = 'left', legend.title = element_blank())
  1215. # ── 12. SUPPLEMENT: GROUP BEHAVIOUR FIGURES --------------------------
  1217. # Experiments 1 & 2: beeswarm accuracy and confidence plots with significance bars
  1218. tidy_acc <- tidyplot(
  1219. av_cor %>% mutate(
  1220. kmed = ifelse(kmed == 'x0.5', 'x0.5', 'x2.0'),
  1221. group = dplyr::recode(group, 'CCD (MRI)'='CCD\n(MRI)', 'NT (Online)'='NT\n(Online)',
  1222. 'CCD (Online)'='CCD\n(Online)', 'NT (MRI)'='NT\n(MRI)')),
  1223. x = group, y = av_cor, colour = kmed) %>%
  1224. add_data_points_beeswarm(alpha=0.1, cex=1) %>% add_sd_errorbar() %>% add_mean_dot(size=2) %>%
  1225. adjust_size(width=75, height=50) %>% adjust_x_axis_title(title='') %>%
  1226. adjust_y_axis_title(title='p(Correct)') %>% adjust_font(fontsize=15) %>%
  1227. adjust_legend_title(title='Coh.') %>% adjust_legend_position(position='top') +
  1228. stat_compare_means(paired=TRUE, hide.ns=TRUE, size=5, show.legend=FALSE,
  1229. label='p.signif', label.y=1.04) +
  1230. coord_cartesian(clip="off")
  1232. tidy_conf <- tidyplot(
  1233. av_co %>% mutate(
  1234. kmed = ifelse(kmed == 'x0.5', 'x0.5', 'x2.0'),
  1235. group = dplyr::recode(group, 'CCD (MRI)'='CCD\n(MRI)', 'NT (Online)'='NT\n(Online)',
  1236. 'CCD (Online)'='CCD\n(Online)', 'NT (MRI)'='NT\n(MRI)'),
  1237. ID = factor(ID)),
  1238. x = group, y = av_co, colour = kmed) %>%
  1239. add_data_points_beeswarm(alpha=0.1, cex=1) %>% add_sd_errorbar() %>% add_mean_dot(size=2) %>%
  1240. adjust_size(width=75, height=50) %>% adjust_x_axis_title(title='') %>%
  1241. adjust_y_axis_title(title='Confidence') %>% adjust_font(fontsize=15) %>%
  1242. adjust_legend_title(title='Coh.') %>% adjust_legend_position(position='top') +
  1243. stat_compare_means(paired=TRUE, hide.ns=TRUE, size=5, show.legend=FALSE,
  1244. label='p.signif', label.y=102) +
  1245. coord_cartesian(clip="off")
  1247. tidy_acc | tidy_conf
  1249. tidyplot_exp3_cor <- av_cor_vr %>%
  1250. mutate(kmed = as.character(kmed),
  1251. Pres = dplyr::recode(PresentationType, 'Monocular'='Mono.', 'Lateralized'='Lat.', 'Binocular'='Bino.'),
  1252. ID = factor(ID))
  1254. tidyplot_exp3_co <- av_co_vr %>%
  1255. mutate(kmed = as.character(kmed),
  1256. Pres = dplyr::recode(PresentationType, 'Monocular'='Mono.', 'Lateralized'='Lat.', 'Binocular'='Bino.'),
  1257. ID = factor(ID))
  1259. Cor <- make_exp3_plot(tidyplot_exp3_cor, 'av_cor', 'p(Correct)', 1.04)
  1260. Conf <- make_exp3_plot(tidyplot_exp3_co, 'av_co', 'Confidence', 0.75)
  1262. Cor / Conf
  1264. # ── 13. SUPPLEMENT: INTELLIGENCE CORRELATION -------------------------------
  1265. # ICAR scores from online platforms linked to individual meta-d' estimates to
  1266. # check whether metacognitive efficiency correlates with fluid intelligence.
  1268. ## NT ------
  1270. nt_join_ICAR_mrat <- read_csv('Data/ICAR_NT_mrat.csv')
  1271. nt_join_ICAR_meta_d <- read_csv('Data/ICAR_NT_meta_d.csv')
  1272. nt_join_ICAR_d1 <- read_csv('Data/ICAR_NT_d1.csv')
  1274. cor.test(nt_join_ICAR_mrat$mean, nt_join_ICAR_mrat$ICAR)
  1275. cor.test(nt_join_ICAR_meta_d$mean, nt_join_ICAR_meta_d$ICAR)
  1276. cor.test(nt_join_ICAR_d1$mean, nt_join_ICAR_d1$ICAR)
  1278. ## CCD ------
  1280. ccd_join_ICAR_mrat <- read_csv('Data/ICAR_CCD_mrat.csv')
  1281. ccd_join_ICAR_meta_d <- read_csv('Data/ICAR_CCD_meta_d.csv')
  1282. ccd_join_ICAR_d1 <- read_csv('Data/ICAR_CCD_d1.csv')
  1284. cor.test(ccd_join_ICAR_mrat$mean, ccd_join_ICAR_mrat$ICAR)
  1285. cor.test(ccd_join_ICAR_meta_d$mean, ccd_join_ICAR_meta_d$ICAR)
  1286. cor.test(ccd_join_ICAR_d1$mean, ccd_join_ICAR_d1$ICAR)
  1288. # ── 14. SUPPLEMENT: EXCLUSION ANALYSIS -------------------------------
  1290. ### Experiment 1 ------
  1292. run_lmer(correct ~ kmed + (1|ID), lmer_data_exc %>% filter(group == 'CCD (Online)'))
  1293. run_lmer(correct ~ kmed + (1|ID), lmer_data %>% filter(group == 'NT (Online)'))
  1294. run_lmer(correct ~ kmed * group + (1|ID), lmer_data_exc %>% filter(group %in% c('NT (Online)', 'CCD (Online)')))
  1296. run_lmer(conf ~ kmed + (1|ID), lmer_data_exc %>% filter(group == 'CCD (Online)'))
  1297. run_lmer(conf ~ kmed + (1|ID), lmer_data %>% filter(group == 'NT (Online)'))
  1298. run_lmer(conf ~ kmed * group + (1|ID), lmer_data_exc %>% filter(group %in% c('NT (Online)', 'CCD (Online)')))
  1300. ### Experiment 3 ---------
  1302. run_lmer(Correct ~ kmed * group + (1|ID), data_vr_tests %>% filter(!ID %in% c('4003', '4004', '4007', '4008')))
  1303. run_lmer(conf ~ kmed * group + (1|ID), data_vr_tests %>% filter(!ID %in% c('4003', '4004', '4007', '4008')))
  1305. #### BINO ----
  1307. run_lmer(Correct ~ kmed + (1|ID), data_vr_tests %>% filter(group=='CCD', PresentationType=='Binocular', !ID %in% c('4003', '4004', '4007', '4008')))
  1308. run_lmer(Correct ~ kmed + (1|ID), data_vr_tests %>% filter(group=='NT', PresentationType=='Binocular'))
  1310. run_lmer(conf ~ kmed + (1|ID), data_vr_tests %>% filter(group=='CCD', PresentationType=='Binocular', !ID %in% c('4003', '4004', '4007', '4008')))
  1311. run_lmer(conf ~ kmed + (1|ID), data_vr_tests %>% filter(group=='NT', PresentationType=='Binocular'))
  1313. #### LAT ----
  1315. run_lmer(Correct ~ kmed + (1|ID), data_vr_tests %>% filter(group=='CCD', PresentationType=='Lateralized', !ID %in% c('4003', '4004', '4007', '4008')))
  1316. run_lmer(Correct ~ kmed + (1|ID), data_vr_tests %>% filter(group=='NT', PresentationType=='Lateralized'))
  1318. run_lmer(conf ~ kmed + (1|ID), data_vr_tests %>% filter(group=='CCD', PresentationType=='Lateralized', !ID %in% c('4003', '4004', '4007', '4008')))
  1319. run_lmer(conf ~ kmed + (1|ID), data_vr_tests %>% filter(group=='NT', PresentationType=='Lateralized'))
  1321. #### MONO ----
  1323. run_lmer(Correct ~ kmed + (1|ID), data_vr_tests %>% filter(group=='CCD', PresentationType=='Monocular', !ID %in% c('4003', '4004', '4007', '4008')))
  1324. run_lmer(Correct ~ kmed + (1|ID), data_vr_tests %>% filter(group=='NT', PresentationType=='Monocular'))
  1326. run_lmer(conf ~ kmed + (1|ID), data_vr_tests %>% filter(group=='CCD', PresentationType=='Monocular', !ID %in% c('4003', '4004', '4007', '4008')))
  1327. run_lmer(conf ~ kmed + (1|ID), data_vr_tests %>% filter(group=='NT', PresentationType=='Monocular'))
  1329. #### Computational -----------------------------------------------------------
  1331. # Experiment 1
  1332. outputONL_CCD_exc <- readRDS('HMetaDFit/Mratio/ONL_CCD_exc.rdata')
  1333. Results_ONL_CCD_exc <- HMeta_post_clean(outputONL_CCD_exc, 'CCD (ONL)', F)
  1335. data_for_permute_exp1_CCD_exc <- checkBothRDK_w_exc %>% filter(group=='CCD (Online)') %>% dplyr::select(ID, Trial, conf, dotDirection, type, kmed, referenceSelection, group)
  1336. onl_CCDvsNTp_exc <- permute_metad_group_mrat(data_for_permute_exp1_CCD_exc, data_for_permute_exp1_NT, cores = 10, nreps = nreps);
  1338. #Experiment 3
  1339. data_for_permute_exp3_CCD_exc <- data_vr %>% filter(BinocularTrial == 1, Group=='CCD', !PID %in% exc_ids)
  1340. data_for_permute_exp3_CCD_L_exc<- data_vr %>% filter(LateralizedTrial == 1, Group=='CCD', !PID %in% exc_ids)
  1341. data_for_permute_exp3_CCD_M_exc<- data_vr %>% filter(MonocularTrial == 1, Group=='CCD', !PID %in% exc_ids)
  1343. VR_CCD_metaDat_exc <- HMetaGroupPrep_vr(data_vr %>% filter(BinocularTrial == 1, !PID %in% exc_ids), 'CCD')
  1344. outputVR_CCD_exc <- metad_group(nR_S1 = list(VR_CCD_metaDat_exc[[1]]), nR_S2 = list(VR_CCD_metaDat_exc[[2]]));
  1346. exc_output_clean <- HMeta_post_clean(outputVR_CCD_exc, 'CCD (VR)', F)
  1347. mu_log_vr_CCD_exc <- exc_output_clean$Fit %>% filter(name == 'mu_logMratio') %>% mutate(mean = exp(mean))
  1348. diffVR_exc = as.numeric(mu_log_exp3_NT[2] - mu_log_vr_CCD_exc[2])
  1350. vr_CCDvsNTp_exc <- permute_metad_group_vr(data_for_permute_exp3_CCD_exc, data_for_permute_exp3_NT,cores = 10, nreps = nreps);
  1351. perm_effect_stats(round(diffVR_exc, 2), vr_CCDvsNTp_exc[,4])

Analysis_Script_Commented_Github.r at commit 8cd58ed, under CC-BY-NC-SA-4.0 · at the source

Overview

Authors: J.M. Barnby1,2, R. Dean3,4, H. Burgess3,4, P. Dayan5,6, L.J. Richards3,4
ORCID iDs: J.M. Barnby, R. Dean, H. Burgess, L.J. Richards
  1. Institute of Psychiatry, Psychology and Neuroscience, King’s College London, UK
  2. Centre for AI and Machine Learning, Edith Cowan University, WA, Australia
  3. Department of Neuroscience, Washington University in St. Louis Medical School, St Louis, MO, USA
  4. The University of Queensland, Queensland Brain Institute, Brisbane, Australia
  5. Max Planck Institute for Biological Cybernetics, Tübingen, Germany
  6. University of Tübingen, Germany
Institutions: King's College London (United Kingdom); Edith Cowan University (Australia); Washington University in St. Louis (United States); The University of Queensland (Australia); Max Planck Institute for Biological Cybernetics (Germany); University of Tübingen (Germany)
Journal: bioRxiv (preprint)
Dates: published online 4 March 2026
Type: Preprint
License: CC BY
Identifiers: DOI 10.64898/2026.03.02.709173 · OpenAlex W7133771409
Open access: green, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), cognitive (subfield)
Methods: Statistics, Connectivity, Physiology & signal measures
Keywords: Agenesis of the corpus callosum, corpus callosum dysgenesis, perception, metacognition, social inference, abstract decision making
Topic: Fetal and Pediatric Neurological Disorders (Pediatrics, Perinatology and Child Health, Medicine), according to OpenAlex
Funding: Wellcome Trust (WT228268/Z/23/Z)
Citations: not cited yet (Europe PMC); 43 references in the paper

Abstract

The corpus callosum is the largest commissure in the mammalian brain and plays a major role in supporting cognitive processes required for adapting to complex environments. Individuals born with Corpus Callosum Dysgenesis (CCD), characterized by malformations of the corpus callosum, commonly exhibit deficits in social navigation, abstract problem-solving, decision-making, and self-awareness. Metacognition is a key cognitive process that supports these functions; however, it has yet to be tested comprehensively in individuals with CCD. Over three experiments, and three CCD cohorts, we tested the impact of this neurodevelopmental disorder on perceptual accuracy, confidence judgements, and metacognitive efficiency using two variants of a Random Dot Kinematogram task within lab, online, and VR conditions. We found that individuals with CCD typically displayed normal perceptual accuracy but failed to adjust their confidence judgements in line with task difficulty. Computational modelling revealed that this difference was explained by lower metacognitive efficiency driven by consistently lower metacognitive sensitivity. Together, these results provide evidence that the corpus callosum plays a crucial role in supporting metacognition.

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

Brain-Development-and-Disorders-Lab/Barnby_etal_2026_ccd_impairs_metacognition

License: CC-BY-NC-SA-4.0
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 8cd58ed054b724da8008732a7ec8978a3bea35d3, 30 July 2026
Languages: MATLAB (34), R (15)
Size: 115 files, 49 scripts
Software Heritage: not archived
Found in: “Data & Code”
Holds: README, license file, 1 notebook
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (14 files), tidyverse (10 files), broom (9 files), ggpubr (9 files), JAGS (9 files), reshape2 (9 files), lmerTest (2 files), brms (1 file), easystats (1 file), lme4 (1 file), Optimization Toolbox (1 file), Parallel Computing Toolbox (1 file), patchwork (1 file), Stan (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
51 files
At the source: github.com/Brain-Development-and-Disorders-Lab/Barnby_etal_2026_ccd_impairs_metacognition

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;
  • 49 scripts, each with its path and the digest of its content;
  • 4 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 & Code

Data and analysis code is freely available here: https://github.com/Brain-Development-and-Disorders-Lab/Barnby_etal_2026_ccd_impairs_metacognition

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

Versions

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

Version 1, 30 September 2026: the first record

Recorded: type, journal, dates, 5 authors, 6 keywords, 1 funder, 39 references.

Cite

This paper

Barnby, J., Dean, R., Burgess, H., Dayan, P., & Richards, L. (2026). Corpus Callosum Dysgenesis impairs metacognition: evidence from multi-modality and multi-cohort replications. bioRxiv (preprint). https://doi.org/10.64898/2026.03.02.709173

BibTeX

@article{barnby2026corpus,
author = {Barnby, J.M. and Dean, R. and Burgess, H. and Dayan, P. and Richards, L.J.},
title = {{Corpus Callosum Dysgenesis impairs metacognition: evidence from multi-modality and multi-cohort replications}},
journal = {bioRxiv (preprint)},
year = {2026},
month = mar,
publisher = {bioRxiv},
issn = {2692-8205},
doi = {10.64898/2026.03.02.709173},
url = {https://doi.org/10.64898/2026.03.02.709173}
}

RIS

TY - JOUR
AU - Barnby, J.M.
AU - Dean, R.
AU - Burgess, H.
AU - Dayan, P.
AU - Richards, L.J.
TI - Corpus Callosum Dysgenesis impairs metacognition: evidence from multi-modality and multi-cohort replications
T2 - bioRxiv (preprint)
J2 - bioRxiv
PY - 2026
DA - 2026/03/04
SN - 2692-8205
PB - bioRxiv
DO - 10.64898/2026.03.02.709173
UR - https://doi.org/10.64898/2026.03.02.709173
ER -

CSL-JSON

{
"id": "10.64898/2026.03.02.709173",
"type": "article",
"title": "Corpus Callosum Dysgenesis impairs metacognition: evidence from multi-modality and multi-cohort replications",
"container-title": "bioRxiv (preprint)",
"author": [
{
"family": "Barnby",
"given": "J.M."
},
{
"family": "Dean",
"given": "R."
},
{
"family": "Burgess",
"given": "H."
},
{
"family": "Dayan",
"given": "P."
},
{
"family": "Richards",
"given": "L.J."
}
],
"container-title-short": "bioRxiv",
"DOI": "10.64898/2026.03.02.709173",
"ISSN": "2692-8205",
"publisher": "bioRxiv",
"URL": "https://doi.org/10.64898/2026.03.02.709173",
"issued": {
"date-parts": [
[
2026,
3,
4
]
]
}
}

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.1162/imag.a.1332 [code]
Decision processes underlying effort avoidance and their relationship with metacognition.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: JAGS, Optimization Toolbox, broom, 5 other tools, cognitive, 3 references
[2] doi:10.1016/j.celrep.2026.117505 [code]
Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice.
Journal: Cell reports
In common: JAGS, Stan, brms, 7 other tools
[3] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Stan, brms, Optimization Toolbox, 7 other tools
[4] doi:10.1162/imag.a.1258 [code]
Non-specific increase in alpha power during a neurofeedback session targeting its downregulation.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Stan, brms, easystats, 7 other tools
[5] doi:10.1038/s44271-026-00431-w [code]
Alpha power increases spontaneously during a neurofeedback session.
Journal: Communications psychology
In common: Stan, brms, easystats, 6 other tools, cognitive
[6] doi:10.1002/hbm.70605 [code]
BrainEnrich: Revealing Biological Insights for Imaging-Derived Phenotypes Through Transcriptomic Enrichment.
Journal: Human brain mapping
In common: JAGS, easystats, broom, 6 other tools
[7] doi:10.1111/psyp.70285 [code]
Spontaneous Modulation of Alpha Power During a Neurofeedback Session Without Instructions.
Journal: Psychophysiology
In common: Stan, brms, easystats, 6 other tools
[8] doi:10.1016/j.isci.2026.116747 [code]
Age and loneliness relate to reduced trust learning and alterations in amygdala function.
Journal: iScience
In common: Stan, brms, easystats, 5 other tools, cognitive
[9] doi:10.1093/braincomms/fcag217 [code]
Investigation of stress hormones across multiday seizure cycles.
Journal: Brain communications
In common: brms, easystats, broom, 6 other tools
[10] doi:10.1186/s13229-026-00730-3 [code]
Predicting emotional valence in autism: a preregistered study in the Bayesian Brain framework.
Journal: Molecular autism
In common: Stan, brms, easystats, 5 other tools, cognitive

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.