Emergent blink rate in early childhood is associated with neural origins of executive function.

1. # ==============================================================================
2. # MainAnalysis.R
3. # Emergent blink rate in early childhood is associated with neural origins
4. # of executive function
5. # Kuwamizu, Yamamoto, Otani & Moriguchi
6. #
7. # Copyright (c) 2026 the authors. Released under the MIT License; see code/LICENSE.
8. # The accompanying data are released under CC BY 4.0; see LICENSE.md.
9. #
10. # Integrated analysis pipeline. Reproduces every statistic reported in the
11. # manuscript and Supplementary Information, and exports figure source data.
12. #
13. # Section 1 Behavioural correlations Fig 2a-c; Supp Fig 1a-d
14. # Section 2 fNIRS linear mixed models Fig 3a-f; Supp Table 3
15. # Section 3 Neural correlations Fig 4a-b; Supp Table 4; Supp Fig 5
16. # Section 4 Robustness and sensitivity Supp Notes 1 and 4; Methods
17. # Section 5 Summary CSVs and session information
18. #
19. # Input : data/SourceData.csv
20. # Output: outputs_summary/results_summary_generated.csv
21. # source_data_figures/*.csv
22. # R >= 4.3
23. # ==============================================================================
26. ## ====== 0. Setup, data, helpers ==============================================
27. if (!require("pacman")) install.packages("pacman")
28. pacman::p_load(tidyverse, lme4, lmerTest, emmeans, ppcor, rstatix,
29. effectsize, performance, mediation)
31. options(scipen = 999)
32. set.seed(365)
34. OUTDIR <- "outputs_summary"
35. FIGDIR <- "source_data_figures"
36. dir.create(OUTDIR, showWarnings = FALSE, recursive = TRUE)
37. dir.create(FIGDIR, showWarnings = FALSE, recursive = TRUE)
39. # ---- single data path (override here if needed) ----
40. candidate_paths <- c(
41. file.path("data", "SourceData.csv"),
42. "SourceData.csv",
43. file.path("..", "data", "SourceData.csv")
44. )
45. data_path <- candidate_paths[file.exists(candidate_paths)][1]
46. if (is.na(data_path)) stop("SourceData.csv not found; set data_path manually.")
47. cat("Loading data from:", data_path, "\n")
49. raw_data <- read.csv(data_path, na.strings = c("NA", "NaN", ""),
50. fileEncoding = "UTF-8", stringsAsFactors = FALSE) %>%
51. mutate(ID = as.factor(ID))
53. fnirs_cols <- grep("^(r|l)_(DL|RL)_(pre|post|mix|sw)$", names(raw_data), value = TRUE)
54. count_cols <- c("sw_corr","sw_rep","Pre_corr","Pre_rep",
55. "Post_corr","Post_rep","Mix_corr","Mix_rep")
56. target_cols <- c("age","sEBR","sw_acc","pre_acc","post_acc","mix_acc",
57. count_cols, fnirs_cols)
58. raw_data <- raw_data %>%
59. mutate(across(any_of(target_cols), ~ as.numeric(as.character(.))),
60. sex = if ("sex" %in% names(.)) factor(sex) else NA,
61. site = if ("site" %in% names(.)) factor(site) else NA)
62. cat("N participants:", nrow(raw_data), "\n")
64. # ---- helper functions ----
65. # Spearman with N and 95% CI.
66. # CI from Fisher's z transform with the Bonett & Wright (2000) standard error
67. # for Spearman's rho: SE_z = sqrt((1 + rho^2/2)/(n - 3)).
68. sp <- function(x, y, conf = 0.95) {
69. d <- na.omit(data.frame(x, y))
70. n <- nrow(d)
71. ct <- cor.test(d$x, d$y, method = "spearman", exact = FALSE)
72. rho <- unname(ct$estimate)
73. if (n > 3 && abs(rho) < 1) {
74. z <- atanh(rho)
75. se <- sqrt((1 + rho^2 / 2) / (n - 3))
76. zc <- qnorm(1 - (1 - conf) / 2)
77. lo <- tanh(z - zc * se); hi <- tanh(z + zc * se)
78. } else { lo <- NA_real_; hi <- NA_real_ }
79. c(rho = rho, p = ct$p.value, n = n, ci_low = lo, ci_high = hi)
80. }
81. # Partial Spearman via ppcor (df = n - 2 - k); covariates are ranked.
82. # CI from Fisher's z transform with SE_z = 1/sqrt(n - k - 3) (k = # covariates).
83. pspear <- function(x, y, Z, conf = 0.95) {
84. Z <- as.data.frame(Z)
85. d <- na.omit(data.frame(x = x, y = y, Z))
86. zc <- setdiff(names(d), c("x", "y"))
87. k <- length(zc)
88. n <- nrow(d)
89. zr <- as.data.frame(lapply(d[zc], function(v) rank(xtfrm(v))))
90. pc <- ppcor::pcor.test(rank(d$x), rank(d$y), zr, method = "pearson")
91. rho <- unname(pc$estimate)
92. if ((n - k - 3) > 0 && abs(rho) < 1) {
93. z <- atanh(rho)
94. se <- 1 / sqrt(n - k - 3)
95. zcrit <- qnorm(1 - (1 - conf) / 2)
96. lo <- tanh(z - zcrit * se); hi <- tanh(z + zcrit * se)
97. } else { lo <- NA_real_; hi <- NA_real_ }
98. c(rho = rho, p = pc$p.value, n = n, ci_low = lo, ci_high = hi)
99. }
100. winsorize <- function(v, p = 0.05) {
101. lo <- quantile(v, p, na.rm = TRUE); hi <- quantile(v, 1-p, na.rm = TRUE)
102. pmin(pmax(v, lo), hi)
103. }
105. # ---- running summary collector (-> results_summary_generated.csv) ----
106. # Every value written to the CSV is also echoed to the console (prefixed "[CSV]"),
107. # so the console log and the CSV contain the same numbers.
108. SUMMARY <- list()
109. # `stat` : NAME of the test statistic or estimate type (e.g. "t", "V", "F")
110. # `stat_value` : VALUE of that test statistic, when it differs from `estimate`
111. # (e.g. estimate = Cohen's d while stat_value = the t statistic).
112. # Left empty when `estimate` already IS the statistic.
113. # `stat_value` is the LAST argument so that every pre-existing positional call
114. # in this script keeps its original meaning.
115. add_row <- function(loc, analysis, estimate, stat = "", df = "",
116. p = "", n = "", model = "", note = "",
117. ci_low = "", ci_high = "", stat_value = "") {
118. row <- data.frame(
119. ms_location = loc, analysis = analysis, estimate = as.character(estimate),
120. stat = stat, stat_value = as.character(stat_value),
121. df = as.character(df), p = as.character(p),
122. n = as.character(n),
123. ci_low = as.character(ci_low), ci_high = as.character(ci_high),
124. model = model, note = note, stringsAsFactors = FALSE)
125. SUMMARY[[length(SUMMARY) + 1]] <<- row
126. # ---- console echo of the same numbers ----
127. e <- as.character(estimate); st <- as.character(stat); d <- as.character(df)
128. sv <- as.character(stat_value)
129. pp <- as.character(p); nn <- as.character(n)
130. cl <- as.character(ci_low); ch <- as.character(ci_high)
131. md <- as.character(model); nt <- as.character(note)
132. part <- function(lab, v) if (length(v) && !is.na(v) && nzchar(v)) paste0(" ", lab, v) else ""
133. ci_str <- if (length(cl) && length(ch) && !is.na(cl) && !is.na(ch) && nzchar(cl) && nzchar(ch))
134. sprintf(" 95%%CI=[%s, %s]", cl, ch) else ""
135. # when a statistic value is present, echo it as e.g. " t=3.214"
136. st_str <- if (length(sv) && !is.na(sv) && nzchar(sv))
137. paste0(" ", st, "=", sv) else part("", st)
138. cat(sprintf("[CSV] %-18s | %-46s | est=%s%s%s%s%s%s%s%s\n",
139. loc, analysis, e,
140. st_str, part("df=", d), part("p=", pp), part("n=", nn), ci_str,
141. part("model=", md), part("note=", nt)))
142. }
143. fmtp <- function(p) ifelse(p < 0.001, "<0.001", signif(p, 3))
144. # 95% CI as a printable string, e.g. "[0.12, 0.34]"
145. fmtci <- function(lo, hi, d = 3) ifelse(is.na(lo) | is.na(hi), "",
146. sprintf(paste0("[%.", d, "f, %.", d, "f]"), lo, hi))
149. ## =============================================================================
150. ## SECTION 1 -- Behavioural correlations [Fig 2a-c; Supp Fig 1a-d]
151. ## =============================================================================
152. cat("\n\n========== SECTION 1: Behaviour (Fig 2; Supp Fig 1) ==========\n")
154. ## 1-1 Friedman + post-hoc [Supp Fig 1a]
155. beh_long <- raw_data %>%
156. dplyr::select(ID, pre_acc, post_acc, mix_acc) %>%
157. pivot_longer(c(pre_acc, post_acc, mix_acc), names_to = "Condition", values_to = "Accuracy") %>%
158. mutate(Condition = factor(Condition, levels = c("pre_acc","post_acc","mix_acc"))) %>%
159. drop_na()
160. fr <- friedman.test(Accuracy ~ Condition | ID, data = beh_long); print(fr)
161. # The omnibus Friedman test is reported in full as chi2(df) and an exact p.
162. # (No separate effect size is added here: the editor's effect-size + CI
163. # requirement targets t-tests and ANOVAs; the rank-based post-hoc comparisons
164. # below carry the effect sizes and CIs.)
165. n_blk <- nrow(raw_data %>% dplyr::select(pre_acc, post_acc, mix_acc) %>% tidyr::drop_na())
166. add_row("SuppFig1a","Pre/Post/Mix accuracy (Friedman)", round(unname(fr$statistic),2),
167. "chi2", fr$parameter, fmtp(fr$p.value), n_blk)
169. # Effect size for the Friedman test: Kendall's W with 95% CI.
170. kw <- as.data.frame(suppressWarnings(
171. effectsize::kendalls_w(Accuracy ~ Condition | ID, data = beh_long,
172. ci = 0.95, alternative = "two.sided")))
173. cat(sprintf("Friedman effect size: Kendall's W = %.3f, 95%% CI [%.3f, %.3f]\n",
174. kw$Kendalls_W, kw$CI_low, kw$CI_high))
175. add_row("SuppFig1a","Friedman effect size (Kendall's W)", round(kw$Kendalls_W,3),
176. "Kendall W", "", "", n_blk, "", "",
177. round(kw$CI_low,3), round(kw$CI_high,3))
179. # Post-hoc paired Wilcoxon signed-rank tests, reported in full and separately.
180. # Statistic V; p from the normal approximation with continuity correction
181. # (exact = FALSE, appropriate here because accuracy values are tied), then
182. # Bonferroni-adjusted across the 2 comparisons; matched-pairs rank-biserial
183. # correlation r (+95% CI) as the effect size.
184. # NOTE: the signed-rank statistic V has no degrees of freedom, so V is
185. # written to `stat_value` and the `df` column is left empty.
186. ph_pairs <- list(c("pre_acc","post_acc","Pre vs Post"),
187. c("pre_acc","mix_acc","Pre vs Mix"))
188. ph_raw_p <- sapply(ph_pairs, function(pr)
189. wilcox.test(raw_data[[pr[1]]], raw_data[[pr[2]]], paired = TRUE, exact = FALSE)$p.value)
190. ph_adj_p <- p.adjust(ph_raw_p, method = "bonferroni")
191. for (i in seq_along(ph_pairs)) {
192. pr <- ph_pairs[[i]]
193. d2 <- na.omit(data.frame(a = raw_data[[pr[1]]], b = raw_data[[pr[2]]]))
194. wt <- wilcox.test(d2$a, d2$b, paired = TRUE, exact = FALSE)
195. rb <- tryCatch(as.data.frame(effectsize::rank_biserial(d2$a, d2$b, paired = TRUE, ci = 0.95)),
196. error = function(e) data.frame(r_rank_biserial = NA, CI_low = NA, CI_high = NA))
197. cat(sprintf("Post-hoc %-11s V = %.0f, p_adj = %s, r_rb = %.3f, 95%% CI %s, n = %d\n",
198. pr[3], unname(wt$statistic), fmtp(ph_adj_p[i]),
199. rb$r_rank_biserial, fmtci(rb$CI_low, rb$CI_high), nrow(d2)))
200. add_row("SuppFig1a", paste0("Post-hoc Wilcoxon ", pr[3]),
201. estimate = round(rb$r_rank_biserial, 3),
202. stat = "V",
203. df = "", # V has no df
204. p = fmtp(ph_adj_p[i]),
205. n = nrow(d2),
206. note = paste0("two-sided, normal approximation; Bonferroni-adjusted p; ",
207. "estimate = matched-pairs rank-biserial r (95% CI)"),
208. ci_low = round(rb$CI_low, 3),
209. ci_high = round(rb$CI_high, 3),
210. stat_value = round(unname(wt$statistic), 1))
211. }
213. ## 1-2 Primary correlations [Fig 2a-c]
214. for (nm in list(c("age","sw_acc","Fig2a","age ~ switching accuracy"),
215. c("age","sEBR","Fig2b","age ~ sEBR"),
216. c("sEBR","sw_acc","Fig2c","sEBR ~ switching accuracy"))) {
217. r <- sp(raw_data[[nm[1]]], raw_data[[nm[2]]])
218. cat(sprintf("%-28s rho=%.2f p=%s 95%% CI %s n=%d\n",
219. nm[4], r["rho"], fmtp(r["p"]), fmtci(r["ci_low"], r["ci_high"]), r["n"]))
220. add_row(nm[3], nm[4], round(r["rho"],3), "Spearman rho", "", fmtp(r["p"]), r["n"],
221. "", "", round(r["ci_low"],3), round(r["ci_high"],3))
222. }
223. r <- pspear(raw_data$sEBR, raw_data$sw_acc, raw_data["age"])
224. cat(sprintf("sEBR ~ switching acc | age rho=%.2f p=%s 95%% CI %s n=%d\n",
225. r["rho"], fmtp(r["p"]), fmtci(r["ci_low"], r["ci_high"]), r["n"]))
226. add_row("Fig2c","sEBR ~ switching accuracy | age", round(r["rho"],3),
227. "partial Spearman", "", fmtp(r["p"]), r["n"],
228. "", "", round(r["ci_low"],3), round(r["ci_high"],3))
230. ## 1-3 Phase-wise sEBR x accuracy + age x phase-accuracy [Supp Fig 1b-d]
231. for (v in c("pre_acc","post_acc","mix_acc")) {
232. rs <- sp(raw_data$sEBR, raw_data[[v]]); rp <- pspear(raw_data$sEBR, raw_data[[v]], raw_data["age"])
233. add_row("SuppFig1", paste0("sEBR ~ ", v), round(rs["rho"],3), "Spearman rho","", fmtp(rs["p"]), rs["n"],
234. "", "", round(rs["ci_low"],3), round(rs["ci_high"],3))
235. add_row("SuppFig1", paste0("sEBR ~ ", v, " | age"), round(rp["rho"],3), "partial Spearman","", fmtp(rp["p"]), rp["n"],
236. "", "", round(rp["ci_low"],3), round(rp["ci_high"],3))
237. }
238. for (v in c("pre_acc","post_acc","mix_acc")) {
239. ra <- sp(raw_data$age, raw_data[[v]])
240. add_row("SuppFig1", paste0("age ~ ", v), round(ra["rho"],3), "Spearman rho","", fmtp(ra["p"]), ra["n"],
241. "", "", round(ra["ci_low"],3), round(ra["ci_high"],3))
242. }
245. ## =============================================================================
246. ## SECTION 2 -- fNIRS LMM [Fig 3a-f; Supp Table 3]
247. ## =============================================================================
248. cat("\n\n========== SECTION 2: fNIRS LMM (Fig 3; Supp Table 3) ==========\n")
250. df_long <- raw_data %>%
251. pivot_longer(matches("^(r|l)_(DL|RL)_(pre|post|mix)$"),
252. names_to = c("Hemisphere","Region","Phase"),
253. names_pattern = "^([rl])_([A-Z]{2})_([a-z]+)", values_to = "OxyHb") %>%
254. mutate(Subject = ID,
255. Hemisphere = factor(Hemisphere, levels = c("l","r"), labels = c("Left","Right")),
256. Region = factor(Region, levels = c("RL","DL"), labels = c("RLPFC","DLPFC")),
257. Phase = factor(Phase, levels = c("pre","mix","post")),
258. PhaseRep = dplyr::case_when(Phase=="pre"~Pre_rep, Phase=="mix"~Mix_rep, Phase=="post"~Post_rep)) %>%
259. drop_na(OxyHb)
260. cat("df_long rows:", nrow(df_long), "\n")
262. ## 2-1 Base model [Fig 3b-d]
263. model_base <- lmer(OxyHb ~ Phase*Hemisphere*Region + (1|Subject), data = df_long)
264. print(anova(model_base))
265. ab <- anova(model_base)
266. for (eff in c("Phase","Hemisphere","Region")) {
267. add_row("Fig3", paste0(eff," main effect"), round(ab[eff,"F value"],2), "F",
268. paste0(ab[eff,"NumDF"], ",", sprintf("%.1f", ab[eff,"DenDF"])), fmtp(ab[eff,"Pr(>F)"]), "", "base")
269. }
270. # Interaction terms reported in full (F, df, exact p) even though non-significant.
271. for (eff in c("Phase:Hemisphere","Phase:Region","Hemisphere:Region","Phase:Hemisphere:Region")) {
272. add_row("Fig3", paste0(eff," interaction"), round(ab[eff,"F value"],2), "F",
273. paste0(ab[eff,"NumDF"], ",", sprintf("%.1f", ab[eff,"DenDF"])), fmtp(ab[eff,"Pr(>F)"]), "", "base")
274. }
275. ## Base-model post-hoc EMM contrasts [Fig 3c,d]
276. ## Reported in full: EMM difference with 95% CI, t(df), p.
277. ## >>> SIGN: contrasts are (first level - second level); all five come out
278. ## >>> NEGATIVE here. The manuscript quotes magnitudes. See the header block.
279. ## `infer = c(TRUE, TRUE)` is added so that the confidence interval is returned
280. ## alongside the test; emmeans applies the SAME multiplicity adjustment to the
281. ## interval as to the p value, so both carry the adjustment named below.
282. ## Phase has three levels -> Tukey adjustment is active.
283. ## Hemisphere and Region have two levels -> a single contrast, no adjustment.
284. ## The corresponding standardised effect sizes (Cohen's d) are written
285. ## separately by report_effsizes() above.
286. emm_contrast_rows <- function(model, fac, adjust, loc, tag) {
287. ct <- as.data.frame(summary(
288. emmeans(model, as.formula(paste0("pairwise ~ ", fac)), adjust = adjust)$contrasts,
289. infer = c(TRUE, TRUE)))
290. cat(sprintf("\n-- %s post-hoc contrasts (adjust = %s) --\n", fac, adjust))
291. print(ct)
292. for (j in seq_len(nrow(ct))) {
293. add_row(loc, paste0(fac, " post-hoc (", ct$contrast[j], ")"),
294. estimate = round(ct$estimate[j], 4),
295. stat = "t",
296. df = sprintf("%.1f", ct$df[j]),
297. p = fmtp(ct$p.value[j]),
298. n = "",
299. model = tag,
300. note = paste0("estimate = UNSTANDARDISED EMM difference in oxy-Hb units ",
301. "(95% CI on that difference, NOT on Cohen's d); ",
302. "sign follows the contrast label (first - second level); ",
303. "p and CI adjustment: ", adjust,
304. " (", nrow(ct), " contrast(s)); ",
305. "the standardised effect size for this contrast is reported ",
306. "separately as a 'Cohen d' row"),
307. ci_low = round(ct$lower.CL[j], 4),
308. ci_high = round(ct$upper.CL[j], 4),
309. stat_value = round(ct$t.ratio[j], 3))
310. }
311. ct
312. }
313. cph <- emm_contrast_rows(model_base, "Phase", "tukey", "Fig3c", "base")
314. chh <- emm_contrast_rows(model_base, "Hemisphere", "none", "Fig3d", "base")
315. crr <- emm_contrast_rows(model_base, "Region", "none", "Fig3d", "base")
317. ## 2-2 Effect sizes for BASE and AGE models
318. report_effsizes <- function(model, label, tag) {
319. cat(sprintf("\n--- Effect sizes [%s] ---\n", label))
320. cat("partial eta^2:\n"); es <- eta_squared(anova(model, type = 3), partial = TRUE, ci = 0.95, alternative = "two.sided"); print(es)
321. cat("R^2 (Nakagawa):\n"); r2 <- r2_nakagawa(model); print(r2)
322. sig <- sigma(model); edf <- df.residual(model)
323. for (fac in c("Phase","Hemisphere","Region")) {
324. cat(sprintf(" Cohen d [%s]:\n", fac))
325. esz <- eff_size(emmeans(model, as.formula(paste0("~ ", fac))), sigma = sig, edf = edf)
326. print(esz)
327. ed <- as.data.frame(esz)
328. for (j in seq_len(nrow(ed))) {
329. add_row(paste0("Fig3 effsize [",tag,"]"),
330. paste0("Cohen d ", fac, ": ", ed$contrast[j]),
331. round(ed$effect.size[j],3), "Cohen d", round(ed$df[j],1), "", "",
332. tag, paste0("post-hoc EMM contrast (sign = first - second level); ",
333. "eff_size() calls contrast(adjust = \"none\"), so this 95% CI is ",
334. "UNADJUSTED for multiplicity even where the corresponding p value ",
335. "is Tukey-adjusted; sigma = residual SD of the LMM, edf = df.residual"),
336. round(ed$lower.CL[j],3), round(ed$upper.CL[j],3))
337. }
338. }
339. es <- as.data.frame(es)
340. for (i in seq_len(nrow(es))) {
341. add_row(paste0("Fig3 effsize [",tag,"]"), paste0("partial eta2 ", es$Parameter[i]),
342. signif(es$Eta2_partial[i],3), "eta2_p","","","", tag, "",
343. round(es$CI_low[i],3), round(es$CI_high[i],3))
344. }
345. add_row(paste0("Fig3 effsize [",tag,"]"),"Marginal R2", round(r2$R2_marginal,3), "R2","","","", tag)
346. add_row(paste0("Fig3 effsize [",tag,"]"),"Conditional R2",round(r2$R2_conditional,3), "R2","","","", tag)
347. }
349. ## 2-3 Age model [Fig 3e,f]
350. model_age <- lmer(OxyHb ~ (Phase*Hemisphere*Region)*age + (1|Subject), data = df_long)
351. cat("\nBase vs Age LRT:\n"); print(anova(model_base, model_age))
352. print(anova(model_age))
353. aa <- anova(model_age)
354. ## Every term of the age-model ANOVA is written to the summary CSV, not only the
355. ## three significant age interactions: the manuscript reports F, df and P for the
356. ## non-significant higher-order interactions as well.
357. for (eff in rownames(aa)) {
358. add_row("Fig3e/f", eff, round(aa[eff,"F value"],2), "F",
359. paste0(aa[eff,"NumDF"], ",", sprintf("%.1f", aa[eff,"DenDF"])),
360. fmtp(aa[eff,"Pr(>F)"]), "", "age",
361. "Type III ANOVA, Satterthwaite df")
362. }
363. ## Post-hoc slope contrasts decomposing the age interactions (Fig 3e,f).
364. ## Reported in full: estimate (difference in age slopes) with 95% CI, t(df), p.
365. ## >>> SIGN: (first level - second level), as for the base-model contrasts above.
366. ## >>> See the sign-convention block in the file header before transcribing.
367. for (fac in c("Hemisphere","Region","Phase")) {
368. ct <- as.data.frame(summary(
369. emtrends(model_age, as.formula(paste0("pairwise ~ ", fac)), var = "age")$contrasts,
370. infer = c(TRUE, TRUE)))
371. cat(sprintf("\n-- age-slope contrasts by %s --\n", fac)); print(ct)
372. for (j in seq_len(nrow(ct))) {
373. # emmeans applies a Tukey adjustment to BOTH the p values and the CIs of
374. # pairwise contrasts by default. Hemisphere and Region have two levels
375. # (a single contrast, so the adjustment is inert); Phase has three levels,
376. # where the adjustment is active. This is recorded in the note.
377. add_row("Fig3e/f", paste0(fac, ":age contrast (", ct$contrast[j], ")"),
378. estimate = round(ct$estimate[j], 4),
379. stat = "t",
380. df = sprintf("%.1f", ct$df[j]),
381. p = fmtp(ct$p.value[j]),
382. n = "",
383. model = "age",
384. note = paste0("estimate = difference in age slopes (95% CI); ",
385. "Tukey-adjusted p and CI across the ", nrow(ct),
386. " contrast(s) of ", fac),
387. ci_low = round(ct$lower.CL[j], 4),
388. ci_high = round(ct$upper.CL[j], 4),
389. stat_value = round(ct$t.ratio[j], 3))
390. }
391. }
393. report_effsizes(model_base, "BASE model: Phase x Hemisphere x Region", "BASE")
394. report_effsizes(model_age, "AGE model: (Phase x Hemisphere x Region) x age", "AGE")
396. ## 2-3b Residual diagnostics for the mixed-effects models [Methods]
397. ## Normality of residuals and of the random intercepts, dependence of the
398. ## residual variance on the fitted values, and equality of residual variance
399. ## across the design cells.
400. lmm_diagnostics <- function(model, tag) {
401. mf <- model@frame
402. r <- residuals(model)
403. f <- fitted(model)
404. sw <- shapiro.test(r)
405. sk <- mean((r - mean(r))^3) / sd(r)^3
406. ku <- mean((r - mean(r))^4) / sd(r)^4 - 3
407. u <- r^2 / mean(r^2)
408. LM <- 0.5 * sum((fitted(lm(u ~ f)) - 1)^2)
409. bp <- pchisq(LM, 1, lower.tail = FALSE)
410. cell <- interaction(mf$Phase, mf$Hemisphere, mf$Region, drop = TRUE)
411. z <- abs(r - ave(r, cell, FUN = median))
412. bf <- anova(lm(z ~ cell))
413. sds <- tapply(r, cell, sd)
414. re <- ranef(model)$Subject[, 1]
415. swr <- shapiro.test(re)
417. cat(sprintf("\n--- Residual diagnostics [%s] ---\n", tag))
418. cat(sprintf("Residual normality: W = %.4f, p = %s (skewness %.3f, excess kurtosis %.3f)\n",
419. sw$statistic, fmtp(sw$p.value), sk, ku))
420. cat(sprintf("Residual variance vs fitted values: LM(1) = %.3f, p = %s\n", LM, fmtp(bp)))
421. cat(sprintf("Residual variance across design cells: F(%d, %d) = %.3f, p = %s\n",
422. bf$Df[1], bf$Df[2], bf$`F value`[1], fmtp(bf$`Pr(>F)`[1])))
423. cat(sprintf("Residual SD by cell: %.4f to %.4f (max/min = %.2f)\n",
424. min(sds), max(sds), max(sds) / min(sds)))
425. cat(sprintf("Random-intercept normality: W = %.4f, p = %s\n",
426. swr$statistic, fmtp(swr$p.value)))
428. add_row("Methods", paste0("Residual normality [", tag, "]"),
429. round(unname(sw$statistic), 4), "Shapiro-Wilk W", "", fmtp(sw$p.value), length(r), tag)
430. add_row("Methods", paste0("Residual skewness / excess kurtosis [", tag, "]"),
431. sprintf("%.3f / %.3f", sk, ku), "desc", "", "", length(r), tag)
432. add_row("Methods", paste0("Residual variance vs fitted [", tag, "]"),
433. round(LM, 3), "Breusch-Pagan LM", 1, fmtp(bp), length(r), tag)
434. add_row("Methods", paste0("Residual variance across cells [", tag, "]"),
435. round(bf$`F value`[1], 3), "Brown-Forsythe F",
436. paste0(bf$Df[1], ",", bf$Df[2]), fmtp(bf$`Pr(>F)`[1]), length(r), tag)
437. add_row("Methods", paste0("Residual SD ratio across cells [", tag, "]"),
438. round(max(sds) / min(sds), 2), "max/min", "", "", length(sds), tag)
439. add_row("Methods", paste0("Random-intercept normality [", tag, "]"),
440. round(unname(swr$statistic), 4), "Shapiro-Wilk W", "", fmtp(swr$p.value), length(re), tag)
441. }
442. lmm_diagnostics(model_base, "BASE")
443. lmm_diagnostics(model_age, "AGE")
445. ## 2-3c Effect size for the base-vs-age model comparison [Fig 3e,f]
446. r2_base <- r2_nakagawa(model_base)
447. r2_age <- r2_nakagawa(model_age)
448. d_marg <- r2_age$R2_marginal - r2_base$R2_marginal
449. d_cond <- r2_age$R2_conditional - r2_base$R2_conditional
450. lrt_ba <- anova(model_base, model_age)
451. aic_diff <- lrt_ba$AIC[2] - lrt_ba$AIC[1]
452. cat(sprintf("\nBase vs Age model comparison: change in marginal R2 = %.3f, ", d_marg))
453. cat(sprintf("change in conditional R2 = %.3f, AIC difference = %.1f\n", d_cond, aic_diff))
454. add_row("Fig3e/f","Base vs Age: likelihood ratio test", round(lrt_ba$Chisq[2],2), "chi2",
455. lrt_ba$Df[2], fmtp(lrt_ba$`Pr(>Chisq)`[2]), "", "age",
456. "LRT of the base model against the model adding age and all its interactions")
457. add_row("Fig3e/f","Base vs Age: change in marginal R2", round(d_marg,3), "dR2","","","","age")
458. add_row("Fig3e/f","Base vs Age: change in conditional R2", round(d_cond,3), "dR2","","","","age")
459. add_row("Fig3e/f","Base vs Age: AIC difference", round(aic_diff,1), "dAIC","","","","age",
460. "AIC(age) - AIC(base); equals 2*Df - Chisq of the LRT above")
462. ## 2-4 One-sample t-tests vs baseline [Supp Table 3]
463. ## Full reporting per cell: t(df), exact p (+Bonferroni), Cohen's d with 95% CI.
464. st3 <- df_long %>% group_by(Phase, Hemisphere, Region) %>%
465. group_modify(~{
466. x <- .x$OxyHb
467. tt <- t.test(x, mu = 0)
468. dd <- tryCatch(as.data.frame(effectsize::cohens_d(x, mu = 0, ci = 0.95)),
469. error = function(e) data.frame(Cohens_d = mean(x)/sd(x),
470. CI_low = NA, CI_high = NA))
471. tibble(N = length(x), Mean_Oxy = mean(x, na.rm = TRUE),
472. Mean_CI_low = unname(tt$conf.int[1]), Mean_CI_high = unname(tt$conf.int[2]),
473. t_value = unname(tt$statistic), df = unname(tt$parameter),
474. p_raw = tt$p.value,
475. Cohen_d = dd$Cohens_d, d_CI_low = dd$CI_low, d_CI_high = dd$CI_high)
476. }) %>%
477. ungroup() %>%
478. mutate(p_bonf = p.adjust(p_raw, method = "bonferroni")) %>%
479. arrange(Phase, Hemisphere, Region)
480. print(as.data.frame(st3))
481. write.csv(st3, file.path(FIGDIR, "source_SuppTable3_onesample.csv"), row.names = FALSE)
482. for (i in seq_len(nrow(st3))) {
483. add_row("SuppTable3",
484. paste0("one-sample t vs 0: ", st3$Phase[i], " ", st3$Hemisphere[i], " ", st3$Region[i]),
485. estimate = round(st3$Cohen_d[i], 3),
486. stat = "t",
487. df = st3$df[i],
488. p = fmtp(st3$p_raw[i]),
489. n = st3$N[i],
490. note = paste0("two-sided; p shown is raw, p_bonf=", fmtp(st3$p_bonf[i]),
491. "; estimate = Cohen d (95% CI)"),
492. ci_low = round(st3$d_CI_low[i], 3),
493. ci_high = round(st3$d_CI_high[i], 3),
494. stat_value = round(st3$t_value[i], 3))
495. ## The Oxy-Hb column of Supplementary Table 3 is the cell mean, which is a
496. ## different quantity from Cohen's d above and is reported in mM*mm. Written
497. ## as its own row so that the table can be reproduced from this file alone.
498. add_row("SuppTable3",
499. paste0("mean oxy-Hb: ", st3$Phase[i], " ", st3$Hemisphere[i], " ", st3$Region[i]),
500. estimate = round(st3$Mean_Oxy[i], 3),
501. stat = "mean",
502. n = st3$N[i],
503. note = "cell mean oxy-Hb in mM*mm with 95% CI; Oxy-Hb column of Supp Table 3",
504. ci_low = round(st3$Mean_CI_low[i], 3),
505. ci_high = round(st3$Mean_CI_high[i], 3))
506. }
508. ## 2-5 Fig 3c-f source-data export (EMMs for Prism) [Fig 3c,d,e,f]
509. write.csv(as.data.frame(summary(emmeans(model_base, ~ Phase))),
510. file.path(FIGDIR, "Fig3c_phase_emm.csv"), row.names = FALSE)
511. write.csv(as.data.frame(summary(emmeans(model_base, ~ Hemisphere*Region))),
512. file.path(FIGDIR, "Fig3d_hemi_region_emm.csv"), row.names = FALSE)
513. age_seq <- seq(floor(min(df_long$age)), ceiling(max(df_long$age)), by = 1)
514. write.csv(as.data.frame(summary(emmeans(model_age, ~ Hemisphere | age, at = list(age = age_seq)), infer = TRUE)),
515. file.path(FIGDIR, "Fig3e_hemisphere_age.csv"), row.names = FALSE)
516. write.csv(as.data.frame(summary(emmeans(model_age, ~ Region | age, at = list(age = age_seq)), infer = TRUE)),
517. file.path(FIGDIR, "Fig3f_region_age.csv"), row.names = FALSE)
520. ## =============================================================================
521. ## SECTION 3 -- Neural correlations [Fig 4a-b; Supp Table 4; Supp Fig 5]
522. ## =============================================================================
523. cat("\n\n========== SECTION 3: Neural correlations (Fig 4; Supp Table 4) ==========\n")
525. rois <- c(r_DL_sw="rDLPFC", l_DL_sw="lDLPFC", r_RL_sw="rRLPFC", l_RL_sw="lRLPFC")
526. outcomes <- c(sEBR="sEBR", sw_acc="Switching accuracy")
527. forest_df <- purrr::map_dfr(names(rois), function(r) purrr::map_dfr(names(outcomes), function(ov) {
528. raw <- sp(raw_data[[r]], raw_data[[ov]])
529. adj <- pspear(raw_data[[r]], raw_data[[ov]], raw_data["age"])
530. tibble(Panel = rois[[r]], Series = outcomes[[ov]],
531. rho_raw = raw["rho"], p_raw = raw["p"], n_raw = raw["n"],
532. ci_low_raw = raw["ci_low"], ci_high_raw = raw["ci_high"],
533. rho_age = adj["rho"], p_age = adj["p"], n_age = adj["n"],
534. ci_low_age = adj["ci_low"], ci_high_age = adj["ci_high"])
535. }))
536. print(as.data.frame(forest_df %>% mutate(across(starts_with("rho"), ~round(.,3)),
537. across(starts_with("p_"), ~signif(.,3)),
538. across(starts_with("ci_"), ~round(.,3)))))
539. write.csv(forest_df, file.path(FIGDIR, "Fig4_SuppTable4_correlations.csv"), row.names = FALSE)
540. for (i in seq_len(nrow(forest_df))) {
541. fr <- forest_df[i,]
542. add_row("Fig4/SuppTable4", paste0(fr$Panel," ~ ", fr$Series), round(fr$rho_raw,3),
543. "Spearman rho","", fmtp(fr$p_raw), fr$n_raw,
544. "", "", round(fr$ci_low_raw,3), round(fr$ci_high_raw,3))
545. if (fr$Panel == "rDLPFC")
546. add_row("Fig4", paste0(fr$Panel," ~ ", fr$Series, " | age"), round(fr$rho_age,3),
547. "partial Spearman","", fmtp(fr$p_age), fr$n_age,
548. "", "", round(fr$ci_low_age,3), round(fr$ci_high_age,3))
549. }
551. ## 3-2 Mediation: sEBR -> rDLPFC -> switching accuracy [Supp Fig 5]
552. med_df <- raw_data %>% dplyr::select(sEBR, sw_acc, r_DL_sw) %>% drop_na()
553. cat("\n-- Mediation (Supp Fig 5): N =", nrow(med_df), "--\n")
554. q75 <- quantile(med_df$sEBR, 0.75); q25 <- quantile(med_df$sEBR, 0.25)
555. ## The sEBR "effect" in this mediation is the interquartile-range increase, so
556. ## the two quartiles define the contrast being estimated and are reported in the
557. ## Supplementary Information. Written to the summary CSV as their own rows.
558. cat(sprintf("sEBR interquartile range used as the treatment contrast: %.2f -> %.2f blinks/min\n",
559. q25, q75))
560. add_row("SuppFig5","sEBR 25th percentile (control value)", round(unname(q25),2), "percentile",
561. "","", nrow(med_df), "", "blinks/min; lower end of the interquartile-range contrast")
562. add_row("SuppFig5","sEBR 75th percentile (treatment value)", round(unname(q75),2), "percentile",
563. "","", nrow(med_df), "", "blinks/min; upper end of the interquartile-range contrast")
564. model.m <- lm(r_DL_sw ~ sEBR, data = med_df)
565. model.y <- lm(sw_acc ~ sEBR + r_DL_sw, data = med_df)
566. set.seed(365)
567. med.out <- mediate(model.m, model.y, treat = "sEBR", mediator = "r_DL_sw",
568. treat.value = q75, control.value = q25,
569. boot = FALSE, sims = 5000, robustSE = TRUE)
570. print(summary(med.out))
571. ms <- summary(med.out)
572. add_row("SuppFig5","ACME (indirect via rDLPFC)", round(ms$d0,4), "quasi-Bayes","", fmtp(ms$d0.p), nrow(med_df), "", "exploratory",
573. round(ms$d0.ci[1],4), round(ms$d0.ci[2],4))
574. add_row("SuppFig5","ADE (direct)", round(ms$z0,4), "quasi-Bayes","", fmtp(ms$z0.p), nrow(med_df), "", "exploratory",
575. round(ms$z0.ci[1],4), round(ms$z0.ci[2],4))
576. add_row("SuppFig5","Proportion mediated", round(ms$n0,4), "quasi-Bayes","", "", nrow(med_df), "",
577. "exploratory; kept to 4 dp because the manuscript quotes this as a percentage to 2 dp",
578. round(ms$n0.ci[1],4), round(ms$n0.ci[2],4))
579. sens.out <- medsens(med.out, rho.by = 0.1, sims = 1000)
580. print(summary(sens.out))
581. ## Sensitivity to violations of sequential ignorability. `err.cr.d` is the value
582. ## of the residual correlation rho at which the ACME is driven to zero -- the
583. ## quantity quoted in the Supplementary Information. It has length 2 only when
584. ## the ACME differs between treatment and control groups, which is not the case
585. ## here; the guard keeps the row well defined either way. Note that the grid is
586. ## stepped by rho.by = 0.1, so this threshold is resolved to that granularity.
587. rho0 <- unname(sens.out$err.cr.d)
588. add_row("SuppFig5","Sensitivity: rho at which ACME = 0",
589. paste(round(rho0, 2), collapse = "; "), "rho",
590. "","", nrow(med_df), "",
591. paste0("residual correlation between the mediator and outcome errors at which ",
592. "the ACME reaches zero; evaluated on a grid with rho.by = 0.1"))
595. ## =============================================================================
596. ## SECTION 4 -- Supp Note 1 (sex/site/trials), Supp Note 4 (outliers),
597. ## Methods (low pre-switch accuracy exclusion)
598. ## =============================================================================
599. cat("\n\n========== SECTION 4: Robustness / sensitivity ==========\n")
601. ## 4-1 sex / site adjusted associations [Supp Note 1]
602. pairs <- list(c("sEBR","r_DL_sw","sEBR-rDLPFC"),
603. c("sEBR","sw_acc","sEBR-sw_acc"),
604. c("r_DL_sw","sw_acc","rDLPFC-sw_acc"))
605. for (cv in c("sex","site")) {
606. if (all(is.na(raw_data[[cv]]))) next
607. cat(sprintf("\n-- adjust for %s --\n", cv))
608. for (pr in pairs) {
609. r <- pspear(raw_data[[pr[1]]], raw_data[[pr[2]]], raw_data[cv])
610. cat(sprintf("%-14s rho=%.2f p=%s n=%d\n", pr[3], r["rho"], fmtp(r["p"]), r["n"]))
611. add_row("SuppNote1", paste0(pr[3]," | ", cv), round(r["rho"],3), "partial Spearman","", fmtp(r["p"]), r["n"],
612. "", "", round(r["ci_low"],3), round(r["ci_high"],3))
613. }
614. }
616. ## 4-2 LME robustness: add sex / site / phase-matched trial count [Supp Note 1]
617. dlx <- df_long %>% mutate(sex = raw_data$sex[match(Subject, raw_data$ID)],
618. site = raw_data$site[match(Subject, raw_data$ID)])
619. fit_lrt <- function(extra, sub) {
620. d <- dlx %>% drop_na(OxyHb, age, all_of(sub))
621. m0 <- lmer(OxyHb ~ (Phase*Hemisphere*Region)*age + (1|Subject), d, REML = FALSE)
622. m1 <- lmer(as.formula(paste0("OxyHb ~ (Phase*Hemisphere*Region)*age + ", extra, " + (1|Subject)")),
623. d, REML = FALSE)
624. an <- anova(m0, m1); list(chi = an$Chisq[2], df = an$Df[2], p = an$`Pr(>Chisq)`[2])
625. }
626. for (s in list(c("sex","sex"), c("site","site"))) {
627. r <- fit_lrt(s[1], s[2]); cat(sprintf("LME +%s: chi2(%d)=%.3f p=%s\n", s[1], r$df, r$chi, fmtp(r$p)))
628. add_row("SuppNote1", paste0("LME age vs age+", s[1]), round(r$chi,3), "chi2", r$df, fmtp(r$p), "", "age")
629. }
630. dl2 <- df_long %>% drop_na(OxyHb, age, PhaseRep)
631. m_a <- lmer(OxyHb ~ (Phase*Hemisphere*Region)*age + (1|Subject), dl2, REML = FALSE)
632. m_r <- lmer(OxyHb ~ (Phase*Hemisphere*Region)*age + PhaseRep + (1|Subject), dl2, REML = FALSE)
633. m_rm <- lmer(OxyHb ~ (Phase*Hemisphere*Region)*age + PhaseRep*age + (1|Subject), dl2, REML = FALSE)
634. # Step 1: does the PhaseRep main effect improve fit over the age model? (1 df)
635. # Step 2: does the PhaseRep x age interaction improve fit over the model that
636. # ALREADY contains the PhaseRep main effect? (1 df)
637. # Step 2 must compare m_r with m_rm, not m_a with m_rm: the latter would test the
638. # main effect and the interaction jointly on 2 df. The comparison below isolates
639. # the interaction on 1 df, as the label states.
640. lr1 <- anova(m_a, m_r) # m_a -> m_a + PhaseRep
641. lr2 <- anova(m_r, m_rm) # m_r -> m_r + PhaseRep:age
642. stopifnot(lr1$Df[2] == 1, lr2$Df[2] == 1) # both steps must be 1 df
643. cat(sprintf("LME +trial: chi2(%d)=%.3f p=%s | +trial x age: chi2(%d)=%.3f p=%s\n",
644. lr1$Df[2], lr1$Chisq[2], fmtp(lr1$`Pr(>Chisq)`[2]),
645. lr2$Df[2], lr2$Chisq[2], fmtp(lr2$`Pr(>Chisq)`[2])))
646. add_row("SuppNote1","LME age vs age+PhaseRep", round(lr1$Chisq[2],3),"chi2",lr1$Df[2],
647. fmtp(lr1$`Pr(>Chisq)`[2]),"","age+rep",
648. "LRT of the PhaseRep main effect (m_a vs m_r)")
649. add_row("SuppNote1","LME PhaseRep x age", round(lr2$Chisq[2],3),"chi2",lr2$Df[2],
650. fmtp(lr2$`Pr(>Chisq)`[2]),"","age+rep",
651. "LRT of the PhaseRep x age interaction over the main-effect model (m_r vs m_rm)")
652. # Reference only: the joint 2-df test of both terms together, retained for
653. # traceability. The 1-df test above is the one reported.
654. lr2_joint <- anova(m_a, m_rm)
655. cat(sprintf(" [reference, joint 2-df test] +trial and trial x age: chi2(%d)=%.3f p=%s\n",
656. lr2_joint$Df[2], lr2_joint$Chisq[2], fmtp(lr2_joint$`Pr(>Chisq)`[2])))
657. add_row("SuppNote1","LME PhaseRep + PhaseRep x age (joint)", round(lr2_joint$Chisq[2],3),
658. "chi2", lr2_joint$Df[2], fmtp(lr2_joint$`Pr(>Chisq)`[2]),"","age+rep",
659. "joint 2-df test of both terms; reference only, superseded by the 1-df test above")
660. arp <- anova(m_r)
661. for (eff in c("Phase:age","Hemisphere:age","Region:age"))
662. add_row("SuppNote1", paste0(eff," (with trial)"), round(arp[eff,"F value"],2),"F",
663. paste0(arp[eff,"NumDF"],",",round(arp[eff,"DenDF"],0)), fmtp(arp[eff,"Pr(>F)"]),"","age+rep")
665. ## Effect sizes for the trial-adjusted model and for the model comparisons.
666. es_rep <- as.data.frame(eta_squared(anova(m_r, type = 3), partial = TRUE,
667. ci = 0.95, alternative = "two.sided"))
668. for (eff in c("Phase:age","Hemisphere:age","Region:age")) {
669. i <- match(eff, es_rep$Parameter)
670. add_row("SuppNote1", paste0("partial eta2 ", eff, " (with trial)"),
671. signif(es_rep$Eta2_partial[i],3), "eta2_p","","","", "age+rep","",
672. round(es_rep$CI_low[i],3), round(es_rep$CI_high[i],3))
673. }
674. for (s in c("sex","site")) {
675. if (all(is.na(raw_data[[s]]))) next
676. d <- dlx %>% drop_na(OxyHb, age, all_of(s))
677. q0 <- lmer(OxyHb ~ (Phase*Hemisphere*Region)*age + (1|Subject), d, REML = FALSE)
678. q1 <- lmer(as.formula(paste0("OxyHb ~ (Phase*Hemisphere*Region)*age + ", s, " + (1|Subject)")),
679. d, REML = FALSE)
680. dm <- r2_nakagawa(q1)$R2_marginal - r2_nakagawa(q0)$R2_marginal
681. cat(sprintf("Change in marginal R2 when adding %s: %.4f\n", s, dm))
682. add_row("SuppNote1", paste0("Change in marginal R2 with ", s), round(dm,4), "dR2","","","","age")
683. }
684. dm_rep <- r2_nakagawa(m_r)$R2_marginal - r2_nakagawa(m_a)$R2_marginal
685. cat(sprintf("Change in marginal R2 when adding trial count: %.4f\n", dm_rep))
686. add_row("SuppNote1","Change in marginal R2 with trial count", round(dm_rep,4),
687. "dR2","","","","age+rep")
689. ## 4-3 rep-corrected partial correlations [Supp Note 1]
690. for (pr in list(c("sEBR","r_DL_sw","sEBR-rDLPFC"), c("sEBR","sw_acc","sEBR-sw_acc"))) {
691. for (z in list(list("sw_rep","| sw_rep"), list(c("age","sw_rep"),"| age+sw_rep"))) {
692. r <- pspear(raw_data[[pr[1]]], raw_data[[pr[2]]], raw_data[z[[1]]])
693. add_row("SuppNote1", paste0(pr[3]," ", z[[2]]), round(r["rho"],3), "partial Spearman","", fmtp(r["p"]), r["n"],
694. "", "", round(r["ci_low"],3), round(r["ci_high"],3))
695. }
696. }
697. ## 4-4 completed-trial descriptives + age x trials [Methods / Supp Note 1]
698. desc <- raw_data %>% summarise(across(c(Pre_rep, Post_rep, Mix_rep),
699. list(m = ~mean(.,na.rm=TRUE), s = ~sd(.,na.rm=TRUE))))
700. print(desc)
701. for (ph in c("Pre_rep","Post_rep","Mix_rep")) {
702. add_row("Methods/SuppNote1", paste0("completed trials ", ph),
703. sprintf("%.1f+/-%.1f", mean(raw_data[[ph]],na.rm=TRUE), sd(raw_data[[ph]],na.rm=TRUE)), "desc","","",nrow(raw_data))
704. ra <- sp(raw_data$age, raw_data[[ph]])
705. add_row("SuppNote1", paste0("age ~ ", ph), round(ra["rho"],3), "Spearman rho","", fmtp(ra["p"]), ra["n"],
706. "", "", round(ra["ci_low"],3), round(ra["ci_high"],3))
707. }
709. ## 4-5 Outlier robustness [Supp Note 4]
710. raw_data$sEBR_w <- winsorize(raw_data$sEBR, 0.05)
711. thr <- mean(raw_data$sEBR, na.rm = TRUE) + 3*sd(raw_data$sEBR, na.rm = TRUE)
712. trim <- raw_data %>% filter(is.na(sEBR) | sEBR <= thr)
713. for (z in list(list("sEBR_w","sw_acc","winsorized","sEBR-sw_acc"),
714. list("sEBR_w","r_DL_sw","winsorized","sEBR-rDLPFC"))) {
715. r <- sp(raw_data[[z[[1]]]], raw_data[[z[[2]]]])
716. add_row("SuppNote4", paste0(z[[4]]," ", z[[3]]), round(r["rho"],3),"Spearman rho","", fmtp(r["p"]), r["n"],
717. "", "", round(r["ci_low"],3), round(r["ci_high"],3))
718. }
719. for (z in list(list("sw_acc","trimmed(n-3)","sEBR-sw_acc"), list("r_DL_sw","trimmed(n-3)","sEBR-rDLPFC"))) {
720. r <- sp(trim$sEBR, trim[[z[[1]]]])
721. add_row("SuppNote4", paste0(z[[3]]," ", z[[2]]), round(r["rho"],3),"Spearman rho","", fmtp(r["p"]), r["n"],
722. "", "", round(r["ci_low"],3), round(r["ci_high"],3))
723. }
725. ## 4-6 Low pre-switch accuracy exclusion [Methods]
726. excl <- raw_data %>% filter(pre_acc <= 0.5)
727. cat("\nExcluded (pre_acc<=0.5): n =", nrow(excl), "| IDs:", paste(excl$ID, collapse=","), "\n")
728. de <- raw_data %>% filter(pre_acc > 0.5)
729. for (z in list(list("sEBR","sw_acc","sEBR-sw_acc"), list("sEBR","r_DL_sw","sEBR-rDLPFC"),
730. list("age","sw_acc","age-sw_acc"), list("age","sEBR","age-sEBR"))) {
731. r <- sp(de[[z[[1]]]], de[[z[[2]]]])
732. add_row("Methods", paste0("after excl: ", z[[3]]), round(r["rho"],3),"Spearman rho","", fmtp(r["p"]), r["n"],
733. "", "", round(r["ci_low"],3), round(r["ci_high"],3))
734. }
737. ## =============================================================================
738. ## SECTION 5 -- Write summary CSVs
739. ## =============================================================================
740. cat("\n\n========== SECTION 5: Writing summaries ==========\n")
741. summary_df <- dplyr::bind_rows(SUMMARY)
742. write.csv(summary_df, file.path(OUTDIR, "results_summary_generated.csv"), row.names = FALSE)
743. cat("Wrote:", file.path(OUTDIR, "results_summary_generated.csv"),
744. "(", nrow(summary_df), "rows )\n")
745. cat("Source-data CSVs in:", FIGDIR, "\n")
748. ## ---- session information -----------------------------------------------
749. ## Written to outputs_summary/sessionInfo.txt as well as echoed to the console.
750. ## The script name below is the single source of truth for that file's header;
751. ## keep it in sync with the filename if this script is revised again.
752. SCRIPT_NAME <- "MainAnalysis.R"
753. si_path <- file.path(OUTDIR, "sessionInfo.txt")
754. si_con <- file(si_path, open = "wt")
755. writeLines(c(
756. paste0(SCRIPT_NAME, " -- session information"),
757. paste0("Generated: ", format(Sys.time(), "%Y-%m-%d %H:%M:%S %Z")),
758. paste0("Data file: ", data_path),
759. ""), si_con)
760. capture.output(sessionInfo(), file = si_con)
761. close(si_con)
762. cat("Wrote:", si_path, "\n")
764. cat("\n=== Session info ===\n"); print(sessionInfo())
765. # ==============================================================================
766. # END
767. # ==============================================================================

MainAnalysis.R, no license · at the source