diff --git a/CHANGELOG.md b/CHANGELOG.md index 98b439c..5363f93 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,122 @@ # Changelog +## Unreleased — 2026-09 audit remediation + +A repo audit (2026-09-04) compared the README/CLAUDE.md against the code and +found 24 defects plus a dozen false README claims. All verified and fixed here. + +### Behaviour changes (read these before re-running a study) + +- **Vertex `absolute` band power is now a density (dB/Hz)**, `10*log10(integral / bandwidth)`, + matching the ROI/electrode definition. Previously `vertex_cluster` / `vertex_specparam` + reported `10*log10(integral)`. Within-band group statistics are unaffected (a per-band + constant shift); absolute values and their plots move by `10*log10(bandwidth)` per band. +- **Vertex modules now honour the top-level `epoch_sampling:` block and per-analysis + overrides** (precedence: global → `vertex.epoch_sampling` → analysis block). Previously only + the `vertex:` block was read. `vertex_cluster` now epoch-samples when the merged config + enables it (it never did before). `n_bootstrap: 0` = full timeseries on the vertex sampler + too (it used to fall through to a single random draw). +- **`--jobs`**: an explicit CLI value (including `1`) wins over the YAML `jobs:`; when omitted + the YAML value is used. `--jobs 0` / `-1` now actually auto-parallelize (they were coerced + to `1`). `roi_connectivity` and `electrode_connectivity` are parallel-capable too. +- **`--force` now removes previous output** (published `tables/` + `figures/` always; the + working `data/` when the `process` step runs). It used to only bypass `--strict-output`. +- **`vertex_spatial` is retired in Python as well**: no subjects are loaded and R is not + called; it writes empty result tables + a note. (R already did this after Python had + processed every subject.) +- **`roi_connectivity` reads a `metrics:` list** from its config block (like + `vertex_connectivity`). Unknown names error. +- **`roi_directed` hypothesis tables are renamed** to the canonical prefix + (`roi_directed_{global,directed_edges,region}_hypotheses.csv`, `roi_directed_omnibus_lmm.csv`, + `roi_directed_global_bar.png`) and cover DTF as well as TE; DTF-only runs no longer skip R. +- **`fcd_comparison` finds its two primaries across paradigm dirs** (they normally live under + `resting` and `vertex`); `sensor_dir` / `source_dir` overrides are accepted. +- Deprecated analysis names print/check the **canonical** output directory. + +### Fixed + +- Evoked R scripts (`roi_evoked`, `electrode_evoked`) looped `config$contrasts`, which is NULL + under `design:`/`hypotheses:`, so their LMM/post-hoc tables came out empty. They now derive + contrasts from the design spec like `roi_psd_analysis.R`. +- `vertex_evoked` joins the hypothesis contract: `--hypothesis` is accepted and + `vertex_evoked_hypotheses.csv` is written via the permutation adapter; `list` shows an + "Evoked Response (Vertex Level)" heading. +- AAC / PPC from `roi_cross_freq` now get hypothesis statistics + (`R/roi_cross_freq_edges_analysis.R`: `roi_cross_freq_{aac,ppc}_{global,directed_edges,region}_hypotheses.csv`); + the old log claimed they ran "via the gating path". +- `R/roi_connectivity_analysis.R` no longer requires coherence columns; it adapts to the + metric columns present (`--metric aec` alone works). +- `aggregate_to_regions()` keeps `delta_ref`, so region-level hypotheses work for that DV. +- The package imports without `mne` (lazy TFR import with a clear ImportError); `matplotlib` + is a declared core dependency; extras are documented. +- R scripts are packaged into wheels/sdists (`share/source-analytics/R`); `find_r_script_dir()` + checks that location and honours `SOURCE_ANALYTICS_R_DIR`. +- `init` writes a config that parses as-is: `design:`/`hypotheses:` (omnibus + pairwise + contrasts), canonical bands, one `paradigms:` block, both discovery layouts; `--output -` + streams YAML to stdout. It no longer emits the legacy `contrasts:` form. +- R-timeout log messages now report the real timeout (they said 600 s for 3600 s runs). +- `ANALYSIS_METADATA` `about` text: aperiodic default window is 12–45 Hz (not 2–50); PSD + `absolute` is described as dB/Hz density; `electrode_comparison` / `fcd_comparison` carry + `supplements` + `requires`. +- Helper-script `source()` calls in the R modules are hard failures again (were silently + swallowed by `tryCatch`). +- Removed the orphaned `R/network_analysis.R` and the dead lookups for non-existent + `roi_network_analysis.R` / `vertex_network_analysis.R`; removed the stale + `run_connectivity_network.sh`; added `scripts/run_study.sh` (dependency-ordered recipe). + +### Docs + +- README synced to the code: `init` behaviour, figures off by default, real analytics/results + trees (paradigm + profile segments), install extras, R package list (`ggsignif`, `optparse`; + `stringr` dropped), `vertex_signature` naming, `fcd_comparison` / `electrode_signature` in + the catalog, `electrode_comparison` needing `roi_psd`, `--jobs` / `--profile` / + `--paradigm` semantics, no signed-STC fallback, no `dB` column, referenced-group subject + filtering, `circos_metrics` passthrough, epoch-sampling defaults. +- CLAUDE.md rewritten (it described the original PSD-only package). + +## v0.6.0 — 2026-08 + +- `fcd_comparison`: source-vs-sensor functional-connectivity-density comparison module. +- Evoked build-out: ERP amplitude/latency measures, induced power, debiased ITC, cycle ramps and + tiled extraction, wired into all three evoked modules; declared hypotheses for `roi_evoked` / + `electrode_evoked` with the measure as the FDR facet. +- `vertex_specparam`: two-fit peak detection, fit-window diagnostic, persisted `offset_centered`. +- Aperiodic default fit window 12–45 Hz (cited in `docs/methods/APERIODIC_FIT_WINDOW.md`); + `vertex_specparam` no longer fits through the line-noise notch. +- `electrode_signature` module; `vertex_mvpa` renamed to `vertex_signature` (multi-model neural + signature, true AUC, valid permutation p, balanced accuracy). +- `delta_ref` (delta-referenced power) DV for `roi_psd` / `electrode_psd` under a profile. +- `--profile` runs via `StudyConfig.for_profile()` writing to `analytics//` and + `results//`; profile-narrowed hypotheses forwarded to R. +- `--jobs` parallel per-subject processing for the vertex modules, `roi_connectivity` and + `electrode_connectivity`; precomputed-connectivity cache for `vertex_graph`. +- Fully-qualified `fdr_family` (member-set identity); canonical low→high band order everywhere. +- ROI PSD `absolute` switched to power density (dB/Hz) with a restricted relative-power range. +- Figures: module figure dir cleared before regeneration; figures regenerable from persisted + data; effect-size mosaics; anatomical-coverage labels for significant clusters; NBS + subnetwork figures; circos polish. +- NBS significant-edge mask is now filled (it was allocated and never written). +- Packaging: `specparam` floor made resolvable (`>=2.0.0rc6`); version single-sourced from git. + +## v0.5.0 — 2026-06 + +- Declarative hypothesis layer: `design:` / `hypotheses:` (kinds `omnibus`, `contrast`, + `regression`, `equivalence`) with declarative FDR family scope and method; emmeans (R) + tabular adapter, Python permutation (map + cluster) adapter, edge/NBS adapter and + directed-edge adapter; `--hypothesis NAME` selection. Auto-gating retired. +- `config.contrasts` derived from the design spec (stored bridge removed); a legacy + `contrasts:` block is lifted into the spec. +- Module renames: `roi_pac` → `roi_cross_freq` (PAC + AAC + PPC), `roi_transfer_entropy` → + `roi_directed` (TE + DTF); connectivity network split into `*_graph` + `*_nbs`. +- New modules: `vertex_cross_freq`, `vertex_directed` (ridge-MVAR DTF), `vertex_evoked`, + `electrode_connectivity` (source-vs-sensor FC comparator). +- New kernels: wPLI, dPLI, AAC, n:m PPC with surrogate significance, Hipp-2012 AEC + (vectorized), directed-aware FCD. +- Per-sub-output selection: `--metric` / `--band` / `--select`. +- `vertex_spatial` GLS statistics retired (R side). +- `ANALYSIS_METADATA` domains + `supplements`; MIT license; README rewritten around the + source-localization handoff. + ## v0.4.0 — 2026-04-25 ### Compatibility diff --git a/CLAUDE.md b/CLAUDE.md index 98542b5..dffaedc 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -1,77 +1,74 @@ # CLAUDE.md — source-analytics -## What This Is +Group-level statistics for source-localized EEG. **Python** orchestrates, loads +reconstructions, does signal processing, and runs the permutation/cluster stats +for vertex and sensor maps. **R** (lme4/lmerTest/emmeans, ggplot2) does the LMM +statistics and figures for the ROI/electrode modules. Python calls `Rscript`. -Statistical analysis toolkit for source-localized EEG data. **Python** handles signal processing (PSD, band power) and data I/O. **R** handles statistics (lme4/lmerTest LMMs, t-tests, effect sizes) and visualization (ggplot2). +`README.md` is the user manual and is kept in sync with the code. For methods, +`docs/methods/` is authoritative. When the two disagree with the code, fix the doc. ## Setup ```bash -# Python -cd /home/edm9fd/sandbox/source-analytics uv venv && source .venv/bin/activate -uv pip install -e . - -# R packages (one-time) -Rscript -e 'install.packages(c("ggplot2","dplyr","tidyr","readr","stringr","forcats","lme4","lmerTest","effectsize","emmeans","yaml","argparse","patchwork","scales"))' +uv pip install -e ".[all]" # mne / scikit-learn / networkx / nibabel extras +Rscript -e 'install.packages(c("ggplot2","dplyr","tidyr","readr","forcats","lme4","lmerTest","effectsize","emmeans","yaml","argparse","optparse","patchwork","scales","ggsignif"))' ``` +Tests: `.venv/bin/python -m pytest -q`. R syntax check: +`for f in R/*.R; do Rscript -e "invisible(parse(file='$f'))"; done`. + ## Usage ```bash -source-analytics validate --study /mnt/d/research/EEG/FORGE/analysis.yaml -source-analytics run --study /mnt/d/research/EEG/FORGE/analysis.yaml --analysis psd -source-analytics list +source-analytics init /path/to/rest_roi --name study # writes rest_roi/analysis/study.yaml +source-analytics validate --study study.yaml +source-analytics list --study study.yaml +source-analytics run --study study.yaml --paradigm resting --analysis roi_psd +source-analytics run ... --steps statistics,figures,summary # figures are OFF by default +scripts/run_study.sh study.yaml # whole study, dependency order ``` -## Architecture +## Layout ``` -Python (orchestration + signal processing) R (statistics + visualization) -───────────────────────────────────────── ───────────────────────────── -1. Discovery (find subjects, load YAML) -2. Load ROI timeseries (pickle/.set) -3. Compute PSD (Welch's, scipy) -4. Extract band power (abs, rel, dB) -5. Export CSVs ──────────────────────────► 6. Read CSVs + YAML config - 7. LMMs (lme4/lmerTest) - 8. t-tests, Hedges' g, FDR - 9. ggplot2 figures - 10. Markdown summary +src/source_analytics/ + cli.py run / validate / list / figure / init + core.py ANALYSIS_REGISTRY (+ deprecated aliases), ANALYSIS_METADATA, StudyAnalyzer + config.py StudyConfig, DesignSpec (design:/hypotheses:), profiles, paradigm scoping + analyses/ one module per analysis, all subclass analyses/base.BaseAnalysis + hypothesis/ Python adapters (tabular / permutation / edge) for declared hypotheses + spectral/, stats/, io/, viz/, atlas/ +R/ _analysis.R scripts, hypothesis.R (emmeans adapter), stats_utils.R +scripts/run_study.sh +docs/methods/ APERIODIC_FIT_WINDOW, CONNECTIVITY_METHODS, HYPOTHESIS, DESIGN_SPEC ``` -### Python modules (`src/source_analytics/`) -- `cli.py` — CLI entry point -- `core.py` — StudyAnalyzer orchestrator -- `config.py` — YAML study config loader -- `io/` — Data loading (pkl/npy/.set readers, subject discovery) -- `spectral/` — PSD and band power extraction -- `analyses/` — BaseAnalysis ABC + PSD module (exports CSV, calls Rscript) - -### R scripts (`R/`) -- `psd_analysis.R` — Main entry point (called by Python) -- `stats_utils.R` — Omnibus LMM (group*roi interaction), emmeans post-hoc, FDR/Holm correction -- `plot_psd.R` — ggplot2 figures (PSD curves, boxplots, heatmaps) -- `report.R` — Markdown summary writer - -## CSV Interface (Python → R) - -Exported to `{output_dir}/psd/data/`: -- `band_power.csv` — subject, group, roi, band, absolute, relative, dB -- `psd_curves.csv` — subject, group, roi, freq_hz, psd -- `study_config.yaml` — copy of study config for R - -R outputs to `{output_dir}/psd/tables/`: -- `psd_omnibus.csv` — Omnibus LMM results (group x ROI interaction, Type III ANOVA) -- `psd_posthoc_roi.csv` — emmeans post-hoc contrasts per ROI (Holm-corrected) - -## R packages required - -ggplot2, dplyr, tidyr, readr, stringr, forcats, lme4, lmerTest, effectsize, emmeans, yaml, argparse, patchwork, scales - -## Adding a New Analysis - -1. Create `src/source_analytics/analyses/my_analysis.py` (Python: data extraction + CSV export) -2. Create `R/my_analysis.R` (R: statistics + visualization) -3. Register in `core.py` ANALYSIS_REGISTRY -4. **Update `README.md`** — add the new module to the "Analysis Modules" section with a description of what it does, its Python/R responsibilities, and its output files +Lifecycle per module: `setup → process(_subject) → aggregate → statistics → +figures → summary`. `DEFAULT_RUN_STEPS` excludes `figures`. + +## Output contract + +- Working tree: `paths.analytics/[/]//` holds + `data/` (per-subject CSVs + `study_config.yaml` snapshot) and `ANALYSIS_SUMMARY.md`. +- Published: `paths.results/[/]{tables,figures}///` + (`BaseAnalysis.tbl_dir` / `fig_dir`). source-lightbox reads these. +- Every inferential module writes `_hypotheses.csv` (one row per + band × spatial cell, or per cluster for map modules) into `tables/`. +- Band power CSVs carry `absolute` (mean power density, dB/Hz) and `relative`. + There is no `dB` column. + +## Conventions + +- Deprecated analysis names (`psd`, `pac`, `vertex_mvpa`, …) are in + `core._DEPRECATED_NAMES`; output always goes under the canonical name. +- R scripts derive contrasts with `contrasts_from_spec(parse_design_spec(config))` + from `R/hypothesis.R`; never read `config$contrasts` directly. +- Optional deps are lazy: mne (evoked/TFR), scikit-learn (signature), networkx + (graph). The package must import without them. +- Adding an analysis: subclass `BaseAnalysis`, register in `ANALYSIS_REGISTRY`, + add an `ANALYSIS_METADATA` entry (`category`, `level`, `domain`, and + `supplements`/`requires` for secondaries), declare `SELECTABLE`, add the R + script if it uses emmeans, then add it to the README catalog and + `scripts/run_study.sh`. Keep the timeout value and its log message in step. diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..96935ab --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,3 @@ +graft R +include README.md LICENSE CHANGELOG.md +global-exclude __pycache__ *.pyc diff --git a/R/electrode_aperiodic_analysis.R b/R/electrode_aperiodic_analysis.R index c6bcdbf..27ad503 100644 --- a/R/electrode_aperiodic_analysis.R +++ b/R/electrode_aperiodic_analysis.R @@ -32,8 +32,8 @@ script_dir <- if (exists("script.dir")) { }, error = function(e) "R") } -tryCatch(source(file.path(script_dir, "stats_utils.R")), error = function(e) NULL) -tryCatch(source(file.path(script_dir, "hypothesis.R")), error = function(e) NULL) +source(file.path(script_dir, "stats_utils.R")) +source(file.path(script_dir, "hypothesis.R")) # --- Argument parsing --- parser <- ArgumentParser(description = "Electrode aperiodic analysis (R)") diff --git a/R/electrode_evoked_analysis.R b/R/electrode_evoked_analysis.R index b0ea49b..ce46240 100644 --- a/R/electrode_evoked_analysis.R +++ b/R/electrode_evoked_analysis.R @@ -33,6 +33,7 @@ script_dir <- if (exists("script.dir")) { } source(file.path(script_dir, "stats_utils.R")) +source(file.path(script_dir, "hypothesis.R")) source(file.path(script_dir, "report.R")) # --- Argument parsing --- @@ -51,6 +52,8 @@ parser$add_argument("--figures-only", action = "store_true", default = FALSE, help = "Skip statistics; regenerate figures from existing data/tables") parser$add_argument("--no-figures", action = "store_true", default = FALSE, help = "Skip all figure generation (stats/tables only)") +parser$add_argument("--hypothesis", default = NULL, + help = "Run only the named hypothesis(es) (comma-separated) from the design spec; default = all") args <- parser$parse_args() data_dir <- args$data_dir @@ -79,6 +82,20 @@ group_colors <- unlist(config$group_colors) group_labels <- unlist(config$groups) group_order <- config$group_order +# Pairwise contrast list for the descriptive group*spatial LMM below. Derived from +# the design:/hypotheses: spec (config$contrasts is NULL for modern configs); +# falls back to a legacy contrasts: block if the spec has none. +contrasts <- tryCatch(contrasts_from_spec(parse_design_spec(config)), + error = function(e) config$contrasts) +if (is.null(contrasts)) contrasts <- list() +if (!is.null(args$hypothesis) && length(contrasts) > 0) { + want <- trimws(strsplit(args$hypothesis, ",")[[1]]) + contrasts <- Filter(function(ct) ct$name %in% want, contrasts) + message("Contrasts narrowed to --hypothesis set: ", + paste(vapply(contrasts, function(ct) ct$name, character(1)), collapse = ", ")) +} +message("Contrasts: ", length(contrasts)) + # Check for electrode_categories in config electrode_categories <- config$electrode_categories @@ -94,7 +111,7 @@ if (!figures_only) { all_omnibus <- list() all_posthoc <- list() -for (contrast in config$contrasts) { +for (contrast in contrasts) { cname <- contrast$name ga <- contrast$group_a gb <- contrast$group_b @@ -194,7 +211,7 @@ if (nrow(omnibus_df) > 0) { # --- Post-hoc emmeans (gated on significant omnibus) --- message("\nRunning post-hoc emmeans...") -for (contrast in config$contrasts) { +for (contrast in contrasts) { cname <- contrast$name ga <- contrast$group_a gb <- contrast$group_b @@ -295,7 +312,7 @@ if (length(electrode_categories) > 0) { all_omnibus_reg <- list() all_posthoc_reg <- list() - for (contrast in config$contrasts) { + for (contrast in contrasts) { cname <- contrast$name ga <- contrast$group_a gb <- contrast$group_b @@ -379,7 +396,7 @@ if (length(electrode_categories) > 0) { } # Region post-hoc - for (contrast in config$contrasts) { + for (contrast in contrasts) { cname <- contrast$name ga <- contrast$group_a gb <- contrast$group_b @@ -459,12 +476,12 @@ for (mname in measure_names) { # Rename channel to roi for stats_utils compatibility gp_data <- measures_df %>% filter(measure_name == mname, - group %in% unlist(lapply(config$contrasts, function(c) c(c$group_a, c$group_b)))) %>% + group %in% unlist(lapply(contrasts, function(c) c(c$group_a, c$group_b)))) %>% rename(roi = channel) if (nrow(gp_data) == 0) next - gp <- run_posthoc_global(gp_data, config$contrasts, spatial_col = "roi", + gp <- run_posthoc_global(gp_data, contrasts, spatial_col = "roi", dv_col = "value", dv_label = mname) if (nrow(gp) > 0) global_posthoc_list[[mname]] <- gp } diff --git a/R/network_analysis.R b/R/network_analysis.R deleted file mode 100644 index c4df3f6..0000000 --- a/R/network_analysis.R +++ /dev/null @@ -1,190 +0,0 @@ -#!/usr/bin/env Rscript -# Network Analysis Report Generator -# Reads graph metrics and NBS results, generates ANALYSIS_SUMMARY.md - -suppressPackageStartupMessages({ - library(optparse) - library(yaml) -}) - -option_list <- list( - make_option("--data-dir", type = "character", help = "Path to data/ directory"), - make_option("--config", type = "character", help = "Path to study_config.yaml"), - make_option("--output-dir", type = "character", help = "Path to output directory"), - make_option("--fig-dir", type = "character", default = NULL, - help = "Directory for figures (default: output-dir/figures)"), - make_option("--tbl-dir", type = "character", default = NULL, - help = "Directory for tables (default: output-dir/tables)"), - make_option("--no-figures", action = "store_true", default = FALSE, - help = "Skip all figure generation") -) -opts <- parse_args(OptionParser(option_list = option_list)) - -no_figures <- isTRUE(opts[["no-figures"]]) - -if (no_figures) { - ggsave <- function(...) invisible(NULL) -} - -data_dir <- opts[["data-dir"]] -config_path <- opts[["config"]] -output_dir <- opts[["output-dir"]] - -config <- read_yaml(config_path) - -fig_dir <- if (!is.null(opts[["fig-dir"]])) opts[["fig-dir"]] else file.path(output_dir, "figures") -tbl_dir <- if (!is.null(opts[["tbl-dir"]])) opts[["tbl-dir"]] else file.path(output_dir, "tables") - -# --- Load data ---------------------------------------------------------------- -global_path <- file.path(data_dir, "network_global_metrics.csv") -nodal_path <- file.path(data_dir, "network_nodal_metrics.csv") - -if (!file.exists(global_path)) { - cat("No network_global_metrics.csv found.\n") - quit(status = 0) -} - -global_metrics <- read.csv(global_path, stringsAsFactors = FALSE) - -# --- Config ------------------------------------------------------------------- -net_cfg <- config$network %||% list() -threshold_method <- net_cfg$threshold_method %||% "proportional" -threshold_value <- net_cfg$threshold_value %||% 0.1 -nbs_threshold <- net_cfg$nbs_threshold %||% 3.0 -nbs_perm <- net_cfg$nbs_permutations %||% 5000 - -# --- Summaries ---------------------------------------------------------------- -n_subjects <- length(unique(global_metrics$subject)) -bands <- unique(global_metrics$band) -groups <- unique(global_metrics$group) - -# Group-level global metric summaries -global_summary <- do.call(rbind, lapply(bands, function(b) { - do.call(rbind, lapply(groups, function(g) { - sub <- global_metrics[global_metrics$band == b & global_metrics$group == g, ] - data.frame( - Band = b, - Group = g, - Mean_Efficiency = round(mean(sub$global_efficiency), 4), - Mean_Modularity = round(mean(sub$modularity), 4), - Mean_SW = round(mean(sub$small_worldness), 3), - Mean_Edges = round(mean(sub$n_edges), 1), - stringsAsFactors = FALSE - ) - })) -})) - -# T-tests on global metrics -global_tests <- do.call(rbind, lapply(bands, function(b) { - sub <- global_metrics[global_metrics$band == b, ] - if (length(groups) < 2) return(NULL) - g1 <- sub[sub$group == groups[1], ] - g2 <- sub[sub$group == groups[2], ] - - do.call(rbind, lapply(c("global_efficiency", "modularity", "small_worldness"), function(m) { - tryCatch({ - tt <- t.test(g1[[m]], g2[[m]]) - data.frame( - Band = b, - Metric = m, - t = round(tt$statistic, 3), - p = round(tt$p.value, 4), - stringsAsFactors = FALSE - ) - }, error = function(e) NULL) - })) -})) - -# --- Write ANALYSIS_SUMMARY.md ----------------------------------------------- -lines <- c( - "# Network Analysis Summary", - "", - sprintf("**Study**: %s", config$name), - "**Analysis**: Graph-theoretic metrics + Network-Based Statistic (NBS)", - sprintf("**Threshold**: %s (%.2f)", threshold_method, threshold_value), - sprintf("**NBS threshold**: t = %.1f", nbs_threshold), - sprintf("**NBS permutations**: %d", nbs_perm), - sprintf("**Subjects**: %d (%s)", n_subjects, paste(groups, collapse = ", ")), - "", - "## Methods", - "", - "Graph metrics (degree, clustering, betweenness, efficiency, modularity,", - "small-worldness) were computed from thresholded imaginary coherence matrices.", - "Group differences in nodal metrics: cluster-based permutation testing.", - "Subnetwork identification: Network-Based Statistic (Zalesky et al., 2010).", - "" -) - -# Epoch info -wb_cfg <- config$vertex %||% list() -epoch_cfg <- wb_cfg$epoch_sampling -if (!is.null(epoch_cfg) && isTRUE(epoch_cfg$enabled)) { - lines <- c(lines, - sprintf("**Epoch sampling**: %d epochs of %.1fs", - epoch_cfg$n_epochs, epoch_cfg$epoch_duration_sec), - "" - ) -} - -lines <- c(lines, - "## Global Metrics Summary", - "", - "| Band | Group | Efficiency | Modularity | Small-World | Edges |", - "|------|-------|------------|------------|-------------|-------|" -) - -for (i in seq_len(nrow(global_summary))) { - r <- global_summary[i, ] - lines <- c(lines, sprintf( - "| %s | %s | %.4f | %.4f | %.3f | %.1f |", - r$Band, r$Group, r$Mean_Efficiency, r$Mean_Modularity, r$Mean_SW, r$Mean_Edges - )) -} - -# Global metric t-tests -if (!is.null(global_tests) && nrow(global_tests) > 0) { - lines <- c(lines, - "", - "## Global Metric Group Comparisons", - "", - "| Band | Metric | t | p |", - "|------|--------|---|---|" - ) - for (i in seq_len(nrow(global_tests))) { - r <- global_tests[i, ] - lines <- c(lines, sprintf("| %s | %s | %.3f | %.4f |", r$Band, r$Metric, r$t, r$p)) - } -} - -# NBS results -nbs_path <- file.path(tbl_dir, "nbs_results.csv") -if (file.exists(nbs_path)) { - nbs <- read.csv(nbs_path, stringsAsFactors = FALSE) - sig_nbs <- nbs[nbs$p_corrected < 0.05, ] - lines <- c(lines, "", "## NBS Results", "") - if (nrow(sig_nbs) > 0) { - lines <- c(lines, sprintf("**%d significant subnetworks** (p < 0.05):", nrow(sig_nbs))) - for (i in seq_len(nrow(sig_nbs))) { - r <- sig_nbs[i, ] - lines <- c(lines, sprintf("- %s component %d: %d edges, p = %.4f", - r$key, r$component, r$n_edges, r$p_corrected)) - } - } else { - lines <- c(lines, "No significant NBS subnetworks at p < 0.05.") - } -} - -lines <- c(lines, - "", - "## Output Files", - "", - "- `data/network_nodal_metrics.csv` — per-vertex graph metrics", - "- `data/network_global_metrics.csv` — global graph metrics", - "- `tables/network_stats.csv` — cluster permutation on nodal metrics", - "- `tables/nbs_results.csv` — NBS subnetwork results", - "- `figures/network_*.png` — nodal metric glass brains", - "" -) - -writeLines(lines, file.path(output_dir, "ANALYSIS_SUMMARY.md")) -cat("Wrote ANALYSIS_SUMMARY.md\n") diff --git a/R/roi_aperiodic_analysis.R b/R/roi_aperiodic_analysis.R index a42a6f2..44b6e98 100644 --- a/R/roi_aperiodic_analysis.R +++ b/R/roi_aperiodic_analysis.R @@ -36,9 +36,9 @@ script_dir <- if (exists("script.dir")) { } # Source helpers -tryCatch(source(file.path(script_dir, "stats_utils.R")), error = function(e) NULL) -tryCatch(source(file.path(script_dir, "hypothesis.R")), error = function(e) NULL) -tryCatch(source(file.path(script_dir, "plot_psd.R")), error = function(e) NULL) +source(file.path(script_dir, "stats_utils.R")) +source(file.path(script_dir, "hypothesis.R")) +source(file.path(script_dir, "plot_psd.R")) if (!exists("theme_pub")) { theme_pub <- function(base_size = 14) { diff --git a/R/roi_connectivity_analysis.R b/R/roi_connectivity_analysis.R index ebaa4d5..424d49e 100644 --- a/R/roi_connectivity_analysis.R +++ b/R/roi_connectivity_analysis.R @@ -33,7 +33,18 @@ script_dir <- if (exists("script.dir")) { }, error = function(e) "R") } -tryCatch(source(file.path(script_dir, "stats_utils.R")), error = function(e) NULL) +source(file.path(script_dir, "stats_utils.R")) + +# Connectivity metrics the Python side can export. Every step below iterates the +# subset actually present in the edge CSV, so a `--metric aec` (or any other +# single-metric) run works end to end. +KNOWN_METRICS <- c("coherence", "imag_coherence", "pli", "wpli", "dwpli", "dpli", + "aec", "partial_corr") +METRIC_LABELS <- c(coherence = "Coherence", imag_coherence = "Imag. Coherence", + pli = "PLI", wpli = "wPLI", dwpli = "dwPLI", dpli = "dPLI", + aec = "AEC", partial_corr = "Partial Corr.") +present_metrics <- function(df) intersect(KNOWN_METRICS, names(df)) +metric_label <- function(m) ifelse(is.na(METRIC_LABELS[m]), m, METRIC_LABELS[m]) # Define sig_stars locally if not sourced if (!exists("sig_stars")) { @@ -64,17 +75,15 @@ theme_pub <- function(base_size = 14) { # =========================================================================== #' Compute global connectivity per subject x band (mean of all edges) -#' @param edges data.frame with columns: subject, group, band, roi1, roi2, coherence, imag_coherence -#' @return data.frame with subject, group, band, mean_coherence, mean_imag_coherence +#' @param edges data.frame with columns: subject, group, band, roi1, roi2, plus one +#' column per connectivity metric present (any of KNOWN_METRICS) +#' @return data.frame with subject, group, band, mean_ per metric present compute_global_connectivity <- function(edges) { + mets <- present_metrics(edges) edges %>% group_by(subject, group, band) %>% summarise( - mean_coherence = mean(coherence, na.rm = TRUE), - mean_imag_coherence = mean(imag_coherence, na.rm = TRUE), - mean_pli = if ("pli" %in% names(.)) mean(pli, na.rm = TRUE) else NA_real_, - mean_dwpli = if ("dwpli" %in% names(.)) mean(dwpli, na.rm = TRUE) else NA_real_, - mean_aec = if ("aec" %in% names(.)) mean(aec, na.rm = TRUE) else NA_real_, + across(all_of(mets), ~ mean(.x, na.rm = TRUE), .names = "mean_{.col}"), n_edges = n(), .groups = "drop" ) @@ -87,21 +96,11 @@ compute_global_connectivity <- function(edges) { #' @param bands named list of band limits #' @return data.frame with t-test results run_global_ttests <- function(global_df, contrasts, bands) { - metrics <- c("mean_coherence", "mean_imag_coherence") - metric_labels <- c("coherence", "imag_coherence") - # Include optional metrics if present in data - if ("mean_pli" %in% names(global_df) && !all(is.na(global_df$mean_pli))) { - metrics <- c(metrics, "mean_pli") - metric_labels <- c(metric_labels, "pli") - } - if ("mean_dwpli" %in% names(global_df) && !all(is.na(global_df$mean_dwpli))) { - metrics <- c(metrics, "mean_dwpli") - metric_labels <- c(metric_labels, "dwpli") - } - if ("mean_aec" %in% names(global_df) && !all(is.na(global_df$mean_aec))) { - metrics <- c(metrics, "mean_aec") - metric_labels <- c(metric_labels, "aec") - } + # Every mean_ column present (and not all-NA) is tested. + metric_labels <- KNOWN_METRICS[paste0("mean_", KNOWN_METRICS) %in% names(global_df)] + metric_labels <- metric_labels[vapply(metric_labels, function(m) + !all(is.na(global_df[[paste0("mean_", m)]])), logical(1))] + metrics <- paste0("mean_", metric_labels) results <- list() for (contrast in contrasts) { @@ -183,7 +182,8 @@ run_global_ttests <- function(global_df, contrasts, bands) { # =========================================================================== #' Map edges to region pairs, average within -#' @param edges data.frame with columns: subject, group, band, roi1, roi2, coherence, imag_coherence +#' @param edges data.frame with columns: subject, group, band, roi1, roi2, plus the +#' metric columns present #' @param roi_categories named list of ROI name vectors #' @return data.frame with region_pair replacing roi1/roi2 aggregate_edges_to_region_pairs <- function(edges, roi_categories) { @@ -210,14 +210,11 @@ aggregate_edges_to_region_pairs <- function(edges, roi_categories) { ) # Average edge values within each subject x band x region_pair + mets <- present_metrics(edges) edges_mapped %>% group_by(subject, group, band, region_pair) %>% summarise( - coherence = mean(coherence, na.rm = TRUE), - imag_coherence = mean(imag_coherence, na.rm = TRUE), - pli = if ("pli" %in% names(.)) mean(pli, na.rm = TRUE) else NA_real_, - dwpli = if ("dwpli" %in% names(.)) mean(dwpli, na.rm = TRUE) else NA_real_, - aec = if ("aec" %in% names(.)) mean(aec, na.rm = TRUE) else NA_real_, + across(all_of(mets), ~ mean(.x, na.rm = TRUE)), n_edges = n(), .groups = "drop" ) @@ -228,7 +225,7 @@ aggregate_edges_to_region_pairs <- function(edges, roi_categories) { #' @param region_pair_df data.frame from aggregate_edges_to_region_pairs() #' @param contrasts list of contrast definitions #' @param bands named list of band limits -#' @param metric character: "coherence" or "imag_coherence" +#' @param metric character: one of the metric columns present (see KNOWN_METRICS) #' @return data.frame with omnibus results run_omnibus_lmm_region_pair <- function(region_pair_df, contrasts, bands, metric = "coherence") { if (!has_lme4) { @@ -332,7 +329,7 @@ run_omnibus_lmm_region_pair <- function(region_pair_df, contrasts, bands, metric #' @param contrasts list of contrast definitions #' @param bands named list of band limits #' @param omnibus_df data.frame from run_omnibus_lmm_region_pair() -#' @param metric character: "coherence" or "imag_coherence" +#' @param metric character: one of the metric columns present (see KNOWN_METRICS) #' @param gate logical: if TRUE, only run for significant omnibus results #' @return data.frame with post-hoc results run_posthoc_emmeans_region_pair <- function(region_pair_df, contrasts, bands, @@ -424,8 +421,8 @@ run_posthoc_emmeans_region_pair <- function(region_pair_df, contrasts, bands, # =========================================================================== #' Plot group-mean connectivity matrices (heatmaps) per band -#' @param edges data.frame with subject, group, band, roi1, roi2, coherence, imag_coherence -#' @param metric_col column name: "coherence" or "imag_coherence" +#' @param edges data.frame with subject, group, band, roi1, roi2, plus metric columns +#' @param metric_col column name: one of the metric columns present #' @param group_colors, group_labels, group_order from config #' @param output_dir figures/ directory plot_connectivity_matrices <- function(edges, metric_col, group_colors, @@ -486,19 +483,18 @@ plot_connectivity_matrices <- function(edges, metric_col, group_colors, plot_global_connectivity_bar <- function(global_df, group_colors, group_labels, group_order, output_dir, sig_df = NULL) { - # Pivot to long format for both metrics + # Pivot to long format over every metric present + mean_cols <- paste0("mean_", KNOWN_METRICS) + mean_cols <- mean_cols[mean_cols %in% names(global_df)] plot_data <- global_df %>% filter(group %in% group_order) %>% pivot_longer( - cols = c(mean_coherence, mean_imag_coherence), + cols = all_of(mean_cols), names_to = "metric", values_to = "value" ) %>% mutate( - metric = case_when( - metric == "mean_coherence" ~ "Coherence", - metric == "mean_imag_coherence" ~ "Imag. Coherence" - ), + metric = unname(metric_label(sub("^mean_", "", metric))), group_label = group_labels[group], group_label = factor(group_label, levels = group_labels[group_order]) ) @@ -536,13 +532,7 @@ plot_global_connectivity_bar <- function(global_df, group_colors, group_labels, if (nrow(sig_hits) > 0) { # Map metric names to facet labels sig_hits <- sig_hits %>% - mutate( - metric_facet = case_when( - metric == "coherence" ~ "Coherence", - metric == "imag_coherence" ~ "Imag. Coherence", - TRUE ~ metric - ) - ) + mutate(metric_facet = unname(metric_label(metric))) # Compute y_max per band x metric for positioning y_maxes <- plot_data %>% @@ -660,7 +650,9 @@ write_connectivity_summary <- function(global_df, global_ttest_df, collapse = ", " ) - add("**Analysis:** Functional Connectivity (Coherence & Imaginary Coherence)") + mets <- sub("^mean_", "", grep("^mean_", names(global_df), value = TRUE)) + mets <- intersect(KNOWN_METRICS, mets) + add("**Analysis:** Functional Connectivity (", paste(metric_label(mets), collapse = ", "), ")") add("") add("**Groups:** ", group_str) add("") @@ -668,7 +660,16 @@ write_connectivity_summary <- function(global_df, global_ttest_df, add("") add("**Frequency Bands:** ", band_str) add("") - add("**Metrics:** Magnitude-squared coherence (MSC) and absolute imaginary coherence (|iCoh|)") + metric_desc <- c( + coherence = "magnitude-squared coherence (MSC)", + imag_coherence = "absolute imaginary coherence (|iCoh|; Nolte 2004)", + pli = "phase lag index (PLI; Stam 2007)", + wpli = "weighted PLI (wPLI; Vinck 2011)", + dwpli = "debiased weighted PLI (dwPLI; Vinck 2011)", + dpli = "directed PLI (dPLI; Stam & van Straaten 2012)", + aec = "orthogonalized amplitude envelope correlation (AEC; Hipp 2012)", + partial_corr = "partial correlation (Marrelec 2006)") + add("**Metrics:** ", paste(ifelse(is.na(metric_desc[mets]), mets, metric_desc[mets]), collapse = "; ")) add("") add("**Timeseries:** Signed (phase-preserving) ROI source timeseries used for all computations.") add("") @@ -881,6 +882,10 @@ dir.create(tbl_dir, showWarnings = FALSE, recursive = TRUE) message("Loading data...") edges <- read_csv(file.path(data_dir, "roi_connectivity_edges.csv"), show_col_types = FALSE) message(" roi_connectivity_edges.csv: ", nrow(edges), " rows") +if (length(present_metrics(edges)) == 0) + stop("roi_connectivity_edges.csv has no known connectivity metric column (", + paste(KNOWN_METRICS, collapse = ", "), ")") +message(" Metrics present: ", paste(present_metrics(edges), collapse = ", ")) # --- Load config --- config <- read_yaml(config_path) @@ -961,8 +966,8 @@ if (!figures_only) { # =========================================================================== message("\nGenerating figures...") -# Connectivity matrices (per metric) -for (mc in c("coherence", "imag_coherence")) { +# Connectivity matrices (one set per metric present in the edge CSV) +for (mc in present_metrics(edges)) { plot_connectivity_matrices(edges, mc, group_colors, group_labels, group_order, fig_dir) } diff --git a/R/roi_cross_freq_edges_analysis.R b/R/roi_cross_freq_edges_analysis.R new file mode 100644 index 0000000..d477aee --- /dev/null +++ b/R/roi_cross_freq_edges_analysis.R @@ -0,0 +1,322 @@ +#!/usr/bin/env Rscript +# roi_cross_freq_edges_analysis.R — AAC / PPC (edge-level cross-frequency) statistics +# +# Called by Python (roi_cross_freq) after roi_pac_analysis.R: +# Rscript R/roi_cross_freq_edges_analysis.R --data-dir ... --config ... --output-dir ... +# +# Reads aac_edges.csv and/or ppc_edges.csv exported by Python (one row per +# subject x freq_pair x roi_x x roi_y; the roi_x slow band is paired with the +# roi_y fast band, so the matrix is ASYMMETRIC and roi_x == roi_y is the local, +# within-ROI coupling). The declared design:/hypotheses: are tested at three +# tiers, mirroring roi_directed: +# 1. Global — per-subject mean over all ROI x ROI cells, one cell per freq_pair +# (directed-edge adapter with a single synthetic "(global)" edge) +# 2. Edges — mass-univariate per ordered roi_x -> roi_y cell +# (directed-edge adapter; FDR across the edge family) +# 3. Region — emmeans-tabular over directed region pairs (roi_categories) +# The band axis is freq_pair (slow-fast), not the study bands, so the explicit +# freq_pair level vector is passed as `bands` (same pattern as roi_pac_analysis.R). +# +# Output tables (per metric m in {aac, ppc}): +# roi_cross_freq__global_hypotheses.csv +# roi_cross_freq__directed_edges_hypotheses.csv +# roi_cross_freq__region_hypotheses.csv (roi_categories only) +# Figure: roi_cross_freq__global_bar.png +# Report: appends an AAC/PPC section to ANALYSIS_SUMMARY.md (written by the PAC +# script), or creates it when the PAC script did not run. + +library(argparse) +library(yaml) +library(readr) +library(dplyr) +library(tidyr) +library(ggplot2) + +script_dir <- if (exists("script.dir")) { + script.dir +} else { + tryCatch({ + args_all <- commandArgs(trailingOnly = FALSE) + file_arg <- grep("^--file=", args_all, value = TRUE) + if (length(file_arg) > 0) { + dirname(normalizePath(sub("^--file=", "", file_arg))) + } else { + "R" + } + }, error = function(e) "R") +} + +source(file.path(script_dir, "stats_utils.R")) +source(file.path(script_dir, "hypothesis.R")) + +has_lme4 <- requireNamespace("lme4", quietly = TRUE) && + requireNamespace("lmerTest", quietly = TRUE) && + requireNamespace("emmeans", quietly = TRUE) + +EDGE_METRICS <- c("aac", "ppc") +METRIC_TITLE <- c(aac = "Amplitude-Amplitude Coupling (AAC)", + ppc = "n:m Phase-Phase Coupling (PPC)") + +theme_pub <- function(base_size = 14) { + theme_minimal(base_size = base_size) + + theme( + panel.grid.minor = element_blank(), + panel.grid.major = element_line(color = "grey92"), + strip.text = element_text(face = "bold", size = base_size), + legend.position = "bottom", + plot.title = element_text(face = "bold", size = base_size + 2) + ) +} + +# --- Helpers --------------------------------------------------------------- + +#' Per-subject mean over all ROI x ROI cells, per freq_pair, per DV. +compute_global_edges <- function(edges, dv_cols) { + edges %>% + group_by(subject, group, freq_pair) %>% + summarise(across(all_of(dv_cols), ~ mean(.x, na.rm = TRUE), .names = "mean_{.col}"), + n_cells = n(), .groups = "drop") +} + +#' Map roi_x -> roi_y cells onto DIRECTED region pairs (slow-region -> fast-region) +#' and average within subject x freq_pair x region_pair. +aggregate_cells_to_region_pairs <- function(edges, roi_categories, dv_cols) { + roi_to_region <- data.frame( + roi = unlist(roi_categories), + region = rep(names(roi_categories), lengths(roi_categories)), + stringsAsFactors = FALSE + ) + edges %>% + inner_join(roi_to_region, by = c("roi_x" = "roi")) %>% + rename(region_x = region) %>% + inner_join(roi_to_region, by = c("roi_y" = "roi")) %>% + rename(region_y = region) %>% + mutate(region_pair = paste(region_x, "->", region_y)) %>% + group_by(subject, group, freq_pair, region_pair) %>% + summarise(across(all_of(dv_cols), ~ mean(.x, na.rm = TRUE)), + n_cells = n(), .groups = "drop") +} + +plot_global_bar <- function(global_df, dv, metric, group_colors, group_labels, + group_order, output_dir) { + mean_col <- paste0("mean_", dv) + if (!mean_col %in% names(global_df)) return(invisible(NULL)) + plot_data <- global_df %>% + filter(group %in% group_order) %>% + mutate(value = .data[[mean_col]], + group_label = factor(group_labels[group], levels = group_labels[group_order])) + color_vals <- group_colors[group_order] + names(color_vals) <- group_labels[group_order] + + p <- ggplot(plot_data, aes(x = freq_pair, y = value, fill = group_label)) + + geom_boxplot(width = 0.6, alpha = 0.7, position = position_dodge(0.8), + outlier.shape = NA) + + geom_jitter(aes(color = group_label), + position = position_jitterdodge(dodge.width = 0.8, jitter.width = 0.1), + size = 1.5, alpha = 0.5, show.legend = FALSE) + + scale_fill_manual(values = color_vals, name = NULL) + + scale_color_manual(values = color_vals, name = NULL) + + labs(x = "Frequency pair (slow-fast)", + y = paste0("Global ", toupper(dv), " (mean of all ROI x ROI cells)"), + title = paste0("Global ", METRIC_TITLE[metric], " by Frequency Pair and Group")) + + theme_pub() + + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + + n_pairs <- length(unique(plot_data$freq_pair)) + fname <- paste0("roi_cross_freq_", metric, "_global_bar", + if (dv == metric) "" else paste0("_", sub(paste0("^", metric, "_"), "", dv)), + ".png") + ggsave(file.path(output_dir, fname), p, width = max(8, 2 * n_pairs), height = 5, dpi = 300) + message(" Saved: ", fname) +} + +# --- Argument parsing -------------------------------------------------------- + +parser <- ArgumentParser(description = "AAC / PPC edge-level cross-frequency statistics (R)") +parser$add_argument("--data-dir", required = TRUE, + help = "Directory containing aac_edges.csv / ppc_edges.csv") +parser$add_argument("--config", required = TRUE, help = "Path to study YAML config") +parser$add_argument("--output-dir", required = TRUE, + help = "Root output directory for this analysis") +parser$add_argument("--fig-dir", default = NULL, + help = "Directory for figures (default: output-dir/figures)") +parser$add_argument("--tbl-dir", default = NULL, + help = "Directory for tables (default: output-dir/tables)") +parser$add_argument("--figures-only", action = "store_true", default = FALSE, + help = "Skip statistics; regenerate figures from existing data/tables") +parser$add_argument("--no-figures", action = "store_true", default = FALSE, + help = "Skip all figure generation (stats/tables only)") +parser$add_argument("--roi-categories", default = NULL, + help = "Path to roi_categories.yaml (atlas ROI groupings)") +parser$add_argument("--hypothesis", default = NULL, + help = "Comma-separated declared hypothesis name(s) to run (default: all)") +parser$add_argument("--metric", default = NULL, + help = "Comma-separated subset of aac,ppc (default: every edge CSV present)") +args <- parser$parse_args() + +data_dir <- args$data_dir +output_dir <- args$output_dir +figures_only <- args$figures_only +no_figures <- args$no_figures +if (no_figures) ggsave <- function(...) invisible(NULL) + +fig_dir <- if (!is.null(args$fig_dir)) args$fig_dir else file.path(output_dir, "figures") +tbl_dir <- if (!is.null(args$tbl_dir)) args$tbl_dir else file.path(output_dir, "tables") +dir.create(fig_dir, showWarnings = FALSE, recursive = TRUE) +dir.create(tbl_dir, showWarnings = FALSE, recursive = TRUE) + +config <- read_yaml(args$config) +group_colors <- unlist(config$group_colors) +group_labels <- unlist(config$groups) +group_order <- config$group_order +message("Study: ", config$name) + +if (!is.null(args$roi_categories) && file.exists(args$roi_categories)) { + rc <- read_yaml(args$roi_categories) + if (length(rc) == 1 && identical(names(rc), "roi_categories")) rc <- rc[["roi_categories"]] + config$roi_categories <- rc + message("Loaded roi_categories from: ", args$roi_categories, + " (", length(config$roi_categories), " regions)") +} + +metrics <- if (!is.null(args$metric)) trimws(strsplit(args$metric, ",")[[1]]) else EDGE_METRICS +metrics <- intersect(EDGE_METRICS, metrics) +metrics <- metrics[file.exists(file.path(data_dir, paste0(metrics, "_edges.csv")))] +if (length(metrics) == 0) { + message("No aac_edges.csv / ppc_edges.csv in ", data_dir, " — nothing to do.") + quit(status = 0) +} + +spec <- tryCatch(parse_design_spec(config), error = function(e) NULL) +if (!figures_only && (is.null(spec) || length(spec$hypotheses) == 0)) + stop("No design:/hypotheses: declared in config — nothing to test.") + +report_lines <- character() +add <- function(...) report_lines <<- c(report_lines, paste0(...)) + +for (metric in metrics) { + message("\n=== ", METRIC_TITLE[metric], " ===") + edges <- read_csv(file.path(data_dir, paste0(metric, "_edges.csv")), show_col_types = FALSE) + message(" ", metric, "_edges.csv: ", nrow(edges), " rows") + dv_cols <- intersect(c(metric, paste0(metric, "_z")), names(edges)) + if (length(dv_cols) == 0) { + message(" no '", metric, "' column — skipping") + next + } + freq_pairs <- unique(as.character(edges$freq_pair)) + prefix <- paste0("roi_cross_freq_", metric) + global_df <- compute_global_edges(edges, dv_cols) + + hyp_global <- data.frame(); hyp_edges <- data.frame(); hyp_region <- data.frame() + if (!figures_only) { + # 1. Global: one synthetic edge, lm(mean_dv ~ group) per freq_pair. + message(" --- Global (hypothesis layer) ---") + global_edges <- global_df + global_edges$roi_x <- "(global)"; global_edges$roi_y <- "(global)" + hyp_global <- bind_rows(lapply(dv_cols, function(dv) { + h <- run_directed_edges(global_edges, names(spec$hypotheses), spec, + dv_col = paste0("mean_", dv), source_col = "roi_x", + target_col = "roi_y", band_col = "freq_pair", + bands = freq_pairs) + if (nrow(h) > 0) h$dv <- dv + h + })) + if (!is.null(args$hypothesis) && nrow(hyp_global) > 0) + hyp_global <- hyp_global[hyp_global$hypothesis %in% + trimws(strsplit(args$hypothesis, ",")[[1]]), , drop = FALSE] + if (nrow(hyp_global) > 0) { + write_csv(hyp_global, file.path(tbl_dir, paste0(prefix, "_global_hypotheses.csv"))) + message(" Saved: ", prefix, "_global_hypotheses.csv (", nrow(hyp_global), " rows)") + } + + # 2. Cells: mass-univariate per ordered roi_x -> roi_y. + message(" --- ROI x ROI cells (hypothesis layer, mass-univariate) ---") + hyp_edges <- write_module_directed_edges( + edges, config, tbl_dir, prefix = prefix, dv_cols = dv_cols, + source_col = "roi_x", target_col = "roi_y", band_col = "freq_pair", + bands = freq_pairs, hypothesis = args$hypothesis) + if (is.null(hyp_edges)) hyp_edges <- data.frame() + + # 3. Region pairs (emmeans-tabular, group * region_pair). + if (length(config$roi_categories) > 0 && has_lme4) { + message(" --- Directed region pairs (hypothesis layer) ---") + region_df <- aggregate_cells_to_region_pairs(edges, config$roi_categories, dv_cols) + message(" Aggregated to ", length(unique(region_df$region_pair)), " directed region pairs") + hyp_region <- write_module_hypotheses( + region_df, config, tbl_dir, prefix = paste0(prefix, "_region"), + dv_cols = dv_cols, spatial_col = "region_pair", band_col = "freq_pair", + bands = freq_pairs, hypothesis = args$hypothesis) + if (is.null(hyp_region)) hyp_region <- data.frame() + } else if (length(config$roi_categories) == 0) { + message(" No roi_categories in config -- skipping region-pair tier") + } else { + message(" lme4/lmerTest/emmeans not available -- skipping region-pair tier") + } + } + + # Figures + for (dv in dv_cols) + plot_global_bar(global_df, dv, metric, group_colors, group_labels, group_order, fig_dir) + + # Report section + if (!figures_only) { + add("## ", METRIC_TITLE[metric], " — edge-level cross-frequency") + add("") + add("Cells are ROI×ROI per slow–fast frequency pair (roi_x carries the slow band, ", + "roi_y the fast band; the matrix is asymmetric and roi_x = roi_y is the local ", + "within-ROI coupling). DV(s): ", paste(dv_cols, collapse = ", "), ".") + add("") + add("**Statistics (declarative hypothesis layer):** (1) global — per-subject mean over all ", + "cells, lm(mean ~ group) per frequency pair; (2) cells — mass-univariate contrast per ", + "ordered roi_x→roi_y cell, FDR across the cell family per frequency pair", + if (nrow(hyp_region) > 0) "; (3) directed region pairs — emmeans contrast over dv ~ group * region_pair + (1|subject)" else "", + ".") + add("") + gc <- if (nrow(hyp_global) > 0) hyp_global[hyp_global$kind == "contrast", , drop = FALSE] else data.frame() + if (nrow(gc) > 0) { + add("| Hypothesis | DV | Freq pair | estimate | t | q | Hedges' g | Sig |") + add("| --- | --- | --- | --- | --- | --- | --- | --- |") + for (i in seq_len(nrow(gc))) { + row <- gc[i, ] + add(sprintf("| %s | %s | %s | %.4f | %.2f | %.4f | %.2f | %s |", + row$label %||% row$hypothesis, row$dv, row$band, + ifelse(is.na(row$estimate), 0, row$estimate), + ifelse(is.na(row$stat), 0, row$stat), + ifelse(is.na(row$q_value), 1, row$q_value), + ifelse(is.na(row$effect_size), 0, row$effect_size), + if (isTRUE(row$significant)) "**Yes**" else "No")) + } + add("") + } else { + add("*No global contrast rows.*") + add("") + } + n_sig_cells <- if (nrow(hyp_edges) > 0) sum(hyp_edges$significant, na.rm = TRUE) else 0 + add("**Cells tested:** ", if (nrow(hyp_edges) > 0) length(unique(hyp_edges$spatial)) else 0, + " per frequency pair; **significant after FDR:** ", n_sig_cells, ".") + add("") + add("Tables: `", prefix, "_global_hypotheses.csv`, `", prefix, "_directed_edges_hypotheses.csv`", + if (nrow(hyp_region) > 0) paste0(", `", prefix, "_region_hypotheses.csv`") else "", ".") + add("") + } +} + +if (!figures_only && length(report_lines) > 0) { + summary_path <- file.path(output_dir, "ANALYSIS_SUMMARY.md") + if (file.exists(summary_path)) { + existing <- readLines(summary_path, warn = FALSE) + # Replace a previous AAC/PPC block (idempotent re-runs) before appending. + marker <- "" + cut <- which(existing == marker) + if (length(cut) > 0) existing <- existing[seq_len(cut[1] - 1)] + writeLines(c(existing, marker, "", report_lines), summary_path) + message(" Report section appended: ", summary_path) + } else { + writeLines(c(paste0("# ROI Cross-Frequency (AAC / PPC) — ", config$name), "", + paste0("**Generated:** ", format(Sys.time(), "%Y-%m-%d %H:%M")), "", + "", "", report_lines), summary_path) + message(" Report written: ", summary_path) + } +} + +message("\nDone. Output: ", output_dir) diff --git a/R/roi_evoked_analysis.R b/R/roi_evoked_analysis.R index 4cfc83a..9598bad 100644 --- a/R/roi_evoked_analysis.R +++ b/R/roi_evoked_analysis.R @@ -33,6 +33,7 @@ script_dir <- if (exists("script.dir")) { } source(file.path(script_dir, "stats_utils.R")) +source(file.path(script_dir, "hypothesis.R")) source(file.path(script_dir, "report.R")) # --- Argument parsing --- @@ -53,6 +54,8 @@ parser$add_argument("--no-figures", action = "store_true", default = FALSE, help = "Skip all figure generation (stats/tables only)") parser$add_argument("--roi-categories", default = NULL, help = "Path to roi_categories.yaml (atlas ROI groupings)") +parser$add_argument("--hypothesis", default = NULL, + help = "Run only the named hypothesis(es) (comma-separated) from the design spec; default = all") args <- parser$parse_args() data_dir <- args$data_dir @@ -89,6 +92,20 @@ group_colors <- unlist(config$group_colors) group_labels <- unlist(config$groups) group_order <- config$group_order +# Pairwise contrast list for the descriptive group*spatial LMM below. Derived from +# the design:/hypotheses: spec (config$contrasts is NULL for modern configs); +# falls back to a legacy contrasts: block if the spec has none. +contrasts <- tryCatch(contrasts_from_spec(parse_design_spec(config)), + error = function(e) config$contrasts) +if (is.null(contrasts)) contrasts <- list() +if (!is.null(args$hypothesis) && length(contrasts) > 0) { + want <- trimws(strsplit(args$hypothesis, ",")[[1]]) + contrasts <- Filter(function(ct) ct$name %in% want, contrasts) + message("Contrasts narrowed to --hypothesis set: ", + paste(vapply(contrasts, function(ct) ct$name, character(1)), collapse = ", ")) +} +message("Contrasts: ", length(contrasts)) + message("Study: ", config$name) message("Groups: ", paste(group_order, collapse = ", ")) @@ -108,7 +125,7 @@ if (!figures_only) { all_omnibus <- list() all_posthoc <- list() -for (contrast in config$contrasts) { +for (contrast in contrasts) { cname <- contrast$name ga <- contrast$group_a gb <- contrast$group_b @@ -208,7 +225,7 @@ if (nrow(omnibus_df) > 0) { # --- Post-hoc emmeans (gated on significant omnibus) --- message("\nRunning post-hoc emmeans...") -for (contrast in config$contrasts) { +for (contrast in contrasts) { cname <- contrast$name ga <- contrast$group_a gb <- contrast$group_b @@ -310,7 +327,7 @@ if (length(config$roi_categories) > 0) { all_omnibus_reg <- list() all_posthoc_reg <- list() - for (contrast in config$contrasts) { + for (contrast in contrasts) { cname <- contrast$name ga <- contrast$group_a gb <- contrast$group_b @@ -394,7 +411,7 @@ if (length(config$roi_categories) > 0) { } # Region post-hoc - for (contrast in config$contrasts) { + for (contrast in contrasts) { cname <- contrast$name ga <- contrast$group_a gb <- contrast$group_b @@ -500,11 +517,11 @@ global_posthoc_list <- list() for (mname in measure_names) { gp_data <- measures_df %>% filter(measure_name == mname, - group %in% unlist(lapply(config$contrasts, function(c) c(c$group_a, c$group_b)))) + group %in% unlist(lapply(contrasts, function(c) c(c$group_a, c$group_b)))) if (nrow(gp_data) == 0) next - gp <- run_posthoc_global(gp_data, config$contrasts, spatial_col = "roi", + gp <- run_posthoc_global(gp_data, contrasts, spatial_col = "roi", dv_col = "value", dv_label = mname) if (nrow(gp) > 0) global_posthoc_list[[mname]] <- gp } diff --git a/R/roi_pac_analysis.R b/R/roi_pac_analysis.R index 17e67e5..41fdfa8 100644 --- a/R/roi_pac_analysis.R +++ b/R/roi_pac_analysis.R @@ -33,7 +33,7 @@ script_dir <- if (exists("script.dir")) { }, error = function(e) "R") } -tryCatch(source(file.path(script_dir, "stats_utils.R")), error = function(e) NULL) +source(file.path(script_dir, "stats_utils.R")) source(file.path(script_dir, "hypothesis.R")) # --- Hypothesis-layer adapters -------------------------------------------- diff --git a/R/roi_transfer_entropy_analysis.R b/R/roi_transfer_entropy_analysis.R index f4ea79c..a61d054 100644 --- a/R/roi_transfer_entropy_analysis.R +++ b/R/roi_transfer_entropy_analysis.R @@ -1,9 +1,13 @@ #!/usr/bin/env Rscript -# roi_transfer_entropy_analysis.R — Transfer entropy statistics and report +# roi_transfer_entropy_analysis.R — roi_directed statistics and report (TE + DTF) +# +# Output tables/figures carry the module's canonical prefix `roi_directed_*`; +# this file keeps its legacy name because Python resolves it by that name. # # Called by Python: Rscript R/roi_transfer_entropy_analysis.R --data-dir ... --config ... --output-dir ... # -# Reads transfer_entropy_edges.csv exported by Python. +# Reads roi_transfer_entropy_edges.csv exported by Python. Directed DVs present +# in the CSV (te, net_te, dtf) are all tested; a DTF-only run works. # Three analysis tiers: # 1. Global TE: mean TE across all directed edges per subject x band, Welch t-test, BH FDR # 2. Directional: paired t-test on TE(X→Y) vs TE(Y→X) within groups (test for net directionality) @@ -37,7 +41,7 @@ script_dir <- if (exists("script.dir")) { } }, error = function(e) "R") } -tryCatch(source(file.path(script_dir, "stats_utils.R")), error = function(e) NULL) +source(file.path(script_dir, "stats_utils.R")) source(file.path(script_dir, "hypothesis.R")) # --- Hypothesis-layer adapters -------------------------------------------- @@ -86,15 +90,27 @@ theme_pub <- function(base_size = 14) { # 1. Global TE analysis: Welch t-tests # =========================================================================== +# Directed DVs the module can export; whichever are present in the edge CSV +# are tested (te + net_te from --metric te, dtf from --metric dtf). +DIRECTED_DVS <- c("te", "net_te", "dtf") + compute_global_te <- function(edges) { - edges %>% + present <- intersect(DIRECTED_DVS, names(edges)) + out <- edges %>% group_by(subject, group, band) %>% summarise( - mean_te = mean(te, na.rm = TRUE), - mean_abs_net_te = mean(abs(net_te), na.rm = TRUE), + across(any_of(setdiff(present, "net_te")), ~ mean(.x, na.rm = TRUE), + .names = "mean_{.col}"), n_edges = n(), .groups = "drop" ) + if ("net_te" %in% present) { + abs_net <- edges %>% + group_by(subject, group, band) %>% + summarise(mean_abs_net_te = mean(abs(net_te), na.rm = TRUE), .groups = "drop") + out <- left_join(out, abs_net, by = c("subject", "group", "band")) + } + out } run_global_ttests <- function(global_df, contrasts, bands) { @@ -253,14 +269,13 @@ aggregate_edges_to_region_pairs <- function(edges, roi_categories) { edges_mapped %>% group_by(subject, group, band, region_pair) %>% summarise( - te = mean(te, na.rm = TRUE), - net_te = mean(net_te, na.rm = TRUE), + across(any_of(DIRECTED_DVS), ~ mean(.x, na.rm = TRUE)), n_edges = n(), .groups = "drop" ) } -run_omnibus_lmm <- function(region_pair_df, contrasts, bands) { +run_omnibus_lmm <- function(region_pair_df, contrasts, bands, dv = "te") { if (!has_lme4) { message(" lme4/lmerTest not available -- skipping region-pair LMM") return(data.frame()) @@ -293,7 +308,8 @@ run_omnibus_lmm <- function(region_pair_df, contrasts, bands) { converged <- TRUE; singular <- FALSE tryCatch({ - fit <- lmer(te ~ group * region_pair + (1 | subject), data = bdata) + bdata$.dv <- bdata[[dv]] + fit <- lmer(.dv ~ group * region_pair + (1 | subject), data = bdata) singular <- isSingular(fit) aov <- anova(fit, type = 3) @@ -438,7 +454,12 @@ run_posthoc_emmeans <- function(region_pair_df, contrasts, bands, omnibus_df, ga # =========================================================================== plot_global_te_bar <- function(global_df, group_colors, group_labels, - group_order, output_dir) { + group_order, output_dir, dv = "te") { + mean_col <- paste0("mean_", dv) + if (!mean_col %in% names(global_df)) return(invisible(NULL)) + dv_label <- c(te = "Transfer Entropy", dtf = "DTF")[dv] + if (is.na(dv_label)) dv_label <- toupper(dv) + global_df$mean_te <- global_df[[mean_col]] plot_data <- global_df %>% filter(group %in% group_order) %>% mutate( @@ -465,15 +486,16 @@ plot_global_te_bar <- function(global_df, group_colors, group_labels, size = 1.5, alpha = 0.5, show.legend = FALSE) + scale_fill_manual(values = color_vals, name = NULL) + scale_color_manual(values = color_vals, name = NULL) + - labs(x = "Frequency Band", y = "Global Transfer Entropy (mean of all directed edges)", - title = "Global Transfer Entropy by Band and Group") + + labs(x = "Frequency Band", y = paste0("Global ", dv_label, " (mean of all directed edges)"), + title = paste0("Global ", dv_label, " by Band and Group")) + theme_pub() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) n_bands <- length(unique(summary_data$band)) - ggsave(file.path(output_dir, "roi_transfer_entropy_global_bar.png"), p, + fname <- if (dv == "te") "roi_directed_global_bar.png" else paste0("roi_directed_global_bar_", dv, ".png") + ggsave(file.path(output_dir, fname), p, width = max(8, 2 * n_bands), height = 5, dpi = 300) - message(" Saved: roi_transfer_entropy_global_bar.png") + message(" Saved: ", fname) } # =========================================================================== @@ -509,7 +531,11 @@ write_te_summary <- function(global_df, hyp_global, hyp_edges, collapse = ", " ) - add("**Analysis:** Binned Transfer Entropy (Schreiber 2000)") + dv_desc <- c(te = "binned transfer entropy (Schreiber 2000)", + net_te = "net transfer entropy TE(X\u2192Y) \u2212 TE(Y\u2192X)", + dtf = "directed transfer function (Kami\u0144ski & Blinowska 1991; ridge-MVAR)") + add("**Analysis:** Directed connectivity \u2014 ", + paste(dv_desc[intersect(names(dv_desc), dv_cols)], collapse = "; ")) add("") add("**Groups:** ", group_str) add("") @@ -517,10 +543,16 @@ write_te_summary <- function(global_df, hyp_global, hyp_edges, add("") add("**Frequency Bands:** ", band_str) add("") - add("**Metric:** TE(X\u2192Y) = H(Y_f, Y_p) + H(Y_p, X_p) \u2212 H(Y_p) \u2212 H(Y_f, Y_p, X_p)") - add("") - add("**Parameters:** lag=1 sample, 5 equal-probability bins (quantile-based)") - add("") + if ("te" %in% dv_cols) { + add("**TE metric:** TE(X\u2192Y) = H(Y_f, Y_p) + H(Y_p, X_p) \u2212 H(Y_p) \u2212 H(Y_f, Y_p, X_p)") + add("") + add("**TE parameters:** lag=1 sample, 5 equal-probability bins (quantile-based)") + add("") + } + if ("dtf" %in% dv_cols) { + add("**DTF:** normalized directed transfer function from a ridge-regularized MVAR fit, band-averaged") + add("") + } add("**Timeseries:** Signed (phase-preserving) ROI source timeseries, band-pass filtered (4th-order Butterworth)") add("") @@ -529,8 +561,9 @@ write_te_summary <- function(global_df, hyp_global, hyp_edges, add("**Statistics (declarative hypothesis layer):**") add("") - add("1. **Global:** between-group contrast on per-subject mean TE (across all ", - n_edges_str, " directed edges) per band \u2014 lm(mean_te ~ group), per-band FDR.") + add("1. **Global:** between-group contrast on the per-subject mean of each directed DV (", + paste(setdiff(dv_cols, "net_te"), collapse = ", "), "; across all ", + n_edges_str, " directed edges) per band \u2014 lm(mean_dv ~ group), per-band FDR.") add("") add("2. **Directed edges:** mass-univariate group contrast per ordered ", "source\u2192target edge (the directed-edge adapter; source\u2192target and ", @@ -549,13 +582,13 @@ write_te_summary <- function(global_df, hyp_global, hyp_edges, gc <- if (!is.null(hyp_global) && nrow(hyp_global) > 0) hyp_global[hyp_global$kind == "contrast", , drop = FALSE] else data.frame() if (nrow(gc) > 0) { - add("| Hypothesis | Band | estimate | t | df | q | Hedges' g | Sig |") - add("| --- | --- | --- | --- | --- | --- | --- | --- |") + add("| Hypothesis | DV | Band | estimate | t | df | q | Hedges' g | Sig |") + add("| --- | --- | --- | --- | --- | --- | --- | --- | --- |") for (i in seq_len(nrow(gc))) { row <- gc[i, ] sig_str <- if (isTRUE(row$significant)) "**Yes**" else "No" - add(sprintf("| %s | %s | %.5f | %.2f | %.1f | %.4f | %.2f | %s |", - row$label %||% row$hypothesis, row$band, + add(sprintf("| %s | %s | %s | %.5f | %.2f | %.1f | %.4f | %.2f | %s |", + row$label %||% row$hypothesis, row$dv %||% "te", row$band, ifelse(is.na(row$estimate), 0, row$estimate), ifelse(is.na(row$stat), 0, row$stat), ifelse(is.na(row$df), 0, row$df), @@ -776,8 +809,13 @@ if (!is.null(args$roi_categories) && file.exists(args$roi_categories)) { # =========================================================================== # 1. Global TE analysis (needed for figures — always compute) # =========================================================================== -message("\n=== Global Transfer Entropy Analysis ===") +message("\n=== Global Directed Connectivity Analysis ===") +dv_cols <- intersect(DIRECTED_DVS, names(edges)) +if (length(dv_cols) == 0) + stop("roi_transfer_entropy_edges.csv has none of the directed DV columns: ", + paste(DIRECTED_DVS, collapse = ", ")) +message("Directed DVs present: ", paste(dv_cols, collapse = ", ")) global_df <- compute_global_te(edges) message(" Global TE computed: ", nrow(global_df), " subject x band rows") @@ -793,17 +831,21 @@ if (!figures_only) { # hypothesis layer has a single cell — handled by the directed-edge adapter # with one synthetic (global) edge (lm(mean_te ~ group), no spatial term). # =========================================================================== - message("\n=== Global TE (hypothesis layer) ===") + message("\n=== Global directed (hypothesis layer) ===") global_edges <- global_df global_edges$source_roi <- "(global)"; global_edges$target_roi <- "(global)" - hyp_global <- run_directed_edges(global_edges, names(spec$hypotheses), spec, - dv_col = "mean_te", band_col = "band") + hyp_global <- bind_rows(lapply(setdiff(dv_cols, "net_te"), function(dv) { + h <- run_directed_edges(global_edges, names(spec$hypotheses), spec, + dv_col = paste0("mean_", dv), band_col = "band") + if (nrow(h) > 0) h$dv <- dv + h + })) if (!is.null(args$hypothesis) && nrow(hyp_global) > 0) hyp_global <- hyp_global[hyp_global$hypothesis %in% trimws(strsplit(args$hypothesis, ",")[[1]]), , drop = FALSE] if (nrow(hyp_global) > 0) { - write_csv(hyp_global, file.path(tbl_dir, "roi_transfer_entropy_global_hypotheses.csv")) - message(" Saved: roi_transfer_entropy_global_hypotheses.csv (", nrow(hyp_global), " rows)") + write_csv(hyp_global, file.path(tbl_dir, "roi_directed_global_hypotheses.csv")) + message(" Saved: roi_directed_global_hypotheses.csv (", nrow(hyp_global), " rows)") for (i in seq_len(nrow(hyp_global))) { row <- hyp_global[i, ] sig_str <- if (isTRUE(row$significant)) " ***" else "" @@ -820,8 +862,8 @@ if (!figures_only) { # =========================================================================== message("\n=== Directed-Edge TE (hypothesis layer, mass-univariate) ===") hyp_edges <- write_module_directed_edges( - edges, config, tbl_dir, prefix = "roi_transfer_entropy", - dv_cols = "te", source_col = "source_roi", target_col = "target_roi", + edges, config, tbl_dir, prefix = "roi_directed", + dv_cols = dv_cols, source_col = "source_roi", target_col = "target_roi", band_col = "band", hypothesis = args$hypothesis) if (is.null(hyp_edges)) hyp_edges <- data.frame() @@ -843,8 +885,8 @@ if (!figures_only) { message(" Aggregated to ", n_region_pairs, " directed region pairs") hyp_region <- write_module_hypotheses( - region_pair_df, config, tbl_dir, prefix = "roi_transfer_entropy_region", - dv_cols = "te", spatial_col = "region_pair", band_col = "band", + region_pair_df, config, tbl_dir, prefix = "roi_directed_region", + dv_cols = dv_cols, spatial_col = "region_pair", band_col = "band", hypothesis = args$hypothesis) if (is.null(hyp_region)) hyp_region <- data.frame() @@ -853,7 +895,11 @@ if (!figures_only) { # Diagnostic omnibus LMM (group x region_pair interaction-F) — NOT a # hypothesis; provides the interaction the marginal hypothesis layer lacks. - omnibus_df <- run_omnibus_lmm(region_pair_df, te_contrasts, config$bands) + omnibus_df <- bind_rows(lapply(setdiff(dv_cols, "net_te"), function(dv) { + o <- run_omnibus_lmm(region_pair_df, te_contrasts, config$bands, dv = dv) + if (nrow(o) > 0) o$dv <- dv + o + })) if (nrow(omnibus_df) > 0) { message("\n === Region-Pair Omnibus (diagnostic) ===") for (i in seq_len(nrow(omnibus_df))) { @@ -865,8 +911,8 @@ if (!figures_only) { row$group_F, row$group_q, grp_sig, row$interaction_F, row$interaction_q, int_sig)) } - write_csv(omnibus_df, file.path(tbl_dir, "roi_transfer_entropy_omnibus_lmm.csv")) - message(" Saved: roi_transfer_entropy_omnibus_lmm.csv (diagnostic)") + write_csv(omnibus_df, file.path(tbl_dir, "roi_directed_omnibus_lmm.csv")) + message(" Saved: roi_directed_omnibus_lmm.csv (diagnostic)") } if (nrow(posthoc_df) > 0) { sig_count <- sum(posthoc_df$significant, na.rm = TRUE) @@ -881,10 +927,10 @@ if (!figures_only) { message("Figures-only mode: loading existing hypothesis tables...") rd <- function(f) tryCatch(read_csv(file.path(tbl_dir, f), show_col_types = FALSE), error = function(e) data.frame()) - hyp_global <- rd("roi_transfer_entropy_global_hypotheses.csv") - hyp_edges <- rd("roi_transfer_entropy_directed_edges_hypotheses.csv") - hyp_region <- rd("roi_transfer_entropy_region_hypotheses.csv") - omnibus_df <- rd("roi_transfer_entropy_omnibus_lmm.csv") + hyp_global <- rd("roi_directed_global_hypotheses.csv") + hyp_edges <- rd("roi_directed_directed_edges_hypotheses.csv") + hyp_region <- rd("roi_directed_region_hypotheses.csv") + omnibus_df <- rd("roi_directed_omnibus_lmm.csv") posthoc_df <- .te_contrast_rows(hyp_region) if (nrow(posthoc_df) > 0) posthoc_df$region_pair <- posthoc_df$spatial } @@ -893,7 +939,8 @@ if (!figures_only) { # Figures # =========================================================================== message("\nGenerating figures...") -plot_global_te_bar(global_df, group_colors, group_labels, group_order, fig_dir) +for (dv in setdiff(dv_cols, "net_te")) + plot_global_te_bar(global_df, group_colors, group_labels, group_order, fig_dir, dv = dv) if (!figures_only) { # =========================================================================== diff --git a/R/stats_utils.R b/R/stats_utils.R index 1c25b39..04220a2 100644 --- a/R/stats_utils.R +++ b/R/stats_utils.R @@ -39,7 +39,7 @@ order_bands <- function(x, ref = NULL) { #' Reports Type III ANOVA F-tests for group, roi, and group:roi interaction. #' FDR (BH) correction applied across bands within each contrast. #' -#' @param band_df data.frame with columns: subject, group, roi, band, absolute, relative, dB +#' @param band_df data.frame with columns: subject, group, roi, band, absolute (dB/Hz), relative [, delta_ref] #' @param contrasts list of lists, each with name, group_a, group_b #' @param bands named list of c(fmin, fmax) — used only for ordering #' @param power_type character — column to use as DV: "relative" or "absolute" @@ -138,7 +138,7 @@ run_omnibus_lmm <- function(band_df, contrasts, bands, power_type = "relative") #' Run emmeans post-hoc contrasts per ROI for significant omnibus results #' -#' @param band_df data.frame with columns: subject, group, roi, band, absolute, relative, dB +#' @param band_df data.frame with columns: subject, group, roi, band, absolute (dB/Hz), relative [, delta_ref] #' @param contrasts list of lists, each with name, group_a, group_b #' @param bands named list of c(fmin, fmax) #' @param omnibus_df data.frame from run_omnibus_lmm() @@ -229,7 +229,7 @@ run_posthoc_emmeans <- function(band_df, contrasts, bands, omnibus_df, #' Aggregate ROI-level data to region-level means #' -#' @param band_df data.frame with columns: subject, group, roi, band, absolute, relative, dB +#' @param band_df data.frame with columns: subject, group, roi, band, absolute (dB/Hz), relative [, delta_ref] #' @param roi_categories named list of ROI name vectors #' @return data.frame with 'region' column replacing 'roi' aggregate_to_regions <- function(band_df, roi_categories) { @@ -239,12 +239,13 @@ aggregate_to_regions <- function(band_df, roi_categories) { stringsAsFactors = FALSE ) + # Average every power DV present (absolute / relative / delta_ref) so a + # config `dvs:` that includes delta_ref survives the region aggregate. band_df %>% inner_join(roi_to_region, by = "roi") %>% group_by(subject, group, region, band) %>% summarise( - absolute = mean(absolute, na.rm = TRUE), - relative = mean(relative, na.rm = TRUE), + across(any_of(c("absolute", "relative", "delta_ref")), ~ mean(.x, na.rm = TRUE)), .groups = "drop" ) } @@ -448,7 +449,7 @@ run_posthoc_emmeans_region <- function(band_df, contrasts, bands, roi_categories #' ROIs/electrodes are retained as replicate observations within regions. #' Model: dv ~ group * region + (1|subject) #' -#' @param band_df data.frame with columns: subject, group, roi, band, absolute, relative, dB +#' @param band_df data.frame with columns: subject, group, roi, band, absolute (dB/Hz), relative [, delta_ref] #' @param contrasts list of contrast definitions #' @param bands named list of frequency band limits #' @param roi_categories named list of ROI/electrode name vectors per region @@ -558,7 +559,7 @@ run_omnibus_lmm_region_nested <- function(band_df, contrasts, bands, roi_categor #' Same model as run_omnibus_lmm_region_nested() — individual ROIs/electrodes #' retained as replicates. emmeans(fit, pairwise ~ group | region). #' -#' @param band_df data.frame with columns: subject, group, roi, band, absolute, relative, dB +#' @param band_df data.frame with columns: subject, group, roi, band, absolute (dB/Hz), relative [, delta_ref] #' @param contrasts list of contrast definitions #' @param bands named list of frequency band limits #' @param roi_categories named list of ROI/electrode name vectors per region diff --git a/R/vertex_spatial_analysis.R b/R/vertex_spatial_analysis.R index 78c3d03..2c004ad 100644 --- a/R/vertex_spatial_analysis.R +++ b/R/vertex_spatial_analysis.R @@ -10,7 +10,6 @@ suppressPackageStartupMessages({ library(optparse) library(yaml) - library(nlme) library(lme4) library(lmerTest) library(dplyr) diff --git a/README.md b/README.md index 54e8b5e..1a40b77 100644 --- a/README.md +++ b/README.md @@ -55,24 +55,35 @@ You have run `source-localization` and have a derivatives tree of per-subject reconstructions. Three commands take you from there to results: ```bash -# 1. Scaffold a study config from the reconstruction directory. -# (--groups-from reuses the group mapping from the source-localization config.) +# 1. Scaffold a study config from the reconstruction directory. It is WRITTEN to +# /analysis/.yaml (status goes to stderr); pass `--output -` to +# print the YAML to stdout instead. --groups-from reuses the subject→group +# mapping from the source-localization config. source-analytics init /path/to/localization/rest_roi \ - --name "My Study" \ - --groups-from /path/to/localization/study_config.yaml > study.yaml + --name study \ + --groups-from /path/to/localization/study_config.yaml +# -> /path/to/localization/rest_roi/analysis/study.yaml -# 2. Edit study.yaml — declare groups, bands, the design/hypotheses, and which -# analyses to run under each paradigm (see "Study configuration"). +# 2. Edit that file — it already has groups, `design:`/`hypotheses:` (an omnibus +# plus every pairwise contrast), the canonical bands, and one `paradigms:` +# block (`resting`, with roi_psd / roi_aperiodic / roi_connectivity). Add +# paradigms and analyses as needed (see "Study configuration"). # 3. Sanity-check config + subject discovery before any long run. -source-analytics validate --study study.yaml +source-analytics validate --study rest_roi/analysis/study.yaml # 4. Run an analysis. --paradigm picks the block under `paradigms:` in the config. -source-analytics run --study study.yaml --paradigm resting --analysis roi_psd +source-analytics run --study rest_roi/analysis/study.yaml --paradigm resting --analysis roi_psd ``` -Each analysis writes a self-contained output directory with a `data/`, `tables/`, -`figures/`, and an `ANALYSIS_SUMMARY.md` (see [Output structure](#output-structure)). +`init` discovers subjects in either layout the toolkit reads: BIDS-style +`derivatives/sub-*/` (groups come from `--groups-from`, otherwise `UNKNOWN` until +you edit them) or `derivatives///` (the folder name is the group). + +Each run writes per-subject data + `ANALYSIS_SUMMARY.md` under `paths.analytics` +and the published `tables/` + `figures/` under `paths.results` (see +[Output structure](#output-structure)). **Figures are not produced by a default +run** — add `--steps …,figures` or use `source-analytics figure`. `source-analytics list` shows every analysis you can run. --- @@ -82,11 +93,25 @@ Each analysis writes a self-contained output directory with a `data/`, `tables/` ### Python (3.10+) ```bash -pip install -e . # or: uv pip install -e . +pip install -e ".[all]" # or: uv pip install -e ".[all]" ``` -Core dependencies: numpy, scipy, pandas, pyyaml, scikit-learn, networkx, -matplotlib, mne. +Core dependencies (always installed): numpy, scipy, pandas, pyyaml, specparam, +joblib, matplotlib. The rest are **extras** — the package imports and the CLI +works without them, and a module that needs one fails with an ImportError naming +the extra to install: + +| Extra | Pulls in | Needed by | +|---|---|---| +| `mne` | mne | `roi_evoked`, `vertex_evoked`, `electrode_evoked` (Morlet TFR) | +| `mvpa` | scikit-learn | `vertex_signature`, `electrode_signature` | +| `network` | networkx | `roi_graph`, `vertex_graph`, `*_nbs`, `*_network` | +| `atlas` | nibabel | atlas readers | +| `all` | all of the above + dev tools | a full study | + +Wheels/sdists ship the R scripts under `/share/source-analytics/R`; an +editable checkout uses `R/` beside `src/`. `SOURCE_ANALYTICS_R_DIR` overrides +the lookup. > **uv users:** run the CLI with `uv run --no-sync source-analytics …`. Plain > `uv run` can trip on the lockfile; `--no-sync` avoids the re-resolve. @@ -111,8 +136,10 @@ uv pip install --python .venv --no-deps "source-analytics @ git+https://github.c ⚠ **`v0.4.0` is a scientific pin, not merely an old version.** It hardcodes `freq_range=(2, 50)` for aperiodic fitting; later releases resolve the window -dynamically and adopt 14–45 Hz, which changes every aperiodic number. Do not -"upgrade" it to reproduce work that cites it. +dynamically and default to **12–45 Hz** (`spectral.aperiodic.DEFAULT_FREQ_RANGE`; +see [`docs/methods/APERIODIC_FIT_WINDOW.md`](docs/methods/APERIODIC_FIT_WINDOW.md)), +which changes every aperiodic number. Do not "upgrade" it to reproduce work that +cites it. ### R @@ -120,12 +147,16 @@ Statistics and most figures are R. Install once: ```r install.packages(c( - "ggplot2", "dplyr", "tidyr", "readr", "stringr", "forcats", + "ggplot2", "dplyr", "tidyr", "readr", "forcats", "ggsignif", "lme4", "lmerTest", "effectsize", "emmeans", - "yaml", "argparse", "patchwork", "scales" + "yaml", "argparse", "optparse", "patchwork", "scales" )) ``` +(`ggsignif` draws the significance brackets in the PSD/aperiodic/evoked figures; +`optparse` is used by the vertex report scripts, `argparse` by the ROI/electrode +scripts.) + --- ## Core concepts @@ -201,11 +232,12 @@ What it looks for depends on the level: | `roi_timeseries_magnitude.set` | EEGLAB | same data + metadata (sfreq) | **Vertex-level** (`vertex_cluster`, `vertex_connectivity`, `vertex_cross_freq`, -`vertex_directed`, `vertex_specparam`, `vertex_mvpa`, `vertex_evoked`): +`vertex_directed`, `vertex_specparam`, `vertex_signature`, `vertex_evoked`): | File | Format | Contents | |------|--------|----------| -| `step5_stc_signed.pkl` | pickle | MNE `SourceEstimate` `(n_vertices, n_times)`, signed (falls back to `step5_stc_magnitude.pkl`, legacy `step5_stc.pkl`) | +| `step5_stc_signed.pkl` | pickle | MNE `SourceEstimate` `(n_vertices, n_times)`, signed. **No fallback**: phase-based modules refuse to read a magnitude file as signed | +| `step5_stc_magnitude.pkl` | pickle | rectified variant; only `magnitude=True` readers use it, and only they fall back to the legacy unsuffixed `step5_stc.pkl` (which is magnitude-only) | | `step3_source_coords_mm.npy` | NumPy | source coordinates `(n_vertices, 3)` in mm | **Electrode-level** (`electrode_psd`, `electrode_aperiodic`, @@ -218,15 +250,28 @@ What it looks for depends on the level: These need a `subject_roster.csv` (`subject_id, group, eeg_filename, eeg_dir`) set via `electrode.subject_roster` in the config. -Expected discovery layout (group folders → subject folders → `data/`): +Two discovery layouts are supported. **Grouped** (the folder name is the group): + +``` +data_dir/ + Group_A/Subject_001//… + Group_A/Subject_002//… + Group_B/Subject_003//… +``` + +**Flat** (BIDS-style `sub-*`, what source-localization writes; groups come from a +`subjects:` map on the paradigm — `init` fills it from `--groups-from`): ``` data_dir/ - Group_A/Subject_001/data/… - Group_A/Subject_002/data/… - Group_B/Subject_003/data/… + sub-001//… + sub-002//… ``` +`data_subdir` defaults to `pipeline/data`. Only subjects whose group is +referenced by a declared hypothesis/contrast are analysed +(`StudyConfig.referenced_groups()`); a group nobody tests is silently skipped. + --- ## Study configuration @@ -272,8 +317,13 @@ bands: # name → [fmin, fmax] Hz High Gamma: [65, 80] circos_metrics: [imag_coherence, dwpli, pli, aec, coherence] # gallery circos chords + # (read by source-lightbox only; source-analytics + # passes it through untouched) + +jobs: -1 # default worker count for --jobs (-1/0 = all but one core) # ── Random epoch sampling (global default; per-analysis override below) ── +# Code defaults when the block is absent: enabled: false, n_bootstrap: 1. epoch_sampling: enabled: true epoch_duration_sec: 2.0 @@ -294,13 +344,17 @@ paradigms: roi_psd: {} roi_aperiodic: {} roi_connectivity: + metrics: [imag_coherence, dwpli, pli, aec, coherence] # subset of the ROI metric set epoch_sampling: {n_bootstrap: 0} # per-analysis override roi_graph: {connectivity_metrics: [imag_coherence, dwpli, pli, aec, coherence]} roi_nbs: {nbs_threshold: 2.5, nbs_permutations: 5000} roi_cross_freq: {} roi_directed: {} electrode_psd: {} - electrode_comparison: {} + electrode_comparison: {} # needs electrode_psd AND roi_psd + electrode_connectivity: {} + fcd_comparison: {} # needs electrode_connectivity (here) + # AND vertex_connectivity (vertex paradigm) vertex: data_dir: ./localization/rest_shell/derivatives @@ -324,7 +378,9 @@ paradigms: | `hypotheses[]` `{name, kind, weights/groups/predictor}` | hypothesis layer | the declarative tests, run by name via `--hypothesis` | | `hypotheses[]` `{label, role}` | figures, gallery | readable labels + grouping tag (no gating) | | `bands` | all spectral/connectivity | frequency bands analysed | -| `epoch_sampling` | spectral/connectivity | random-epoch resampling (`n_bootstrap: 0` = full timeseries) | +| `epoch_sampling` | spectral/connectivity, all levels | random-epoch resampling (`n_bootstrap: 0` = full timeseries). Precedence: global → `vertex.epoch_sampling` → per-analysis block | +| `jobs` | `run --jobs` default | worker count when `--jobs` is not given | +| `.{include_analyses, include_hypotheses, bands, rois}` | `run --profile` | a narrowed study written to its own tree (see below) | | `paths.{analytics, results}` | I/O + gallery | working vs published output trees | | `paradigms.

