From e2766219e84536a9000e4e5fc56ef74009ef1cca Mon Sep 17 00:00:00 2001
From: Kevin Sitek <sitek.kevin@gmail.com>
Date: Wed, 23 Sep 2026 22:26:50 -0500
Subject: [PATCH 01/10] Add opt-in --filter-method iir and --reject-mode abs
 preprocessing

Some published FFR pipelines band-pass with zero-phase first-order
Butterworth filters (butter + filtfilt) and reject trials on absolute
amplitude (max|x| >= threshold) rather than peak-to-peak. Both are now
opt-in (defaults unchanged), plumbed through the CLI, nipype workflow
nodes and sidecar (RejectionMode). Tests compare the IIR output with
scipy butter+sosfiltfilt and check abs vs ptp on constructed epochs.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 CHANGES.rst                      |  10 ++
 docs/source/pipeline_details.rst |   7 +-
 ffrprep/ffrprep_cli.py           |  32 ++++++-
 ffrprep/preproc.py               |  97 +++++++++++++++++--
 ffrprep/tests/test_cli.py        |  53 +++++++++++
 ffrprep/tests/test_preproc.py    | 158 +++++++++++++++++++++++++++++++
 6 files changed, 346 insertions(+), 11 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 2b715b0..d668269 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -19,6 +19,16 @@ CLI
 - New flag ``--difference-pairs A:B [C:D …]`` to compute difference
   evokeds across explicit pairs. For 2-type datasets the difference
   is auto-emitted; for 3+ types this flag is required to opt in.
+- New flag ``--filter-method {fir,iir}`` (default ``fir``): ``iir`` applies
+  zero-phase (forward-backward) first-order Butterworth high-pass and
+  low-pass filters (12 dB/octave overall), matching the ``butter`` /
+  ``filtfilt`` band-pass used in some published FFR pipelines.
+- New flag ``--reject-mode {ptp,abs}`` (default ``ptp``): ``abs`` drops
+  epochs whose absolute amplitude reaches ``--reject-eeg`` at any sample
+  (``max|x| >= threshold``, drop reason ``ABS_AMP``) instead of MNE's
+  peak-to-peak criterion. The preprocessing sidecar records
+  ``RejectionMode`` (``peak-to-peak`` / ``absolute-amplitude``) next to
+  ``RejectionThresholds``.
 - New flag ``--no-report``: skip HTML report generation (preprocessing
   and analysis) while still writing all derivatives. Intended for bulk
   runs over many subjects, where report figures dominate runtime.
diff --git a/docs/source/pipeline_details.rst b/docs/source/pipeline_details.rst
index faea4e0..8d69e01 100644
--- a/docs/source/pipeline_details.rst
+++ b/docs/source/pipeline_details.rst
@@ -132,7 +132,9 @@ Preprocessing Workflow Nodes
 *Default Parameters:*
    - **High-pass:** 1.0 Hz (removes slow drifts, preserves FFR frequencies)
    - **Low-pass:** 40.0 Hz (removes EMG and high-frequency noise)
-   - **Filter design:** Zero-phase FIR with automatic transition bandwidth
+   - **Filter design:** Zero-phase FIR with automatic transition bandwidth.
+     ``--filter-method iir`` switches to zero-phase first-order
+     Butterworth high-pass + low-pass passes (12 dB/octave overall).
 
 *Outputs:* Filtered EEG data
 
@@ -158,7 +160,8 @@ amplitude-based rejection.
      around stimulus onset
    - **Baseline:** ``--baseline -0.2 0`` seconds (pre-stimulus)
    - **Rejection:** ``--reject-eeg 75e-6`` (75 µV peak-to-peak); pass
-     ``--no-auto-reject`` to disable
+     ``--no-auto-reject`` to disable, or ``--reject-mode abs`` to drop
+     epochs whose absolute amplitude reaches the threshold at any sample
 
 *Outputs:* Epoched EEG data and the post-rejection drop log.
 
