OSCR

Decoding everyday levels of musical training from subcortical white-matter architecture.

Code ↔ Paper

7 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 7 matches · 2 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Results › Sample characteristics and analytical approach ↔ scripts/preprocessing/unpack_connectome_log_transformation.m, lines 51–91 · score 0.65 · log transformed, outcome variable, eITV, sex, structural connectivity, age
  2. [2] § Methods › Mediation and moderation models ↔ scripts/mediation/mediation_moderation.R, lines 154–227 · score 0.61 · full model, BCa, delta, bootstrapped, mediation, indirect
  3. [3] § Methods › Network-level inference ↔ scripts/analysis/network_statistics.py, lines 1903–1988 · score 0.56 · confidence intervals, prediction strength, network pairs, brain networks, FDR, permutation
  4. [4] § Methods › dMRI preprocessing and structural brain network construction ↔ mrtrix_pipeline_QC.sh, the whole file · a weak match · score 0.55 · mrtrix3, orientation, intensity, diffusion
  5. [5] § Methods › Mediation and moderation models ↔ scripts/mediation/mediation_moderation.R, lines 1–43 · score 0.55 · SubC, FIML, exogenous, SEM, mediation, moderation
  6. [6] § Methods › dMRI preprocessing and structural brain network construction ↔ mrtrix_pipeline_step_3.sh, the whole file · a weak match · score 0.54 · backtracking, cutoff, dynamic, FOD, seeding, algorithm
  7. [7] § Results › Sample characteristics and analytical approach ↔ scripts/analysis/network_statistics.py, lines 2147–2286 · score 0.51 · NBS Predict, brain networks, Destrieux, pipeline, structural connectivity, log

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 · 629 lines · 23 KB · no license · 2 matches

  1. # ============================================================
  2. # Subcortical (between) connectivity): Mediation vs Moderation
  3. # ------------------------------------------------------------
  4. # Model 1 (Mediation): F3 → SubC → MET (with IQ, Age, Gender covariates)
  5. # Model 1R (No-mediator): Remove MET ~ SubC for a true nested test (Δdf=1)
  6. # Model 2 (Moderation): MET ~ F3 + SubC + F3×SubC + IQ + Age + Gender
  7. # ------------------------------------------------------------
  8. # Key changes vs previous:
  9. # - F3 treated as EXOGENOUS (no regression on covariates)
  10. # - True nested mediation comparison (remove the b-path)
  11. # - Consistent preprocessing (z-score continuous; Gender as 0/1)
  12. # - Proper MI pooling with semTools::runMI (SEM) and mice::with/pool (lm)
  13. # - Report bootstrapped CIs for indirect when using FIML; pooled CIs otherwise
  14. # ============================================================
  16. suppressPackageStartupMessages({
  17. library(lavaan)
  18. library(semPlot)
  19. library(semTools) # runMI, lavTestLRT.mi, parameter pooling
  20. library(psych)
  21. library(ggplot2)
  22. library(dplyr)
  23. library(tidyr)
  24. library(corrplot)
  25. library(ggpubr)
  26. library(VIM) # missing data visualization
  27. library(mice) # multiple imputation
  28. library(interactions) # interact_plot (OLS moderation)
  29. library(lmtest) # coeftest
  30. library(sandwich) # vcovHC
  31. })
  33. # ------------------------
  34. # 0) Paths & I/O
  35. # ------------------------
  36. setwd(".")
  37. output_dir <- "results/mediation"
  38. if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)
  40. # ------------------------
  41. # 1) Load data
  42. # ------------------------
  43. dat <- read.csv("results/connectivity_features/connectivity_mediation_features.csv")
  45. # ------------------------
  46. # 2) Select & rename
  47. # ------------------------
  48. dat0 <- dat %>%
  49. select(Age, Gender, TotalMETScore, F3, IQ, subcortical_between_mean_conn) %>%
  50. rename(
  51. MET = TotalMETScore,
  52. F3_S = F3, # lifetime training
  53. SubC = subcortical_between_mean_conn
  54. )
  56. # ------------------------
  57. # 3) Gender handling (binary 0/1)
  58. # ------------------------
  59. if (!is.numeric(dat0$Gender)) {
  60. if (is.factor(dat0$Gender) || is.character(dat0$Gender)) {
  61. lev <- unique(as.character(dat0$Gender))
  62. if (length(lev) == 2) {
  63. lev <- sort(lev)
  64. map <- setNames(c(0, 1), lev)
  65. dat0$Gender_num <- as.integer(map[as.character(dat0$Gender)])
  66. cat("Gender mapping (0/1):\n")
  67. print(map)
  68. } else {
  69. stop("Gender has more than 2 levels. Please recode to binary (0/1) before running.")
  70. }
  71. } else {
  72. stop("Gender is not numeric/factor/character. Please provide binary 0/1.")
  73. }
  74. } else {
  75. # Assume already 0/1 (prints summary for sanity)
  76. dat0$Gender_num <- dat0$Gender
  77. cat("Gender is numeric. Summary:\n")
  78. print(summary(dat0$Gender_num))
  79. }
  81. # ------------------------
  82. # 4) Missingness overview
  83. # ------------------------
  84. miss_tbl <- data.frame(
  85. Variable = c("Age", "Gender_num", "MET", "F3_S", "IQ", "SubC"),
  86. N_Missing = c(
  87. sum(is.na(dat0$Age)),
  88. sum(is.na(dat0$Gender_num)),
  89. sum(is.na(dat0$MET)),
  90. sum(is.na(dat0$F3_S)),
  91. sum(is.na(dat0$IQ)),
  92. sum(is.na(dat0$SubC))
  93. )
  94. )
  95. miss_tbl$Pct_Missing <- round(100 * miss_tbl$N_Missing / nrow(dat0), 1)
  96. write.csv(miss_tbl, file.path(output_dir, "missingness_summary.csv"), row.names = FALSE)
  97. print(miss_tbl)
  99. pdf(file.path(output_dir, "missing_data_pattern_SubC.pdf"), width = 8, height = 6)
  100. aggr(dat0[, c("Age","MET","F3_S","IQ","SubC")],
  101. col = c('navyblue','red'), numbers = TRUE, sortVars = TRUE,
  102. main = "Missing Data Pattern - SC-SM Analysis")
  103. dev.off()
  105. # Decision: impute if IQ missingness > 20%
  106. iq_missing_pct <- mean(is.na(dat0$IQ))*100
  107. use_imputation <- iq_missing_pct > 20
  108. cat("IQ missing %:", round(iq_missing_pct,1), "| Using multiple imputation:", use_imputation, "\n")
  110. # ------------------------
  111. # 5) Scale continuous variables (z-scores). Do NOT scale Gender.
  112. # ------------------------
  113. scale01 <- function(x) as.numeric(scale(x))
  114. datZ <- dat0 %>%
  115. mutate(
  116. Age_z = scale01(Age),
  117. MET_z = scale01(MET),
  118. F3_z = scale01(F3_S),
  119. IQ_z = scale01(IQ),
  120. SubC_z= scale01(SubC)
  121. )
  123. # ------------------------
  124. # 6) Descriptives & correlations
  125. # ------------------------
  126. desc <- psych::describe(datZ[, c("Age_z","MET_z","F3_z","IQ_z","SubC_z")])
  127. write.csv(round(desc, 3), file.path(output_dir, "descriptives.csv"))
  128. cor_mat <- cor(datZ[, c("MET","F3_S","IQ","Age","SubC")], use = "pairwise.complete.obs")
  129. write.csv(round(cor_mat, 3), file.path(output_dir, "correlation_matrix.csv"))
  130. pdf(file.path(output_dir, "correlation_matrix_SubC.pdf"), width = 7, height = 7)
  131. corrplot(cor_mat, method = "color", type = "upper", order = "hclust",
  132. tl.col = "black", tl.srt = 45, addCoef.col = "black", number.cex = .7)
  133. dev.off()
  135. # ------------------------
  136. # 7) Mediation models (F3 EXOGENOUS)
  137. # Full: SubC_z ~ a*F3_z + covs ; MET_z ~ b*SubC_z + cprime*F3_z + covs
  138. # No-M: SubC_z ~ a*F3_z + covs ; MET_z ~ ctot*F3_z + covs
  139. # Exogenous covariances among F3_z, IQ_z, Age_z, Gender_num
  140. # ------------------------
  142. # --- Model specifications ---
  143. # CRITICAL: Constrain MET_z ~~ 0*SubC_z in BOTH models to prevent
  144. # the reduced model from regaining the b-path via residual covariance
  146. model_med_full <- '
  147. # mediator & outcome regressions
  148. SubC_z ~ a*F3_z + IQ_z + Age_z + Gender_num
  149. MET_z ~ b*SubC_z + cprime*F3_z + IQ_z + Age_z + Gender_num
  151. # Constrain residual covariance to zero
  152. MET_z ~~ 0*SubC_z
  154. # effects
  155. indirect := a*b
  156. direct := cprime
  157. total := indirect + direct
  158. '
  160. model_med_nomedi <- '
  161. # mediator regression (same as full model)
  162. SubC_z ~ a*F3_z + IQ_z + Age_z + Gender_num
  164. # outcome regression WITHOUT b-path (SubC_z removed)
  165. MET_z ~ ctot*F3_z + IQ_z + Age_z + Gender_num
  167. # CRITICAL: Constrain residual covariance to zero (same as full model)
  168. MET_z ~~ 0*SubC_z
  169. '
  171. # ------------------------
  172. # 8) Fit mediation (FIML or MI pooled)
  173. # ------------------------
  174. set.seed(123)
  176. if (!use_imputation) {
  177. # ---- FIML + bootstrap CIs for indirect
  178. fit_full <- sem(
  179. model_med_full, data = datZ,
  180. missing = "fiml",
  181. meanstructure = TRUE,
  182. fixed.x = TRUE, # Treat exogenous vars as fixed
  183. auto.cov.y = FALSE, # Don't auto-add endogenous covariances
  184. se = "bootstrap", bootstrap = 5000
  185. )
  186. fit_nom <- sem(
  187. model_med_nomedi, data = datZ,
  188. missing = "fiml",
  189. meanstructure = TRUE,
  190. fixed.x = TRUE, # Treat exogenous vars as fixed
  191. auto.cov.y = FALSE, # Don't auto-add endogenous covariances
  192. se = "bootstrap", bootstrap = 5000
  193. )
  195. # LRT (true nested: removes b-path => Δdf=1)
  196. med_LRT <- anova(fit_full, fit_nom)
  197. write.csv(as.data.frame(med_LRT), file.path(output_dir, "mediation_LRT.csv"), row.names = FALSE)
  199. # Parameter estimates with bootstrapped CIs
  200. pe_full <- parameterEstimates(fit_full, standardized = TRUE,
  201. boot.ci.type = "bca.simple", ci = TRUE)
  202. write.csv(pe_full, file.path(output_dir, "mediation_param_estimates_full.csv"), row.names = FALSE)
  204. } else {
  205. # ---- Multiple Imputation
  206. # Build a mids object for variables used in SEM
  207. imp_vars <- datZ %>% select(MET_z, F3_z, SubC_z, IQ_z, Age_z, Gender_num)
  208. # Basic imputation spec (predictors default); m=20 for stability
  209. imp <- mice(imp_vars, m = 20, maxit = 20, printFlag = FALSE, seed = 123)
  211. fit_full <- runMI(model_med_full, fun = "sem", data = imp, se = "standard",
  212. estimator = "MLR", meanstructure = TRUE,
  213. fixed.x = TRUE, auto.cov.y = FALSE)
  214. fit_nom <- runMI(model_med_nomedi, fun = "sem", data = imp, se = "standard",
  215. estimator = "MLR", meanstructure = TRUE,
  216. fixed.x = TRUE, auto.cov.y = FALSE)
  218. # MI-pooled LRT (D3 method)
  219. med_LRT <- lavTestLRT.mi(fit_full, fit_nom)
  220. capture.output(med_LRT, file = file.path(output_dir, "mediation_LRT_MI.txt"))
  222. # MI-pooled parameter estimates (delta-method CIs)
  223. pe_full <- summary(fit_full, standardized = TRUE, ci = TRUE)
  224. capture.output(pe_full, file = file.path(output_dir, "mediation_param_estimates_full_MI.txt"))
  226. # NOTE: Bootstrapped CIs under MI require bootstrapLavaan.mi; omitted for brevity.
  227. }
  229. # ------------------------
  230. # 9) Extract mediation effects
  231. # ------------------------
  232. get_effect <- function(pe, label) {
  233. if (is.data.frame(pe)) {
  234. row <- pe[pe$label == label, , drop = FALSE]
  235. if (nrow(row) == 0) return(c(NA, NA, NA, NA))
  236. return(c(row$est, row$ci.lower, row$ci.upper, row$pvalue))
  237. } else {
  238. # pe is a "summary" print under MI; re-extract via parameterEstimates if possible
  239. # For MI, grab parameterEstimates from the fitted object:
  240. pe_tab <- tryCatch(parameterEstimates(fit_full, standardized = TRUE, ci = TRUE),
  241. error = function(e) NULL)
  242. if (is.null(pe_tab)) return(c(NA, NA, NA, NA))
  243. row <- pe_tab[pe_tab$label == label, , drop = FALSE]
  244. if (nrow(row) == 0) return(c(NA, NA, NA, NA))
  245. return(c(row$est, row$ci.lower, row$ci.upper, row$pvalue))
  246. }
  247. }
  248. indirect_vec <- get_effect(pe_full, "indirect")
  249. direct_vec <- get_effect(pe_full, "direct")
  250. total_vec <- get_effect(pe_full, "total")
  252. # ------------------------
  253. # 10) Moderation (OLS)
  254. # Compare base vs interaction to get ΔR² and test
  255. # If MI: pool coefficients and average ΔR² across imputations.
  256. # ------------------------
  257. form_base <- as.formula(MET_z ~ F3_z + IQ_z + SubC_z + Age_z + Gender_num)
  258. form_int <- as.formula(MET_z ~ F3_z + IQ_z + SubC_z + F3_z:SubC_z + Age_z + Gender_num)
  260. if (!use_imputation) {
  261. mod_data <- datZ[, all.vars(form_int)]
  262. cc <- complete.cases(mod_data)
  263. mod_data <- mod_data[cc, ]
  265. fit_base <- lm(form_base, data = mod_data)
  266. fit_int <- lm(form_int, data = mod_data)
  268. # ΔR² and ANOVA
  269. R2_base <- summary(fit_base)$r.squared
  270. R2_int <- summary(fit_int)$r.squared
  271. dR2 <- R2_int - R2_base
  272. anv <- anova(fit_base, fit_int)
  274. # Robust SEs (HC3) for interaction term
  275. robust_coefs <- coeftest(fit_int, vcov. = vcovHC(fit_int, type = "HC3"))
  276. int_row <- robust_coefs[rownames(robust_coefs) == "F3_z:SubC_z", , drop = FALSE]
  278. write.csv(broom::tidy(fit_int), file.path(output_dir, "moderation_lm_tidy.csv"), row.names = FALSE)
  279. write.csv(as.data.frame(anv), file.path(output_dir, "moderation_lm_anova_base_vs_int.csv"))
  280. robust_df <- data.frame(
  281. term = rownames(robust_coefs),
  282. estimate = robust_coefs[, "Estimate"],
  283. std.error = robust_coefs[, "Std. Error"],
  284. statistic = robust_coefs[, "t value"],
  285. p.value = robust_coefs[, "Pr(>|t|)"],
  286. row.names = NULL,
  287. check.names = FALSE
  288. )
  289. write.csv(robust_df, file.path(output_dir, "moderation_lm_robustHC3.csv"), row.names = FALSE)
  292. } else {
  293. # Moderation with MI pooling
  294. # Use the same mids (imp); fit base and interaction in each imputed dataset
  295. fit_base_mi <- with(data = imp, expr = lm(form_base))
  296. fit_int_mi <- with(data = imp, expr = lm(form_int))
  297. pool_int <- pool(fit_int_mi)
  299. # Extract pooled interaction coefficient
  300. summ_pool <- summary(pool_int, conf.int = TRUE)
  301. int_pool <- summ_pool[summ_pool$term == "F3_z:SubC_z", , drop = FALSE]
  303. # ΔR²: compute per-imputation and average
  304. R2_base_vec <- sapply(fit_base_mi$analyses, function(m) summary(m)$r.squared)
  305. R2_int_vec <- sapply(fit_int_mi$analyses, function(m) summary(m)$r.squared)
  306. dR2 <- mean(R2_int_vec - R2_base_vec, na.rm = TRUE)
  308. # Save
  309. write.csv(summ_pool, file.path(output_dir, "moderation_pooled_coefficients.csv"), row.names = FALSE)
  310. write.csv(data.frame(R2_base = R2_base_vec, R2_int = R2_int_vec,
  311. dR2 = R2_int_vec - R2_base_vec),
  312. file.path(output_dir, "moderation_R2_by_imputation.csv"), row.names = FALSE)
  313. }
  315. # ------------------------
  316. # 10.5) Conditional effects (simple slopes) of F3 at SC–SM = -1, 0, +1 SD
  317. # ------------------------
  319. simple_slopes_out <- file.path(output_dir, "moderation_simple_slopes.csv")
  321. if (!use_imputation) {
  322. # --- Single-model case (use robust HC3 variance) ---
  323. stopifnot(exists("fit_int"), inherits(fit_int, "lm"))
  325. b <- coef(fit_int)
  326. V <- sandwich::vcovHC(fit_int, type = "HC3") # robust variance-covariance
  328. # Helper: slope of F3 at SC level z, its SE using delta method, t, p
  329. slope_at <- function(z) {
  330. # beta_F3(z) = b_F3 + z * b_F3xSC
  331. beta <- unname(b["F3_z"] + z * b["F3_z:SubC_z"])
  332. var <- V["F3_z","F3_z"] +
  333. (z^2)*V["F3_z:SubC_z","F3_z:SubC_z"] +
  334. 2*z*V["F3_z","F3_z:SubC_z"]
  335. se <- sqrt(var)
  336. tval <- beta / se
  337. df <- fit_int$df.residual
  338. pval <- 2*pt(abs(tval) * -1, df = df) # two-sided
  339. c(beta = beta, SE = se, t = tval, p = pval)
  340. }
  342. levels <- c(-1, 0, 1)
  343. labs <- c("Low SC–SM (−1 SD)", "Mean SC–SM (0)", "High SC–SM (+1 SD)")
  344. M <- t(vapply(levels, slope_at, numeric(4L)))
  345. simple_df <- data.frame(
  346. Level = labs,
  347. `β (std.)` = round(M[, "beta"], 3),
  348. SE = round(M[, "SE"], 3),
  349. t = round(M[, "t"], 2),
  350. `p-value` = signif(M[, "p"], 3),
  351. check.names = FALSE
  352. )
  353. write.csv(simple_df, simple_slopes_out, row.names = FALSE)
  354. print(simple_df)
  356. } else {
  357. # --- MI case: compute per-imputation slopes then pool with Rubin's rules ---
  358. stopifnot(exists("fit_int_mi"))
  360. # Extract per-imputation coefficients and (model-based) vcov
  361. analyses <- fit_int_mi$analyses
  362. m <- length(analyses)
  363. if (m < 2) stop("Not enough imputations to pool simple slopes.")
  365. # Compute per-imputation Q (slope) and U (its variance) at each z
  366. levels <- c(-1, 0, 1)
  367. labs <- c("Low SC–SM (−1 SD)", "Mean SC–SM (0)", "High SC–SM (+1 SD)")
  369. pool_one_level <- function(z) {
  370. Q <- numeric(m) # estimates
  371. U <- numeric(m) # variances
  373. for (i in seq_len(m)) {
  374. bi <- coef(analyses[[i]])
  375. Vi <- vcov(analyses[[i]]) # model-based vcov; robust + MI is non-trivial
  377. beta <- unname(bi["F3_z"] + z * bi["F3_z:SubC_z"])
  378. var <- Vi["F3_z","F3_z"] +
  379. (z^2)*Vi["F3_z:SubC_z","F3_z:SubC_z"] +
  380. 2*z*Vi["F3_z","F3_z:SubC_z"]
  382. Q[i] <- beta
  383. U[i] <- var
  384. }
  386. # Rubin's rules
  387. Qbar <- mean(Q)
  388. Ubar <- mean(U)
  389. B <- var(Q) # between-imputation variance
  390. Tvar <- Ubar + (1 + 1/m) * B
  391. SE <- sqrt(Tvar)
  393. # df (Barnard & Rubin)
  394. lam <- (1 + 1/m) * B / Tvar
  395. df_old <- (m - 1) / (lam^2)
  396. # approximate df for Ubar (complete-data) – use average residual df
  397. df_complete <- mean(vapply(analyses, function(mod) mod$df.residual, numeric(1)))
  398. df <- 1 / ( (1/df_old) + (1/df_complete) )
  400. tval <- Qbar / SE
  401. pval <- 2 * pt(abs(tval) * -1, df = df)
  403. c(beta = Qbar, SE = SE, t = tval, p = pval)
  404. }
  406. M <- t(vapply(levels, pool_one_level, numeric(4L)))
  407. simple_df <- data.frame(
  408. Level = labs,
  409. `β (std.)` = round(M[, "beta"], 3),
  410. SE = round(M[, "SE"], 3),
  411. t = round(M[, "t"], 2),
  412. `p-value` = signif(M[, "p"], 3),
  413. check.names = FALSE
  414. )
  415. write.csv(simple_df, simple_slopes_out, row.names = FALSE)
  416. print(simple_df)
  417. }
  419. # Also write a short summary row with R²_base, R²_int, and ΔR² (if available)
  420. r2_summary_path <- file.path(output_dir, "moderation_model_comparison.csv")
  421. if (exists("R2_base") && exists("R2_int")) {
  422. write.csv(
  423. data.frame(
  424. R2_main_effects = round(R2_base, 3),
  425. R2_with_interaction = round(R2_int, 3),
  426. delta_R2 = round(R2_int - R2_base, 3)
  427. ),
  428. r2_summary_path,
  429. row.names = FALSE
  430. )
  431. } else if (exists("R2_base_vec") && exists("R2_int_vec")) {
  432. write.csv(
  433. data.frame(
  434. R2_main_effects_mean = round(mean(R2_base_vec), 3),
  435. R2_with_interaction_mean = round(mean(R2_int_vec), 3),
  436. delta_R2_mean = round(mean(R2_int_vec - R2_base_vec), 3)
  437. ),
  438. r2_summary_path,
  439. row.names = FALSE
  440. )
  441. }
  443. # ------------------------
  444. # 11) Visuals
  445. # ------------------------
  446. # Mediation path diagram
  447. pdf(file.path(output_dir, "sem_paths_mediation_full.pdf"), width = 10, height = 7)
  448. semPaths(fit_full, what = "std", layout = "tree2",
  449. edge.label.cex = .7, sizeMan = 7, sizeLat = 9,
  450. residuals = FALSE, nCharNodes = 0,
  451. title = TRUE, title.cex = 1.1)
  452. dev.off()
  454. # ---- Mediation effects plot: ALWAYS create a page ----
  455. # ---- Mediation effects plot: ALWAYS create a page ----
  456. pdf(file.path(output_dir, "mediation_effects_CI.pdf"), width = 8, height = 5)
  458. pe_ok <- (exists("pe_full") &&
  459. is.data.frame(pe_full) &&
  460. any(pe_full$label %in% c("indirect","direct","total")))
  462. if (pe_ok) {
  463. eff_tab <- subset(pe_full, label %in% c("indirect","direct","total"),
  464. select = c(label, est, ci.lower, ci.upper))
  465. # coerce to numeric
  466. eff_tab$est <- as.numeric(eff_tab$est)
  467. eff_tab$ci.lower <- as.numeric(eff_tab$ci.lower)
  468. eff_tab$ci.upper <- as.numeric(eff_tab$ci.upper)
  470. # ORDER: Direct, Indirect, Total
  471. eff_tab$label <- factor(eff_tab$label,
  472. levels = c("direct","indirect","total"),
  473. labels = c("Direct (c′)", "Indirect (a×b)", "Total"))
  475. g <- ggplot(eff_tab, aes(x = label, y = est, fill = label)) +
  476. geom_col() +
  477. geom_errorbar(aes(ymin = ci.lower, ymax = ci.upper), width = .2) +
  478. geom_hline(yintercept = 0, linetype = "dashed") +
  479. scale_fill_manual(values = c(
  480. "Direct (c′)" = "#d9d9d9", # light gray
  481. "Indirect (a×b)" = "#f2f2f2", # lightest gray
  482. "Total" = "#7a7a7a" # darker gray
  483. ), guide = "none") +
  484. theme_minimal() +
  485. labs(title = "Mediation effects with 95% CI",
  486. x = NULL, y = "Standardized effect (β)")
  487. print(g)
  488. } else {
  489. plot.new()
  490. text(.5, .6, "Mediation effects not available", cex = 1.2)
  491. text(.5, .45, "Expected labels in pe_full: indirect, direct, total.", cex = .9)
  492. }
  493. dev.off()
  495. # ---- Moderation interaction plot: ALWAYS create a page ----
  496. # ---- Moderation interaction plot: ALWAYS create a page ----
  497. pdf(file.path(output_dir, "moderation_interaction_plot.pdf"), width = 8, height = 6)
  499. if (!use_imputation && exists("fit_int") && inherits(fit_int, "lm")) {
  500. # ensure mod_data exists
  501. if (!exists("mod_data")) {
  502. mod_data <- tryCatch(model.frame(fit_int), error = function(e) NULL)
  503. }
  505. if (!is.null(mod_data)) {
  506. at_sc <- c(-1, 0, 1)
  507. x_min <- min(mod_data$F3_z, na.rm = TRUE)
  508. x_max <- max(mod_data$F3_z, na.rm = TRUE)
  510. # prediction grid for lines
  511. grid <- expand.grid(
  512. F3_z = seq(x_min, x_max, length.out = 120),
  513. SubC_z = at_sc,
  514. IQ_z = 0,
  515. Age_z = 0,
  516. Gender_num = 0
  517. )
  518. grid$MET_hat <- predict(fit_int, newdata = grid)
  519. grid$SC_lab <- factor(grid$SubC_z, levels = at_sc,
  520. labels = c("SC–SM: -1 SD", "SC–SM: 0", "SC–SM: +1 SD"))
  522. # group points to match the 3 SC–SM reference levels
  523. mod_data$SC_lab <- cut(mod_data$SubC_z,
  524. breaks = c(-Inf, -1, 1, Inf),
  525. labels = c("SC–SM: -1 SD", "SC–SM: 0", "SC–SM: +1 SD"),
  526. right = TRUE)
  528. g2 <- ggplot() +
  529. # observed points
  530. geom_point(data = mod_data, aes(x = F3_z, y = MET_z, shape = SC_lab),
  531. alpha = 0.35, size = 1, na.rm = TRUE) +
  532. # predicted lines
  533. geom_line(data = grid, aes(x = F3_z, y = MET_hat, linetype = SC_lab)) +
  534. theme_minimal() +
  535. labs(title = "Moderation: Simple slopes of F3 at SC–SM levels",
  536. x = "F3 (z)", y = "MET (z)",
  537. linetype = NULL, shape = NULL)
  538. print(g2)
  539. } else {
  540. plot.new()
  541. text(.5, .6, "Interaction model not available for plotting", cex = 1.2)
  542. text(.5, .45, "Could not reconstruct model.frame(fit_int).", cex = .9)
  543. }
  544. } else {
  545. plot.new()
  546. text(.5, .6, "Interaction model not available for plotting", cex = 1.2)
  547. text(.5, .45, "Check fit_int exists and use_imputation == FALSE.", cex = .9)
  548. }
  549. dev.off()
  553. # ------------------------
  554. # 12) Summaries & comparison table
  555. # ------------------------
  556. # LRT p-value
  557. if (!use_imputation) {
  558. LRT_p <- as.numeric(med_LRT$`Pr(>Chisq)`[2])
  559. } else {
  560. # lavTestLRT.mi prints an object; extract p-value if available, else NA
  561. LRT_p <- tryCatch({
  562. as.numeric(med_LRT["Pr(>Chisq)"][2, 1])
  563. }, error = function(e) NA)
  564. }
  566. # Moderation stats
  567. if (!use_imputation) {
  568. int_beta <- as.numeric(coef(fit_int)["F3_z:SubC_z"])
  569. int_p <- summary(fit_int)$coefficients["F3_z:SubC_z", "Pr(>|t|)"]
  570. # robust p (optional):
  571. int_p_HC3 <- tryCatch({
  572. as.numeric(coeftest(fit_int, vcov. = vcovHC(fit_int, type="HC3"))["F3_z:SubC_z","Pr(>|t|)"])
  573. }, error = function(e) NA)
  574. R2b <- R2_base; R2i <- R2_int
  575. } else {
  576. int_beta <- if (nrow(int_pool) == 0) NA else int_pool$estimate
  577. int_p <- if (nrow(int_pool) == 0) NA else int_pool$p.value
  578. int_p_HC3 <- NA
  579. R2b <- mean(R2_base_vec); R2i <- mean(R2_int_vec)
  580. }
  582. comp <- data.frame(
  583. Model = c("Mediation (full vs no-mediator)", "Moderation (interaction)"),
  584. Key_Test = c("LRT for b-path (Δdf=1)", "β(F3×SubC) & ΔR²"),
  585. Estimate_or_Δ = c(
  586. paste0("Indirect (a*b) = ", round(indirect_vec[1], 3),
  587. " [", round(indirect_vec[2], 3), ", ", round(indirect_vec[3], 3), "]"),
  588. paste0("β = ", round(int_beta, 3),
  589. ", ΔR² = ", round(R2i - R2b, 3))
  590. ),
  591. p_value = c(ifelse(is.na(LRT_p), "NA", formatC(LRT_p, format="f", digits=3)),
  592. ifelse(is.na(int_p), "NA", formatC(int_p, format="f", digits=3)))
  593. )
  594. write.csv(comp, file.path(output_dir, "comparison_table.csv"), row.names = FALSE)
  595. print(comp)
  597. # ------------------------
  598. # 13) Text report
  599. # ------------------------
  600. sink(file.path(output_dir, "report_SubC_reimpl.txt"))
  601. cat("SC–SM: Mediation vs Moderation (Re-implementation)\n")
  602. cat("===============================================\n\n")
  603. cat("F3 treated as EXOGENOUS. Covariates (IQ, Age, Gender) enter SubC and MET.\n")
  604. cat("Mediation comparison is TRUE NESTED (remove MET ~ SubC; Δdf=1).\n\n")
  606. cat("Missingness summary (%, N=", nrow(dat0), "):\n", sep = "")
  607. print(miss_tbl); cat("\n")
  608. cat("IQ missing %:", round(iq_missing_pct,1), " | Multiple Imputation:", use_imputation, "\n\n")
  610. cat("Mediation LRT p-value (full vs no-mediator): ", ifelse(is.na(LRT_p), "NA", LRT_p), "\n")
  611. cat("Indirect (a*b): est=", round(indirect_vec[1],3),
  612. " CI=[", round(indirect_vec[2],3), ", ", round(indirect_vec[3],3), "]",
  613. " p=", round(indirect_vec[4],3), "\n\n")
  615. cat("Moderation:\n")
  616. cat(" β(F3×SubC) =", round(int_beta,3),
  617. " | p =", ifelse(is.na(int_p), "NA", round(int_p,3)),
  618. " | ΔR² =", round(R2i - R2b, 3), "\n")
  619. if (!is.na(int_p_HC3)) cat(" (Robust HC3 p ≈", round(int_p_HC3,3), ")\n")
  620. cat("\nNotes:\n")
  621. cat("- Global SEM fit indices are not relied upon because the models are near-saturated.\n")
  622. cat("- Inference for mediation comes from indirect effect CI and nested LRT (Δdf=1).\n")
  623. cat("- Inference for moderation comes from the interaction term and ΔR².\n")
  624. cat("- MI pooling uses semTools::runMI for SEM and mice::with/pool for OLS.\n")
  625. sink()
  627. cat("\n=== DONE ===\nOutput in: ", normalizePath(output_dir), "\n", sep = "")