.data_dir` / `data_subdir` | discovery | where subject reconstructions live | | `paradigms.

.analyses.` | that analysis | enables it + sets its parameters | @@ -364,25 +420,39 @@ source-analytics run --study study.yaml --paradigm resting --analysis roi_psd [o | Flag | Meaning | |---|---| | `--study PATH` | study YAML (required) | -| `--paradigm NAME` | paradigm block under `paradigms:` (required for multi-paradigm configs) | -| `--analysis NAME` | analysis to run (see [catalog](#analysis-catalog--what-exists)) | +| `--paradigm NAME` | paradigm block under `paradigms:`. Omit it (and `--analysis`) on a multi-paradigm config to run **every** listed analysis of every paradigm; `--analysis` without `--paradigm` is an error | +| `--analysis NAME` | analysis to run (see [catalog](#analysis-catalog--what-exists)). Omit it with `--paradigm` to run everything listed for that paradigm | | `--steps a,b,…` | lifecycle steps to run. Valid: `setup, process, aggregate, statistics, figures, summary` | +| `--jobs N`, `-j N` | worker processes for the per-subject `process` step. `0`/`-1` = all but one core. Explicit `N` wins over the YAML `jobs:`; omitted = YAML value, else serial. Used by the vertex modules, `roi_connectivity`, `electrode_connectivity`; results are identical to serial | +| `--profile NAME` | run under the top-level `:` profile block (narrowed bands / ROIs / hypotheses / analyses) and write to a separate tree, `analytics//…` + `results//…`. Narrowing ROIs changes the FDR family, so profile q-values are not comparable to the default run's | | `--metric m,…` | restrict a module's metrics (shorthand for `--select metric=…`) | | `--band b,…` | restrict bands, case/format-insensitive (shorthand for `--select band=…`) | | `--hypothesis n,…` | test only these declared hypotheses (shorthand for `--select hypothesis=…`) | | `--select DIM=v,…` | generic sub-output selection, repeatable (see `list` for a module's dims) | -| `--force` | overwrite the output directory if it exists | -| `--strict-output` | error if the output directory exists (unless `--force`) | +| `--force` | remove the analysis's previous output first: its published `tables/` + `figures/` always, and its working dir (`data/` + summary) when the `process` step runs. Also overrides `--strict-output` | +| `--strict-output` | error if the working dir already holds output (unless `--force`) | -**Lifecycle steps.** A run is `setup → process → aggregate → statistics → figures → -summary`. `--steps` re-runs a subset against existing on-disk data — e.g. recompute -only the statistics and report after a config change, without reprocessing subjects: +**Lifecycle steps.** The full lifecycle is `setup → process → aggregate → +statistics → figures → summary`. **A default `run` executes everything except +`figures`** (`DEFAULT_RUN_STEPS` in `analyses/base.py`), so tables and the summary +appear but no images do. Render figures with an explicit step list, or with +`source-analytics figure` for the on-demand summary figures: ```bash +# full run including figures source-analytics run --study study.yaml --paradigm resting --analysis roi_psd \ - --steps statistics,summary,figures + --steps setup,process,aggregate,statistics,figures,summary +# recompute only statistics + figures + report from persisted data/ (no reprocessing) +source-analytics run --study study.yaml --paradigm resting --analysis roi_psd \ + --steps statistics,figures,summary ``` +`--steps` re-runs a subset against the on-disk `data/` of an earlier run; the +`figures` step clears the module's figure dir before regenerating so stale images +never linger. Deprecated analysis names (`psd`, `pac`, `vertex_mvpa`, …) still +resolve, and their output always lands under the **canonical** name (`roi_psd/`, +`roi_cross_freq/`, `vertex_signature/`). + ### `validate`, `list`, `figure`, `init` ```bash @@ -391,9 +461,14 @@ source-analytics list [--study study.yaml] # paradigm-aware when --stud source-analytics figure --study study.yaml --paradigm resting --analysis roi_psd --list source-analytics figure --study study.yaml --paradigm resting --analysis roi_psd \ --type effect_heatmap [--contrast disease_effect --band low_gamma] -source-analytics init /path/to/reconstruction_dir --name "Study" --groups-from sl_config.yaml +source-analytics init /path/to/reconstruction_dir --name study --groups-from sl_config.yaml \ + [--paradigm resting] [--analyses roi_psd,roi_aperiodic] [--output PATH | -] ``` +`init` writes `/analysis/.yaml` (or stdout with +`--output -`) and parses it back to prove the scaffold loads. `list` groups the +catalog by paradigm category and level, tagging each module's `--select` dims. + --- ## Analysis catalog — what exists @@ -408,17 +483,17 @@ directed families is tracked, equation-checked, in | Analysis | Level | Computes | Reference | |---|---|---|---| -| `roi_psd`, `electrode_psd` | ROI, elec | band power (Welch PSD: absolute/relative/dB) | Welch 1967 | -| `roi_aperiodic`, `electrode_aperiodic`, `vertex_specparam` | ROI, elec, vtx | 1/f aperiodic (offset, exponent) + oscillatory peaks | Donoghue 2020 (specparam) | -| `vertex_cluster` | vtx | per-vertex band power / fALFF / slope / peak-α, cluster-corrected maps | Maris & Oostenveld 2007 | -| `vertex_mvpa` | vtx | whole-brain pattern decoding (linear SVM, LOOCV, permutation) | — (linear SVM) | +| `roi_psd`, `electrode_psd` | ROI, elec | band power (Welch PSD). CSV columns: `absolute` = mean power density in dB/Hz, `relative` = fraction of total; optional `delta_ref` under a profile. There is no separate `dB` column | Welch 1967 | +| `roi_aperiodic`, `electrode_aperiodic`, `vertex_specparam` | ROI, elec, vtx | 1/f aperiodic (offset, exponent) + oscillatory peaks; default fit window **12–45 Hz** | Donoghue 2020 (specparam) | +| `vertex_cluster` | vtx | per-vertex band power (same dB/Hz `absolute` as `roi_psd`) / fALFF / slope / peak-α, cluster-corrected maps; honours `epoch_sampling` | Maris & Oostenveld 2007 | +| `vertex_signature` (alias `vertex_mvpa`) | vtx | whole-brain neural signature: multi-model decoding (PCA-reduced, back-projected), permutation p | — | | `vertex_spatial` *(RETIRED)* | vtx | was: spatial-covariance GLS robustness check | — | > `vertex_spatial` is **retired** (it produced a spatial-covariance robustness > table, never a manuscript result, and did not survive the design-spec migration). > Spatially-resolved vertex inference is `vertex_cluster` (glass-brain clusters) + -> `vertex_nbs` (network-based statistic). The module exits cleanly with empty -> frames + a note. +> `vertex_nbs` (network-based statistic). The module processes no subjects and +> calls no R: it writes empty result tables + a retirement note and exits. ### Connectivity (same-frequency functional connectivity) @@ -430,41 +505,54 @@ directed families is tracked, equation-checked, in | `electrode_connectivity` | elec | FC-six all-pairs + per-channel FCD — the **source-vs-sensor comparator** | as above | > `roi_network` / `vertex_network` are **combined aliases** that run graph + NBS -> together; the split modules (`*_graph`, `*_nbs`) are preferred for the gallery. -> `dpli` is directed and is auto-excluded from the undirected graph/NBS layer. +> together and write a Python summary (there is no R report for them); the split +> modules (`*_graph`, `*_nbs`) are preferred for the gallery. `dpli` is directed +> and is auto-excluded from the undirected graph/NBS layer. `roi_connectivity` +> takes an optional `metrics:` list in its config block (like +> `vertex_connectivity`); the R report adapts to whichever metric columns the +> edge CSV carries, so `--metric aec` alone is fine. ### Cross-frequency | Analysis | Level | Computes | Reference | |---|---|---|---| -| `roi_cross_freq`, `vertex_cross_freq` | ROI, vtx | PAC (Modulation Index, surrogate-z); cross-frequency AAC; n:m PPC | Tort 2010; Bruns 2000 / Masimore 2004; Tass 1998 / Palva 2005 | +| `roi_cross_freq`, `vertex_cross_freq` | ROI, vtx | PAC (Modulation Index, surrogate-z); cross-frequency AAC; n:m PPC. ROI: PAC hypotheses via `roi_pac_analysis.R`; AAC/PPC via `roi_cross_freq_edges_analysis.R`, three tiers each: `roi_cross_freq_{aac,ppc}_{global,directed_edges,region}_hypotheses.csv` (PPC has DVs `ppc` + `ppc_z`) | Tort 2010; Bruns 2000 / Masimore 2004; Tass 1998 / Palva 2005 | ### Directed | Analysis | Level | Computes | Reference | |---|---|---|---| -| `roi_directed` | ROI | transfer entropy (`te`, `net_te`); DTF (`dtf`, ridge-MVAR) | Schreiber 2000; Kamiński & Blinowska 1991 | -| `vertex_directed` | vtx | DTF outflow / inflow / netflow (ridge-MVAR), cluster-corrected | Kamiński & Blinowska 1991 | +| `roi_directed` | ROI | transfer entropy (`te`, `net_te`); DTF (`dtf`, ridge-MVAR). `--metric te,dtf`; hypothesis tables `roi_directed_{global,directed_edges,region}_hypotheses.csv` carry a `dv` column covering every exported DV (`te`, `net_te`, `dtf`) | Schreiber 2000; Kamiński & Blinowska 1991 | +| `vertex_directed` | vtx | DTF outflow / inflow / netflow (ridge-MVAR), cluster-corrected. Filter with `--select measure=outflow,…` (not `--metric`) | Kamiński & Blinowska 1991 | > Source ROIs/vertices are strongly collinear (mean inter-node |corr| ≈ 0.64), so > DTF uses a **ridge-regularized** MVAR — plain LS-MVAR is non-stationary; the > module warns if a fit is unstable. -### Sensor-level (validation) +### Source vs sensor (validation) -| Analysis | Level | Computes | -|---|---|---| -| `electrode_comparison` *(suppl. of `electrode_psd`)* | elec | source-vs-electrode effect-size validation | +| Analysis | Level | Reads | Computes | +|---|---|---|---| +| `electrode_comparison` *(suppl.)* | elec | `electrode_psd` **and** `roi_psd` (same paradigm) | source-vs-electrode band-power concordance + effect-size validation | +| `fcd_comparison` *(suppl.)* | elec | `electrode_connectivity` **and** `vertex_connectivity` — normally in *different* paradigms; sibling paradigm dirs are searched, or set `fcd_comparison.{sensor_dir,source_dir}` | source-vs-sensor FCD comparison (mean + spatial CV) per band × metric | +| `electrode_signature` *(suppl. of `electrode_psd`)* | elec | `electrode_psd` | sensor-level neural signature (decoding on electrode band power) — the source-vs-sensor counterpart of `vertex_signature` | + +`ANALYSIS_METADATA` records these as `supplements` (the primary the gallery nests +them under) plus `requires` (every upstream module, for run ordering). ### Evoked (trial-based paradigms only) | Analysis | Level | Computes | |---|---|---| -| `roi_evoked`, `vertex_evoked`, `electrode_evoked` | ROI, vtx, elec | ITC, ERSP, single-trial power | +| `roi_evoked`, `vertex_evoked`, `electrode_evoked` | ROI, vtx, elec | ITC (raw + debiased), ERSP, single-trial power, induced power, ERP amplitude/latency. ROI/electrode: descriptive `group × unit` LMM in R plus declared hypotheses (`_evoked_hypotheses.py`, measure as facet); vertex: cluster permutation plus declared hypotheses via the permutation adapter (`vertex_evoked_hypotheses.csv`, band = measure name, dv = measure type) | **Renames (2026-06).** `roi_pac` → `roi_cross_freq` (now also AAC + PPC); `roi_transfer_entropy` → `roi_directed`. Old names still work as deprecated aliases -(`psd`/`aperiodic`/`pac`/`mvpa`/`wholebrain`/… also map to the canonical names). +(`psd`/`aperiodic`/`pac`/`roi_pac`/`mvpa`/`vertex_mvpa`/`wholebrain`/`spatial_lmm`/ +`specparam_vertex`/`transfer_entropy`/`roi_transfer_entropy`/`evoked`/`electrode` +map to the canonical names; output always lands under the canonical directory). +Two R scripts keep their legacy filenames on purpose: `roi_pac_analysis.R` (for +`roi_cross_freq`) and `roi_transfer_entropy_analysis.R` (for `roi_directed`). --- @@ -518,10 +606,16 @@ source-analytics run --study study.yaml --paradigm resting \ Output is **additive**: a `_hypotheses.csv` written alongside the module's other tables, with one tidy row per band × spatial cell (estimate, SE, CI, stat, p, -q, effect size, `fdr_family`, plus legacy-named aliases for existing figure -consumers). Modules with multiple spatial tiers emit one table per tier — e.g. -`roi_directed` writes `…_global_hypotheses.csv`, `…_directed_edges_hypotheses.csv`, -and `…_region_hypotheses.csv`. +q, effect size, `fdr_family`). Modules with multiple spatial tiers emit one table +per tier — `roi_directed` writes `roi_directed_global_hypotheses.csv`, +`roi_directed_directed_edges_hypotheses.csv`, and `roi_directed_region_hypotheses.csv`. + +Two limits worth knowing: `kind: regression` is accepted by the config and the +emmeans (R) adapter, but the Python permutation (map) and edge (NBS) adapters +return **no rows** for it yet — a continuous predictor is not wired into those +paths. And the R evoked scripts and every emmeans module derive their contrast +list from `design:`/`hypotheses:`; a legacy `contrasts:` block is lifted into the +same spec, so both config styles work. --- @@ -549,23 +643,30 @@ Selectable dimensions vary by module (`metric` / `band` / `hypothesis` / ## Output structure -Each analysis writes a self-contained directory under `paths.analytics` -(working tree). The published `tables/` + `figures/` are mirrored to -`paths.results`, which `source-lightbox` reads. +Two trees, both segmented by paradigm (and by profile when `--profile` is used). +The **working tree** under `paths.analytics` holds the per-subject data and the +narrative; the **published tree** under `paths.results` holds the tables and +figures that `source-lightbox` reads. `tables/` and `figures/` do **not** live +inside the analysis working directory. ``` -// +/[/]// ANALYSIS_SUMMARY.md # methods + results narrative (markdown) data/ _*.csv # the computed per-subject measures (the inputs to stats) study_config.yaml # the resolved config snapshot used for this run - tables/ + +/[/]tables/// _hypotheses.csv # the hypothesis-layer result (one row per band×cell) … # any module-specific diagnostic tables - figures/ - *.png # ggplot2 / glass-brain / matplotlib figures +/[/]figures/// + *.png # ggplot2 / glass-brain / matplotlib figures (figures step only) ``` +`config.for_paradigm_analysis()` sets the working dir to `analytics/`; +`BaseAnalysis.tbl_dir` / `fig_dir` resolve the published dirs. A legacy +single-paradigm config (no `paradigms:`) omits the paradigm segment. + The `_hypotheses.csv` is the canonical statistical contract across all emmeans/permutation modules; figure and gallery consumers read it (plus legacy column aliases during the migration). @@ -579,7 +680,8 @@ column aliases during the migration). ## Running a full study, in order A study run is just the analyses invoked in dependency order (primaries before -their supplements). The canonical recipes live in `scripts/`; the essential order: +their supplements). [`scripts/run_study.sh study.yaml [run flags]`](scripts/run_study.sh) +is the canonical recipe; the essential order: ```bash SA="source-analytics run --study study.yaml --paradigm" @@ -594,8 +696,9 @@ $SA resting --analysis roi_cross_freq # PAC + AAC + PPC (--metric to p $SA resting --analysis roi_directed # transfer entropy + DTF (--metric te|dtf) $SA resting --analysis electrode_psd # PRIMARY $SA resting --analysis electrode_aperiodic -$SA resting --analysis electrode_comparison # ↳ after electrode_psd +$SA resting --analysis electrode_comparison # ↳ after electrode_psd AND roi_psd $SA resting --analysis electrode_connectivity # sensor FC comparator +$SA resting --analysis electrode_signature # ↳ after electrode_psd # Vertex paradigm — whole-brain $SA vertex --analysis vertex_connectivity # PRIMARY (slow; computes matrices) @@ -603,10 +706,13 @@ $SA vertex --analysis vertex_graph # ↳ after vertex_connectivity $SA vertex --analysis vertex_nbs # ↳ after vertex_connectivity $SA vertex --analysis vertex_cluster $SA vertex --analysis vertex_specparam -$SA vertex --analysis vertex_mvpa +$SA vertex --analysis vertex_signature $SA vertex --analysis vertex_cross_freq # local PAC + AAC + PPC $SA vertex --analysis vertex_directed # vertex DTF (outflow/inflow/netflow) +# Source-vs-sensor FCD: reads electrode_connectivity (resting) + vertex_connectivity (vertex) +$SA resting --analysis fcd_comparison + # Evoked paradigm (trial-based data only) $SA evoked --analysis roi_evoked $SA evoked --analysis vertex_evoked diff --git a/docs/methods/DESIGN_SPEC.md b/docs/methods/DESIGN_SPEC.md index 9073415..ee9555d 100644 --- a/docs/methods/DESIGN_SPEC.md +++ b/docs/methods/DESIGN_SPEC.md @@ -58,7 +58,7 @@ study.yaml: design: + hypotheses: ▼ hypothesis layer ───────────────────────────────────────────────────────────── R/hypothesis.R load_design_spec() + run_hypothesis() + EMMEANS adapter - (sourced by roi_psd_analysis.R, electrode_psd_analysis.R, …) + (sourced by roi_psd_analysis.R, electrode_analysis.R, …) src/source_analytics/ load_design_spec() + run_hypothesis() + PERMUTATION adapter hypothesis/ (built on stats/cluster_permutation.py, stats/graph_metrics.py; __init__.py imported by vertex_network_analysis.py, …) diff --git a/docs/methods/HYPOTHESIS.md b/docs/methods/HYPOTHESIS.md index c3069ff..fcff561 100644 --- a/docs/methods/HYPOTHESIS.md +++ b/docs/methods/HYPOTHESIS.md @@ -192,6 +192,14 @@ This fits modules with a modest spatial cardinality (≤ ~32 ROIs / ~30 channels and connectivity-**matrix** subnetworks (roi_nbs / vertex_nbs) belong to the edge/NBS adapter — a `group × edge` LMM is infeasible and wrong for either (~496 edges → ~2500 params). +**Long-DV modules take the Python tabular adapter instead.** When a module exports one `value` +column plus a label column saying what it holds (the evoked pair: `value` + `measure_name`), +there is no `dv_cols` vector to pass. Hand the label column to `write_module_hypotheses_tabular` +as a **facet** — `facet_cols=("measure_name",)` — which runs an independent FDR family per facet +across the band × spatial grid. Facet ≡ DV: R gives each `dv_col` its own `run_hypothesis()` call +and so its own family, and a facet is that same family boundary expressed in a long table. +`analyses/_evoked_hypotheses.py` is the live example, shared by both evoked modules. + ## 9. Status - ✅ **emmeans adapter** built + verified (bit-exact vs legacy on real data); wired into @@ -209,9 +217,19 @@ and connectivity-**matrix** subnetworks (roi_nbs / vertex_nbs) belong to the edg `vertex_nbs`, and the combined `roi_network`/`vertex_network` aliases** via `write_module_hypotheses_edge()` (additive `_hypotheses.csv`). Verified bit-exact vs the legacy NBS on real FORGE data (Low Gamma / imag_coherence / KO_VEH vs WT_VEH). -- ⏳ **deferred:** `roi_evoked` / `electrode_evoked` (long-format DV; no data in the resting - study). **specials:** `vertex_mvpa` (decoding), `vertex_spatial` (GLS), `electrode_comparison` - (agreement — may not take hypotheses). +- ✅ **evoked (tabular adapter, measure as facet):** wired into `roi_evoked` (spatial `roi`) and + `electrode_evoked` (spatial `channel`) via the shared `analyses/_evoked_hypotheses.py`. These + two export a long DV — one `value` column faceted by `measure_name`, not one column per + measure — so the measure is passed as a **facet**, giving an independent FDR family per measure + across the spatial grid. That is the only defensible family: the measures are on incomparable + scales (ITC 0-1, ERSP dB, ERP amplitude in signal units, latency in seconds). There is **no band + axis** — each measure definition already fixes its own band and time window — so the band + coordinate is null. Additive: the descriptive `group * roi` LMM in the R scripts is unchanged + and still runs. Verified on planted synthetic signal in `tests/test_evoked_hypotheses.py`. +- ✅ **vertex_evoked (permutation adapter):** per-measure vertex maps through the same map+cluster + contract as `vertex_cluster` (band coordinate = measure name, dv = measure type). +- ⏳ **specials:** `vertex_signature` (decoding), `electrode_comparison` / `fcd_comparison` + (agreement — may not take hypotheses). `vertex_spatial` is retired (no inference). - ⏳ **migration / retirement:** move `study_treatment.yaml` to `design:`/`hypotheses:` and delete the retired gating code (`apply_hypothesis_gating`, `build_rescue_verdicts`, `gate_on`). diff --git a/pyproject.toml b/pyproject.toml index 683685a..faa32b1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,8 +20,17 @@ dependencies = [ # that exactly, and this floor merely lets the declaration resolve. "specparam>=2.0.0rc6", "joblib>=1.3", + # Figures (matplotlib glass brains, comparison plots) are part of the core + # lifecycle, so matplotlib is a core dependency, not an extra. + "matplotlib>=3.7", ] +# Optional extras. The package imports without any of them; the modules that +# need one raise a clear ImportError naming the extra to install: +# mne -> the evoked / TFR modules (roi_evoked, vertex_evoked, electrode_evoked) +# mvpa -> vertex_signature / electrode_signature (scikit-learn) +# network -> roi_graph / vertex_graph / *_network / *_nbs (networkx) +# atlas -> nibabel atlas readers [project.optional-dependencies] mne = ["mne>=1.5"] atlas = ["nibabel>=5.0"] @@ -37,6 +46,13 @@ source-analytics = "source_analytics.cli:main" [tool.setuptools.packages.find] where = ["src"] +# Ship the R scripts in wheels/sdists. They live outside the Python package +# (R/ beside src/), so they go in as data files under +# /share/source-analytics/R, where analyses.base.find_r_script_dir() +# looks after the repo checkout. MANIFEST.in grafts R/ into the sdist. +[tool.setuptools.data-files] +"share/source-analytics/R" = ["R/*.R"] + [tool.black] line-length = 99 diff --git a/run_connectivity_network.sh b/run_connectivity_network.sh deleted file mode 100755 index 5c58b90..0000000 --- a/run_connectivity_network.sh +++ /dev/null @@ -1,12 +0,0 @@ -#!/bin/bash -source /home/edm9fd/sandbox/source-analytics/.venv/bin/activate - -echo "=== Starting vertex_connectivity at $(date) ===" >> /home/edm9fd/sandbox/source-analytics/bg_analyses.log 2>&1 -source-analytics run --study /mnt/d/research/EEG/FORGE/analysis_wholebrain.yaml --analysis vertex_connectivity >> /home/edm9fd/sandbox/source-analytics/bg_analyses.log 2>&1 -echo "=== vertex_connectivity finished at $(date) ===" >> /home/edm9fd/sandbox/source-analytics/bg_analyses.log 2>&1 - -echo "=== Starting network at $(date) ===" >> /home/edm9fd/sandbox/source-analytics/bg_analyses.log 2>&1 -source-analytics run --study /mnt/d/research/EEG/FORGE/analysis_wholebrain.yaml --analysis network >> /home/edm9fd/sandbox/source-analytics/bg_analyses.log 2>&1 -echo "=== network finished at $(date) ===" >> /home/edm9fd/sandbox/source-analytics/bg_analyses.log 2>&1 - -echo "=== ALL DONE at $(date) ===" >> /home/edm9fd/sandbox/source-analytics/bg_analyses.log 2>&1 diff --git a/scripts/run_study.sh b/scripts/run_study.sh new file mode 100755 index 0000000..34ef1be --- /dev/null +++ b/scripts/run_study.sh @@ -0,0 +1,50 @@ +#!/usr/bin/env bash +# Run a full source-analytics study in dependency order (primaries before the +# supplements that read their output). This is the canonical recipe the README +# refers to; trim the lists to the analyses your config declares. +# +# scripts/run_study.sh study.yaml [extra `run` flags, e.g. --jobs -1 --steps ...] +# +# Paradigm names (resting / vertex / evoked) must match the `paradigms:` keys in +# your study YAML. fcd_comparison reads electrode_connectivity (resting) AND +# vertex_connectivity (vertex), so it runs last. Figures are OFF by default; add +# --steps setup,process,aggregate,statistics,figures,summary (or run +# `source-analytics figure ...`) to render them. +set -euo pipefail + +STUDY="${1:?usage: $0 study.yaml [run flags...]}"; shift || true +SA=(source-analytics run --study "$STUDY" "$@") + +run() { local paradigm="$1" analysis="$2"; echo "=== $paradigm / $analysis ==="; "${SA[@]}" --paradigm "$paradigm" --analysis "$analysis"; } + +# Resting paradigm — ROI + electrode +run resting roi_psd +run resting roi_aperiodic +run resting roi_connectivity # PRIMARY +run resting roi_graph # after roi_connectivity +run resting roi_nbs # after roi_connectivity +run resting roi_cross_freq # PAC + AAC + PPC (--metric to pick one) +run resting roi_directed # transfer entropy + DTF (--metric te|dtf) +run resting electrode_psd # PRIMARY +run resting electrode_aperiodic +run resting electrode_comparison # after electrode_psd AND roi_psd +run resting electrode_connectivity # sensor FC comparator +run resting electrode_signature # after electrode_psd + +# Vertex paradigm — whole-brain +run vertex vertex_connectivity # PRIMARY (slow; computes matrices) +run vertex vertex_graph # after vertex_connectivity +run vertex vertex_nbs # after vertex_connectivity +run vertex vertex_cluster +run vertex vertex_specparam +run vertex vertex_signature +run vertex vertex_cross_freq # local PAC + AAC + PPC +run vertex vertex_directed # vertex DTF (outflow/inflow/netflow) + +# Source-vs-sensor FCD comparison (cross-paradigm: reads resting + vertex output) +run resting fcd_comparison + +# Evoked paradigm (trial-based data only) — uncomment if your study has one +# run evoked roi_evoked +# run evoked vertex_evoked +# run evoked electrode_evoked diff --git a/src/source_analytics/analyses/_evoked_hypotheses.py b/src/source_analytics/analyses/_evoked_hypotheses.py new file mode 100644 index 0000000..d20a7ad --- /dev/null +++ b/src/source_analytics/analyses/_evoked_hypotheses.py @@ -0,0 +1,132 @@ +"""Declared-hypothesis wiring shared by the evoked modules. + +``roi_evoked`` and ``electrode_evoked`` are near-twins: the same measure +extraction over a different spatial unit (ROI vs channel), writing the same +row schema. They therefore need the same hypothesis wiring, and it lives here +so the two cannot drift — the same reason ``_network_base`` exists for the +graph/NBS pair. + +**Why these two were deferred.** ``docs/methods/HYPOTHESIS.md`` §9 lists +``roi_evoked`` / ``electrode_evoked`` as deferred with the reason "long-format +DV". Every other emmeans-tabular module exports one column per dependent +variable (``roi_psd`` has a column per band-power), so the module hands +``write_module_hypotheses`` a ``dv_cols`` vector. The evoked modules instead +export a single ``value`` column with ``measure_name`` telling you what it +holds. There is exactly one DV and 20-odd measures inside it. + +**The resolution is that the measure is a FACET, not a DV.** The tabular adapter +already carries ``facet_cols``, which runs an independent FDR family per facet +combination across the band x spatial grid. Pointing that at ``measure_name`` +gives one family per measure across the ROI/channel grid, which is the only +defensible family here: the measures are on incomparable scales (ITC is 0-1, +ERSP is dB, ERP amplitude is signal units, ERP latency is seconds), so a family +spanning them would pool quantities with no common null. The declared ``fdr:`` +scope still modulates *within* a measure. + +**There is no band axis.** Each measure definition already fixes its own band +and time window — ``band_lo``/``band_hi`` are properties of the measure, not a +factor to test across — so no ``band`` column is exported and the band +coordinate stays null. That is a genuine absence, not a dropped axis. + +This is **additive**. It does not replace the descriptive ``group * roi`` LMM +the R scripts run: that is a whole-brain omnibus over every spatial unit, while +this is the a-priori path, one declared contrast/omnibus/equivalence per +``hypotheses:`` entry. Both tables are written. +""" + +from __future__ import annotations + +import logging +from pathlib import Path + +import pandas as pd + +logger = logging.getLogger(__name__) + +# Every evoked measure row carries the group under this key, whatever the study +# calls its design factor. +_ROW_GROUP_KEY = "group" + + +def evoked_hypothesis_frame(analysis, measures_csv: str) -> pd.DataFrame | None: + """The measure table to test, from memory or from the exported CSV. + + Prefers the in-memory rows, so a full run needs no round-trip. Falls back to + the module's exported measures CSV so ``--steps statistics`` works on its own + against an earlier run's export. + + Returns ``None`` (having logged why) when there is nothing testable. + """ + rows = getattr(analysis, "_measure_rows", None) + if rows: + df = pd.DataFrame(rows) + else: + csv = Path(analysis.output_dir) / "data" / measures_csv + if not csv.exists(): + logger.warning( + "No measure rows in memory and no %s — skipping %s hypotheses.", + csv.name, analysis.name, + ) + return None + df = pd.read_csv(csv) + + if df.empty: + return None + + # The design factor is whatever the spec names it; rows always carry the + # group under "group", so alias when a study renames the factor. + factor = analysis.config.design_spec.factor + if factor not in df.columns: + if _ROW_GROUP_KEY not in df.columns: + logger.warning( + "%s measure table has neither '%s' nor '%s' — skipping hypotheses.", + analysis.name, factor, _ROW_GROUP_KEY, + ) + return None + df[factor] = df[_ROW_GROUP_KEY] + return df + + +def write_evoked_hypotheses(analysis, *, spatial_col: str, measures_csv: str): + """Run every declared hypothesis over an evoked measure table. + + Writes the additive ``_hypotheses.csv``. No-op when nothing is + declared. See the module docstring for why the measure is the facet. + + Parameters + ---------- + analysis + The calling :class:`BaseAnalysis`; supplies config, ``tbl_dir``, ``name`` + and the ``--hypothesis`` selection. + spatial_col + The spatial unit column — ``"roi"`` or ``"channel"``. + measures_csv + Basename of the module's exported measures CSV, used as the fallback + source when no rows are in memory. + """ + spec = analysis.config.design_spec + if spec is None or not spec.hypotheses: + return None + + df = evoked_hypothesis_frame(analysis, measures_csv) + if df is None: + return None + if spatial_col not in df.columns: + logger.warning( + "%s measure table has no '%s' column — skipping hypotheses.", + analysis.name, spatial_col, + ) + return None + + from ..hypothesis import write_module_hypotheses_tabular + + wanted = analysis._selection.get("hypothesis") + return write_module_hypotheses_tabular( + df, analysis.config, analysis.tbl_dir, prefix=analysis.name, + value_col="value", spatial_col=spatial_col, + facet_cols=("measure_name",), + # No band axis: absent from the frame, so the adapter runs a single + # null-band pass rather than iterating bands. + band_col="band", + hypothesis=",".join(sorted(wanted)) if wanted else None, + ) diff --git a/src/source_analytics/analyses/base.py b/src/source_analytics/analyses/base.py index 019f973..5b2c4aa 100644 --- a/src/source_analytics/analyses/base.py +++ b/src/source_analytics/analyses/base.py @@ -17,20 +17,32 @@ def find_r_script_dir() -> Path: - """Locate the R/ directory relative to this package. + """Locate the R/ scripts directory. - Searches upward from the analyses/ directory to find the R/ scripts - directory that lives at the package root (sibling to src/). + Checked in order: ``SOURCE_ANALYTICS_R_DIR`` (env override), the repo + checkout (``R/`` beside ``src/``, the editable-install case), the wheel's + data files (``/share/source-analytics/R``, see + ``[tool.setuptools.data-files]`` in pyproject.toml), then ``./R``. """ + import os + import sys + + env = os.environ.get("SOURCE_ANALYTICS_R_DIR") + candidates: list[Path] = [Path(env)] if env else [] pkg_root = Path(__file__).resolve().parent.parent.parent.parent # src/../.. - r_dir = pkg_root / "R" - if r_dir.is_dir(): - return r_dir - for candidate in [Path.cwd() / "R", Path(__file__).parent.parent.parent / "R"]: + candidates += [ + pkg_root / "R", + Path(sys.prefix) / "share" / "source-analytics" / "R", + Path(__file__).parent.parent.parent / "R", + Path.cwd() / "R", + ] + for candidate in candidates: if candidate.is_dir(): return candidate raise FileNotFoundError( - "Cannot find R/ scripts directory. Expected at: " + str(pkg_root / "R") + "Cannot find the R/ scripts directory. Looked in: " + + ", ".join(str(c) for c in candidates) + + ". Set SOURCE_ANALYTICS_R_DIR to point at it." ) @@ -295,21 +307,42 @@ def _parallel_capable(self) -> bool: """True when this module overrides the parallel per-subject compute.""" return type(self)._compute_subject is not BaseAnalysis._compute_subject - def _resolve_jobs(self, jobs: int) -> int: - """Effective worker count. CLI ``jobs`` wins; else ``jobs:`` in the study - config; ``-1``/``0`` → all but one core. Clamped to [1, #subjects].""" + def _resolve_jobs(self, jobs: int | None) -> int: + """Effective worker count. + + An explicit CLI ``--jobs N`` (any integer, including ``1``) wins. When + the CLI gave nothing (``None``), ``jobs:`` in the study config is used, + else serial. ``0``/``-1`` from either source mean "all but one core". + """ import os n = jobs - if n in (None, 1): - cfg_n = self.config.raw.get("jobs") - if cfg_n is not None: - n = int(cfg_n) if n is None: - n = 1 + cfg_n = self.config.raw.get("jobs") + n = int(cfg_n) if cfg_n is not None else 1 + n = int(n) if n <= 0: # -1 / 0 → auto: leave one core free n = max(1, (os.cpu_count() or 1) - 1) - return max(1, int(n)) + return max(1, n) + + def _vertex_epoch_config(self) -> dict | None: + """Merged epoch-sampling config for the vertex modules, or None if disabled. + + Precedence (lowest → highest): top-level ``epoch_sampling:`` → the + ``vertex:`` block's ``epoch_sampling:`` → this analysis's own + ``epoch_sampling:`` (under ``paradigms.