diff --git a/ffrprep/ffrprep_cli.py b/ffrprep/ffrprep_cli.py
index 9bf318c..5861d30 100644
--- a/ffrprep/ffrprep_cli.py
+++ b/ffrprep/ffrprep_cli.py
@@ -221,6 +221,8 @@ def _make_preproc_payload(args_snap, deriv_snap, subject, task_label, run_label,
         "ref_channels": ref_channels,
         "high_pass": effective_l,
         "low_pass": effective_h,
+        "filter_method": args_snap.get("filter_method", "fir"),
+        "reject_mode": args_snap.get("reject_mode", "ptp"),
         "baseline": baseline,
         "tmin": args_snap.get("tmin"),
         "tmax": args_snap.get("tmax"),
@@ -308,6 +310,8 @@ def _make_concat_payload(args_snap, deriv_snap, subject, task_label, runs,
         "ref_channels": ref_channels,
         "high_pass": effective_l,
         "low_pass": effective_h,
+        "filter_method": args_snap.get("filter_method", "fir"),
+        "reject_mode": args_snap.get("reject_mode", "ptp"),
         "baseline": baseline,
         "tmin": args_snap.get("tmin"),
         "tmax": args_snap.get("tmax"),
@@ -524,6 +528,8 @@ def _build_preproc_workflow(payload):
     target.ref_channels = payload["ref_channels"]
     target.high_pass = payload["high_pass"]
     target.low_pass = payload["low_pass"]
+    target.filter_method = payload.get("filter_method", "fir")
+    target.reject_mode = payload.get("reject_mode", "ptp")
     target.baseline = payload["baseline"]
     target.tmin = payload["tmin"]
     target.tmax = payload["tmax"]
@@ -1663,11 +1669,35 @@ def get_parser():
         "--reject-eeg",
         type=float,
         help=(
-            "Peak-to-peak rejection threshold for EEG channels in Volts. "
+            "Rejection threshold for EEG channels in Volts (peak-to-peak "
+            "by default; see --reject-mode). "
             "Set to 0 to disable automatic rejection."
         ),
         default=75e-6,
     )
+    preproc_group.add_argument(
+        "--reject-mode",
+        choices=["ptp", "abs"],
+        default="ptp",
+        help=(
+            "How --reject-eeg is applied: 'ptp' = MNE peak-to-peak "
+            "criterion (default); 'abs' = drop epochs whose absolute "
+            "amplitude reaches the threshold at any sample (max|x| >= "
+            "threshold), e.g. --reject-mode abs --reject-eeg 35e-6."
+        ),
+    )
+    preproc_group.add_argument(
+        "--filter-method",
+        choices=["fir", "iir"],
+        default="fir",
+        help=(
+            "Band-pass filter design: 'fir' = MNE zero-phase FIR with "
+            "automatic transition bandwidth (default); 'iir' = zero-phase "
+            "first-order Butterworth high-pass then low-pass (12 dB/octave "
+            "overall, forward-backward), as used in some published FFR "
+            "pipelines."
+        ),
+    )
     preproc_group.add_argument(
         "--no-auto-reject",
         action="store_true",
diff --git a/ffrprep/preproc.py b/ffrprep/preproc.py
index 5382145..8c18f3a 100644
--- a/ffrprep/preproc.py
+++ b/ffrprep/preproc.py
@@ -441,7 +441,8 @@ def reference_raw_file(input_raw_path, ref_channels, output_dir, subject=None, t
 
 
 def filter_raw_file(
-    input_raw_path, high_pass, low_pass, output_dir, subject=None, task=None, session=None, run=None
+    input_raw_path, high_pass, low_pass, output_dir, subject=None, task=None, session=None, run=None,
+    filter_method="fir",
 ):
     """
     Load a raw FIF file, apply filtering, save filtered FIF. Returns path.
@@ -455,7 +456,9 @@ def filter_raw_file(
     from ffrprep.preproc import filter_data as _filter_data
     from ffrprep.preproc import save_raw_helper as _save_raw_helper
 
-    filtered = _filter_data(raw, high_pass=high_pass, low_pass=low_pass)
+    filtered = _filter_data(
+        raw, high_pass=high_pass, low_pass=low_pass, filter_method=filter_method,
+    )
 
     # Derive identifiers from filename as in reference_raw_file
     import re
@@ -625,7 +628,7 @@ def reference_data(eeg_data=None, ref_channels=None):
     return eeg_data
 
 
-def filter_data(eeg_data=None, high_pass=None, low_pass=None):
+def filter_data(eeg_data=None, high_pass=None, low_pass=None, filter_method="fir"):
     """
     Apply frequency filters to the provided EEG data object.
 
@@ -639,6 +642,13 @@ def filter_data(eeg_data=None, high_pass=None, low_pass=None):
     low_pass : float
         Upper passband edge (in Hertz).
         Default = None.
+    filter_method : {"fir", "iir"}
+        ``"fir"`` (default) uses MNE's zero-phase FIR design with automatic
+        transition bandwidths. ``"iir"`` applies zero-phase (forward-
+        backward) first-order Butterworth filters, a high-pass then a
+        low-pass, so the band-pass rolls off at 12 dB/octave overall.
+        This mirrors the ``butter``/``filtfilt`` pipeline used in some
+        published FFR analyses.
 
     Returns
     -------
@@ -664,7 +674,25 @@ def filter_data(eeg_data=None, high_pass=None, low_pass=None):
 
     # Apply filtering using MNE's filter method
     # l_freq is low-frequency cutoff (high-pass), h_freq is high-frequency cutoff (low-pass)  # noqa: E501
-    filtered_eeg.filter(l_freq=high_pass, h_freq=low_pass)
+    if filter_method == "fir":
+        filtered_eeg.filter(l_freq=high_pass, h_freq=low_pass)
+    elif filter_method == "iir":
+        # Fresh iir_params per call: MNE writes the designed filter back
+        # into the dict it is given.
+        if high_pass is not None:
+            filtered_eeg.filter(
+                l_freq=high_pass, h_freq=None, method="iir",
+                iir_params=dict(order=1, ftype="butter"), phase="zero",
+            )
+        if low_pass is not None:
+            filtered_eeg.filter(
+                l_freq=None, h_freq=low_pass, method="iir",
+                iir_params=dict(order=1, ftype="butter"), phase="zero",
+            )
+    else:
+        raise ValueError(
+            f"filter_method must be 'fir' or 'iir', got {filter_method!r}"
+        )
 
     # Return the filtered data
     return filtered_eeg
@@ -686,6 +714,7 @@ def epoch_data(
     run=None,
     trial_types=None,
     verbose=True,
+    reject_mode="ptp",
 ):
     """
     Epoch the provided EEG data object based on events.
@@ -706,6 +735,13 @@ def epoch_data(
         Start time before event (in seconds).
     tmax : float, optional
         End time after event (in seconds).
+    reject : dict, optional
+        Rejection thresholds in Volts, e.g. ``{"eeg": 75e-6}``.
+    reject_mode : {"ptp", "abs"}
+        How ``reject["eeg"]`` is interpreted. ``"ptp"`` (default) is MNE's
+        peak-to-peak criterion. ``"abs"`` drops any epoch whose absolute
+        amplitude reaches the threshold at any EEG sample (``max|x| >=
+        threshold``); dropped epochs are logged with reason ``ABS_AMP``.
     verbose : bool
         Whether to print verbose output.
 
@@ -917,6 +953,16 @@ def _coerce_int(val):
             )
         chosen_event_id = {k: chosen_event_id[k] for k in requested}
 
+    if reject_mode not in ("ptp", "abs"):
+        raise ValueError(
+            f"reject_mode must be 'ptp' or 'abs', got {reject_mode!r}"
+        )
+    abs_threshold = None
+    epochs_reject = reject
+    if reject_mode == "abs" and reject:
+        abs_threshold = float(reject["eeg"])
+        epochs_reject = None
+
     epoched_data = Epochs(
         eeg_data,
         events=events,
@@ -927,13 +973,27 @@ def _coerce_int(val):
         tmax=tmax,
         baseline=baseline_window,
         verbose=verbose,
-        reject=reject,
+        reject=epochs_reject,
     )
 
     # drop_bad() will apply rejection heuristics and return the filtered
     # Epochs object; keep this separate for clarity and easier debugging.
     epoched_data = epoched_data.drop_bad()
 
+    if abs_threshold is not None:
+        import numpy as np
+        from mne import pick_types
+
+        eeg_picks = pick_types(epoched_data.info, eeg=True)
+        if len(eeg_picks) and len(epoched_data):
+            amplitude = np.abs(epoched_data.get_data(picks=eeg_picks))
+            too_large = np.flatnonzero(amplitude.max(axis=(1, 2)) >= abs_threshold)
+            if too_large.size:
+                epoched_data.drop(too_large, reason="ABS_AMP", verbose=False)
+        # Read back by _save_one_preproc_epochs for the JSON sidecar
+        # (MNE's own ``reject`` attribute is None in this mode).
+        epoched_data.ffrprep_reject_abs = abs_threshold
+
     return epoched_data, (tmin, tmax)
 
 
@@ -1126,6 +1186,8 @@ def create_preprocessing_workflow(name="ffrprep_preproc", disk_backed=False):
                 "ref_channels",
                 "high_pass",
                 "low_pass",
+                "filter_method",
+                "reject_mode",
                 "baseline",
                 "tmin",
                 "tmax",
@@ -1199,6 +1261,7 @@ def create_preprocessing_workflow(name="ffrprep_preproc", disk_backed=False):
                     "task",
                     "session",
                     "run",
+                    "filter_method",
                 ],
                 output_names=["filtered_raw_path"],
                 function=filter_raw_file,
@@ -1208,7 +1271,7 @@ def create_preprocessing_workflow(name="ffrprep_preproc", disk_backed=False):
     else:
         filter_node = Node(
             Function(
-                input_names=["eeg_data", "high_pass", "low_pass"],
+                input_names=["eeg_data", "high_pass", "low_pass", "filter_method"],
                 output_names=["filtered_data"],
                 function=filter_data,
             ),
@@ -1244,6 +1307,7 @@ def create_preprocessing_workflow(name="ffrprep_preproc", disk_backed=False):
                     "task",
                     "run",
                     "trial_types",
+                    "reject_mode",
                 ],
                 output_names=["epochs", "time_window"],
                 function=epoch_data,
@@ -1268,6 +1332,7 @@ def create_preprocessing_workflow(name="ffrprep_preproc", disk_backed=False):
                     "task",
                     "run",
                     "trial_types",
+                    "reject_mode",
                 ],
                 output_names=["epochs", "time_window"],
                 function=epoch_data,
@@ -1339,6 +1404,7 @@ def create_preprocessing_workflow(name="ffrprep_preproc", disk_backed=False):
                     [
                         ("high_pass", "high_pass"),
                         ("low_pass", "low_pass"),
+                        ("filter_method", "filter_method"),
                         ("output_dir", "output_dir"),
                         ("sub_label", "subject"),
                         ("task_label", "task"),
@@ -1363,6 +1429,7 @@ def create_preprocessing_workflow(name="ffrprep_preproc", disk_backed=False):
                         ("tmin", "tmin"),
                         ("tmax", "tmax"),
                         ("reject", "reject"),
+                        ("reject_mode", "reject_mode"),
                         ("picks", "picks"),
                         ("on_missing", "on_missing"),
                         ("event_id", "event_id"),
@@ -1423,7 +1490,15 @@ def create_preprocessing_workflow(name="ffrprep_preproc", disk_backed=False):
                 (inputnode, reference_node, [("ref_channels", "ref_channels")]),
                 # Connect reference to filter
                 (reference_node, filter_node, [("referenced_data", "eeg_data")]),
-                (inputnode, filter_node, [("high_pass", "high_pass"), ("low_pass", "low_pass")]),
+                (
+                    inputnode,
+                    filter_node,
+                    [
+                        ("high_pass", "high_pass"),
+                        ("low_pass", "low_pass"),
+                        ("filter_method", "filter_method"),
+                    ],
+                ),
                 # Connect filter to epoch
                 (filter_node, epoch_node, [("filtered_data", "eeg_data")]),
                 (
@@ -1434,6 +1509,7 @@ def create_preprocessing_workflow(name="ffrprep_preproc", disk_backed=False):
                         ("tmin", "tmin"),
                         ("tmax", "tmax"),
                         ("reject", "reject"),
+                        ("reject_mode", "reject_mode"),
                         ("picks", "picks"),
                         ("on_missing", "on_missing"),
                         ("event_id", "event_id"),
@@ -1822,7 +1898,12 @@ def _save_one_preproc_epochs(
             ),
         },
     }
-    if reject_thresholds:
+    abs_reject = getattr(epochs, "ffrprep_reject_abs", None)
+    if abs_reject is not None:
+        sidecar["RejectionMode"] = "absolute-amplitude"
+        sidecar["RejectionThresholds"] = {"eeg": float(abs_reject)}
+    elif reject_thresholds:
+        sidecar["RejectionMode"] = "peak-to-peak"
         sidecar["RejectionThresholds"] = {
             k: float(v) for k, v in reject_thresholds.items()
         }
diff --git a/ffrprep/tests/test_cli.py b/ffrprep/tests/test_cli.py
index a9a360a..815a013 100644
--- a/ffrprep/tests/test_cli.py
+++ b/ffrprep/tests/test_cli.py
@@ -372,6 +372,59 @@ def test_run_ffrprep_both_stages(
         mock_analysis_workflow.run.assert_called_once()
 
 
+def test_paper_exact_preprocessing_flags_default_to_mne_behaviour():
+    parser = get_parser()
+    base = ["/path/to/bids", "/path/to/output", "participant"]
+    args = parser.parse_args(base)
+    assert args.filter_method == "fir"
+    assert args.reject_mode == "ptp"
+
+    args = parser.parse_args(
+        base + ["--filter-method", "iir", "--reject-mode", "abs", "--reject-eeg", "35e-6"]
+    )
+    assert args.filter_method == "iir"
+    assert args.reject_mode == "abs"
+    assert args.reject_eeg == pytest.approx(35e-6)
+
+
+@pytest.mark.parametrize("flag, bad", [("--filter-method", "butter"), ("--reject-mode", "mean")])
+def test_paper_exact_preprocessing_flags_reject_unknown_choices(flag, bad):
+    parser = get_parser()
+    with pytest.raises(SystemExit):
+        parser.parse_args(["/b", "/o", "participant", flag, bad])
+
+
+def test_workflow_inputs_carry_filter_method_and_reject_mode(tmp_path):
+    """The CLI payload reaches the nipype inputnode (both node layouts)."""
+    from ffrprep.ffrprep_cli import _build_preproc_workflow, _make_preproc_payload
+
+    args_snap = {
+        "bids_dir": str(tmp_path / "bids"), "tmin": -0.04, "tmax": 0.213,
+        "picks": None, "on_missing": "warn", "event_id": None,
+        "events_file": None, "save_each_node": False, "work_dir": None,
+        "filter_method": "iir", "reject_mode": "abs",
+    }
+    deriv_snap = {
+        "derivatives_root": str(tmp_path / "deriv"),
+        "preprocessing_dir": str(tmp_path / "deriv" / "ffrprep-preprocessing"),
+        "preprocessing_subject_dir": str(tmp_path / "deriv" / "p" / "sub-01" / "eeg"),
+        "analysis_dir": str(tmp_path / "deriv" / "ffrprep-analysis"),
+        "analysis_subject_dir": str(tmp_path / "deriv" / "a" / "sub-01"),
+    }
+    for save_each_node in (False, True):
+        args_snap["save_each_node"] = save_each_node
+        payload = _make_preproc_payload(
+            args_snap, deriv_snap, "01", "da", "01", [], 70.0, 2000.0,
+            (-0.04, 0.0), {"eeg": 35e-6},
+        )
+        assert payload["filter_method"] == "iir"
+        assert payload["reject_mode"] == "abs"
+        wf = _build_preproc_workflow(payload)
+        inputs = wf.get_node("inputnode").inputs
+        assert inputs.filter_method == "iir"
+        assert inputs.reject_mode == "abs"
+
+
 def test_no_report_flag_defaults_off_and_can_be_enabled():
     parser = get_parser()
     base = ["/path/to/bids", "/path/to/output", "participant"]
diff --git a/ffrprep/tests/test_preproc.py b/ffrprep/tests/test_preproc.py
index 635cfac..55ba571 100644
--- a/ffrprep/tests/test_preproc.py
+++ b/ffrprep/tests/test_preproc.py
@@ -590,6 +590,164 @@ def test_epoch_data_trial_types_unknown_raises(tmp_path):
         )
 
 
+def _single_channel_raw(sfreq=16384.0, seconds=8.0, seed=0):
+    import numpy as np
+    from mne import create_info
+    from mne.io import RawArray
+
+    rng = np.random.default_rng(seed)
+    n = int(sfreq * seconds)
+    t = np.arange(n) / sfreq
+    data = (
+        rng.normal(0, 1e-6, n)
+        + 5e-6 * np.sin(2 * np.pi * 100.0 * t)
+        + 3e-6 * np.sin(2 * np.pi * 20.0 * t)
+    )[np.newaxis, :]
+    return RawArray(data, create_info(["Cz"], sfreq, ch_types=["eeg"]), verbose=False)
+
+
+def test_filter_data_default_is_fir_and_matches_explicit_fir():
+    import numpy as np
+
+    raw = _single_channel_raw(sfreq=2000.0, seconds=4.0)
+    default = filter_data(raw, high_pass=70.0, low_pass=500.0)
+    explicit = filter_data(raw, high_pass=70.0, low_pass=500.0, filter_method="fir")
+    np.testing.assert_array_equal(default.get_data(), explicit.get_data())
+
+
+def test_filter_data_iir_matches_scipy_butter_filtfilt():
+    """iir = order-1 Butterworth HP then LP, each zero-phase (filtfilt)."""
+    import numpy as np
+    from scipy import signal
+
+    sfreq = 16384.0
+    raw = _single_channel_raw(sfreq=sfreq, seconds=6.0)
+    filtered = filter_data(raw, high_pass=70.0, low_pass=2000.0, filter_method="iir")
+
+    x = raw.get_data()[0]
+    hp = signal.butter(1, 70.0 / (sfreq / 2), "high", output="sos")
+    lp = signal.butter(1, 2000.0 / (sfreq / 2), "low", output="sos")
+    expected = signal.sosfiltfilt(lp, signal.sosfiltfilt(hp, x))
+
+    interior = slice(int(sfreq), -int(sfreq))  # skip edge-padding differences
+    np.testing.assert_allclose(
+        filtered.get_data()[0][interior], expected[interior], rtol=0, atol=1e-12,
+    )
+    # And it actually filtered: the 20 Hz component is strongly attenuated.
+    assert np.std(filtered.get_data()[0][interior]) < np.std(x[interior])
+
+
+def test_filter_data_iir_supports_single_sided_and_leaves_input_untouched():
+    import numpy as np
+
+    raw = _single_channel_raw(sfreq=2000.0, seconds=4.0)
+    before = raw.get_data().copy()
+    high_only = filter_data(raw, high_pass=70.0, low_pass=None, filter_method="iir")
+    low_only = filter_data(raw, high_pass=None, low_pass=200.0, filter_method="iir")
+    np.testing.assert_array_equal(raw.get_data(), before)
+    assert not np.allclose(high_only.get_data(), before)
+    assert not np.allclose(low_only.get_data(), before)
+
+
+def test_filter_data_unknown_method_raises():
+    raw = _single_channel_raw(sfreq=1000.0, seconds=2.0)
+    with pytest.raises(ValueError, match="filter_method"):
+        filter_data(raw, high_pass=1.0, low_pass=40.0, filter_method="butter")
+
+
+def _raw_and_events_for_rejection(tmp_path):
+    """Three-epoch recording: clean, spike (+45 uV), +/-20 uV square wave."""
+    import numpy as np
+    import pandas as pd
+    from mne import create_info
+    from mne.io import RawArray
+
+    sfreq = 1000.0
+    n = int(sfreq * 6)
+    rng = np.random.default_rng(1)
+    data = rng.normal(0, 0.5e-6, n)
+    # epoch 2 (onset 2.0 s): single +45 uV sample after stimulus onset
+    data[int(2.0 * sfreq) + 50] += 45e-6
+    # epoch 3 (onset 3.0 s): +/-20 uV square wave -> ptp = 40 uV, max|x| = 20 uV
+    seg = slice(int(3.0 * sfreq) + 10, int(3.0 * sfreq) + 110)
+    data[seg] += 20e-6 * np.sign(np.sin(np.linspace(0, 8 * np.pi, 100)))
+    raw = RawArray(data[np.newaxis, :], create_info(["Cz"], sfreq, ["eeg"]), verbose=False)
+    events_path = tmp_path / "events.tsv"
+    pd.DataFrame({
+        "onset": [1.0, 2.0, 3.0],
+        "duration": [0.0] * 3,
+        "trial_type": ["a", "a", "a"],
+    }).to_csv(events_path, sep="\t", index=False)
+    return raw, events_path
+
+
+def test_epoch_data_ptp_rejection_default_rejects_square_wave(tmp_path):
+    raw, events_path = _raw_and_events_for_rejection(tmp_path)
+    epochs, _ = epoch_data(
+        raw, baseline=[-0.05, 0.0], events_file=str(events_path),
+        tmin=-0.05, tmax=0.2, reject={"eeg": 35e-6}, verbose=False,
+    )
+    assert len(epochs) == 1  # spike (ptp 45) and square wave (ptp 40) both dropped
+
+
+def test_epoch_data_abs_rejection_keeps_epoch_whose_ptp_only_exceeds(tmp_path):
+    raw, events_path = _raw_and_events_for_rejection(tmp_path)
+    epochs, _ = epoch_data(
+        raw, baseline=[-0.05, 0.0], events_file=str(events_path),
+        tmin=-0.05, tmax=0.2, reject={"eeg": 35e-6}, reject_mode="abs",
+        verbose=False,
+    )
+    # Only the +45 uV spike reaches |x| >= 35 uV; the +/-20 uV wave stays.
+    assert len(epochs) == 2
+    reasons = [r for r in epochs.drop_log if r]
+    assert reasons == [("ABS_AMP",)]
+    assert epochs.ffrprep_reject_abs == pytest.approx(35e-6)
+
+
+def test_epoch_data_abs_mode_without_reject_is_a_noop(tmp_path):
+    raw, events_path = _raw_and_events_for_rejection(tmp_path)
+    epochs, _ = epoch_data(
+        raw, baseline=[-0.05, 0.0], events_file=str(events_path),
+        tmin=-0.05, tmax=0.2, reject=None, reject_mode="abs", verbose=False,
+    )
+    assert len(epochs) == 3
+
+
+def test_epoch_data_unknown_reject_mode_raises(tmp_path):
+    raw, events_path = _raw_and_events_for_rejection(tmp_path)
+    with pytest.raises(ValueError, match="reject_mode"):
+        epoch_data(
+            raw, baseline=[-0.05, 0.0], events_file=str(events_path),
+            tmin=-0.05, tmax=0.2, reject={"eeg": 35e-6}, reject_mode="mean",
+            verbose=False,
+        )
+
+
+@pytest.mark.parametrize(
+    "reject_mode, reject, expected_mode",
+    [
+        ("abs", {"eeg": 35e-6}, "absolute-amplitude"),
+        ("ptp", {"eeg": 75e-6}, "peak-to-peak"),
+    ],
+)
+def test_sidecar_records_rejection_mode(tmp_path, reject_mode, reject, expected_mode):
+    raw, events_path = _raw_and_events_for_rejection(tmp_path)
+    epochs, _ = epoch_data(
+        raw, baseline=[-0.05, 0.0], events_file=str(events_path),
+        tmin=-0.05, tmax=0.2, reject=reject, reject_mode=reject_mode, verbose=False,
+    )
+    bids_root = tmp_path / "bids"
+    bids_root.mkdir()
+    out = save_preprocessing_outputs(
+        epochs, bids_root, subject="01", task="da", run="01",
+    )
+    out = out if isinstance(out, list) else [out]
+    with open(out[0].with_suffix(".json")) as f:
+        sidecar = json.load(f)
+    assert sidecar["RejectionMode"] == expected_mode
+    assert sidecar["RejectionThresholds"] == {"eeg": pytest.approx(reject["eeg"])}
+
+
 def test_create_preprocessing_workflow(tmp_path):
     """Test the create_preprocessing_workflow function."""
     # Test that the workflow can be created without errors

From 92af1e10a36c8372ea354575c913b3fe50f331b6 Mon Sep 17 00:00:00 2001
From: Kevin Sitek <sitek.kevin@gmail.com>
Date: Wed, 23 Sep 2026 22:29:14 -0500
Subject: [PATCH 02/10] Add generic FFR harmonic-amplitude and
 stimulus-to-response metrics

harmonic_amplitudes computes zero-padded, band-averaged FFT amplitudes at
the harmonics of f0 (fundamental and summed upper harmonics, in uV) and
stim_to_resp_xcorr the max cross-correlation in xcorr 'coeff' form with an
optional lag window and Fisher z. load_wav_mono / resample_signal bring a
stimulus to the EEG rate. Defaults follow Whiteford et al. (2025).

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 CHANGES.rst                              |  11 ++
 docs/source/api_ref.rst                  |   4 +
 ffrprep/analysis.py                      | 169 +++++++++++++++++++++++
 ffrprep/tests/test_analysis_numerical.py | 164 ++++++++++++++++++++++
 4 files changed, 348 insertions(+)

diff --git a/CHANGES.rst b/CHANGES.rst
index d668269..4cdb084 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -143,6 +143,17 @@ Analysis
     Python loop of ``scipy.stats.pearsonr`` calls (same values and
     pair ordering; ~0.3 s vs ~9 min for 3000 epochs x 4147 samples).
   - ``compute_fft``: amplitude spectrum helper.
+  - ``harmonic_amplitudes``: FFT amplitude (zero-padded, single-sided,
+    microvolts) of a response window, averaged in a ``bin_hz``-wide
+    band around each harmonic of ``f0``; returns the per-harmonic values,
+    the fundamental, and the summed upper harmonics. Defaults follow the
+    /da/ measure of Whiteford et al. (2025) (60-180 ms, 100 Hz, 10
+    harmonics, 60 Hz bins).
+  - ``stim_to_resp_xcorr``: maximum stimulus-to-response correlation in
+    MATLAB ``xcorr(..., 'coeff')`` form (no mean removal) with an
+    optional lag window, returning ``(r, Fisher z, lag_ms)``.
+  - ``load_wav_mono`` / ``resample_signal``: WAV reader and polyphase
+    resampler used to bring a stimulus to the EEG sampling rate.
 
 - Analysis worker granularity changed from per-file to
   per-(task, run) group. ``_collect_analysis_groups`` stitches
diff --git a/docs/source/api_ref.rst b/docs/source/api_ref.rst
index 8364c19..bdfa560 100644
--- a/docs/source/api_ref.rst
+++ b/docs/source/api_ref.rst
@@ -62,6 +62,10 @@ Reference API
    corr_resp_to_resp
    response_consistency
    compute_fft
+   harmonic_amplitudes
+   stim_to_resp_xcorr
+   load_wav_mono
+   resample_signal
 
 
 .. _ref_reports:
diff --git a/ffrprep/analysis.py b/ffrprep/analysis.py
index c7a3b88..71d2d56 100644
--- a/ffrprep/analysis.py
+++ b/ffrprep/analysis.py
@@ -1060,3 +1060,172 @@ def compute_fft(signal_data, sfreq, fmin=None, fmax=None):
         amplitude = amplitude[mask]
 
     return freqs, amplitude
+
+
+def harmonic_amplitudes(
+    evoked,
+    f0=100.0,
+    n_harmonics=10,
+    bin_hz=60.0,
+    tmin=0.060,
+    tmax=0.180,
+    pad_pow=2,
+    ch_name=None,
+    unit_scale=1e6,
+):
+    """
+    Spectral amplitude of the fundamental and its harmonics in an FFR.
+
+    The steady-state portion of the response (``tmin``-``tmax``) is
+    Fourier transformed with zero padding to
+    ``2 ** (ceil(log2(L)) + pad_pow)`` points (``L`` = window length in
+    samples), converted to single-sided amplitude ``2 * |FFT| / L`` and,
+    for harmonic ``k`` (``k = 1 .. n_harmonics``, harmonic 1 = ``f0``),
+    averaged over all FFT points within ``k * f0 +/- bin_hz / 2``. The
+    defaults reproduce the /da/ (100 Hz F0) F0 / upper-harmonic measure of
+    Whiteford et al. (2025), which is computed on the *sum* of the two
+    polarity averages.
+
+    Parameters
+    ----------
+    evoked : mne.Evoked
+        Response to analyze (in Volts), e.g. the sum of the per-polarity
+        averages.
+    f0 : float, default=100.0
+        Fundamental frequency of the stimulus in Hz.
+    n_harmonics : int, default=10
+        Number of harmonics (including the fundamental) to evaluate.
+    bin_hz : float, default=60.0
+        Full width of the averaging band centred on each harmonic, in Hz.
+    tmin, tmax : float
+        Analysis window in seconds relative to the Evoked time axis
+        (inclusive at both ends).
+    pad_pow : int, default=2
+        Extra powers of two of zero padding beyond the next power of two.
+    ch_name : str, optional
+        Channel to analyze. Defaults to the first channel.
+    unit_scale : float, default=1e6
+        Multiplier applied to the data before the FFT (1e6 -> microvolts).
+
+    Returns
+    -------
+    result : dict
+        ``"harmonics"`` : ndarray of shape (n_harmonics,) with the band-
+        averaged amplitude of each harmonic; ``"f0"`` : amplitude of the
+        first harmonic; ``"upper_harmonics"`` : sum of harmonics 2 to
+        ``n_harmonics`` (0.0 if ``n_harmonics`` is 1).
+    """
+    sfreq = float(evoked.info["sfreq"])
+    if (n_harmonics * f0 + bin_hz / 2.0) >= sfreq / 2.0:
+        raise ValueError(
+            f"Harmonic band up to {n_harmonics * f0 + bin_hz / 2.0:g} Hz "
+            f"reaches the Nyquist frequency ({sfreq / 2.0:g} Hz)."
+        )
+
+    picks = 0 if ch_name is None else evoked.ch_names.index(ch_name)
+    segment = evoked.copy().crop(tmin=tmin, tmax=tmax).data[picks]
+    n_samples = segment.size
+    if n_samples < 2:
+        raise ValueError("Analysis window contains fewer than 2 samples.")
+
+    n_fft = 2 ** (int(np.ceil(np.log2(n_samples))) + int(pad_pow))
+    amplitude = 2.0 * np.abs(np.fft.rfft(segment * unit_scale, n_fft)) / n_samples
+    freqs = np.fft.rfftfreq(n_fft, d=1.0 / sfreq)
+
+    half_bin = bin_hz / 2.0 + 1e-9
+    harmonics = np.array([
+        amplitude[np.abs(freqs - k * f0) <= half_bin].mean()
+        for k in range(1, n_harmonics + 1)
+    ])
+    return {
+        "harmonics": harmonics,
+        "f0": float(harmonics[0]),
+        "upper_harmonics": float(harmonics[1:].sum()),
+    }
+
+
+def stim_to_resp_xcorr(stim, resp, sfreq, lag_range_ms=None):
+    """
+    Maximum stimulus-to-response cross-correlation (coefficient form).
+
+    Both inputs must share the sampling rate ``sfreq`` and be already
+    cropped to the analysis windows (e.g. stimulus vowel 50-170 ms vs.
+    response 60-180 ms). The cross-correlation is normalized like MATLAB's
+    ``xcorr(..., 'coeff')``: divided by ``sqrt(sum(stim**2) *
+    sum(resp**2))`` with **no** mean removal. The reported value is the
+    maximum over lags (optionally restricted to ``lag_range_ms``).
+
+    Parameters
+    ----------
+    stim, resp : array-like
+        1-D stimulus and response waveforms sampled at ``sfreq``.
+    sfreq : float
+        Sampling frequency in Hz.
+    lag_range_ms : tuple of (float, float), optional
+        Inclusive lag window in ms in which to search for the maximum.
+        Positive lag means the response is delayed relative to the
+        stimulus. ``None`` searches all lags.
+
+    Returns
+    -------
+    r : float
+        Maximum correlation coefficient within the searched lags.
+    z : float
+        Fisher r-to-z transform of ``r`` (``arctanh``).
+    lag_ms : float
+        Lag in ms at which the maximum occurs.
+    """
+    stim = np.asarray(stim, dtype=float).ravel()
+    resp = np.asarray(resp, dtype=float).ravel()
+    if stim.size == 0 or resp.size == 0:
+        raise ValueError("stim and resp must be non-empty.")
+
+    denom = np.sqrt(np.sum(stim ** 2) * np.sum(resp ** 2))
+    if denom == 0:
+        return float("nan"), float("nan"), float("nan")
+
+    corr = signal.correlate(resp, stim, mode="full", method="fft") / denom
+    lags_ms = signal.correlation_lags(resp.size, stim.size, mode="full") / sfreq * 1000.0
+
+    mask = np.ones(corr.shape, dtype=bool)
+    if lag_range_ms is not None:
+        lo, hi = lag_range_ms
+        mask = (lags_ms >= lo - 1e-9) & (lags_ms <= hi + 1e-9)
+        if not mask.any():
+            return float("nan"), float("nan"), float("nan")
+
+    masked_idx = np.flatnonzero(mask)
+    best = masked_idx[np.argmax(corr[masked_idx])]
+    r = float(corr[best])
+    with np.errstate(divide="ignore"):
+        z = float(np.arctanh(r))
+    return r, z, float(lags_ms[best])
+
+
+def load_wav_mono(path):
+    """Read a WAV file as a float, mono (channel-averaged) waveform.
+
+    Returns
+    -------
+    data : ndarray
+        1-D float waveform (arbitrary amplitude scale).
+    sfreq : float
+        Sampling rate of the file in Hz.
+    """
+    from scipy.io import wavfile
+
+    sfreq, data = wavfile.read(str(path))
+    data = np.asarray(data, dtype=float)
+    if data.ndim > 1:
+        data = data.mean(axis=1)
+    return data, float(sfreq)
+
+
+def resample_signal(data, sfreq_in, sfreq_out):
+    """Polyphase-resample a 1-D signal (no-op if the rates match)."""
+    from fractions import Fraction
+
+    if int(round(sfreq_in)) == int(round(sfreq_out)):
+        return np.asarray(data, dtype=float)
+    ratio = Fraction(sfreq_out / sfreq_in).limit_denominator(100000)
+    return signal.resample_poly(np.asarray(data, dtype=float), ratio.numerator, ratio.denominator)
diff --git a/ffrprep/tests/test_analysis_numerical.py b/ffrprep/tests/test_analysis_numerical.py
index 99f46f3..6c6ff47 100644
--- a/ffrprep/tests/test_analysis_numerical.py
+++ b/ffrprep/tests/test_analysis_numerical.py
@@ -214,3 +214,167 @@ def test_response_consistency_single_epoch_returns_nan_and_empty():
         mean_r, r_vals = response_consistency(epochs)
     assert r_vals.size == 0
     assert np.isnan(mean_r)
+
+
+# ---------------------------------------------------------------------------
+# harmonic_amplitudes / stim_to_resp_xcorr / wav + resampling helpers
+# ---------------------------------------------------------------------------
+
+def _harmonic_evoked(amps_uv, f0=100.0, sfreq=16384.0, tmin=-0.04, tmax=0.213, phase=0.3):
+    """Evoked (Volts) made of sines at k*f0 with given amplitudes in microvolts."""
+    n = int(round((tmax - tmin) * sfreq)) + 1
+    times = tmin + np.arange(n) / sfreq
+    data = np.zeros(n)
+    for k, amp in enumerate(amps_uv, start=1):
+        data += amp * 1e-6 * np.sin(2 * np.pi * k * f0 * times + phase * k)
+    info = mne.create_info(["Cz"], sfreq, "eeg")
+    return mne.EvokedArray(data[np.newaxis, :], info, tmin=tmin, nave=1, verbose="ERROR")
+
+
+def _reference_harmonic_amplitudes(evoked, f0, n_harmonics, bin_hz, tmin, tmax, pad_pow):
+    """Straight-line restatement of the published MATLAB recipe (dataFFT.m)."""
+    sfreq = evoked.info["sfreq"]
+    # Nearest-sample, inclusive window (MATLAB-style index lookup).
+    first = int(round((tmin - evoked.tmin) * sfreq))
+    last = int(round((tmax - evoked.tmin) * sfreq))
+    x = evoked.data[0][first:last + 1] * 1e6
+    length = x.size
+    nfft = 2 ** (int(np.ceil(np.log2(length))) + pad_pow)
+    amp = 2 * np.abs(np.fft.fft(x, nfft) / length)
+    freqs = np.arange(nfft) * sfreq / nfft
+    out = []
+    for k in range(1, n_harmonics + 1):
+        band = (freqs >= k * f0 - bin_hz / 2) & (freqs <= k * f0 + bin_hz / 2)
+        out.append(amp[band].mean())
+    return np.array(out)
+
+
+def test_harmonic_amplitudes_matches_reference_recipe():
+    from ffrprep.analysis import harmonic_amplitudes
+
+    evoked = _harmonic_evoked([1.0, 0.5, 0.25, 0.1, 0.05])
+    result = harmonic_amplitudes(evoked)
+    expected = _reference_harmonic_amplitudes(
+        evoked, f0=100.0, n_harmonics=10, bin_hz=60.0, tmin=0.06, tmax=0.18, pad_pow=2,
+    )
+    np.testing.assert_allclose(result["harmonics"], expected, rtol=1e-12)
+    assert result["f0"] == pytest.approx(expected[0])
+    assert result["upper_harmonics"] == pytest.approx(expected[1:].sum())
+
+
+def test_harmonic_amplitudes_scales_linearly_and_reports_microvolts():
+    from ffrprep.analysis import harmonic_amplitudes
+
+    base = harmonic_amplitudes(_harmonic_evoked([1.0, 0.5]))
+    doubled = harmonic_amplitudes(_harmonic_evoked([2.0, 1.0]))
+    assert doubled["f0"] == pytest.approx(2 * base["f0"])
+    assert doubled["upper_harmonics"] == pytest.approx(2 * base["upper_harmonics"])
+    # A 1 uV tone yields a microvolt-scale (not volt-scale) amplitude.
+    assert 1e-3 < base["f0"] < 1.0
+
+    single_uv = harmonic_amplitudes(_harmonic_evoked([1.0]))
+    single_volts = harmonic_amplitudes(_harmonic_evoked([1.0]), unit_scale=1.0)
+    assert single_volts["f0"] == pytest.approx(single_uv["f0"] * 1e-6)
+
+
+def test_harmonic_amplitudes_separates_fundamental_from_upper_harmonics():
+    from ffrprep.analysis import harmonic_amplitudes
+
+    fundamental_only = harmonic_amplitudes(_harmonic_evoked([1.0]))
+    harmonics_only = harmonic_amplitudes(_harmonic_evoked([0.0, 1.0, 1.0, 1.0]))
+    assert fundamental_only["f0"] > fundamental_only["upper_harmonics"] * 5
+    assert harmonics_only["upper_harmonics"] > harmonics_only["f0"] * 5
+
+
+def test_harmonic_amplitudes_validates_inputs():
+    from ffrprep.analysis import harmonic_amplitudes
+
+    low_rate = _harmonic_evoked([1.0], sfreq=1000.0)
+    with pytest.raises(ValueError, match="Nyquist"):
+        harmonic_amplitudes(low_rate)
+    with pytest.raises(ValueError, match="fewer than 2 samples"):
+        harmonic_amplitudes(_harmonic_evoked([1.0]), tmin=0.1, tmax=0.1)
+
+
+def test_stim_to_resp_xcorr_recovers_delay_and_matches_numpy_reference():
+    from ffrprep.analysis import stim_to_resp_xcorr
+
+    rng = np.random.default_rng(4)
+    sfreq = 16384.0
+    stim = rng.normal(size=1967)
+    delay = 164  # samples (~10 ms)
+    resp = np.concatenate([np.zeros(delay), stim])[: stim.size] + 0.05 * rng.normal(size=stim.size)
+
+    r, z, lag_ms = stim_to_resp_xcorr(stim, resp, sfreq)
+
+    full = np.correlate(resp, stim, mode="full") / np.sqrt(np.sum(stim**2) * np.sum(resp**2))
+    assert r == pytest.approx(full.max())
+    assert z == pytest.approx(np.arctanh(r))
+    assert lag_ms == pytest.approx(delay / sfreq * 1000.0)
+
+
+def test_stim_to_resp_xcorr_does_not_remove_the_mean():
+    """'coeff' normalization (MATLAB xcorr) keeps DC: not a Pearson r."""
+    from ffrprep.analysis import stim_to_resp_xcorr
+
+    stim = np.ones(200)
+    resp = np.ones(200)
+    r, _, lag_ms = stim_to_resp_xcorr(stim, resp, 1000.0)
+    assert r == pytest.approx(1.0)
+    assert lag_ms == pytest.approx(0.0)
+
+
+def test_stim_to_resp_xcorr_lag_range_restricts_search():
+    from ffrprep.analysis import stim_to_resp_xcorr
+
+    rng = np.random.default_rng(6)
+    sfreq = 16384.0
+    stim = rng.normal(size=1500)
+    resp = np.concatenate([np.zeros(164), stim])[:1500]  # true lag ~ 10 ms
+    r_all, _, lag_all = stim_to_resp_xcorr(stim, resp, sfreq)
+    r_win, _, lag_win = stim_to_resp_xcorr(stim, resp, sfreq, lag_range_ms=(0.0, 5.0))
+    assert lag_all == pytest.approx(164 / sfreq * 1000)
+    assert 0.0 <= lag_win <= 5.0 + 1e-6
+    assert r_win < r_all
+    # Window containing the true lag finds it again.
+    r_in, _, lag_in = stim_to_resp_xcorr(stim, resp, sfreq, lag_range_ms=(6.9, 10.9))
+    assert lag_in == pytest.approx(lag_all)
+    assert r_in == pytest.approx(r_all)
+
+
+def test_stim_to_resp_xcorr_degenerate_inputs():
+    from ffrprep.analysis import stim_to_resp_xcorr
+
+    r, z, lag = stim_to_resp_xcorr(np.zeros(50), np.ones(50), 1000.0)
+    assert np.isnan(r) and np.isnan(z) and np.isnan(lag)
+    with pytest.raises(ValueError):
+        stim_to_resp_xcorr([], [1.0], 1000.0)
+    r, _, _ = stim_to_resp_xcorr(np.ones(10), np.ones(10), 1000.0, lag_range_ms=(500.0, 600.0))
+    assert np.isnan(r)
+
+
+def test_load_wav_mono_and_resample_signal(tmp_path):
+    from scipy.io import wavfile
+
+    from ffrprep.analysis import load_wav_mono, resample_signal
+
+    sfreq_in = 24414
+    t = np.arange(4151) / sfreq_in
+    left = (np.sin(2 * np.pi * 100 * t) * 20000).astype(np.int16)
+    right = (np.sin(2 * np.pi * 100 * t) * 10000).astype(np.int16)
+    path = tmp_path / "stim.wav"
+    wavfile.write(path, sfreq_in, np.column_stack([left, right]))
+
+    data, sfreq = load_wav_mono(path)
+    assert sfreq == sfreq_in
+    assert data.ndim == 1
+    np.testing.assert_allclose(data, (left.astype(float) + right.astype(float)) / 2)
+
+    out = resample_signal(data, sfreq_in, 16384.0)
+    assert out.size == pytest.approx(data.size * 16384 / sfreq_in, abs=2)
+    same = resample_signal(data, 16384.0, 16384.0)
+    np.testing.assert_array_equal(same, np.asarray(data, dtype=float))
+    # A 100 Hz tone stays a 100 Hz tone after resampling.
+    spectrum = np.abs(np.fft.rfft(out))
+    peak_hz = np.fft.rfftfreq(out.size, 1 / 16384.0)[np.argmax(spectrum)]
+    assert peak_hz == pytest.approx(100.0, abs=3.0)

From deba371b45418185be51113ae7f833d64afd0952 Mon Sep 17 00:00:00 2001
From: Kevin Sitek <sitek.kevin@gmail.com>
Date: Wed, 23 Sep 2026 22:34:36 -0500
Subject: [PATCH 03/10] Add opt-in FFR metrics and covariates to the group
 metrics table

With --f0 the group table gains rms_snr_polarity_sum, f0_uv and
upper_harmonics_uv computed on the sum of the two per-trial-type averages;
--stimulus adds stimulus-to-response r/z/lag (and a lag-limited variant);
--n-trials-presented adds usable_pct. participants.tsv and --covariates are
joined onto the saved TSV. discover_group_inputs returns per-type evokeds,
and the report section shows histograms instead of an unreadably tall table
for large cohorts. Defaults leave existing group outputs unchanged.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 CHANGES.rst                 |  12 ++
 docs/source/api_ref.rst     |   1 +
 docs/source/usage.rst       |  18 +++
 ffrprep/ffrprep_cli.py      | 119 ++++++++++++++++
 ffrprep/group.py            | 266 ++++++++++++++++++++++++++++++++++--
 ffrprep/reports.py          |  85 ++++++++----
 ffrprep/tests/test_cli.py   |  35 +++++
 ffrprep/tests/test_group.py | 253 ++++++++++++++++++++++++++++++++++
 8 files changed, 746 insertions(+), 43 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 4cdb084..3413ef5 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -63,6 +63,18 @@ Group-level analysis
   - Deliberately aggregation-only, not inferential: no group-level
     statistics are computed, matching the scope other BIDS Apps
     (e.g. MRIQC) use for their own "group" level.
+- Group-level FFR metrics (opt-in): ``--f0`` adds ``rms_snr_polarity_sum``,
+  ``f0_uv`` and ``upper_harmonics_uv`` (band-averaged FFT amplitude at the
+  harmonics of F0, computed on the sum of the two per-trial-type
+  averages; see ``--n-harmonics``, ``--harmonic-bin-hz``,
+  ``--harmonic-window``); ``--stimulus`` adds ``stim2resp_r/z/lag_ms``
+  (and ``stim2resp_lim_*`` with ``--xcorr-lag-range``); ``--n-trials-presented``
+  adds ``usable_pct``. ``discover_group_inputs`` now also returns the
+  per-trial-type evoked files (``"by_type"``).
+- The group metrics TSV joins subject-level covariates from
+  ``<bids_dir>/participants.tsv`` and ``--covariates`` (``merge_covariates``);
+  the report section switches to per-metric histograms for cohorts above
+  40 subjects (the per-subject values stay in the TSV).
 - ``ffrprep.reports`` gains ``build_group_report`` and
   ``build_metrics_table_section``; ``build_subject_report`` /
   ``build_analysis_report`` pass new ``entity_label`` / ``meta_label``
diff --git a/docs/source/api_ref.rst b/docs/source/api_ref.rst
index bdfa560..2817923 100644
--- a/docs/source/api_ref.rst
+++ b/docs/source/api_ref.rst
@@ -112,6 +112,7 @@ Reference API
    discover_group_inputs
    compute_grand_average
    compute_subject_metrics
+   merge_covariates
    save_group_outputs
    save_group_metrics
    run_group_level
diff --git a/docs/source/usage.rst b/docs/source/usage.rst
index 0866486..64e433a 100644
--- a/docs/source/usage.rst
+++ b/docs/source/usage.rst
@@ -265,6 +265,24 @@ recomputed from each subject's saved evoked/epochs), and a single
 ``--task`` restrict which subjects/tasks are aggregated, same as at
 the participant level.
 
+Optional metrics for FFR studies are switched on with flags. With
+``--f0`` the metrics table gains ``f0_uv`` and ``upper_harmonics_uv`` (FFT
+amplitude, in microvolts, averaged in a 60 Hz band around each harmonic of
+the stimulus F0 over 60-180 ms) and ``rms_snr_polarity_sum``, all computed
+on the **sum** of the two per-trial-type averages. With ``--stimulus`` (a
+WAV file) it also gains ``stim2resp_r`` / ``stim2resp_z`` /
+``stim2resp_lag_ms``, the maximum stimulus-to-response cross-correlation
+(``xcorr`` "coeff" normalization, Fisher z). ``participants.tsv`` from the
+BIDS directory (plus any ``--covariates`` TSVs) is joined onto the saved
+table so downstream statistics have group, age, etc. alongside the metrics:
+
+.. code-block:: bash
+
+    ffrprep /data /data/derivatives group \
+      --f0 100 --n-harmonics 10 --harmonic-window 0.06 0.18 \
+      --stimulus /path/to/da.wav --xcorr-stim-window 0.05 0.17 \
+      --xcorr-resp-window 0.06 0.18 --n-trials-presented 6000
+
 This step only aggregates outputs participant-level ffrprep has
 already computed — it does not perform any group-level statistics
 (no hypothesis tests, no GLM). Downstream statistical analysis is
diff --git a/ffrprep/ffrprep_cli.py b/ffrprep/ffrprep_cli.py
index 5861d30..4fda5e2 100644
--- a/ffrprep/ffrprep_cli.py
+++ b/ffrprep/ffrprep_cli.py
@@ -1807,6 +1807,125 @@ def get_parser():
         default=[0.100, 0.200],
     )
 
+    # Group-level options
+    group_group = parser.add_argument_group(
+        "Group-level options",
+        "Used with the 'group' analysis level. They add optional columns to "
+        "the group metrics table; nothing here performs statistical inference.",
+    )
+    group_group.add_argument(
+        "--f0",
+        type=float,
+        default=None,
+        help=(
+            "Stimulus fundamental frequency in Hz. When set, the group metrics "
+            "table gains rms_snr_polarity_sum, f0_uv and upper_harmonics_uv "
+            "(band-averaged FFT amplitude at each harmonic, in microvolts) "
+            "computed on the SUM of the two per-trial-type averages."
+        ),
+    )
+    group_group.add_argument(
+        "--n-harmonics",
+        dest="n_harmonics",
+        type=int,
+        default=10,
+        help="Number of harmonics (including the fundamental) for --f0.",
+    )
+    group_group.add_argument(
+        "--harmonic-bin-hz",
+        dest="harmonic_bin_hz",
+        type=float,
+        default=60.0,
+        help="Full width in Hz of the band averaged around each harmonic.",
+    )
+    group_group.add_argument(
+        "--harmonic-window",
+        dest="harmonic_window",
+        nargs=2,
+        type=float,
+        metavar=("START", "END"),
+        default=[0.060, 0.180],
+        help="Response window in seconds for the harmonic FFT.",
+    )
+    group_group.add_argument(
+        "--stimulus",
+        type=Path,
+        default=None,
+        help=(
+            "WAV file of the stimulus (t=0 at stimulus onset). When set, the "
+            "group metrics table gains stim2resp_r / stim2resp_z / "
+            "stim2resp_lag_ms: the maximum stimulus-to-response "
+            "cross-correlation ('coeff' normalization, Fisher z) on the "
+            "polarity-summed response. The stimulus is resampled to the EEG "
+            "sampling rate."
+        ),
+    )
+    group_group.add_argument(
+        "--xcorr-stim-window",
+        dest="xcorr_stim_window",
+        nargs=2,
+        type=float,
+        metavar=("START", "END"),
+        default=[0.050, 0.170],
+        help="Stimulus segment in seconds for --stimulus.",
+    )
+    group_group.add_argument(
+        "--xcorr-resp-window",
+        dest="xcorr_resp_window",
+        nargs=2,
+        type=float,
+        metavar=("START", "END"),
+        default=[0.060, 0.180],
+        help="Response segment in seconds for --stimulus.",
+    )
+    group_group.add_argument(
+        "--xcorr-lag-range",
+        dest="xcorr_lag_range",
+        nargs=2,
+        type=float,
+        metavar=("MIN_MS", "MAX_MS"),
+        default=None,
+        help=(
+            "Additionally report the maximum correlation restricted to this "
+            "lag window in ms (columns stim2resp_lim_*). Positive lag = "
+            "response delayed relative to the stimulus, with both segments "
+            "aligned at their start."
+        ),
+    )
+    group_group.add_argument(
+        "--xcorr-lag-resp-window",
+        dest="xcorr_lag_resp_window",
+        nargs=2,
+        type=float,
+        metavar=("START", "END"),
+        default=None,
+        help=(
+            "Response segment in seconds for the lag-restricted correlation "
+            "(defaults to --xcorr-resp-window)."
+        ),
+    )
+    group_group.add_argument(
+        "--n-trials-presented",
+        dest="n_trials_presented",
+        type=int,
+        default=None,
+        help=(
+            "Trials presented per recording; adds usable_pct (kept trials / "
+            "presented) to the group metrics table."
+        ),
+    )
+    group_group.add_argument(
+        "--covariates",
+        nargs="+",
+        type=Path,
+        default=None,
+        help=(
+            "Extra subject-level TSV files (with a participant_id column) to "
+            "join onto the group metrics table, in addition to "
+            "<bids_dir>/participants.tsv."
+        ),
+    )
+
     # General options
     parser.add_argument(
         "--skip_bids_validation",
diff --git a/ffrprep/group.py b/ffrprep/group.py
index 2ceaede..c8a7007 100644
--- a/ffrprep/group.py
+++ b/ffrprep/group.py
@@ -20,7 +20,15 @@
 import pandas as pd
 
 import ffrprep.reports as reports
-from ffrprep.analysis import compute_power, response_consistency, rms_snr
+from ffrprep.analysis import (
+    compute_power,
+    harmonic_amplitudes,
+    load_wav_mono,
+    resample_signal,
+    response_consistency,
+    rms_snr,
+    stim_to_resp_xcorr,
+)
 
 # Matches combined/diff evoked filenames written by
 # preproc.save_analysis_outputs(), e.g.
@@ -34,6 +42,7 @@
     r"_desc-(?P<desc>[A-Za-z0-9]+)\.fif$"
 )
 _DIFF_DESC_RE = re.compile(r"^evokedDiff([A-Za-z0-9]+)Vs([A-Za-z0-9]+)$")
+_PER_TYPE_DESC_RE = re.compile(r"^evoked(?P<label>[A-Za-z0-9]+)$")
 
 # Matches preprocessing epochs filenames written by
 # preproc.save_preprocessing_outputs(), e.g.
@@ -67,13 +76,16 @@ def discover_group_inputs(derivatives_root, participant_label=None, task=None):
     Returns
     -------
     inputs : dict
-        Dictionary with keys ``"evoked"``, ``"diff"``, and
-        ``"epochs"``. ``"evoked"`` and ``"epochs"`` map a
+        Dictionary with keys ``"evoked"``, ``"diff"``, ``"by_type"``
+        and ``"epochs"``. ``"evoked"`` and ``"epochs"`` map a
         ``(task, session, run)`` tuple to ``{subject: Path}``
         (``{subject: [Path, ...]}`` for ``"epochs"``); ``"diff"``
         maps a ``(task, session, run, condition_a, condition_b)``
-        tuple to ``{subject: Path}``. ``session`` and ``run`` are
-        None when not present in a filename.
+        tuple to ``{subject: Path}``; ``"by_type"`` maps a
+        ``(task, session, run)`` tuple to
+        ``{subject: {condition_label: Path}}`` for the per-trial-type
+        evoked files. ``session`` and ``run`` are None when not present
+        in a filename.
     """
     derivatives_root = Path(derivatives_root)
     analysis_dir = derivatives_root / "ffrprep-analysis"
@@ -81,6 +93,7 @@ def discover_group_inputs(derivatives_root, participant_label=None, task=None):
 
     evoked = defaultdict(dict)
     diff = defaultdict(dict)
+    by_type = defaultdict(lambda: defaultdict(dict))
     for path in sorted(analysis_dir.glob("sub-*/*.fif")):
         match = _EVOKED_RE.match(path.name)
         if not match:
@@ -101,6 +114,10 @@ def discover_group_inputs(derivatives_root, participant_label=None, task=None):
         if diff_match:
             a, b = diff_match.groups()
             diff[key + (a, b)][subject] = path
+            continue
+        type_match = _PER_TYPE_DESC_RE.match(desc)
+        if type_match:
+            by_type[key][subject][type_match.group("label")] = path
 
     epochs = defaultdict(lambda: defaultdict(list))
     for path in sorted(preproc_dir.glob("sub-*/eeg/*.fif")):
@@ -120,6 +137,10 @@ def discover_group_inputs(derivatives_root, participant_label=None, task=None):
     return {
         "evoked": dict(evoked),
         "diff": dict(diff),
+        "by_type": {
+            key: {subject: dict(conds) for subject, conds in subjects.items()}
+            for key, subjects in by_type.items()
+        },
         "epochs": {key: dict(value) for key, value in epochs.items()},
     }
 
@@ -164,7 +185,49 @@ def compute_grand_average(evoked_paths):
     return grand_average, subjects
 
 
-def compute_subject_metrics(evoked_paths, epochs_paths=None, response_window=(0.100, 0.200)):
+def _load_polarity_sum(cond_paths):
+    """Sum of a subject's two per-trial-type Evokeds (or None).
+
+    FFR analyses that suppress the cochlear microphonic / stimulus
+    artifact add the responses to the two stimulus polarities (a **sum**
+    of the per-polarity averages, not a trial-weighted mean). Returns
+    None unless exactly two per-type files are available.
+    """
+    if not cond_paths or len(cond_paths) != 2:
+        return None
+    first, second = (_read_evoked_with_baseline(cond_paths[c]) for c in sorted(cond_paths))
+    summed = mne.combine_evoked([first, second], weights=[1, 1])
+    summed.baseline = first.baseline
+    return summed
+
+
+def _window(array, sfreq, tmin, tmax):
+    """Nearest-sample, inclusive ``tmin``-``tmax`` slice of a signal that starts at t=0."""
+    return array[int(round(tmin * sfreq)): int(round(tmax * sfreq)) + 1]
+
+
+def _stim_to_resp_columns(summed, stim, spec, prefix="stim2resp"):
+    """r / z / lag columns for one stimulus-to-response window definition."""
+    sfreq = float(summed.info["sfreq"])
+    stim_tmin, stim_tmax = spec["stim_window"]
+    resp_tmin, resp_tmax = spec["resp_window"]
+    stim_seg = _window(stim, sfreq, stim_tmin, stim_tmax)
+    resp_seg = summed.copy().crop(tmin=resp_tmin, tmax=resp_tmax).data[0]
+    r, z, lag_ms = stim_to_resp_xcorr(
+        stim_seg, resp_seg, sfreq, lag_range_ms=spec.get("lag_range_ms"),
+    )
+    return {f"{prefix}_r": r, f"{prefix}_z": z, f"{prefix}_lag_ms": lag_ms}
+
+
+def compute_subject_metrics(
+    evoked_paths,
+    epochs_paths=None,
+    response_window=(0.100, 0.200),
+    by_type_paths=None,
+    harmonics=None,
+    stimulus=None,
+    n_trials_presented=None,
+):
     """
     Compute per-subject scalar FFR metrics from saved derivatives.
 
@@ -175,6 +238,10 @@ def compute_subject_metrics(evoked_paths, epochs_paths=None, response_window=(0.
     Nothing here is computed from raw data; this only aggregates
     already-computed participant-level derivatives.
 
+    Optionally (``harmonics`` / ``stimulus``) it also derives spectral and
+    stimulus-following measures from the **sum of the two per-trial-type
+    (polarity) averages**, the response commonly analyzed in FFR studies.
+
     Parameters
     ----------
     evoked_paths : dict[str, pathlib.Path]
@@ -187,27 +254,49 @@ def compute_subject_metrics(evoked_paths, epochs_paths=None, response_window=(0.
     response_window : tuple of (float, float)
         Response window for RMS SNR / band power, matching the
         participant-level ``--response-window`` default.
+    by_type_paths : dict[str, dict[str, pathlib.Path]], optional
+        ``{subject: {condition_label: per-type Evoked path}}``. Required
+        for the polarity-sum measures; a subject needs exactly two
+        conditions, otherwise those columns are NaN.
+    harmonics : dict, optional
+        Keyword arguments for :func:`ffrprep.analysis.harmonic_amplitudes`
+        (``f0``, ``n_harmonics``, ``bin_hz``, ``tmin``, ``tmax``). Adds
+        ``rms_snr_polarity_sum``, ``f0_uv`` and ``upper_harmonics_uv``.
+    stimulus : dict, optional
+        ``{"path": wav, "stim_window": (t0, t1), "resp_window": (t0, t1)}``
+        plus optionally ``"lag_range_ms": (lo, hi)`` and
+        ``"lag_resp_window": (t0, t1)``. Adds ``stim2resp_r``,
+        ``stim2resp_z``, ``stim2resp_lag_ms`` (and ``stim2resp_lim_*`` when
+        a lag range is given). The stimulus is resampled to the EEG rate.
+    n_trials_presented : int, optional
+        Trials presented per recording; when given, ``usable_pct`` =
+        100 * kept trials / presented is added.
 
     Returns
     -------
     metrics : pandas.DataFrame
         One row per subject with columns ``subject``, ``n_avg``,
-        ``rms_snr``, ``band_power_90_110hz``, ``response_consistency``.
+        ``rms_snr``, ``band_power_90_110hz``, ``response_consistency`` and
+        the optional columns described above.
     """
     epochs_paths = epochs_paths or {}
+    by_type_paths = by_type_paths or {}
     resp_lower, resp_upper = response_window
+    use_sum = harmonics is not None or stimulus is not None
+    stim_cache = {}
     rows = []
     for subject in sorted(evoked_paths):
         evoked = _read_evoked_with_baseline(evoked_paths[subject])
-        row = {
-            "subject": subject,
-            "n_avg": int(evoked.nave),
+        row = {"subject": subject, "n_avg": int(evoked.nave)}
+        if n_trials_presented:
+            row["usable_pct"] = 100.0 * int(evoked.nave) / float(n_trials_presented)
+        row.update({
             "rms_snr": None,
             "band_power_90_110hz": compute_power(
                 evoked, f_low=90, f_high=110, t_low=resp_lower, t_high=resp_upper,
             ),
             "response_consistency": None,
-        }
+        })
         if evoked.baseline is not None:
             row["rms_snr"] = rms_snr(evoked, response_lower=resp_lower, response_upper=resp_upper)
 
@@ -224,10 +313,94 @@ def compute_subject_metrics(evoked_paths, epochs_paths=None, response_window=(0.
         if consistencies:
             row["response_consistency"] = sum(consistencies) / len(consistencies)
 
+        if use_sum:
+            summed = _load_polarity_sum(by_type_paths.get(subject))
+            if summed is None:
+                print(
+                    f"Warning: sub-{subject} does not have exactly two per-trial-type "
+                    "evoked files; polarity-sum metrics are left empty."
+                )
+            row.update(_polarity_sum_metrics(
+                summed, harmonics, stimulus, response_window, stim_cache,
+            ))
+
         rows.append(row)
     return pd.DataFrame(rows)
 
 
+def _polarity_sum_metrics(summed, harmonics, stimulus, response_window, stim_cache):
+    """Columns derived from the polarity-summed response (NaN when unavailable)."""
+    nan = float("nan")
+    columns = {"rms_snr_polarity_sum": nan}
+    if harmonics is not None:
+        columns.update({"f0_uv": nan, "upper_harmonics_uv": nan})
+    if stimulus is not None:
+        columns.update({"stim2resp_r": nan, "stim2resp_z": nan, "stim2resp_lag_ms": nan})
+        if stimulus.get("lag_range_ms") is not None:
+            columns.update({
+                "stim2resp_lim_r": nan, "stim2resp_lim_z": nan, "stim2resp_lim_lag_ms": nan,
+            })
+    if summed is None:
+        return columns
+
+    if summed.baseline is not None:
+        columns["rms_snr_polarity_sum"] = rms_snr(
+            summed, response_lower=response_window[0], response_upper=response_window[1],
+        )
+    if harmonics is not None:
+        result = harmonic_amplitudes(summed, **harmonics)
+        columns["f0_uv"] = result["f0"]
+        columns["upper_harmonics_uv"] = result["upper_harmonics"]
+    if stimulus is not None:
+        sfreq = float(summed.info["sfreq"])
+        if sfreq not in stim_cache:
+            wav, wav_sfreq = load_wav_mono(stimulus["path"])
+            stim_cache[sfreq] = resample_signal(wav, wav_sfreq, sfreq)
+        stim = stim_cache[sfreq]
+        columns.update(_stim_to_resp_columns(summed, stim, stimulus))
+        if stimulus.get("lag_range_ms") is not None:
+            lim_spec = {
+                "stim_window": stimulus["stim_window"],
+                "resp_window": stimulus.get("lag_resp_window", stimulus["resp_window"]),
+                "lag_range_ms": stimulus["lag_range_ms"],
+            }
+            columns.update(_stim_to_resp_columns(summed, stim, lim_spec, prefix="stim2resp_lim"))
+    return columns
+
+
+def merge_covariates(metrics, covariate_paths):
+    """Left-join subject-level covariate TSVs onto a metrics table.
+
+    Each TSV needs a ``participant_id`` column (``sub-<label>``, as in a
+    BIDS ``participants.tsv`` or ``phenotype/*.tsv``). Columns already
+    present in ``metrics`` are not overwritten. Missing files are skipped
+    with a warning.
+
+    Returns
+    -------
+    pandas.DataFrame
+        ``metrics`` with the covariate columns appended.
+    """
+    merged = metrics.copy()
+    key = "sub-" + merged["subject"].astype(str)
+    for path in covariate_paths or []:
+        path = Path(path)
+        if not path.exists():
+            print(f"Warning: covariate file not found, skipping: {path}")
+            continue
+        table = pd.read_csv(path, sep="\t", dtype={"participant_id": str})
+        if "participant_id" not in table.columns:
+            print(f"Warning: {path} has no 'participant_id' column, skipping.")
+            continue
+        new_columns = [c for c in table.columns if c != "participant_id" and c not in merged.columns]
+        if not new_columns:
+            continue
+        lookup = table.drop_duplicates("participant_id").set_index("participant_id")[new_columns]
+        for column in new_columns:
+            merged[column] = key.map(lookup[column])
+    return merged
+
+
 def _sort_key(key):
     """Sort key for (task, session, run, ...) tuples that may mix None
     with strings — plain tuple comparison raises TypeError in that case.
@@ -321,6 +494,59 @@ def _write_group_dataset_description(group_dir):
         json.dump(dataset_desc, f, indent=2)
 
 
+def _metric_options_from_args(args):
+    """Translate group-level CLI args into ``compute_subject_metrics`` options.
+
+    Every attribute is read with ``getattr`` so callers that build a bare
+    namespace (tests, scripts) keep working; unset options leave the
+    corresponding metrics off.
+    """
+    harmonics = None
+    f0 = getattr(args, "f0", None)
+    if f0:
+        tmin, tmax = getattr(args, "harmonic_window", None) or (0.060, 0.180)
+        harmonics = {
+            "f0": float(f0),
+            "n_harmonics": int(getattr(args, "n_harmonics", None) or 10),
+            "bin_hz": float(getattr(args, "harmonic_bin_hz", None) or 60.0),
+            "tmin": float(tmin),
+            "tmax": float(tmax),
+        }
+
+    stimulus = None
+    stimulus_path = getattr(args, "stimulus", None)
+    if stimulus_path:
+        stimulus = {
+            "path": Path(stimulus_path),
+            "stim_window": tuple(getattr(args, "xcorr_stim_window", None) or (0.050, 0.170)),
+            "resp_window": tuple(getattr(args, "xcorr_resp_window", None) or (0.060, 0.180)),
+        }
+        lag_range = getattr(args, "xcorr_lag_range", None)
+        if lag_range:
+            stimulus["lag_range_ms"] = tuple(lag_range)
+            lag_resp_window = getattr(args, "xcorr_lag_resp_window", None)
+            if lag_resp_window:
+                stimulus["lag_resp_window"] = tuple(lag_resp_window)
+
+    return {
+        "harmonics": harmonics,
+        "stimulus": stimulus,
+        "n_trials_presented": getattr(args, "n_trials_presented", None),
+    }
+
+
+def _covariate_paths_from_args(args):
+    """``participants.tsv`` (if present next to bids_dir) plus ``--covariates``."""
+    paths = []
+    bids_dir = getattr(args, "bids_dir", None)
+    if bids_dir:
+        participants = Path(bids_dir) / "participants.tsv"
+        if participants.exists():
+            paths.append(participants)
+    paths.extend(Path(p) for p in (getattr(args, "covariates", None) or []))
+    return paths
+
+
 def run_group_level(args):
     """
     Run the group-level aggregation step (BIDS-App ``analysis_level=group``).
@@ -333,12 +559,19 @@ def run_group_level(args):
     returns the list of output paths written.
 
     This step only aggregates outputs participant-level ffrprep has
-    already computed; it does not run any group-level statistics.
+    already computed; it does not run any group-level statistics. When
+    ``args.f0`` / ``args.stimulus`` are set, the metrics table also gets
+    harmonic-amplitude and stimulus-to-response columns computed from the
+    sum of the per-trial-type averages, and subject-level covariates from
+    ``<bids_dir>/participants.tsv`` and ``args.covariates`` are joined onto
+    the saved TSV (never used for inference here).
     """
     output_dir = Path(args.output_dir)
     response_window = tuple(getattr(args, "response_window", None) or (0.100, 0.200))
     participant_label = getattr(args, "participant_label", None)
     task_filter = getattr(args, "task", None)
+    metric_options = _metric_options_from_args(args)
+    covariate_paths = _covariate_paths_from_args(args)
 
     inputs = discover_group_inputs(output_dir, participant_label=participant_label, task=task_filter)
 
@@ -363,8 +596,13 @@ def run_group_level(args):
         written.append(evoked_path)
 
         epochs_paths = inputs["epochs"].get((task, session, run), {})
-        metrics = compute_subject_metrics(evoked_paths, epochs_paths, response_window=response_window)
-        metrics_path = save_group_metrics(output_dir, task, session, run, metrics)
+        metrics = compute_subject_metrics(
+            evoked_paths, epochs_paths, response_window=response_window,
+            by_type_paths=inputs["by_type"].get((task, session, run), {}),
+            **metric_options,
+        )
+        metrics_with_covariates = merge_covariates(metrics, covariate_paths)
+        metrics_path = save_group_metrics(output_dir, task, session, run, metrics_with_covariates)
         written.append(metrics_path)
 
         sections = [
diff --git a/ffrprep/reports.py b/ffrprep/reports.py
index 4298a2e..ffe44af 100644
--- a/ffrprep/reports.py
+++ b/ffrprep/reports.py
@@ -359,48 +359,75 @@ def build_evoked_section(
     }
 
 
-def build_metrics_table_section(metrics, *, section_id, title):
+def build_metrics_table_section(metrics, *, section_id, title, max_table_rows=40):
     """Build a section descriptor summarizing a group-level metrics table.
 
     ``metrics`` is a per-(task, run) DataFrame such as
     :func:`ffrprep.group.compute_subject_metrics` returns: one row per
     subject, columns are scalar FFR metrics recomputed from saved
     participant-level derivatives. The summary table reports the group
-    mean +/- SD for each metric; the full per-subject table is rendered
-    as a figure so any number of subjects fits without changing the
-    report layout.
+    mean +/- SD for each metric. Up to ``max_table_rows`` subjects the
+    full per-subject table is rendered as a figure; for larger cohorts
+    (where a table would be unreadably tall) the figure shows one
+    histogram per metric instead and the per-subject values live in the
+    accompanying ``_metrics.tsv``.
     """
+    numeric_columns = [
+        column for column in metrics.select_dtypes("number").columns if column != "n_avg"
+    ]
     summary = {"Subjects": str(len(metrics))}
-    for column in metrics.select_dtypes("number").columns:
-        if column == "n_avg":
-            continue
+    for column in numeric_columns:
         values = metrics[column].dropna()
         if values.empty:
             continue
-        summary[column.replace("_", " ")] = f"{values.mean():.3g} ± {values.std():.3g}"
-
-    # Column-wise ".3g" formatting (rather than a fixed .round(4)) so
-    # small-magnitude columns like band power (~1e-12) don't collapse
-    # to "0.0" next to O(1) columns like RMS SNR.
-    formatted = metrics.astype(object)
-    for column in metrics.select_dtypes("number").columns:
-        formatted[column] = metrics[column].map(lambda v: f"{v:.3g}")
+        summary[column.replace("_", " ")] = f"{values.mean():.3g} \u00b1 {values.std():.3g}"
+
+    if len(metrics) > max_table_rows and numeric_columns:
+        n_panels = len(numeric_columns)
+        n_cols = min(3, n_panels)
+        n_rows = -(-n_panels // n_cols)
+        fig, axes = plt.subplots(
+            n_rows, n_cols, figsize=(4.2 * n_cols, 3.0 * n_rows), squeeze=False,
+        )
+        for ax, column in zip(axes.ravel(), numeric_columns):
+            values = metrics[column].dropna()
+            ax.hist(values, bins=min(30, max(5, len(values) // 4)), color="#6B4E8C")
+            ax.set_title(column.replace("_", " "), fontsize=10)
+            ax.tick_params(labelsize=8)
+        for ax in axes.ravel()[n_panels:]:
+            ax.axis("off")
+        fig.tight_layout()
+        figure_title = "Metric distributions"
+        caption = (
+            f"Distribution of each metric across {len(metrics)} subjects; the "
+            "per-subject values are in the accompanying _metrics.tsv."
+        )
+    else:
+        # Column-wise ".3g" formatting (rather than a fixed .round(4)) so
+        # small-magnitude columns like band power (~1e-12) don't collapse
+        # to "0.0" next to O(1) columns like RMS SNR.
+        formatted = metrics.astype(object)
+        for column in metrics.select_dtypes("number").columns:
+            formatted[column] = metrics[column].map(lambda v: f"{v:.3g}")
+
+        fig, ax = plt.subplots(
+            figsize=(max(6, 1.4 * len(metrics.columns)), 0.4 * (len(metrics) + 1) + 0.5)
+        )
+        ax.axis("off")
+        table = ax.table(
+            cellText=formatted.values,
+            colLabels=list(metrics.columns),
+            loc="center",
+        )
+        table.auto_set_font_size(False)
+        table.set_fontsize(9)
+        fig.tight_layout()
+        figure_title = "Per-subject metrics"
+        caption = "Recomputed from saved participant-level derivatives; not persisted elsewhere."
 
-    fig, ax = plt.subplots(
-        figsize=(max(6, 1.4 * len(metrics.columns)), 0.4 * (len(metrics) + 1) + 0.5)
-    )
-    ax.axis("off")
-    table = ax.table(
-        cellText=formatted.values,
-        colLabels=list(metrics.columns),
-        loc="center",
-    )
-    table.auto_set_font_size(False)
-    table.set_fontsize(9)
-    fig.tight_layout()
     figures = [{
-        "title": "Per-subject metrics",
-        "caption": "Recomputed from saved participant-level derivatives; not persisted elsewhere.",
+        "title": figure_title,
+        "caption": caption,
         "data_uri": _fig_to_data_uri(fig),
     }]
     plt.close(fig)
diff --git a/ffrprep/tests/test_cli.py b/ffrprep/tests/test_cli.py
index 815a013..d474ad5 100644
--- a/ffrprep/tests/test_cli.py
+++ b/ffrprep/tests/test_cli.py
@@ -425,6 +425,41 @@ def test_workflow_inputs_carry_filter_method_and_reject_mode(tmp_path):
         assert inputs.reject_mode == "abs"
 
 
+def test_group_level_option_defaults_leave_extra_metrics_off():
+    args = get_parser().parse_args(["/b", "/o", "group"])
+    assert args.f0 is None
+    assert args.stimulus is None
+    assert args.n_harmonics == 10
+    assert args.harmonic_bin_hz == 60.0
+    assert args.harmonic_window == [0.06, 0.18]
+    assert args.xcorr_stim_window == [0.05, 0.17]
+    assert args.xcorr_resp_window == [0.06, 0.18]
+    assert args.xcorr_lag_range is None
+    assert args.n_trials_presented is None
+    assert args.covariates is None
+
+
+def test_group_level_options_parse_and_reach_group_module():
+    from ffrprep.group import _covariate_paths_from_args, _metric_options_from_args
+
+    args = get_parser().parse_args([
+        "/bids", "/out", "group",
+        "--f0", "100", "--harmonic-window", "0.06", "0.18", "--n-harmonics", "8",
+        "--stimulus", "/stim/da.wav", "--xcorr-lag-range", "6.9", "10.9",
+        "--xcorr-lag-resp-window", "0.05", "0.2", "--n-trials-presented", "6000",
+        "--covariates", "/pheno/a.tsv", "/pheno/b.tsv",
+    ])
+    options = _metric_options_from_args(args)
+    assert options["harmonics"] == {
+        "f0": 100.0, "n_harmonics": 8, "bin_hz": 60.0, "tmin": 0.06, "tmax": 0.18,
+    }
+    assert options["stimulus"]["path"] == Path("/stim/da.wav")
+    assert options["stimulus"]["lag_range_ms"] == (6.9, 10.9)
+    assert options["stimulus"]["lag_resp_window"] == (0.05, 0.2)
+    assert options["n_trials_presented"] == 6000
+    assert _covariate_paths_from_args(args) == [Path("/pheno/a.tsv"), Path("/pheno/b.tsv")]
+
+
 def test_no_report_flag_defaults_off_and_can_be_enabled():
     parser = get_parser()
     base = ["/path/to/bids", "/path/to/output", "participant"]
diff --git a/ffrprep/tests/test_group.py b/ffrprep/tests/test_group.py
index b688908..7f3a8b5 100644
--- a/ffrprep/tests/test_group.py
+++ b/ffrprep/tests/test_group.py
@@ -9,6 +9,7 @@
 import matplotlib
 import mne
 import numpy as np
+import pytest
 
 from ffrprep.group import (
     compute_grand_average,
@@ -252,3 +253,255 @@ def test_run_group_level_with_no_derivatives_prints_message_and_returns_empty(tm
 
     assert written == []
     assert "No participant-level ffrprep-analysis derivatives found" in capsys.readouterr().out
+
+
+# ---------------------------------------------------------------------------
+# per-type evokeds, polarity sum, harmonic / stimulus metrics, covariates
+# ---------------------------------------------------------------------------
+
+DA_SFREQ = 8192.0
+DA_TMIN, DA_TMAX = -0.04, 0.213
+
+
+def _chirp(t):
+    """A 170 ms 100->120 Hz glide with a smooth envelope (unique xcorr peak)."""
+    inside = (t >= 0) & (t <= 0.17)
+    envelope = np.sin(np.pi * np.clip(t, 0, 0.17) / 0.17) ** 2
+    phase = 2 * np.pi * (100.0 * t + 0.5 * (20.0 / 0.17) * t ** 2)
+    return np.where(inside, np.sin(phase) * envelope, 0.0)
+
+
+def _da_evoked(data_v, nave=3000):
+    info = mne.create_info(["Cz"], DA_SFREQ, ch_types=["eeg"])
+    evoked = mne.EvokedArray(
+        data_v[np.newaxis, :], info, tmin=DA_TMIN, nave=nave, verbose=False,
+    )
+    evoked.baseline = (DA_TMIN, 0.0)
+    return evoked
+
+
+def _da_times():
+    n = int(round((DA_TMAX - DA_TMIN) * DA_SFREQ)) + 1
+    return DA_TMIN + np.arange(n) / DA_SFREQ
+
+
+def _write_da_subject(root, subject, amplitude_uv=1.0, delay_s=0.012, seed=0):
+    """Combined + two per-polarity evokeds; polarity artifact cancels in the sum."""
+    rng = np.random.default_rng(seed)
+    times = _da_times()
+    signal = amplitude_uv * 1e-6 * _chirp(times - delay_s)
+    artifact = 5e-6 * rng.normal(size=times.size)
+    pol1, pol2 = signal + artifact, signal - artifact
+
+    subject_dir = root / "ffrprep-analysis" / f"sub-{subject}"
+    subject_dir.mkdir(parents=True, exist_ok=True)
+    base = f"sub-{subject}_task-da_run-01"
+    _da_evoked(pol1).save(subject_dir / f"{base}_desc-evokedDaPol1.fif", overwrite=True)
+    _da_evoked(pol2).save(subject_dir / f"{base}_desc-evokedDaPol2.fif", overwrite=True)
+    _da_evoked((pol1 + pol2) / 2, nave=6000).save(
+        subject_dir / f"{base}_desc-evoked.fif", overwrite=True,
+    )
+    return signal
+
+
+def _write_stimulus_wav(path, sfreq=24414):
+    from scipy.io import wavfile
+
+    t = np.arange(int(round(0.213 * sfreq))) / sfreq
+    wavfile.write(path, sfreq, (_chirp(t) * 20000).astype(np.int16))
+
+
+def test_discover_group_inputs_collects_per_type_evoked(tmp_path):
+    for subject in ("01", "02"):
+        _write_da_subject(tmp_path, subject, seed=int(subject))
+    inputs = discover_group_inputs(tmp_path)
+
+    key = ("da", None, "01")
+    assert set(inputs["by_type"][key]["01"]) == {"DaPol1", "DaPol2"}
+    assert set(inputs["by_type"][key].keys()) == {"01", "02"}
+    # the combined file is still the "evoked" entry, not a per-type one
+    assert set(inputs["evoked"][key].keys()) == {"01", "02"}
+    assert inputs["diff"] == {}
+
+
+def test_load_polarity_sum_is_the_plain_sum_and_needs_two_types(tmp_path):
+    from ffrprep.group import _load_polarity_sum
+
+    _write_da_subject(tmp_path, "01", seed=3)
+    paths = discover_group_inputs(tmp_path)["by_type"][("da", None, "01")]["01"]
+    summed = _load_polarity_sum(paths)
+
+    first = mne.read_evokeds(str(paths["DaPol1"]), condition=0, verbose=False)
+    second = mne.read_evokeds(str(paths["DaPol2"]), condition=0, verbose=False)
+    np.testing.assert_allclose(summed.data, first.data + second.data, atol=1e-15)
+    assert summed.baseline is not None
+    assert _load_polarity_sum({"DaPol1": paths["DaPol1"]}) is None
+    assert _load_polarity_sum(None) is None
+
+
+def test_compute_subject_metrics_polarity_sum_harmonics_and_stimulus(tmp_path):
+    from ffrprep.analysis import harmonic_amplitudes
+
+    signal = _write_da_subject(tmp_path, "01", amplitude_uv=1.0, delay_s=0.012, seed=1)
+    _write_da_subject(tmp_path, "02", amplitude_uv=2.0, delay_s=0.012, seed=2)
+    stim_path = tmp_path / "stim.wav"
+    _write_stimulus_wav(stim_path)
+
+    inputs = discover_group_inputs(tmp_path)
+    key = ("da", None, "01")
+    metrics = compute_subject_metrics(
+        inputs["evoked"][key],
+        by_type_paths=inputs["by_type"][key],
+        response_window=(0.0, 0.213),
+        harmonics={"f0": 100.0, "n_harmonics": 10, "bin_hz": 60.0, "tmin": 0.06, "tmax": 0.18},
+        stimulus={
+            "path": stim_path, "stim_window": (0.05, 0.17), "resp_window": (0.06, 0.18),
+        },
+        n_trials_presented=6000,
+    )
+
+    assert list(metrics["subject"]) == ["01", "02"]
+    assert list(metrics["n_avg"]) == [6000, 6000]
+    assert metrics["usable_pct"].tolist() == [100.0, 100.0]
+    for column in ("rms_snr_polarity_sum", "f0_uv", "upper_harmonics_uv",
+                   "stim2resp_r", "stim2resp_z", "stim2resp_lag_ms"):
+        assert metrics[column].notna().all(), column
+
+    # harmonics come from 2 * signal (the sum of the two polarity averages)
+    expected = harmonic_amplitudes(_da_evoked(2 * signal))
+    assert metrics.loc[0, "f0_uv"] == pytest.approx(expected["f0"], rel=1e-6)
+    assert metrics.loc[0, "upper_harmonics_uv"] == pytest.approx(expected["upper_harmonics"], rel=1e-6)
+    # a 2x larger response gives 2x larger spectral amplitudes
+    assert metrics.loc[1, "f0_uv"] == pytest.approx(2 * metrics.loc[0, "f0_uv"], rel=1e-6)
+    # stimulus following: strong correlation, z = atanh(r), lag = delay - 10 ms window offset
+    assert metrics.loc[0, "stim2resp_r"] > 0.9
+    assert metrics.loc[0, "stim2resp_z"] == pytest.approx(np.arctanh(metrics.loc[0, "stim2resp_r"]))
+    assert metrics.loc[0, "stim2resp_lag_ms"] == pytest.approx(2.0, abs=0.3)
+
+
+def test_compute_subject_metrics_lag_limited_variant(tmp_path):
+    _write_da_subject(tmp_path, "01", delay_s=0.012, seed=1)
+    stim_path = tmp_path / "stim.wav"
+    _write_stimulus_wav(stim_path)
+    inputs = discover_group_inputs(tmp_path)
+    key = ("da", None, "01")
+
+    metrics = compute_subject_metrics(
+        {"01": inputs["evoked"][key]["01"]},
+        by_type_paths=inputs["by_type"][key],
+        stimulus={
+            "path": stim_path, "stim_window": (0.05, 0.17), "resp_window": (0.06, 0.18),
+            "lag_range_ms": (6.9, 10.9), "lag_resp_window": (0.05, 0.20),
+        },
+    )
+    # both segments start at 50 ms, so the lag is the true 12 ms delay's neighbourhood
+    assert 6.9 <= metrics.loc[0, "stim2resp_lim_lag_ms"] <= 10.9
+    assert metrics.loc[0, "stim2resp_lim_r"] <= metrics.loc[0, "stim2resp_r"] + 1e-9
+
+
+def test_polarity_sum_metrics_are_nan_without_two_types(tmp_path, capsys):
+    _write_da_subject(tmp_path, "01", seed=1)
+    inputs = discover_group_inputs(tmp_path)
+    key = ("da", None, "01")
+    only_one = {"01": {"DaPol1": inputs["by_type"][key]["01"]["DaPol1"]}}
+
+    metrics = compute_subject_metrics(
+        {"01": inputs["evoked"][key]["01"]},
+        by_type_paths=only_one,
+        harmonics={"f0": 100.0},
+    )
+    assert np.isnan(metrics.loc[0, "f0_uv"])
+    assert np.isnan(metrics.loc[0, "rms_snr_polarity_sum"])
+    assert "exactly two per-trial-type" in capsys.readouterr().out
+
+
+def test_metrics_without_new_options_keep_original_columns(tmp_path):
+    _write_da_subject(tmp_path, "01", seed=1)
+    inputs = discover_group_inputs(tmp_path)
+    key = ("da", None, "01")
+    metrics = compute_subject_metrics({"01": inputs["evoked"][key]["01"]})
+    assert list(metrics.columns) == [
+        "subject", "n_avg", "rms_snr", "band_power_90_110hz", "response_consistency",
+    ]
+
+
+def test_merge_covariates_joins_on_participant_id_without_overwriting(tmp_path, capsys):
+    import pandas as pd
+
+    from ffrprep.group import merge_covariates
+
+    metrics = pd.DataFrame({"subject": ["bu001", "bu002", "zz999"], "f0_uv": [0.1, 0.2, 0.3]})
+    participants = tmp_path / "participants.tsv"
+    participants.write_text(
+        "participant_id\tsite\tgroup\tage\nsub-bu001\tBU\tMus\t20\nsub-bu002\tBU\tNMus\t33\n"
+    )
+    phenotype = tmp_path / "pheno.tsv"
+    phenotype.write_text("participant_id\tage\tdp\nsub-bu001\t99\t1.5\nsub-bu002\t99\t0.5\n")
+
+    merged = merge_covariates(metrics, [participants, phenotype, tmp_path / "missing.tsv"])
+
+    assert list(merged["group"][:2]) == ["Mus", "NMus"]
+    assert pd.isna(merged.loc[2, "group"])  # subject without a covariate row
+    assert list(merged["age"][:2]) == [20, 33]  # first file wins, not overwritten by pheno
+    assert list(merged["dp"][:2]) == [1.5, 0.5]
+    assert "covariate file not found" in capsys.readouterr().out
+
+
+def test_run_group_level_end_to_end_with_harmonics_stimulus_and_covariates(tmp_path):
+    import pandas as pd
+
+    bids_dir = tmp_path / "bids"
+    bids_dir.mkdir()
+    (bids_dir / "participants.tsv").write_text(
+        "participant_id\tsite\tgroup\n"
+        "sub-01\tBU\tMus\nsub-02\tBU\tNMus\nsub-03\tCMU\tVar\n"
+    )
+    out = tmp_path / "out"
+    for i, subject in enumerate(["01", "02", "03"], start=1):
+        _write_da_subject(out, subject, amplitude_uv=float(i), seed=i)
+    stim_path = tmp_path / "stim.wav"
+    _write_stimulus_wav(stim_path)
+
+    args = SimpleNamespace(
+        output_dir=out, bids_dir=bids_dir, participant_label=None, task=None,
+        response_window=(0.0, 0.213), f0=100.0, n_harmonics=10, harmonic_bin_hz=60.0,
+        harmonic_window=(0.06, 0.18), stimulus=stim_path,
+        xcorr_stim_window=(0.05, 0.17), xcorr_resp_window=(0.06, 0.18),
+        xcorr_lag_range=None, xcorr_lag_resp_window=None,
+        n_trials_presented=6000, covariates=None,
+    )
+    written = run_group_level(args)
+
+    group_dir = out / "ffrprep-group"
+    table = pd.read_csv(group_dir / "task-da_run-01_metrics.tsv", sep="\t", dtype={"subject": str})
+    assert list(table["subject"]) == ["01", "02", "03"]
+    assert list(table["group"]) == ["Mus", "NMus", "Var"]
+    assert list(table["site"]) == ["BU", "BU", "CMU"]
+    for column in ("f0_uv", "upper_harmonics_uv", "stim2resp_z", "usable_pct"):
+        assert table[column].notna().all(), column
+    assert table["f0_uv"].is_monotonic_increasing  # amplitudes 1, 2, 3 uV
+    assert (group_dir / "group_report.html").exists()
+    assert all(path.exists() for path in written)
+
+
+def test_metrics_section_uses_histograms_for_large_cohorts():
+    import pandas as pd
+
+    from ffrprep.reports import build_metrics_table_section
+
+    rng = np.random.default_rng(0)
+    big = pd.DataFrame({
+        "subject": [f"{i:03d}" for i in range(80)],
+        "n_avg": 6000,
+        "f0_uv": rng.normal(0.14, 0.03, 80),
+        "upper_harmonics_uv": rng.normal(0.2, 0.05, 80),
+    })
+    small = big.head(5)
+
+    assert build_metrics_table_section(big, section_id="a", title="t")["figures"][0]["title"] == (
+        "Metric distributions"
+    )
+    assert build_metrics_table_section(small, section_id="b", title="t")["figures"][0]["title"] == (
+        "Per-subject metrics"
+    )
+    assert "f0 uv" in build_metrics_table_section(big, section_id="a", title="t")["summary"]

From 78b74406fa3eb8743e2859640f2b0a208366b3e0 Mon Sep 17 00:00:00 2001
From: Kevin Sitek <sitek.kevin@gmail.com>
Date: Wed, 23 Sep 2026 22:45:22 -0500
Subject: [PATCH 04/10] Do not lag-restrict the primary stimulus-to-response
 correlation

The stim2resp_r/z/lag_ms columns must search all lags; --xcorr-lag-range
only defines the separate stim2resp_lim_* variant. The primary columns were
accidentally computed with the restriction whenever a lag range was given.
Validated on real data: max-over-all-lags r = 0.1607 vs 0.1610 published.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 ffrprep/group.py            | 8 +++++++-
 ffrprep/tests/test_group.py | 8 ++++++--
 2 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/ffrprep/group.py b/ffrprep/group.py
index c8a7007..3d1e6df 100644
--- a/ffrprep/group.py
+++ b/ffrprep/group.py
@@ -357,7 +357,13 @@ def _polarity_sum_metrics(summed, harmonics, stimulus, response_window, stim_cac
             wav, wav_sfreq = load_wav_mono(stimulus["path"])
             stim_cache[sfreq] = resample_signal(wav, wav_sfreq, sfreq)
         stim = stim_cache[sfreq]
-        columns.update(_stim_to_resp_columns(summed, stim, stimulus))
+        # The primary measure searches ALL lags; ``lag_range_ms`` belongs to
+        # the separate lag-limited ("lim") variant below.
+        primary_spec = {
+            "stim_window": stimulus["stim_window"],
+            "resp_window": stimulus["resp_window"],
+        }
+        columns.update(_stim_to_resp_columns(summed, stim, primary_spec))
         if stimulus.get("lag_range_ms") is not None:
             lim_spec = {
                 "stim_window": stimulus["stim_window"],
diff --git a/ffrprep/tests/test_group.py b/ffrprep/tests/test_group.py
index 7f3a8b5..a5e8ac0 100644
--- a/ffrprep/tests/test_group.py
+++ b/ffrprep/tests/test_group.py
@@ -394,9 +394,13 @@ def test_compute_subject_metrics_lag_limited_variant(tmp_path):
             "lag_range_ms": (6.9, 10.9), "lag_resp_window": (0.05, 0.20),
         },
     )
-    # both segments start at 50 ms, so the lag is the true 12 ms delay's neighbourhood
+    # The lag-limited variant is confined to its window (the true 12 ms delay lies outside)...
     assert 6.9 <= metrics.loc[0, "stim2resp_lim_lag_ms"] <= 10.9
-    assert metrics.loc[0, "stim2resp_lim_r"] <= metrics.loc[0, "stim2resp_r"] + 1e-9
+    # ...while the primary measure searches ALL lags: its windows are offset by 10 ms, so
+    # the 12 ms delay appears at a 2 ms array lag, outside the restriction, with a
+    # strictly higher correlation than anything inside the 6.9-10.9 ms window.
+    assert metrics.loc[0, "stim2resp_lag_ms"] == pytest.approx(2.0, abs=0.3)
+    assert metrics.loc[0, "stim2resp_r"] > metrics.loc[0, "stim2resp_lim_r"] + 0.05
 
 
 def test_polarity_sum_metrics_are_nan_without_two_types(tmp_path, capsys):

From 6c10251b589b20ac2dd21765f80c4d646ffa06d6 Mon Sep 17 00:00:00 2001
From: Kevin Sitek <sitek.kevin@gmail.com>
Date: Wed, 23 Sep 2026 23:17:15 -0500
Subject: [PATCH 05/10] Grand-average single-channel evokeds with differing
 channel names

Multi-site cohorts can label the one recorded channel differently (A32 vs
Cz) and mne.grand_average refuses to combine them. Rename single-channel
evokeds to the most common name (with a note); mismatched multi-channel
sets raise an explicit error.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 CHANGES.rst                 |  3 +++
 ffrprep/group.py            | 36 ++++++++++++++++++++++++++++++++++++
 ffrprep/tests/test_group.py | 32 ++++++++++++++++++++++++++++++++
 3 files changed, 71 insertions(+)

diff --git a/CHANGES.rst b/CHANGES.rst
index 3413ef5..1354442 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -71,6 +71,9 @@ Group-level analysis
   (and ``stim2resp_lim_*`` with ``--xcorr-lag-range``); ``--n-trials-presented``
   adds ``usable_pct``. ``discover_group_inputs`` now also returns the
   per-trial-type evoked files (``"by_type"``).
+- ``compute_grand_average`` renames single-channel evokeds that carry different
+  channel names across sites (e.g. ``A32`` vs ``Cz``) to a common name before
+  averaging; mismatched multi-channel sets raise a clear ``ValueError``.
 - The group metrics TSV joins subject-level covariates from
   ``<bids_dir>/participants.tsv`` and ``--covariates`` (``merge_covariates``);
   the report section switches to per-metric histograms for cohorts above
diff --git a/ffrprep/group.py b/ffrprep/group.py
index 3d1e6df..8266d76 100644
--- a/ffrprep/group.py
+++ b/ffrprep/group.py
@@ -181,10 +181,46 @@ def compute_grand_average(evoked_paths):
     """
     subjects = sorted(evoked_paths)
     evokeds = [_read_evoked_with_baseline(evoked_paths[subject]) for subject in subjects]
+    _harmonize_channel_names(evokeds, subjects)
     grand_average = mne.grand_average(evokeds)
     return grand_average, subjects
 
 
+def _harmonize_channel_names(evokeds, subjects):
+    """Give single-channel Evokeds a common channel name (in place).
+
+    Multi-site FFR cohorts often record one channel that different sites
+    label differently (e.g. ``A32`` vs ``Cz``); ``mne.grand_average``
+    refuses to combine those. When every Evoked has exactly one channel
+    the name carries no information, so all are renamed to the most
+    common one. Differing names across multi-channel Evokeds cannot be
+    aligned automatically and raise a ``ValueError``.
+    """
+    names = [tuple(ev.ch_names) for ev in evokeds]
+    if len(set(names)) <= 1:
+        return
+    if all(len(n) == 1 for n in names):
+        counts = {}
+        for (name,) in names:
+            counts[name] = counts.get(name, 0) + 1
+        target = max(sorted(counts), key=counts.get)
+        print(
+            f"Note: single-channel evokeds use different channel names {sorted(counts)}; "
+            f"treating them as the same channel and renaming to '{target}'."
+        )
+        for ev in evokeds:
+            if ev.ch_names[0] != target:
+                ev.rename_channels({ev.ch_names[0]: target})
+        return
+    first = names[0]
+    for subject, name in zip(subjects, names):
+        if name != first:
+            raise ValueError(
+                f"sub-{subject} has channels {list(name)} but sub-{subjects[0]} has "
+                f"{list(first)}; grand averaging needs identical channel sets."
+            )
+
+
 def _load_polarity_sum(cond_paths):
     """Sum of a subject's two per-trial-type Evokeds (or None).
 
diff --git a/ffrprep/tests/test_group.py b/ffrprep/tests/test_group.py
index a5e8ac0..6bddaeb 100644
--- a/ffrprep/tests/test_group.py
+++ b/ffrprep/tests/test_group.py
@@ -509,3 +509,35 @@ def test_metrics_section_uses_histograms_for_large_cohorts():
         "Per-subject metrics"
     )
     assert "f0 uv" in build_metrics_table_section(big, section_id="a", title="t")["summary"]
+
+
+def test_compute_grand_average_harmonizes_single_channel_names(tmp_path, capsys):
+    """Sites often name the one recorded channel differently (A32 vs Cz)."""
+    paths = {}
+    for subject, ch_name in (("01", "A32"), ("02", "A32"), ("03", "Cz")):
+        info = mne.create_info([ch_name], 1000.0, ch_types=["eeg"])
+        data = np.random.default_rng(int(subject)).normal(0, 1e-6, (1, 300))
+        evoked = mne.EvokedArray(data, info, tmin=-0.04, nave=10, verbose=False)
+        evoked.baseline = (-0.04, 0.0)
+        path = tmp_path / f"sub-{subject}_ave.fif"
+        evoked.save(path, overwrite=True)
+        paths[subject] = path
+
+    grand_average, subjects = compute_grand_average(paths)
+
+    assert subjects == ["01", "02", "03"]
+    assert grand_average.ch_names == ["A32"]  # most common name wins
+    assert "different channel names" in capsys.readouterr().out
+
+
+def test_compute_grand_average_rejects_mismatched_multichannel_names(tmp_path):
+    paths = {}
+    for subject, names in (("01", ["Cz", "Pz"]), ("02", ["Cz", "Fz"])):
+        info = mne.create_info(names, 1000.0, ch_types=["eeg"] * 2)
+        evoked = mne.EvokedArray(np.zeros((2, 100)), info, tmin=0.0, nave=5, verbose=False)
+        path = tmp_path / f"sub-{subject}_ave.fif"
+        evoked.save(path, overwrite=True)
+        paths[subject] = path
+
+    with pytest.raises(ValueError, match="identical channel sets"):
+        compute_grand_average(paths)

From 710352d8f3a2597f3a85edd415c3c28ebc21ab87 Mon Sep 17 00:00:00 2001
From: Kevin Sitek <sitek.kevin@gmail.com>
Date: Thu, 24 Sep 2026 16:09:31 -0500
Subject: [PATCH 06/10] Add --skip-existing, --clean-work-dir and
 --no-keep-epochs

Bulk runs over hundreds of subjects were limited by disk space (2.5-4 GB
of Nipype working files per subject) and needed external wrapper scripts
to resume after interruptions.

- --skip-existing skips (task, run) iterations whose JSON sidecars already
  exist; completion is judged from sidecars so it survives --no-keep-epochs.
  A finished subject builds no report and no longer errors when the epochs
  were removed.
- --clean-work-dir removes an iteration's Nipype working files after it
  succeeds (logs are kept; failed iterations keep theirs for debugging).
- --keep-epochs/--no-keep-epochs (default keep) deletes *_epo.fif after the
  analysis stage; ignored with a warning for --stage preprocessing.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 CHANGES.rst                             |   6 +
 docs/source/usage.rst                   |  27 ++++
 ffrprep/ffrprep_cli.py                  | 181 ++++++++++++++++++++---
 ffrprep/tests/test_cli.py               | 182 ++++++++++++++++++++++++
 ffrprep/tests/test_parallel_dispatch.py | 129 +++++++++++++++++
 5 files changed, 506 insertions(+), 19 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 5faa023..6a86b83 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -36,6 +36,12 @@ CLI
 - New flag ``--no-report``: skip HTML report generation (preprocessing
   and analysis) while still writing all derivatives. Intended for bulk
   runs over many subjects, where report figures dominate runtime.
+- New flags for resumable, disk-friendly bulk runs: ``--skip-existing``
+  (skip (task, run) iterations whose sidecar outputs already exist),
+  ``--clean-work-dir`` (delete each iteration's Nipype working files after
+  it succeeds; logs are kept) and ``--keep-epochs`` /
+  ``--no-keep-epochs`` (default keep; ``--no-keep-epochs`` deletes the
+  ``*_epo.fif`` files after analysis and keeps the JSON sidecars).
 - New flag ``--with-stimuli`` on ``ffrprep-download example``
   (``BooleanOptionalAction``, default ``False``): additionally
   fetches the BIDS ``/stimuli/`` directory needed by stimulus-aware
diff --git a/docs/source/usage.rst b/docs/source/usage.rst
index 64e433a..41e8676 100644
--- a/docs/source/usage.rst
+++ b/docs/source/usage.rst
@@ -353,6 +353,33 @@ Cluster (GNU parallel on a single beefy box)
          --participant_label {} --n_procs 4" \
       ::: 01 02 03 04 05 06 07 08
 
+Resuming and saving disk space
+------------------------------
+
+Large cohorts are limited by disk space more than by CPU: each subject's
+Nipype working directory can hold several GB, and the epochs files are
+much larger than the evokeds. Three flags make bulk runs resumable and
+lean:
+
+.. code-block:: bash
+
+    ffrprep /data /data/derivatives participant \
+        --participant_label 01 02 03 \
+        --skip-existing --clean-work-dir --no-keep-epochs
+
+* ``--skip-existing`` skips (task, run) iterations whose JSON sidecar
+  outputs already exist, so an interrupted run can simply be started again.
+  Completion is judged from the sidecars, so it still works after
+  ``--no-keep-epochs``. Iterations that are already complete build no
+  report.
+* ``--clean-work-dir`` deletes each iteration's Nipype working files once
+  it has succeeded. Worker logs are kept, and the working files of a failed
+  iteration are kept for debugging.
+* ``--no-keep-epochs`` deletes the ``*_epo.fif`` files after the analysis
+  stage (evokeds and all JSON sidecars stay). It is ignored, with a
+  warning, for ``--stage preprocessing`` because the analysis stage still
+  needs the epochs.
+
 Failure handling
 ----------------
 
diff --git a/ffrprep/ffrprep_cli.py b/ffrprep/ffrprep_cli.py
index 440da7b..94b64cf 100644
--- a/ffrprep/ffrprep_cli.py
+++ b/ffrprep/ffrprep_cli.py
@@ -2,6 +2,7 @@
 import logging
 import multiprocessing
 import os
+import shutil
 import sys
 import time
 from concurrent.futures import ProcessPoolExecutor, as_completed
@@ -223,6 +224,7 @@ def _make_preproc_payload(args_snap, deriv_snap, subject, task_label, run_label,
         "low_pass": effective_h,
         "filter_method": args_snap.get("filter_method", "fir"),
         "reject_mode": args_snap.get("reject_mode", "ptp"),
+        "clean_work_dir": bool(args_snap.get("clean_work_dir", False)),
         "baseline": baseline,
         "tmin": args_snap.get("tmin"),
         "tmax": args_snap.get("tmax"),
@@ -312,6 +314,7 @@ def _make_concat_payload(args_snap, deriv_snap, subject, task_label, runs,
         "low_pass": effective_h,
         "filter_method": args_snap.get("filter_method", "fir"),
         "reject_mode": args_snap.get("reject_mode", "ptp"),
+        "clean_work_dir": bool(args_snap.get("clean_work_dir", False)),
         "baseline": baseline,
         "tmin": args_snap.get("tmin"),
         "tmax": args_snap.get("tmax"),
@@ -505,6 +508,7 @@ def _make_analysis_payload(args_snap, deriv_snap, subject, group):
         ),
         "analysis_subject_dir": str(deriv_snap["analysis_subject_dir"]),
         "derivatives_root": str(deriv_snap.get("derivatives_root", "")),
+        "clean_work_dir": bool(args_snap.get("clean_work_dir", False)),
     }
 
 
@@ -1303,6 +1307,75 @@ def _build_analysis_report(args, derivatives_info, subject):
     print(f"Analysis report written to: {out_path}")
 
 
+def _find_preproc_outputs(payload, suffix=".fif"):
+    """Preprocessing outputs (``_desc-preproc*_epo<suffix>``) for one iteration.
+
+    ``suffix=".json"`` finds the sidecars instead of the epochs; those are
+    written with every output and survive ``--no-keep-epochs``, so they are
+    what ``--skip-existing`` keys on.
+    """
+    expected_dir = Path(payload["output_dir"])
+    subject = payload["subject"]
+    task_label = payload["task_label"]
+    candidates = sorted(expected_dir.glob(
+        f"sub-{subject}_*task-{task_label}_*desc-preproc*_epo{suffix}"
+    ))
+    if payload["kind"] == "per_run":
+        run_label = payload["run_label"]
+        if run_label is not None:
+            candidates = [c for c in candidates if f"run-{run_label}" in c.name]
+    else:  # concat
+        candidates = [c for c in candidates if "_run-" not in c.name]
+    return candidates
+
+
+def _preproc_done(payload):
+    """True when this preprocessing iteration's outputs already exist."""
+    return bool(_find_preproc_outputs(payload, suffix=".json"))
+
+
+def _analysis_done(payload):
+    """True when this (task, run) group's combined evoked sidecar exists."""
+    base_stem = Path(payload["preproc_files"][0]).stem.partition("_desc-preproc")[0]
+    return (Path(payload["analysis_subject_dir"]) / f"{base_stem}_desc-evoked.json").exists()
+
+
+def _analysis_outputs_complete(analysis_subject_dir):
+    """True when any combined-evoked sidecar exists in the subject's analysis dir."""
+    return any(Path(analysis_subject_dir).glob("*_desc-evoked.json"))
+
+
+def _select_pending(payloads, done_fn, skip_existing, label):
+    """Drop already-complete iterations when ``--skip-existing`` is set."""
+    if not skip_existing:
+        return payloads
+    pending = [p for p in payloads if not done_fn(p)]
+    print(
+        f"--skip-existing: {len(payloads) - len(pending)}/{len(payloads)} "
+        f"{label} iteration(s) already complete"
+    )
+    return pending
+
+
+def _remove_workflow_files(work_dir, workflow_name):
+    """Delete a finished iteration's nipype node directories, keeping ``*.log``.
+
+    nipype writes ``<work_dir>/<workflow name>/<node>/...``; the worker log
+    lives directly in ``work_dir`` and is not touched.
+    """
+    shutil.rmtree(Path(work_dir) / workflow_name, ignore_errors=True)
+
+
+def _remove_epoch_files(preproc_subject_dir):
+    """Delete ``*_epo.fif`` files (sidecars stay); returns (count, bytes freed)."""
+    count, freed = 0, 0
+    for path in Path(preproc_subject_dir).glob("*_epo.fif"):
+        freed += path.stat().st_size
+        path.unlink()
+        count += 1
+    return count, freed
+
+
 def _check_preproc_output_or_raise(payload):
     """Confirm at least one preprocessing output landed on disk.
 
@@ -1321,16 +1394,11 @@ def _check_preproc_output_or_raise(payload):
     expected_dir = Path(payload["output_dir"])
     subject = payload["subject"]
     task_label = payload["task_label"]
-    candidates = sorted(expected_dir.glob(
-        f"sub-{subject}_*task-{task_label}_*desc-preproc*_epo.fif"
-    ))
+    candidates = _find_preproc_outputs(payload)
     if payload["kind"] == "per_run":
         run_label = payload["run_label"]
-        if run_label is not None:
-            candidates = [c for c in candidates if f"run-{run_label}" in c.name]
         marker = f"run-{run_label}" if run_label is not None else "single"
     else:  # concat
-        candidates = [c for c in candidates if "_run-" not in c.name]
         marker = "concat"
     if not candidates:
         raise FileNotFoundError(
@@ -1357,6 +1425,8 @@ def _preproc_iteration(payload):
     wf = _build_preproc_workflow(payload)
     wf.run(plugin="Linear")
     outs = _check_preproc_output_or_raise(payload)
+    if payload.get("clean_work_dir"):
+        _remove_workflow_files(work_dir, wf.name)
     return IterationResult(
         identifier=payload["identifier"],
         output_files=[str(p) for p in outs],
@@ -1381,6 +1451,8 @@ def _concat_iteration(payload):
     wf = _build_preproc_workflow(payload)
     wf.run(plugin="Linear")
     outs = _check_preproc_output_or_raise(payload)
+    if payload.get("clean_work_dir"):
+        _remove_workflow_files(work_dir, wf.name)
     return IterationResult(
         identifier=payload["identifier"],
         output_files=[str(p) for p in outs],
@@ -1444,6 +1516,8 @@ def _analysis_iteration(payload):
     # derivatives root and constructs the per-subject path internally.
     analysis_wf.inputs.inputnode.derivatives_root = payload["derivatives_root"]
     analysis_wf.run(plugin="Linear")
+    if payload.get("clean_work_dir"):
+        _remove_workflow_files(work_dir, analysis_wf.name)
 
     # Propagate provenance from each per-condition source file's sidecar
     # into the corresponding analysis-output sidecars. The first file
@@ -1943,6 +2017,40 @@ def get_parser():
             "many subjects, where report figures dominate runtime."
         ),
     )
+    parser.add_argument(
+        "--skip-existing",
+        dest="skip_existing",
+        action="store_true",
+        default=False,
+        help=(
+            "Skip (task, run) iterations whose outputs already exist (matched "
+            "on the JSON sidecars, so this still works after --no-keep-epochs). "
+            "Makes an interrupted run resumable by re-issuing the same command."
+        ),
+    )
+    parser.add_argument(
+        "--clean-work-dir",
+        dest="clean_work_dir",
+        action="store_true",
+        default=False,
+        help=(
+            "Delete each iteration's nipype working files once it succeeds "
+            "(they can be several GB per subject); worker logs are kept. "
+            "Work dirs of failed iterations are left in place for debugging."
+        ),
+    )
+    parser.add_argument(
+        "--keep-epochs",
+        dest="keep_epochs",
+        action=argparse.BooleanOptionalAction,
+        default=True,
+        help=(
+            "Keep the large preprocessed *_epo.fif files after the analysis "
+            "stage. --no-keep-epochs deletes them (JSON sidecars stay) once "
+            "analysis and its report are done; ignored with --stage "
+            "preprocessing because analysis still needs them."
+        ),
+    )
     parser.add_argument(
         "--n_procs", type=int, default=1,
         help=(
@@ -2323,10 +2431,14 @@ def run_ffrprep():
                     )
                     for task_label in tasks
                 ]
-                _dispatch(_concat_iteration, payloads,
-                          n_procs=args.n_procs, kind_label="concat")
-                if not args.no_report:
-                    _build_preproc_report(args, derivatives_info, subject)
+                payloads = _select_pending(
+                    payloads, _preproc_done, args.skip_existing, "preprocessing",
+                )
+                if payloads:
+                    _dispatch(_concat_iteration, payloads,
+                              n_procs=args.n_procs, kind_label="concat")
+                    if not args.no_report:
+                        _build_preproc_report(args, derivatives_info, subject)
             else:
                 # One worker per (task, run) — fully independent iterations.
                 payloads = [
@@ -2338,10 +2450,14 @@ def run_ffrprep():
                     for task_label in tasks
                     for run in runs
                 ]
-                _dispatch(_preproc_iteration, payloads,
-                          n_procs=args.n_procs, kind_label="preproc")
-                if not args.no_report:
-                    _build_preproc_report(args, derivatives_info, subject)
+                payloads = _select_pending(
+                    payloads, _preproc_done, args.skip_existing, "preprocessing",
+                )
+                if payloads:
+                    _dispatch(_preproc_iteration, payloads,
+                              n_procs=args.n_procs, kind_label="preproc")
+                    if not args.no_report:
+                        _build_preproc_report(args, derivatives_info, subject)
 
         if args.stage in ["analysis", "both"]:
             print("\n" + "=" * 60)
@@ -2355,6 +2471,14 @@ def run_ffrprep():
             preproc_subject_dir = Path(derivatives_info["preprocessing_subject_dir"])
             groups = _collect_analysis_groups(preproc_subject_dir)
 
+            if not groups and args.skip_existing and _analysis_outputs_complete(
+                derivatives_info["analysis_subject_dir"]
+            ):
+                print(
+                    "--skip-existing: analysis outputs already complete "
+                    "(preprocessed epochs were removed); nothing to do"
+                )
+                continue
             if not groups:
                 print(f"ERROR: No preprocessing outputs found for subject {subject}.")
                 print(f"Expected location: {preproc_subject_dir}")
@@ -2371,11 +2495,30 @@ def run_ffrprep():
                 _make_analysis_payload(args_snap, deriv_snap, subject, g)
                 for g in groups
             ]
-            _dispatch(_analysis_iteration, payloads,
-                      n_procs=args.n_procs, kind_label="analysis")
-            print(f"Analysis completed. Outputs saved to: {derivatives_info['analysis_subject_dir']}")
-            if not args.no_report:
-                _build_analysis_report(args, derivatives_info, subject)
+            payloads = _select_pending(
+                payloads, _analysis_done, args.skip_existing, "analysis",
+            )
+            if payloads:
+                _dispatch(_analysis_iteration, payloads,
+                          n_procs=args.n_procs, kind_label="analysis")
+                print(f"Analysis completed. Outputs saved to: {derivatives_info['analysis_subject_dir']}")
+                if not args.no_report:
+                    _build_analysis_report(args, derivatives_info, subject)
+
+        if not args.keep_epochs:
+            if args.stage in ("analysis", "both"):
+                n_removed, freed = _remove_epoch_files(
+                    derivatives_info["preprocessing_subject_dir"],
+                )
+                print(
+                    f"--no-keep-epochs: removed {n_removed} epochs file(s) "
+                    f"({freed / 1e6:.0f} MB); JSON sidecars kept"
+                )
+            else:
+                print(
+                    "Warning: --no-keep-epochs is ignored with --stage "
+                    "preprocessing (the analysis stage still needs the epochs)."
+                )
 
     print("\n" + "=" * 60)
     print("ffrprep processing completed successfully!")
diff --git a/ffrprep/tests/test_cli.py b/ffrprep/tests/test_cli.py
index cddb143..fd20b59 100644
--- a/ffrprep/tests/test_cli.py
+++ b/ffrprep/tests/test_cli.py
@@ -347,6 +347,9 @@ def test_run_ffrprep_both_stages(
         mock_args.h_freq = None
         mock_args.n_procs = 1
         mock_args.no_report = False
+        mock_args.skip_existing = False
+        mock_args.clean_work_dir = False
+        mock_args.keep_epochs = True
 
         mock_parser.return_value.parse_args.return_value = mock_args
         mock_get_participants.return_value = ["01"]
@@ -537,6 +540,9 @@ def test_run_ffrprep_no_report_skips_html_reports(
         mock_args.h_freq = None
         mock_args.n_procs = 1
         mock_args.no_report = no_report
+        mock_args.skip_existing = False
+        mock_args.clean_work_dir = False
+        mock_args.keep_epochs = True
 
         mock_parser.return_value.parse_args.return_value = mock_args
         mock_get_participants.return_value = ["01"]
@@ -564,6 +570,182 @@ def test_run_ffrprep_no_report_skips_html_reports(
         assert mock_analysis_report.call_count == expected_calls
 
 
+def test_resume_and_cleanup_flags_default_to_previous_behaviour():
+    parser = get_parser()
+    args = parser.parse_args(["/path/to/bids", "/path/to/output", "participant"])
+    assert args.skip_existing is False
+    assert args.clean_work_dir is False
+    assert args.keep_epochs is True
+
+
+def test_resume_and_cleanup_flags_parse():
+    parser = get_parser()
+    args = parser.parse_args([
+        "/path/to/bids", "/path/to/output", "participant",
+        "--skip-existing", "--clean-work-dir", "--no-keep-epochs",
+    ])
+    assert args.skip_existing is True
+    assert args.clean_work_dir is True
+    assert args.keep_epochs is False
+
+
+def _resume_flow_args(tmp_dir, **overrides):
+    mock_args = MagicMock()
+    mock_args.bids_dir = Path(tmp_dir) / "bids"
+    mock_args.output_dir = Path(tmp_dir) / "output"
+    mock_args.analysis_level = "participant"
+    mock_args.stage = "both"
+    mock_args.skip_bids_validation = True
+    mock_args.participant_label = None
+    mock_args.baseline = "-0.2,0"
+    mock_args.ref_channels = "average"
+    mock_args.high_pass = 1.0
+    mock_args.low_pass = 40.0
+    mock_args.tmin = -0.2
+    mock_args.tmax = 0.6
+    mock_args.work_dir = None
+    mock_args.task = None
+    mock_args.run = None
+    mock_args.concat_runs = False
+    mock_args.no_filter = False
+    mock_args.no_auto_reject = False
+    mock_args.save_each_node = False
+    mock_args.events_file = None
+    mock_args.picks = None
+    mock_args.event_id = None
+    mock_args.reject_eeg = 75e-6
+    mock_args.on_missing = "warn"
+    mock_args.l_freq = None
+    mock_args.h_freq = None
+    mock_args.n_procs = 1
+    mock_args.no_report = False
+    mock_args.skip_existing = False
+    mock_args.clean_work_dir = False
+    mock_args.keep_epochs = True
+    for key, value in overrides.items():
+        setattr(mock_args, key, value)
+    return mock_args
+
+
+def _resume_flow_dirs(tmp_dir):
+    derivatives_root = Path(tmp_dir) / "derivatives"
+    preproc_dir = derivatives_root / "ffrprep-preprocessing" / "sub-01"
+    analysis_dir = derivatives_root / "ffrprep-analysis" / "sub-01"
+    preproc_dir.mkdir(parents=True)
+    analysis_dir.mkdir(parents=True)
+    return {
+        "derivatives_root": derivatives_root,
+        "preprocessing_dir": derivatives_root / "ffrprep-preprocessing",
+        "preprocessing_subject_dir": preproc_dir,
+        "analysis_dir": derivatives_root / "ffrprep-analysis",
+        "analysis_subject_dir": analysis_dir,
+    }
+
+
+@patch("ffrprep.ffrprep_cli._build_analysis_report")
+@patch("ffrprep.ffrprep_cli._build_preproc_report")
+@patch("ffrprep.ffrprep_cli._dispatch")
+@patch("ffrprep.ffrprep_cli.get_parser")
+@patch("ffrprep.ffrprep_cli.validate_input_dir")
+@patch("ffrprep.ffrprep_cli.get_participants")
+@patch("ffrprep.ffrprep_cli.get_sessions_tasks_runs")
+@patch("ffrprep.ffrprep_cli.setup_derivatives_directories")
+def test_run_ffrprep_skip_existing_with_removed_epochs_does_nothing(
+    mock_setup_dirs, mock_get_sessions, mock_get_participants, mock_validate,
+    mock_parser, mock_dispatch, mock_preproc_report, mock_analysis_report,
+    capsys,
+):
+    """Sidecars from an earlier --no-keep-epochs run mean nothing is left to do."""
+    with tempfile.TemporaryDirectory() as tmp_dir:
+        mock_parser.return_value.parse_args.return_value = _resume_flow_args(
+            tmp_dir, skip_existing=True,
+        )
+        mock_get_participants.return_value = ["01"]
+        mock_get_sessions.return_value = {
+            "sessions": [None], "tasks": ["passive"], "runs": [1],
+        }
+        dirs = _resume_flow_dirs(tmp_dir)
+        mock_setup_dirs.return_value = dirs
+        # Finished earlier: sidecars only (epochs were deleted).
+        (dirs["preprocessing_subject_dir"]
+         / "sub-01_task-passive_run-1_desc-preproc_epo.json").write_text("{}")
+        (dirs["analysis_subject_dir"]
+         / "sub-01_task-passive_run-1_desc-evoked.json").write_text("{}")
+
+        run_ffrprep()
+
+        mock_dispatch.assert_not_called()
+        mock_preproc_report.assert_not_called()
+        mock_analysis_report.assert_not_called()
+        out = capsys.readouterr().out
+        assert "1/1 preprocessing iteration(s) already complete" in out
+        assert "analysis outputs already complete" in out
+
+
+@patch("ffrprep.ffrprep_cli._build_analysis_report")
+@patch("ffrprep.ffrprep_cli._build_preproc_report")
+@patch("ffrprep.ffrprep_cli._dispatch")
+@patch("ffrprep.ffrprep_cli.get_parser")
+@patch("ffrprep.ffrprep_cli.validate_input_dir")
+@patch("ffrprep.ffrprep_cli.get_participants")
+@patch("ffrprep.ffrprep_cli.get_sessions_tasks_runs")
+@patch("ffrprep.ffrprep_cli.setup_derivatives_directories")
+def test_run_ffrprep_no_keep_epochs_removes_fif_after_analysis(
+    mock_setup_dirs, mock_get_sessions, mock_get_participants, mock_validate,
+    mock_parser, mock_dispatch, mock_preproc_report, mock_analysis_report,
+    capsys,
+):
+    with tempfile.TemporaryDirectory() as tmp_dir:
+        mock_parser.return_value.parse_args.return_value = _resume_flow_args(
+            tmp_dir, keep_epochs=False,
+        )
+        mock_get_participants.return_value = ["01"]
+        mock_get_sessions.return_value = {
+            "sessions": [None], "tasks": ["passive"], "runs": [1],
+        }
+        dirs = _resume_flow_dirs(tmp_dir)
+        mock_setup_dirs.return_value = dirs
+        fif = dirs["preprocessing_subject_dir"] / "sub-01_task-passive_run-1_desc-preproc_epo.fif"
+        sidecar = fif.with_suffix(".json")
+        fif.write_bytes(b"x" * 100)
+        sidecar.write_text("{}")
+
+        run_ffrprep()
+
+        assert not fif.exists()
+        assert sidecar.exists()
+        assert "--no-keep-epochs: removed 1 epochs file(s)" in capsys.readouterr().out
+
+
+@patch("ffrprep.ffrprep_cli._dispatch")
+@patch("ffrprep.ffrprep_cli.get_parser")
+@patch("ffrprep.ffrprep_cli.validate_input_dir")
+@patch("ffrprep.ffrprep_cli.get_participants")
+@patch("ffrprep.ffrprep_cli.get_sessions_tasks_runs")
+@patch("ffrprep.ffrprep_cli.setup_derivatives_directories")
+def test_run_ffrprep_no_keep_epochs_is_ignored_for_preprocessing_stage(
+    mock_setup_dirs, mock_get_sessions, mock_get_participants, mock_validate,
+    mock_parser, mock_dispatch, capsys,
+):
+    with tempfile.TemporaryDirectory() as tmp_dir:
+        mock_parser.return_value.parse_args.return_value = _resume_flow_args(
+            tmp_dir, stage="preprocessing", keep_epochs=False, no_report=True,
+        )
+        mock_get_participants.return_value = ["01"]
+        mock_get_sessions.return_value = {
+            "sessions": [None], "tasks": ["passive"], "runs": [1],
+        }
+        dirs = _resume_flow_dirs(tmp_dir)
+        mock_setup_dirs.return_value = dirs
+        fif = dirs["preprocessing_subject_dir"] / "sub-01_task-passive_run-1_desc-preproc_epo.fif"
+        fif.write_bytes(b"x")
+
+        run_ffrprep()
+
+        assert fif.exists()
+        assert "--no-keep-epochs is ignored with --stage preprocessing" in capsys.readouterr().out
+
+
 @patch("ffrprep.ffrprep_cli.get_parser")
 @patch("ffrprep.ffrprep_cli.get_participants")
 def test_run_ffrprep_no_participants(mock_get_participants, mock_parser):
diff --git a/ffrprep/tests/test_parallel_dispatch.py b/ffrprep/tests/test_parallel_dispatch.py
index af73e6b..1446b80 100644
--- a/ffrprep/tests/test_parallel_dispatch.py
+++ b/ffrprep/tests/test_parallel_dispatch.py
@@ -17,12 +17,22 @@
 
 import pytest
 
+from unittest.mock import MagicMock, patch
+
 from ffrprep.ffrprep_cli import (
     IterationResult,
+    _analysis_done,
+    _analysis_outputs_complete,
     _dispatch,
+    _find_preproc_outputs,
     _make_analysis_payload,
     _make_concat_payload,
     _make_preproc_payload,
+    _preproc_done,
+    _preproc_iteration,
+    _remove_epoch_files,
+    _remove_workflow_files,
+    _select_pending,
     _setup_worker_log,
     _snapshot_args,
     _snapshot_deriv,
@@ -266,3 +276,122 @@ def test_dispatch_propagates_worker_exception_across_processes():
     ]
     with pytest.raises(RuntimeError, match="injected failure for task-X"):
         _dispatch(_fail_task_X, payloads, n_procs=2, kind_label="test")
+
+
+# ----- resume (--skip-existing) and cleanup (--clean-work-dir / --no-keep-epochs) helpers
+
+
+def _preproc_payload(out_dir, kind="per_run", run="1", task="da"):
+    return {"kind": kind, "subject": "01", "task_label": task, "run_label": run,
+            "output_dir": str(out_dir), "work_dir": str(out_dir / "work"),
+            "identifier": f"task-{task}_run-{run}"}
+
+
+def _touch(path, size=0):
+    path.parent.mkdir(parents=True, exist_ok=True)
+    path.write_bytes(b"x" * size)
+    return path
+
+
+def test_find_preproc_outputs_filters_by_run_and_suffix(tmp_path):
+    _touch(tmp_path / "sub-01_task-da_run-1_desc-preprocA_epo.fif")
+    _touch(tmp_path / "sub-01_task-da_run-1_desc-preprocA_epo.json")
+    _touch(tmp_path / "sub-01_task-da_run-2_desc-preprocA_epo.json")
+    _touch(tmp_path / "sub-01_task-other_run-1_desc-preprocA_epo.json")
+
+    payload = _preproc_payload(tmp_path, run="1")
+    assert [p.name for p in _find_preproc_outputs(payload)] == ["sub-01_task-da_run-1_desc-preprocA_epo.fif"]
+    assert [p.name for p in _find_preproc_outputs(payload, suffix=".json")] == [
+        "sub-01_task-da_run-1_desc-preprocA_epo.json",
+    ]
+
+
+def test_preproc_done_uses_sidecars_so_it_survives_removed_epochs(tmp_path):
+    payload = _preproc_payload(tmp_path, run="1")
+    assert not _preproc_done(payload)
+    _touch(tmp_path / "sub-01_task-da_run-1_desc-preprocA_epo.json")   # no .fif: epochs were deleted
+    assert _preproc_done(payload)
+    assert not _preproc_done(_preproc_payload(tmp_path, run="2"))
+
+
+def test_preproc_done_concat_ignores_per_run_outputs(tmp_path):
+    payload = _preproc_payload(tmp_path, kind="concat", run=["1", "2"])
+    _touch(tmp_path / "sub-01_task-da_run-1_desc-preprocA_epo.json")
+    assert not _preproc_done(payload)
+    _touch(tmp_path / "sub-01_task-da_desc-preprocA_epo.json")         # no run token = concatenated
+    assert _preproc_done(payload)
+
+
+def test_analysis_done_and_outputs_complete(tmp_path):
+    analysis_dir = tmp_path / "analysis" / "sub-01"
+    payload = {
+        "preproc_files": [str(tmp_path / "sub-01_task-da_run-1_desc-preprocA_epo.fif")],
+        "analysis_subject_dir": str(analysis_dir),
+    }
+    assert not _analysis_done(payload)
+    assert not _analysis_outputs_complete(analysis_dir)
+    _touch(analysis_dir / "sub-01_task-da_run-1_desc-evokedA.json")    # per-type only: not "complete"
+    assert not _analysis_done(payload)
+    _touch(analysis_dir / "sub-01_task-da_run-1_desc-evoked.json")
+    assert _analysis_done(payload)
+    assert _analysis_outputs_complete(analysis_dir)
+
+
+def test_select_pending_only_filters_when_skip_existing(tmp_path, capsys):
+    done = _preproc_payload(tmp_path, run="1")
+    todo = _preproc_payload(tmp_path, run="2")
+    _touch(tmp_path / "sub-01_task-da_run-1_desc-preprocA_epo.json")
+
+    assert _select_pending([done, todo], _preproc_done, False, "preprocessing") == [done, todo]
+    assert capsys.readouterr().out == ""
+    assert _select_pending([done, todo], _preproc_done, True, "preprocessing") == [todo]
+    assert "1/2 preprocessing iteration(s) already complete" in capsys.readouterr().out
+
+
+def test_remove_epoch_files_deletes_fif_but_keeps_sidecars(tmp_path):
+    _touch(tmp_path / "a_epo.fif", size=1000)
+    _touch(tmp_path / "b_epo.fif", size=500)
+    _touch(tmp_path / "a_epo.json")
+    _touch(tmp_path / "keep.log")
+
+    assert _remove_epoch_files(tmp_path) == (2, 1500)
+    assert sorted(p.name for p in tmp_path.iterdir()) == ["a_epo.json", "keep.log"]
+    assert _remove_epoch_files(tmp_path) == (0, 0)
+
+
+def test_remove_workflow_files_keeps_logs_and_tolerates_missing_dirs(tmp_path):
+    _touch(tmp_path / "ffrprep_preproc" / "filter_data" / "result.pklz", size=10)
+    _touch(tmp_path / "iteration.log")
+    _remove_workflow_files(tmp_path, "ffrprep_preproc")
+    assert not (tmp_path / "ffrprep_preproc").exists()
+    assert (tmp_path / "iteration.log").exists()
+    _remove_workflow_files(tmp_path, "ffrprep_preproc")    # already gone: no error
+
+
+def _run_preproc_iteration(tmp_path, clean, check_raises=None):
+    payload = _preproc_payload(tmp_path)
+    payload["clean_work_dir"] = clean
+    work_dir = Path(payload["work_dir"])
+    _touch(work_dir / "ffrprep_preproc" / "epoch_data" / "result.pklz", size=10)
+    wf = MagicMock()
+    wf.name = "ffrprep_preproc"
+    with patch("ffrprep.ffrprep_cli._build_preproc_workflow", return_value=wf), \
+            patch("ffrprep.ffrprep_cli._check_preproc_output_or_raise",
+                  side_effect=check_raises, return_value=[]):
+        result = _preproc_iteration(payload)
+    return work_dir, result
+
+
+def test_preproc_iteration_cleans_work_dir_on_success_only_when_asked(tmp_path):
+    work_dir, _ = _run_preproc_iteration(tmp_path / "keep", clean=False)
+    assert (work_dir / "ffrprep_preproc").exists()
+
+    work_dir, result = _run_preproc_iteration(tmp_path / "clean", clean=True)
+    assert not (work_dir / "ffrprep_preproc").exists()
+    assert Path(result.log_path).exists()          # worker log survives the cleanup
+
+
+def test_preproc_iteration_keeps_work_dir_when_it_fails(tmp_path):
+    with pytest.raises(FileNotFoundError):
+        _run_preproc_iteration(tmp_path, clean=True, check_raises=FileNotFoundError("no outputs"))
+    assert (tmp_path / "work" / "ffrprep_preproc").exists()   # left for debugging

From 88bf7c267b34461b97a4938a98fbd530a131686c Mon Sep 17 00:00:00 2001
From: Kevin Sitek <sitek.kevin@gmail.com>
Date: Thu, 24 Sep 2026 16:11:08 -0500
Subject: [PATCH 07/10] Add QC flag columns to the group metrics table

--min-usable-pct and --min-snr add qc_usable_pct_ok / qc_snr_ok,
qc_include and qc_reason columns to the group metrics TSV. Recordings are
flagged, never dropped: the table and grand averages keep every subject so
downstream analyses choose whether to filter, and the report summary shows
the flagged count. A NaN value fails its check; --min-usable-pct without
--n-trials-presented is a clear error.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 CHANGES.rst                 |  5 +++
 docs/source/api_ref.rst     |  1 +
 docs/source/usage.rst       |  8 ++++
 ffrprep/ffrprep_cli.py      | 22 ++++++++++
 ffrprep/group.py            | 62 +++++++++++++++++++++++++-
 ffrprep/reports.py          |  6 ++-
 ffrprep/tests/test_cli.py   |  9 ++++
 ffrprep/tests/test_group.py | 86 +++++++++++++++++++++++++++++++++++++
 8 files changed, 195 insertions(+), 4 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 6a86b83..39c6ed2 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -81,6 +81,11 @@ Group-level analysis
   (and ``stim2resp_lim_*`` with ``--xcorr-lag-range``); ``--n-trials-presented``
   adds ``usable_pct``. ``discover_group_inputs`` now also returns the
   per-trial-type evoked files (``"by_type"``).
+- QC flags on the group metrics table: ``--min-usable-pct`` (needs
+  ``--n-trials-presented``) and ``--min-snr`` add ``qc_usable_pct_ok`` /
+  ``qc_snr_ok``, ``qc_include`` and ``qc_reason`` columns (``add_qc_flags``).
+  Subjects are flagged, never dropped: the table and the grand averages keep
+  every subject, and the report summary shows how many were flagged.
 - ``compute_grand_average`` renames single-channel evokeds that carry different
   channel names across sites (e.g. ``A32`` vs ``Cz``) to a common name before
   averaging; mismatched multi-channel sets raise a clear ``ValueError``.
diff --git a/docs/source/api_ref.rst b/docs/source/api_ref.rst
index 2817923..f07b1ae 100644
--- a/docs/source/api_ref.rst
+++ b/docs/source/api_ref.rst
@@ -112,6 +112,7 @@ Reference API
    discover_group_inputs
    compute_grand_average
    compute_subject_metrics
+   add_qc_flags
    merge_covariates
    save_group_outputs
    save_group_metrics
diff --git a/docs/source/usage.rst b/docs/source/usage.rst
index 41e8676..9ac53c0 100644
--- a/docs/source/usage.rst
+++ b/docs/source/usage.rst
@@ -283,6 +283,14 @@ table so downstream statistics have group, age, etc. alongside the metrics:
       --stimulus /path/to/da.wav --xcorr-stim-window 0.05 0.17 \
       --xcorr-resp-window 0.06 0.18 --n-trials-presented 6000
 
+``--min-usable-pct`` and ``--min-snr`` add quality-control flag columns
+(``qc_usable_pct_ok``, ``qc_snr_ok``, ``qc_include`` and a ``qc_reason`` text
+column) to the table. Subjects below a threshold are **flagged, not removed**:
+the table and the grand averages still contain them, and you decide
+downstream whether to filter on ``qc_include``. ``--min-usable-pct`` needs
+``--n-trials-presented``; ``--min-snr`` uses ``rms_snr_polarity_sum`` when
+``--f0`` is set and the combined-evoked ``rms_snr`` otherwise.
+
 This step only aggregates outputs participant-level ffrprep has
 already computed — it does not perform any group-level statistics
 (no hypothesis tests, no GLM). Downstream statistical analysis is
diff --git a/ffrprep/ffrprep_cli.py b/ffrprep/ffrprep_cli.py
index 94b64cf..3ef8f56 100644
--- a/ffrprep/ffrprep_cli.py
+++ b/ffrprep/ffrprep_cli.py
@@ -1988,6 +1988,28 @@ def get_parser():
             "presented) to the group metrics table."
         ),
     )
+    group_group.add_argument(
+        "--min-usable-pct",
+        dest="min_usable_pct",
+        type=float,
+        default=None,
+        help=(
+            "QC threshold on usable_pct (needs --n-trials-presented). Adds "
+            "qc_usable_pct_ok, qc_include and qc_reason columns; flagged "
+            "subjects are kept in the table and grand averages."
+        ),
+    )
+    group_group.add_argument(
+        "--min-snr",
+        dest="min_snr",
+        type=float,
+        default=None,
+        help=(
+            "QC threshold on RMS SNR (polarity-sum when available). Adds "
+            "qc_snr_ok, qc_include and qc_reason columns; flagged subjects "
+            "are kept in the table and grand averages."
+        ),
+    )
     group_group.add_argument(
         "--covariates",
         nargs="+",
diff --git a/ffrprep/group.py b/ffrprep/group.py
index 8266d76..3ec2561 100644
--- a/ffrprep/group.py
+++ b/ffrprep/group.py
@@ -410,6 +410,48 @@ def _polarity_sum_metrics(summed, harmonics, stimulus, response_window, stim_cac
     return columns
 
 
+def add_qc_flags(metrics, min_usable_pct=None, min_snr=None):
+    """Flag recordings that fail quality thresholds; never drop rows.
+
+    Adds boolean ``qc_usable_pct_ok`` (``usable_pct >= min_usable_pct``)
+    and/or ``qc_snr_ok`` (polarity-sum RMS SNR, else the combined-evoked
+    ``rms_snr``, ``>= min_snr``) columns for each threshold given, plus
+    ``qc_include`` (all checks pass) and ``qc_reason`` (the failed checks,
+    ``n/a`` when none). A missing value (NaN) fails its check. The metrics
+    table keeps every subject and the grand averages are unchanged, so
+    downstream analyses choose whether to filter on ``qc_include``.
+
+    Returns ``metrics`` unchanged when no threshold is given.
+    """
+    if min_usable_pct is None and min_snr is None:
+        return metrics
+
+    flagged = metrics.copy()
+    ok_columns = []
+    reasons = [[] for _ in range(len(flagged))]
+
+    def _check(column, threshold, ok_name, label):
+        values = flagged[column]
+        passed = (values >= threshold).fillna(False).astype(bool)
+        flagged[ok_name] = passed
+        ok_columns.append(ok_name)
+        for i, ok in enumerate(passed):
+            if not ok:
+                reasons[i].append(f"{label}<{threshold:g}")
+
+    if min_usable_pct is not None:
+        if "usable_pct" not in flagged:
+            raise ValueError("--min-usable-pct needs --n-trials-presented (no usable_pct column)")
+        _check("usable_pct", min_usable_pct, "qc_usable_pct_ok", "usable_pct")
+    if min_snr is not None:
+        snr_column = "rms_snr_polarity_sum" if "rms_snr_polarity_sum" in flagged else "rms_snr"
+        _check(snr_column, min_snr, "qc_snr_ok", snr_column)
+
+    flagged["qc_include"] = flagged[ok_columns].all(axis=1)
+    flagged["qc_reason"] = ["; ".join(r) or "n/a" for r in reasons]
+    return flagged
+
+
 def merge_covariates(metrics, covariate_paths):
     """Left-join subject-level covariate TSVs onto a metrics table.
 
@@ -536,6 +578,14 @@ def _write_group_dataset_description(group_dir):
         json.dump(dataset_desc, f, indent=2)
 
 
+def _qc_summary(metrics):
+    """Report summary entry for the QC flags (empty when no threshold was set)."""
+    if "qc_include" not in metrics:
+        return {}
+    n_flagged = int((~metrics["qc_include"]).sum())
+    return {"QC flagged": f"{n_flagged} of {len(metrics)}"}
+
+
 def _metric_options_from_args(args):
     """Translate group-level CLI args into ``compute_subject_metrics`` options.
 
@@ -606,7 +656,9 @@ def run_group_level(args):
     harmonic-amplitude and stimulus-to-response columns computed from the
     sum of the per-trial-type averages, and subject-level covariates from
     ``<bids_dir>/participants.tsv`` and ``args.covariates`` are joined onto
-    the saved TSV (never used for inference here).
+    the saved TSV (never used for inference here). ``args.min_usable_pct`` /
+    ``args.min_snr`` add ``qc_*`` flag columns (see :func:`add_qc_flags`);
+    flagged subjects stay in the table and the grand averages.
     """
     output_dir = Path(args.output_dir)
     response_window = tuple(getattr(args, "response_window", None) or (0.100, 0.200))
@@ -614,6 +666,10 @@ def run_group_level(args):
     task_filter = getattr(args, "task", None)
     metric_options = _metric_options_from_args(args)
     covariate_paths = _covariate_paths_from_args(args)
+    min_usable_pct = getattr(args, "min_usable_pct", None)
+    min_snr = getattr(args, "min_snr", None)
+    if min_usable_pct is not None and not metric_options["n_trials_presented"]:
+        raise ValueError("--min-usable-pct needs --n-trials-presented")
 
     inputs = discover_group_inputs(output_dir, participant_label=participant_label, task=task_filter)
 
@@ -643,6 +699,7 @@ def run_group_level(args):
             by_type_paths=inputs["by_type"].get((task, session, run), {}),
             **metric_options,
         )
+        metrics = add_qc_flags(metrics, min_usable_pct=min_usable_pct, min_snr=min_snr)
         metrics_with_covariates = merge_covariates(metrics, covariate_paths)
         metrics_path = save_group_metrics(output_dir, task, session, run, metrics_with_covariates)
         written.append(metrics_path)
@@ -656,9 +713,10 @@ def run_group_level(args):
                 response_window=response_window,
             ),
             reports.build_metrics_table_section(
-                metrics,
+                metrics.drop(columns=[c for c in metrics.columns if c.startswith("qc_")]),
                 section_id=f"{label}-metrics",
                 title="Per-subject metrics",
+                extra_summary=_qc_summary(metrics),
             ),
         ]
 
diff --git a/ffrprep/reports.py b/ffrprep/reports.py
index ffe44af..292614f 100644
--- a/ffrprep/reports.py
+++ b/ffrprep/reports.py
@@ -359,7 +359,7 @@ def build_evoked_section(
     }
 
 
-def build_metrics_table_section(metrics, *, section_id, title, max_table_rows=40):
+def build_metrics_table_section(metrics, *, section_id, title, max_table_rows=40, extra_summary=None):
     """Build a section descriptor summarizing a group-level metrics table.
 
     ``metrics`` is a per-(task, run) DataFrame such as
@@ -370,7 +370,8 @@ def build_metrics_table_section(metrics, *, section_id, title, max_table_rows=40
     full per-subject table is rendered as a figure; for larger cohorts
     (where a table would be unreadably tall) the figure shows one
     histogram per metric instead and the per-subject values live in the
-    accompanying ``_metrics.tsv``.
+    accompanying ``_metrics.tsv``. ``extra_summary`` (a dict of label to
+    text) adds entries to the summary table, e.g. the QC-flag count.
     """
     numeric_columns = [
         column for column in metrics.select_dtypes("number").columns if column != "n_avg"
@@ -381,6 +382,7 @@ def build_metrics_table_section(metrics, *, section_id, title, max_table_rows=40
         if values.empty:
             continue
         summary[column.replace("_", " ")] = f"{values.mean():.3g} \u00b1 {values.std():.3g}"
+    summary.update(extra_summary or {})
 
     if len(metrics) > max_table_rows and numeric_columns:
         n_panels = len(numeric_columns)
diff --git a/ffrprep/tests/test_cli.py b/ffrprep/tests/test_cli.py
index fd20b59..ffbd673 100644
--- a/ffrprep/tests/test_cli.py
+++ b/ffrprep/tests/test_cli.py
@@ -467,6 +467,15 @@ def test_group_level_option_defaults_leave_extra_metrics_off():
     assert args.covariates is None
 
 
+def test_qc_threshold_flags_default_off_and_parse():
+    parser = get_parser()
+    base = ["/path/to/bids", "/path/to/output", "group"]
+    args = parser.parse_args(base)
+    assert args.min_usable_pct is None and args.min_snr is None
+    args = parser.parse_args(base + ["--min-usable-pct", "60", "--min-snr", "1.5"])
+    assert args.min_usable_pct == 60.0 and args.min_snr == 1.5
+
+
 def test_group_level_options_parse_and_reach_group_module():
     from ffrprep.group import _covariate_paths_from_args, _metric_options_from_args
 
diff --git a/ffrprep/tests/test_group.py b/ffrprep/tests/test_group.py
index 6bddaeb..a0f91a4 100644
--- a/ffrprep/tests/test_group.py
+++ b/ffrprep/tests/test_group.py
@@ -12,6 +12,7 @@
 import pytest
 
 from ffrprep.group import (
+    add_qc_flags,
     compute_grand_average,
     compute_subject_metrics,
     discover_group_inputs,
@@ -541,3 +542,88 @@ def test_compute_grand_average_rejects_mismatched_multichannel_names(tmp_path):
 
     with pytest.raises(ValueError, match="identical channel sets"):
         compute_grand_average(paths)
+
+
+# ---------------------------------------------------------------------------
+# QC flag columns
+# ---------------------------------------------------------------------------
+
+def _qc_frame():
+    import pandas as pd
+
+    return pd.DataFrame({
+        "subject": ["a", "b", "c", "d"],
+        "usable_pct": [90.0, 60.0, 59.9, np.nan],
+        "rms_snr_polarity_sum": [3.0, 1.5, 1.49, 2.0],
+        "rms_snr": [9.0, 9.0, 9.0, 9.0],
+    })
+
+
+def test_add_qc_flags_is_a_no_op_without_thresholds():
+    frame = _qc_frame()
+    assert add_qc_flags(frame) is frame
+
+
+def test_add_qc_flags_thresholds_are_inclusive_and_nan_fails():
+    flagged = add_qc_flags(_qc_frame(), min_usable_pct=60, min_snr=1.5)
+
+    assert list(flagged["qc_usable_pct_ok"]) == [True, True, False, False]
+    assert list(flagged["qc_snr_ok"]) == [True, True, False, True]
+    assert list(flagged["qc_include"]) == [True, True, False, False]
+    assert list(flagged["qc_reason"]) == [
+        "n/a",
+        "n/a",
+        "usable_pct<60; rms_snr_polarity_sum<1.5",
+        "usable_pct<60",
+    ]
+    assert len(flagged) == 4  # rows are flagged, never dropped
+
+
+def test_add_qc_flags_single_threshold_adds_only_its_column():
+    flagged = add_qc_flags(_qc_frame(), min_snr=1.5)
+    assert "qc_snr_ok" in flagged and "qc_usable_pct_ok" not in flagged
+    assert list(flagged["qc_include"]) == [True, True, False, True]
+
+
+def test_add_qc_flags_falls_back_to_combined_snr_column():
+    frame = _qc_frame().drop(columns="rms_snr_polarity_sum")
+    flagged = add_qc_flags(frame, min_snr=10)
+    assert not flagged["qc_include"].any()
+    assert flagged["qc_reason"].iloc[0] == "rms_snr<10"
+
+
+def test_add_qc_flags_min_usable_pct_needs_usable_pct_column():
+    frame = _qc_frame().drop(columns="usable_pct")
+    with pytest.raises(ValueError, match="n-trials-presented"):
+        add_qc_flags(frame, min_usable_pct=60)
+
+
+def test_run_group_level_writes_qc_flag_columns_and_keeps_all_subjects(tmp_path):
+    import pandas as pd
+
+    out = tmp_path / "out"
+    for i, subject in enumerate(["01", "02", "03"], start=1):
+        _write_da_subject(out, subject, amplitude_uv=float(i), seed=i)
+    args = SimpleNamespace(
+        output_dir=out, participant_label=None, task=None, response_window=(0.0, 0.213),
+        n_trials_presented=6000, min_usable_pct=101, min_snr=None,  # nave=6000 -> usable 100%
+    )
+    run_group_level(args)
+
+    table = pd.read_csv(out / "ffrprep-group" / "task-da_run-01_metrics.tsv", sep="\t", dtype={"subject": str})
+    assert list(table["subject"]) == ["01", "02", "03"]
+    assert not table["qc_include"].any()
+    assert set(table["qc_reason"]) == {"usable_pct<101"}
+    assert "QC flagged" in (out / "ffrprep-group" / "group_report.html").read_text()
+
+
+def test_run_group_level_min_usable_pct_without_trials_presented_is_an_error(tmp_path):
+    out = tmp_path / "out"
+    for subject in ["01", "02"]:
+        _write_da_subject(out, subject, amplitude_uv=1.0, seed=int(subject))
+    args = SimpleNamespace(
+        output_dir=out, participant_label=None, task=None, response_window=(0.0, 0.213),
+        min_usable_pct=60,
+    )
+    with pytest.raises(ValueError, match="n-trials-presented"):
+        run_group_level(args)

From 7e6024eb76bd673e1732361cad678fa6ccc5a3d2 Mon Sep 17 00:00:00 2001
From: Kevin Sitek <sitek.kevin@gmail.com>
Date: Thu, 24 Sep 2026 16:12:24 -0500
Subject: [PATCH 08/10] Write a data dictionary beside the group metrics TSV

The group step now writes task-<task>[_run-<run>]_metrics.json with a
description and units for every metrics column (BIDS-style), embedding the
windows, harmonic settings and QC thresholds the run used. Covariate columns
take their description/levels/units from the JSON sidecar next to the
covariate TSV when present. merge_covariates gains return_sources and
save_group_metrics a dictionary argument.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 CHANGES.rst                 |   8 ++
 docs/source/api_ref.rst     |   1 +
 docs/source/usage.rst       |   6 ++
 ffrprep/group.py            | 174 ++++++++++++++++++++++++++++++++++--
 ffrprep/tests/test_group.py | 125 ++++++++++++++++++++++++++
 5 files changed, 306 insertions(+), 8 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 39c6ed2..db8e2fa 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -86,6 +86,14 @@ Group-level analysis
   ``qc_snr_ok``, ``qc_include`` and ``qc_reason`` columns (``add_qc_flags``).
   Subjects are flagged, never dropped: the table and the grand averages keep
   every subject, and the report summary shows how many were flagged.
+- The group metrics TSV now has a BIDS-style data dictionary,
+  ``task-<task>[_run-<run>]_metrics.json``, beside it (``build_metrics_dictionary``):
+  a description and units for every column, with the actual windows and
+  thresholds of the run embedded. Covariate columns take their description,
+  levels and units from the JSON sidecar next to the covariate TSV
+  (``participants.json``, ``phenotype/*.json``) when one exists.
+  ``merge_covariates`` gains ``return_sources`` and ``save_group_metrics``
+  a ``dictionary`` argument.
 - ``compute_grand_average`` renames single-channel evokeds that carry different
   channel names across sites (e.g. ``A32`` vs ``Cz``) to a common name before
   averaging; mismatched multi-channel sets raise a clear ``ValueError``.
diff --git a/docs/source/api_ref.rst b/docs/source/api_ref.rst
index f07b1ae..d0503c3 100644
--- a/docs/source/api_ref.rst
+++ b/docs/source/api_ref.rst
@@ -113,6 +113,7 @@ Reference API
    compute_grand_average
    compute_subject_metrics
    add_qc_flags
+   build_metrics_dictionary
    merge_covariates
    save_group_outputs
    save_group_metrics
diff --git a/docs/source/usage.rst b/docs/source/usage.rst
index 9ac53c0..10e4d20 100644
--- a/docs/source/usage.rst
+++ b/docs/source/usage.rst
@@ -291,6 +291,12 @@ downstream whether to filter on ``qc_include``. ``--min-usable-pct`` needs
 ``--n-trials-presented``; ``--min-snr`` uses ``rms_snr_polarity_sum`` when
 ``--f0`` is set and the combined-evoked ``rms_snr`` otherwise.
 
+Next to ``*_metrics.tsv`` the group step writes ``*_metrics.json``, a
+BIDS-style data dictionary with a description and units for every column. It
+records the windows, harmonic settings and QC thresholds the run used, and
+copies the descriptions of joined covariates from the ``.json`` next to their
+TSV (for example ``participants.json``) when there is one.
+
 This step only aggregates outputs participant-level ffrprep has
 already computed — it does not perform any group-level statistics
 (no hypothesis tests, no GLM). Downstream statistical analysis is
diff --git a/ffrprep/group.py b/ffrprep/group.py
index 3ec2561..da06f1d 100644
--- a/ffrprep/group.py
+++ b/ffrprep/group.py
@@ -452,7 +452,7 @@ def _check(column, threshold, ok_name, label):
     return flagged
 
 
-def merge_covariates(metrics, covariate_paths):
+def merge_covariates(metrics, covariate_paths, return_sources=False):
     """Left-join subject-level covariate TSVs onto a metrics table.
 
     Each TSV needs a ``participant_id`` column (``sub-<label>``, as in a
@@ -463,9 +463,12 @@ def merge_covariates(metrics, covariate_paths):
     Returns
     -------
     pandas.DataFrame
-        ``metrics`` with the covariate columns appended.
+        ``metrics`` with the covariate columns appended. With
+        ``return_sources=True`` a ``(table, sources)`` tuple, where
+        ``sources`` maps each appended column to the file it came from.
     """
     merged = metrics.copy()
+    sources = {}
     key = "sub-" + merged["subject"].astype(str)
     for path in covariate_paths or []:
         path = Path(path)
@@ -482,7 +485,144 @@ def merge_covariates(metrics, covariate_paths):
         lookup = table.drop_duplicates("participant_id").set_index("participant_id")[new_columns]
         for column in new_columns:
             merged[column] = key.map(lookup[column])
-    return merged
+            sources[column] = path
+    return (merged, sources) if return_sources else merged
+
+
+def _covariate_sidecar_entry(path, column):
+    """The ``column`` entry of the JSON sidecar next to a covariate TSV, or ``None``."""
+    sidecar = Path(path).with_suffix(".json")
+    if not sidecar.exists():
+        return None
+    try:
+        entry = json.loads(sidecar.read_text()).get(column)
+    except (OSError, ValueError):
+        return None
+    return entry if isinstance(entry, dict) else None
+
+
+def build_metrics_dictionary(
+    columns,
+    *,
+    response_window=(0.100, 0.200),
+    harmonics=None,
+    stimulus=None,
+    n_trials_presented=None,
+    min_usable_pct=None,
+    min_snr=None,
+    covariate_sources=None,
+):
+    """Build a BIDS-style data dictionary for the group metrics TSV.
+
+    Parameters
+    ----------
+    columns : sequence of str
+        Column names of the saved metrics table, in order.
+    response_window, harmonics, stimulus, n_trials_presented, min_usable_pct, min_snr
+        The options the metrics were computed with (as returned by
+        ``_metric_options_from_args``); descriptions embed the actual
+        windows and thresholds so the dictionary documents the run.
+    covariate_sources : dict, optional
+        ``{column: path}`` from ``merge_covariates(..., return_sources=True)``.
+        A covariate's description (and levels/units) is copied from the JSON
+        sidecar next to its TSV (``participants.json``, ``phenotype/*.json``)
+        when present.
+
+    Returns
+    -------
+    dict
+        ``{column: {"Description": ..., "Units": ...}}`` with one entry per
+        column, in the order of ``columns``. Unknown columns get a stub entry.
+    """
+    resp = f"{response_window[0] * 1e3:g}-{response_window[1] * 1e3:g} ms"
+    info = {
+        "subject": {"Description": "Participant label (without the 'sub-' prefix)"},
+        "n_avg": {"Description": "Number of trials averaged in the combined evoked (kept after rejection)",
+                  "Units": "trials"},
+        "usable_pct": {
+            "Description": f"100 * n_avg / {n_trials_presented} trials presented"
+            if n_trials_presented else "100 * n_avg / trials presented",
+            "Units": "%",
+        },
+        "rms_snr": {
+            "Description": f"RMS of the response window ({resp}) divided by RMS of the baseline, "
+                           "combined evoked response",
+            "Units": "ratio",
+        },
+        "band_power_90_110hz": {
+            "Description": f"Mean 90-110 Hz power in the response window ({resp}), combined evoked response",
+            "Units": "V^2",
+        },
+        "response_consistency": {
+            "Description": "Mean pairwise trial-to-trial correlation of the preprocessed epochs "
+                           "(averaged over runs; empty when a recording has fewer than 10 epochs)",
+            "Units": "r",
+        },
+        "rms_snr_polarity_sum": {
+            "Description": f"RMS of the response window ({resp}) divided by RMS of the baseline, "
+                           "sum of the two per-trial-type averages",
+            "Units": "ratio",
+        },
+        "qc_usable_pct_ok": {
+            "Description": f"usable_pct >= {min_usable_pct:g}" if min_usable_pct is not None
+            else "usable_pct passes the threshold",
+            "Levels": {"True": "passes", "False": "fails or missing"},
+        },
+        "qc_snr_ok": {
+            "Description": f"SNR (rms_snr_polarity_sum when present, else rms_snr) >= {min_snr:g}"
+            if min_snr is not None else "SNR passes the threshold",
+            "Levels": {"True": "passes", "False": "fails or missing"},
+        },
+        "qc_include": {
+            "Description": "All requested QC checks pass; subjects are flagged, never removed",
+            "Levels": {"True": "include", "False": "flagged"},
+        },
+        "qc_reason": {"Description": "Failed QC checks ('n/a' when none)"},
+    }
+    if harmonics:
+        window = f"{harmonics['tmin'] * 1e3:g}-{harmonics['tmax'] * 1e3:g} ms"
+        band = f"+/-{harmonics['bin_hz'] / 2:g} Hz"
+        info["f0_uv"] = {
+            "Description": f"Mean FFT amplitude (2|FFT|/L) within {band} of {harmonics['f0']:g} Hz, "
+                           f"{window}, sum of the two per-trial-type averages",
+            "Units": "uV",
+        }
+        info["upper_harmonics_uv"] = {
+            "Description": f"Sum over harmonics 2..{harmonics['n_harmonics']} of the mean FFT amplitude "
+                           f"within {band} of k*{harmonics['f0']:g} Hz, {window}, sum of the two "
+                           "per-trial-type averages",
+            "Units": "uV",
+        }
+    if stimulus:
+        stim_win = f"{stimulus['stim_window'][0] * 1e3:g}-{stimulus['stim_window'][1] * 1e3:g} ms"
+        resp_win = f"{stimulus['resp_window'][0] * 1e3:g}-{stimulus['resp_window'][1] * 1e3:g} ms"
+        base = (f"cross-correlation ('coeff' normalization) of the stimulus ({stim_win}) with the "
+                f"polarity-summed response ({resp_win})")
+        info["stim2resp_r"] = {"Description": f"Maximum {base} over all lags", "Units": "r"}
+        info["stim2resp_z"] = {"Description": "Fisher z transform (arctanh) of stim2resp_r", "Units": "z"}
+        info["stim2resp_lag_ms"] = {
+            "Description": "Lag of the maximum stimulus-to-response correlation "
+                           "(positive = response delayed)", "Units": "ms",
+        }
+        lag_range = stimulus.get("lag_range_ms")
+        if lag_range:
+            lim = f"restricted to lags of {lag_range[0]:g} to {lag_range[1]:g} ms"
+            info["stim2resp_lim_r"] = {"Description": f"Maximum {base}, {lim}", "Units": "r"}
+            info["stim2resp_lim_z"] = {"Description": "Fisher z transform of stim2resp_lim_r", "Units": "z"}
+            info["stim2resp_lim_lag_ms"] = {
+                "Description": f"Lag of the maximum lag-restricted correlation ({lim})", "Units": "ms",
+            }
+
+    dictionary = {}
+    for column in columns:
+        entry = info.get(column)
+        if entry is None and column in (covariate_sources or {}):
+            path = covariate_sources[column]
+            entry = _covariate_sidecar_entry(path, column) or {
+                "Description": f"Joined from {Path(path).name}",
+            }
+        dictionary[column] = entry or {"Description": "No description available"}
+    return dictionary
 
 
 def _sort_key(key):
@@ -550,8 +690,12 @@ def save_group_outputs(derivatives_root, task, session, run, grand_average, subj
     return evoked_path
 
 
-def save_group_metrics(derivatives_root, task, session, run, metrics):
-    """Save a per-subject metrics table to a group-level TSV file."""
+def save_group_metrics(derivatives_root, task, session, run, metrics, dictionary=None):
+    """Save a per-subject metrics table to a group-level TSV file.
+
+    When ``dictionary`` (see :func:`build_metrics_dictionary`) is given it
+    is written as ``<base>_metrics.json`` next to the TSV.
+    """
     derivatives_root = Path(derivatives_root)
     group_dir = derivatives_root / "ffrprep-group"
     group_dir.mkdir(parents=True, exist_ok=True)
@@ -560,6 +704,8 @@ def save_group_metrics(derivatives_root, task, session, run, metrics):
     base = _basename(task, session, run)
     metrics_path = group_dir / f"{base}_metrics.tsv"
     metrics.to_csv(metrics_path, sep="\t", index=False)
+    if dictionary is not None:
+        metrics_path.with_suffix(".json").write_text(json.dumps(dictionary, indent=2) + "\n")
     return metrics_path
 
 
@@ -700,9 +846,21 @@ def run_group_level(args):
             **metric_options,
         )
         metrics = add_qc_flags(metrics, min_usable_pct=min_usable_pct, min_snr=min_snr)
-        metrics_with_covariates = merge_covariates(metrics, covariate_paths)
-        metrics_path = save_group_metrics(output_dir, task, session, run, metrics_with_covariates)
-        written.append(metrics_path)
+        metrics_with_covariates, covariate_sources = merge_covariates(
+            metrics, covariate_paths, return_sources=True,
+        )
+        dictionary = build_metrics_dictionary(
+            metrics_with_covariates.columns,
+            response_window=response_window,
+            min_usable_pct=min_usable_pct,
+            min_snr=min_snr,
+            covariate_sources=covariate_sources,
+            **metric_options,
+        )
+        metrics_path = save_group_metrics(
+            output_dir, task, session, run, metrics_with_covariates, dictionary=dictionary,
+        )
+        written.extend([metrics_path, metrics_path.with_suffix(".json")])
 
         sections = [
             reports.build_evoked_section(
diff --git a/ffrprep/tests/test_group.py b/ffrprep/tests/test_group.py
index a0f91a4..d84c9d5 100644
--- a/ffrprep/tests/test_group.py
+++ b/ffrprep/tests/test_group.py
@@ -13,6 +13,7 @@
 
 from ffrprep.group import (
     add_qc_flags,
+    build_metrics_dictionary,
     compute_grand_average,
     compute_subject_metrics,
     discover_group_inputs,
@@ -627,3 +628,127 @@ def test_run_group_level_min_usable_pct_without_trials_presented_is_an_error(tmp
     )
     with pytest.raises(ValueError, match="n-trials-presented"):
         run_group_level(args)
+
+
+# ---------------------------------------------------------------------------
+# metrics data dictionary
+# ---------------------------------------------------------------------------
+
+_HARMONICS = {"f0": 100.0, "n_harmonics": 10, "bin_hz": 60.0, "tmin": 0.06, "tmax": 0.18}
+_STIMULUS = {"path": "da.wav", "stim_window": (0.05, 0.17), "resp_window": (0.06, 0.18)}
+
+
+def test_metrics_dictionary_has_one_entry_per_column_in_order():
+    columns = ["subject", "n_avg", "usable_pct", "rms_snr", "band_power_90_110hz", "response_consistency",
+               "rms_snr_polarity_sum", "f0_uv", "upper_harmonics_uv", "stim2resp_r", "stim2resp_z",
+               "stim2resp_lag_ms", "qc_usable_pct_ok", "qc_include", "qc_reason", "mystery"]
+    dictionary = build_metrics_dictionary(
+        columns, response_window=(0.0, 0.213), harmonics=_HARMONICS, stimulus=_STIMULUS,
+        n_trials_presented=6000, min_usable_pct=60,
+    )
+    assert list(dictionary) == columns
+    assert all("Description" in entry for entry in dictionary.values())
+    assert dictionary["f0_uv"]["Units"] == "uV"
+    assert dictionary["usable_pct"]["Units"] == "%"
+    assert dictionary["stim2resp_lag_ms"]["Units"] == "ms"
+    assert dictionary["mystery"]["Description"] == "No description available"
+
+
+def test_metrics_dictionary_descriptions_embed_the_run_parameters():
+    dictionary = build_metrics_dictionary(
+        ["usable_pct", "rms_snr", "f0_uv", "upper_harmonics_uv", "stim2resp_r", "stim2resp_lim_r",
+         "qc_usable_pct_ok", "qc_snr_ok"],
+        response_window=(0.0, 0.213), harmonics=_HARMONICS,
+        stimulus={**_STIMULUS, "lag_range_ms": (6.9, 10.9)},
+        n_trials_presented=6000, min_usable_pct=60, min_snr=1.5,
+    )
+    assert "6000 trials presented" in dictionary["usable_pct"]["Description"]
+    assert "0-213 ms" in dictionary["rms_snr"]["Description"]
+    assert "+/-30 Hz of 100 Hz" in dictionary["f0_uv"]["Description"]
+    assert "60-180 ms" in dictionary["f0_uv"]["Description"]
+    assert "2..10" in dictionary["upper_harmonics_uv"]["Description"]
+    assert "50-170 ms" in dictionary["stim2resp_r"]["Description"]
+    assert "6.9 to 10.9 ms" in dictionary["stim2resp_lim_r"]["Description"]
+    assert dictionary["qc_usable_pct_ok"]["Description"] == "usable_pct >= 60"
+    assert dictionary["qc_snr_ok"]["Description"].endswith(">= 1.5")
+
+
+def test_metrics_dictionary_omits_optional_metrics_that_were_off():
+    dictionary = build_metrics_dictionary(["subject", "f0_uv"], response_window=(0.1, 0.2))
+    assert dictionary["f0_uv"]["Description"] == "No description available"
+
+
+def test_metrics_dictionary_copies_covariate_descriptions_from_sibling_json(tmp_path):
+    (tmp_path / "participants.tsv").write_text("participant_id\tgroup\tage\n")
+    (tmp_path / "participants.json").write_text(json.dumps({
+        "group": {"Description": "Musicianship group", "Levels": {"Mus": "musician"}},
+    }))
+    dictionary = build_metrics_dictionary(
+        ["group", "age"],
+        covariate_sources={"group": tmp_path / "participants.tsv", "age": tmp_path / "participants.tsv"},
+    )
+    assert dictionary["group"] == {"Description": "Musicianship group", "Levels": {"Mus": "musician"}}
+    assert dictionary["age"] == {"Description": "Joined from participants.tsv"}
+
+
+def test_merge_covariates_reports_the_source_file_of_each_new_column(tmp_path):
+    import pandas as pd
+
+    from ffrprep.group import merge_covariates
+
+    metrics = pd.DataFrame({"subject": ["01", "02"], "age": [1, 2]})
+    first = tmp_path / "participants.tsv"
+    first.write_text("participant_id\tage\tgroup\nsub-01\t9\tMus\nsub-02\t9\tNMus\n")
+    second = tmp_path / "pheno.tsv"
+    second.write_text("participant_id\tgroup\tdp\nsub-01\tX\t1.5\nsub-02\tX\t0.5\n")
+
+    merged, sources = merge_covariates(metrics, [first, second], return_sources=True)
+
+    assert sources == {"group": first, "dp": second}
+    assert list(merged.columns) == ["subject", "age", "group", "dp"]
+    assert merge_covariates(metrics, [first]).equals(merged[["subject", "age", "group"]])
+
+
+def test_save_group_metrics_writes_dictionary_beside_the_tsv(tmp_path):
+    import pandas as pd
+
+    metrics = pd.DataFrame({"subject": ["01"], "n_avg": [5]})
+    dictionary = {"subject": {"Description": "s"}, "n_avg": {"Description": "n", "Units": "trials"}}
+    path = save_group_metrics(tmp_path, "da", None, "01", metrics, dictionary=dictionary)
+
+    assert path.name == "task-da_run-01_metrics.tsv"
+    assert json.loads(path.with_suffix(".json").read_text()) == dictionary
+    plain = save_group_metrics(tmp_path / "other", "da", None, "01", metrics)
+    assert not plain.with_suffix(".json").exists()
+
+
+def test_run_group_level_writes_dictionary_covering_every_tsv_column(tmp_path):
+    import pandas as pd
+
+    bids_dir = tmp_path / "bids"
+    bids_dir.mkdir()
+    (bids_dir / "participants.tsv").write_text(
+        "participant_id\tgroup\nsub-01\tMus\nsub-02\tNMus\nsub-03\tVar\n"
+    )
+    (bids_dir / "participants.json").write_text(
+        json.dumps({"group": {"Description": "Musicianship group"}})
+    )
+    out = tmp_path / "out"
+    for i, subject in enumerate(["01", "02", "03"], start=1):
+        _write_da_subject(out, subject, amplitude_uv=float(i), seed=i)
+    stim_path = tmp_path / "stim.wav"
+    _write_stimulus_wav(stim_path)
+    args = SimpleNamespace(
+        output_dir=out, bids_dir=bids_dir, participant_label=None, task=None,
+        response_window=(0.0, 0.213), f0=100.0, stimulus=stim_path,
+        n_trials_presented=6000, min_usable_pct=60, min_snr=1.5,
+    )
+    written = run_group_level(args)
+
+    group_dir = out / "ffrprep-group"
+    table = pd.read_csv(group_dir / "task-da_run-01_metrics.tsv", sep="\t")
+    dictionary = json.loads((group_dir / "task-da_run-01_metrics.json").read_text())
+    assert list(dictionary) == list(table.columns)
+    assert dictionary["group"]["Description"] == "Musicianship group"
+    assert "No description available" not in {v["Description"] for v in dictionary.values()}
+    assert group_dir / "task-da_run-01_metrics.json" in written

From 55365ef00a7713ba5785309b5c941e8186402a37 Mon Sep 17 00:00:00 2001
From: Kevin Sitek <sitek.kevin@gmail.com>
Date: Thu, 24 Sep 2026 16:20:05 -0500
Subject: [PATCH 09/10] Avoid an RST emphasis warning in the --keep-epochs help

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 ffrprep/ffrprep_cli.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ffrprep/ffrprep_cli.py b/ffrprep/ffrprep_cli.py
index 3ef8f56..bfbeeb7 100644
--- a/ffrprep/ffrprep_cli.py
+++ b/ffrprep/ffrprep_cli.py
@@ -2067,7 +2067,7 @@ def get_parser():
         action=argparse.BooleanOptionalAction,
         default=True,
         help=(
-            "Keep the large preprocessed *_epo.fif files after the analysis "
+            "Keep the large preprocessed epochs (_epo.fif) files after the analysis "
             "stage. --no-keep-epochs deletes them (JSON sidecars stay) once "
             "analysis and its report are done; ignored with --stage "
             "preprocessing because analysis still needs them."

From de611bc33fc01b81ddc84cb8b8c007288873dcda Mon Sep 17 00:00:00 2001
From: Kevin Sitek <sitek.kevin@gmail.com>
Date: Thu, 24 Sep 2026 16:57:57 -0500
Subject: [PATCH 10/10] Update docs for the group step, resume/cleanup flags
 and sidecar fields

pipeline_details.rst: the group section now covers the polarity-sum and
stimulus metrics, covariates, QC flags, the metrics dictionary, channel-name
harmonization and the ffrprep-group output files; the sidecar text lists
RejectionMode and Sources/RawSources and explains the epoch counts; the
workflow section describes --skip-existing, --clean-work-dir and
--no-keep-epochs. The walkthrough, index and README point to the group step.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 CHANGES.rst                      |  4 ++
 README.rst                       |  5 +-
 docs/source/index.rst            |  4 +-
 docs/source/pipeline_details.rst | 81 +++++++++++++++++++++++++++-----
 docs/source/walkthrough.rst      |  3 ++
 5 files changed, 84 insertions(+), 13 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index db8e2fa..481ca9a 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -86,6 +86,10 @@ Group-level analysis
   ``qc_snr_ok``, ``qc_include`` and ``qc_reason`` columns (``add_qc_flags``).
   Subjects are flagged, never dropped: the table and the grand averages keep
   every subject, and the report summary shows how many were flagged.
+- Docs: the pipeline details page now covers the group-level metrics, covariates,
+  QC flags and output files, the ``--skip-existing`` / ``--clean-work-dir`` /
+  ``--no-keep-epochs`` behaviour, and the ``RejectionMode`` and ``Sources``
+  sidecar fields.
 - The group metrics TSV now has a BIDS-style data dictionary,
   ``task-<task>[_run-<run>]_metrics.json``, beside it (``build_metrics_dictionary``):
   a description and units for every column, with the actual windows and
diff --git a/README.rst b/README.rst
index fdb4e14..b4e49ce 100644
--- a/README.rst
+++ b/README.rst
@@ -51,7 +51,10 @@ under your scheduler of choice.
 Running the ``group`` analysis level against an ``output_dir`` that
 already holds participant-level derivatives aggregates them across
 subjects: a grand-average evoked response, a per-subject metrics
-table, and a group HTML report per (task, run). This is aggregation
+table, and a group HTML report per (task, run). The table can also
+hold F0 and harmonic amplitudes, stimulus-to-response correlation,
+covariates from ``participants.tsv`` and QC flag columns, and comes
+with a ``_metrics.json`` data dictionary. This is aggregation
 only, in the spirit of a "prep" tool — no group-level statistics are
 computed; see the `Usage guide
 <https://sitek.github.io/ffrprep/usage.html>`_ for details.
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 67580a0..a55f69a 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -67,7 +67,9 @@ parallel workers per subject) with the published Docker image:
 Outputs land under
 ``/path/to/your/bids_dataset/derivatives/ffrprep-preprocessing/``
 and ``…/ffrprep-analysis/``, each carrying a
-``sub-XX_<stage>_report.html`` per subject. See
+``sub-XX_<stage>_report.html`` per subject. Running the ``group``
+analysis level on the same output directory afterwards writes a grand
+average and a per-subject metrics table to ``…/ffrprep-group/``. See
 :ref:`installation` for image pinning and a Singularity build, and
 :ref:`usage` for the full list of CLI flags plus worked examples
 (per-trial-type analysis, explicit difference pairs, cluster
diff --git a/docs/source/pipeline_details.rst b/docs/source/pipeline_details.rst
index ab2c1a1..9ede345 100644
--- a/docs/source/pipeline_details.rst
+++ b/docs/source/pipeline_details.rst
@@ -188,9 +188,11 @@ combined file under ``--no-split-by-trial-type``.
      doesn't leak through ``.save()``).
    - Write each sibling ``.json`` sidecar with ``EpochCount`` /
      ``EpochCountTotal`` / ``EpochCountRejected`` /
-     ``RejectionThresholds`` / ``Filtering`` / ``SamplingFrequency``
-     / ``EpochTmin`` / ``EpochTmax`` / ``Channels`` plus run /
-     session / ``ConcatenatedRuns`` / ``Condition`` provenance.
+     ``RejectionMode`` / ``RejectionThresholds`` / ``Filtering`` /
+     ``SamplingFrequency`` / ``EpochTmin`` / ``EpochTmax`` /
+     ``Channels`` plus run / session / ``ConcatenatedRuns`` /
+     ``Condition`` provenance. ``Sources`` and ``RawSources`` list the
+     raw EEG file(s) the epochs came from, when they can be found.
    - Initialize the per-derivatives ``dataset_description.json`` if
      missing.
 
@@ -213,10 +215,16 @@ Under ``--no-split-by-trial-type`` the per-trial-type files are
 replaced by a single ``_desc-preproc_epo.fif`` + its sidecar.
 
 The sidecar JSON carries provenance, ``EpochCount`` /
-``EpochCountTotal`` / ``EpochCountRejected``, ``RejectionThresholds``,
-``Filtering`` (high-pass and low-pass cut-offs), sampling frequency,
-run / session identifiers, and ``Condition`` (for per-trial-type
-files only).
+``EpochCountTotal`` / ``EpochCountRejected``, ``RejectionMode``
+(``peak-to-peak`` or ``absolute-amplitude``) with
+``RejectionThresholds``, ``Filtering`` (high-pass and low-pass
+cut-offs), sampling frequency, run / session identifiers, the raw
+file(s) in ``Sources`` / ``RawSources``, and ``Condition`` (for
+per-trial-type files only). The epoch counts cover only the trial
+types being written: in a per-trial-type file they are that trial
+type's counts, and events outside ``event_id`` (for example markers)
+are not counted. Epochs dropped by ``--reject-mode abs`` count as
+rejected.
 
 *Outputs:* File paths, processing metadata.
 
@@ -410,6 +418,20 @@ Pipeline Integration and Quality Control
 - nipype caches per-iteration intermediates under ``work/`` so
   re-runs that already have a saved ``_desc-preproc_epo.fif`` skip
   the workflow re-execution
+- ``--skip-existing`` skips a (task, run) iteration whose sidecar
+  outputs already exist, so an interrupted run can be repeated with
+  the same command. The check uses the ``.json`` sidecars rather than
+  the epochs files, so it also works after ``--no-keep-epochs``. An
+  analysis iteration counts as done when its ``_desc-evoked.json``
+  exists, and no report is rebuilt for iterations that are skipped.
+- ``--clean-work-dir`` deletes an iteration's Nipype working files
+  once it has succeeded. Log files are kept, and the files of a
+  failed iteration are left in place.
+- ``--no-keep-epochs`` deletes the ``_epo.fif`` files after the
+  analysis stage and its report have finished. The JSON sidecars
+  stay. With ``--stage preprocessing`` the flag is ignored and a
+  warning is printed, because the analysis stage still needs the
+  epochs.
 
 **Quality Control Checkpoints:**
 
@@ -449,18 +471,55 @@ session, run) with derivatives from at least two subjects, it:
   present, difference (``_desc-evokedDiff{A}Vs{B}.fif``) Evoked files
   via :py:func:`ffrprep.group.discover_group_inputs`, and computes a
   grand average with ``mne.grand_average``
-  (:py:func:`ffrprep.group.compute_grand_average`);
+  (:py:func:`ffrprep.group.compute_grand_average`). If the subjects'
+  single channel has different names (for example ``A32`` and ``Cz``
+  at different sites), it is renamed to a common name first.
+  Multi-channel sets that differ raise an error;
 - recomputes the same scalar FFR metrics the participant-level
   report shows (RMS SNR, band power) directly from each subject's
   saved combined Evoked, plus trial-to-trial response consistency
   from that subject's saved preprocessing Epochs
   (:py:func:`ffrprep.group.compute_subject_metrics`) — nothing is
   recomputed from raw data;
-- writes the grand-average Evoked(s), a per-subject metrics TSV, and
-  a single ``group_report.html`` (reusing the same section builders
-  and template as the participant-level reports) under
+- adds further columns when asked. ``--f0`` adds
+  ``rms_snr_polarity_sum``, ``f0_uv`` and ``upper_harmonics_uv``;
+  ``--stimulus`` adds ``stim2resp_r``, ``stim2resp_z`` and
+  ``stim2resp_lag_ms`` (plus ``stim2resp_lim_*`` with
+  ``--xcorr-lag-range``); ``--n-trials-presented`` adds
+  ``usable_pct``. These are computed from the sum of the two
+  per-trial-type evoked files. ``f0_uv`` and ``upper_harmonics_uv``
+  are FFT amplitudes averaged in a band around each harmonic of F0
+  (:py:func:`ffrprep.analysis.harmonic_amplitudes`), and
+  ``stim2resp_r`` is the maximum normalized cross-correlation between
+  the stimulus and the response
+  (:py:func:`ffrprep.analysis.stim_to_resp_xcorr`);
+- joins subject-level covariates onto the table by ``participant_id``,
+  from ``<bids_dir>/participants.tsv`` and any ``--covariates`` files
+  (:py:func:`ffrprep.group.merge_covariates`). A column already in
+  the table is not overwritten;
+- with ``--min-usable-pct`` and/or ``--min-snr``, adds ``qc_*``
+  columns that mark recordings below a threshold
+  (:py:func:`ffrprep.group.add_qc_flags`). Rows are flagged, not
+  removed, and the grand averages still include every subject;
+- writes the grand-average Evoked(s), a per-subject metrics TSV with a
+  ``_metrics.json`` data dictionary next to it
+  (:py:func:`ffrprep.group.build_metrics_dictionary`), and a single
+  ``group_report.html`` (reusing the same section builders and
+  template as the participant-level reports) under
   ``output_dir/ffrprep-group/``.
 
+*Output Structure:* ::
+
+    derivatives/ffrprep-group/
+    ├── dataset_description.json
+    ├── task-YY_run-ZZ_desc-grandAverage_ave.fif
+    ├── task-YY_run-ZZ_desc-grandAverage_ave.json
+    ├── task-YY_run-ZZ_desc-grandAverageDiffPositiveVsNegative_ave.fif
+    ├── task-YY_run-ZZ_desc-grandAverageDiffPositiveVsNegative_ave.json
+    ├── task-YY_run-ZZ_metrics.tsv
+    ├── task-YY_run-ZZ_metrics.json
+    └── group_report.html
+
 This step is deliberately scoped to aggregation, not inference: it
 summarizes what participant-level ffrprep already computed and
 performs no group-level statistics (no hypothesis tests, no GLM).
diff --git a/docs/source/walkthrough.rst b/docs/source/walkthrough.rst
index 2232fd4..0be20d5 100644
--- a/docs/source/walkthrough.rst
+++ b/docs/source/walkthrough.rst
@@ -455,6 +455,9 @@ After completing this tutorial, you can:
 - Process your own FFR datasets using the same workflow
 - Modify preprocessing parameters for your specific experimental setup
 - Use the generated outputs for further statistical analysis
+- Aggregate several subjects with the ``group`` analysis level, which
+  writes a grand average and a per-subject metrics table to
+  ``derivatives/ffrprep-group/`` (Example 5 in :ref:`usage`)
 - Integrate ``ffrprep`` into automated processing pipelines
 
 For more advanced usage, see the :ref:`usage` documentation and `API reference <api_ref.html>`_.