mediation_moderation.R at commit f21d2ab, no license · at the source

Overview

Authors: Massimo Lumaca1, Marcus T. Pearce1,2, Peter E. Keller1,3, Peter Vuust1, Elvira Brattico1,4, Giosuè Baggio5, Katarzyna Hat6, Ole A. Heggli1, Kristian Sandberg7,8
ORCID iDs: Massimo Lumaca, Marcus T. Pearce, Peter E. Keller, Elvira Brattico, Giosuè Baggio, Katarzyna Hat, Ole A. Heggli, Kristian Sandberg
  1. Center for Music in the Brain, Department of Clinical Medicine, Aarhus University & The Royal Academy of Music Aarhus/Aalborg, Aarhus, Denmark
  2. School of Electronic Engineering and Computer Science, Queen Mary University of London, London, United Kingdom
  3. The MARCS Institute for Brain, Behaviour and Development, Western Sydney University, Penrith, Australia
  4. Department of Education, Psychology, Communication, University of Bari, Bari, Italy
  5. Language Acquisition and Language Processing Lab, Norwegian University of Science and Technology, Trondheim, Norway
  6. Consciousness Lab, Institute of Psychology, Jagiellonian University, Kraków, Poland
  7. Center of Functionally Integrative Neuroscience, Department of Clinical Medicine, Aarhus University, Aarhus, Denmark
  8. Neurobiology Research Unit, Copenhagen University Hospital Rigshospitalet, Copenhagen, Denmark