.analyses.``). The ROI + and electrode modules already merge global + per-analysis in + ``__init__``; the vertex modules used to read only the ``vertex:`` + block, so the documented global default never reached them. + """ + raw = self.config.raw or {} + global_cfg = raw.get("epoch_sampling") or {} + vertex_cfg = (self.config.vertex or {}).get("epoch_sampling") or {} + analysis_cfg = (raw.get(self.name) or {}).get("epoch_sampling") or {} + merged = {**global_cfg, **vertex_cfg, **analysis_cfg} + if not merged.get("enabled", False): + return None + return merged def _process_subjects_parallel(self, subjects: list[SubjectInfo], n_jobs: int) -> None: """Compute per-subject payloads in worker processes, merge in the parent. @@ -546,7 +579,7 @@ def run( subjects: list[SubjectInfo], steps: set[str] | None = None, select: dict[str, frozenset[str]] | None = None, - jobs: int = 1, + jobs: int | None = None, ) -> None: """Execute the analysis lifecycle. diff --git a/src/source_analytics/analyses/electrode_evoked_analysis.py b/src/source_analytics/analyses/electrode_evoked_analysis.py index 46ad7d5..aa94c11 100644 --- a/src/source_analytics/analyses/electrode_evoked_analysis.py +++ b/src/source_analytics/analyses/electrode_evoked_analysis.py @@ -27,6 +27,7 @@ extract_measure_in_tiles, resolve_n_cycles, ) +from ._evoked_hypotheses import write_evoked_hypotheses from .base import BaseAnalysis, find_r_script_dir logger = logging.getLogger(__name__) @@ -46,6 +47,8 @@ class ElectrodeEvokedAnalysis(BaseAnalysis): name = "electrode_evoked" + SELECTABLE = {"hypothesis": "declared hypothesis"} + def __init__(self, config: StudyConfig, output_dir: Path): super().__init__(config, output_dir) self._measure_rows: list[dict] = [] @@ -358,8 +361,17 @@ def aggregate(self) -> None: ) def statistics(self) -> None: - """Delegated to R.""" - pass + """Declared hypotheses over the channel measure table (additive). + + Delegates to the shared evoked wiring — see + :mod:`._evoked_hypotheses` for why the measure is the facet and the + band coordinate is null. The descriptive LMM in + ``electrode_evoked_analysis.R`` is unaffected and still runs. + """ + write_evoked_hypotheses( + self, spatial_col="channel", + measures_csv="electrode_evoked_measures.csv", + ) def figures(self) -> None: """Regenerate R figures from existing data/tables.""" @@ -405,6 +417,10 @@ def summary(self) -> None: "--tbl-dir", str(self.tbl_dir), ] cmd.extend(self._r_no_figures_flags()) + # --hypothesis NAME[,NAME] narrows the R-side pairwise contrast list too. + wanted_hyp = self._selection.get("hypothesis") + if wanted_hyp: + cmd.extend(["--hypothesis", ",".join(sorted(wanted_hyp))]) logger.info("Calling R: %s", " ".join(cmd)) try: diff --git a/src/source_analytics/analyses/fcd_comparison_analysis.py b/src/source_analytics/analyses/fcd_comparison_analysis.py index 4fddc3e..1ebd5b2 100644 --- a/src/source_analytics/analyses/fcd_comparison_analysis.py +++ b/src/source_analytics/analyses/fcd_comparison_analysis.py @@ -82,18 +82,58 @@ def __init__(self, config: StudyConfig, output_dir: Path): cfg = config.raw.get(self.name, {}) self._metric_filter = cfg.get("metrics") - def setup(self) -> None: - base = self.config.output_dir - sensor_csv = base / "electrode_connectivity" / "data" / "electrode_fcd.csv" - source_csv = base / "vertex_connectivity" / "data" / "vertex_fcd.csv" - if not sensor_csv.exists(): - raise FileNotFoundError( - f"{sensor_csv} not found — run 'electrode_connectivity' first." - ) - if not source_csv.exists(): + def _find_upstream_csv(self, module: str, filename: str, override_key: str) -> Path: + """Locate a primary module's data CSV, in this paradigm or a sibling one. + + The two primaries usually live in *different* paradigms (the canonical + layout runs ``electrode_connectivity`` under ``resting`` and + ``vertex_connectivity`` under ``vertex``), so a same-paradigm lookup is + not enough. Search order: + + 1. an explicit ``: `` in this module's config block + (a directory containing ``data/`` or the CSV itself); + 2. this paradigm's working dir (``analytics///data``); + 3. every sibling paradigm dir under the same analytics root. + """ + cfg = self.config.raw.get(self.name, {}) or {} + override = cfg.get(override_key) + if override: + p = Path(override) + if not p.is_absolute(): + p = (self.config.output_dir / p).resolve() + candidate = p if p.suffix == ".csv" else p / "data" / filename + if candidate.exists(): + return candidate raise FileNotFoundError( - f"{source_csv} not found — run 'vertex_connectivity' first." + f"{self.name}: {override_key}={override!r} does not contain {filename}" ) + + base = self.config.output_dir + same = base / module / "data" / filename + if same.exists(): + return same + + searched = [same] + root = base.parent if self.config.paradigm_name else base + for sibling in sorted(p for p in root.iterdir() if p.is_dir()) if root.is_dir() else []: + candidate = sibling / module / "data" / filename + searched.append(candidate) + if candidate.exists(): + logger.info("%s: using %s from paradigm dir %s", self.name, module, sibling.name) + return candidate + raise FileNotFoundError( + f"{filename} not found — run '{module}' first (any paradigm). Searched: " + + ", ".join(str(s) for s in searched) + + f". Or set {self.name}.{override_key} to its output directory." + ) + + def setup(self) -> None: + sensor_csv = self._find_upstream_csv( + "electrode_connectivity", "electrode_fcd.csv", "sensor_dir", + ) + source_csv = self._find_upstream_csv( + "vertex_connectivity", "vertex_fcd.csv", "source_dir", + ) self._sensor_df = pd.read_csv(sensor_csv) self._source_df = pd.read_csv(source_csv) logger.info( diff --git a/src/source_analytics/analyses/roi_aperiodic_analysis.py b/src/source_analytics/analyses/roi_aperiodic_analysis.py index cfaac1c..2b26008 100644 --- a/src/source_analytics/analyses/roi_aperiodic_analysis.py +++ b/src/source_analytics/analyses/roi_aperiodic_analysis.py @@ -209,7 +209,7 @@ def summary(self) -> None: except FileNotFoundError: logger.error("Rscript not found. Install R to enable statistics and visualization.") except subprocess.TimeoutExpired: - logger.error("R script timed out after 600 seconds") + logger.error("R script timed out after 3600 seconds") # Render brain mosaics from posthoc effect sizes if self._generate_figures: diff --git a/src/source_analytics/analyses/roi_connectivity_analysis.py b/src/source_analytics/analyses/roi_connectivity_analysis.py index 51596f5..1e2a9ef 100644 --- a/src/source_analytics/analyses/roi_connectivity_analysis.py +++ b/src/source_analytics/analyses/roi_connectivity_analysis.py @@ -63,8 +63,22 @@ def __init__(self, config: StudyConfig, output_dir: Path): self._metrics: list[str] = list(self._ROI_METRICS) def setup(self) -> None: - # Restrict emitted metrics to --metric / --select metric=... if given. - self._metrics = self._select("metric", self._ROI_METRICS) + # Configured set: the module's `metrics:` list (under + # paradigms.

.analyses.roi_connectivity), else every ROI metric — the + # same contract vertex_connectivity honours. --metric / --select then + # narrows within that set. + cfg_metrics = (self.config.raw.get(self.name) or {}).get("metrics") + if cfg_metrics: + unknown = [m for m in cfg_metrics if m not in self._ROI_METRICS] + if unknown: + raise ValueError( + f"roi_connectivity: unknown metrics in config {unknown}; " + f"known: {list(self._ROI_METRICS)}" + ) + configured = [m for m in self._ROI_METRICS if m in set(cfg_metrics)] + else: + configured = list(self._ROI_METRICS) + self._metrics = self._select("metric", configured) self._edge_rows.clear() def _compute_subject(self, subject: SubjectInfo): @@ -345,4 +359,4 @@ def summary(self) -> None: "Rscript not found. Install R to enable statistics and visualization." ) except subprocess.TimeoutExpired: - logger.error("R script timed out after 600 seconds") + logger.error("R script timed out after 3600 seconds") diff --git a/src/source_analytics/analyses/roi_cross_freq_analysis.py b/src/source_analytics/analyses/roi_cross_freq_analysis.py index 00cad66..60cb423 100644 --- a/src/source_analytics/analyses/roi_cross_freq_analysis.py +++ b/src/source_analytics/analyses/roi_cross_freq_analysis.py @@ -10,8 +10,10 @@ All three share the signed (phase-preserving) ROI front-end and the same valid cross-frequency band pairs (``get_valid_pac_pairs``). PAC keeps its R statistics -+ brain mosaics; AAC/PPC emit edge CSVs (group statistics ride the connectivity -/ gating R path). ++ brain mosaics (``roi_pac_analysis.R``); AAC/PPC emit ROI×ROI edge CSVs whose +declared hypotheses are tested by ``roi_cross_freq_edges_analysis.R`` at the +global, cell and directed-region-pair tiers +(``roi_cross_freq_{aac,ppc}_*_hypotheses.csv``). """ from __future__ import annotations @@ -230,39 +232,67 @@ def aggregate(self) -> None: logger.info("Exported %s_edges.csv (%d rows)", metric, len(df)) def statistics(self) -> None: - """Delegated to R (PAC) / gating engine (AAC, PPC).""" + """Delegated to R: PAC via ``roi_pac_analysis.R``, AAC/PPC via + ``roi_cross_freq_edges_analysis.R`` (both run from :meth:`summary`).""" pass + def _edge_metrics_on_disk(self) -> list[str]: + """The AAC/PPC metrics selected for this run whose edge CSV exists.""" + data_dir = self.output_dir / "data" + return [ + m for m in ("aac", "ppc") + if m in self._metrics and (data_dir / f"{m}_edges.csv").exists() + ] + def figures(self) -> None: - """Regenerate PAC R figures from existing data/tables.""" + """Regenerate the R figures from existing data/tables.""" if "pac" in self._metrics: self._call_r_figures_only("roi_pac_analysis.R", "pac_values.csv") + for m in self._edge_metrics_on_disk(): + self._call_r_figures_only("roi_cross_freq_edges_analysis.R", f"{m}_edges.csv") + break # one call handles every edge CSV present # -------------------------------------------------------------- summary def summary(self) -> None: - """PAC statistics + figures via R; AAC/PPC emit matrices only (for now).""" + """Statistics + figures + report via R (PAC first, then AAC/PPC).""" if "pac" in self._metrics: self._run_pac_r() - if self._aac_rows or self._ppc_rows: - logger.info( - "AAC/PPC edge CSVs written; group statistics run via the " - "connectivity/gating R path (no dedicated R script yet)." - ) + edge_metrics = self._edge_metrics_on_disk() + if edge_metrics: + self._run_edges_r(edge_metrics) + + def _run_edges_r(self, metrics: list[str]) -> None: + """AAC/PPC hypotheses (global / cells / region pairs) via R.""" + self._run_r_script( + "roi_cross_freq_edges_analysis.R", + extra_args=["--metric", ",".join(metrics)], + label="AAC/PPC", + ) def _run_pac_r(self) -> None: data_dir = self.output_dir / "data" if not (data_dir / "pac_values.csv").exists(): logger.error("pac_values.csv not found -- skipping PAC R analysis") return + if self._run_r_script("roi_pac_analysis.R", label="PAC") and self._generate_figures: + self._render_brain_mosaics() + + def _run_r_script(self, script_name: str, *, extra_args: list[str] | None = None, + label: str = "R", timeout: int = 3600) -> bool: + """Run one of this module's R scripts with the standard argument set. + + Returns True when the script exited 0. + """ + data_dir = self.output_dir / "data" try: r_dir = _find_r_script_dir() except FileNotFoundError as e: logger.error(str(e)) - return - r_script = r_dir / "roi_pac_analysis.R" + return False + r_script = r_dir / script_name if not r_script.exists(): logger.error("R script not found: %s", r_script) - return + return False config_path = data_dir / "study_config.yaml" config_data = dict(self.config.raw) @@ -287,25 +317,26 @@ def _run_pac_r(self) -> None: wanted_hyp = self._selection.get("hypothesis") if wanted_hyp: cmd.extend(["--hypothesis", ",".join(sorted(wanted_hyp))]) + if extra_args: + cmd.extend(extra_args) - logger.info("Calling R: %s", " ".join(cmd)) + logger.info("Calling R (%s): %s", label, " ".join(cmd)) try: - result = subprocess.run(cmd, capture_output=True, text=True, timeout=3600) + result = subprocess.run(cmd, capture_output=True, text=True, timeout=timeout) for stream in (result.stdout, result.stderr): if stream: for line in stream.strip().split("\n"): if line.strip(): logger.info("[R] %s", line) if result.returncode != 0: - logger.error("R script failed with exit code %d", result.returncode) + logger.error("%s R script failed with exit code %d", label, result.returncode) + return False + return True except FileNotFoundError: - logger.error("Rscript not found. Install R to enable PAC statistics.") + logger.error("Rscript not found. Install R to enable %s statistics.", label) except subprocess.TimeoutExpired: - logger.error("PAC R script timed out") - return - - if self._generate_figures: - self._render_brain_mosaics() + logger.error("%s R script timed out after %d seconds", label, timeout) + return False def _render_brain_mosaics(self) -> None: """Render brain ROI mosaics from PAC region-level posthoc CSVs.""" diff --git a/src/source_analytics/analyses/roi_directed_analysis.py b/src/source_analytics/analyses/roi_directed_analysis.py index 446dc47..258095e 100644 --- a/src/source_analytics/analyses/roi_directed_analysis.py +++ b/src/source_analytics/analyses/roi_directed_analysis.py @@ -41,15 +41,17 @@ def _find_r_script_dir() -> Path: class ROIDirectedAnalysis(BaseAnalysis): - """ROI-level directed connectivity (transfer entropy; DTF planned). + """ROI-level directed connectivity (transfer entropy + DTF). Uses **signed** (phase-preserving) ROI timeseries to compute binned transfer entropy for all n*(n-1) directed ROI pairs (40 brain ROIs → 1,560 directed pairs; 6 corpus callosum white matter tracts excluded). - Python computes directed matrices and exports directed edge-level CSV. - R (lme4, ggplot2) handles global t-tests, directional paired t-tests, - region-pair LMM, and summary report. + Python computes directed matrices and exports a directed edge-level CSV + (``roi_transfer_entropy_edges.csv``, one column per measure). R runs the + declared hypotheses over every directed DV present (``te``, ``net_te``, + ``dtf``) at the global, directed-edge and region-pair tiers, writing + ``roi_directed_*_hypotheses.csv`` tables and the summary report. ``--metric`` selects which directed measure(s) to compute: ``te`` (transfer entropy, which also yields ``net_te``) and/or ``dtf`` @@ -184,18 +186,12 @@ def figures(self) -> None: def summary(self) -> None: """Call Rscript for statistics and summary report. - The R stats currently cover transfer entropy; a DTF-only run skips R - (the ``dtf`` column is still exported in the edge CSV for downstream use). + The R script tests whichever directed DV columns the edge CSV carries + (``te``/``net_te`` and/or ``dtf``), so a ``--metric dtf`` run gets the + same hypothesis tables as a TE run. """ data_dir = self.output_dir / "data" - if "te" not in self._metrics: - logger.info( - "roi_directed: metrics=%s — DTF has no R stats yet; edge CSV written, " - "skipping R.", self._metrics, - ) - return - if not (data_dir / "roi_transfer_entropy_edges.csv").exists(): logger.error("roi_transfer_entropy_edges.csv not found -- skipping R analysis") return @@ -238,13 +234,14 @@ def summary(self) -> None: if wanted_hyp: cmd.extend(["--hypothesis", ",".join(sorted(wanted_hyp))]) + r_timeout = 3600 logger.info("Calling R: %s", " ".join(cmd)) try: result = subprocess.run( cmd, capture_output=True, text=True, - timeout=3600, + timeout=r_timeout, ) if result.stdout: for line in result.stdout.strip().split("\n"): @@ -260,4 +257,4 @@ def summary(self) -> None: "Rscript not found. Install R to enable statistics and visualization." ) except subprocess.TimeoutExpired: - logger.error("R script timed out after 600 seconds") + logger.error("R script timed out after %d seconds", r_timeout) diff --git a/src/source_analytics/analyses/roi_evoked_analysis.py b/src/source_analytics/analyses/roi_evoked_analysis.py index d28f9a7..938579c 100644 --- a/src/source_analytics/analyses/roi_evoked_analysis.py +++ b/src/source_analytics/analyses/roi_evoked_analysis.py @@ -22,6 +22,7 @@ resolve_n_cycles, extract_measure_in_tiles, ) +from ._evoked_hypotheses import write_evoked_hypotheses from .base import BaseAnalysis, find_r_script_dir logger = logging.getLogger(__name__) @@ -39,6 +40,8 @@ class ROIEvokedAnalysis(BaseAnalysis): name = "roi_evoked" + SELECTABLE = {"hypothesis": "declared hypothesis"} + def __init__(self, config: StudyConfig, output_dir: Path): super().__init__(config, output_dir) self._measure_rows: list[dict] = [] @@ -255,8 +258,16 @@ def aggregate(self) -> None: logger.info("Exported roi_evoked_tfr.csv (%d rows)", len(tfr_df)) def statistics(self) -> None: - """Delegated to R.""" - pass + """Declared hypotheses over the ROI measure table (additive). + + Delegates to the shared evoked wiring — see + :mod:`._evoked_hypotheses` for why the measure is the facet and the + band coordinate is null. The descriptive ``group * roi`` LMM in + ``roi_evoked_analysis.R`` is unaffected and still runs. + """ + write_evoked_hypotheses( + self, spatial_col="roi", measures_csv="roi_evoked_measures.csv" + ) def figures(self) -> None: """Regenerate R figures from existing data/tables.""" @@ -299,6 +310,10 @@ def summary(self) -> None: ] cmd.extend(self._r_no_figures_flags()) cmd.extend(self._r_roi_categories_flags()) + # --hypothesis NAME[,NAME] narrows the R-side pairwise contrast list too. + wanted_hyp = self._selection.get("hypothesis") + if wanted_hyp: + cmd.extend(["--hypothesis", ",".join(sorted(wanted_hyp))]) logger.info("Calling R: %s", " ".join(cmd)) try: diff --git a/src/source_analytics/analyses/roi_network_analysis.py b/src/source_analytics/analyses/roi_network_analysis.py index 1beec7f..262736c 100644 --- a/src/source_analytics/analyses/roi_network_analysis.py +++ b/src/source_analytics/analyses/roi_network_analysis.py @@ -15,7 +15,6 @@ from __future__ import annotations import logging -import subprocess from pathlib import Path import numpy as np @@ -27,7 +26,6 @@ from ..stats.graph_metrics import compute_graph_metrics from ..viz.constants import CC_ROIS, METRIC_LABELS from ._network_base import NetworkAnalysisBase -from .base import find_r_script_dir logger = logging.getLogger(__name__) @@ -496,21 +494,10 @@ def figures(self) -> None: self._graph_figures() def summary(self) -> None: - # Keep the richer R report when available; else the python summary. + # The combined alias has no R report of its own (graph + NBS statistics + # run in Python via the hypothesis layer); write the Python summary. data_dir = self.output_dir / "data" config_path = data_dir / "study_config.yaml" with open(config_path, "w") as f: yaml.dump(dict(self.config.raw), f, default_flow_style=False) - try: - r_script = find_r_script_dir() / "roi_network_analysis.R" - if r_script.exists(): - cmd = ["Rscript", str(r_script), "--data-dir", str(data_dir), - "--config", str(config_path), "--output-dir", str(self.output_dir), - "--fig-dir", str(self.fig_dir), "--tbl-dir", str(self.tbl_dir)] - cmd.extend(self._r_no_figures_flags()) - cmd.extend(self._r_roi_categories_flags()) - if subprocess.run(cmd, capture_output=True, text=True, timeout=3600).returncode == 0: - return - except (FileNotFoundError, subprocess.TimeoutExpired): - pass self._write_summary(graph=True, nbs=True) diff --git a/src/source_analytics/analyses/roi_psd_analysis.py b/src/source_analytics/analyses/roi_psd_analysis.py index fc32d3f..c2b016b 100644 --- a/src/source_analytics/analyses/roi_psd_analysis.py +++ b/src/source_analytics/analyses/roi_psd_analysis.py @@ -245,7 +245,7 @@ def summary(self) -> None: except FileNotFoundError: logger.error("Rscript not found. Install R to enable statistics and visualization.") except subprocess.TimeoutExpired: - logger.error("R script timed out after 600 seconds") + logger.error("R script timed out after 3600 seconds") # Render brain mosaics from posthoc effect sizes if self._generate_figures: diff --git a/src/source_analytics/analyses/vertex_cluster_analysis.py b/src/source_analytics/analyses/vertex_cluster_analysis.py index 63604f1..73fd613 100644 --- a/src/source_analytics/analyses/vertex_cluster_analysis.py +++ b/src/source_analytics/analyses/vertex_cluster_analysis.py @@ -29,6 +29,7 @@ compute_spectral_slope, compute_peak_frequency, ) +from ..spectral.epoch_sampler import sample_epochs from ..stats.cluster_permutation import ( cluster_permutation_test, has_significant_cluster as _has_significant_cluster, @@ -105,6 +106,10 @@ def __init__(self, config: StudyConfig, output_dir: Path): ) self._correction_method = "cluster" + # Global epoch_sampling → vertex: block → per-analysis block (see base). + # None = full continuous PSD (the historical vertex_cluster behaviour). + self._epoch_config = self._vertex_epoch_config() + def setup(self) -> None: self._band_power_rows.clear() self._feature_rows.clear() @@ -131,7 +136,21 @@ def _compute_subject(self, subject: SubjectInfo): # Compute PSD for all vertices: (n_vertices, n_freqs) fmax = max(hi for _, hi in self.config.bands.values()) + 10 - freqs, psd = compute_psd_vertices(stc_data, sfreq, fmax=fmax) + if self._epoch_config is not None: + epochs = sample_epochs( + stc_data, sfreq, + epoch_duration_sec=self._epoch_config.get("epoch_duration_sec", 2.0), + n_epochs=self._epoch_config.get("n_epochs", 80), + seed=self._epoch_config.get("seed", 42), + n_bootstrap=self._epoch_config.get("n_bootstrap", 1), + ) + all_psd = [] + for ep in epochs: + freqs, p = compute_psd_vertices(ep, sfreq, fmax=fmax) + all_psd.append(p) + psd = np.mean(all_psd, axis=0) + else: + freqs, psd = compute_psd_vertices(stc_data, sfreq, fmax=fmax) # Extract band power metrics band_power = extract_band_power_vertices( @@ -694,7 +713,7 @@ def summary(self) -> None: logger.warning("Rscript not found — writing Python summary") self._write_python_summary() except subprocess.TimeoutExpired: - logger.error("R script timed out after 600 seconds") + logger.error("R script timed out after 3600 seconds") self._write_python_summary() def _write_python_summary(self) -> None: diff --git a/src/source_analytics/analyses/vertex_connectivity_analysis.py b/src/source_analytics/analyses/vertex_connectivity_analysis.py index 3a8b338..1e5d8ff 100644 --- a/src/source_analytics/analyses/vertex_connectivity_analysis.py +++ b/src/source_analytics/analyses/vertex_connectivity_analysis.py @@ -27,7 +27,7 @@ compute_fcd, FCD_CENTER, ) -from ..spectral.epoch_sampler import sample_epochs, get_epoch_config +from ..spectral.epoch_sampler import sample_epochs from ..stats.cluster_permutation import ( cluster_permutation_test, has_significant_cluster as _has_significant_cluster, @@ -73,7 +73,8 @@ def __init__(self, config: StudyConfig, output_dir: Path): self._adjacency_distance = float(wb_cfg.get("adjacency_distance_mm", 5.0)) self._cluster_threshold = float(wb_cfg.get("cluster_threshold", 2.0)) - self._epoch_config = get_epoch_config(wb_cfg) + # Global epoch_sampling → vertex: block → per-analysis block (see base). + self._epoch_config = self._vertex_epoch_config() self._cluster_results: dict = {} diff --git a/src/source_analytics/analyses/vertex_evoked_analysis.py b/src/source_analytics/analyses/vertex_evoked_analysis.py index e9664f1..24a315d 100644 --- a/src/source_analytics/analyses/vertex_evoked_analysis.py +++ b/src/source_analytics/analyses/vertex_evoked_analysis.py @@ -46,6 +46,7 @@ class VertexEvokedAnalysis(BaseAnalysis): """ name = "vertex_evoked" + SELECTABLE = {"hypothesis": "declared hypothesis"} def __init__(self, config: StudyConfig, output_dir: Path): super().__init__(config, output_dir) @@ -289,6 +290,36 @@ def statistics(self) -> None: ) logger.info("Exported vertex_evoked_stats.csv (%d rows)", len(all_stats)) + # --- Declared hypotheses (hypothesis layer; additive, map+cluster contract) --- + # Same wiring as vertex_cluster: every declared hypothesis runs over the + # per-subject vertex maps through the permutation adapter and lands in + # tables/vertex_evoked_hypotheses.csv. The measure is the cell's "band" + # coordinate (each measure fixes its own band + time window, so there is + # no separate band axis — see analyses/_evoked_hypotheses.py) and the + # dv is the measure type (itc / ersp / stp / ...). --hypothesis narrows. + from ..hypothesis import write_module_hypotheses_perm + + measure_types = { + m["measure_name"]: m["measure_type"] for m in self._measure_rows + } + maps_by_cell = { + (mname, measure_types.get(mname, "value")): { + uid: self._subject_measures[uid][mname] + for uid in self._subject_groups + if mname in self._subject_measures.get(uid, {}) + } + for mname in measure_names + } + wanted_hyp = self._selection.get("hypothesis") + write_module_hypotheses_perm( + maps_by_cell, self._subject_groups, coords, self.config, self.tbl_dir, + prefix="vertex_evoked", + n_perms=self._n_permutations, threshold=self._cluster_threshold, + distance_mm=self._adjacency_distance, + hypothesis=",".join(sorted(wanted_hyp)) if wanted_hyp else None, + atlas_dir=self._atlas_dir, + ) + self._save_cluster_state() def figures(self) -> None: @@ -325,7 +356,9 @@ def summary(self) -> None: "", "- `data/vertex_evoked_measures.csv` — per-subject per-vertex measures", "- `data/source_coords.csv` — vertex coordinates (mm)", - "- `tables/vertex_evoked_stats.csv` — cluster-permutation statistics", + "- `tables/vertex_evoked_stats.csv` — cluster-permutation statistics (pairwise contrasts)", + "- `tables/vertex_evoked_hypotheses.csv` — declared hypotheses (permutation adapter; " + "band = measure name, dv = measure type)", "- `figures/evoked_*.png` — glass-brain measure maps", "", ] diff --git a/src/source_analytics/analyses/vertex_network_analysis.py b/src/source_analytics/analyses/vertex_network_analysis.py index db1b8ff..75180e5 100644 --- a/src/source_analytics/analyses/vertex_network_analysis.py +++ b/src/source_analytics/analyses/vertex_network_analysis.py @@ -17,7 +17,6 @@ import logging import pickle -import subprocess from functools import lru_cache from pathlib import Path @@ -32,7 +31,7 @@ compute_vertex_connectivity_matrix, compute_vertex_connectivity_matrix_epochs, ) -from ..spectral.epoch_sampler import sample_epochs, get_epoch_config +from ..spectral.epoch_sampler import sample_epochs from ..stats.graph_metrics import ( GLOBAL_METRIC_NAMES, compute_auc, @@ -44,7 +43,6 @@ plot_nbs_roi_circos, ) from ._network_base import NetworkAnalysisBase -from .base import find_r_script_dir logger = logging.getLogger(__name__) @@ -111,7 +109,8 @@ def __init__(self, config: StudyConfig, output_dir: Path): self._density_min = float(cfg.get("density_min", 0.05)) self._density_max = float(cfg.get("density_max", 0.40)) self._density_step = float(cfg.get("density_step", 0.01)) - self._epoch_config = get_epoch_config(config.vertex) + # Global epoch_sampling → vertex: block → per-analysis block (see base). + self._epoch_config = self._vertex_epoch_config() def setup(self) -> None: self._auc_rows.clear() @@ -583,15 +582,6 @@ def summary(self) -> None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: yaml.dump(config_data, f, default_flow_style=False) - try: - r_script = find_r_script_dir() / "vertex_network_analysis.R" - if r_script.exists(): - cmd = ["Rscript", str(r_script), "--data-dir", str(data_dir), - "--config", str(config_path), "--output-dir", str(self.output_dir), - "--fig-dir", str(self.fig_dir), "--tbl-dir", str(self.tbl_dir)] - cmd.extend(self._r_no_figures_flags()) - if subprocess.run(cmd, capture_output=True, text=True, timeout=3600).returncode == 0: - return - except (FileNotFoundError, subprocess.TimeoutExpired): - pass + # The combined alias has no R report of its own (graph + NBS statistics + # run in Python via the hypothesis layer); write the Python summary. self._write_summary(graph=True, nbs=True) diff --git a/src/source_analytics/analyses/vertex_signature_analysis.py b/src/source_analytics/analyses/vertex_signature_analysis.py index 3e52b13..b93ab24 100644 --- a/src/source_analytics/analyses/vertex_signature_analysis.py +++ b/src/source_analytics/analyses/vertex_signature_analysis.py @@ -26,7 +26,7 @@ from ..io.discovery import SubjectInfo from ..io.loader import SubjectLoader from ..spectral.vertex import compute_psd_vertices, extract_band_power_vertices -from ..spectral.epoch_sampler import sample_epochs, get_epoch_config +from ..spectral.epoch_sampler import sample_epochs from ..stats.signature import ( SignatureResult, classifier_label, @@ -87,7 +87,8 @@ def __init__(self, config: StudyConfig, output_dir: Path): if self._noise_exclude is not None: self._noise_exclude = tuple(self._noise_exclude) - self._epoch_config = get_epoch_config(wb_cfg) + # Global epoch_sampling → vertex: block → per-analysis block (see base). + self._epoch_config = self._vertex_epoch_config() self._signature_results: dict[str, object] = {} def setup(self) -> None: diff --git a/src/source_analytics/analyses/vertex_spatial_analysis.py b/src/source_analytics/analyses/vertex_spatial_analysis.py index 227dc9a..ecfde3e 100644 --- a/src/source_analytics/analyses/vertex_spatial_analysis.py +++ b/src/source_analytics/analyses/vertex_spatial_analysis.py @@ -1,340 +1,93 @@ -"""Vertex-level spatial analysis. - -Fits a single model per band using nlme::gls with exponential spatial -correlation structure, accounting for spatial autocorrelation and avoiding -the multiple comparison problem inherent in vertex-wise testing. - -The heavy lifting is done in R (nlme::gls). Python handles data preparation, -figure generation, and orchestration. +"""Vertex-level spatial analysis — RETIRED. + +This module used to fit a per-contrast spatial-covariance GLS (nlme::gls, +corExp + nugget) on per-vertex band power as a robustness check on the vertex +group difference. It iterated the legacy ``config$contrasts`` block, which the +declarative ``design:``/``hypotheses:`` spec no longer populates, and the +spatial-covariance table was never a manuscript result. It was retired rather +than migrated (2026-06): spatially-resolved vertex inference is delivered by +``vertex_cluster`` (cluster-based permutation glass-brain maps) and +``vertex_nbs`` (network-based statistic). + +The module is kept in the registry so old configs and ``--analysis +vertex_spatial`` (and its ``spatial_lmm`` alias) still resolve, but it does no +work: it neither loads source estimates nor calls R. It writes well-formed +empty result tables and a retirement note so downstream consumers (the +gallery) find the expected files, then exits cleanly. ``R/vertex_spatial_analysis.R`` +is retained for reference only and is not invoked. """ from __future__ import annotations import logging -import subprocess from pathlib import Path -import numpy as np import pandas as pd -import yaml from ..config import StudyConfig from ..io.discovery import SubjectInfo -from ..io.loader import SubjectLoader -from ..spectral.vertex import compute_psd_vertices, extract_band_power_vertices -from ..spectral.epoch_sampler import sample_epochs, get_epoch_config -from ..viz.glass_brain import plot_glass_brain, plot_anatomical_glass_brain from .base import BaseAnalysis logger = logging.getLogger(__name__) - -def _find_r_script_dir() -> Path: - pkg_root = Path(__file__).resolve().parent.parent.parent.parent - r_dir = pkg_root / "R" - if r_dir.is_dir(): - return r_dir - for candidate in [Path.cwd() / "R", Path(__file__).parent.parent.parent / "R"]: - if candidate.is_dir(): - return candidate - raise FileNotFoundError("Cannot find R/ scripts directory") +RETIRE_NOTE = ( + "vertex_spatial is RETIRED (design-spec migration, 2026-06). The per-contrast " + "GLS spatial-covariance robustness model iterated config$contrasts, which the " + "declarative design:/hypotheses: spec no longer populates. Spatially-resolved " + "vertex inference is provided by vertex_cluster (cluster-permutation glass-brain " + "maps) and vertex_nbs (network-based statistic). No data was processed." +) class VertexSpatialAnalysis(BaseAnalysis): - """Vertex-level spatial analysis (R-primary).""" + """RETIRED — writes empty result tables + a note; processes no subjects.""" name = "vertex_spatial" SELECTABLE = {"band": "frequency band"} def __init__(self, config: StudyConfig, output_dir: Path): super().__init__(config, output_dir) - self._power_rows: list[dict] = [] - self._source_coords: np.ndarray | None = None - self._vertex_indices: np.ndarray | None = None - self._sfreq: float | None = None - - # Config - slmm_cfg = config.raw.get("vertex_spatial", config.raw.get("spatial_lmm", {})) - self._stat_method = slmm_cfg.get("stat_method", "gls") - self._correlation_structure = slmm_cfg.get("correlation_structure", "exponential") - self._spatial_range_mm = float(slmm_cfg.get("spatial_range_mm", 3.0)) - - wb_cfg = config.vertex - self._noise_exclude = wb_cfg.get("noise_exclude_hz") - if self._noise_exclude is not None: - self._noise_exclude = tuple(self._noise_exclude) + self._warned = False - self._epoch_config = get_epoch_config(wb_cfg) + def _warn_once(self) -> None: + if not self._warned: + logger.warning(RETIRE_NOTE) + self._warned = True def setup(self) -> None: - self._power_rows.clear() - self._source_coords = None - self._vertex_indices = None - - def process_subject(self, subject: SubjectInfo) -> None: - loader = SubjectLoader(subject.data_dir) - uid = f"{subject.group}_{subject.subject_id}" - - stc_data = loader.load_source_timecourses() - sfreq = loader.load_sfreq() - coords = loader.load_source_coords() - - if self._sfreq is None: - self._sfreq = sfreq - - # Apply vertex filter (compute mask once from first subject) - if self._vertex_indices is None: - mask = self.config.get_vertex_mask(coords) - self._vertex_indices = np.where(mask)[0] - self._source_coords = coords[mask] - if self.config.has_vertex_filter: - logger.info( - "Vertex filter: %d/%d vertices retained", - len(self._vertex_indices), len(coords), - ) - - stc_data = stc_data[self._vertex_indices] - coords = self._source_coords + self._warn_once() - # Compute PSD - fmax = max(hi for _, hi in self.config.bands.values()) + 10 - if self._epoch_config is not None: - epochs = sample_epochs( - stc_data, sfreq, - epoch_duration_sec=self._epoch_config.get("epoch_duration_sec", 2.0), - n_epochs=self._epoch_config.get("n_epochs", 80), - seed=self._epoch_config.get("seed", 42), - n_bootstrap=self._epoch_config.get("n_bootstrap", 1), - ) - all_psd = [] - for ep in epochs: - f, p = compute_psd_vertices(ep, sfreq, fmax=fmax) - all_psd.append(p) - freqs = f - psd = np.mean(all_psd, axis=0) - else: - freqs, psd = compute_psd_vertices(stc_data, sfreq, fmax=fmax) + def process_subject(self, subject: SubjectInfo) -> None: # noqa: D401 + """No-op: the retired module loads nothing.""" - band_power = extract_band_power_vertices( - freqs, psd, self._selected_bands(), noise_exclude=self._noise_exclude, - ) - - n_vertices = stc_data.shape[0] - for band_name, bp in band_power.items(): - for vi in range(n_vertices): - self._power_rows.append({ - "subject": uid, - "group": subject.group, - "vertex_idx": int(self._vertex_indices[vi]), - "x": float(coords[vi, 0]), - "y": float(coords[vi, 1]), - "z": float(coords[vi, 2]), - "band": band_name, - "relative": float(bp["relative"][vi]), - "absolute": float(bp["absolute"][vi]), - }) - - def aggregate(self) -> None: - data_dir = self.output_dir / "data" - - power_df = pd.DataFrame(self._power_rows) - if power_df.empty: - logger.warning("No vertex spatial data collected") - return - power_df.to_csv(data_dir / "vertex_spatial_data.csv", index=False) - logger.info("Exported vertex_spatial_data.csv (%d rows)", len(power_df)) - - if self._source_coords is not None: - coords_df = pd.DataFrame(self._source_coords, columns=["x", "y", "z"]) - coords_df.index.name = "vertex_idx" - coords_df.to_csv(data_dir / "source_coords.csv") + def aggregate(self) -> None: # noqa: D401 + """No-op: nothing was computed.""" def statistics(self) -> None: - """Delegate all statistical analysis to R.""" - pass - - def figures(self) -> None: - """Generate residual maps and per-vertex t-score brain figures.""" - from scipy import stats as sp_stats - + """Write the well-formed empty tables downstream consumers expect.""" tbl_dir = self.tbl_dir - fig_dir = self.fig_dir - - residuals_csv = tbl_dir / "spatial_residuals.csv" - if residuals_csv.exists() and self._source_coords is not None: - resid_df = pd.read_csv(residuals_csv) - for band in resid_df["band"].unique(): - sub = resid_df[resid_df["band"] == band] - mean_resid = sub.groupby("vertex_idx")["residual"].mean().values - if len(mean_resid) == len(self._source_coords): - safe_name = band.lower().replace(" ", "_") - plot_glass_brain( - coords=self._source_coords, - values=mean_resid, - title=f"Spatial Residuals — {band}", - output_path=fig_dir / f"spatial_residuals_{safe_name}.png", - cmap="RdBu_r", - ) - - # Per-vertex t-score anatomical brain figures for gamma bands - data_csv = self.output_dir / "data" / "vertex_spatial_data.csv" - coords_csv = self.output_dir / "data" / "source_coords.csv" - if not data_csv.exists() or not coords_csv.exists(): - return - - df = pd.read_csv(data_csv) - coords_df = pd.read_csv(coords_csv) - coords = coords_df[["x", "y", "z"]].values - n_verts = len(coords) - vert_map = dict(zip(coords_df["vertex_idx"], range(n_verts))) - - all_bands = list(df["band"].unique()) - if not all_bands: - return + tbl_dir.mkdir(parents=True, exist_ok=True) + for name in ("vertex_spatial_results.csv", "vertex_spatial_residuals.csv"): + pd.DataFrame().to_csv(tbl_dir / name, index=False) + logger.info("vertex_spatial (retired): wrote empty result tables to %s", tbl_dir) - contrasts = self._pairwise_contrasts() - for contrast in contrasts: - ga, gb = contrast.group_a, contrast.group_b - contrast_name = f"{ga}_vs_{gb}" - - band_data = {} - for band in all_bands: - bdata = df[df["band"] == band] - t_values = np.zeros(n_verts) - for vidx in coords_df["vertex_idx"]: - vdata = bdata[bdata["vertex_idx"] == vidx] - a_vals = vdata[vdata["group"] == ga]["absolute"].values - b_vals = vdata[vdata["group"] == gb]["absolute"].values - if len(a_vals) >= 2 and len(b_vals) >= 2: - t_stat, _ = sp_stats.ttest_ind( - a_vals, b_vals, equal_var=False) - t_values[vert_map[vidx]] = t_stat - - n_sig = int((np.abs(t_values) > 2.0).sum()) - band_data[band] = { - "values": t_values, - "n_sig": n_sig, - "n_total": n_verts, - } - logger.info(" %s %s: %d/%d vertices |t|>2.0", - contrast_name, band, n_sig, n_verts) - - # Symmetric vlim for diverging colormap - all_t = np.concatenate([bd["values"] for bd in band_data.values()]) - vmax = float(np.ceil(np.nanmax(np.abs(all_t)) * 10) / 10) - if vmax == 0: - continue - mean_dir = "lower" if np.mean(all_t) < 0 else "higher" - - plot_anatomical_glass_brain( - coords=coords, - band_data=band_data, - output_path=fig_dir / f"band_tscores_{contrast_name}.png", - title=(f"Dorsal Vertex-Level Band Power: " - f"{ga} vs {gb} (vertex spatial posthoc)"), - subtitle=(f"Large circles = uncorrected |t| > 2.0; " - f"Mean direction: {ga} {mean_dir} than {gb}"), - cmap="RdBu_r", - vlim=(-vmax, vmax), - sig_threshold=2.0, - colorbar_label=f"t-statistic ({ga} vs {gb})", - dpi=300, - ) + def figures(self) -> None: # noqa: D401 + """No figures: there are no results.""" def summary(self) -> None: - """Run R script for spatial GLS fitting + report generation.""" - data_dir = self.output_dir / "data" - - config_path = data_dir / "study_config.yaml" - config_data = dict(self.config.raw) - if self._sfreq is not None: - config_data["sfreq"] = self._sfreq - with open(config_path, "w") as f: - yaml.dump(config_data, f, default_flow_style=False) - - try: - r_dir = _find_r_script_dir() - except FileNotFoundError as e: - logger.warning(str(e)) - self._write_python_summary() - return - - r_script = r_dir / "vertex_spatial_analysis.R" - if not r_script.exists(): - logger.warning("R script not found: %s", r_script) - self._write_python_summary() - return - - cmd = [ - "Rscript", str(r_script), - "--data-dir", str(data_dir), - "--config", str(config_path), - "--output-dir", str(self.output_dir), - "--fig-dir", str(self.fig_dir), - "--tbl-dir", str(self.tbl_dir), - ] - cmd.extend(self._r_no_figures_flags()) - - logger.info("Calling R: %s", " ".join(cmd)) - try: - result = subprocess.run( - cmd, capture_output=True, text=True, timeout=7200, # 2 hr for large multi-group studies - ) - if result.stdout: - for line in result.stdout.strip().split("\n"): - logger.info("[R] %s", line) - if result.stderr: - for line in result.stderr.strip().split("\n"): - if line.strip(): - logger.info("[R] %s", line) - if result.returncode != 0: - logger.error("R script failed (exit %d)", result.returncode) - self._write_python_summary() - except FileNotFoundError: - logger.warning("Rscript not found — writing Python summary") - self._write_python_summary() - except subprocess.TimeoutExpired: - logger.error("R script timed out after 1200 seconds") - self._write_python_summary() - - def _write_python_summary(self) -> None: lines = [ - "# Vertex Spatial Analysis Summary", + "# Vertex Spatial Analysis — RETIRED", "", f"**Study**: {self.config.name}", - "**Analysis**: Vertex Spatial Analysis", - f"**Correlation structure**: {self._correlation_structure}", - f"**Spatial range**: {self._spatial_range_mm} mm", - "", - "## Methods", - "", - "Spatial generalized least squares (nlme::gls) was used to model vertex-level " - "band power as a function of group, with an exponential spatial correlation " - "structure (`corExp(form = ~x+y+z | subject)`). This single-model approach " - "accounts for spatial autocorrelation and avoids the multiple comparison " - "problem inherent in vertex-wise testing.", "", - "## Status", + RETIRE_NOTE, "", - "R analysis not available. Data exported to `data/vertex_spatial_data.csv` " - "for manual R analysis.", - "", - ] - - if self._epoch_config is not None: - lines.insert(-2, - f"**Epoch sampling**: {self._epoch_config.get('n_epochs', 80)} epochs " - f"of {self._epoch_config.get('epoch_duration_sec', 2.0)}s" - ) - - lines.extend([ "## Output Files", "", - "- `data/vertex_spatial_data.csv` — per-subject per-vertex band power with coordinates", - "- `data/source_coords.csv` — vertex coordinates (mm)", + "- `tables/vertex_spatial_results.csv` — empty (retired)", + "- `tables/vertex_spatial_residuals.csv` — empty (retired)", "", - ]) - - summary_path = self.output_dir / "ANALYSIS_SUMMARY.md" - summary_path.write_text("\n".join(lines)) - logger.info("Wrote %s", summary_path) + ] + path = self.output_dir / "ANALYSIS_SUMMARY.md" + path.write_text("\n".join(lines)) + logger.info("Wrote %s", path) diff --git a/src/source_analytics/analyses/vertex_specparam_analysis.py b/src/source_analytics/analyses/vertex_specparam_analysis.py index caa0216..9eb030d 100644 --- a/src/source_analytics/analyses/vertex_specparam_analysis.py +++ b/src/source_analytics/analyses/vertex_specparam_analysis.py @@ -25,7 +25,7 @@ from ..spectral.aperiodic import band_peak_reachability, resolve_freq_range from ..spectral.vertex import compute_psd_vertices from ..spectral.vertex_aperiodic import fit_aperiodic_vertices -from ..spectral.epoch_sampler import sample_epochs, get_epoch_config +from ..spectral.epoch_sampler import sample_epochs from ..stats.cluster_permutation import ( cluster_permutation_test, has_significant_cluster as _has_significant_cluster, @@ -102,7 +102,8 @@ def __init__(self, config: StudyConfig, output_dir: Path): if self._noise_exclude is not None: self._noise_exclude = tuple(self._noise_exclude) - self._epoch_config = get_epoch_config(wb_cfg) + # Global epoch_sampling → vertex: block → per-analysis block (see base). + self._epoch_config = self._vertex_epoch_config() self._cluster_results: dict = {} def setup(self) -> None: diff --git a/src/source_analytics/cli.py b/src/source_analytics/cli.py index da28b7d..cf22262 100644 --- a/src/source_analytics/cli.py +++ b/src/source_analytics/cli.py @@ -10,7 +10,7 @@ import yaml from .config import StudyConfig -from .core import StudyAnalyzer, ANALYSIS_REGISTRY, ANALYSIS_METADATA +from .core import StudyAnalyzer, ANALYSIS_REGISTRY, ANALYSIS_METADATA, canonical_analysis_name from .analyses.base import VALID_STEPS, BaseAnalysis @@ -41,13 +41,14 @@ def _run_single( analysis_name: str, steps: set[str] | None = None, select: dict[str, frozenset[str]] | None = None, - jobs: int = 1, + jobs: int | None = None, ): """Run one analysis on a (possibly paradigm-scoped) config.""" analyzer = StudyAnalyzer(config) _print_study_summary(config, analyzer) analyzer.run_analysis(analysis_name, steps=steps, select=select, jobs=jobs) - print(f"\nDone. Output: {config.output_dir / analysis_name}") + # run_analysis writes under the CANONICAL name (deprecated aliases resolve). + print(f"\nDone. Output: {config.output_dir / canonical_analysis_name(analysis_name)}") def _parse_selection(args) -> dict[str, frozenset[str]] | None: @@ -110,22 +111,51 @@ def _add(dim: str, raw: str): return {d: frozenset(v) for d, v in sel.items()} -def _check_output_clean(out_dir, analysis_name, *, strict, force): - """Enforce --strict-output: error if `out_dir / analysis_name` exists. +def _prepare_output(config, analysis_name, *, strict, force, steps=None): + """Apply --strict-output / --force to the analysis's output directories. - --force overrides (re-process anyway). + All paths use the CANONICAL analysis name (a deprecated alias such as + ``psd`` writes to ``roi_psd/``), matching what ``run_analysis`` does. + + --strict-output: error if the working dir (``analytics//``) + already holds output, unless --force. + + --force: actually remove the previous output before running — the + published ``tables/`` and ``figures/`` dirs always, and the working dir + (``data/`` + summary) only when the ``process`` step will re-create it. + With ``--steps`` that excludes ``process`` the persisted data is kept, since + the requested steps read it back. """ - if not strict or force: - return - target = out_dir / analysis_name - if target.exists() and any(target.iterdir()): + import shutil + + canonical = canonical_analysis_name(analysis_name) + work = config.output_dir / canonical + paradigm = config.paradigm_name or "" + published = [ + config.results_dir / "tables" / paradigm / canonical, + config.results_dir / "figures" / paradigm / canonical, + ] + has_output = work.exists() and any(work.iterdir()) + + if strict and has_output and not force: print( - f"ERROR: --strict-output set and analysis output already exists: {target}", + f"ERROR: --strict-output set and analysis output already exists: {work}", file=sys.stderr, ) print("Pass --force to overwrite, or remove the directory first.", file=sys.stderr) sys.exit(1) + if force: + reprocess = steps is None or "process" in steps + targets = list(published) + ([work] if reprocess else []) + for t in targets: + if t.exists(): + shutil.rmtree(t) + print(f"--force: removed {t}") + if not reprocess and work.exists(): + print(f"--force: kept {work} (data/) because --steps excludes 'process'") + return work + def cmd_run(args): """Run an analysis module.""" @@ -176,7 +206,7 @@ def cmd_run(args): # Parse --metric / --band / --select sub-output selection select = _parse_selection(args) - jobs = getattr(args, "jobs", 1) or 1 + jobs = getattr(args, "jobs", None) # None → YAML `jobs:` → serial; 0/-1 → auto strict = getattr(args, "strict_output", False) force = getattr(args, "force", False) @@ -187,7 +217,7 @@ def cmd_run(args): if args.analysis: # Scope to one paradigm + one analysis aconfig = config.for_paradigm_analysis(args.paradigm, args.analysis) - _check_output_clean(aconfig.output_dir, args.analysis, strict=strict, force=force) + _prepare_output(aconfig, args.analysis, strict=strict, force=force, steps=steps) _run_single(aconfig, args.analysis, steps=steps, select=select, jobs=jobs) else: # Run all analyses listed for this paradigm @@ -200,7 +230,7 @@ def cmd_run(args): print(f"Paradigm: {args.paradigm} | Analysis: {analysis_name}") print(f"{'='*60}") aconfig = config.for_paradigm_analysis(args.paradigm, analysis_name) - _check_output_clean(aconfig.output_dir, analysis_name, strict=strict, force=force) + _prepare_output(aconfig, analysis_name, strict=strict, force=force, steps=steps) _run_single(aconfig, analysis_name, steps=steps, select=select, jobs=jobs) print() else: @@ -219,7 +249,7 @@ def cmd_run(args): print(f"Paradigm: {pname} | Analysis: {analysis_name}") print(f"{'='*60}") aconfig = config.for_paradigm_analysis(pname, analysis_name) - _check_output_clean(aconfig.output_dir, analysis_name, strict=strict, force=force) + _prepare_output(aconfig, analysis_name, strict=strict, force=force, steps=steps) _run_single(aconfig, analysis_name, steps=steps, select=select, jobs=jobs) print() else: @@ -227,11 +257,11 @@ def cmd_run(args): if not args.analysis: print("ERROR: --analysis is required for single-paradigm configs.") sys.exit(1) - _check_output_clean(config.output_dir, args.analysis, strict=strict, force=force) + _prepare_output(config, args.analysis, strict=strict, force=force, steps=steps) analyzer = StudyAnalyzer(config) _print_study_summary(config, analyzer) analyzer.run_analysis(args.analysis, steps=steps, select=select, jobs=jobs) - print(f"\nDone. Output: {config.output_dir / args.analysis}") + print(f"\nDone. Output: {config.output_dir / canonical_analysis_name(args.analysis)}") def _validate_single(config: StudyConfig, paradigm_name: str | None = None): @@ -388,6 +418,7 @@ def cmd_list(args): "resting|vertex": "Resting State (Vertex Level)", "resting|electrode": "Resting State (Electrode Level)", "evoked|roi": "Evoked Response (ROI Level)", + "evoked|vertex": "Evoked Response (Vertex Level)", "evoked|electrode": "Evoked Response (Electrode Level)", } @@ -493,103 +524,172 @@ def cmd_figure(args): print("\nNo figures generated (no data matched filters or no significant results).") +def _discover_init_subjects(derivatives_dir: Path, data_subdir: str): + """Find subjects under a reconstruction root, in either supported layout. + + Returns ``(layout, subject_groups)`` where ``layout`` is ``"flat"`` (BIDS + ``sub-*`` dirs directly under ``derivatives/``; groups unknown unless + ``--groups-from`` supplies them) or ``"grouped"`` (``//`` + dirs; the group is the folder name). + """ + flat = sorted( + d.name for d in derivatives_dir.iterdir() + if d.is_dir() and d.name.startswith("sub-") + ) + if flat: + return "flat", {name: "UNKNOWN" for name in flat} + + grouped: dict[str, str] = {} + for group_dir in sorted(d for d in derivatives_dir.iterdir() if d.is_dir()): + for subj_dir in sorted(d for d in group_dir.iterdir() if d.is_dir()): + if (subj_dir / data_subdir).is_dir() or (subj_dir / "data").is_dir(): + grouped[subj_dir.name] = group_dir.name + if grouped: + return "grouped", grouped + return "flat", {} + + +# Starter bands for a scaffolded config — the README's canonical set. +_INIT_BANDS = { + "Delta": [1, 4], + "Theta": [4, 10], + "Alpha": [10, 13], + "Beta": [13, 30], + "Low Gamma": [30, 55], + "High Gamma": [65, 80], +} + +_INIT_DEFAULT_ANALYSES = ["roi_psd", "roi_aperiodic", "roi_connectivity"] + + def cmd_init(args): - """Scaffold a source-analytics YAML config from a paradigm directory.""" + """Scaffold a study YAML (design:/hypotheses: + paradigms:) from a reconstruction dir. + + Writes ``{paradigm_dir}/analysis/{name}.yaml`` by default; ``--output -`` + prints the YAML to stdout (status goes to stderr) so it can be redirected. + The file parses with ``StudyConfig.from_yaml`` as-is; edit groups/hypotheses + and the analysis list, then ``validate`` it. + """ paradigm_dir = Path(args.paradigm_dir).resolve() derivatives_dir = paradigm_dir / "derivatives" + data_subdir = args.data_subdir + err = sys.stderr if not derivatives_dir.is_dir(): - print(f"ERROR: derivatives directory not found: {derivatives_dir}") + print(f"ERROR: derivatives directory not found: {derivatives_dir}", file=err) sys.exit(1) - # Discover sub-* directories - subject_dirs = sorted( - d.name for d in derivatives_dir.iterdir() - if d.is_dir() and d.name.startswith("sub-") - ) - if not subject_dirs: - print(f"ERROR: no sub-* directories found in {derivatives_dir}") + layout, subject_groups = _discover_init_subjects(derivatives_dir, data_subdir) + if not subject_groups: + print( + f"ERROR: no subjects found in {derivatives_dir} (expected sub-* dirs, " + "or // dirs).", file=err, + ) sys.exit(1) - # Build subject_groups mapping - subject_groups = {} + # --groups-from: source-localization study_config has subjects[].id/.group if args.groups_from: - groups_path = Path(args.groups_from).resolve() - with open(groups_path) as f: - src_config = yaml.safe_load(f) - # source-localization study_config has subjects[].id and subjects[].group + with open(Path(args.groups_from).resolve()) as f: + src_config = yaml.safe_load(f) or {} id_to_group = {} - for s in src_config.get("subjects", []): - sid = s.get("id", "") - group = s.get("group") + for s_ in src_config.get("subjects", []) or []: + sid, group = str(s_.get("id", "")), s_.get("group") if sid and group: + id_to_group[sid] = group id_to_group[f"sub-{sid}"] = group - for subj in subject_dirs: + for subj in subject_groups: if subj in id_to_group: subject_groups[subj] = id_to_group[subj] - else: - subject_groups[subj] = "UNKNOWN" - else: - for subj in subject_dirs: - subject_groups[subj] = "UNKNOWN" - - # Determine unique groups (excluding UNKNOWN) - unique_groups = sorted(set( - g for g in subject_groups.values() if g != "UNKNOWN" - )) + elif subj.startswith("sub-") and subj[4:] in id_to_group: + subject_groups[subj] = id_to_group[subj[4:]] + + unique_groups = sorted({g for g in subject_groups.values() if g != "UNKNOWN"}) + levels = unique_groups or ["Group1", "Group2"] + + # Declared hypotheses: an omnibus when there are >2 groups, plus every + # pairwise contrast (weights: later level minus earlier level). + hypotheses: list[dict] = [] + if len(levels) > 2: + hypotheses.append({"name": "group_omnibus", "kind": "omnibus", "role": "phenotype"}) + for i, ga in enumerate(levels): + for gb in levels[i + 1:]: + hypotheses.append({ + "name": f"{gb}_vs_{ga}", + "kind": "contrast", + "label": f"{gb} vs {ga}", + "weights": {gb: 1, ga: -1}, + "role": "phenotype", + }) - # Build YAML structure config_name = args.name or paradigm_dir.name + out_dir = paradigm_dir / "analysis" + analyses = [a.strip() for a in args.analyses.split(",") if a.strip()] + unknown = [a for a in analyses if a not in ANALYSIS_REGISTRY] + if unknown: + print(f"ERROR: unknown analyses: {', '.join(unknown)} (see `source-analytics list`)", file=err) + sys.exit(1) + + paradigm_block: dict = { + "data_dir": str(derivatives_dir), + "data_subdir": data_subdir, + "analyses": {a: {} for a in analyses}, + } + if layout == "flat": + # Flat sub-* layout needs the explicit subject → group map. + paradigm_block["subjects"] = dict(subject_groups) + config = { "name": config_name, - "groups": {g: g for g in unique_groups} if unique_groups else {"Group1": "Group 1"}, - "group_order": unique_groups if unique_groups else ["Group1"], + "groups": {g: g for g in levels}, + "group_order": list(levels), "group_colors": {}, - "contrasts": [], - "bands": { - "delta": [2, 3.5], - "theta": [3.5, 7.5], - "alpha_1": [8, 10], - "alpha_2": [10.5, 12.5], - "beta": [13, 30], - "gamma_1": [30, 55], - "gamma_2": [65, 80], - "epsilon": [81, 120], - }, - "roi_categories": {}, - "discovery": { - "data_subdir": "pipeline/data", - "subject_groups": subject_groups, + "design": {"factor": "group", "reference": levels[0], "levels": list(levels)}, + "hypotheses": hypotheses, + "bands": dict(_INIT_BANDS), + "epoch_sampling": { + "enabled": False, "epoch_duration_sec": 2.0, "n_epochs": 80, "n_bootstrap": 1, }, + # Absolute so the file works from any working directory. + "output_dir": str(out_dir / "analytics"), + "results_dir": str(out_dir / "results"), + "paradigms": {args.paradigm: paradigm_block}, } - # Generate pairwise contrasts - if len(unique_groups) >= 2: - for i, ga in enumerate(unique_groups): - for gb in unique_groups[i + 1:]: - config["contrasts"].append({ - "name": f"{ga}_vs_{gb}", - "group_a": ga, - "group_b": gb, - }) - - # Write config - analysis_dir = paradigm_dir / "analysis" - analysis_dir.mkdir(exist_ok=True) - out_path = analysis_dir / f"{config_name}.yaml" - - with open(out_path, "w") as f: - yaml.dump(config, f, default_flow_style=False, sort_keys=False) - - print(f"Config written: {out_path}") - print(f"Subjects: {len(subject_dirs)}") - if unique_groups: - for g in unique_groups: - n = sum(1 for v in subject_groups.values() if v == g) - print(f" {g}: n={n}") + header = ( + f"# source-analytics study config, scaffolded by `source-analytics init` from\n" + f"# {paradigm_dir}\n" + f"# Layout: {layout}. Edit groups / hypotheses / analyses, then:\n" + f"# source-analytics validate --study \n" + f"# source-analytics run --study --paradigm {args.paradigm} " + f"--analysis {analyses[0] if analyses else 'roi_psd'}\n" + ) + text = header + yaml.dump(config, default_flow_style=False, sort_keys=False, allow_unicode=True) + + if args.output == "-": + sys.stdout.write(text) + out_path = None + else: + out_path = Path(args.output).resolve() if args.output else out_dir / f"{config_name}.yaml" + out_path.parent.mkdir(parents=True, exist_ok=True) + out_path.write_text(text) + print(f"Config written: {out_path}", file=err) + + print(f"Layout: {layout}; subjects: {len(subject_groups)}", file=err) + for g in unique_groups: + n = sum(1 for v in subject_groups.values() if v == g) + print(f" {g}: n={n}", file=err) n_unknown = sum(1 for v in subject_groups.values() if v == "UNKNOWN") if n_unknown: - print(f" UNKNOWN: n={n_unknown} (edit config to assign groups)") + print( + f" UNKNOWN: n={n_unknown} (pass --groups-from, or edit " + f"paradigms.{args.paradigm}.subjects to assign groups)", file=err, + ) + if out_path is not None: + # Prove the scaffold parses before the user edits it. + try: + StudyConfig.from_yaml(out_path) + except Exception as exc: # noqa: BLE001 + print(f"WARNING: scaffolded config failed to parse: {exc}", file=err) def main(): @@ -619,11 +719,13 @@ def main(): f"Valid: {', '.join(sorted(VALID_STEPS))}", ) p_run.add_argument( - "--jobs", "-j", type=int, default=1, metavar="N", - help="Parallel worker processes for the per-subject process step " - "(default 1 = serial). N<=0 uses all-but-one core. Only parallel-capable " - "modules (the vertex analyses) use it; others run serially. Results are " - "identical to serial regardless of N.", + "--jobs", "-j", type=int, default=None, metavar="N", + help="Parallel worker processes for the per-subject process step. " + "An explicit N wins over the study YAML's top-level `jobs:`; when omitted " + "that YAML value is used, else serial. N<=0 uses all-but-one core. Only " + "parallel-capable modules use it (the vertex analyses, roi_connectivity, " + "electrode_connectivity); others run serially. Results are identical to " + "serial regardless of N.", ) p_run.add_argument( "--metric", @@ -662,7 +764,10 @@ def main(): p_run.add_argument( "--force", action="store_true", - help="Overwrite the analysis output directory if it already exists", + help="Remove the analysis's previous output before running: its published " + "tables/ and figures/ dirs always, and its analytics working dir " + "(data/ + summary) when the process step runs. Also overrides " + "--strict-output.", ) p_run.set_defaults(func=cmd_run) @@ -690,10 +795,24 @@ def main(): p_fig.set_defaults(func=cmd_figure) # init - p_init = subparsers.add_parser("init", help="Scaffold analysis config from paradigm directory") - p_init.add_argument("paradigm_dir", type=Path, help="Paradigm directory (contains derivatives/)") + p_init = subparsers.add_parser( + "init", + help="Scaffold a study YAML (design/hypotheses + one paradigm) from a " + "reconstruction directory; writes

