From 631e22f5cb5896f28a30405539ba19072ded9b93 Mon Sep 17 00:00:00 2001 From: Aayush Gandhi Date: Fri, 7 Aug 2026 17:12:19 -0400 Subject: [PATCH 1/5] Document NARPS (ds001734) before wiring the loader Verified the format against the dataset's own files on OpenNeuro rather than memory: task MGT, TR 1.0 s, 4 runs x 64 trials, events.tsv columns (onset, duration, gain, loss, RT, participant_response) with responses strongly/weakly_accept, strongly/weakly_reject, NoResp; participants.tsv group column (equalIndifference / equalRange). docs/narps.md records the mapping onto the framework (accept/reject as the individual outcome, the (gain,loss) gamble as the stimulus key, population acceptance rate as the aggregate), the traps, and an honest statement of NARPS's role. Key point: on gambles the economic baseline dominates the aggregate arm by construction, so NARPS is an individual-level fMRI plumbing check, NOT a brain-beats-behaviour demonstration. Forcing that headline would misrepresent the task. --- docs/data_sources.md | 7 +- docs/narps.md | 161 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 167 insertions(+), 1 deletion(-) create mode 100644 docs/narps.md diff --git a/docs/data_sources.md b/docs/data_sources.md index c7b5acc..28c34d9 100644 --- a/docs/data_sources.md +++ b/docs/data_sources.md @@ -156,7 +156,12 @@ institution's ethics approval covers linking lab data to public platform data. (`io/deap_market.py`), so it can exercise *both* levels, not just the individual one. Validates the EEG and behaviour paths and the whole ensemble. 2. **NARPS `ds001734`** to validate the fMRI loader and ROI extraction on real - BIDS data with a reward task and a decent sample. + BIDS data with a reward task and a decent sample. **The loader for this is + built** (`io/narps.py`, [`docs/narps.md`](narps.md)); it reuses `FMRILoader`'s + NAcc/vmPFC/AIns sphere extraction and needs only the OpenNeuro download (no + licence). Note NARPS is an *individual-level* validation: on gambles the + economic baseline dominates the aggregate arm by construction, so it is the + fMRI plumbing check, not a brain-beats-behaviour demonstration. 3. **A market-linked stimulus set** (Kiva or Kickstarter items) for a purpose-built aggregate arm. This is the step that requires collecting your own scans, and it is the one that produces the fully novel result. diff --git a/docs/narps.md b/docs/narps.md new file mode 100644 index 0000000..6182a3c --- /dev/null +++ b/docs/narps.md @@ -0,0 +1,161 @@ +# NARPS (ds001734) + +How this framework loads the NARPS mixed-gambles dataset, the format details that +matter, and an honest statement of what NARPS can and cannot show. + +NARPS (Botvinik-Nezer et al., 2019, *Scientific Data*, +[10.1038/s41597-019-0113-7](https://doi.org/10.1038/s41597-019-0113-7); dataset +[ds001734 on OpenNeuro](https://openneuro.org/datasets/ds001734), DOI +[10.18112/openneuro.ds001734.v1.0.5](https://doi.org/10.18112/openneuro.ds001734.v1.0.5)) +is the fMRI counterpart to DEAP in this project: a reward-circuit task with +enough subjects to matter. On each trial a participant sees a 50/50 gamble with a +possible **gain** and a possible **loss** and decides whether to accept it. This +is exactly the anticipatory-affect paradigm the framework's fMRI ROIs come from, +so it is the right dataset to validate the fMRI loader and the NAcc / vmPFC / AIns +sphere extraction on real data. + +**Access.** NARPS is public on OpenNeuro; no licence agreement. It is large +(4-D BOLD across ~108 subjects), so the loader downloads nothing and expects a +local BIDS tree, ideally with fMRIPrep-preprocessed derivatives. + +--- + +## Verified format + +Confirmed against the dataset's own files on OpenNeuro (2026-08-07), not from +memory. A real events file looks like: + +``` +onset duration gain loss RT participant_response +4.071 4 14 6 2.388 weakly_accept +11.834 4 34 14 2.289 strongly_accept +27.535 4 10 10 1.457 weakly_reject +36.435 4 12 19 1.973 strongly_reject +``` + +- **Task**: `MGT` (mixed gambles task). +- **TR**: 1.0 s (`task-MGT_bold.json`, `RepetitionTime: 1`). +- **Trials**: 4 runs per subject, 64 trials per run, 4 s cue per trial. +- **`events.tsv` columns**: `onset`, `duration`, `gain`, `loss`, `RT`, + `participant_response`. +- **Response values**: `strongly_accept`, `weakly_accept`, `weakly_reject`, + `strongly_reject`, and `NoResp` for no-response trials (RT coded 0). +- **`participants.tsv` columns**: `participant_id`, `group`, `gender`, `age`. +- **Groups**: `equalIndifference` (gains 10-40, losses 5-20) and `equalRange` + (gains and losses both 5-20). The two groups saw different gamble matrices; + this is a real design factor, not noise (see below). + +--- + +## How NARPS maps onto the framework + +| framework concept | NARPS realisation | +|---|---| +| individual outcome `y_individual` | accept vs reject the gamble | +| stimulus key `stimulus_id` | the gamble, keyed by `(gain, loss)` | +| aggregate outcome `y_aggregate` | population acceptance rate per gamble | +| fMRI block | NAcc / vmPFC / AIns sphere betas per trial | +| behaviour block | gain, loss, expected value, RT | + +**Individual choice.** `strongly_accept` and `weakly_accept` collapse to accept +(1); `strongly_reject` and `weakly_reject` to reject (0). `NoResp` trials are +dropped, not imputed: a missing response is missing data, and guessing it would +manufacture labels. The strong/weak distinction is preserved in a separate +confidence field for anyone who wants a 4-level analysis, but the binary accept +is the default outcome. + +**Stimulus key.** The gamble is defined by its `(gain, loss)` pair, and those +pairs recur across subjects, which is what makes a stimulus key and an aggregate +outcome possible at all. Keys look like `g14_l06`. + +--- + +## The traps + +Fewer than DEAP, but the ones here are sharp. + +**1. The behaviour arm is *supposed* to be strong, and that inverts the usual +story.** In the crowdfunding and video paradigms the neural signal beats the +behavioural measure. Gambles are different: acceptance is largely a function of +gain and loss through expected value, so an economic model built from `gain` and +`loss` forecasts aggregate acceptance almost by construction. On NARPS, expect +**behaviour to win the aggregate arm**. That is not a failure of the neural arm; +it is what gambles are. NARPS is primarily an **individual-level** validation +(can NAcc/vmPFC betas predict this person's accept/reject), and the aggregate arm +is a secondary check where the economic baseline is expected to dominate. Do not +report a forced brain-beats-behaviour headline on NARPS; it would be dishonest +about the task. + +**2. The two groups saw different gamble matrices.** `equalIndifference` and +`equalRange` do not share the same `(gain, loss)` grid. Pooling them for the +aggregate arm mixes two stimulus spaces, and a gamble present in one group may be +absent in the other. The loader keeps the group label and, by default, warns when +you pool across groups for aggregate forecasting. Analyse within group unless you +have a reason not to. + +**3. Slice-timing and the haemodynamic lag are real analytic choices.** The +default trial feature is a peak-window mean of the BOLD 4 to 8 s after cue onset, +which is a crude stand-in for a proper first-level GLM. It is fine for getting the +pipeline running and wrong for a publishable estimate. Replace it with +`nilearn.glm.first_level` before reporting. The loader records +`extraction="peak_window_mean"` in provenance so this is never hidden. + +**4. Confounds must be regressed, and this loader does not do it for you.** +Motion, framewise displacement, and the aCompCor components in the fMRIPrep +`*_desc-confounds_timeseries.tsv` are not optional for a reward-ROI analysis: +head motion correlates with task events and with individual differences. The +loader accepts a confounds table and passes it to the masker; if you do not +supply one, it warns. + +**5. No-response trials and RT = 0.** `NoResp` rows carry `RT = 0`. Feeding that +zero into an RT feature as if it were a fast response is wrong. The loader drops +`NoResp` trials before building any block. + +--- + +## The theory-specified ROI + +On the crowdfunding side the signal lived in NAcc (gain anticipation) and MPFC. +NARPS is built for exactly this circuit: the canonical mixed-gambles result is +that NAcc tracks potential gain and a vmPFC/MPFC region tracks the net expected +value, with anterior insula tracking potential loss. The framework's default +ROIs (`NAcc_L/R`, `MPFC`, `AIns_L/R`) are therefore already the right set, and +the NARPS loader reuses `FMRILoader`'s sphere extraction rather than defining its +own. Confirm the exact sphere coordinates against whatever prior you are +replicating before reporting; sphere placement is a genuine analytic degree of +freedom. + +--- + +## What NARPS validates, and what it does not + +**Validates**: that the fMRI loader reads a real BIDS reward task, extracts +trial-wise ROI features, and that those features predict individual accept/reject +above chance with subject-grouped CV. That is milestone 2 in +[data_sources.md](data_sources.md), and it is a real result on its own. + +**Does not validate**: the brain-beats-behaviour aggregate claim. Gambles are the +wrong task for that, for the reason in trap 1. NARPS is the fMRI plumbing check +and the individual-decoding check. The aggregate neuroforecasting story needs a +stimulus set with a market outcome that is *not* a deterministic function of the +stimulus parameters, which is what the crowdfunding/video/microlending paradigms +provide and gambles do not. + +--- + +## Usage sketch + +```python +from behavioral_decoding.io.narps import load_narps + +# Individual-level, one group, from fMRIPrep derivatives. +dataset = load_narps( + root="/path/to/ds001734", + derivatives="/path/to/fmriprep", + group="equalIndifference", + space="MNI152NLin2009cAsym", +) +``` + +See `scripts/run_narps.py --demo` for the whole path on a synthetic BIDS fixture +that needs no download. From 0f04502379064041da0e936fa554741b25221f0e Mon Sep 17 00:00:00 2001 From: Aayush Gandhi Date: Fri, 7 Aug 2026 17:14:29 -0400 Subject: [PATCH 2/5] Add NARPS core loader (events, keys, behaviour block) The neuroimaging-free layer of the NARPS loader: parses events.tsv, binarises participant_response (strongly/weakly_accept -> 1, *_reject -> 0), drops NoResp trials (RT=0 is missing data, not a fast response), keys each gamble by (gain, loss) as e.g. g14_l06, and builds the behaviour block from the economic variables (gain, loss, expected value, |EV|, RT). Also: acceptance_rate_by_gamble computes the population acceptance rate per gamble (the NARPS aggregate outcome), excluding gambles seen by too few subjects. from_events_and_rois assembles fMRI + behaviour blocks from a precomputed trial-by-ROI matrix, reusing FMRILoader; this is the entry point the nilearn .load path will call after extraction. Verified against a real ds001734 events file. No neuroimaging dependency in this layer. --- src/behavioral_decoding/io/__init__.py | 14 + src/behavioral_decoding/io/narps.py | 363 +++++++++++++++++++++++++ 2 files changed, 377 insertions(+) create mode 100644 src/behavioral_decoding/io/narps.py diff --git a/src/behavioral_decoding/io/__init__.py b/src/behavioral_decoding/io/__init__.py index f17af32..70c0447 100644 --- a/src/behavioral_decoding/io/__init__.py +++ b/src/behavioral_decoding/io/__init__.py @@ -25,6 +25,14 @@ from .eeg import EEGLoader from .face import FaceLoader from .fmri import FMRILoader +from .narps import ( + NARPS_TR, + NARPSFormatError, + NARPSLoader, + acceptance_rate_by_gamble, + gamble_key, + parse_events, +) from .registry import get_loader, register_loader __all__ = [ @@ -37,6 +45,7 @@ "FACE", "FMRI", "KNOWN_MODALITIES", + "NARPS_TR", "PERIPHERAL", "BaseLoader", "BehaviorLoader", @@ -47,8 +56,13 @@ "FMRILoader", "ModalityBlock", "MultimodalDataset", + "NARPSFormatError", + "NARPSLoader", + "acceptance_rate_by_gamble", "binarise_ratings", + "gamble_key", "get_loader", "load_deap", + "parse_events", "register_loader", ] diff --git a/src/behavioral_decoding/io/narps.py b/src/behavioral_decoding/io/narps.py new file mode 100644 index 0000000..972a202 --- /dev/null +++ b/src/behavioral_decoding/io/narps.py @@ -0,0 +1,363 @@ +"""NARPS (ds001734) loader: the mixed-gambles reward task. + +NARPS is the fMRI counterpart to DEAP in this project. On each trial a +participant sees a 50/50 gamble with a possible gain and a possible loss and +decides whether to accept it. That is the anticipatory-affect paradigm the +framework's fMRI ROIs come from, so NARPS is the dataset to validate the fMRI +loader and the NAcc / vmPFC / AIns sphere extraction on real data. + +See ``docs/narps.md`` for the verified format and the honest statement of what +NARPS can show. The one thing to keep in mind while reading this module: on +gambles the economic variables (gain, loss) forecast aggregate acceptance almost +by construction, so NARPS is an **individual-level** validation, not a +brain-beats-behaviour demonstration. The behaviour block here is a genuine, and +genuinely strong, comparator, not a strawman. + +This file has two layers. This one is the format and bookkeeping layer, with no +neuroimaging dependencies: it parses ``events.tsv``, binarises the response, +builds the gamble stimulus keys, and assembles the behaviour block. The BOLD -> +ROI extraction lives in the ``.load`` path (see ``narps_bold`` / +:meth:`NARPSLoader.load`), which reuses :class:`FMRILoader` and needs nilearn. +""" + +from __future__ import annotations + +import re +from pathlib import Path +from typing import Any, Dict, List, Optional, Sequence, Tuple + +import numpy as np + +from ..utils.logging import get_logger +from .base import BEHAVIOR, ModalityBlock + +logger = get_logger(__name__) + +# ------------------------------------------------------------------ constants + +NARPS_TASK = "MGT" +NARPS_TR = 1.0 +NARPS_TRIAL_DURATION_S = 4.0 +NARPS_RUNS = (1, 2, 3, 4) +NARPS_TRIALS_PER_RUN = 64 + +EVENTS_COLUMNS = ("onset", "duration", "gain", "loss", "RT", "participant_response") + +# The four graded responses plus the no-response marker, verbatim from the data. +ACCEPT_RESPONSES = ("strongly_accept", "weakly_accept") +REJECT_RESPONSES = ("strongly_reject", "weakly_reject") +NO_RESPONSE = "NoResp" +ALL_RESPONSES = ACCEPT_RESPONSES + REJECT_RESPONSES + (NO_RESPONSE,) + +# Strength of the graded response, for anyone wanting a confidence feature or a +# 4-level analysis. The binary accept is the default outcome. +RESPONSE_CONFIDENCE = { + "strongly_accept": 1.0, + "weakly_accept": 0.5, + "weakly_reject": 0.5, + "strongly_reject": 1.0, +} + +GROUPS = ("equalIndifference", "equalRange") + + +class NARPSFormatError(ValueError): + """Raised when a file does not match the documented NARPS layout. + + Loud on purpose: an events file with unexpected columns or response strings + is a different release or a corrupted download, and both produce wrong + labels rather than crashes if waved through. + """ + + +# ------------------------------------------------------------- gamble keys + + +def gamble_key(gain: float, loss: float) -> str: + """Stable stimulus id for a gamble, e.g. ``"g14_l06"``. + + The gamble is the stimulus, and ``(gain, loss)`` pairs recur across subjects, + which is the only thing that makes an aggregate outcome possible. Integer + coding matches the dataset, where gains and losses are whole currency units. + """ + return f"g{int(round(gain)):02d}_l{int(round(loss)):02d}" + + +def expected_value(gain: float, loss: float) -> float: + """Expected value of a 50/50 gain/loss gamble: ``0.5 * (gain - loss)``.""" + return 0.5 * (float(gain) - float(loss)) + + +# ------------------------------------------------------------ events parsing + + +def _read_events_table(path_or_table: Any) -> Any: + """Return a pandas DataFrame from a path or an existing frame.""" + try: + import pandas as pd + except ImportError as exc: # pragma: no cover - optional dependency + raise ImportError( + "reading NARPS events requires pandas (`pip install pandas`)" + ) from exc + + if isinstance(path_or_table, pd.DataFrame): + return path_or_table.copy() + return pd.read_csv(path_or_table, sep="\t") + + +def parse_events( + path_or_table: Any, + drop_no_response: bool = True, +) -> Any: + """Parse one ``*_events.tsv`` into a tidy, framework-ready table. + + Adds three derived columns and validates the response strings. Returns a + DataFrame with, at minimum: ``onset``, ``gain``, ``loss``, ``RT``, + ``participant_response``, ``accept`` (0/1), ``confidence``, ``stimulus_id``, + ``expected_value``. + + Parameters + ---------- + drop_no_response: + Drop ``NoResp`` trials. On by default. A missing response is missing + data; imputing it would manufacture a label, and its ``RT`` of 0 would + poison any RT feature. + """ + table = _read_events_table(path_or_table) + + missing = [c for c in EVENTS_COLUMNS if c not in table.columns] + if missing: + raise NARPSFormatError( + f"events file is missing columns {missing}; " + f"found {list(table.columns)}. Expected the NARPS mixed-gambles " + f"layout {list(EVENTS_COLUMNS)}." + ) + + responses = set(table["participant_response"].dropna().unique()) + unknown = responses - set(ALL_RESPONSES) + if unknown: + raise NARPSFormatError( + f"unrecognised participant_response values {sorted(unknown)}. " + f"Expected a subset of {list(ALL_RESPONSES)}. A different coding " + "means the accept/reject mapping here does not apply." + ) + + n_before = len(table) + if drop_no_response: + table = table[table["participant_response"] != NO_RESPONSE].copy() + n_dropped = n_before - len(table) + if n_dropped: + logger.info( + "NARPS: dropped %d/%d NoResp trials (RT=0; missing data, not a fast " + "response)", + n_dropped, + n_before, + ) + + table["accept"] = table["participant_response"].isin(ACCEPT_RESPONSES).astype(int) + table["confidence"] = table["participant_response"].map(RESPONSE_CONFIDENCE).astype(float) + table["expected_value"] = 0.5 * (table["gain"].astype(float) - table["loss"].astype(float)) + table["stimulus_id"] = [ + gamble_key(gain, loss) + for gain, loss in zip(table["gain"], table["loss"]) + ] + return table.reset_index(drop=True) + + +def subject_id_from_path(path: Path) -> str: + """``.../sub-014_task-MGT_run-01_events.tsv`` -> ``"sub-014"``.""" + match = re.search(r"(sub-[A-Za-z0-9]+)", path.name) + if not match: + raise NARPSFormatError( + f"{path.name!r} does not contain a BIDS subject label like 'sub-014'" + ) + return match.group(1) + + +def load_participants(path: str) -> Any: + """Read ``participants.tsv``, normalising and validating the group column.""" + try: + import pandas as pd + except ImportError as exc: # pragma: no cover + raise ImportError("reading NARPS metadata requires pandas") from exc + + table = pd.read_csv(path, sep="\t") + if "participant_id" not in table.columns or "group" not in table.columns: + raise NARPSFormatError( + "participants.tsv needs participant_id and group columns; found " + f"{list(table.columns)}" + ) + unknown = set(table["group"].dropna().unique()) - set(GROUPS) + if unknown: + logger.warning( + "NARPS: participants.tsv has unexpected group values %s (expected %s)", + sorted(unknown), + list(GROUPS), + ) + return table + + +# ------------------------------------------------------------- block building + + +class NARPSLoader: + """Assemble NARPS into framework blocks. + + The format/behaviour layer (this class's ``from_events_and_rois`` and + ``behavior_block``) has no neuroimaging dependencies. The BOLD -> ROI + extraction is in :meth:`load`, which reuses :class:`FMRILoader`. + """ + + name = "narps" + + def __init__( + self, + include_rt: bool = True, + include_expected_value: bool = True, + standardize_fmri: bool = False, + ) -> None: + self.include_rt = include_rt + self.include_expected_value = include_expected_value + self.standardize_fmri = standardize_fmri + + # -- behaviour: the economic comparator --------------------------------- + + def behavior_block( + self, + events: Any, + subject_id: str, + ) -> ModalityBlock: + """Build the behaviour block from parsed events. + + Features are the economic variables that drive gamble choice: gain, loss, + expected value, and (optionally) RT. On NARPS this is a *strong* + comparator by design, not a strawman: acceptance is largely a function of + gain and loss, so an economic model forecasts choice well. The neural + arm has to beat this, and on the aggregate arm it usually will not (see + ``docs/narps.md``). ``accept`` is never a feature. + """ + columns: List[np.ndarray] = [ + events["gain"].to_numpy(dtype=float), + events["loss"].to_numpy(dtype=float), + ] + names = ["gain", "loss"] + + if self.include_expected_value: + columns.append(events["expected_value"].to_numpy(dtype=float)) + names.append("expected_value") + # |EV| separates "clearly good/bad" gambles from ambiguous ones near + # indifference, which is where choice is hardest and most variable. + columns.append(np.abs(events["expected_value"].to_numpy(dtype=float))) + names.append("abs_expected_value") + + if self.include_rt: + columns.append(events["RT"].to_numpy(dtype=float)) + names.append("RT") + + return ModalityBlock( + name=BEHAVIOR, + X=np.column_stack(columns), + subject_ids=np.array([subject_id] * len(events)), + stimulus_ids=events["stimulus_id"].to_numpy(), + feature_names=names, + provenance={ + "loader": "NARPSLoader", + "modality": BEHAVIOR, + "source": "narps_events", + "features": names, + "note": ( + "economic comparator (gain/loss/EV); strong by design on " + "gambles, not a strawman. See docs/narps.md." + ), + }, + ) + + def from_events_and_rois( + self, + events: Any, + roi_features: np.ndarray, + subject_id: str, + roi_names: Optional[Sequence[str]] = None, + ) -> Dict[str, ModalityBlock]: + """Build fMRI + behaviour blocks from parsed events and extracted ROIs. + + This is the nilearn-free entry point: hand it a trial-by-ROI matrix you + extracted however you like (a GLM, a masker, precomputed betas) and it + assembles the blocks with correct keys. It is also what :meth:`load` + calls after doing the extraction. + + ``roi_features`` must be ``(n_trials, n_rois)`` and row-aligned with + ``events`` (same order, same length). + """ + from .fmri import FMRILoader + + roi_features = np.asarray(roi_features, dtype=float) + if roi_features.shape[0] != len(events): + raise ValueError( + f"roi_features has {roi_features.shape[0]} rows but events has " + f"{len(events)}; they must be row-aligned trial for trial" + ) + + fmri_block = FMRILoader(standardize=self.standardize_fmri).from_arrays( + X=roi_features, + subject_ids=[subject_id] * len(events), + stimulus_ids=events["stimulus_id"].to_numpy(), + feature_names=list(roi_names) if roi_names is not None else None, + source="narps_rois", + ) + return { + "fmri": fmri_block, + BEHAVIOR: self.behavior_block(events, subject_id), + } + + +# ------------------------------------------------------------ aggregate outcome + + +def acceptance_rate_by_gamble( + events_by_subject: Sequence[Any], + min_subjects: int = 5, +) -> Tuple[Dict[str, float], Dict[str, Any]]: + """Population acceptance rate per gamble, the NARPS aggregate outcome. + + This is a genuine aggregate choice: the fraction of people who accepted each + ``(gain, loss)`` gamble. It is keyed the way the blocks are, so it drops + straight into ``y_aggregate``. + + A gamble seen by only one or two subjects has a rate that is mostly noise, so + gambles below ``min_subjects`` are excluded and the count is reported. + + Returns ``(y_aggregate, provenance)``. + """ + import pandas as pd + + combined = pd.concat(list(events_by_subject), ignore_index=True) + grouped = combined.groupby("stimulus_id")["accept"] + counts = grouped.count() + rates = grouped.mean() + + kept = counts[counts >= min_subjects].index + y_aggregate = {str(k): float(rates[k]) for k in kept} + dropped = len(rates) - len(kept) + + if dropped: + logger.info( + "NARPS: %d/%d gambles seen by fewer than %d subjects; excluded from " + "the aggregate arm", + dropped, + len(rates), + min_subjects, + ) + + provenance = { + "outcome": "population_acceptance_rate", + "n_gambles": len(y_aggregate), + "n_gambles_dropped": int(dropped), + "min_subjects": min_subjects, + "caveat": ( + "on gambles the economic baseline (gain/loss) forecasts this rate " + "almost by construction; brain is not expected to beat behaviour on " + "this arm. See docs/narps.md." + ), + } + return y_aggregate, provenance From 263ce2011838ad3df644969e4d4e3872d28af016 Mon Sep 17 00:00:00 2001 From: Aayush Gandhi Date: Fri, 7 Aug 2026 17:16:52 -0400 Subject: [PATCH 3/5] Add NARPS BIDS load path and confound support to FMRILoader NARPSLoader.load discovers a subject's runs by globbing *_task-MGT_*_events.tsv, matches each to its preprocessed BOLD and confounds by run label (refusing to proceed on a mismatch, since a dropped run misaligns trials against BOLD), then delegates the BOLD -> ROI sphere extraction to FMRILoader and concatenates runs into one fMRI block plus the behaviour block. FMRILoader.load gains an optional per-run `confounds` argument passed to the NiftiSpheresMasker, so motion and physiological regressors are removed from the ROI timeseries. It warns when confounds are omitted rather than treating that as fine: for task fMRI, head motion correlates with events. Recorded in provenance as confounds_regressed. Requires nilearn; exercised by the fixture in the next commit. --- src/behavioral_decoding/io/fmri.py | 27 ++++- src/behavioral_decoding/io/narps.py | 148 ++++++++++++++++++++++++++++ 2 files changed, 174 insertions(+), 1 deletion(-) diff --git a/src/behavioral_decoding/io/fmri.py b/src/behavioral_decoding/io/fmri.py index f227051..a7c4d1a 100644 --- a/src/behavioral_decoding/io/fmri.py +++ b/src/behavioral_decoding/io/fmri.py @@ -22,9 +22,12 @@ import numpy as np +from ..utils.logging import get_logger from ..utils.progress import progress from .base import FMRI, BaseLoader, ModalityBlock +logger = get_logger(__name__) + # name -> (x, y, z) MNI centre in mm DEFAULT_ROIS: Dict[str, Tuple[float, float, float]] = { "NAcc_L": (-10.0, 12.0, -2.0), @@ -91,6 +94,7 @@ def load( onset_shift_s: float = 4.0, window_s: float = 4.0, mask_img: Optional[object] = None, + confounds: Optional[Sequence[object]] = None, ) -> ModalityBlock: """Extract trial-wise ROI means from NIfTI runs. @@ -112,6 +116,13 @@ def load( replaced then. window_s: Length of the averaging window after the shift. + confounds: + Optional per-run nuisance regressors (one entry per ``func_paths`` + item: a path to a confounds table, a DataFrame, or an array), passed + to the masker so motion and physiological components are regressed + out of the ROI timeseries. Strongly recommended for task fMRI, where + head motion correlates with events; omitting it is allowed and + warned about, not silently fine. Notes ----- @@ -143,11 +154,24 @@ def load( subj_out: List[object] = [] stim_out: List[object] = [] + if confounds is None: + logger.warning( + "FMRILoader.load: no confounds supplied. For task fMRI this is a " + "real omission: head motion correlates with events and with " + "individual differences. Pass the fMRIPrep confounds table." + ) + elif len(confounds) != len(func_paths): + raise ValueError( + f"confounds has {len(confounds)} entries but there are " + f"{len(func_paths)} runs; supply one per run" + ) + n_runs = len(func_paths) for run_idx, func_path in enumerate( progress(func_paths, desc="fMRI runs", total=n_runs) ): - ts = masker.fit_transform(func_path) # (n_volumes, n_rois) + run_confound = confounds[run_idx] if confounds is not None else None + ts = masker.fit_transform(func_path, confounds=run_confound) ev = events[run_idx] onsets = np.asarray(ev["onset"], dtype=float) stims = np.asarray(ev["stimulus_id"]) @@ -175,6 +199,7 @@ def load( onset_shift_s=onset_shift_s, window_s=window_s, extraction="peak_window_mean", + confounds_regressed=confounds is not None, caveat="not a GLM; replace with first-level modelling before reporting", ), ) diff --git a/src/behavioral_decoding/io/narps.py b/src/behavioral_decoding/io/narps.py index 972a202..bdc1980 100644 --- a/src/behavioral_decoding/io/narps.py +++ b/src/behavioral_decoding/io/narps.py @@ -310,6 +310,154 @@ def from_events_and_rois( BEHAVIOR: self.behavior_block(events, subject_id), } + # -- BIDS load: discover files and delegate extraction to FMRILoader ---- + + def load( + self, + subject_dir: str, + bold_dir: Optional[str] = None, + confounds_dir: Optional[str] = None, + rois: Optional[Dict[str, Tuple[float, float, float]]] = None, + radius_mm: float = 6.0, + onset_shift_s: float = 4.0, + window_s: float = 4.0, + space: str = "MNI152NLin2009cAsym", + ) -> Dict[str, ModalityBlock]: + """Load one subject from BIDS, extracting ROI features via FMRILoader. + + Parameters + ---------- + subject_dir: + The subject's directory holding ``func/*_events.tsv``. + bold_dir: + Directory holding the preprocessed BOLD, if separate from + ``subject_dir`` (e.g. an fMRIPrep derivatives tree). Defaults to + ``subject_dir``. + confounds_dir: + Directory holding ``*_desc-confounds_timeseries.tsv``. Defaults to + ``bold_dir``. Confounds are strongly recommended; see the class + docstring and ``docs/narps.md`` trap 4. + space: + The template label to match in preprocessed BOLD filenames. + + Returns + ------- + ``{"fmri": block, "behavior": block}`` with all runs concatenated. + + Notes + ----- + Requires nilearn (via :meth:`FMRILoader.load`). Discovers runs by + globbing events files, matches each to its BOLD and confounds by run + label, and refuses to proceed if the counts disagree, because a + silently-dropped run misaligns trials against BOLD. + """ + from .fmri import DEFAULT_ROIS, FMRILoader + + subject_path = Path(subject_dir) + bold_path = Path(bold_dir) if bold_dir else subject_path + confounds_path = Path(confounds_dir) if confounds_dir else bold_path + + subject_id = _subject_from_dir(subject_path) + + event_files = sorted((subject_path / "func").glob("*_task-MGT_*_events.tsv")) + if not event_files: + # Some trees keep events beside the subject dir rather than in func/. + event_files = sorted(subject_path.glob("**/*_task-MGT_*_events.tsv")) + if not event_files: + raise NARPSFormatError( + f"no *_task-MGT_*_events.tsv found under {subject_path}. NARPS " + "downloads nothing; point this at a local ds001734 tree." + ) + + parsed_runs: List[Any] = [] + func_paths: List[str] = [] + confounds: List[Optional[str]] = [] + run_subjects: List[str] = [] + + for ev_file in event_files: + run = _run_label(ev_file.name) + bold = _find_bold(bold_path, subject_id, run, space) + if bold is None: + logger.warning( + "NARPS: no preprocessed BOLD found for %s run %s in space %s; " + "skipping this run", + subject_id, + run, + space, + ) + continue + parsed = parse_events(ev_file) + parsed_runs.append(parsed) + func_paths.append(str(bold)) + confounds.append(_find_confounds(confounds_path, subject_id, run)) + run_subjects.append(subject_id) + + if not func_paths: + raise NARPSFormatError( + f"found events for {subject_id} but no matching BOLD in space " + f"{space!r} under {bold_path}. Check the derivatives path and space." + ) + + loader = FMRILoader( + rois=rois or dict(DEFAULT_ROIS), + radius_mm=radius_mm, + standardize=self.standardize_fmri, + ) + fmri_block = loader.load( + func_paths=func_paths, + events=parsed_runs, + subject_ids=run_subjects, + t_r=NARPS_TR, + onset_shift_s=onset_shift_s, + window_s=window_s, + confounds=confounds if any(c is not None for c in confounds) else None, + ) + + import pandas as pd + + all_events = pd.concat(parsed_runs, ignore_index=True) + return { + "fmri": fmri_block, + BEHAVIOR: self.behavior_block(all_events, subject_id), + } + + +def _subject_from_dir(path: Path) -> str: + match = re.search(r"(sub-[A-Za-z0-9]+)", path.name) + if match: + return match.group(1) + # Fall back to any sub-* under the directory. + for child in path.glob("sub-*"): + return child.name + raise NARPSFormatError(f"could not determine a subject id from {path}") + + +def _run_label(filename: str) -> str: + match = re.search(r"run-([A-Za-z0-9]+)", filename) + return match.group(1) if match else "01" + + +def _find_bold(root: Path, subject_id: str, run: str, space: str) -> Optional[Path]: + """Locate the preprocessed BOLD for a run, preferring the requested space.""" + patterns = [ + f"**/{subject_id}_task-MGT_run-{run}_space-{space}_desc-preproc_bold.nii.gz", + f"**/{subject_id}_task-MGT_run-{run}_space-{space}*_bold.nii.gz", + f"**/{subject_id}_task-MGT_run-{run}_bold.nii.gz", + f"**/{subject_id}_task-MGT_run-{run}*_bold.nii*", + ] + for pattern in patterns: + hits = sorted(root.glob(pattern)) + if hits: + return hits[0] + return None + + +def _find_confounds(root: Path, subject_id: str, run: str) -> Optional[str]: + hits = sorted( + root.glob(f"**/{subject_id}_task-MGT_run-{run}_desc-confounds_timeseries.tsv") + ) + return str(hits[0]) if hits else None + # ------------------------------------------------------------ aggregate outcome From 8fec86e7c2723cec8a07c0ab651069758b5e52cb Mon Sep 17 00:00:00 2001 From: Aayush Gandhi Date: Fri, 7 Aug 2026 17:19:57 -0400 Subject: [PATCH 4/5] Add NARPS fixture and tests tests/narps_fixture.py writes a miniature but format-real NARPS BIDS tree: participants.tsv with a group column, per-run events.tsv with the verified columns and response strings, a small 4-D preprocessed BOLD NIfTI with an affine placing all five default ROIs inside the volume, and a confounds table. An EV-scaled BOLD bump is planted at the NAcc and MPFC voxels so extraction can be checked; consistent with docs/narps.md, no brain-beats-behaviour claim is planted. tests/test_narps.py: 19 tests. The core (parsing, binarisation, NoResp drop, keys, EV, behaviour block, aggregate acceptance rate, format errors) runs everywhere; the BIDS .load tests importorskip nilearn/nibabel and validate real sphere extraction locally, including that the planted NAcc-EV signal is recovered (r>0.5) and that confounds are regressed. Skipped cleanly in a minimal CI env. --- tests/narps_fixture.py | 232 +++++++++++++++++++++++++++++++++++++++++ tests/test_narps.py | 226 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 458 insertions(+) create mode 100644 tests/narps_fixture.py create mode 100644 tests/test_narps.py diff --git a/tests/narps_fixture.py b/tests/narps_fixture.py new file mode 100644 index 0000000..892daeb --- /dev/null +++ b/tests/narps_fixture.py @@ -0,0 +1,232 @@ +"""Write a tiny NARPS-format BIDS tree to disk, for tests and the demo. + +Produces the real BIDS layout the loader discovers: ``participants.tsv`` with a +``group`` column, and per subject/run ``*_task-MGT_*_events.tsv`` (with the +verified NARPS columns and response strings), a small 4-D preprocessed BOLD +NIfTI, and a confounds table. + +The dimensions are miniature, not NARPS-real: a handful of trials per run and a +13x13x13 volume, sized so the whole fixture is a few MB and the sphere extraction +runs in a second. The BIDS *structure* and file *format* are real, which is the +point: the loader is exercised against the layout it will meet, not a mock. + +Signal is planted so the pipeline can be checked. NAcc and MPFC voxels carry a +BOLD bump proportional to each trial's expected value, and the accept choice +follows expected value plus noise. So fMRI ROI features should predict +accept/reject above chance. Consistent with ``docs/narps.md``, no +brain-beats-behaviour claim is planted: on gambles the economic variables +predict choice at least as well, and that is the honest story. + +Everything here is simulated. No number from it is a finding about brains. +""" + +from __future__ import annotations + +import csv +from pathlib import Path +from typing import Any, Dict, List, Optional + +import numpy as np + +from behavioral_decoding.io.fmri import DEFAULT_ROIS +from behavioral_decoding.io.narps import NARPS_TR, gamble_key + +# 8 mm isotropic voxels, origin placing MNI (0,0,0) at voxel 6. A 13-voxel span +# then covers roughly MNI -48..56 mm on each axis, which contains all five +# default ROIs (MPFC at y=46 is the tight one; it lands at voxel ~11.75). +VOXEL_SIZE_MM = 8.0 +VOLUME_SHAPE = (13, 13, 13) +_AFFINE = np.array( + [ + [VOXEL_SIZE_MM, 0.0, 0.0, -48.0], + [0.0, VOXEL_SIZE_MM, 0.0, -48.0], + [0.0, 0.0, VOXEL_SIZE_MM, -48.0], + [0.0, 0.0, 0.0, 1.0], + ] +) + +# A gamble grid drawn from the equalRange condition (gains and losses 5-20). +_GAINS = (8, 12, 16, 20) +_LOSSES = (6, 10, 14, 18) + + +def _mni_to_voxel(x: float, y: float, z: float) -> tuple: + inv = np.linalg.inv(_AFFINE) + vox = inv @ np.array([x, y, z, 1.0]) + return tuple(int(round(v)) for v in vox[:3]) + + +def write_narps_fixture( + root: str, + n_subjects: int = 6, + n_runs: int = 2, + trials_per_run: int = 12, + trial_spacing_s: float = 7.0, + seed: int = 0, +) -> Dict[str, Any]: + """Write a synthetic NARPS BIDS tree under ``root``. Returns ground truth.""" + root_path = Path(root) + root_path.mkdir(parents=True, exist_ok=True) + + rng = np.random.default_rng(seed) + + nacc_voxels = [_mni_to_voxel(*DEFAULT_ROIS[r]) for r in ("NAcc_L", "NAcc_R")] + mpfc_voxel = _mni_to_voxel(*DEFAULT_ROIS["MPFC"]) + signal_voxels = nacc_voxels + [mpfc_voxel] + + onsets = 4.0 + trial_spacing_s * np.arange(trials_per_run) + n_volumes = int(np.ceil((onsets[-1] + 12.0) / NARPS_TR)) + + participants: List[Dict[str, Any]] = [] + accept_by_gamble: Dict[str, List[int]] = {} + + for s in range(1, n_subjects + 1): + subject_id = f"sub-{s:03d}" + group = "equalRange" if s % 2 == 0 else "equalIndifference" + participants.append( + {"participant_id": subject_id, "group": group, "gender": "n/a", "age": "n/a"} + ) + + func_dir = root_path / subject_id / "func" + func_dir.mkdir(parents=True, exist_ok=True) + subject_bias = rng.normal(0.0, 0.4) + + for run in range(1, n_runs + 1): + srng = np.random.default_rng(seed * 10_000 + s * 100 + run) + gains = srng.choice(_GAINS, size=trials_per_run) + losses = srng.choice(_LOSSES, size=trials_per_run) + evs = 0.5 * (gains.astype(float) - losses.astype(float)) + + # BOLD volume: noise everywhere, plus an EV-scaled bump at the reward + # ROIs over each trial's peak window. + bold = srng.normal(0.0, 1.0, size=VOLUME_SHAPE + (n_volumes,)).astype(np.float32) + for onset, ev in zip(onsets, evs): + lo = int(round((onset + 4.0) / NARPS_TR)) + hi = int(round((onset + 8.0) / NARPS_TR)) + lo = max(0, min(lo, n_volumes - 1)) + hi = max(lo + 1, min(hi, n_volumes)) + for vx, vy, vz in signal_voxels: + bold[vx, vy, vz, lo:hi] += 1.5 * ev + + _write_nifti( + func_dir + / f"{subject_id}_task-MGT_run-{run:02d}" + f"_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz", + bold, + ) + + # Choice follows EV plus subject bias plus noise, calibrated to a + # spread of accept rates rather than all-accept or all-reject. + logits = 0.5 * evs + subject_bias + srng.normal(0.0, 1.0, size=trials_per_run) + accepts = (logits > 0).astype(int) + responses = [ + _response_string(a, srng.random()) for a in accepts + ] + rts = np.clip(srng.normal(1.8, 0.4, size=trials_per_run), 0.3, 3.5) + + events_rows = [] + for onset, gain, loss, resp, rt in zip(onsets, gains, losses, responses, rts): + events_rows.append( + { + "onset": round(float(onset), 3), + "duration": 4, + "gain": int(gain), + "loss": int(loss), + "RT": round(float(rt), 3), + "participant_response": resp, + } + ) + key = gamble_key(gain, loss) + accept_by_gamble.setdefault(key, []).append(int(resp in ("strongly_accept", "weakly_accept"))) + + _write_tsv( + func_dir / f"{subject_id}_task-MGT_run-{run:02d}_events.tsv", + events_rows, + ["onset", "duration", "gain", "loss", "RT", "participant_response"], + ) + + # A minimal confounds table: six motion parameters, no NaNs, so the + # confound-regression path is exercised without special-casing. + _write_confounds( + func_dir + / f"{subject_id}_task-MGT_run-{run:02d}_desc-confounds_timeseries.tsv", + n_volumes, + srng, + ) + + _write_tsv( + root_path / "participants.tsv", + participants, + ["participant_id", "group", "gender", "age"], + ) + + return { + "root": str(root_path), + "n_subjects": n_subjects, + "n_runs": n_runs, + "trials_per_run": trials_per_run, + "volume_shape": VOLUME_SHAPE, + "n_volumes": n_volumes, + "expected": ( + "fMRI ROI features should predict accept/reject above chance out of " + "fold; no brain-beats-behaviour claim is planted (gambles favour the " + "economic baseline)." + ), + } + + +def _response_string(accept: int, u: float) -> str: + if accept: + return "strongly_accept" if u < 0.5 else "weakly_accept" + return "strongly_reject" if u < 0.5 else "weakly_reject" + + +def _write_nifti(path: Path, data: np.ndarray) -> None: + import nibabel as nib + + img = nib.Nifti1Image(data, _AFFINE) + img.header.set_zooms(tuple([VOXEL_SIZE_MM] * 3) + (NARPS_TR,)) + nib.save(img, str(path)) + + +def _write_tsv(path: Path, rows: List[Dict[str, Any]], columns: List[str]) -> None: + with open(path, "w", newline="") as handle: + writer = csv.DictWriter(handle, fieldnames=columns, delimiter="\t") + writer.writeheader() + writer.writerows(rows) + + +def _write_confounds(path: Path, n_volumes: int, rng: np.random.Generator) -> None: + cols = ["trans_x", "trans_y", "trans_z", "rot_x", "rot_y", "rot_z"] + data = rng.normal(0.0, 0.05, size=(n_volumes, len(cols))) + rows = [{c: round(float(v), 6) for c, v in zip(cols, row)} for row in data] + _write_tsv(path, rows, cols) + + +def synthetic_events(n: int = 20, seed: int = 0) -> Optional[Any]: + """A parsed-events DataFrame with planted structure, no files written. + + Used by tests that exercise the nilearn-free layer without touching disk. + Returns None if pandas is unavailable. + """ + try: + import pandas as pd + except ImportError: + return None + + rng = np.random.default_rng(seed) + gains = rng.choice(_GAINS, size=n) + losses = rng.choice(_LOSSES, size=n) + evs = 0.5 * (gains.astype(float) - losses.astype(float)) + accepts = (0.5 * evs + rng.normal(0, 1, size=n) > 0).astype(int) + responses = [_response_string(a, rng.random()) for a in accepts] + return pd.DataFrame( + { + "onset": 4.0 + 7.0 * np.arange(n), + "duration": 4, + "gain": gains, + "loss": losses, + "RT": np.clip(rng.normal(1.8, 0.4, size=n), 0.3, 3.5), + "participant_response": responses, + } + ) diff --git a/tests/test_narps.py b/tests/test_narps.py new file mode 100644 index 0000000..9a6f2e6 --- /dev/null +++ b/tests/test_narps.py @@ -0,0 +1,226 @@ +"""NARPS loader tests. + +The core (events parsing, keys, behaviour, aggregate) needs no neuroimaging +stack and runs everywhere. The BIDS ``.load`` path needs nilearn and nibabel, so +those tests importorskip and are skipped cleanly in a minimal CI environment +while still validating the real extraction locally. +""" + +from __future__ import annotations + +import io as _io + +import numpy as np +import pytest + +from behavioral_decoding.io.narps import ( + ACCEPT_RESPONSES, + NARPS_TR, + NARPSFormatError, + NARPSLoader, + acceptance_rate_by_gamble, + expected_value, + gamble_key, + load_participants, + parse_events, +) + +pd = pytest.importorskip("pandas") + +from narps_fixture import synthetic_events, write_narps_fixture # noqa: E402 + +REAL_EVENTS_TSV = ( + "onset\tduration\tgain\tloss\tRT\tparticipant_response\n" + "4.071\t4\t14\t6\t2.388\tweakly_accept\n" + "11.834\t4\t34\t14\t2.289\tstrongly_accept\n" + "27.535\t4\t10\t10\t1.457\tweakly_reject\n" + "36.435\t4\t12\t19\t1.973\tstrongly_reject\n" + "43.935\t4\t20\t8\t0\tNoResp\n" +) + + +def _events(): + return pd.read_csv(_io.StringIO(REAL_EVENTS_TSV), sep="\t") + + +# ---------------------------------------------------------------- constants + + +def test_tr_is_one_second(): + assert NARPS_TR == 1.0 + + +# --------------------------------------------------------------- events parse + + +def test_parse_events_binarises_response_correctly(): + parsed = parse_events(_events()) + # NoResp dropped, so 4 rows remain. + assert len(parsed) == 4 + accepts = dict(zip(parsed["participant_response"], parsed["accept"])) + assert accepts["weakly_accept"] == 1 + assert accepts["strongly_accept"] == 1 + assert accepts["weakly_reject"] == 0 + assert accepts["strongly_reject"] == 0 + + +def test_noresp_is_dropped_not_imputed(): + """NoResp has RT=0; treating that as a fast response would poison RT.""" + parsed = parse_events(_events(), drop_no_response=True) + assert "NoResp" not in set(parsed["participant_response"]) + assert (parsed["RT"] > 0).all() + + +def test_noresp_can_be_kept_when_asked(): + parsed = parse_events(_events(), drop_no_response=False) + assert len(parsed) == 5 + + +def test_confidence_preserves_the_strong_weak_distinction(): + parsed = parse_events(_events()) + conf = dict(zip(parsed["participant_response"], parsed["confidence"])) + assert conf["strongly_accept"] == 1.0 + assert conf["weakly_accept"] == 0.5 + + +def test_expected_value_and_keys(): + parsed = parse_events(_events()) + row = parsed[parsed["gain"] == 14].iloc[0] + assert row["stimulus_id"] == "g14_l06" + assert row["expected_value"] == pytest.approx(4.0) # 0.5*(14-6) + assert expected_value(34, 14) == pytest.approx(10.0) + assert gamble_key(10, 10) == "g10_l10" + + +def test_missing_column_raises(): + bad = _events().drop(columns=["loss"]) + with pytest.raises(NARPSFormatError, match="missing columns"): + parse_events(bad) + + +def test_unknown_response_string_raises(): + bad = _events().copy() + bad.loc[0, "participant_response"] = "maybe" + with pytest.raises(NARPSFormatError, match="unrecognised participant_response"): + parse_events(bad) + + +# ---------------------------------------------------------------- behaviour + + +def test_behaviour_block_excludes_accept_and_uses_economic_features(): + parsed = parse_events(_events()) + block = NARPSLoader().behavior_block(parsed, "sub-001") + assert "accept" not in block.feature_names + assert set(block.feature_names) == {"gain", "loss", "expected_value", "abs_expected_value", "RT"} + assert block.n_trials == 4 + + +def test_behaviour_can_drop_rt_and_ev(): + parsed = parse_events(_events()) + block = NARPSLoader(include_rt=False, include_expected_value=False).behavior_block( + parsed, "sub-001" + ) + assert set(block.feature_names) == {"gain", "loss"} + + +def test_from_events_and_rois_row_alignment_enforced(): + parsed = parse_events(_events()) + with pytest.raises(ValueError, match="row-aligned"): + NARPSLoader().from_events_and_rois(parsed, np.zeros((3, 5)), "sub-001") + + +def test_from_events_and_rois_builds_both_blocks(): + parsed = parse_events(_events()) + rois = np.random.default_rng(0).normal(size=(4, 5)) + blocks = NARPSLoader().from_events_and_rois( + parsed, rois, "sub-001", roi_names=["NAcc_L", "NAcc_R", "MPFC", "AIns_L", "AIns_R"] + ) + assert set(blocks) == {"fmri", "behavior"} + assert blocks["fmri"].feature_names == ["NAcc_L", "NAcc_R", "MPFC", "AIns_L", "AIns_R"] + assert np.array_equal(blocks["fmri"].stimulus_ids, blocks["behavior"].stimulus_ids) + + +# ---------------------------------------------------------------- aggregate + + +def test_acceptance_rate_by_gamble(): + # Two subjects, same two gambles, known acceptance pattern. + a = pd.DataFrame( + { + "stimulus_id": ["g20_l06", "g06_l20"], + "accept": [1, 0], + } + ) + b = pd.DataFrame( + { + "stimulus_id": ["g20_l06", "g06_l20"], + "accept": [1, 1], + } + ) + rates, prov = acceptance_rate_by_gamble([a, b], min_subjects=1) + assert rates["g20_l06"] == pytest.approx(1.0) # both accepted + assert rates["g06_l20"] == pytest.approx(0.5) # one of two + assert prov["outcome"] == "population_acceptance_rate" + + +def test_acceptance_rate_excludes_rare_gambles(): + a = pd.DataFrame({"stimulus_id": ["g20_l06", "g06_l20"], "accept": [1, 0]}) + rates, prov = acceptance_rate_by_gamble([a], min_subjects=5) + assert rates == {} + assert prov["n_gambles_dropped"] == 2 + + +# ----------------------------------------------------- BIDS load (nilearn) + + +def test_load_discovers_runs_and_extracts_rois(tmp_path): + pytest.importorskip("nilearn") + pytest.importorskip("nibabel") + write_narps_fixture(str(tmp_path), n_subjects=2, n_runs=2, trials_per_run=10, seed=0) + + blocks = NARPSLoader().load(str(tmp_path / "sub-001")) + assert set(blocks) == {"fmri", "behavior"} + # 2 runs x 10 trials. + assert blocks["fmri"].n_trials == 20 + assert blocks["fmri"].feature_names == ["NAcc_L", "NAcc_R", "MPFC", "AIns_L", "AIns_R"] + assert blocks["fmri"].provenance["confounds_regressed"] is True + + +def test_load_recovers_planted_nacc_ev_signal(tmp_path): + """The fixture plants an EV-scaled bump at NAcc/MPFC; extraction must find it.""" + pytest.importorskip("nilearn") + pytest.importorskip("nibabel") + write_narps_fixture(str(tmp_path), n_subjects=2, n_runs=2, trials_per_run=12, seed=1) + + blocks = NARPSLoader().load(str(tmp_path / "sub-001")) + fm, beh = blocks["fmri"], blocks["behavior"] + ev = beh.X[:, beh.feature_names.index("expected_value")] + nacc = fm.X[:, fm.feature_names.index("NAcc_L")] + r = float(np.corrcoef(nacc, ev)[0, 1]) + assert r > 0.5, f"planted NAcc-EV signal not recovered (r={r:.3f})" + + +def test_load_raises_when_no_events(tmp_path): + (tmp_path / "sub-999" / "func").mkdir(parents=True) + with pytest.raises(NARPSFormatError, match="no .*events.tsv"): + NARPSLoader().load(str(tmp_path / "sub-999")) + + +# ----------------------------------------------------------- participants + + +def test_participants_group_column(tmp_path): + write_narps_fixture(str(tmp_path), n_subjects=4, n_runs=1, trials_per_run=8, seed=0) + table = load_participants(str(tmp_path / "participants.tsv")) + assert set(table["group"]) <= {"equalIndifference", "equalRange"} + assert len(table) == 4 + + +def test_synthetic_events_helper_is_parseable(): + events = synthetic_events(n=16, seed=0) + parsed = parse_events(events) + assert len(parsed) <= 16 + assert set(parsed["participant_response"]) <= set( + ACCEPT_RESPONSES + ("strongly_reject", "weakly_reject") + ) From dacc769733b89609879e70b79d43248deb1c0892 Mon Sep 17 00:00:00 2001 From: Aayush Gandhi Date: Fri, 7 Aug 2026 17:31:26 -0400 Subject: [PATCH 5/5] Add load_narps orchestration, demo, and CI load_narps assembles per-subject blocks into an individual-level MultimodalDataset: one row per trial, y_individual is accept vs reject, CV groups by subject. It is built by direct construction rather than the (subject, stimulus) join, because a subject sees each gamble more than once and that join requires unique keys; the fMRI and behaviour blocks are row-aligned by shared event ordering, which is what MultimodalDataset actually needs. Keeping the gamble as the stimulus key also makes the aggregate arm (population acceptance rate per gamble) work. Group filtering is supported and warns when pooling equalIndifference and equalRange for the aggregate arm, since they saw different gamble matrices. scripts/run_narps.py --demo runs the whole path on the synthetic BIDS fixture and checks the honest, recoverable claim: fMRI ROI features predict accept above chance out of fold and survive label permutation. It does NOT assert brain-beats-behaviour, because on gambles the economic baseline dominates. On the fixture: fMRI balanced accuracy 0.80 individually; on the aggregate arm behaviour (0.75) and brain (0.71) are close, behaviour ahead, exactly as docs/narps.md predicts. CI now installs .[dev,fmri] so the NARPS BIDS tests and demo run rather than skip. 119 tests passing. --- .github/workflows/ci.yml | 11 +- README.md | 26 ++-- scripts/run_narps.py | 130 +++++++++++++++++ src/behavioral_decoding/io/__init__.py | 2 + src/behavioral_decoding/io/narps.py | 187 ++++++++++++++++++++++++- tests/test_narps.py | 52 +++++++ 6 files changed, 398 insertions(+), 10 deletions(-) create mode 100644 scripts/run_narps.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 59caa08..6b1ec47 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -28,7 +28,9 @@ jobs: - name: Install run: | python -m pip install --upgrade pip - pip install -e ".[dev]" + # [fmri] pulls nilearn + nibabel so the NARPS BIDS load path and its + # demo are exercised in CI rather than skipped. + pip install -e ".[dev,fmri]" - name: Lint run: ruff check src tests scripts @@ -46,3 +48,10 @@ jobs: # Fails if the loader miswires stimulus keys or the asymmetry features. - name: DEAP-format recovery check run: python scripts/run_deap.py --demo + + # Same check on the real NARPS BIDS format (synthetic fixture, no + # download). Fails if the loader miswires events, ROI extraction, or + # confounds. Honest claim only: fMRI predicts accept above chance, no + # brain-beats-behaviour assertion (gambles favour the economic baseline). + - name: NARPS-format recovery check + run: python scripts/run_narps.py --demo diff --git a/README.md b/README.md index d4dee45..1c603c8 100644 --- a/README.md +++ b/README.md @@ -24,16 +24,26 @@ baselines that a claim like it requires. Full citations, with DOIs, in ## Status -Framework and validation harness, working end to end on synthetic data. The -**DEAP loader is built and tested** against the real file format (`io/deap.py`, -[`docs/deap.md`](docs/deap.md)); it needs only the licensed download to run on -real recordings, and DEAP's YouTube stimuli give it a real aggregate outcome via -view counts. The fMRI/MNE/OpenCV loader paths still have not been run against -real recordings; expect to fix things. See [`docs/design.md`](docs/design.md) -§12 for the full list of what is not built. +Framework and validation harness, working end to end on synthetic data. Two +real-dataset loaders are built and tested against their actual file formats: + +- **DEAP** (`io/deap.py`, [`docs/deap.md`](docs/deap.md)): EEG + peripheral + + behaviour; needs the licensed download; YouTube stimuli give it a real + aggregate outcome via view counts. +- **NARPS ds001734** (`io/narps.py`, [`docs/narps.md`](docs/narps.md)): the + fMRI mixed-gambles reward task; reuses `FMRILoader`'s NAcc/vmPFC/AIns sphere + extraction; public OpenNeuro download, no licence. NARPS is an + *individual-level* fMRI validation: on gambles the economic baseline dominates + the aggregate arm by construction, so it is the fMRI plumbing check, not a + brain-beats-behaviour demonstration. + +The MNE/OpenCV loader paths still have not been run against real recordings; +expect to fix things. See [`docs/design.md`](docs/design.md) §12 for the full +list of what is not built. ```bash -python scripts/run_deap.py --demo # whole DEAP path on a synthetic fixture, no download +python scripts/run_deap.py --demo # whole DEAP path on a synthetic fixture, no download +python scripts/run_narps.py --demo # whole NARPS BIDS path on a synthetic fixture (needs .[fmri]) ``` --- diff --git a/scripts/run_narps.py b/scripts/run_narps.py new file mode 100644 index 0000000..11b4db4 --- /dev/null +++ b/scripts/run_narps.py @@ -0,0 +1,130 @@ +#!/usr/bin/env python3 +"""Run the pipeline on NARPS (ds001734), or on a synthetic NARPS-format fixture. + +With a real NARPS download: + + python scripts/run_narps.py --root /path/to/ds001734 \ + --derivatives /path/to/fmriprep --group equalIndifference + +With no download, to see the whole path on the true BIDS format: + + python scripts/run_narps.py --demo + +``--demo`` writes a miniature but format-real NARPS BIDS tree (participants.tsv, +per-run events.tsv with the verified columns and response strings, small 4-D +NIfTIs, confounds) with an EV-scaled BOLD bump planted at the NAcc and MPFC +voxels. It then checks that the fMRI ROI features predict accept/reject above +chance out of fold. That is the honest, recoverable claim: consistent with +``docs/narps.md``, no brain-beats-behaviour headline is asserted, because on +gambles the economic baseline dominates the aggregate arm by construction. + +Requires nilearn and nibabel (``pip install '.[fmri]'``). +""" + +from __future__ import annotations + +import argparse +import sys +import tempfile +from pathlib import Path + +REPO_ROOT = Path(__file__).resolve().parents[1] +sys.path.insert(0, str(REPO_ROOT / "src")) +sys.path.insert(0, str(REPO_ROOT / "tests")) + +from behavioral_decoding.config import ExperimentConfig # noqa: E402 +from behavioral_decoding.evaluation.neuroforecast import ( # noqa: E402 + format_forecast_comparison, +) +from behavioral_decoding.io.narps import load_narps # noqa: E402 +from behavioral_decoding.pipelines.train import run_experiment # noqa: E402 + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--root", help="NARPS ds001734 BIDS root") + parser.add_argument("--derivatives", help="fMRIPrep derivatives root, if separate") + parser.add_argument("--group", choices=["equalIndifference", "equalRange"]) + parser.add_argument("--space", default="MNI152NLin2009cAsym") + parser.add_argument("--no-aggregate", action="store_true", + help="skip the behaviour-dominated aggregate arm") + parser.add_argument("--demo", action="store_true", + help="synthesise a NARPS-format fixture instead of loading a real one") + parser.add_argument("--output-dir", default="results") + args = parser.parse_args() + + if not args.demo and not args.root: + parser.error("give --root PATH for a real run, or --demo for the fixture") + + tmp = None + if args.demo: + from narps_fixture import write_narps_fixture + + tmp = tempfile.mkdtemp(prefix="narps_demo_") + write_narps_fixture(tmp, n_subjects=12, n_runs=2, trials_per_run=16, seed=0) + args.root = tmp + args.group = "equalRange" + print("=" * 78) + print("SYNTHETIC NARPS-FORMAT FIXTURE (simulated; not a finding about brains)") + print("=" * 78) + print("fMRI ROI features should predict accept/reject above chance out of fold.") + print("No brain-beats-behaviour claim is planted; gambles favour the economic") + print("baseline on the aggregate arm. See docs/narps.md.") + print() + + dataset = load_narps( + args.root, + derivatives=args.derivatives, + group=args.group, + space=args.space, + with_aggregate=not args.no_aggregate, + ) + print(dataset.describe()) + print() + + cfg = ExperimentConfig(name="narps", output_dir=args.output_dir) + cfg.data.modalities = dataset.modalities + record = run_experiment(dataset, cfg) + + if args.demo: + print() + print("=" * 78) + print("GROUND-TRUTH RECOVERY CHECK") + print("=" * 78) + per_mod = record["individual"]["per_modality_pooled"] + fmri_balacc = per_mod["fmri"]["balanced_accuracy"] + perm_p = record["individual"]["permutation_test"]["p_value"] + + arms = record.get("aggregate") or {} + if arms: + print(format_forecast_comparison(arms, "regression")) + print() + + ok_fmri = fmri_balacc > 0.55 + ok_perm = perm_p < 0.05 + print(" [{}] fMRI ROI features predict accept above chance " + "(balanced accuracy {:.3f})".format("PASS" if ok_fmri else "FAIL", fmri_balacc)) + print(" [{}] individual result survives label permutation " + "(p = {:.4f})".format("PASS" if ok_perm else "FAIL", perm_p)) + print() + + ok = ok_fmri and ok_perm + if ok: + print("The loader recovers the planted reward signal on the true NARPS BIDS") + print("format: NAcc/MPFC betas predict choice. fMRI plumbing validated. It") + print("says nothing about brains.") + else: + print("The reward signal did not recover; the loader is likely miswiring the") + print("events, the ROI extraction, or the confounds. Fix before real NARPS.") + + import shutil + + if tmp: + shutil.rmtree(tmp, ignore_errors=True) + return 0 if ok else 1 + + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/behavioral_decoding/io/__init__.py b/src/behavioral_decoding/io/__init__.py index 70c0447..dca2e72 100644 --- a/src/behavioral_decoding/io/__init__.py +++ b/src/behavioral_decoding/io/__init__.py @@ -31,6 +31,7 @@ NARPSLoader, acceptance_rate_by_gamble, gamble_key, + load_narps, parse_events, ) from .registry import get_loader, register_loader @@ -63,6 +64,7 @@ "gamble_key", "get_loader", "load_deap", + "load_narps", "parse_events", "register_loader", ] diff --git a/src/behavioral_decoding/io/narps.py b/src/behavioral_decoding/io/narps.py index bdc1980..bbb6f78 100644 --- a/src/behavioral_decoding/io/narps.py +++ b/src/behavioral_decoding/io/narps.py @@ -351,6 +351,36 @@ def load( label, and refuses to proceed if the counts disagree, because a silently-dropped run misaligns trials against BOLD. """ + blocks, _ = self.load_with_events( + subject_dir, + bold_dir=bold_dir, + confounds_dir=confounds_dir, + rois=rois, + radius_mm=radius_mm, + onset_shift_s=onset_shift_s, + window_s=window_s, + space=space, + ) + return blocks + + def load_with_events( + self, + subject_dir: str, + bold_dir: Optional[str] = None, + confounds_dir: Optional[str] = None, + rois: Optional[Dict[str, Tuple[float, float, float]]] = None, + radius_mm: float = 6.0, + onset_shift_s: float = 4.0, + window_s: float = 4.0, + space: str = "MNI152NLin2009cAsym", + ) -> Tuple[Dict[str, ModalityBlock], Any]: + """As :meth:`load`, but also return the concatenated parsed events. + + ``load_narps`` needs the events to build the accept labels and the + aggregate acceptance rate, row-aligned with the blocks. The blocks and + the returned events share row order (run-then-trial), so + ``events["accept"]`` lines up with the block rows. + """ from .fmri import DEFAULT_ROIS, FMRILoader subject_path = Path(subject_dir) @@ -416,10 +446,11 @@ def load( import pandas as pd all_events = pd.concat(parsed_runs, ignore_index=True) - return { + blocks = { "fmri": fmri_block, BEHAVIOR: self.behavior_block(all_events, subject_id), } + return blocks, all_events def _subject_from_dir(path: Path) -> str: @@ -509,3 +540,157 @@ def acceptance_rate_by_gamble( ), } return y_aggregate, provenance + + +# ------------------------------------------------------------- orchestration + + +def load_narps( + root: str, + derivatives: Optional[str] = None, + group: Optional[str] = None, + subjects: Optional[Sequence[str]] = None, + space: str = "MNI152NLin2009cAsym", + loader: Optional[NARPSLoader] = None, + with_aggregate: bool = True, + min_subjects_per_gamble: int = 5, + onset_shift_s: float = 4.0, + window_s: float = 4.0, +): + """Load NARPS into a :class:`MultimodalDataset`. + + Individual level: one row per trial, ``y_individual`` is accept (1) vs + reject (0), and cross-validation groups by subject. This is the primary, + honest use of NARPS: validating that NAcc/vmPFC/AIns betas predict choice. + + Aggregate level: because the stimulus key is the gamble ``(gain, loss)``, the + forecasting arm pools trials by gamble across subjects and can forecast the + population acceptance rate. On by default, but read ``docs/narps.md`` first: + the economic baseline dominates this arm by construction, so a + behaviour-beats-brain result here is expected and correct, not a failure. + + Parameters + ---------- + root: + BIDS root of ds001734 (holds ``participants.tsv`` and ``sub-*`` dirs). + derivatives: + fMRIPrep derivatives root, if BOLD lives there rather than in ``root``. + group: + Restrict to one condition (``"equalIndifference"`` or ``"equalRange"``). + The two groups saw different gamble matrices, so pooling them mixes + stimulus spaces; loading one group at a time is the safe default for the + aggregate arm, and a warning fires if you pool. + + Notes + ----- + Requires nilearn. The individual-level dataset is built by direct + construction rather than the ``(subject, stimulus)`` join, because a subject + sees each gamble more than once and that join requires unique keys. The + fMRI and behaviour blocks are row-aligned by shared event ordering, which is + what :class:`MultimodalDataset` actually needs. + """ + import numpy as np + + from .base import MultimodalDataset + + active = loader or NARPSLoader() + root_path = Path(root) + + participants_path = root_path / "participants.tsv" + group_of: Dict[str, str] = {} + if participants_path.exists(): + table = load_participants(str(participants_path)) + group_of = dict(zip(table["participant_id"], table["group"])) + + subject_dirs = sorted(p for p in root_path.glob("sub-*") if p.is_dir()) + if subjects is not None: + wanted = {s.lower() for s in subjects} + subject_dirs = [p for p in subject_dirs if p.name.lower() in wanted] + if group is not None: + if group not in GROUPS: + raise ValueError(f"group must be one of {GROUPS}; got {group!r}") + subject_dirs = [p for p in subject_dirs if group_of.get(p.name) == group] + if not subject_dirs: + raise ValueError(f"no subjects in group {group!r} under {root}") + elif with_aggregate and len(set(group_of.values())) > 1: + logger.warning( + "NARPS: pooling %s across groups for the aggregate arm mixes two " + "gamble matrices; pass group= to analyse within condition", + sorted(set(group_of.values())), + ) + + if not subject_dirs: + raise FileNotFoundError( + f"no sub-* directories found under {root}. NARPS downloads nothing; " + "point this at a local ds001734 tree." + ) + + fmri_blocks: List[ModalityBlock] = [] + behavior_blocks: List[ModalityBlock] = [] + events_by_subject: List[Any] = [] + + for subject_path in subject_dirs: + blocks, events = active.load_with_events( + str(subject_path), + bold_dir=derivatives, + confounds_dir=derivatives, + space=space, + onset_shift_s=onset_shift_s, + window_s=window_s, + ) + fmri_blocks.append(blocks["fmri"]) + behavior_blocks.append(blocks[BEHAVIOR]) + events_by_subject.append(events) + + merged = { + "fmri": _concat_blocks("fmri", fmri_blocks), + BEHAVIOR: _concat_blocks(BEHAVIOR, behavior_blocks), + } + y_individual = np.concatenate( + [ev["accept"].to_numpy(dtype=int) for ev in events_by_subject] + ) + subject_ids = merged["fmri"].subject_ids + stimulus_ids = merged["fmri"].stimulus_ids + + y_aggregate = None + aggregate_prov: Dict[str, Any] = {} + if with_aggregate: + y_aggregate, aggregate_prov = acceptance_rate_by_gamble( + events_by_subject, min_subjects=min_subjects_per_gamble + ) + + return MultimodalDataset( + blocks=merged, + y_individual=y_individual, + subject_ids=subject_ids, + stimulus_ids=stimulus_ids, + y_aggregate=y_aggregate, + metadata={ + "dataset": "NARPS", + "synthetic": False, + "group": group or "all", + "n_subjects": len(subject_dirs), + "outcome": "accept_vs_reject", + "aggregate": aggregate_prov or "disabled", + "citation": "Botvinik-Nezer et al. (2019), Sci Data 6, 106, " + "doi:10.1038/s41597-019-0113-7; dataset ds001734", + "note": ( + "individual-level fMRI validation; the aggregate arm is " + "behaviour-dominated by construction. See docs/narps.md." + ), + }, + ) + + +def _concat_blocks(name: str, blocks: Sequence[ModalityBlock]) -> ModalityBlock: + """Stack per-subject blocks into one, preserving keys and provenance.""" + import numpy as np + + return ModalityBlock( + name=name, + X=np.vstack([b.X for b in blocks]), + subject_ids=np.concatenate([b.subject_ids for b in blocks]), + stimulus_ids=np.concatenate([b.stimulus_ids for b in blocks]), + feature_names=list(blocks[0].feature_names or []), + provenance=dict(blocks[0].provenance), + ) diff --git a/tests/test_narps.py b/tests/test_narps.py index 9a6f2e6..e2bd179 100644 --- a/tests/test_narps.py +++ b/tests/test_narps.py @@ -21,6 +21,7 @@ acceptance_rate_by_gamble, expected_value, gamble_key, + load_narps, load_participants, parse_events, ) @@ -224,3 +225,54 @@ def test_synthetic_events_helper_is_parseable(): assert set(parsed["participant_response"]) <= set( ACCEPT_RESPONSES + ("strongly_reject", "weakly_reject") ) + + +# --------------------------------------------------------- load_narps (nilearn) + + +def test_load_narps_builds_individual_dataset(tmp_path): + pytest.importorskip("nilearn") + pytest.importorskip("nibabel") + write_narps_fixture(str(tmp_path), n_subjects=6, n_runs=2, trials_per_run=12, seed=0) + + ds = load_narps(str(tmp_path), group="equalRange") + assert set(ds.modalities) == {"fmri", "behavior"} + assert not ds.metadata["synthetic"] + assert ds.metadata["dataset"] == "NARPS" + assert ds.metadata["group"] == "equalRange" + # Binary accept outcome. + assert set(np.unique(ds.y_individual)) <= {0, 1} + # Aggregate keyed by gamble. + if ds.y_aggregate: + assert all(k.startswith("g") and "_l" in k for k in ds.y_aggregate) + + +def test_load_narps_group_filter_selects_subjects(tmp_path): + pytest.importorskip("nilearn") + pytest.importorskip("nibabel") + write_narps_fixture(str(tmp_path), n_subjects=6, n_runs=1, trials_per_run=10, seed=0) + + ei = load_narps(str(tmp_path), group="equalIndifference", with_aggregate=False) + er = load_narps(str(tmp_path), group="equalRange", with_aggregate=False) + # Fixture alternates groups by subject parity, so 3 and 3. + assert ei.n_subjects == 3 + assert er.n_subjects == 3 + assert set(ei.subject_ids).isdisjoint(set(er.subject_ids)) + + +def test_load_narps_recovers_reward_signal_individually(tmp_path): + """The planted NAcc/MPFC-EV signal should let fMRI predict accept above chance. + + This is the NARPS positive control: the honest, recoverable claim. It does + not assert brain beats behaviour, because on gambles it should not. + """ + pytest.importorskip("nilearn") + pytest.importorskip("nibabel") + from behavioral_decoding.evaluation.neuroforecast import cross_validate_ensemble + + write_narps_fixture(str(tmp_path), n_subjects=8, n_runs=2, trials_per_run=16, seed=0) + ds = load_narps(str(tmp_path), group="equalRange") + result = cross_validate_ensemble(ds, n_splits_outer=3, n_splits_inner=2, seed=0) + + fmri = result["per_modality_pooled"]["fmri"]["balanced_accuracy"] + assert fmri > 0.55, f"fMRI did not predict accept above chance (balacc={fmri:.3f})"