Institutions: Aarhus University (Denmark); Royal Academy of Music (Denmark); Queen Mary University of London (United Kingdom); Western Sydney University (Australia); University of Bari Aldo Moro (Italy); Norwegian University of Science and Technology (Norway); Jagiellonian University (Poland); Rigshospitalet (Denmark)
Journal: Imaging neuroscience (Cambridge, Mass.), volume 4, article IMAG.a.1325
Dates: received 25 February 2026; accepted 8 July 2026; published online 4 August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1162/imag.a.1325 · PMID 42559177 · PMCID PMC13440155 · OpenAlex W7168242511
Open access: diamond, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), systems (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, fMRI & imaging, Physiology & signal measures
Keywords: network neuroscience, machine learning, musical training, neuroplasticity, subcortical circuits, non-musicians
Topic: Neuroscience and Music Perception (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: European Cooperation in Science and Technology (COST) (CA18106); Center for Music in the Brain (MIB) is funded by The Lundbeck Foundation (R469-2024-1573); Center for Music in the Brain (MIB) is funded by Købmand Herman Sallings Fond
Citations: not cited yet (Europe PMC); 176 references in the paper
Research resources: RRID:SCR_001847, 2017 RRID:SCR_002438, 2001 RRID:SCR_002823, RRID:SCR_008796, 2019 RRID:SCR_016216

Abstract

Most people engage with music informally or through the standard school curriculum rather than through extensive professional practice, resulting in graded levels of accumulated musical training. In this study, we used cross-validated machine learning on whole-brain structural connectomes to investigate whether individual differences in brain wiring predict variation in musical training across 225 adults with none-to-moderate training levels. Connectomes were weighted by fiber bundle capacity (FBC), a quantitative diffusion MRI metric of structural connectivity. Musical training was quantified using the Gold-MSI Musical Training subscale (F3), a self-reported measure encompassing both formal and informal musical practice. Subcortical white-matter pathways linking thalamus, putamen, and pallidum with sensorimotor cortex carried the strongest predictive signal (r = 0.261; R2 = 0.068; p_perm = 0.003). Furthermore, this connectivity predicted the association between musical training and auditory perception skill, which was present in individuals with stronger subcortical–sensorimotor connectivity (+1 SD: β = 0.192, p = 0.004), but absent in those with weaker connectivity (−1 SD: β = −0.066, p = 0.488). Subcortical–sensorimotor connectivity is, therefore, the strongest structural correlate of graded musical training in non-specialists within the general population, identifying the structural neural substrate on which the training–skill association depends. This work warrants renewed targeted focus on subcortical circuits in models of musical expertise.

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

Repositories

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

MassimoLumaca/neuroARC

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: cc7dd46694fc0d44559bd9f2eacc16550512d4a3, 28 February 2024
Languages: Shell (8)
Size: 12 files, 8 scripts
Software Heritage: not archived
Found in: “Data and Code Availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: MRtrix3 (7 files), FreeSurfer (2 files), FSL (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
9 files
At the source: github.com/MassimoLumaca/neuroARC

MassimoLumaca/neuroTraining

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: f21d2abe75d365b4b22f3b18df0ef6e1ade7aab6, 27 January 2026
Languages: Python (3), R (1), MATLAB (1)
Size: 22 files, 5 scripts
Software Heritage: not archived
Found in: “Data and Code Availability”
Holds: README, environment (requirements.txt)
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: Matplotlib (3 files), NumPy (3 files), pandas (3 files), SciPy (3 files), seaborn (3 files), statsmodels (3 files), broom (1 file), ggplot2 (1 file), ggpubr (1 file), lavaan (1 file), Statistics and Machine Learning Toolbox (1 file), psych (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
6 files
At the source: github.com/MassimoLumaca/neuroTraining

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:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 13 scripts, each with its path and the digest of its content;
  • 7 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 and Code Availability

Data cannot be shared publicly as it is part of an ongoing study, and thus considered unanonymized under Danish law, even if pseudonymized. Researchers who wish to access the data may contact Dr Kristian Sandberg () at The Center of Functionally Integrative Neuroscience and/or The Technology Transfer Office () at Aarhus University, Denmark, to establish a data-sharing agreement. After permission has been given by the relevant ethics committee, data will be made available to the researchers for replication purposes. As the project is ongoing, sharing requests for other purposes will be evaluated on a case-by-case basis. Codes for reproducing diffusion-MRI analyses and results are available on GitHub: https://github.com/MassimoLumaca/neuroARC and https://github.com/MassimoLumaca/neuroTraining.

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, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 9 authors, 6 keywords, 3 funders, 173 references, 5 RRIDs.

Cite

This paper

Lumaca, M., Pearce, M. T., Keller, P. E., Vuust, P., Brattico, E., Baggio, G., Hat, K., Heggli, O. A., & Sandberg, K. (2026). Decoding everyday levels of musical training from subcortical white-matter architecture. Imaging neuroscience (Cambridge, Mass.), 4, IMAG.a.1325. https://doi.org/10.1162/imag.a.1325

BibTeX

@article{lumaca2026decoding,
author = {Lumaca, Massimo and Pearce, Marcus T. and Keller, Peter E. and Vuust, Peter and Brattico, Elvira and Baggio, Giosuè and Hat, Katarzyna and Heggli, Ole A. and Sandberg, Kristian},
title = {{Decoding everyday levels of musical training from subcortical white-matter architecture}},
journal = {Imaging neuroscience (Cambridge, Mass.)},
year = {2026},
month = aug,
volume = {4},
pages = {IMAG.a.1325},
publisher = {MIT Press},
issn = {2837-6056},
doi = {10.1162/imag.a.1325},
url = {https://doi.org/10.1162/imag.a.1325},
pmid = {42559177},
pmcid = {PMC13440155}
}

RIS

TY - JOUR
AU - Lumaca, Massimo
AU - Pearce, Marcus T.
AU - Keller, Peter E.
AU - Vuust, Peter
AU - Brattico, Elvira
AU - Baggio, Giosuè
AU - Hat, Katarzyna
AU - Heggli, Ole A.
AU - Sandberg, Kristian
TI - Decoding everyday levels of musical training from subcortical white-matter architecture
T2 - Imaging neuroscience (Cambridge, Mass.)
J2 - Imaging Neurosci (Camb)
PY - 2026
DA - 2026/08/04
VL - 4
SP - IMAG.a.1325
SN - 2837-6056
PB - MIT Press
DO - 10.1162/imag.a.1325
UR - https://doi.org/10.1162/imag.a.1325
LA - en
ER -

CSL-JSON

{
"id": "10.1162/imag.a.1325",
"type": "article-journal",
"title": "Decoding everyday levels of musical training from subcortical white-matter architecture",
"container-title": "Imaging neuroscience (Cambridge, Mass.)",
"author": [
{
"family": "Lumaca",
"given": "Massimo"
},
{
"family": "Pearce",
"given": "Marcus T."
},
{
"family": "Keller",
"given": "Peter E."
},
{
"family": "Vuust",
"given": "Peter"
},
{
"family": "Brattico",
"given": "Elvira"
},
{
"family": "Baggio",
"given": "Giosuè"
},
{
"family": "Hat",
"given": "Katarzyna"
},
{
"family": "Heggli",
"given": "Ole A."
},
{
"family": "Sandberg",
"given": "Kristian"
}
],
"container-title-short": "Imaging Neurosci (Camb)",
"volume": "4",
"page": "IMAG.a.1325",
"DOI": "10.1162/imag.a.1325",
"PMID": "42559177",
"PMCID": "PMC13440155",
"ISSN": "2837-6056",
"publisher": "MIT Press",
"URL": "https://doi.org/10.1162/imag.a.1325",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
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.1038/s41467-026-73366-9 [code]
Cortical and white matter myelination proceed in concert during early infancy.
Journal: Nature communications
In common: MRtrix3, FreeSurfer, statsmodels, 5 other tools, 11 references
[2] doi:10.1126/sciadv.aec2348 [code]
Congenital blindness reduces myelination in human visual cortex.
Journal: Science advances
In common: lavaan, FreeSurfer, ggplot2, 3 other tools, structural MRI / diffusion, 13 references
[3] doi:10.1038/s41467-026-73072-6 [code]
Mapping the spatiotemporal continuum of structural connectivity development across the human connectome in youth.
Journal: Nature communications
In common: MRtrix3, psych, FreeSurfer, 5 other tools, structural MRI / diffusion, 10 references
[4] doi:10.1162/imag.a.1153 [code]
TRAMFIX: TRavelling Across Melbourne for FIXel-based analysis (a reproducibility and reliability study).
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: MRtrix3, FSL, statsmodels, 5 other tools, structural MRI / diffusion, 10 references
[5] doi:10.1038/s41467-026-71151-2 [code]
Common and distinct neural correlates of social interaction processing and theory of mind in narratives.
Journal: Nature communications
In common: MRtrix3, FreeSurfer, FSL, 7 other tools, 8 references
[6] doi:10.1016/j.nicl.2026.104012 [code]
Structural-functional multilayer brain network properties and outcome of combined repetitive transcranial magnetic stimulation and psychotherapy for obsessive-compulsive disorder.
Journal: NeuroImage. Clinical
In common: psych, FreeSurfer, ggpubr, 9 other tools, structural MRI / diffusion, 5 references
[7] doi:10.7554/elife.103097 [code]
Canonical neurodevelopmental trajectories of structural and functional manifolds.
Journal: eLife
In common: ggpubr, statsmodels, Statistics and Machine Learning Toolbox, 5 other tools, structural MRI / diffusion, 10 references
[8] doi:10.1038/s41467-026-71918-7 [code]
Developmental disinhibition gates language lateralization in childhood.
Journal: Nature communications
In common: MRtrix3, FreeSurfer, FSL, 5 other tools, 8 references
[9] doi:10.7554/elife.108109 [code]
Multimodal MRI marker of cognition explains the association between cognition and mental health in the UK Biobank.
Journal: eLife
In common: lavaan, psych, ggplot2, 5 other tools, structural MRI / diffusion, 8 references
[10] doi:10.1002/hbm.70562 [code]
Distinct Physiological Mechanisms Drive Grey Matter Plasticity in Complex Versus Simple Sequence Learning.
Journal: Human brain mapping
In common: MRtrix3, statsmodels, seaborn, 4 other tools, structural MRI / diffusion, 6 references

Contribute

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

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

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.