Reproduce Whiteford paper results with new preproc and analysis

CHANGES.rst

@@ -19,13 +19,29 @@ 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``.
 - Fixed: ``--reject-eeg 0`` now disables automatic rejection as its help
   text documents (``parse_reject``). It previously built a ``{"eeg": 0.0}``
   threshold, which rejects every epoch. Negative thresholds are now an
   argument error.
 - 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
@@ -57,6 +73,38 @@ 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"``).
+- 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.
+- 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
+  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``.
+- 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``
@@ -149,6 +197,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

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.

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:
@@ -108,6 +112,9 @@ Reference API
    discover_group_inputs
    compute_grand_average
    compute_subject_metrics
+   add_qc_flags
+   build_metrics_dictionary
+   merge_covariates
    save_group_outputs
    save_group_metrics
    run_group_level

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

docs/source/pipeline_details.rst

@@ -132,7 +132,9 @@ Preprocessing Workflow Nodes
 *Default Parameters:*
    - **High-pass:** 70.0 Hz (removes slow drifts and low-frequency noise, preserves FFR frequencies)
    - **Low-pass:** 1000.0 Hz (preserves FFR harmonics while removing 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.
 
@@ -185,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.
 
@@ -210,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.
 
@@ -407,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:**
 
@@ -446,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).

docs/source/usage.rst

@@ -265,6 +265,38 @@ 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
+
+``--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.
+
+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
@@ -335,6 +367,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
 ----------------
 

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>`_.