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/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. 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 f17af32..dca2e72 100644 --- a/src/behavioral_decoding/io/__init__.py +++ b/src/behavioral_decoding/io/__init__.py @@ -25,6 +25,15 @@ 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, + load_narps, + parse_events, +) from .registry import get_loader, register_loader __all__ = [ @@ -37,6 +46,7 @@ "FACE", "FMRI", "KNOWN_MODALITIES", + "NARPS_TR", "PERIPHERAL", "BaseLoader", "BehaviorLoader", @@ -47,8 +57,14 @@ "FMRILoader", "ModalityBlock", "MultimodalDataset", + "NARPSFormatError", + "NARPSLoader", + "acceptance_rate_by_gamble", "binarise_ratings", + "gamble_key", "get_loader", "load_deap", + "load_narps", + "parse_events", "register_loader", ] 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 new file mode 100644 index 0000000..bbb6f78 --- /dev/null +++ b/src/behavioral_decoding/io/narps.py @@ -0,0 +1,696 @@ +"""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), + } + + # -- 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. + """ + 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) + 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) + blocks = { + "fmri": fmri_block, + BEHAVIOR: self.behavior_block(all_events, subject_id), + } + return blocks, all_events + + +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 + + +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 + + +# ------------------------------------------------------------- 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/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..e2bd179 --- /dev/null +++ b/tests/test_narps.py @@ -0,0 +1,278 @@ +"""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_narps, + 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") + ) + + +# --------------------------------------------------------- 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})"