Trial-to-trial reliability metrics for cue and lever epochs

CHANGELOG.md

@@ -4,6 +4,10 @@ Versions follow [Semantic Versioning](https://semver.org) (`<major>.<minor>.<pat
 
 ## Unreleased
 
+### Added
+
+- Trial-to-trial reliability metrics in the `<stem>_metrics-session.csv` written by `process metrics` and `process peri-event --with-metrics`, for go/no-go sessions: `epoch_correlation`, `response_fraction`, `variance_quench`, `signal_fraction` and `reliability_n`, for all trials and per `sdt_type` (`<metric>_sdt-<type>`) and go/no-go stimulus (`<metric>_stim-go`, `<metric>_stim-nogo`). Cue-aligned epochs end at each trial's response (`--no-mask-response` to disable), and lever-aligned epochs run from the cue to the push, with the baseline before the cue (`--cue-baseline`). Options `--min-rt`, `--response-sd` and `--time-warp`. The existing columns are unchanged. The peri-event report's session table shows the new columns with a trial-group selector ([#173](https://github.com/DuguidLab/mesoscopy/issues/173)).
+
 ### Fixed
 
 - `process peri-event --event response` wrongly aligned windows to `start_time + response_time`, which falls before the cue since `response_time` is measured from the cue and every trial starts with the ITI. The lever push is now taken at `cue_onset + response_time` ([#176](https://github.com/DuguidLab/mesoscopy/issues/176)).

docs/how-to/metrics.md

@@ -45,6 +45,9 @@ Pearson correlation between trial traces over the response window) and, for each
 `<metric>_sd` and `<metric>_cv` across trials, ignoring trials where the metric is undefined. `onset_n`, `decay_n`
 and `offset_n` count the trials where each was found.
 
+For go/no-go sessions aligned to the cue, trial start or lever push, the session table also has the
+[reliability](#reliability) columns.
+
 Any metric is empty when it is undefined for that trial, for example a decay when the trace never falls back to
 the fraction, or an onset when the trace never crosses the threshold.
 
@@ -85,6 +88,47 @@ estimate that the `sd` onset threshold depends on.
 - Session CV is unstable for any metric whose mean sits near zero, which the extrapolated onset often does.
   Read the SD in that case.
 
+## Reliability
+
+How consistent a region's response is from trial to trial, over the part of each trial before the lever push, so
+that movement does not inflate it. Taken for go/no-go sessions only, and needs the `cue_onset`, `response_time` and
+`sdt_type` trials columns (carried over by `process peri-event`, or joined with `--trials`). Without them, or for
+reward-aligned windows, the columns are skipped with a warning.
+
+**Epoch.** The event the windows are aligned to is read from which trials column matches `event_time`.
+
+- Cue- or trial-start-aligned: the response window, ending at each trial's response. Misses and correct rejections
+  end at the median response time of the trials with a response. `--no-mask-response` uses the whole response
+  window.
+- Lever-aligned (`--event response`): from the cue to the lever push. The baseline is `--cue-baseline START END`
+  (default `-1 0`) before each trial's cue, so extract these windows with a `--pre` long enough to reach it (e.g.
+  `--pre 5`); trials whose baseline falls outside the window have no response fraction or variance quench
+  baseline, and a warning gives their count.
+
+Trials with a response time below `--min-rt` (default 0.2 s) are left out.
+
+**Metrics.**
+
+| Column | Meaning |
+| --- | --- |
+| `epoch_correlation` | Mean zero-lag Pearson correlation between every pair of trials, over the samples both have. How alike single trials are. |
+| `response_fraction` | Fraction of trials whose mean epoch value exceeds their baseline mean by `--response-sd` (default 2) baseline SDs. |
+| `variance_quench` | Across-trial variance over the epoch divided by that over the baseline. Below 1 when variability drops after the event. |
+| `signal_fraction` | Fraction of single-trial variance explained by the trial-mean trace. |
+| `reliability_n` | Trials used. |
+
+**Groups.** Each metric is taken over all trials (`<metric>`), per trial type (`<metric>_sdt-hit`,
+`_sdt-miss`, `_sdt-false_alarm`, `_sdt-correct_rejection`) and per stimulus (`<metric>_stim-go` for hits and
+misses, `<metric>_stim-nogo` for false alarms and correct rejections). Every column is always written; a group with
+fewer than 10 trials is empty.
+
+**Trial counts.** None of the metrics depends on the number of trials, but all are noisier with fewer, so read
+them alongside `reliability_n` when comparing groups.
+
+**Time warping.** Epochs end at different times. By default, trials are compared sample by sample on the shared
+time axis. `--time-warp` instead resamples each epoch onto a common grid from its start to its end, comparing the
+shape of the response regardless of response time.
+
 ## As a library
 
 ```python
@@ -97,4 +141,6 @@ per_trial, per_session = metrics.metrics_tables(perievent, smoothing=5)
 
 `metrics.trial_metrics` works on a `(n_trials, n_samples)` array for one region, and the per-metric functions
 (`peak`, `auc`, `onset_time`, `extrapolated_onset`, `offset_time`, `decay_time`, `trace_correlation`) are
-available on their own.
+available on their own. `metrics.reliability_metrics` takes the reliability columns for one region, with options
+in a `metrics.ReliabilityOptions`; `epoch_correlation`, `response_fraction`, `variance_quench` and
+`signal_fraction` work on masked `(n_trials, n_samples)` arrays from `reliability_epochs` and `epoch_traces`.

docs/how-to/reports.md

@@ -33,5 +33,6 @@ panel in the grid.
 When the `*_metrics.csv` and `*_metrics-session.csv` from `process metrics` sit next to the input, the
 report adds a metrics section: onset, peak and offset markers on the traces with the across-trial spread of
 each, per-trial box plots for the selected region, the atlas coloured by a session metric, and the session
-table. Peri-event files written without the trials columns have no trial types; pass the trials CSV with
+table. For go/no-go sessions the table also shows the reliability metrics, with a selector for the trial group
+they are taken over. Peri-event files written without the trials columns have no trial types; pass the trials CSV with
 `-t/--trials` to recover them.

docs/typical-workflow.md

@@ -96,7 +96,9 @@ mesoscopy process metrics /path/to/example-recording_regions_event-cueonset_peri
 writes `example-recording_regions_event-cueonset_metrics.csv` with the onset time, peak time, amplitude, area
 under the curve, decay time, offset time and duration of every trial in every region, with the trials columns
 carried over from the peri-event file, and `example-recording_regions_event-cueonset_metrics-session.csv` with the mean,
-SD and CV of each metric across trials per region, and the mean pairwise correlation between trial traces.
+SD and CV of each metric across trials per region, and the mean pairwise correlation between trial traces. For
+go/no-go sessions the session table also has trial-to-trial reliability metrics over the cue-to-lever epoch, for all
+trials, per trial type and per stimulus; see [Response metrics](how-to/metrics.md#reliability).
 
 Traces are baseline-subtracted with the mean over `--baseline START END` (default: all pre-event samples), and
 metrics are taken over `--response START END` (default: all post-event samples). Onset is the first

src/mesoscopy/process/__init__.py

@@ -23,8 +23,11 @@
 
 from __future__ import annotations
 
+import dataclasses
+import functools
 import os
 import typing
+import warnings
 from pathlib import Path
 
 import click
@@ -490,6 +493,74 @@ def _metrics_options(command: Callable[..., None]) -> Callable[..., None]:
     return command
 
 
+def _reliability_options(command: Callable[..., None]) -> Callable[..., None]:
+    """Reliability metric options, passed to the callback as one `reliability` argument.
+
+    Args:
+        command (Callable[..., None]): Command callback to decorate; takes a `reliability` keyword argument.
+
+    Returns:
+        Callable[..., None]: The callback with the options attached.
+    """
+    names = [field.name for field in dataclasses.fields(pm.ReliabilityOptions) if field.name != "min_trials"]
+
+    @functools.wraps(command)
+    def wrapper(*args: typing.Any, **kwargs: typing.Any) -> None:
+        kwargs["reliability"] = pm.ReliabilityOptions(**{name: kwargs.pop(name) for name in names})
+        command(*args, **kwargs)
+
+    options = [
+        click.option(
+            "--mask-response/--no-mask-response",
+            default=True,
+            show_default=True,
+            help="End cue-aligned reliability epochs at each trial's response; trials without one end at the"
+            " median response time.",
+        ),
+        click.option(
+            "--min-rt",
+            type=click.FloatRange(min=0),
+            default=0.2,
+            show_default=True,
+            help="Drop trials with a response time below this, in seconds, from the reliability metrics.",
+        ),
+        click.option(
+            "--response-sd",
+            type=float,
+            default=2.0,
+            show_default=True,
+            help="Baseline SD multiple a trial's mean epoch response must exceed for response_fraction.",
+        ),
+        click.option(
+            "--time-warp",
+            is_flag=True,
+            default=False,
+            help="Resample each reliability epoch onto a common 0-1 grid.",
+        ),
+        click.option(
+            "--cue-baseline",
+            type=(float, float),
+            default=(-1.0, 0.0),
+            show_default=True,
+            help="Baseline window START END relative to each trial's cue, for lever-aligned reliability metrics.",
+        ),
+    ]
+    for option in reversed(options):
+        wrapper = option(wrapper)
+    return wrapper
+
+
+def _echo_warnings(caught: list[warnings.WarningMessage]) -> None:
+    """Echo the user warnings recorded while computing metrics.
+
+    Args:
+        caught (list[warnings.WarningMessage]): Warnings recorded by `warnings.catch_warnings(record=True)`.
+    """
+    for warning in caught:
+        if issubclass(warning.category, UserWarning):
+            click.echo(f"Warning: {warning.message}")
+
+
 def _validate_metrics_options(smoothing: int, extrapolate_range: tuple[float, float]) -> None:
     """Reject metric options that `trial_metrics` would refuse.
 
@@ -589,6 +660,7 @@ def _write_metrics(trial_table: pd.DataFrame, session_table: pd.DataFrame, out_d
     help="Also write the response metrics of `process metrics`, joined with TRIALS_PATH. CSV input only.",
 )
 @_metrics_options
+@_reliability_options
 def perievent_cmd(
     recording_path: str,
     trials_path: str,
@@ -608,6 +680,7 @@ def perievent_cmd(
     extrapolate_range: tuple[float, float],
     smoothing: int,
     decay_fraction: float,
+    reliability: pm.ReliabilityOptions,
 ) -> None:
     """Extract per-trial windows around a behavioural event from a behaviour-aligned recording.
 
@@ -660,7 +733,8 @@ def perievent_cmd(
     if with_metrics and n_kept == 0:
         click.echo("No trials kept, skipping metrics.")
     elif with_metrics:
-        with timer.Timer(message="Extracting metrics"):
+        with timer.Timer(message="Extracting metrics"), warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("always", UserWarning)
             tables = pm.metrics_tables(
                 windows,
                 baseline=baseline,
@@ -672,7 +746,10 @@ def perievent_cmd(
                 extrapolate_range=extrapolate_range,
                 decay_fraction=decay_fraction,
                 smoothing=smoothing,
+                event=event,
+                reliability=reliability,
             )
+        _echo_warnings(caught)
         _write_metrics(*tables, out_dir, Path(outpath).stem.removesuffix("_perievent"))
 
 
@@ -839,6 +916,7 @@ def _perievent_csv(
     help="Trials CSV to join onto the per-trial table by trial_index, for peri-event files without trials columns.",
 )
 @_metrics_options
+@_reliability_options
 def metrics_cmd(
     path: str,
     out_dir: str,
@@ -852,13 +930,16 @@ def metrics_cmd(
     extrapolate_range: tuple[float, float],
     smoothing: int,
     decay_fraction: float,
+    reliability: pm.ReliabilityOptions,
 ) -> None:
     """Extract per-trial response metrics and per-session variability from peri-event traces.
 
     PATH is the long-format *_perievent.csv written by `process peri-event`. Writes <stem>_metrics.csv with onset
     time, peak time, amplitude, area under the curve, decay time, offset time and duration per trial per region,
     and <stem>_metrics-session.csv with the mean, SD and CV of each across trials per region, plus the mean
-    pairwise trial-trace correlation.
+    pairwise trial-trace correlation. For go/no-go sessions aligned to the cue, trial start or lever push, the
+    session table also has trial-to-trial reliability metrics over the cue-to-lever epoch, for all trials, per
+    sdt_type and per go/no-go stimulus.
     """  # noqa: DOC501
     import pandas as pd
 
@@ -878,7 +959,8 @@ def metrics_cmd(
         click.echo(f"Loading trials from {trials_path}...")
         trials = pd.read_csv(trials_path)
 
-    with timer.Timer(message="Extracting metrics"):
+    with timer.Timer(message="Extracting metrics"), warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter("always", UserWarning)
         try:
             tables = pm.metrics_tables(
                 perievent,
@@ -892,9 +974,11 @@ def metrics_cmd(
                 extrapolate_range=extrapolate_range,
                 decay_fraction=decay_fraction,
                 smoothing=smoothing,
+                reliability=reliability,
             )
         except ValueError as error:
             msg = f"{path}: {error} Expected the columns written by `process peri-event`."
             raise click.ClickException(msg) from error
+    _echo_warnings(caught)
 
     _write_metrics(*tables, out_dir, Path(path).stem.removesuffix("_perievent"))