/analysis/.yaml", + ) + p_init.add_argument("paradigm_dir", type=Path, help="Reconstruction directory (contains derivatives/)") p_init.add_argument("--name", help="Study name (default: directory name)") p_init.add_argument("--groups-from", type=Path, help="source-localization study_config.yaml for group mappings") + p_init.add_argument("--paradigm", default="resting", help="Name of the paradigm block to emit (default: resting)") + p_init.add_argument( + "--analyses", default=",".join(_INIT_DEFAULT_ANALYSES), metavar="a,b,...", + help=f"Analyses to list under the paradigm (default: {','.join(_INIT_DEFAULT_ANALYSES)})", + ) + p_init.add_argument("--data-subdir", default="pipeline/data", help="Per-subject data subdir (default: pipeline/data)") + p_init.add_argument( + "--output", "-o", metavar="PATH", + help="Where to write the YAML (default: /analysis/.yaml); '-' prints to stdout", + ) p_init.set_defaults(func=cmd_init) args = parser.parse_args() diff --git a/src/source_analytics/core.py b/src/source_analytics/core.py index 337495c..4d2b0cb 100644 --- a/src/source_analytics/core.py +++ b/src/source_analytics/core.py @@ -96,6 +96,11 @@ ANALYSIS_REGISTRY[_old] = ANALYSIS_REGISTRY[_new] +def canonical_analysis_name(name: str) -> str: + """Canonical registry name for ``name`` (alias-resolved), without warning.""" + return _DEPRECATED_NAMES.get(name, name) + + def resolve_analysis_name(name: str) -> str: """Resolve a possibly-deprecated analysis name to the canonical name. @@ -116,11 +121,15 @@ def resolve_analysis_name(name: str) -> str: # domain = how analyses are grouped/listed (by the data they use) # supplements = a SECONDARY analysis that can only run after the named primary # (it consumes the primary's output). Absent = primary. -ANALYSIS_METADATA: dict[str, dict[str, str]] = { +# requires = every upstream module whose output the analysis reads (a +# superset of `supplements` for modules with >1 primary). The +# lightbox groups by `supplements`; `requires` is the honest +# dependency list for run ordering. +ANALYSIS_METADATA: dict[str, dict] = { "roi_psd": {"category": "resting", "level": "roi", "domain": "Spectral", "description": "PSD (power spectral density)", - "about": "Resting-state oscillatory power in each source-localized ROI. Per subject, ROI, and band, power is integrated from the ROI power spectrum and reported two ways: absolute power in decibels (10*log10 of the band's integrated power) and relative power (the band's fraction of total power across the analyzed spectrum). Groups are compared per ROI and band with a linear mixed model (subject as a random effect, via lmerTest) and contrasts estimated with emmeans; effect sizes are Hedges' g, with p-values FDR-corrected (Benjamini-Hochberg) within each band. Read it as: which ROIs differ in oscillatory power, in which bands, and in which direction (an up arrow means the first-listed group of the pair is higher). Relative power controls for overall spectral amplitude and isolates spectral shape; absolute (dB) is the band's log power."}, + "about": "Resting-state oscillatory power in each source-localized ROI. Per subject, ROI, and band, power is integrated from the ROI power spectrum and reported two ways: absolute power as the mean power density over the band in dB/Hz (10*log10 of the band's integrated power divided by its bandwidth, so the across-band 1/f shape is visible rather than width-confounded) and relative power (the band's fraction of total power across the analyzed spectrum). Groups are compared per ROI and band with a linear mixed model (subject as a random effect, via lmerTest) and contrasts estimated with emmeans; effect sizes are Hedges' g, with p-values FDR-corrected (Benjamini-Hochberg) within each band. Read it as: which ROIs differ in oscillatory power, in which bands, and in which direction (an up arrow means the first-listed group of the pair is higher). Relative power controls for overall spectral amplitude and isolates spectral shape; absolute (dB/Hz) is the band's log power density."}, "roi_aperiodic": {"category": "resting", "level": "roi", "domain": "Spectral", "description": "1/f aperiodic decomposition", - "about": "The aperiodic (1/f) component of each ROI's resting power spectrum, fit with specparam/FOOOF over 2-50 Hz: the exponent (steepness of the 1/f decay, a proxy for excitation/inhibition balance -- steeper usually read as more inhibition) and the offset (broadband power level). Groups are compared per ROI with a linear mixed model (dv ~ group * ROI, subject as a random effect; Type-III ANOVA, Satterthwaite df) followed by emmeans pairwise contrasts gated on the group-by-ROI omnibus; effect sizes are Hedges' g and per-ROI post-hoc p-values are Holm-corrected. There is no frequency-band axis -- exponent and offset are single broadband parameters. Read it as: which ROIs show a group shift in spectral slope (E/I tilt) or broadband power, and in which direction (an up arrow means the first-listed group is higher)."}, + "about": "The aperiodic (1/f) component of each ROI's resting power spectrum, fit with specparam/FOOOF over 12-45 Hz by default (spectral.aperiodic.DEFAULT_FREQ_RANGE; see docs/methods/APERIODIC_FIT_WINDOW.md): the exponent (steepness of the 1/f decay, a proxy for excitation/inhibition balance -- steeper usually read as more inhibition) and the offset (broadband power level). Groups are compared per ROI with a linear mixed model (dv ~ group * ROI, subject as a random effect; Type-III ANOVA, Satterthwaite df) followed by emmeans pairwise contrasts gated on the group-by-ROI omnibus; effect sizes are Hedges' g and per-ROI post-hoc p-values are Holm-corrected. There is no frequency-band axis -- exponent and offset are single broadband parameters. Read it as: which ROIs show a group shift in spectral slope (E/I tilt) or broadband power, and in which direction (an up arrow means the first-listed group is higher)."}, "roi_connectivity": {"category": "resting", "level": "roi", "domain": "Connectivity", "description": "ROI pairwise connectivity matrices (descriptive; inference in roi_nbs / roi_graph)", "about": "Resting functional connectivity between every pair of source-localized ROIs, per band. For each pair it estimates one or more coupling measures from the cross-spectrum: coherence, imaginary coherence (Nolte et al. 2004), the phase-lag index and its weighted/debiased forms (PLI Stam et al. 2007; wPLI and dwPLI Vinck et al. 2011), the directed phase-lag index (dPLI Stam & van Straaten 2012), amplitude-envelope correlation with leakage orthogonalization (AEC Hipp et al. 2012), and shrinkage partial correlation (Ledoit-Wolf). Corpus-callosum tracts are excluded, leaving the cortical ROI pairs. This module is descriptive: it builds the group-mean ROI-by-ROI connectivity matrix per band and metric and renders it (circos + heatmap, with the between-group difference), but carries no per-edge inferential test -- the legacy per-pair Welch t-test and region-pair linear mixed model were retired because edge-by-edge testing is underpowered at this ROI count. Group-level inference for the connectivity family is provided by roi_nbs (the Network-Based Statistic, for connected sub-networks that differ between groups) and roi_graph (graph-theoretic network organization). Read it as: the connectivity structure per band and its group differences (a warm difference / up arrow means the first-listed group is higher); imaginary-coherence/PLI-family metrics suppress volume-conduction/zero-lag artifacts."}, "roi_cross_freq": {"category": "resting", "level": "roi", "domain": "Cross-frequency", "description": "Cross-frequency coupling (PAC, AAC, n:m PPC)", @@ -142,8 +151,8 @@ def resolve_analysis_name(name: str) -> str: "electrode_signature": {"category": "resting", "level": "electrode", "domain": "Source vs Sensor", "display_name": "Neural signature", "description": "Sensor-level neural signature (classification/decoding on electrode band power) — the source-vs-sensor counterpart of vertex_signature"}, "vertex_signature": {"category": "resting", "level": "vertex", "domain": "Multivariate", "display_name": "Neural signature", "description": "Multivariate/ML neural signature (classification, decoding; PCA-reduced with back-projection)"}, "vertex_cluster": {"category": "resting", "level": "vertex", "domain": "Spectral", "description": "Vertex-level cluster permutation", - "about": "Whole-brain resting spectral maps on the dorsal source surface: per vertex it computes band power (absolute in dB and relative), the 1/f spectral slope, and the peak alpha frequency, then tests where the groups differ. Inference is a cluster-based permutation test -- per-vertex t-statistics are threshold-clustered over neighbouring vertices and each cluster's extent is compared to a permutation null, giving family-wise (FWE) control (Maris & Oostenveld 2007); a threshold-free TFCE variant (Smith & Nichols 2009) is available. Effect sizes are per-vertex Hedges' g. Read it as: spatially-contiguous clusters where the groups differ in a spectral measure -- a cluster with p_corrected < 0.05 marks a region of difference, and the sign of its t-values gives the direction. This is the whole-brain, unparcellated counterpart to the ROI spectral analyses."}, - "vertex_spatial": {"category": "resting", "level": "vertex", "domain": "Spectral", "description": "Spatial GLS (vertex-level generalized least squares)"}, + "about": "Whole-brain resting spectral maps on the dorsal source surface: per vertex it computes band power (absolute as mean density in dB/Hz, the same definition as roi_psd, and relative), the 1/f spectral slope, and the peak alpha frequency, then tests where the groups differ. Inference is a cluster-based permutation test -- per-vertex t-statistics are threshold-clustered over neighbouring vertices and each cluster's extent is compared to a permutation null, giving family-wise (FWE) control (Maris & Oostenveld 2007); a threshold-free TFCE variant (Smith & Nichols 2009) is available. Effect sizes are per-vertex Hedges' g. Read it as: spatially-contiguous clusters where the groups differ in a spectral measure -- a cluster with p_corrected < 0.05 marks a region of difference, and the sign of its t-values gives the direction. This is the whole-brain, unparcellated counterpart to the ROI spectral analyses."}, + "vertex_spatial": {"category": "resting", "level": "vertex", "domain": "Spectral", "description": "RETIRED — was: spatial GLS robustness check; exits with empty tables (use vertex_cluster / vertex_nbs)"}, "vertex_specparam": {"category": "resting", "level": "vertex", "domain": "Spectral", "description": "Vertex-level spectral parameterization", "about": "The aperiodic (1/f) spectrum fit per vertex across the dorsal source surface with specparam/FOOOF -- exponent, offset, and per-band oscillatory peaks (presence, frequency, power). Group differences in the exponent and offset maps, and in per-band peak power, are tested with a cluster-based permutation test (threshold-clustered vertex t-statistics with cluster-extent FWE correction by permutation); band peak presence is compared with a per-vertex chi-square test. Effect sizes are per-vertex Hedges' g. Read it as: spatially-contiguous clusters of vertices where the groups differ in spectral slope, broadband power, or an oscillatory peak -- a cluster with p_corrected < 0.05 marks a region of difference, the sign of its t-values gives direction. This is the whole-brain, unparcellated version of roi_aperiodic (plus peaks)."}, "vertex_connectivity": {"category": "resting", "level": "vertex", "domain": "Connectivity", "description": "Vertex pairwise connectivity", @@ -151,14 +160,14 @@ def resolve_analysis_name(name: str) -> str: "vertex_cross_freq": {"category": "resting", "level": "vertex", "domain": "Cross-frequency", "description": "Vertex cross-frequency coupling (local PAC, AAC, n:m PPC)", "about": "Whole-brain cross-frequency coupling on the dorsal source surface -- the unparcellated counterpart of roi_cross_freq, using the same kernels for each valid slow-phase x fast-amplitude band pair. Phase-amplitude coupling (PAC, Tort et al. 2010 modulation index, surrogate z-scored) is computed locally -- the slow phase and fast amplitude come from the same vertex -- yielding a whole-brain coupling map. Amplitude-amplitude coupling (AAC, power-envelope correlation) and n:m phase-phase coupling (PPC, Palva et al. 2005 phase-locking factor) are computed all-to-all across vertices and summarized to a per-vertex node strength (mean off-diagonal coupling). Group differences in these maps are tested with a cluster-based permutation test (per-vertex t-statistics clustered by spatial adjacency, cluster-extent FWE from a permutation null; Maris & Oostenveld 2007), with per-vertex Hedges' g. Read it as: spatially-contiguous clusters where the groups differ in cross-frequency coupling -- a cluster with p_corrected < 0.05 marks a region of difference, the sign of its t-values gives direction. PAC here is the primary source-spatial-advantage measure (local, no leakage between nodes)."}, "electrode_psd": {"category": "resting", "level": "electrode", "domain": "Sensor-level", "description": "Sensor-level PSD analysis", - "about": "Resting band power at each scalp electrode -- the sensor-space counterpart of roi_psd. Per channel and band, power is reported as absolute power in decibels (10*log10 of the band's integrated power) and relative power (the band's fraction of total 1-100 Hz power). Groups are compared per channel and band with a linear mixed model (dv ~ group * channel, subject as a random effect), with an optional region-nested model over the configured electrode regions (channels as replicates); per-contrast effects come from the hypothesis layer as Hedges' g with band-wise Benjamini-Hochberg FDR. Read it as: which electrodes/regions differ in band power, in which bands and direction (up = the first-listed group is higher) -- the sensor-level check against the source (ROI) result."}, + "about": "Resting band power at each scalp electrode -- the sensor-space counterpart of roi_psd. Per channel and band, power is reported as absolute power density in dB/Hz (10*log10 of the band's integrated power divided by its bandwidth, as in roi_psd) and relative power (the band's fraction of total 1-100 Hz power). Groups are compared per channel and band with a linear mixed model (dv ~ group * channel, subject as a random effect), with an optional region-nested model over the configured electrode regions (channels as replicates); per-contrast effects come from the hypothesis layer as Hedges' g with band-wise Benjamini-Hochberg FDR. Read it as: which electrodes/regions differ in band power, in which bands and direction (up = the first-listed group is higher) -- the sensor-level check against the source (ROI) result."}, "electrode_aperiodic": {"category": "resting", "level": "electrode", "domain": "Sensor-level", "description": "Sensor-level aperiodic (1/f) analysis", - "about": "The aperiodic (1/f) exponent and offset at each scalp electrode, fit with specparam/FOOOF over 2-50 Hz -- the sensor-space counterpart of roi_aperiodic. Groups are compared per channel with a linear mixed model (dv ~ group * channel, subject as a random effect) plus a region-nested model (dv ~ group * region with (1|subject/channel), treating channels as replicates within scalp regions); effect sizes are Hedges' g and post-hoc p-values are Benjamini-Hochberg FDR-corrected across channels. Read it as: which electrodes/regions show a group shift in spectral slope (E/I proxy) or broadband offset, and in which direction (up = the first-listed group is higher)."}, - "electrode_comparison": {"category": "resting", "level": "electrode", "domain": "Source vs Sensor", "display_name": "PSD", "description": "Source vs electrode comparison", + "about": "The aperiodic (1/f) exponent and offset at each scalp electrode, fit with specparam/FOOOF over 12-45 Hz by default (spectral.aperiodic.DEFAULT_FREQ_RANGE) -- the sensor-space counterpart of roi_aperiodic. Groups are compared per channel with a linear mixed model (dv ~ group * channel, subject as a random effect) plus a region-nested model (dv ~ group * region with (1|subject/channel), treating channels as replicates within scalp regions); effect sizes are Hedges' g and post-hoc p-values are Benjamini-Hochberg FDR-corrected across channels. Read it as: which electrodes/regions show a group shift in spectral slope (E/I proxy) or broadband offset, and in which direction (up = the first-listed group is higher)."}, + "electrode_comparison": {"category": "resting", "level": "electrode", "domain": "Source vs Sensor", "display_name": "PSD", "supplements": "electrode_psd", "requires": ["electrode_psd", "roi_psd"], "description": "Source vs electrode comparison (needs electrode_psd AND roi_psd)", "about": "A source-versus-sensor check on resting band power: for each subject and band, electrode power (averaged over channels) is compared to source power (averaged over ROIs). It reports (1) the cross-subject concordance between sensor and source power (Pearson r per band) and (2) whether the group effect agrees at both levels -- per contrast, Hedges' g with 95% CIs at the electrode level and the source (ROI/region) level, plus an 'exceeds_electrode' flag where a region's effect is larger than the global sensor effect. There is no cluster/FWE correction here; significance is read from whether the 95% CI excludes zero. Read it as: does the source reconstruction recover the same spectral group effect the scalp shows, and does it localize it more sharply than the sensor average?"}, "electrode_connectivity": {"category": "resting", "level": "electrode", "domain": "Sensor-level", "description": "Sensor pairwise connectivity + FCD (source-vs-sensor comparator)", "about": "Resting connectivity between the 30 scalp electrodes -- the sensor-space comparator for vertex_connectivity. Per band it computes all-to-all channel coupling with the leakage/volume-conduction-robust subset of the same kernels (AEC, imaginary coherence, PLI, wPLI, dwPLI, dPLI) and the per-channel functional connectivity density (FCD; degree/(n-1) above threshold, Tomasi & Volkow 2010). Groups are compared per channel with a Welch t-test and Benjamini-Hochberg FDR across the 30 channels (effect sizes Hedges' g), plus a hypothesis-layer cluster-permutation test over the sensor montage (adjacency from channel coordinates). Read it as: which electrodes differ in connectivity/FCD, in which band and direction (up = the first-listed group is higher) -- the scalp-level check on whether the source FCD effect is also visible without source reconstruction. Because sensor space is blurred by volume conduction, the volume-conduction-sensitive metrics are omitted here."}, - "fcd_comparison": {"category": "resting", "level": "electrode", "domain": "Source vs Sensor", "display_name": "Connectivity", "description": "Source vs sensor FCD comparison (mean + spatial CV)", + "fcd_comparison": {"category": "resting", "level": "electrode", "domain": "Source vs Sensor", "display_name": "Connectivity", "supplements": "electrode_connectivity", "requires": ["electrode_connectivity", "vertex_connectivity"], "description": "Source vs sensor FCD comparison (mean + spatial CV; needs electrode_connectivity AND vertex_connectivity, cross-paradigm)", "about": "A source-versus-sensor check on functional connectivity density (FCD), pairing the whole-brain vertex FCD maps (vertex_connectivity) against the scalp channel FCD (electrode_connectivity). FCD is each node's fraction of supra-threshold connections (degree/(n-1), threshold 0.05; for dPLI the deviation from its 0.5 no-lag center; Tomasi & Volkow 2010). Per subject and band it summarizes each map two ways -- mean FCD (overall coupling density) and spatial coefficient of variation (CV = SD/mean, how heterogeneous the map is) -- then reports (1) cross-subject concordance between source and sensor summaries (Pearson r per band) and (2) whether the group effect agrees at both levels: per contrast, Hedges' g with 95% CIs at each level plus a sign-concordance flag. There is no cluster/FWE correction here; significance is read from whether a 95% CI excludes zero. Read it as: does source-space recover the same connectivity-density group effect the scalp shows, and is the spatial pattern preserved?"}, "roi_evoked": {"category": "evoked", "level": "roi", "domain": "Evoked", "description": "ITC, ERSP, STP for trial-based paradigms"}, "vertex_evoked": {"category": "evoked", "level": "vertex", "domain": "Evoked", "description": "Vertex-level ITC, ERSP, STP (cluster-corrected) for trial-based paradigms"}, @@ -231,7 +240,7 @@ def run_analysis( analysis_name: str, steps: set[str] | None = None, select: dict[str, frozenset[str]] | None = None, - jobs: int = 1, + jobs: int | None = None, ) -> None: """Run a single named analysis. diff --git a/src/source_analytics/hypothesis/tabular.py b/src/source_analytics/hypothesis/tabular.py index e11a9e3..de61097 100644 --- a/src/source_analytics/hypothesis/tabular.py +++ b/src/source_analytics/hypothesis/tabular.py @@ -363,10 +363,22 @@ def write_module_hypotheses_tabular( if fam_rows: method, scope = _resolve_fdr(hyp, spec) facet_map = dict(zip(facet_cols, fkey)) + # The facet combination is what makes this family distinct from + # the next, so it is the `dv` coordinate of the label: it plays + # exactly the role R's `dv_col` plays, where each DV gets its own + # run_hypothesis() call and so its own family. Without it, two + # facets over the same band x spatial grid produce BYTE-IDENTICAL + # labels for genuinely different families -- the member hash + # cannot separate them, because the member cells really are the + # same. Label-only: q-values are corrected per facet either way. _apply_fdr( fam_rows, method, scope, hypothesis=hyp.name, - dv=str(facet_map.get("dv", "NA")), + dv=( + "|".join(f"{c}={facet_map[c]}" for c in facet_cols) + if facet_cols + else str(facet_map.get("dv", "NA")) + ), spatial_name=spatial_col or "spatial", ) all_rows.extend(fam_rows) diff --git a/src/source_analytics/spectral/epoch_sampler.py b/src/source_analytics/spectral/epoch_sampler.py index 657988a..c953ed0 100644 --- a/src/source_analytics/spectral/epoch_sampler.py +++ b/src/source_analytics/spectral/epoch_sampler.py @@ -57,7 +57,10 @@ def sample_epochs( Sampled epochs. ``effective_n_epochs`` equals ``n_epochs`` when ``n_bootstrap == 1``, or up to ``n_bootstrap * n_epochs`` otherwise. """ - if n_epochs is None or n_epochs == 0: + # n_epochs=0 or n_bootstrap<=0 both mean "no sampling": return the full + # timeseries as a single epoch. The ROI path (sample_roi_epochs) has + # honoured n_bootstrap=0 this way all along; the vertex path now matches. + if n_epochs is None or n_epochs == 0 or (n_bootstrap is not None and n_bootstrap <= 0): return data[np.newaxis, :, :] # (1, n_channels, n_times) if n_bootstrap > 1: diff --git a/src/source_analytics/spectral/tfr.py b/src/source_analytics/spectral/tfr.py index 16cac25..a938098 100644 --- a/src/source_analytics/spectral/tfr.py +++ b/src/source_analytics/spectral/tfr.py @@ -10,7 +10,24 @@ import logging import numpy as np -from mne.time_frequency import tfr_array_morlet + + +def tfr_array_morlet(*args, **kwargs): + """Lazy proxy for :func:`mne.time_frequency.tfr_array_morlet`. + + MNE is an optional extra (``pip install source-analytics[mne]``); importing + it here at module level would make every spectral import — and therefore + the whole CLI — depend on it. Resolve it on first use instead, with a clear + message when it is missing. + """ + try: + from mne.time_frequency import tfr_array_morlet as _impl + except ImportError as exc: # pragma: no cover - exercised only without mne + raise ImportError( + "The evoked (TFR) analyses need MNE-Python. Install it with " + "'pip install \"source-analytics[mne]\"' (or [all])." + ) from exc + return _impl(*args, **kwargs) logger = logging.getLogger(__name__) diff --git a/src/source_analytics/spectral/vertex.py b/src/source_analytics/spectral/vertex.py index d905fbb..56fd5da 100644 --- a/src/source_analytics/spectral/vertex.py +++ b/src/source_analytics/spectral/vertex.py @@ -95,7 +95,8 @@ def extract_band_power_vertices( ------- dict[str, dict[str, ndarray]] band_name -> {"absolute": (n_vertices,), "relative": (n_vertices,)} - where absolute is 10*log10(integrated power) in dB. + where absolute is the mean band power density, 10*log10(integrated + power / bandwidth), in dB/Hz — the same definition as the ROI path. """ # Total power (optionally excluding noise band) if noise_exclude is not None: @@ -129,7 +130,15 @@ def extract_band_power_vertices( abs_power = trapezoid(band_psd, band_freqs, axis=-1) rel_power = abs_power / total_power - db_power = np.where(abs_power > 0, 10.0 * np.log10(abs_power), -np.inf) + + # "absolute" is the mean power DENSITY over the band (integrated power + # / bandwidth) in dB/Hz — the same definition as the ROI/electrode path + # in band_power.py, so the column is comparable across levels. (A + # per-band constant shift relative to the old integrated-power form; + # within-band group statistics are unaffected.) + bandwidth = fmax - fmin + density = abs_power / bandwidth if bandwidth > 0 else abs_power + db_power = np.where(density > 0, 10.0 * np.log10(density), -np.inf) result[band_name] = { "absolute": db_power, diff --git a/src/source_analytics/stats/graph_metrics.py b/src/source_analytics/stats/graph_metrics.py index e04890b..2e3e3a9 100644 --- a/src/source_analytics/stats/graph_metrics.py +++ b/src/source_analytics/stats/graph_metrics.py @@ -611,8 +611,20 @@ def nbs_permutation_test( for size in observed_components ] - # Build significant edge mask + # Build significant edge mask. + # + # A component IS a connected component of the suprathreshold graph, and the + # components partition the nodes, so a component's edges are exactly the + # suprathreshold edges whose endpoints are both in its node set. Without + # this loop the mask was allocated and returned unfilled — always all-False, + # even for a 99-edge component at p = 0.0. sig_edges = np.zeros((n, n), dtype=bool) + for pval, nodes in zip(component_pvalues, observed_nodes): + if pval < 0.05 and len(nodes) >= 2: + idx = np.asarray(nodes, dtype=int) + within = np.zeros((n, n), dtype=bool) + within[np.ix_(idx, idx)] = True + sig_edges |= suprathresh & within n_sig = sum(1 for p in component_pvalues if p < 0.05) logger.info("NBS: %d/%d components significant (p<0.05)", n_sig, len(observed_components)) diff --git a/tests/test_audit_fixes.py b/tests/test_audit_fixes.py new file mode 100644 index 0000000..4abe400 --- /dev/null +++ b/tests/test_audit_fixes.py @@ -0,0 +1,291 @@ +"""Regression tests for the 2026-09 audit fixes (see CHANGELOG, Unreleased).""" + +from __future__ import annotations + +import os +import subprocess +import sys +from pathlib import Path + +import numpy as np +import pytest +import yaml + +from source_analytics.config import StudyConfig +from source_analytics.core import canonical_analysis_name, ANALYSIS_METADATA +from source_analytics.spectral.band_power import extract_band_power +from source_analytics.spectral.vertex import extract_band_power_vertices +from source_analytics.spectral.epoch_sampler import sample_epochs + + +# ---- #6: the package must import without the optional extras --------------- +def test_package_imports_without_mne(): + code = ( + "import sys; sys.modules['mne'] = None; sys.modules['mne.time_frequency'] = None\n" + "import source_analytics.core, source_analytics.cli, source_analytics.spectral\n" + "print('ok')" + ) + out = subprocess.run([sys.executable, "-c", code], capture_output=True, text=True) + assert out.returncode == 0, out.stderr + assert "ok" in out.stdout + + +def test_tfr_raises_clear_error_without_mne(monkeypatch): + import source_analytics.spectral.tfr as tfr + monkeypatch.setitem(sys.modules, "mne", None) + monkeypatch.setitem(sys.modules, "mne.time_frequency", None) + with pytest.raises(ImportError, match=r"source-analytics\[mne\]"): + tfr.tfr_array_morlet(np.zeros((1, 1, 10)), sfreq=100.0, freqs=[5.0], n_cycles=2) + + +# ---- #9: ROI and vertex `absolute` share one definition (dB/Hz density) ----- +def test_vertex_absolute_matches_roi_density(): + freqs = np.linspace(1, 100, 397) + psd_1d = 1e-12 / freqs # 1/f + bands = {"Alpha": (8, 13), "Beta": (13, 30)} + roi = extract_band_power(freqs, psd_1d, bands) + vtx = extract_band_power_vertices(freqs, psd_1d[np.newaxis, :], bands, noise_exclude=None) + for b in bands: + assert vtx[b]["absolute"][0] == pytest.approx(roi[b]["absolute"], rel=1e-9) + assert vtx[b]["relative"][0] == pytest.approx(roi[b]["relative"], rel=1e-9) + + +# ---- #12: n_bootstrap: 0 means "full timeseries" on the vertex sampler too -- +def test_sample_epochs_n_bootstrap_zero_returns_full_data(): + data = np.random.default_rng(0).standard_normal((3, 5000)) + out = sample_epochs(data, 500.0, epoch_duration_sec=1.0, n_epochs=4, seed=1, n_bootstrap=0) + assert out.shape == (1, 3, 5000) + np.testing.assert_array_equal(out[0], data) + sampled = sample_epochs(data, 500.0, epoch_duration_sec=1.0, n_epochs=4, seed=1, n_bootstrap=1) + assert sampled.shape == (4, 3, 500) + + +# ---- #19: vertex modules see global + per-analysis epoch_sampling ---------- +def test_vertex_epoch_config_merges_global_vertex_and_analysis(): + from source_analytics.analyses.vertex_cluster_analysis import VertexClusterAnalysis + + class _Cfg: + raw = { + "epoch_sampling": {"enabled": True, "n_epochs": 40, "n_bootstrap": 5}, + "vertex_cluster": {"epoch_sampling": {"n_bootstrap": 0}}, + } + vertex = {"epoch_sampling": {"n_epochs": 60}} + + a = VertexClusterAnalysis.__new__(VertexClusterAnalysis) + a.config = _Cfg() + merged = a._vertex_epoch_config() + assert merged == {"enabled": True, "n_epochs": 60, "n_bootstrap": 0} + + class _Off: + raw = {"epoch_sampling": {"n_epochs": 40}} + vertex = {} + + a.config = _Off() + assert a._vertex_epoch_config() is None + + +# ---- #17: deprecated names resolve to the canonical output dir ------------- +def test_canonical_analysis_name(): + assert canonical_analysis_name("psd") == "roi_psd" + assert canonical_analysis_name("vertex_mvpa") == "vertex_signature" + assert canonical_analysis_name("roi_psd") == "roi_psd" + + +# ---- #4/#5: dependency metadata is honest -------------------------------- +def test_comparison_modules_declare_requires(): + assert ANALYSIS_METADATA["fcd_comparison"]["requires"] == [ + "electrode_connectivity", "vertex_connectivity"] + assert ANALYSIS_METADATA["fcd_comparison"]["supplements"] == "electrode_connectivity" + assert ANALYSIS_METADATA["electrode_comparison"]["requires"] == ["electrode_psd", "roi_psd"] + + +# ---- #4: fcd_comparison finds its primaries across paradigm dirs ---------- +def test_fcd_comparison_cross_paradigm_lookup(tmp_path): + from source_analytics.analyses.fcd_comparison_analysis import FCDComparisonAnalysis + + analytics = tmp_path / "analytics" + (analytics / "resting" / "electrode_connectivity" / "data").mkdir(parents=True) + (analytics / "vertex" / "vertex_connectivity" / "data").mkdir(parents=True) + rows = "subject,group,band,metric,fcd\ns1,A,Alpha,pli,0.5\ns2,B,Alpha,pli,0.4\n" + (analytics / "resting" / "electrode_connectivity" / "data" / "electrode_fcd.csv").write_text(rows) + (analytics / "vertex" / "vertex_connectivity" / "data" / "vertex_fcd.csv").write_text(rows) + + class _Cfg: + raw = {} + output_dir = analytics / "resting" + paradigm_name = "resting" + results_dir = tmp_path / "results" + vertex = {} + rois = None + roi_categories = {} + atlas_dir = None + + a = FCDComparisonAnalysis.__new__(FCDComparisonAnalysis) + a.config = _Cfg() + a._sensor_df = a._source_df = None + a._selection = {} + src = a._find_upstream_csv("vertex_connectivity", "vertex_fcd.csv", "source_dir") + assert src == analytics / "vertex" / "vertex_connectivity" / "data" / "vertex_fcd.csv" + sen = a._find_upstream_csv("electrode_connectivity", "electrode_fcd.csv", "sensor_dir") + assert sen.parent.parent.parent.name == "resting" + with pytest.raises(FileNotFoundError, match="run 'vertex_connectivity' first"): + a._find_upstream_csv("vertex_connectivity", "nope.csv", "source_dir") + + +# ---- #1: init writes a parseable design/hypotheses/paradigms config -------- +def _fake_reconstruction(root: Path, flat: bool): + deriv = root / "derivatives" + if flat: + for sid in ("sub-01", "sub-02", "sub-03"): + (deriv / sid / "pipeline" / "data").mkdir(parents=True) + else: + for g, sid in (("WT", "s1"), ("WT", "s2"), ("KO", "s3")): + (deriv / g / sid / "pipeline" / "data").mkdir(parents=True) + return deriv + + +def _run_cli(*args, cwd=None): + return subprocess.run( + [sys.executable, "-m", "source_analytics.cli", *args], + capture_output=True, text=True, cwd=cwd, + ) + + +def test_init_grouped_layout_writes_parseable_config(tmp_path): + root = tmp_path / "rest_roi" + _fake_reconstruction(root, flat=False) + out = _run_cli("init", str(root), "--name", "demo") + assert out.returncode == 0, out.stderr + path = root / "analysis" / "demo.yaml" + assert path.exists(), out.stderr + data = yaml.safe_load(path.read_text()) + assert data["design"]["levels"] == ["KO", "WT"] + assert [h["name"] for h in data["hypotheses"]] == ["WT_vs_KO"] + assert data["hypotheses"][0]["weights"] == {"WT": 1, "KO": -1} + assert list(data["paradigms"]) == ["resting"] + assert "subjects" not in data["paradigms"]["resting"] # grouped layout: dirs give groups + cfg = StudyConfig.from_yaml(path) + assert cfg.has_paradigms + assert [c.name for c in cfg.contrasts] == ["WT_vs_KO"] + scoped = cfg.for_paradigm_analysis("resting", "roi_psd") + assert scoped.output_dir == root / "analysis" / "analytics" / "resting" + + +def test_init_flat_layout_with_groups_from_and_stdout(tmp_path): + root = tmp_path / "rest_roi" + _fake_reconstruction(root, flat=True) + sl = tmp_path / "sl.yaml" + sl.write_text(yaml.safe_dump({"subjects": [ + {"id": "01", "group": "A"}, {"id": "02", "group": "B"}, {"id": "03", "group": "C"}]})) + out = _run_cli("init", str(root), "--groups-from", str(sl), "--output", "-", + "--analyses", "roi_psd,roi_connectivity") + assert out.returncode == 0, out.stderr + data = yaml.safe_load(out.stdout) # stdout is pure YAML + assert data["paradigms"]["resting"]["subjects"] == {"sub-01": "A", "sub-02": "B", "sub-03": "C"} + names = [h["name"] for h in data["hypotheses"]] + assert names[0] == "group_omnibus" and len(names) == 4 + assert list(data["paradigms"]["resting"]["analyses"]) == ["roi_psd", "roi_connectivity"] + assert not (root / "analysis").exists() + + +# ---- #10: --jobs 0 / -1 reach the resolver; #16: --force wipes ------------- +def test_cli_jobs_default_is_none(): + from source_analytics import cli + parser_ns = cli.main.__globals__ # sanity: module loaded + assert "canonical_analysis_name" in parser_ns + import argparse + # Build the parser the same way main() does by invoking with --help is noisy; + # instead check the argparse default via a dry parse of `run`. + out = _run_cli("run", "--help") + assert "--jobs" in out.stdout and "roi_connectivity" in out.stdout + + +def test_prepare_output_force_wipes_published_and_working_dirs(tmp_path): + from source_analytics.cli import _prepare_output + + class _Cfg: + output_dir = tmp_path / "analytics" / "resting" + results_dir = tmp_path / "results" + paradigm_name = "resting" + + work = _Cfg.output_dir / "roi_psd" / "data" + tbl = _Cfg.results_dir / "tables" / "resting" / "roi_psd" + fig = _Cfg.results_dir / "figures" / "resting" / "roi_psd" + for d in (work, tbl, fig): + d.mkdir(parents=True) + (d / "x.csv").write_text("stale") + + # deprecated alias resolves to the canonical dir + target = _prepare_output(_Cfg, "psd", strict=False, force=True, steps=None) + assert target == _Cfg.output_dir / "roi_psd" + assert not work.exists() and not tbl.exists() and not fig.exists() + + # --steps without process keeps data/, clears published dirs + for d in (work, tbl): + d.mkdir(parents=True) + (d / "x.csv").write_text("stale") + _prepare_output(_Cfg, "roi_psd", strict=False, force=True, steps={"statistics"}) + assert work.exists() and not tbl.exists() + + # --strict-output without --force errors on existing output + with pytest.raises(SystemExit): + _prepare_output(_Cfg, "roi_psd", strict=True, force=False, steps=None) + + +# ---- #8: vertex_spatial is retired end to end ------------------------------ +def test_vertex_spatial_processes_nothing(tmp_path): + from source_analytics.analyses.vertex_spatial_analysis import VertexSpatialAnalysis + + class _Cfg: + raw = {} + vertex = {} + name = "t" + results_dir = tmp_path / "results" + paradigm_name = None + roi_categories = {} + atlas_dir = None + + a = VertexSpatialAnalysis.__new__(VertexSpatialAnalysis) + a.config = _Cfg() + a.output_dir = tmp_path / "vertex_spatial" + a.output_dir.mkdir() + a._warned = False + a.setup() + a.process_subject(object()) # must not touch a loader + a.statistics() + a.summary() + assert (tmp_path / "results" / "tables" / "vertex_spatial" / "vertex_spatial_results.csv").exists() + assert "RETIRED" in (a.output_dir / "ANALYSIS_SUMMARY.md").read_text() + + +# ---- #13: vertex_evoked is on the hypothesis contract ---------------------- +def test_vertex_evoked_selectable_hypothesis(): + from source_analytics.analyses.vertex_evoked_analysis import VertexEvokedAnalysis + assert "hypothesis" in VertexEvokedAnalysis.SELECTABLE + + +# ---- #7: R scripts are discoverable from an installed prefix -------------- +def test_find_r_script_dir_env_override(tmp_path, monkeypatch): + from source_analytics.analyses.base import find_r_script_dir + d = tmp_path / "Rdir"; d.mkdir() + monkeypatch.setenv("SOURCE_ANALYTICS_R_DIR", str(d)) + assert find_r_script_dir() == d + + +# ---- roi_connectivity honours the YAML metrics: list ---------------------- +def test_roi_connectivity_reads_config_metrics(): + from source_analytics.analyses.roi_connectivity_analysis import ConnectivityAnalysis + + class _Cfg: + raw = {"roi_connectivity": {"metrics": ["aec", "pli"]}} + + a = ConnectivityAnalysis.__new__(ConnectivityAnalysis) + a.config = _Cfg() + a._selection = {} + a._edge_rows = [] + a.setup() + assert a._metrics == [m for m in ConnectivityAnalysis._ROI_METRICS if m in ("aec", "pli")] + + a.config.raw = {"roi_connectivity": {"metrics": ["bogus"]}} + with pytest.raises(ValueError, match="unknown metrics"): + a.setup() diff --git a/tests/test_edge.py b/tests/test_edge.py index 2d5cfd4..9c48396 100644 --- a/tests/test_edge.py +++ b/tests/test_edge.py @@ -140,3 +140,52 @@ def test_nbs_permutation_test_exposes_component_nodes(): # node membership is aligned 1:1 with component sizes/pvalues assert len(res.component_nodes) == len(res.component_sizes) assert all(isinstance(nodes, list) and len(nodes) >= 2 for nodes in res.component_nodes) + + +def test_nbs_significant_edges_mask_is_populated(): + """A significant component must be reflected in ``significant_edges``. + + Regression: the mask was allocated under a "Build significant edge mask" + comment and then never filled, so ``nbs_permutation_test`` returned an + all-False mask even for a large component at p = 0.0. Nothing in the shipped + pipeline reads the field (``_network_base`` uses ``component_nodes`` / + ``component_sizes``), so it was dormant — but it is a public field of the + returned dataclass and any new consumer silently gets "no edges". + + The mask must equal the suprathreshold edges of the significant components: + components partition the nodes, so a component's edges are the + suprathreshold edges with both endpoints in its node set. + """ + import numpy as np + + from source_analytics.stats.graph_metrics import nbs_permutation_test + + rng = np.random.default_rng(4) + n, n_sub = 10, 16 + + def net(shift): + m = np.abs(rng.normal(0.5, 0.1, size=(n, n))) + m = (m + m.T) / 2 + # a dense block of genuinely different edges between groups + m[:5, :5] += shift + m = (m + m.T) / 2 + np.fill_diagonal(m, 0.0) + return m + + A = [net(0.6) for _ in range(n_sub)] + B = [net(0.0) for _ in range(n_sub)] + res = nbs_permutation_test(A, B, nbs_threshold=2.0, n_permutations=500, seed=7) + + assert res.n_significant_components > 0, "fixture produced no component" + mask = np.triu(res.significant_edges, 1) + assert mask.sum() > 0, "significant component but an empty edge mask" + + # the mask must be a subset of the suprathreshold edges, and must exactly + # account for the significant components' sizes + supra = np.triu(np.abs(res.t_matrix) > 2.0, 1) + assert not (mask & ~supra).any(), "mask contains sub-threshold edges" + expected = sum(sz for sz, p in zip(res.component_sizes, res.component_pvalues) + if p < 0.05) + assert int(mask.sum()) == expected, ( + f"mask has {int(mask.sum())} edges, significant components total {expected}" + ) diff --git a/tests/test_evoked_hypotheses.py b/tests/test_evoked_hypotheses.py new file mode 100644 index 0000000..c8044a5 --- /dev/null +++ b/tests/test_evoked_hypotheses.py @@ -0,0 +1,243 @@ +"""Verification of the declared-hypothesis wiring on the evoked modules. + +``roi_evoked`` / ``electrode_evoked`` were the two modules +``docs/methods/HYPOTHESIS.md`` §9 listed as deferred, for "long-format DV": they +export one ``value`` column faceted by ``measure_name`` instead of one column +per DV. The wiring resolves that by making the measure a FACET, so each measure +gets its own FDR family across the spatial grid. + +Checks (synthetic measure table, 3 groups x 2 measures x 3 ROIs, no band axis): + 1. the additive ``_hypotheses.csv`` is written, with a row per + (hypothesis x measure x spatial) cell; + 2. spatial AND measure specificity — the planted cell is significant, and + neither a null ROI nor the null measure is; + 3. the band coordinate is null, because evoked measures carry no band axis; + 4. each measure is its own FDR family (a null measure cannot dilute a + signal measure, nor vice versa); + 5. the CSV fallback drives ``--steps statistics`` with no rows in memory; + 6. a study that renames the design factor still resolves; + 7. nothing is written when no hypotheses are declared; + 8. ``electrode_evoked`` behaves identically on ``channel``. + +Run: uv run --extra dev pytest tests/test_evoked_hypotheses.py -q +""" + +from __future__ import annotations + +from types import SimpleNamespace + +import numpy as np +import pandas as pd + +from source_analytics.analyses._evoked_hypotheses import write_evoked_hypotheses +from source_analytics.config import DesignSpec, Hypothesis + +GROUPS = ["Vehicle", "AUT00206", "AUT00201"] +ROIS = ["Auditory_L", "Auditory_R", "Visual_L"] +MEASURES = ["itc_onset", "erp_n1_amp"] + +# The planted effect: AUT00206 > Vehicle, on itc_onset, in Auditory_L only. +SIGNAL = ("itc_onset", "Auditory_L", "AUT00206", 1.4) + + +def _measure_rows(spatial_col: str = "roi", seed: int = 7) -> list[dict]: + """Long evoked rows; one measure x one spatial unit x one group carries signal.""" + rng = np.random.default_rng(seed) + m_sig, sp_sig, g_sig, amp = SIGNAL + rows = [] + sid = 0 + for g in GROUPS: + for _ in range(12): + sid += 1 + uid = f"m{sid:03d}" + for sp in ROIS: + for m in MEASURES: + val = 1.0 + rng.normal(0, 0.5) + if m == m_sig and sp == sp_sig and g == g_sig: + val += amp + rows.append({ + "subject": uid, "group": g, spatial_col: sp, + "measure_name": m, "measure_type": "itc", + "band_lo": np.nan, "band_hi": np.nan, + "time_lo": 0.0, "time_hi": 0.1, + "value": val, "n_epochs": 40, + }) + return rows + + +def _hyp() -> Hypothesis: + return Hypothesis(name="drug_effect", kind="contrast", + weights={"AUT00206": 1.0, "Vehicle": -1.0}) + + +def _analysis(tmp_path, *, rows, spatial_col="roi", name="roi_evoked", + hyps=None, factor="group", selection=None): + """Minimal stand-in for the BaseAnalysis surface the wiring touches.""" + spec = DesignSpec(factor=factor, reference="Vehicle", levels=GROUPS, + hypotheses=_hyp_list(hyps)) + tbl = tmp_path / "tables" + tbl.mkdir(parents=True, exist_ok=True) + return SimpleNamespace( + name=name, + config=SimpleNamespace(design_spec=spec, bands={}), + tbl_dir=tbl, + output_dir=tmp_path, + _selection=selection or {}, + _measure_rows=rows, + ) + + +def _hyp_list(hyps): + if hyps is None: + return [_hyp()] + return hyps + + +def _run(tmp_path, **kw): + spatial_col = kw.pop("spatial_col", "roi") + name = kw.get("name", "roi_evoked") + an = _analysis(tmp_path, spatial_col=spatial_col, **kw) + write_evoked_hypotheses( + an, spatial_col=spatial_col, measures_csv=f"{name}_measures.csv" + ) + return an.tbl_dir / f"{name}_hypotheses.csv" + + +# --------------------------------------------------------------------------- # + +def test_writes_table_with_cell_per_measure_and_spatial(tmp_path): + out = _run(tmp_path, rows=_measure_rows()) + assert out.exists(), "additive hypotheses CSV was not written" + df = pd.read_csv(out) + # one hypothesis x 2 measures x 3 ROIs + assert len(df) == len(MEASURES) * len(ROIS) + assert set(df["measure_name"]) == set(MEASURES) + assert set(df["spatial"]) == set(ROIS) + assert set(df["hypothesis"]) == {"drug_effect"} + + +def test_measure_and_spatial_specificity(tmp_path): + df = pd.read_csv(_run(tmp_path, rows=_measure_rows())) + m_sig, sp_sig, _, _ = SIGNAL + + def cell(m, sp): + r = df[(df.measure_name == m) & (df.spatial == sp)] + assert len(r) == 1 + return r.iloc[0] + + # planted cell recovered, and in the planted direction + hit = cell(m_sig, sp_sig) + assert bool(hit["significant"]), "planted cell not recovered" + assert hit["estimate"] > 0 + + # same measure, other ROIs: null + for sp in ROIS: + if sp != sp_sig: + assert not bool(cell(m_sig, sp)["significant"]) + # the other measure: null everywhere, including the signal ROI + for sp in ROIS: + assert not bool(cell("erp_n1_amp", sp)["significant"]) + + +def test_band_coordinate_is_null(tmp_path): + """Evoked measures fix their own band/time window, so there is no band axis.""" + df = pd.read_csv(_run(tmp_path, rows=_measure_rows())) + assert df["band"].isna().all() + + +def test_each_measure_is_its_own_fdr_family(tmp_path): + """A null measure must not dilute a signal measure's correction, or vice versa.""" + df = pd.read_csv(_run(tmp_path, rows=_measure_rows())) + fam_sizes = df.groupby("measure_name")["fdr_family"].nunique() + assert (fam_sizes == 1).all(), "a measure's cells split across FDR families" + + # The family is the spatial grid within one measure: q == BH over 3 ROIs. + for m in MEASURES: + sub = df[df.measure_name == m].sort_values("spatial") + p = sub["p_value"].to_numpy() + n = len(p) + order = np.argsort(p) + ranked = p[order] * n / (np.arange(n) + 1) + expect = np.minimum.accumulate(ranked[::-1])[::-1] + got = sub["q_value"].to_numpy()[order] + assert np.allclose(got, expect, atol=1e-12), f"{m}: q is not BH over its ROIs" + + +def test_csv_fallback_when_no_rows_in_memory(tmp_path): + """--steps statistics on its own reads the exported measures CSV.""" + data = tmp_path / "data" + data.mkdir() + pd.DataFrame(_measure_rows()).to_csv(data / "roi_evoked_measures.csv", index=False) + + out = _run(tmp_path, rows=[]) # nothing in memory + assert out.exists() + df = pd.read_csv(out) + m_sig, sp_sig, _, _ = SIGNAL + hit = df[(df.measure_name == m_sig) & (df.spatial == sp_sig)].iloc[0] + assert bool(hit["significant"]) + + +def test_missing_rows_and_missing_csv_is_a_noop(tmp_path): + out = _run(tmp_path, rows=[]) + assert not out.exists() + + +def test_renamed_design_factor_resolves(tmp_path): + """Rows always carry 'group'; a study naming its factor otherwise still works.""" + out = _run(tmp_path, rows=_measure_rows(), factor="treatment") + df = pd.read_csv(out) + m_sig, sp_sig, _, _ = SIGNAL + assert bool(df[(df.measure_name == m_sig) & (df.spatial == sp_sig)].iloc[0]["significant"]) + + +def test_no_declared_hypotheses_writes_nothing(tmp_path): + out = _run(tmp_path, rows=_measure_rows(), hyps=[]) + assert not out.exists() + + +def test_hypothesis_selection_filters(tmp_path): + other = Hypothesis(name="other_effect", kind="contrast", + weights={"AUT00201": 1.0, "Vehicle": -1.0}) + out = _run(tmp_path, rows=_measure_rows(), hyps=[_hyp(), other], + selection={"hypothesis": frozenset({"drug_effect"})}) + df = pd.read_csv(out) + assert set(df["hypothesis"]) == {"drug_effect"} + + +def test_electrode_arm_behaves_identically_on_channel(tmp_path): + """The electrode comparator is the same wiring over 'channel'.""" + out = _run(tmp_path, rows=_measure_rows(spatial_col="channel"), + spatial_col="channel", name="electrode_evoked") + df = pd.read_csv(out) + m_sig, sp_sig, _, _ = SIGNAL + assert len(df) == len(MEASURES) * len(ROIS) + assert bool(df[(df.measure_name == m_sig) & (df.spatial == sp_sig)].iloc[0]["significant"]) + + +def test_missing_spatial_column_is_a_noop(tmp_path): + rows = _measure_rows() + for r in rows: + r.pop("roi") + out = _run(tmp_path, rows=rows) + assert not out.exists() + + +def test_fdr_family_label_separates_the_measures(tmp_path): + """Two measures are two FDR families, and must not share one label. + + Regression for the facet-identity defect: ``fdr_family`` encodes family + IDENTITY so q-values from different families are provably non-comparable + (REPORT_PLAN §10b). The member hash cannot do it here — both families span + the same ROI cells — so the facet has to appear in the label. Before the fix + both measures emitted ``key=all|drug_effect|NA members=cell[3] hash=...``, + byte-identical across genuinely different families. + """ + df = pd.read_csv(_run(tmp_path, rows=_measure_rows())) + per_measure = df.groupby("measure_name")["fdr_family"].agg(set) + assert all(len(s) == 1 for s in per_measure), "a measure split across labels" + labels = {m: next(iter(s)) for m, s in per_measure.items()} + assert len(set(labels.values())) == len(MEASURES), ( + f"distinct FDR families share one label: {labels}" + ) + for m, lab in labels.items(): + assert f"measure_name={m}" in lab, f"facet identity missing from {lab!r}" diff --git a/tests/test_parallel_subjects.py b/tests/test_parallel_subjects.py index 38719f4..55a435d 100644 --- a/tests/test_parallel_subjects.py +++ b/tests/test_parallel_subjects.py @@ -79,6 +79,8 @@ def test_resolve_jobs_auto_and_config(): class _Cfg: raw = {"jobs": 3} a.config = _Cfg() - assert a._resolve_jobs(1) == 3 # config fallback when CLI=1 + assert a._resolve_jobs(None) == 3 # CLI not given -> config `jobs:` + assert a._resolve_jobs(1) == 1 # explicit CLI 1 wins over config assert a._resolve_jobs(2) == 2 # CLI wins over config assert a._resolve_jobs(-1) >= 1 # auto (all-but-one core) + assert a._resolve_jobs(0) >= 1 # auto (all-but-one core) diff --git a/tests/test_r_scripts_smoke.py b/tests/test_r_scripts_smoke.py new file mode 100644 index 0000000..e1778d2 --- /dev/null +++ b/tests/test_r_scripts_smoke.py @@ -0,0 +1,187 @@ +"""End-to-end smoke tests for the R entry points against synthetic CSVs. + +These pin the audit fixes on the R side: + +- the evoked scripts derive their contrast list from ``design:``/``hypotheses:`` + (they used to iterate ``config$contrasts``, which is NULL for modern configs); +- ``roi_connectivity_analysis.R`` runs to completion on a single-metric CSV + (it used to require ``coherence`` + ``imag_coherence`` columns); +- ``roi_cross_freq_edges_analysis.R`` gives AAC/PPC a hypothesis-layer path; +- ``roi_transfer_entropy_analysis.R`` writes ``roi_directed_*`` tables and + tests a DTF-only edge CSV. + +Skipped when ``Rscript`` is not on PATH. +""" + +from __future__ import annotations + +import itertools +import shutil +import subprocess +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest +import yaml + +R_DIR = Path(__file__).resolve().parent.parent / "R" + +pytestmark = pytest.mark.skipif(shutil.which("Rscript") is None, reason="Rscript not installed") + +GROUPS = {"WT_VEH": 6, "KO_VEH": 6} +ROIS = ["Motor_L", "Motor_R", "Hipp_L", "Hipp_R", "Thal_L", "Thal_R"] +BANDS = {"Theta": [4, 10], "Alpha": [10, 13], "Low Gamma": [30, 55]} + + +def _subjects(): + return [(f"{g}_s{i:02d}", g) for g, n in GROUPS.items() for i in range(n)] + + +def _run(script: str, data_dir: Path, config: Path, out: Path, *extra: str) -> str: + cmd = [ + "Rscript", str(R_DIR / script), + "--data-dir", str(data_dir), "--config", str(config), + "--output-dir", str(out), "--fig-dir", str(out / "figures"), + "--tbl-dir", str(out / "tables"), "--no-figures", *extra, + ] + res = subprocess.run(cmd, capture_output=True, text=True, timeout=600) + assert res.returncode == 0, f"{script} failed:\n{res.stdout}\n{res.stderr}" + return res.stdout + res.stderr + + +@pytest.fixture +def design_config(tmp_path) -> Path: + """A modern config: design:/hypotheses: only, no legacy contrasts: block.""" + cfg = { + "name": "Smoke", + "groups": {"WT_VEH": "WT", "KO_VEH": "KO"}, + "group_order": ["WT_VEH", "KO_VEH"], + "group_colors": {"WT_VEH": "#3498DB", "KO_VEH": "#E74C3C"}, + "bands": BANDS, + "roi_categories": {"Motor": ROIS[:2], "Hipp": ROIS[2:4], "Thal": ROIS[4:]}, + "sfreq": 500, + "design": {"factor": "group", "reference": "WT_VEH", "levels": ["WT_VEH", "KO_VEH"]}, + "hypotheses": [ + {"name": "group_omnibus", "kind": "omnibus"}, + {"name": "disease_effect", "kind": "contrast", + "weights": {"KO_VEH": 1, "WT_VEH": -1}, "label": "KO vs WT"}, + ], + } + p = tmp_path / "config.yaml" + p.write_text(yaml.safe_dump(cfg)) + return p + + +def test_evoked_scripts_derive_contrasts_from_design_spec(tmp_path, design_config): + rng = np.random.default_rng(0) + for level, col, units in (("roi", "roi", ROIS), ("electrode", "channel", ["Fz", "Cz", "Pz"])): + rows = [] + for uid, g in _subjects(): + for m, mt in (("itc_theta", "itc"), ("ersp_gamma", "ersp")): + for u in units: + rows.append({ + "subject": uid, "group": g, col: u, "measure_name": m, + "measure_type": mt, "band_lo": 4, "band_hi": 10, + "time_lo": 0.0, "time_hi": 0.3, + "value": rng.normal(0.5 + (0.2 if g == "KO_VEH" else 0), 0.1), + "n_epochs": 50, + }) + out = tmp_path / level + data = out / "data" + data.mkdir(parents=True) + pd.DataFrame(rows).to_csv(data / f"{level}_evoked_measures.csv", index=False) + + log = _run(f"{level}_evoked_analysis.R", data, design_config, out) + assert "Contrasts: 1" in log + omnibus = pd.read_csv(out / "tables" / f"{level}_evoked_omnibus.csv") + assert len(omnibus) > 0, "descriptive LMM table is empty — contrast loop did not run" + assert set(omnibus["contrast"]) == {"disease_effect"} + + +def test_connectivity_script_runs_on_single_metric_csv(tmp_path, design_config): + rng = np.random.default_rng(1) + rows = [] + for uid, g in _subjects(): + for b in BANDS: + for r1, r2 in itertools.combinations(ROIS, 2): + rows.append({"subject": uid, "group": g, "band": b, "roi1": r1, "roi2": r2, + "aec": rng.normal(0.2, 0.03)}) + out = tmp_path / "conn" + data = out / "data" + data.mkdir(parents=True) + pd.DataFrame(rows).to_csv(data / "roi_connectivity_edges.csv", index=False) + + log = _run("roi_connectivity_analysis.R", data, design_config, out) + assert "Metrics present: aec" in log + summary = (out / "ANALYSIS_SUMMARY.md").read_text() + assert "AEC" in summary + assert "Imaginary Coherence" not in summary + + +def test_cross_freq_edges_script_writes_hypothesis_tables(tmp_path, design_config): + rng = np.random.default_rng(2) + pairs = [("Theta", "Low Gamma"), ("Alpha", "Low Gamma")] + data = tmp_path / "xf" / "data" + data.mkdir(parents=True) + for metric in ("aac", "ppc"): + rows = [] + for uid, g in _subjects(): + for pb, ab in pairs: + for rx in ROIS: + for ry in ROIS: + row = {"subject": uid, "group": g, "phase_band": pb, "amp_band": ab, + "freq_pair": f"{pb}-{ab}", "n": 4, "m": 1, + "roi_x": rx, "roi_y": ry, + metric: rng.normal(0.3 + (0.15 if g == "KO_VEH" else 0), 0.05)} + if metric == "ppc": + row["ppc_z"] = rng.normal(1, 1) + rows.append(row) + pd.DataFrame(rows).to_csv(data / f"{metric}_edges.csv", index=False) + + out = tmp_path / "xf" + _run("roi_cross_freq_edges_analysis.R", data, design_config, out, "--metric", "aac,ppc") + tbl = out / "tables" + for m in ("aac", "ppc"): + for tier in ("global", "directed_edges", "region"): + path = tbl / f"roi_cross_freq_{m}_{tier}_hypotheses.csv" + assert path.exists(), path.name + df = pd.read_csv(path) + assert {"hypothesis", "kind", "band", "q_value", "significant"} <= set(df.columns) + assert set(df["band"]) <= {f"{pb}-{ab}" for pb, ab in pairs} + # PPC carries the surrogate z as a second DV + ppc_edges = pd.read_csv(tbl / "roi_cross_freq_ppc_directed_edges_hypotheses.csv") + assert set(ppc_edges["dv"]) == {"ppc", "ppc_z"} + assert "Amplitude-Amplitude Coupling" in (out / "ANALYSIS_SUMMARY.md").read_text() + + +@pytest.mark.parametrize("cols", [["dtf"], ["te", "net_te", "dtf"]]) +def test_directed_script_uses_canonical_prefix_and_tests_every_dv(tmp_path, design_config, cols): + rng = np.random.default_rng(3) + rows = [] + for uid, g in _subjects(): + for b in BANDS: + for r1 in ROIS: + for r2 in ROIS: + if r1 == r2: + continue + row = {"subject": uid, "group": g, "band": b, + "source_roi": r1, "target_roi": r2} + for c in cols: + row[c] = rng.normal(0.1 + (0.03 if g == "KO_VEH" else 0), 0.02) + rows.append(row) + out = tmp_path / "dir" + data = out / "data" + data.mkdir(parents=True) + pd.DataFrame(rows).to_csv(data / "roi_transfer_entropy_edges.csv", index=False) + + _run("roi_transfer_entropy_analysis.R", data, design_config, out) + tbl = out / "tables" + assert not list(tbl.glob("roi_transfer_entropy_*")) + for name in ("roi_directed_global_hypotheses.csv", + "roi_directed_directed_edges_hypotheses.csv", + "roi_directed_region_hypotheses.csv", + "roi_directed_omnibus_lmm.csv"): + assert (tbl / name).exists(), name + edges = pd.read_csv(tbl / "roi_directed_directed_edges_hypotheses.csv") + assert set(edges["dv"]) == set(cols) diff --git a/uv.lock b/uv.lock index a8063f6..24ca65a 100644 --- a/uv.lock +++ b/uv.lock @@ -2042,6 +2042,8 @@ version = "0.6.0" source = { editable = "." } dependencies = [ { name = "joblib" }, + { name = "matplotlib", version = "3.10.9", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, + { name = "matplotlib", version = "3.11.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, { name = "numpy", version = "2.4.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version == '3.11.*'" }, { name = "numpy", version = "2.5.2", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.12'" }, @@ -2095,6 +2097,7 @@ requires-dist = [ { name = "black", marker = "extra == 'dev'" }, { name = "joblib", specifier = ">=1.3" }, { name = "markdown", marker = "extra == 'render'", specifier = ">=3.4" }, + { name = "matplotlib", specifier = ">=3.7" }, { name = "mne", marker = "extra == 'mne'", specifier = ">=1.5" }, { name = "networkx", marker = "extra == 'network'", specifier = ">=3.0" }, { name = "nibabel", marker = "extra == 'atlas'", specifier = ">=5.0" },