From c89dcb097f148f5e82eb772f8b5529c709a80e15 Mon Sep 17 00:00:00 2001 From: gbradham Date: Thu, 29 Oct 2026 13:16:53 -0500 Subject: [PATCH 1/4] feat(review): assign sessions to groups from the session pool Select sessions in Session Review's pool (Ctrl/Shift-click), right-click, and move them to a group, a new group, or out of one. A cohort opened as a folder saves the groups to glider_project.json, keeping its subjects and other session fields; a broken manifest is never saved over. Sessions opened as files are grouped for the review only. --- src/glider/gui/behavior/analysis_window.py | 58 +++++++++++- src/glider/gui/review/pool.py | 66 ++++++++++++- .../unit/gui/behavior/test_analysis_window.py | 92 +++++++++++++++++++ tests/unit/gui/review/test_pool.py | 48 ++++++++++ 4 files changed, 261 insertions(+), 3 deletions(-) diff --git a/src/glider/gui/behavior/analysis_window.py b/src/glider/gui/behavior/analysis_window.py index b04285cc..7cb1ecdb 100644 --- a/src/glider/gui/behavior/analysis_window.py +++ b/src/glider/gui/behavior/analysis_window.py @@ -332,6 +332,8 @@ def __init__(self, parent=None): self._heatmap_grid: tuple[np.ndarray, np.ndarray, np.ndarray] | None = None self._cohort: list[tuple[Path, SessionView]] = [] self._ids: list[str] = [] + # The ids before _unique: what glider_project.json knows them by. + self._manifest_ids: list[str] = [] self._groups: list[str] = [] self._timelines: list = [] self._recordings_of: list = [] @@ -375,6 +377,7 @@ def __init__(self, parent=None): # focus after a click -- those keys are the whole point of the window. self._pool.tree.setFocusPolicy(Qt.FocusPolicy.NoFocus) self._pool.session_picked.connect(self._on_session_picked) + self._pool.group_assigned.connect(self.set_groups) self._canvas = KeypointCanvas() self._transport = Transport() viewer = QFrame() @@ -986,7 +989,7 @@ def _set_cohort(self, loaded, *, ids=None, groups=None) -> None: loaded = list(loaded) ids = list(ids) if ids is not None else [session_id_for(p) for p, _ in loaded] groups = list(groups) if groups is not None else [""] * len(loaded) - ids = _unique(ids) + manifest_ids, ids = ids, _unique(ids) built = [self._timeline_for(path, view) for path, view in loaded] # Before any assignment too: an unreadable file pops a modal dialog, # whose event loop must still see the old cohort whole. @@ -996,6 +999,7 @@ def _set_cohort(self, loaded, *, ids=None, groups=None) -> None: else self._open_store(self._cohort_root / mk.COHORT_FILE, t0=None) ) self._cohort, self._ids, self._groups = loaded, ids, groups + self._manifest_ids = manifest_ids self._invalidate_cohort_cache() self._marker_stores = {} self._retrieval_state, self._retrievals = {}, {} @@ -1005,6 +1009,58 @@ def _set_cohort(self, loaded, *, ids=None, groups=None) -> None: self._fill_pool(0) self._show_session(0) + def set_groups(self, indices: list[int], group: str) -> None: + """Put the sessions at ``indices`` in ``group`` ("" for none). + + A cohort opened as a folder keeps it in ``glider_project.json``, so the + next open, the MCP tools and anything else reading the manifest agree. + Sessions opened as files have no folder to hold one: the grouping + lasts until they are closed, and the status bar says so. + """ + group = group.strip() + chosen = {i for i in indices if 0 <= i < len(self._cohort)} + if not chosen: + return + # Two recordings sharing a name share one manifest record: grouping + # one groups both, on disk and so on screen. + names = {self._manifest_ids[i] for i in chosen} + chosen |= {i for i, sid in enumerate(self._manifest_ids) if sid in names} + if self._cohort_root is not None: + if not self._save_groups(sorted(names), group): + return + else: + self.statusBar().showMessage( + "Groups set for this review only: open the sessions as a folder to save them.", + 8000, + ) + for i in chosen: + self._groups[i] = group + # Every range's rows carry the group. + self._invalidate_cohort_cache() + self._fill_pool(self._shown) + self._fill_epochs() + + def _save_groups(self, session_ids: list[str], group: str) -> bool: + """Write the groups into the cohort's manifest; False if it could not be.""" + from glider.core.project import Project, ProjectError + + try: + project = Project.load(self._cohort_root) + for sid in session_ids: + project.set_session(sid, group=group) + path = project.save() + except ProjectError as e: + # Never saved over: the manifest holds the subject mapping too. + QMessageBox.warning( + self, "Set group", f"{e}\n\nFix or remove it, then set the groups again." + ) + return False + except OSError as e: + QMessageBox.warning(self, "Set group", f"Could not save the groups: {e}") + return False + self.statusBar().showMessage(f"Groups saved to {path.name}", 5000) + return True + def _fill_pool(self, selected: int) -> None: order = behavior_order(label for _, view in self._cohort for label in view.labels) self._pool.set_entries( diff --git a/src/glider/gui/review/pool.py b/src/glider/gui/review/pool.py index 09899a31..3ed4d516 100644 --- a/src/glider/gui/review/pool.py +++ b/src/glider/gui/review/pool.py @@ -3,7 +3,8 @@ A cohort is the unit of analysis, so every loaded session stays visible, grouped by the treatment recorded in ``glider_project.json``. Each row says which evidence the session has -- video, poses, hardware -- before you click -it, and carries its ethogram in miniature. +it, and carries its ethogram in miniature. Right-clicking a selection puts it +in a group; the window decides where that is written. """ from __future__ import annotations @@ -11,14 +12,17 @@ from dataclasses import dataclass import numpy as np -from PyQt6.QtCore import QSize, Qt, pyqtSignal +from PyQt6.QtCore import QPoint, QSize, Qt, pyqtSignal from PyQt6.QtGui import QColor, QImage, QPixmap from PyQt6.QtWidgets import ( + QAbstractItemView, QFrame, QGridLayout, QHBoxLayout, + QInputDialog, QLabel, QLineEdit, + QMenu, QTreeWidget, QTreeWidgetItem, QVBoxLayout, @@ -113,6 +117,8 @@ class SessionPool(QFrame): """Every loaded session, grouped by treatment; clicking one shows it.""" session_picked = pyqtSignal(int) + #: Session indices and the group to put them in; "" takes them out of one. + group_assigned = pyqtSignal(list, str) def __init__(self, parent=None): super().__init__(parent) @@ -143,9 +149,15 @@ def __init__(self, parent=None): self.tree.setHeaderHidden(True) self.tree.setRootIsDecorated(False) self.tree.setIndentation(0) + # Ctrl- and Shift-click gather animals to group; the current row is + # still the one shown. + self.tree.setSelectionMode(QAbstractItemView.SelectionMode.ExtendedSelection) + self.tree.setContextMenuPolicy(Qt.ContextMenuPolicy.CustomContextMenu) + self.tree.customContextMenuRequested.connect(self._context_menu) column.addWidget(self.tree, 1) self._items: list[QTreeWidgetItem] = [] self._names: list[str] = [] + self._entry_groups: list[str] = [] self._groups: list[QTreeWidgetItem] = [] self.filter.textChanged.connect(self._apply_filter) self.tree.currentItemChanged.connect(self._on_current) @@ -154,6 +166,7 @@ def set_entries(self, entries: list[PoolEntry]) -> None: self.tree.blockSignals(True) self.tree.clear() self._items, self._names, self._groups = [], [], [] + self._entry_groups = [e.group for e in entries] # PoolEntry holds a QImage, which is unhashable, so group names are # collected, not entries. grouped = bool(entries) and list(dict.fromkeys(e.group for e in entries)) != [""] @@ -194,6 +207,55 @@ def current(self) -> int: index = None if item is None else item.data(0, _INDEX) return -1 if index is None else int(index) + def selected_indices(self) -> list[int]: + """The selected sessions, in pool order; group headers are not sessions.""" + picked = {item.data(0, _INDEX) for item in self.tree.selectedItems()} + return sorted(int(i) for i in picked if i is not None) + + def group_names(self) -> list[str]: + """Every group a loaded session is in, in the order the pool shows them.""" + return [g for g in dict.fromkeys(self._entry_groups) if g] + + def group_menu(self, indices: list[int]) -> QMenu: + """Put ``indices`` in an existing group, a new one, or none.""" + menu = QMenu(self) + count = len(indices) + title = menu.addAction(f"{count} session{'s' if count != 1 else ''}") + title.setEnabled(False) + current = {self._entry_groups[i] for i in indices} + for name in self.group_names(): + action = menu.addAction(f"Move to {name}") + # Already there: nothing to do. + action.setEnabled(current != {name}) + action.triggered.connect(lambda _=False, g=name: self.group_assigned.emit(indices, g)) + menu.addSeparator() + new = menu.addAction("New group…") + new.triggered.connect(lambda: self._new_group(indices)) + clear = menu.addAction("Remove from group") + clear.setEnabled(current != {""}) + clear.triggered.connect(lambda: self.group_assigned.emit(indices, "")) + return menu + + def _ask_group_name(self) -> str: + name, ok = QInputDialog.getText(self, "New group", "Group name:") + return name.strip() if ok else "" + + def _new_group(self, indices: list[int]) -> None: + name = self._ask_group_name() + if name: + self.group_assigned.emit(indices, name) + + def _context_menu(self, pos: QPoint) -> None: + item = self.tree.itemAt(pos) + index = None if item is None else item.data(0, _INDEX) + if index is None: + return + indices = self.selected_indices() + # A right-click outside the selection means that row alone. + if int(index) not in indices: + indices = [int(index)] + self.group_menu(indices).exec(self.tree.viewport().mapToGlobal(pos)) + def _on_current(self, item, _previous) -> None: index = None if item is None else item.data(0, _INDEX) if index is not None: diff --git a/tests/unit/gui/behavior/test_analysis_window.py b/tests/unit/gui/behavior/test_analysis_window.py index 4ef27a6f..ac95d66e 100644 --- a/tests/unit/gui/behavior/test_analysis_window.py +++ b/tests/unit/gui/behavior/test_analysis_window.py @@ -1419,6 +1419,98 @@ def test_a_broken_manifest_warns_and_still_loads(self, qtbot, tmp_path, monkeypa assert warned and "ungrouped" in warned[0] +class TestSettingGroups: + def _cohort(self, tmp_path, manifest=None): + for sid in ("s1", "s2", "s3"): + _ethogram(tmp_path / "sessions" / sid / "analysis", ["groom"] * 300) + if manifest is not None: + (tmp_path / "glider_project.json").write_text(manifest, encoding="utf-8") + + def _open(self, qtbot, tmp_path): + win = AnalysisWindow() + qtbot.addWidget(win) + win.load_folder(tmp_path) + return win + + def test_a_folder_without_a_manifest_gets_one(self, qtbot, tmp_path): + from glider.core.project import Project + + self._cohort(tmp_path) + win = self._open(qtbot, tmp_path) + win.set_groups([0, 2], "ChR2") + assert win._groups == ["ChR2", "", "ChR2"] + assert Project.load(tmp_path).groups() == {"ChR2": ["s1", "s3"]} + assert win._pool.tree.topLevelItemCount() == 2 # ChR2, Ungrouped + + def test_the_rest_of_the_manifest_is_kept(self, qtbot, tmp_path): + from glider.core.project import Project + + self._cohort( + tmp_path, + '{"subjects": {"m1": {"sex": "F"}}, ' + '"sessions": {"s1": {"subject": "m1", "group": "eYFP", "day": "1"}}}', + ) + win = self._open(qtbot, tmp_path) + win.set_groups([0], "ChR2") + project = Project.load(tmp_path) + record = project.record("s1") + assert (record.group, record.subject, record.day) == ("ChR2", "m1", "1") + assert project.subjects["m1"].sex == "F" + + def test_reopening_shows_the_groups_set(self, qtbot, tmp_path): + self._cohort(tmp_path) + self._open(qtbot, tmp_path).set_groups([1], "Saline") + assert self._open(qtbot, tmp_path)._groups == ["", "Saline", ""] + + def test_removing_from_a_group(self, qtbot, tmp_path): + from glider.core.project import Project + + self._cohort(tmp_path, '{"sessions": {"s1": {"group": "eYFP"}}}') + win = self._open(qtbot, tmp_path) + win.set_groups([0], "") + assert win._groups == ["", "", ""] + assert Project.load(tmp_path).groups() == {} + + def test_the_epoch_rows_carry_the_new_group(self, qtbot, tmp_path): + self._cohort(tmp_path) + win = self._open(qtbot, tmp_path) + before = win.range_rows(0.0, 5.0) + win.set_groups([0], "ChR2") + after = win.range_rows(0.0, 5.0) + assert [r["group"] for r in before] == ["", "", ""] + assert [r["group"] for r in after] == ["ChR2", "", ""] + + def test_a_broken_manifest_is_never_saved_over(self, qtbot, tmp_path, monkeypatch): + self._cohort(tmp_path) + win = self._open(qtbot, tmp_path) + (tmp_path / "glider_project.json").write_text("{not json", encoding="utf-8") + warned = [] + monkeypatch.setattr( + "glider.gui.behavior.analysis_window.QMessageBox.warning", + lambda *a, **k: warned.append(a[-1]), + ) + win.set_groups([0], "ChR2") + assert warned and "not valid JSON" in warned[0] + assert win._groups == ["", "", ""] + assert (tmp_path / "glider_project.json").read_text(encoding="utf-8") == "{not json" + + def test_sessions_opened_as_files_are_grouped_for_now(self, qtbot, tmp_path): + self._cohort(tmp_path) + win = AnalysisWindow() + qtbot.addWidget(win) + win.load_many(sorted(tmp_path.rglob("ethogram_raw.csv"))) + win.set_groups([1], "ChR2") + assert win._groups[1] == "ChR2" + assert not (tmp_path / "glider_project.json").exists() + assert "this review only" in win.statusBar().currentMessage() + + def test_the_pool_menu_reaches_the_window(self, qtbot, tmp_path): + self._cohort(tmp_path) + win = self._open(qtbot, tmp_path) + win._pool.group_assigned.emit([0, 1], "ChR2") + assert win._groups == ["ChR2", "ChR2", ""] + + class TestEditingKeys: def _win(self, qtbot, tmp_path): win = AnalysisWindow() diff --git a/tests/unit/gui/review/test_pool.py b/tests/unit/gui/review/test_pool.py index 1ae3f6bc..be501e45 100644 --- a/tests/unit/gui/review/test_pool.py +++ b/tests/unit/gui/review/test_pool.py @@ -93,3 +93,51 @@ def test_a_narrow_pool_shows_every_badge_whole(qtbot): for badge in badges: assert badge.isVisible() assert badge.sizeHint().width() <= badge.width(), badge.text() + + +class TestGrouping: + def _actions(self, menu): + return {a.text(): a for a in menu.actions() if a.text()} + + def test_several_rows_can_be_selected(self, pool): + pool._items[0].setSelected(True) + pool._items[2].setSelected(True) + assert pool.selected_indices() == [0, 2] + + def test_group_headers_are_not_sessions(self, pool): + pool.tree.topLevelItem(0).setSelected(True) + assert pool.selected_indices() == [] + + def test_the_menu_offers_every_group_and_says_how_many(self, pool): + actions = self._actions(pool.group_menu([0, 2])) + assert "2 sessions" in actions + assert {"Move to ChR2", "Move to eYFP", "New group…", "Remove from group"} <= set(actions) + + def test_moving_to_the_group_they_are_in_is_disabled(self, pool): + actions = self._actions(pool.group_menu([0, 1])) + assert not actions["Move to ChR2"].isEnabled() + assert actions["Move to eYFP"].isEnabled() + + def test_choosing_a_group_reports_the_sessions_and_the_group(self, pool): + seen = [] + pool.group_assigned.connect(lambda i, g: seen.append((i, g))) + self._actions(pool.group_menu([0, 2]))["Move to eYFP"].trigger() + self._actions(pool.group_menu([2]))["Remove from group"].trigger() + assert seen == [([0, 2], "eYFP"), ([2], "")] + + def test_a_new_group_asks_for_its_name(self, pool, monkeypatch): + seen = [] + pool.group_assigned.connect(lambda i, g: seen.append((i, g))) + monkeypatch.setattr(pool, "_ask_group_name", lambda: "Saline") + self._actions(pool.group_menu([1]))["New group…"].trigger() + monkeypatch.setattr(pool, "_ask_group_name", lambda: "") # cancelled + self._actions(pool.group_menu([1]))["New group…"].trigger() + assert seen == [([1], "Saline")] + + def test_an_ungrouped_cohort_can_remove_nothing(self, qtbot): + pool = SessionPool() + qtbot.addWidget(pool) + pool.set_entries([PoolEntry("a", "", 1.0, False, False, False)]) + actions = self._actions(pool.group_menu([0])) + assert not actions["Remove from group"].isEnabled() + assert not any(t.startswith("Move to") for t in actions) From 735b4df7e5e46ddb687ee65c78b7f04388bed5d0 Mon Sep 17 00:00:00 2001 From: gbradham Date: Thu, 29 Oct 2026 13:21:12 -0500 Subject: [PATCH 2/4] feat(review): export time bins and show keypoint distances Export > Time bins (CSV) splits each session, on its own zero, into bins of 1, 3, 5, 10 or a custom number of minutes, tidy or wide, over the whole session or the selected range. The viewer's Distances chip draws the pixel distance between keypoint pairs; right-click it to choose the pairs, which are remembered per keypoint set. --- src/glider/analysis/epochs.py | 29 +++ src/glider/gui/behavior/analysis_window.py | 166 ++++++++++++++++++ src/glider/gui/review/time_bins_dialog.py | 80 +++++++++ src/glider/gui/review/viewer.py | 141 ++++++++++++++- tests/unit/analysis/test_epochs.py | 25 +++ .../unit/gui/behavior/test_analysis_window.py | 106 ++++++++++- tests/unit/gui/review/test_viewer.py | 64 +++++++ 7 files changed, 609 insertions(+), 2 deletions(-) create mode 100644 src/glider/gui/review/time_bins_dialog.py diff --git a/src/glider/analysis/epochs.py b/src/glider/analysis/epochs.py index 571329f0..284a9983 100644 --- a/src/glider/analysis/epochs.py +++ b/src/glider/analysis/epochs.py @@ -29,6 +29,7 @@ "summarize", "threshold_text", "tidy_frame", + "time_bins", "wide_frame", ] @@ -232,6 +233,34 @@ def tidy_frame( return pd.DataFrame.from_records(records, columns=columns) +def _bin_edge(seconds: float, unit: float) -> str: + """``90`` s as ``1.5`` in minutes: whole numbers without a trailing ``.0``.""" + return f"{seconds / unit:g}" + + +def time_bins(start_s: float, end_s: float, bin_s: float) -> list[Epoch]: + """``[start_s, end_s)`` cut into consecutive bins of ``bin_s``, as epochs. + + Bin edges sit on whole multiples of ``bin_s`` from zero (flow start, or the + first video frame), so minute 3 is the same stretch in every export + whatever span was asked for. The last bin is clipped to ``end_s`` and named + for what it covers. Labelled in minutes when the bin is a whole number of + them, in seconds otherwise. + """ + if not bin_s > 0 or not end_s > start_s: + return [] + unit, suffix = (60.0, "min") if float(bin_s) % 60 == 0 else (1.0, "s") + first = math.floor(start_s / bin_s + 1e-9) + bins = [] + k = first + while k * bin_s < end_s - 1e-9: + lo, hi = max(k * bin_s, start_s), min((k + 1) * bin_s, end_s) + name = f"{_bin_edge(lo, unit)}–{_bin_edge(hi, unit)} {suffix}" + bins.append(Epoch(f"bin_{k - first + 1}", name, lo, hi)) + k += 1 + return bins + + def wide_frame(tidy: pd.DataFrame, metrics: list[str]) -> pd.DataFrame: """One row per session: ``session, group, __…``. diff --git a/src/glider/gui/behavior/analysis_window.py b/src/glider/gui/behavior/analysis_window.py index b04285cc..2d2848a4 100644 --- a/src/glider/gui/behavior/analysis_window.py +++ b/src/glider/gui/behavior/analysis_window.py @@ -15,6 +15,7 @@ import hashlib import logging +import math from collections import Counter from dataclasses import dataclass, replace from pathlib import Path @@ -71,6 +72,7 @@ metric_catalog, threshold_text, tidy_frame, + time_bins, wide_frame, ) from glider.analysis.timeline import ( @@ -586,6 +588,7 @@ def _build_top_bar(self) -> QFrame: ("Range stats (CSV)…", self._export_window), ("Epoch table, tidy (CSV)…", lambda: self._export_epochs("tidy")), ("Epoch table, wide (CSV)…", lambda: self._export_epochs("wide")), + ("Time bins (CSV)…", self._export_bins), ("Markers (CSV)…", self._export_markers), ("Retrieval tables (CSV)…", self._export_retrieval), ("Nest distance per frame (CSV)…", self._export_nest_distance), @@ -636,6 +639,10 @@ def _wire_transport(self) -> None: ) t.video_on.toggled.connect(self._canvas.set_show_video) t.poses_on.toggled.connect(self._canvas.set_show_poses) + t.distances_on.toggled.connect(self._canvas.set_show_distances) + t.distances_on.customContextMenuRequested.connect( + lambda pos: self._pairs_menu().exec(t.distances_on.mapToGlobal(pos)) + ) t.pups_on.setVisible(False) # shown for a pup session (_publish_retrieval) t.pups_on.toggled.connect(self._canvas.set_show_pups) t.trail_on.toggled.connect(self._apply_trail) @@ -1133,6 +1140,7 @@ def _adopt(self, path: Path, view: SessionView, timeline) -> None: self._canvas.set_hud_names(self._bar.behavior_order()) self._bar.set_hidden(_settings().value(self._hidden_key(), [], type=list) or []) self._canvas.set_view(view) + self._restore_pairs() # The overlay belongs to the session that just left. self._canvas.set_heatmap(None) self._refresh_heatmap_export_state() @@ -1661,6 +1669,59 @@ def _hidden_key(self) -> str: path = str(Path(self._cohort[self._shown][0]).resolve()) return f"review/hidden/{hashlib.sha1(path.encode('utf-8')).hexdigest()[:16]}" + def _pairs_key(self) -> str: + """Per keypoint set: the same model's sessions share their chosen pairs.""" + names = "|".join(self._view.keypoint_names) if self._view is not None else "" + return f"review/kp_pairs/{hashlib.sha1(names.encode('utf-8')).hexdigest()[:16]}" + + def _restore_pairs(self) -> None: + saved = _settings().value(self._pairs_key(), None, type=list) + if not saved: + self._canvas.set_pairs(None) + return + try: + pairs = [tuple(int(v) for v in str(item).split("-")) for item in saved] + except ValueError: + pairs = None + self._canvas.set_pairs(pairs) + + def _set_pairs(self, pairs) -> None: + self._canvas.set_pairs(pairs) + if pairs is None: + _settings().remove(self._pairs_key()) + else: + _settings().setValue(self._pairs_key(), [f"{i}-{j}" for i, j in pairs]) + # Choosing pairs is asking to see them. + self._transport.distances_on.setChecked(True) + + def _pairs_menu(self) -> QMenu: + """Which keypoint pairs the Distances overlay measures.""" + from glider.gui.review.viewer import all_pairs + + menu = QMenu(self) + names = list(self._view.keypoint_names) if self._view is not None else [] + if len(names) < 2: + menu.addAction("Needs two or more keypoints").setEnabled(False) + return menu + current = set(self._canvas.pairs()) + menu.addAction("Default pairs", lambda: self._set_pairs(None)) + menu.addAction("All pairs", lambda: self._set_pairs(all_pairs(names))) + menu.addAction("None", lambda: self._set_pairs([])) + menu.addSeparator() + for pair in all_pairs(names): + i, j = pair + action = menu.addAction(f"{names[i]} – {names[j]}") + action.setCheckable(True) + action.setChecked(pair in current) + action.toggled.connect( + lambda on, p=pair: self._set_pairs( + sorted(set(self._canvas.pairs()) | {p}) + if on + else [q for q in self._canvas.pairs() if q != p] + ) + ) + return menu + def _remember_hidden(self, keys: list) -> None: # Only on a change: _adopt's own set_hidden echoes the saved list back. key = self._hidden_key() @@ -2973,6 +3034,111 @@ def _export_epochs(self, shape: str) -> None: return self.statusBar().showMessage(f"Wrote {path}", 10000) + def _cohort_span(self) -> tuple[float, float] | None: + """Seconds from the earliest scored frame to the end of the latest, cohort-wide. + + Each session on its own zero (see :meth:`range_rows`), so the span is + what a range would have to cover to reach every frame of every animal. + """ + spans = [] + for index, (_path, view) in enumerate(self._cohort): + bounds = self._scored_bounds(view) + if bounds is None: + continue + timeline = self._timeline_at(index) + spans.append( + ( + mk.seconds_at(timeline, view.fps, bounds[0]), + mk.seconds_at(timeline, view.fps, bounds[1] + 1), + ) + ) + if not spans: + return None + return min(s for s, _ in spans), max(e for _, e in spans) + + def _run_bins_dialog(self, minutes: float, *, has_selection: bool): + """``(bin_s, whole_session, shape)`` as chosen, or None if cancelled.""" + from glider.gui.review.time_bins_dialog import TimeBinsDialog + + dialog = TimeBinsDialog(minutes, has_selection=has_selection, parent=self) + try: + if dialog.exec() != QDialog.DialogCode.Accepted: + return None + return dialog.bin_seconds(), dialog.whole_session(), dialog.layout_shape() + finally: + dialog.deleteLater() + + def bin_frame(self, start_s: float, end_s: float, bin_s: float, shape: str): + """Every metric over consecutive ``bin_s`` bins of ``[start_s, end_s)``. + + The epoch table's arithmetic with bins as the epochs: tidy is a row per + session per bin, wide a row per session with ``__`` + columns. Every metric the cohort can answer goes in, as the range + stats export writes every column rather than the table's picks. + """ + bins = time_bins(start_s, end_s, bin_s) + if not bins: + return None + rows = _zero_fill_across_epochs({b.key: self.range_rows(b.start_s, b.end_s) for b in bins}) + metrics = [m.key for m in metric_catalog(r for block in rows.values() for r in block)] + tidy = tidy_frame(bins, rows, metrics) + if shape == "wide": + return wide_frame(tidy, metrics) + edges = {b.name: (b.start_s, b.end_s) for b in bins} + tidy = tidy.rename(columns={"range": "bin"}) + tidy.insert(3, "bin_start_s", [edges[name][0] for name in tidy["bin"]]) + tidy.insert(4, "bin_end_s", [edges[name][1] for name in tidy["bin"]]) + return tidy + + def _export_bins(self) -> None: + """The cohort cut into fixed-width time bins: every metric, per bin.""" + if not self._cohort: + return + settings = _settings() + try: + minutes = float(settings.value("review/bin_minutes", 1.0)) + except (TypeError, ValueError): + minutes = 1.0 + picked = self._run_bins_dialog(minutes, has_selection=self._current_span is not None) + if picked is None: + return + bin_s, whole, shape = picked + settings.setValue("review/bin_minutes", bin_s / 60.0) + span = self._cohort_span() if whole else self._current_span + if span is None: + QMessageBox.information( + self, "Export time bins", "Nothing to export: no session has scored frames." + ) + return + start_s, end_s = span + if whole: + # From zero (flow start), or earlier for frames before it, on a + # bin edge so every bin is whole but the last. + start_s = math.floor(min(0.0, start_s) / bin_s) * bin_s + root = self._cohort_root or self._session_folder(max(0, self._shown)) + path, _ = QFileDialog.getSaveFileName( + self, + "Export time bins", + str(root / f"time_bins_{bin_s / 60.0:g}min_{shape}.csv"), + "CSV Files (*.csv)", + ) + if not path: + return + QApplication.setOverrideCursor(Qt.CursorShape.WaitCursor) + try: + frame = self.bin_frame(start_s, end_s, bin_s, shape) + finally: + QApplication.restoreOverrideCursor() + if frame is None: + QMessageBox.information(self, "Export time bins", "The span holds no whole frame.") + return + try: + frame.to_csv(path, index=False) + except OSError as e: + QMessageBox.critical(self, "Export time bins", f"Could not write {path}: {e}") + return + self.statusBar().showMessage(f"Wrote {path}", 10000) + def _export_heatmap(self) -> None: """Write the heatmap on screen as a PNG figure plus a CSV of its grid.""" selection = self._bar.selection() diff --git a/src/glider/gui/review/time_bins_dialog.py b/src/glider/gui/review/time_bins_dialog.py new file mode 100644 index 00000000..9023161c --- /dev/null +++ b/src/glider/gui/review/time_bins_dialog.py @@ -0,0 +1,80 @@ +"""Choose how to cut the session into time bins for the binned export.""" + +from __future__ import annotations + +from PyQt6.QtWidgets import ( + QComboBox, + QDialog, + QDialogButtonBox, + QDoubleSpinBox, + QFormLayout, +) + +__all__ = ["BIN_PRESETS_MIN", "TimeBinsDialog"] + +#: Offered in the bin box's drop-down; any other width can be typed. +BIN_PRESETS_MIN = (1.0, 3.0, 5.0, 10.0) + +WHOLE, SELECTION = "whole", "selection" + + +class TimeBinsDialog(QDialog): + """Bin width, the span to bin, and the CSV's shape.""" + + def __init__(self, minutes: float, *, has_selection: bool, parent=None): + super().__init__(parent) + self.setWindowTitle("Export time bins") + form = QFormLayout(self) + self.preset = QComboBox() + for value in BIN_PRESETS_MIN: + self.preset.addItem(f"{value:g} min", value) + self.preset.addItem("Custom", None) + self.bin_width = QDoubleSpinBox() + self.bin_width.setRange(0.05, 24 * 60.0) + self.bin_width.setDecimals(2) + self.bin_width.setSingleStep(1.0) + self.bin_width.setSuffix(" min") + self.bin_width.setValue(minutes) + self.span = QComboBox() + self.span.addItem("Whole session", WHOLE) + self.span.addItem("Selected range", SELECTION) + if has_selection: + self.span.setCurrentIndex(1) + else: + # Disabled rather than hidden: the option is the hint that a + # selection narrows the export. + self.span.model().item(1).setEnabled(False) + self.shape = QComboBox() + self.shape.addItem("Tidy: a row per session per bin", "tidy") + self.shape.addItem("Wide: a row per session", "wide") + self._sync_preset(minutes) + self.preset.activated.connect(self._preset_picked) + self.bin_width.valueChanged.connect(self._sync_preset) + form.addRow("Bin width", self.preset) + form.addRow("", self.bin_width) + form.addRow("Span", self.span) + form.addRow("Layout", self.shape) + buttons = QDialogButtonBox( + QDialogButtonBox.StandardButton.Ok | QDialogButtonBox.StandardButton.Cancel + ) + buttons.accepted.connect(self.accept) + buttons.rejected.connect(self.reject) + form.addRow(buttons) + + def _preset_picked(self, index: int) -> None: + value = self.preset.itemData(index) + if value is not None: + self.bin_width.setValue(value) + + def _sync_preset(self, minutes: float) -> None: + index = self.preset.findData(round(float(minutes), 2)) + self.preset.setCurrentIndex(index if index >= 0 else self.preset.count() - 1) + + def bin_seconds(self) -> float: + return self.bin_width.value() * 60.0 + + def whole_session(self) -> bool: + return self.span.currentData() == WHOLE + + def layout_shape(self) -> str: + return self.shape.currentData() diff --git a/src/glider/gui/review/viewer.py b/src/glider/gui/review/viewer.py index ae6006a0..8f70eec9 100644 --- a/src/glider/gui/review/viewer.py +++ b/src/glider/gui/review/viewer.py @@ -39,12 +39,64 @@ from glider.gui.widgets.pastel_glyphs import lucide_icon from glider.gui.widgets.tool_ui import ElidedLabel, data_font, set_button_role, set_text_role -__all__ = ["TRAIL_DEFAULT_S", "KeypointCanvas", "Transport"] +__all__ = [ + "TRAIL_DEFAULT_S", + "KeypointCanvas", + "Transport", + "all_pairs", + "default_pairs", + "pair_distances", +] logger = logging.getLogger(__name__) TRAIL_DEFAULT_S = 5.0 +# The mouse skeleton the lab's pose models share: head and flanks to the body +# centre, ears to the nose. Any spelling of the centre the pipeline accepts. +_CENTRES = ("center", "center_body", "body_center", "centroid") +_MOUSE_EDGES = ( + ("nose", "{c}"), + ("l_ear", "nose"), + ("r_ear", "nose"), + ("l_body", "{c}"), + ("r_body", "{c}"), + ("tail_base", "{c}"), +) + + +def all_pairs(names: list[str]) -> list[tuple[int, int]]: + """Every unordered pair of keypoints, in list order.""" + return [(i, j) for i in range(len(names)) for j in range(i + 1, len(names))] + + +def default_pairs(names: list[str]) -> list[tuple[int, int]]: + """The pairs measured until the reviewer picks others. + + The mouse skeleton's edges where the keypoints are the lab's mouse set; + otherwise each keypoint to the next, which for a model's own order is + usually a chain along the body. + """ + index = {name.lower(): i for i, name in enumerate(names)} + centre = next((c for c in _CENTRES if c in index), None) + if centre is not None: + edges = [(a, b.format(c=centre)) for a, b in _MOUSE_EDGES] + found = [(index[a], index[b]) for a, b in edges if a in index and b in index] + if len(found) >= 2: + return found + return [(i, i + 1) for i in range(len(names) - 1)] + + +def pair_distances(points: np.ndarray, pairs) -> list[float]: + """Euclidean distance in pixels for each ``(i, j)``; NaN where either is missing.""" + out = [] + for i, j in pairs: + if max(i, j) >= len(points): + out.append(float("nan")) + continue + out.append(float(np.hypot(*(points[i] - points[j])))) + return out + class KeypointCanvas(QWidget): """The animal drawn from its poses, with a trailing centroid track. @@ -70,6 +122,10 @@ def __init__(self, parent=None): self._video_failed = False # this view's video would not open; retried on set_view self._cached: tuple[int, QImage] | None = None self._show_poses = True + self._show_distances = False + # Keypoint index pairs to measure, per keypoint set: a cohort can mix + # pose models, and indices only mean something within one. + self._pairs: dict[tuple[str, ...], list[tuple[int, int]]] = {} self._show_zones = True self._show_hud = True self._hud: tuple[tuple[str, QColor] | None, list[tuple[str, QColor]]] | None = None @@ -282,6 +338,10 @@ def to_widget(point): ) painter.setPen(QPen(QColor(colors.CANVAS), 1)) painter.drawEllipse(to_widget(point), 5, 5) + if self._show_distances: + painter.save() + self._paint_distances(painter, to_widget, points) + painter.restore() if self._show_pups and self._pups is not None: painter.save() @@ -437,6 +497,77 @@ def set_show_poses(self, enabled: bool) -> None: self._show_poses = bool(enabled) self.update() + def set_show_distances(self, enabled: bool) -> None: + """Label each measured pair with its length in source-video pixels.""" + self._show_distances = bool(enabled) + self.update() + + def _names(self) -> tuple[str, ...]: + return tuple(self._view.keypoint_names) if self._view is not None else () + + def set_pairs(self, pairs) -> None: + """The keypoint pairs to measure for this session's keypoint set, as + index pairs; None goes back to :func:`default_pairs`.""" + names = self._names() + if pairs is None: + self._pairs.pop(names, None) + else: + self._pairs[names] = [ + (int(i), int(j)) for i, j in pairs if max(i, j) < len(names) and i != j + ] + self.update() + + def pairs(self) -> list[tuple[int, int]]: + names = self._names() + if names in self._pairs: + return list(self._pairs[names]) + return default_pairs(list(names)) + + def _paint_distances(self, painter, to_widget, points: np.ndarray) -> None: + """A dashed line per pair, its length in px beside the line's middle. + + Measured in the video's own pixels, not the widget's: the number must + not change with the window size. Each label sits off its line, on the + side away from the body unless an earlier label is already there: + mirrored pairs (both ears to the nose) would otherwise stack. + """ + font = QFont(self.font()) + size = font.pointSizeF() if font.pointSizeF() > 0 else 10.0 + font.setPointSizeF(size * 0.9) + painter.setFont(font) + metrics = QFontMetricsF(font) + line = QPen(colors.qcolor_with_alpha(QColor(colors.TEXT_PRIMARY), 0.7), 1) + line.setStyle(Qt.PenStyle.DashLine) + finite = points[np.isfinite(points).all(axis=1)] + body = to_widget(finite.mean(axis=0)) if len(finite) else QPointF() + placed: list[QRectF] = [] + pairs = self.pairs() + for (i, j), px in zip(pairs, pair_distances(points, pairs), strict=True): + if not np.isfinite(px): + continue + a, b = to_widget(points[i]), to_widget(points[j]) + painter.setPen(line) + painter.drawLine(a, b) + text = f"{px:.0f} px" + mid = (a + b) / 2 + dx, dy = b.x() - a.x(), b.y() - a.y() + length = max((dx * dx + dy * dy) ** 0.5, 1e-9) + normal = QPointF(-dy / length, dx / length) + away = mid - body + if normal.x() * away.x() + normal.y() * away.y() < 0: + normal = -normal + half_w, half_h = metrics.horizontalAdvance(text) / 2, metrics.height() / 2 + # Far enough along the normal that the text box clears the line. + reach = abs(normal.x()) * half_w + abs(normal.y()) * half_h + 3 + boxes = [ + QRectF(c.x() - half_w, c.y() - half_h, 2 * half_w, 2 * half_h) + for c in (mid + normal * reach, mid - normal * reach) + ] + box = next((r for r in boxes if not any(r.intersects(q) for q in placed)), boxes[0]) + placed.append(box) + at = QPointF(box.left(), box.top() + metrics.ascent()) + self._paint_label(painter, at, text, QColor(colors.TEXT_PRIMARY)) + def set_retrieval(self, pups, result, nest_at) -> None: """Draw a pup session: its ``PupTracks``, its ``SessionRetrieval`` (or None, unscored) and ``nest_at(frame)``, the nest ``Zone`` in force on a frame. @@ -607,6 +738,13 @@ def __init__(self, parent=None): row.addWidget(self.bout, 1) self.video_on = _toggle("Video", "Draw the session's video", checked=True) self.poses_on = _toggle("Poses", "Draw the tracked keypoints", checked=True) + self.distances_on = _toggle( + "Distances", + "Label keypoint pairs with their distance in video pixels. " + "Right-click to choose the pairs.", + checked=False, + ) + self.distances_on.setContextMenuPolicy(Qt.ContextMenuPolicy.CustomContextMenu) self.pups_on = _toggle( "Pups", "Draw the nest, each pup and the dam's nose line", checked=True ) @@ -621,6 +759,7 @@ def __init__(self, parent=None): for toggle in ( self.video_on, self.poses_on, + self.distances_on, self.pups_on, self.trail_on, self.heatmap_on, diff --git a/tests/unit/analysis/test_epochs.py b/tests/unit/analysis/test_epochs.py index 0ec88c64..8c6bc909 100644 --- a/tests/unit/analysis/test_epochs.py +++ b/tests/unit/analysis/test_epochs.py @@ -18,6 +18,7 @@ summarize, threshold_text, tidy_frame, + time_bins, wide_frame, ) @@ -176,3 +177,27 @@ def test_retrieval_metrics_join_the_catalog_only_when_present(): keys = [m.key for m in metric_catalog([row])] assert {"pups_retrieved", "dam_in_nest_frac", "nose_nest_mean_cm"} <= set(keys) assert "pups_retrieved" not in [m.key for m in metric_catalog([{}])] + + +class TestTimeBins: + def test_whole_minutes_are_named_in_minutes(self): + bins = time_bins(0.0, 180.0, 60.0) + assert [b.name for b in bins] == ["0–1 min", "1–2 min", "2–3 min"] + assert [(b.start_s, b.end_s) for b in bins] == [(0, 60), (60, 120), (120, 180)] + assert len({b.key for b in bins}) == 3 + + def test_the_last_bin_is_clipped_to_the_end(self): + bins = time_bins(0.0, 150.0, 60.0) + assert bins[-1].name == "2–2.5 min" + assert bins[-1].end_s == 150.0 + + def test_other_widths_are_named_in_seconds(self): + assert [b.name for b in time_bins(0.0, 60.0, 30.0)] == ["0–30 s", "30–60 s"] + + def test_edges_stay_on_multiples_of_the_width(self): + bins = time_bins(90.0, 200.0, 60.0) + assert [(b.start_s, b.end_s) for b in bins] == [(90, 120), (120, 180), (180, 200)] + + @pytest.mark.parametrize("start, end, width", [(0, 10, 0), (10, 10, 5), (10, 0, 5)]) + def test_an_empty_span_or_width_gives_no_bins(self, start, end, width): + assert time_bins(start, end, width) == [] diff --git a/tests/unit/gui/behavior/test_analysis_window.py b/tests/unit/gui/behavior/test_analysis_window.py index 4ef27a6f..f70d3f4b 100644 --- a/tests/unit/gui/behavior/test_analysis_window.py +++ b/tests/unit/gui/behavior/test_analysis_window.py @@ -11,7 +11,7 @@ pytest.importorskip("PyQt6") -from PyQt6.QtCore import Qt # noqa: E402 +from PyQt6.QtCore import QSettings, Qt # noqa: E402 from PyQt6.QtGui import QColor, QKeyEvent # noqa: E402 from glider.analysis.behavior.session_view import SessionView # noqa: E402 @@ -2167,6 +2167,109 @@ def test_the_status_bar_counts_markers(self, qtbot, tmp_path): assert "1 marker" in win._status.text() +class TestTimeBinExport: + """Fixed-width bins of every session: the epoch table's maths, bins as epochs.""" + + def _win(self, qtbot, tmp_path): + # 300 frames at 30 fps: 5 s of groom, then 5 s of rear. + paths = [ + _ethogram(tmp_path / sid / "analysis", ["groom"] * 150 + ["rear"] * 150) + for sid in ("s1", "s2") + ] + win = AnalysisWindow() + qtbot.addWidget(win) + win.load_many(paths) + return win + + def test_tidy_is_a_row_per_session_per_bin(self, qtbot, tmp_path): + win = self._win(qtbot, tmp_path) + frame = win.bin_frame(0.0, 10.0, 5.0, "tidy") + assert len(frame) == 4 + assert list(frame.columns[:5]) == ["session", "group", "bin", "bin_start_s", "bin_end_s"] + first = frame[frame["bin"] == "0–5 s"] + second = frame[frame["bin"] == "5–10 s"] + assert first["groom_s"].tolist() == pytest.approx([5.0, 5.0]) + assert first["rear_s"].tolist() == pytest.approx([0.0, 0.0]) + assert second["rear_s"].tolist() == pytest.approx([5.0, 5.0]) + + def test_wide_is_a_row_per_session(self, qtbot, tmp_path): + win = self._win(qtbot, tmp_path) + frame = win.bin_frame(0.0, 10.0, 5.0, "wide") + assert len(frame) == 2 + assert "5–10 s__rear_s" in frame.columns + + def test_the_whole_session_span_reaches_the_last_frame(self, qtbot, tmp_path): + win = self._win(qtbot, tmp_path) + assert win._cohort_span() == pytest.approx((0.0, 10.0)) + + def test_export_writes_the_bins(self, qtbot, tmp_path, monkeypatch): + win = self._win(qtbot, tmp_path) + out = tmp_path / "bins.csv" + monkeypatch.setattr(win, "_run_bins_dialog", lambda *a, **k: (4.0, True, "tidy")) + monkeypatch.setattr( + "glider.gui.behavior.analysis_window._settings", + lambda: QSettings(str(tmp_path / "s.ini"), QSettings.Format.IniFormat), + ) + monkeypatch.setattr( + "glider.gui.behavior.analysis_window.QFileDialog.getSaveFileName", + lambda *a, **k: (str(out), ""), + ) + win._export_bins() + written = pd.read_csv(out) + # 10 s in 4 s bins: two whole and one clipped, for each of two sessions. + assert sorted(set(written["bin"])) == ["0–4 s", "4–8 s", "8–10 s"] + assert len(written) == 6 + + def test_cancelling_the_dialog_writes_nothing(self, qtbot, tmp_path, monkeypatch): + win = self._win(qtbot, tmp_path) + monkeypatch.setattr(win, "_run_bins_dialog", lambda *a, **k: None) + monkeypatch.setattr( + "glider.gui.behavior.analysis_window._settings", + lambda: QSettings(str(tmp_path / "s.ini"), QSettings.Format.IniFormat), + ) + asked = [] + monkeypatch.setattr( + "glider.gui.behavior.analysis_window.QFileDialog.getSaveFileName", + lambda *a, **k: asked.append(a) or ("", ""), + ) + win._export_bins() + assert asked == [] + + +class TestKeypointDistances: + def _win(self, qtbot, tmp_path, monkeypatch): + monkeypatch.setattr( + "glider.gui.behavior.analysis_window._settings", + lambda: QSettings(str(tmp_path / "s.ini"), QSettings.Format.IniFormat), + ) + win = AnalysisWindow() + qtbot.addWidget(win) + win.load(_session(tmp_path / "v")) + return win + + def test_the_toggle_drives_the_canvas(self, qtbot, tmp_path, monkeypatch): + win = self._win(qtbot, tmp_path, monkeypatch) + assert win._canvas._show_distances is False + win._transport.distances_on.setChecked(True) + assert win._canvas._show_distances is True + + def test_chosen_pairs_survive_reopening(self, qtbot, tmp_path, monkeypatch): + win = self._win(qtbot, tmp_path, monkeypatch) + assert len(win._view.keypoint_names) >= 2 + win._set_pairs([(0, 1)]) + assert win._transport.distances_on.isChecked() + again = AnalysisWindow() + qtbot.addWidget(again) + again.load(_session(tmp_path / "w")) + assert again._canvas.pairs() == [(0, 1)] + + def test_the_menu_lists_every_pair(self, qtbot, tmp_path, monkeypatch): + win = self._win(qtbot, tmp_path, monkeypatch) + n = len(win._view.keypoint_names) + checkable = [a for a in win._pairs_menu().actions() if a.isCheckable()] + assert len(checkable) == n * (n - 1) // 2 + + class TestTheEpochTableInTheWindow: """Cohort range markers become columns; a row is an animal; exports match.""" @@ -2262,6 +2365,7 @@ def test_the_export_menu(self, qtbot, tmp_path): "Range stats (CSV)…", "Epoch table, tidy (CSV)…", "Epoch table, wide (CSV)…", + "Time bins (CSV)…", "Markers (CSV)…", "Retrieval tables (CSV)…", "Nest distance per frame (CSV)…", diff --git a/tests/unit/gui/review/test_viewer.py b/tests/unit/gui/review/test_viewer.py index ed64be65..fd67555c 100644 --- a/tests/unit/gui/review/test_viewer.py +++ b/tests/unit/gui/review/test_viewer.py @@ -17,6 +17,9 @@ CHIP_H, KeypointCanvas, Transport, + all_pairs, + default_pairs, + pair_distances, ) from glider.gui.styles import colors # noqa: E402 @@ -252,3 +255,64 @@ def test_a_video_that_will_not_open_is_not_forgotten_by_the_session(canvas, monk assert not canvas.has_video() canvas.set_view(view) # shown again: tried again assert canvas.has_video() + + +MOUSE = ["nose", "l_ear", "r_ear", "l_body", "center", "r_body", "tail_base"] + + +class TestKeypointPairs: + def test_the_mouse_set_gets_its_skeleton(self): + pairs = default_pairs(MOUSE) + assert (0, 4) in pairs and (1, 0) in pairs and (6, 4) in pairs + assert len(pairs) == 6 + + def test_any_spelling_of_the_centre(self): + assert default_pairs(["nose", "body_center", "tail_base"]) == [(0, 1), (2, 1)] + + def test_an_unknown_set_is_chained_in_order(self): + assert default_pairs(["a", "b", "c"]) == [(0, 1), (1, 2)] + assert default_pairs(["a"]) == [] + + def test_all_pairs(self): + assert all_pairs(["a", "b", "c"]) == [(0, 1), (0, 2), (1, 2)] + + def test_distances_are_in_pixels_and_nan_when_missing(self): + points = np.array([[0.0, 0.0], [3.0, 4.0], [np.nan, np.nan]]) + got = pair_distances(points, [(0, 1), (0, 2), (0, 9)]) + assert got[0] == pytest.approx(5.0) + assert np.isnan(got[1]) and np.isnan(got[2]) + + +class TestDistanceOverlay: + def _mouse_view(self): + xy = np.zeros((10, 7, 2)) + xy[:, :, 0] = np.arange(7) * 20 + 100 + xy[:, :, 1] = 200.0 + return _view(labels=["groom"] * 10, frames=np.arange(10), xy=xy, keypoint_names=MOUSE) + + def test_pairs_default_until_chosen(self, canvas): + canvas.set_view(self._mouse_view()) + assert canvas.pairs() == default_pairs(MOUSE) + canvas.set_pairs([(0, 6), (0, 99)]) # an index past the set is dropped + assert canvas.pairs() == [(0, 6)] + canvas.set_pairs(None) + assert canvas.pairs() == default_pairs(MOUSE) + + def test_pairs_belong_to_their_keypoint_set(self, canvas): + canvas.set_view(self._mouse_view()) + canvas.set_pairs([(0, 6)]) + canvas.set_view(_view()) # one keypoint: nothing to measure + assert canvas.pairs() == [] + canvas.set_view(self._mouse_view()) + assert canvas.pairs() == [(0, 6)] + + def test_drawing_distances_paints_without_error(self, canvas): + canvas.set_view(self._mouse_view()) + canvas.set_show_distances(True) + canvas.grab() # runs paintEvent + + def test_the_transport_has_the_toggle_off(self, qtbot): + transport = Transport() + qtbot.addWidget(transport) + assert transport.distances_on.text() == "Distances" + assert not transport.distances_on.isChecked() From 1b75b1bd4d2fd9b272625d36d83855e506ea8755 Mon Sep 17 00:00:00 2001 From: gbradham Date: Thu, 29 Oct 2026 13:21:12 -0500 Subject: [PATCH 3/4] feat(stats): in-house statistics engine for Session Review Qt-free, scipy-only tests on the epoch and time-bin tables: group vs group, epoch vs epoch, group x epoch (mixed ANOVA with sphericity corrections) and trends across bins (per-animal slopes, polynomial contrasts, first vs last). Auto test selection explains its choice; effect sizes carry CIs. Reference values pinned against pingouin and scikit-posthocs. --- src/glider/analysis/stats/__init__.py | 89 +++ src/glider/analysis/stats/compare.py | 704 +++++++++++++++++++++ src/glider/analysis/stats/corrections.py | 46 ++ src/glider/analysis/stats/frames.py | 98 +++ src/glider/analysis/stats/repeated.py | 483 ++++++++++++++ src/glider/analysis/stats/result.py | 87 +++ src/glider/analysis/stats/trend.py | 358 +++++++++++ tests/unit/analysis/test_stats.py | 353 +++++++++++ tests/unit/analysis/test_stats_repeated.py | 407 ++++++++++++ tests/unit/analysis/test_stats_trend.py | 227 +++++++ 10 files changed, 2852 insertions(+) create mode 100644 src/glider/analysis/stats/__init__.py create mode 100644 src/glider/analysis/stats/compare.py create mode 100644 src/glider/analysis/stats/corrections.py create mode 100644 src/glider/analysis/stats/frames.py create mode 100644 src/glider/analysis/stats/repeated.py create mode 100644 src/glider/analysis/stats/result.py create mode 100644 src/glider/analysis/stats/trend.py create mode 100644 tests/unit/analysis/test_stats.py create mode 100644 tests/unit/analysis/test_stats_repeated.py create mode 100644 tests/unit/analysis/test_stats_trend.py diff --git a/src/glider/analysis/stats/__init__.py b/src/glider/analysis/stats/__init__.py new file mode 100644 index 00000000..812af1cf --- /dev/null +++ b/src/glider/analysis/stats/__init__.py @@ -0,0 +1,89 @@ +"""Statistics on the cohort's own numbers: the tests a paper needs. Qt-free. + +The input is the epoch table's tidy frame (:func:`glider.analysis.epochs.tidy_frame`): +one row per animal per epoch, ``session, group, range, …``. The animal is +always the unit -- one value per animal per cell, never frames. + +Every test returns a :class:`TestResult` that carries what a reader needs to +trust it: the statistic and its degrees of freedom, the exact p, an effect size +with a 95% CI, n per cell, the assumption checks, which animals were left out +and why, and -- when the test was picked automatically -- the reason it was. + +* :func:`compare_groups` -- between subjects, one epoch: two groups or several. +* :func:`compare_epochs` -- within subjects: the same animals in two or more epochs. +* :func:`compare_mixed` -- group x epoch: the mixed (split-plot) ANOVA. +* :func:`compare_trend` -- change across time bins: slopes, contrasts, first vs last. +* :func:`run_between` / :func:`run_paired` -- the same on plain arrays. +* :mod:`.corrections` -- Bonferroni, Holm and Benjamini-Hochberg. +""" + +from glider.analysis.stats.compare import ( + BETWEEN_TESTS, + ONE_SAMPLE_TESTS, + PAIRED_TESTS, + compare_epochs, + compare_groups, + dunn, + games_howell, + run_between, + run_one_sample, + run_paired, + tukey, +) +from glider.analysis.stats.corrections import adjust +from glider.analysis.stats.frames import UNGROUPED, Dropped, between_cells, within_matrix +from glider.analysis.stats.repeated import ( + REPEATED_TESTS, + MixedResult, + Sphericity, + compare_mixed, + pairwise_within, + run_mixed, + run_repeated, + sphericity, +) +from glider.analysis.stats.result import TestResult +from glider.analysis.stats.trend import ( + GroupTrend, + TrendResult, + animal_slopes, + compare_trend, + polynomial_contrasts, + run_trend, + trend_contrasts, +) + +__all__ = [ + "BETWEEN_TESTS", + "ONE_SAMPLE_TESTS", + "PAIRED_TESTS", + "REPEATED_TESTS", + "MixedResult", + "Sphericity", + "UNGROUPED", + "Dropped", + "GroupTrend", + "TestResult", + "TrendResult", + "adjust", + "animal_slopes", + "between_cells", + "compare_epochs", + "compare_groups", + "compare_mixed", + "compare_trend", + "dunn", + "games_howell", + "pairwise_within", + "polynomial_contrasts", + "run_between", + "run_mixed", + "run_one_sample", + "run_paired", + "run_repeated", + "run_trend", + "sphericity", + "trend_contrasts", + "tukey", + "within_matrix", +] diff --git a/src/glider/analysis/stats/compare.py b/src/glider/analysis/stats/compare.py new file mode 100644 index 00000000..75d404c9 --- /dev/null +++ b/src/glider/analysis/stats/compare.py @@ -0,0 +1,704 @@ +"""Between- and within-subject comparisons, with the follow-ups they need. + +scipy supplies the distributions and the standard tests. Welch's ANOVA, +Games-Howell and Dunn's test are written out here -- scipy has none of them +before 1.15 and GLIDER supports 1.10 -- and are pinned against pingouin and +scikit-posthocs reference values in the tests. + +"auto" picks the test from the data and says why in ``TestResult.reason``: + +* two groups -- Welch's t when every group passes Shapiro-Wilk, else + Mann-Whitney U. Welch rather than Student because it costs nothing when the + variances are equal and is right when they are not; +* several groups -- one-way ANOVA + Tukey when normal with equal variances + (Levene), Welch's ANOVA + Games-Howell when normal with unequal variances, + Kruskal-Wallis + Dunn otherwise; +* two epochs, same animals -- paired t when the differences pass + Shapiro-Wilk, else Wilcoxon signed-rank. +""" + +from __future__ import annotations + +import inspect +import itertools +import math +import warnings +from functools import partial + +import numpy as np +import pandas as pd +from scipy import stats + +from glider.analysis.stats.corrections import METHODS, adjust +from glider.analysis.stats.frames import Dropped, between_cells, within_matrix +from glider.analysis.stats.result import TestResult, format_p + +__all__ = [ + "BETWEEN_TESTS", + "ONE_SAMPLE_TESTS", + "PAIRED_TESTS", + "compare_epochs", + "compare_groups", + "dunn", + "games_howell", + "run_between", + "run_one_sample", + "run_paired", + "tukey", +] + +#: The between-subject tests by key, and the name each reports under. +BETWEEN_TESTS = { + "auto": "Auto", + "welch": "Welch's t-test", + "student": "Student's t-test", + "mannwhitney": "Mann-Whitney U", + "permutation": "Permutation test (mean difference)", + "anova": "One-way ANOVA", + "welch_anova": "Welch's ANOVA", + "kruskal": "Kruskal-Wallis", +} +#: Tests that compare exactly two groups. +_TWO_ONLY = {"welch", "student", "mannwhitney", "permutation"} + +#: The within-subject (two epoch) tests by key. +PAIRED_TESTS = { + "auto": "Auto", + "paired_t": "Paired t-test", + "wilcoxon": "Wilcoxon signed-rank", + "permutation": "Permutation test (paired)", +} + +#: Below this, Shapiro-Wilk has too little power for "normal" to mean much. +_SMALL_N = 5 +_ALPHA = 0.05 +_RESAMPLES = 9999 +_SEED = 0 + + +# --------------------------------------------------------------------------- +# helpers +# --------------------------------------------------------------------------- + + +def _seeded() -> dict: + """A fixed generator for scipy's resampling, under whichever name it takes. + + Resampled p values and CIs must not change between two clicks on the + same data. ``rng`` arrived in scipy 1.15; ``random_state`` before it. + """ + generator = np.random.default_rng(_SEED) + params = inspect.signature(stats.bootstrap).parameters + return {"rng": generator} if "rng" in params else {"random_state": generator} + + +def _shapiro(values: np.ndarray) -> float: + """Shapiro-Wilk p, NaN under three values; 1.0 for constant data.""" + if len(values) < 3: + return float("nan") + if np.ptp(values) == 0: + return 1.0 + return float(stats.shapiro(values).pvalue) + + +def _levene(cells: list[np.ndarray]) -> float: + """Levene's test (median-centred, i.e. Brown-Forsythe) p, NaN if it cannot run.""" + if len(cells) < 2 or any(len(c) < 2 for c in cells): + return float("nan") + if all(np.ptp(c) == 0 for c in cells): + return 1.0 + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + return float(stats.levene(*cells, center="median").pvalue) + + +def _t_crit(dof: float) -> float: + return float(abs(stats.t.ppf(0.025, dof))) + + +def hedges_g(a: np.ndarray, b: np.ndarray) -> tuple[float, tuple[float, float]]: + """Hedges' g for ``a - b`` and its 95% CI (the R ``cohen.d`` / pingouin interval).""" + n1, n2 = len(a), len(b) + pooled = math.sqrt( + ((n1 - 1) * np.var(a, ddof=1) + (n2 - 1) * np.var(b, ddof=1)) / (n1 + n2 - 2) + ) + if pooled == 0: + return float("nan"), (float("nan"), float("nan")) + g = (np.mean(a) - np.mean(b)) / pooled * (1 - 3 / (4 * (n1 + n2) - 9)) + se = math.sqrt((n1 + n2) / (n1 * n2) + g**2 / (2 * (n1 + n2))) + crit = _t_crit(n1 + n2 - 2) + return float(g), (float(g - crit * se), float(g + crit * se)) + + +def cohen_dz(diff: np.ndarray) -> tuple[float, tuple[float, float]]: + """Cohen's d_z (mean difference over its SD) and its 95% CI.""" + n = len(diff) + sd = np.std(diff, ddof=1) + if sd == 0: + return float("nan"), (float("nan"), float("nan")) + dz = float(np.mean(diff) / sd) + se = math.sqrt(1 / n + dz**2 / (2 * n)) + crit = _t_crit(n - 1) + return dz, (dz - crit * se, dz + crit * se) + + +def _difference_ci(a, b, *, paired: bool) -> tuple[float, float] | None: + """Percentile bootstrap 95% CI of ``mean(a) - mean(b)``; None if it cannot run.""" + if len(a) < 2 or len(b) < 2: + return None + + def mean_diff(x, y, axis=-1): + return np.mean(x, axis=axis) - np.mean(y, axis=axis) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + ci = stats.bootstrap( + (a, b), + mean_diff, + paired=paired, + vectorized=True, + n_resamples=_RESAMPLES, + method="percentile", + **_seeded(), + ).confidence_interval + low, high = float(ci.low), float(ci.high) + return None if math.isnan(low) or math.isnan(high) else (low, high) + + +def _cells_warnings(n: dict[str, int]) -> list[str]: + out = [] + tiny = [name for name, count in n.items() if count < 3] + if tiny: + out.append(f"n < 3 in {', '.join(tiny)}: too few animals for a dependable test") + return out + + +def _failed(test: str, metric: str, n: dict[str, int], why: str, **extra) -> TestResult: + return TestResult( + test=test, + metric=metric, + statistic_name="", + statistic=float("nan"), + df=(), + p=float("nan"), + effect_name="", + effect=float("nan"), + effect_ci=None, + n=n, + warnings=[why], + **extra, + ) + + +# --------------------------------------------------------------------------- +# post-hocs: one row per pair, ``a, b, diff, statistic, p, p_adj`` +# --------------------------------------------------------------------------- + + +def _posthoc_frame(rows: list[dict]) -> pd.DataFrame: + return pd.DataFrame(rows, columns=["a", "b", "diff", "statistic", "p", "p_adj"]) + + +def tukey(cells: dict[str, np.ndarray]) -> pd.DataFrame: + """Tukey-Kramer HSD. ``statistic`` is t; p is already family-wise (p_adj = p).""" + names = list(cells) + k = len(names) + n_total = sum(len(v) for v in cells.values()) + dof = n_total - k + mse = sum(((v - v.mean()) ** 2).sum() for v in cells.values()) / dof + rows = [] + for a, b in itertools.combinations(names, 2): + x, y = cells[a], cells[b] + diff = float(x.mean() - y.mean()) + se = math.sqrt(mse * (1 / len(x) + 1 / len(y))) + t = diff / se if se else float("nan") + p = float(stats.studentized_range.sf(abs(t) * math.sqrt(2), k, dof)) + rows.append({"a": a, "b": b, "diff": diff, "statistic": t, "p": p, "p_adj": p}) + return _posthoc_frame(rows) + + +def games_howell(cells: dict[str, np.ndarray]) -> pd.DataFrame: + """Games-Howell: Tukey's idea with each pair's own variances and Welch df.""" + k = len(cells) + rows = [] + for a, b in itertools.combinations(list(cells), 2): + x, y = cells[a], cells[b] + va, vb = np.var(x, ddof=1) / len(x), np.var(y, ddof=1) / len(y) + diff = float(x.mean() - y.mean()) + se = math.sqrt(va + vb) + if se == 0: + t = p = float("nan") + else: + t = diff / se + dof = (va + vb) ** 2 / (va**2 / (len(x) - 1) + vb**2 / (len(y) - 1)) + p = float(stats.studentized_range.sf(abs(t) * math.sqrt(2), k, dof)) + rows.append({"a": a, "b": b, "diff": diff, "statistic": t, "p": p, "p_adj": p}) + return _posthoc_frame(rows) + + +def dunn(cells: dict[str, np.ndarray], correction: str = "holm") -> pd.DataFrame: + """Dunn's test on the Kruskal-Wallis ranks, tie-corrected; ``diff`` is in mean ranks.""" + names = list(cells) + pooled = np.concatenate([cells[name] for name in names]) + n_total = len(pooled) + ranks = stats.rankdata(pooled) + _, ties = np.unique(pooled, return_counts=True) + tie_term = float(((ties**3) - ties).sum()) / (12 * (n_total - 1)) if n_total > 1 else 0.0 + mean_rank, start = {}, 0 + for name in names: + size = len(cells[name]) + mean_rank[name] = float(ranks[start : start + size].mean()) + start += size + rows = [] + for a, b in itertools.combinations(names, 2): + diff = mean_rank[a] - mean_rank[b] + var = (n_total * (n_total + 1) / 12 - tie_term) * (1 / len(cells[a]) + 1 / len(cells[b])) + z = diff / math.sqrt(var) if var > 0 else float("nan") + p = float(2 * stats.norm.sf(abs(z))) if not math.isnan(z) else float("nan") + rows.append({"a": a, "b": b, "diff": diff, "statistic": z, "p": p}) + frame = pd.DataFrame(rows, columns=["a", "b", "diff", "statistic", "p"]) + frame["p_adj"] = adjust(frame["p"], correction) + return frame + + +# --------------------------------------------------------------------------- +# between subjects +# --------------------------------------------------------------------------- + + +def _welch_anova(cells: list[np.ndarray]) -> tuple[float, float, float]: + """Welch's F, its two df, for several groups with unequal variances.""" + k = len(cells) + n = np.array([len(c) for c in cells], dtype=float) + means = np.array([c.mean() for c in cells]) + variances = np.array([np.var(c, ddof=1) for c in cells]) + weights = n / variances + total = weights.sum() + grand = (weights * means).sum() / total + numerator = (weights * (means - grand) ** 2).sum() / (k - 1) + lam = ((1 - weights / total) ** 2 / (n - 1)).sum() + f = numerator / (1 + 2 * (k - 2) / (k**2 - 1) * lam) + return float(f), float(k - 1), float((k**2 - 1) / (3 * lam)) + + +def _eta_squared(cells: list[np.ndarray]) -> float: + pooled = np.concatenate(cells) + total = ((pooled - pooled.mean()) ** 2).sum() + between = sum(len(c) * (c.mean() - pooled.mean()) ** 2 for c in cells) + return float(between / total) if total else float("nan") + + +def _choose_between(cells: dict[str, np.ndarray], checks: dict[str, float]) -> tuple[str, str]: + """The test "auto" runs, and the sentence that says why.""" + shapiro = {name: checks.get(f"Shapiro {name}", float("nan")) for name in cells} + failing = {name: p for name, p in shapiro.items() if p < _ALPHA} + untested = [name for name, p in shapiro.items() if math.isnan(p)] + note = f" ({', '.join(untested)} too small to check)" if untested else "" + two = len(cells) == 2 + if failing: + worst = min(failing, key=failing.get) + why = f"Shapiro-Wilk {format_p(failing[worst])} in {worst}: not normal" + return ("mannwhitney" if two else "kruskal"), why + if two: + return "welch", f"Shapiro-Wilk p ≥ .05 in every group{note}" + levene = checks.get("Levene", float("nan")) + if levene < _ALPHA: + return "welch_anova", f"normal{note}, but Levene {format_p(levene)}: unequal variances" + return "anova", f"normal{note}, equal variances (Levene {format_p(levene)})" + + +def run_between( + cells: dict[str, np.ndarray], + metric: str = "", + test: str = "auto", + *, + correction: str = "holm", + dropped: list[Dropped] | None = None, +) -> TestResult: + """Compare independent groups: ``cells`` maps group -> values, one per animal.""" + if test not in BETWEEN_TESTS: + raise ValueError(f"unknown test {test!r}; one of {sorted(BETWEEN_TESTS)}") + if correction not in METHODS: + raise ValueError(f"unknown correction {correction!r}; one of {sorted(METHODS)}") + cells = {name: np.asarray(v, dtype=float) for name, v in cells.items() if len(v)} + n = {name: len(v) for name, v in cells.items()} + extra = {"dropped": list(dropped or [])} + label = BETWEEN_TESTS[test] + if len(cells) < 2: + return _failed(label, metric, n, "needs at least two groups with data", **extra) + if test in _TWO_ONLY and len(cells) != 2: + return _failed(label, metric, n, f"{label} compares two groups, not {len(cells)}", **extra) + if any(count < 2 for count in n.values()): + small = ", ".join(name for name, count in n.items() if count < 2) + return _failed(label, metric, n, f"needs two or more animals per group ({small})", **extra) + + checks = {f"Shapiro {name}": _shapiro(v) for name, v in cells.items()} + checks["Levene"] = _levene(list(cells.values())) + reason = "" + if test == "auto": + test, reason = _choose_between(cells, checks) + result = _between(test, cells, metric, n, correction) + result.reason = reason + result.checks = checks + result.dropped = extra["dropped"] + result.warnings = _cells_warnings(n) + result.warnings + if reason and test in ("welch", "anova", "welch_anova") and min(n.values()) < _SMALL_N: + result.warnings.append( + f"normality was checked on fewer than {_SMALL_N} animals per group, " + "where Shapiro-Wilk rarely rejects; a rank test is the safer read" + ) + return result + + +def _between(test, cells, metric, n, correction) -> TestResult: + names = list(cells) + values = list(cells.values()) + common = {"test": BETWEEN_TESTS[test], "metric": metric, "n": n} + if test in _TWO_ONLY: + a, b = values + g, g_ci = hedges_g(a, b) + difference = { + "difference": float(a.mean() - b.mean()), + "difference_ci": _difference_ci(a, b, paired=False), + } + if test in ("welch", "student"): + equal = test == "student" + res = stats.ttest_ind(a, b, equal_var=equal) + if equal: + dof = len(a) + len(b) - 2 + else: + va, vb = np.var(a, ddof=1) / len(a), np.var(b, ddof=1) / len(b) + dof = (va + vb) ** 2 / (va**2 / (len(a) - 1) + vb**2 / (len(b) - 1)) + return TestResult( + **common, + statistic_name="t", + statistic=float(res.statistic), + df=(float(dof),), + p=float(res.pvalue), + effect_name="Hedges' g", + effect=g, + effect_ci=g_ci, + **difference, + ) + if test == "mannwhitney": + res = stats.mannwhitneyu(a, b, alternative="two-sided") + rbc = 2 * float(res.statistic) / (len(a) * len(b)) - 1 + return TestResult( + **common, + statistic_name="U", + statistic=float(res.statistic), + df=(), + p=float(res.pvalue), + effect_name="rank-biserial r", + effect=rbc, + effect_ci=None, + **difference, + ) + res = stats.permutation_test( + (a, b), + lambda x, y, axis=-1: np.mean(x, axis=axis) - np.mean(y, axis=axis), + permutation_type="independent", + vectorized=True, + n_resamples=_RESAMPLES, + **_seeded(), + ) + return TestResult( + **common, + statistic_name="Δmean", + statistic=float(res.statistic), + df=(), + p=float(res.pvalue), + effect_name="Hedges' g", + effect=g, + effect_ci=g_ci, + **difference, + ) + + k, n_total = len(values), sum(n.values()) + if test == "anova": + res = stats.f_oneway(*values) + stat, dfs, posthoc, name = float(res.statistic), (k - 1, n_total - k), tukey, "Tukey HSD" + effect_name, effect = "η²", _eta_squared(values) + elif test == "welch_anova": + f, df1, df2 = _welch_anova(values) + stat, dfs = f, (df1, df2) + posthoc, name = games_howell, "Games-Howell" + effect_name, effect = "η²", _eta_squared(values) + else: # kruskal + res = stats.kruskal(*values) + stat, dfs = float(res.statistic), (k - 1,) + posthoc = partial(dunn, correction=correction) + name = f"Dunn ({METHODS[correction]})" + effect_name, effect = "ε²", stat / (n_total - 1) + dist = stats.chi2(*dfs) if test == "kruskal" else stats.f(*dfs) + p = float(dist.sf(stat)) if not math.isnan(stat) else float("nan") + result = TestResult( + **common, + statistic_name="H" if test == "kruskal" else "F", + statistic=stat, + df=tuple(float(d) for d in dfs), + p=p, + effect_name=effect_name, + effect=float(effect), + effect_ci=None, + ) + if k > 2: + result.posthoc = posthoc(dict(zip(names, values, strict=True))) + result.posthoc_name = name + return result + + +def compare_groups( + tidy: pd.DataFrame, + metric: str, + epoch: str | None = None, + test: str = "auto", + *, + correction: str = "holm", +) -> TestResult: + """Groups against each other on ``metric`` within one epoch of the tidy frame.""" + cells, _sessions, dropped = between_cells(tidy, metric, epoch) + return run_between(cells, metric, test, correction=correction, dropped=dropped) + + +# --------------------------------------------------------------------------- +# within subjects: the same animals in two epochs +# --------------------------------------------------------------------------- + + +def _matched_rank_biserial(diff: np.ndarray) -> float: + nonzero = diff[diff != 0] + if not len(nonzero): + return float("nan") + ranks = stats.rankdata(np.abs(nonzero)) + plus, minus = ranks[nonzero > 0].sum(), ranks[nonzero < 0].sum() + return float((plus - minus) / (plus + minus)) + + +def _location( + values: np.ndarray, + test: str, + *, + names: dict[str, str], + metric: str, + n: dict[str, int], + dropped: list[Dropped], + difference_ci, +) -> TestResult: + """A paired or one-sample test: both are a test of ``values`` (the differences, + or the values less ``mu``) against zero. + + ``names`` gives the test, effect and "differences" wording for the kind. + """ + checks = {f"Shapiro {names['values']}": _shapiro(values)} + reason = "" + if test == "auto": + p = checks[f"Shapiro {names['values']}"] + if p < _ALPHA: + test, reason = "wilcoxon", f"Shapiro-Wilk {format_p(p)} on the {names['values']}" + else: + test = "t" + reason = ( + f"{names['values']} too few to check for normality" + if math.isnan(p) + else f"Shapiro-Wilk p ≥ .05 on the {names['values']}" + ) + common = { + "test": names[test], + "metric": metric, + "n": n, + "difference": float(values.mean()), + "difference_ci": difference_ci, + "reason": reason, + "checks": checks, + "dropped": dropped, + } + d, d_ci = cohen_dz(values) + if test == "t": + res = stats.ttest_1samp(values, 0.0) + result = TestResult( + **common, + statistic_name="t", + statistic=float(res.statistic), + df=(float(len(values) - 1),), + p=float(res.pvalue), + effect_name=names["effect"], + effect=d, + effect_ci=d_ci, + ) + elif test == "wilcoxon": + if not np.any(values): + why = f"every {names['values'][:-1]} is zero" + return _failed(names[test], metric, n, why, dropped=dropped, checks=checks) + res = stats.wilcoxon(values) + result = TestResult( + **common, + statistic_name="W", + statistic=float(res.statistic), + df=(), + p=float(res.pvalue), + effect_name="rank-biserial r", + effect=_matched_rank_biserial(values), + effect_ci=None, + ) + else: + # One sample with permutation_type="samples" flips signs: the same as + # swapping each animal's two epochs. + res = stats.permutation_test( + (values,), + lambda x, axis=-1: np.mean(x, axis=axis), + permutation_type="samples", + vectorized=True, + n_resamples=_RESAMPLES, + **_seeded(), + ) + result = TestResult( + **common, + statistic_name="Δmean", + statistic=float(res.statistic), + df=(), + p=float(res.pvalue), + effect_name=names["effect"], + effect=d, + effect_ci=d_ci, + ) + result.warnings = _cells_warnings({"animals": len(values)}) + result.warnings + return result + + +_PAIRED_KIND = { + "paired_t": "t", + "wilcoxon": "wilcoxon", + "permutation": "permutation", + "auto": "auto", +} + + +def run_paired( + a, + b, + names: tuple[str, str] = ("A", "B"), + metric: str = "", + test: str = "auto", + *, + dropped: list[Dropped] | None = None, +) -> TestResult: + """Compare two epochs of the same animals: ``a[i]`` and ``b[i]`` are one animal.""" + if test not in PAIRED_TESTS: + raise ValueError(f"unknown test {test!r}; one of {sorted(PAIRED_TESTS)}") + a, b = np.asarray(a, dtype=float), np.asarray(b, dtype=float) + if a.shape != b.shape: + raise ValueError("paired samples must be the same length") + n = {names[0]: len(a), names[1]: len(b)} + if len(a) < 2: + return _failed( + PAIRED_TESTS[test], metric, n, "needs two or more animals", dropped=list(dropped or []) + ) + return _location( + a - b, + _PAIRED_KIND[test], + names={ + "values": "differences", + "effect": "Cohen's d_z", + "t": PAIRED_TESTS["paired_t"], + "wilcoxon": PAIRED_TESTS["wilcoxon"], + "permutation": PAIRED_TESTS["permutation"], + }, + metric=metric, + n=n, + dropped=list(dropped or []), + difference_ci=_difference_ci(a, b, paired=True), + ) + + +#: One-sample tests against a fixed value (a slope against 0, a ratio against 1). +ONE_SAMPLE_TESTS = { + "auto": "Auto", + "t": "One-sample t-test", + "wilcoxon": "Wilcoxon signed-rank", + "permutation": "Permutation test (sign-flip)", +} + + +def run_one_sample( + values, + mu: float = 0.0, + metric: str = "", + test: str = "auto", + *, + name: str = "", + dropped: list[Dropped] | None = None, +) -> TestResult: + """Whether ``values`` -- one per animal -- differ from ``mu``. + + ``difference`` is the mean less ``mu``, with its bootstrap 95% CI. + """ + if test not in ONE_SAMPLE_TESTS: + raise ValueError(f"unknown test {test!r}; one of {sorted(ONE_SAMPLE_TESTS)}") + values = np.asarray(values, dtype=float) + n = {name or metric or "animals": len(values)} + if len(values) < 2: + return _failed( + ONE_SAMPLE_TESTS[test], + metric, + n, + "needs two or more animals", + dropped=list(dropped or []), + ) + shifted = values - mu + return _location( + shifted, + test, + names={ + "values": "values", + "effect": "Cohen's d", + "t": ONE_SAMPLE_TESTS["t"], + "wilcoxon": ONE_SAMPLE_TESTS["wilcoxon"], + "permutation": ONE_SAMPLE_TESTS["permutation"], + }, + metric=metric, + n=n, + dropped=list(dropped or []), + difference_ci=_mean_ci(shifted), + ) + + +def _mean_ci(values: np.ndarray) -> tuple[float, float] | None: + """Percentile bootstrap 95% CI of the mean; None if it cannot run.""" + if len(values) < 2: + return None + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + ci = stats.bootstrap( + (values,), + np.mean, + vectorized=True, + n_resamples=_RESAMPLES, + method="percentile", + **_seeded(), + ).confidence_interval + low, high = float(ci.low), float(ci.high) + return None if math.isnan(low) or math.isnan(high) else (low, high) + + +def compare_epochs( + tidy: pd.DataFrame, + metric: str, + epochs: list[str], + test: str = "auto", + *, + group: str | None = None, +) -> TestResult: + """The same animals across epochs of the tidy frame, optionally one group only. + + Two epochs are a paired test (:data:`PAIRED_TESTS`); three or more a + repeated-measures one (:data:`~glider.analysis.stats.repeated.REPEATED_TESTS`). + """ + if len(epochs) < 2: + raise ValueError("compare_epochs needs two or more epochs") + matrix, _sessions, dropped = within_matrix(tidy, metric, list(epochs), group) + if len(epochs) > 2: + from glider.analysis.stats.repeated import run_repeated + + return run_repeated(matrix, list(epochs), metric, test, dropped=dropped) + return run_paired(matrix[:, 0], matrix[:, 1], tuple(epochs), metric, test, dropped=dropped) diff --git a/src/glider/analysis/stats/corrections.py b/src/glider/analysis/stats/corrections.py new file mode 100644 index 00000000..359394de --- /dev/null +++ b/src/glider/analysis/stats/corrections.py @@ -0,0 +1,46 @@ +"""Multiple-comparison corrections of a family of p values. + +NaN p values (a comparison that could not be run) pass through as NaN and do +not count towards the family size. +""" + +from __future__ import annotations + +import numpy as np + +__all__ = ["METHODS", "adjust"] + +#: Accepted ``method`` names, and how a methods line names each. +METHODS = { + "bonferroni": "Bonferroni", + "holm": "Holm", + "fdr_bh": "Benjamini-Hochberg", +} + + +def adjust(p_values, method: str = "holm") -> np.ndarray: + """Adjusted p values, in the input order, capped at 1.""" + p = np.asarray(p_values, dtype=float) + out = np.full(p.shape, np.nan) + ok = ~np.isnan(p) + m = int(ok.sum()) + if m == 0: + return out + if method not in METHODS: + raise ValueError(f"unknown correction {method!r}; one of {sorted(METHODS)}") + pv = p[ok] + if method == "bonferroni": + adjusted = pv * m + else: + order = np.argsort(pv, kind="stable") + ranked = pv[order] + if method == "holm": + # Step-down: the k-th smallest times (m - k), never below an earlier one. + stepped = np.maximum.accumulate(ranked * (m - np.arange(m))) + else: + # Step-up: the k-th smallest times m / k, never above a later one. + stepped = np.minimum.accumulate((ranked * m / np.arange(1, m + 1))[::-1])[::-1] + adjusted = np.empty(m) + adjusted[order] = stepped + out[ok] = np.minimum(adjusted, 1.0) + return out diff --git a/src/glider/analysis/stats/frames.py b/src/glider/analysis/stats/frames.py new file mode 100644 index 00000000..1bc5e5c7 --- /dev/null +++ b/src/glider/analysis/stats/frames.py @@ -0,0 +1,98 @@ +"""From the epoch table's tidy frame to the cells a test compares. + +Two shapes: *between* (one epoch, the animals split by group) and *within* +(the same animals across epochs, as an animals x epochs matrix). Every animal +that cannot contribute is named with the reason, never silently dropped: a +reviewer reading "n = 7" when the cohort has 8 needs to know which one and why. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass + +import numpy as np +import pandas as pd + +__all__ = ["UNGROUPED", "Dropped", "between_cells", "within_matrix"] + +#: The group of an animal with none in ``glider_project.json``. +UNGROUPED = "Ungrouped" + + +@dataclass(frozen=True) +class Dropped: + """An animal left out of a test, and why.""" + + session: str + reason: str + + +def _missing(value) -> bool: + try: + return value is None or math.isnan(float(value)) + except (TypeError, ValueError): + return True + + +def _group(value) -> str: + return UNGROUPED if value is None or str(value) in ("", "nan") else str(value) + + +def between_cells( + tidy: pd.DataFrame, metric: str, epoch: str | None = None +) -> tuple[dict[str, np.ndarray], dict[str, list[str]], list[Dropped]]: + """The metric per group within one epoch. + + Returns ``(values, sessions, dropped)``: each group's values and the + sessions they came from, in frame order, groups sorted by name with + Ungrouped last. ``epoch`` None takes the frame as a single epoch. + """ + rows = tidy if epoch is None else tidy[tidy["range"] == epoch] + values: dict[str, list[float]] = {} + sessions: dict[str, list[str]] = {} + dropped = [] + for _, row in rows.iterrows(): + value = row.get(metric) + if _missing(value): + where = f" in {epoch}" if epoch is not None else "" + dropped.append(Dropped(str(row["session"]), f"no {metric}{where}")) + continue + group = _group(row.get("group")) + values.setdefault(group, []).append(float(value)) + sessions.setdefault(group, []).append(str(row["session"])) + order = sorted(values, key=lambda g: (g == UNGROUPED, g)) + return ( + {g: np.asarray(values[g], dtype=float) for g in order}, + {g: sessions[g] for g in order}, + dropped, + ) + + +def within_matrix( + tidy: pd.DataFrame, metric: str, epochs: list[str], group: str | None = None +) -> tuple[np.ndarray, list[str], list[Dropped]]: + """The metric as an animals x epochs matrix, complete animals only. + + A within-subject test needs every animal in every epoch, so an animal + missing any one is left out entirely (listwise), and named. ``group`` + restricts to one group's animals. + """ + rows = tidy[tidy["range"].isin(epochs)] + if group is not None: + rows = rows[rows["group"].map(_group) == group] + table: dict[str, dict[str, float]] = {} + for _, row in rows.iterrows(): + cells = table.setdefault(str(row["session"]), {}) + value = row.get(metric) + cells[str(row["range"])] = float("nan") if _missing(value) else float(value) + kept, matrix, dropped = [], [], [] + for session, cells in table.items(): + missing = [e for e in epochs if _missing(cells.get(e))] + if missing: + dropped.append(Dropped(session, f"no {metric} in {', '.join(missing)}")) + continue + kept.append(session) + matrix.append([cells[e] for e in epochs]) + out = np.asarray(matrix, dtype=float).reshape(len(kept), len(epochs)) + return out, kept, dropped diff --git a/src/glider/analysis/stats/repeated.py b/src/glider/analysis/stats/repeated.py new file mode 100644 index 00000000..512797f4 --- /dev/null +++ b/src/glider/analysis/stats/repeated.py @@ -0,0 +1,483 @@ +"""Several epochs of the same animals, alone or split by group. + +* :func:`run_repeated` -- one-way repeated-measures ANOVA, or Friedman, over + an animals x epochs matrix, with pairwise follow-ups. +* :func:`run_mixed` -- the split-plot (mixed) ANOVA: group between animals, + epoch within them. Group, epoch, and the interaction -- the effect a + ChR2-vs-eYFP stimulation experiment is usually about. + +Sphericity follows R (``mauchly.test``, ``car::Anova``, afex) and SPSS: +Mauchly's test and the Greenhouse-Geisser and Huynh-Feldt epsilons come from +the covariance of orthonormal contrasts -- in a mixed design, the *pooled +within-group* covariance, so a difference between groups is not mistaken for +a departure from sphericity. When Mauchly rejects, the within-subject p +values reported are Greenhouse-Geisser corrected, and the result says so. +""" + +from __future__ import annotations + +import itertools +import math +from dataclasses import dataclass, field + +import numpy as np +import pandas as pd +from scipy import stats + +from glider.analysis.stats.compare import run_between, run_paired +from glider.analysis.stats.corrections import METHODS, adjust +from glider.analysis.stats.frames import Dropped, between_cells, within_matrix +from glider.analysis.stats.result import TestResult, format_p + +__all__ = [ + "REPEATED_TESTS", + "MixedResult", + "Sphericity", + "pairwise_within", + "run_mixed", + "compare_mixed", + "run_repeated", + "sphericity", +] + +REPEATED_TESTS = { + "auto": "Auto", + "rm_anova": "Repeated-measures ANOVA", + "friedman": "Friedman test", +} + +_ALPHA = 0.05 + + +# --------------------------------------------------------------------------- +# sphericity +# --------------------------------------------------------------------------- + + +def _contrasts(k: int) -> np.ndarray: + """A (k, k-1) orthonormal basis orthogonal to the grand mean.""" + centred = np.eye(k) - 1.0 / k + q, _ = np.linalg.qr(centred) + return q[:, : k - 1] + + +@dataclass(frozen=True) +class Sphericity: + """Mauchly's test and the two epsilon corrections for one within factor.""" + + w: float + chi2: float + df: float + p: float + gg: float # Greenhouse-Geisser epsilon + hf: float # Huynh-Feldt epsilon (Lecoutre's correction for mixed designs) + + @property + def violated(self) -> bool: + return not math.isnan(self.p) and self.p < _ALPHA + + +def sphericity(residuals: np.ndarray, df_resid: int) -> Sphericity: + """Sphericity of a (subjects, k) matrix of within-group residuals. + + ``df_resid`` is n - 1 for a single group and N - groups in a mixed design. + With two levels sphericity cannot fail: W = 1, epsilon = 1. + """ + k = residuals.shape[1] + d = k - 1 + if d <= 1: + return Sphericity(1.0, 0.0, 0.0, 1.0, 1.0, 1.0) + contrasts = residuals @ _contrasts(k) + cov = np.cov(contrasts, rowvar=False, ddof=1) + trace = float(np.trace(cov)) + gg = min(trace**2 / (d * float(np.trace(cov @ cov))), 1.0) if trace > 0 else 1.0 + n_plus = df_resid + 1 # n for one group; N - groups + 1 in general (Lecoutre) + hf_denominator = d * (df_resid - d * gg) + hf = (n_plus * d * gg - 2) / hf_denominator if hf_denominator > 0 else 1.0 + hf = float(min(hf, 1.0)) + # Mauchly's W, and its chi-square approximation with R's second-order term. + if df_resid < d: + # Fewer animals than contrasts: the covariance is singular. + return Sphericity(float("nan"), float("nan"), float("nan"), float("nan"), gg, hf) + sign, logdet = np.linalg.slogdet(cov) + log_w = logdet - d * math.log(trace / d) if sign > 0 and trace > 0 else -math.inf + w = math.exp(log_w) + dof = d * (d + 1) / 2 - 1 + f = 1 - (2 * d**2 + d + 2) / (6 * d * df_resid) + w2 = (d + 2) * (d - 1) * (d - 2) * (2 * d**3 + 6 * d**2 + 3 * k + 2) + w2 /= 288 * (df_resid * d * f) ** 2 + chi = -df_resid * f * log_w + p1, p2 = stats.chi2.sf(chi, dof), stats.chi2.sf(chi, dof + 4) + return Sphericity(w, float(chi), float(dof), float(p1 + w2 * (p2 - p1)), gg, hf) + + +def _corrected_p(f: float, df1: float, df2: float, eps: float) -> float: + """The F test's p on epsilon-scaled df (floored at 1, as R and pingouin do).""" + return float(stats.f.sf(f, max(df1 * eps, 1.0), max(df2 * eps, 1.0))) + + +# --------------------------------------------------------------------------- +# pairwise follow-ups +# --------------------------------------------------------------------------- + + +def pairwise_within( + matrix: np.ndarray, epochs: list[str], *, rank: bool = False, correction: str = "holm" +) -> pd.DataFrame: + """Every pair of epochs, paired t (or Wilcoxon when ``rank``), adjusted together.""" + rows = [] + for i, j in itertools.combinations(range(len(epochs)), 2): + result = run_paired( + matrix[:, i], + matrix[:, j], + (epochs[i], epochs[j]), + test="wilcoxon" if rank else "paired_t", + ) + rows.append( + { + "a": epochs[i], + "b": epochs[j], + "diff": float(np.mean(matrix[:, i] - matrix[:, j])), + "statistic": result.statistic, + "p": result.p, + "effect": result.effect, + } + ) + frame = pd.DataFrame(rows, columns=["a", "b", "diff", "statistic", "p", "effect"]) + frame["p_adj"] = adjust(frame["p"], correction) + return frame + + +# --------------------------------------------------------------------------- +# one group, several epochs +# --------------------------------------------------------------------------- + + +def _rm_sums(matrix: np.ndarray) -> tuple[float, float, float]: + """``(SS_epoch, SS_error, SS_subjects)`` of a one-way repeated-measures design.""" + n, k = matrix.shape + grand = matrix.mean() + ss_epoch = n * float(((matrix.mean(axis=0) - grand) ** 2).sum()) + ss_subjects = k * float(((matrix.mean(axis=1) - grand) ** 2).sum()) + ss_total = float(((matrix - grand) ** 2).sum()) + return ss_epoch, ss_total - ss_epoch - ss_subjects, ss_subjects + + +def _residuals(matrix: np.ndarray) -> np.ndarray: + """What the additive subject + epoch model leaves: what normality is about.""" + return matrix - matrix.mean(axis=1, keepdims=True) - matrix.mean(axis=0) + matrix.mean() + + +def _shapiro(values: np.ndarray) -> float: + values = np.ravel(values) + if len(values) < 3 or np.ptp(values) < 1e-12 * max(1.0, float(np.abs(values).max())): + return float("nan") if len(values) < 3 else 1.0 + return float(stats.shapiro(values).pvalue) + + +def run_repeated( + matrix, + epochs: list[str], + metric: str = "", + test: str = "auto", + *, + correction: str = "holm", + dropped: list[Dropped] | None = None, +) -> TestResult: + """Several epochs of the same animals: row = animal, column = epoch, complete.""" + if test not in REPEATED_TESTS: + raise ValueError(f"unknown test {test!r}; one of {sorted(REPEATED_TESTS)}") + if correction not in METHODS: + raise ValueError(f"unknown correction {correction!r}; one of {sorted(METHODS)}") + matrix = np.asarray(matrix, dtype=float) + n, k = matrix.shape + counts = dict.fromkeys(epochs, n) + base = {"metric": metric, "n": counts, "dropped": list(dropped or [])} + if k < 3: + raise ValueError("run_repeated wants three or more epochs; two are a paired test") + if n < 2: + return TestResult( + test=REPEATED_TESTS[test], + statistic_name="", + statistic=float("nan"), + df=(), + p=float("nan"), + effect_name="", + effect=float("nan"), + effect_ci=None, + warnings=["needs two or more animals with every epoch"], + **base, + ) + ss_epoch, ss_error, _ = _rm_sums(matrix) + df1, df2 = k - 1, (n - 1) * (k - 1) + spher = sphericity(matrix - matrix.mean(axis=0), n - 1) + checks = {"Shapiro residuals": _shapiro(_residuals(matrix)), "Mauchly": spher.p} + reason = "" + if test == "auto": + p = checks["Shapiro residuals"] + if p < _ALPHA: + test, reason = "friedman", f"Shapiro-Wilk {format_p(p)} on the residuals" + else: + test, reason = "rm_anova", "Shapiro-Wilk p ≥ .05 on the residuals" + warnings = [] + if n < 3: + warnings.append("n < 3 animals: too few for a dependable test") + if test == "rm_anova": + f = (ss_epoch / df1) / (ss_error / df2) if ss_error > 0 else float("nan") + p_unc = float(stats.f.sf(f, df1, df2)) if not math.isnan(f) else float("nan") + p, dfs = p_unc, (float(df1), float(df2)) + if spher.violated: + p = _corrected_p(f, df1, df2, spher.gg) + dfs = (df1 * spher.gg, df2 * spher.gg) + reason = (reason + "; " if reason else "") + ( + f"Mauchly {format_p(spher.p)}: Greenhouse-Geisser corrected (ε = {spher.gg:.3f})" + ) + elif math.isnan(spher.p): + warnings.append("too few animals to test sphericity; the uncorrected p assumes it") + result = TestResult( + test=REPEATED_TESTS[test], + statistic_name="F", + statistic=float(f), + df=dfs, + p=p, + effect_name="partial η²", + effect=ss_epoch / (ss_epoch + ss_error) if ss_epoch + ss_error else float("nan"), + effect_ci=None, + reason=reason, + checks=checks, + warnings=warnings, + **base, + ) + result.details = { + "p uncorrected": p_unc, + "p GG": _corrected_p(f, df1, df2, spher.gg), + "p HF": _corrected_p(f, df1, df2, spher.hf), + "GG ε": spher.gg, + "HF ε": spher.hf, + "Mauchly W": spher.w, + } + result.posthoc = pairwise_within(matrix, epochs, correction=correction) + result.posthoc_name = f"Paired t-tests ({METHODS[correction]})" + return result + res = stats.friedmanchisquare(*matrix.T) + chi = float(res.statistic) + result = TestResult( + test=REPEATED_TESTS[test], + statistic_name="χ²", + statistic=chi, + df=(float(k - 1),), + p=float(res.pvalue), + effect_name="Kendall's W", + effect=chi / (n * (k - 1)), + effect_ci=None, + reason=reason, + checks=checks, + warnings=warnings, + **base, + ) + result.posthoc = pairwise_within(matrix, epochs, rank=True, correction=correction) + result.posthoc_name = f"Wilcoxon signed-rank ({METHODS[correction]})" + return result + + +# --------------------------------------------------------------------------- +# group x epoch +# --------------------------------------------------------------------------- + + +@dataclass +class MixedResult: + """A mixed ANOVA: the three effects, sphericity, and simple effects. + + ``effects`` maps ``"group"``, ``factor`` and ``"group × "`` to a + :class:`TestResult` each; ``factor`` is ``"epoch"`` unless the caller + names it (``"bin"`` for time bins). ``simple_effects`` is filled only when the + interaction is significant: the epoch effect within each group, then the + group effect within each epoch, Holm-adjusted together. + """ + + metric: str + groups: list[str] + epochs: list[str] + n: dict[str, int] + effects: dict[str, TestResult] + sphericity: Sphericity + simple_effects: pd.DataFrame | None = None + warnings: list[str] = field(default_factory=list) + dropped: list[Dropped] = field(default_factory=list) + factor: str = "epoch" + + def table(self) -> pd.DataFrame: + """The ANOVA table, one row per effect.""" + rows = [] + for name, r in self.effects.items(): + rows.append( + { + "effect": name, + "F": r.statistic, + "df1": r.df[0] if r.df else float("nan"), + "df2": r.df[1] if len(r.df) > 1 else float("nan"), + "p": r.p, + "p_uncorrected": r.details.get("p uncorrected", r.p), + "partial_eta2": r.effect, + } + ) + return pd.DataFrame(rows) + + def methods(self) -> str: + cells = ", ".join(f"{g} (n = {n})" for g, n in self.n.items()) + parts = [ + f"{name}: {r.statistic_text()}, {format_p(r.p)}, partial η² = {r.effect:.2f}" + for name, r in self.effects.items() + ] + corrected = ( + f" Greenhouse-Geisser corrected (Mauchly {format_p(self.sphericity.p)}," + f" ε = {self.sphericity.gg:.3f})." + if self.sphericity.violated + else "" + ) + return ( + f"Mixed ANOVA on {self.metric}, group × {self.factor} ({', '.join(self.epochs)}), " + f"{cells}: " + "; ".join(parts) + "." + corrected + ) + + +def run_mixed( + cells: dict[str, np.ndarray], + epochs: list[str], + metric: str = "", + *, + correction: str = "holm", + dropped: list[Dropped] | None = None, + factor: str = "epoch", +) -> MixedResult: + """Group (between) x epoch (within). ``cells`` maps group -> animals x epochs matrix. + + ``factor`` names the within-subject factor in the effects and the text. + """ + cells = {g: np.asarray(m, dtype=float) for g, m in cells.items() if len(m)} + groups = list(cells) + n = {g: m.shape[0] for g, m in cells.items()} + k = len(epochs) + if len(groups) < 2 or k < 2 or any(m.shape[1] != k for m in cells.values()): + raise ValueError("run_mixed wants two or more groups, each an animals x epochs matrix") + data = np.vstack([cells[g] for g in groups]) + n_total, n_groups = data.shape[0], len(groups) + grand = data.mean() + ss_total = float(((data - grand) ** 2).sum()) + ss_subjects = k * float(((data.mean(axis=1) - grand) ** 2).sum()) + ss_group = k * sum(n[g] * (cells[g].mean() - grand) ** 2 for g in groups) + ss_epoch = n_total * float(((data.mean(axis=0) - grand) ** 2).sum()) + ss_cells = sum(n[g] * float(((cells[g].mean(axis=0) - grand) ** 2).sum()) for g in groups) + ss_inter = ss_cells - ss_group - ss_epoch + ss_err_between = ss_subjects - ss_group + ss_err_within = ss_total - ss_subjects - ss_epoch - ss_inter + df_group, df_err_between = n_groups - 1, n_total - n_groups + df_epoch, df_inter = k - 1, (n_groups - 1) * (k - 1) + df_err_within = (n_total - n_groups) * (k - 1) + + residuals = np.vstack([cells[g] - cells[g].mean(axis=0) for g in groups]) + spher = sphericity(residuals, df_err_between) + warnings = [] + if min(n.values()) < 3: + warnings.append("n < 3 in a group: too few animals for a dependable test") + if math.isnan(spher.p) and k > 2: + warnings.append("too few animals to test sphericity; the uncorrected p assumes it") + pieces = [_residuals(m).ravel() for m in cells.values() if m.shape[0] > 1] + p_resid = _shapiro(np.concatenate(pieces)) if pieces else float("nan") + if p_resid < _ALPHA: + warnings.append( + f"residuals are not normal (Shapiro-Wilk {format_p(p_resid)}); there is no " + "standard rank-based mixed ANOVA, so read the effects with care" + ) + + def effect(name, ss, df1, ss_err, df2, *, within: bool) -> TestResult: + f = (ss / df1) / (ss_err / df2) if ss_err > 0 else float("nan") + p_unc = float(stats.f.sf(f, df1, df2)) if not math.isnan(f) else float("nan") + p, dfs = p_unc, (float(df1), float(df2)) + if within and spher.violated: + p, dfs = _corrected_p(f, df1, df2, spher.gg), (df1 * spher.gg, df2 * spher.gg) + result = TestResult( + test="Mixed ANOVA", + metric=f"{metric} ({name})", + statistic_name="F", + statistic=float(f), + df=dfs, + p=p, + effect_name="partial η²", + effect=ss / (ss + ss_err) if ss + ss_err else float("nan"), + effect_ci=None, + n=dict(n), + checks={"Shapiro residuals": p_resid, "Mauchly": spher.p} if within else {}, + ) + result.details = {"p uncorrected": p_unc} + if within: + result.details.update( + { + "p GG": _corrected_p(f, df1, df2, spher.gg), + "p HF": _corrected_p(f, df1, df2, spher.hf), + } + ) + return result + + effects = { + "group": effect("group", ss_group, df_group, ss_err_between, df_err_between, within=False), + factor: effect(factor, ss_epoch, df_epoch, ss_err_within, df_err_within, within=True), + f"group × {factor}": effect( + f"group × {factor}", ss_inter, df_inter, ss_err_within, df_err_within, within=True + ), + } + result = MixedResult( + metric=metric, + groups=groups, + epochs=list(epochs), + n=n, + effects=effects, + sphericity=spher, + warnings=warnings, + dropped=list(dropped or []), + factor=factor, + ) + if effects[f"group × {factor}"].p < _ALPHA: + result.simple_effects = _simple_effects(cells, epochs, metric, correction, factor) + return result + + +def _simple_effects(cells, epochs, metric, correction, factor) -> pd.DataFrame: + """What the interaction is made of: epoch within each group, group within each epoch.""" + rows = [] + for g, matrix in cells.items(): + if len(epochs) == 2: + r = run_paired(matrix[:, 0], matrix[:, 1], tuple(epochs), metric, "paired_t") + else: + r = run_repeated(matrix, epochs, metric, "rm_anova", correction=correction) + rows.append({"effect": factor, "within": g, "test": r.test, "result": r}) + for j, epoch in enumerate(epochs): + column = {g: m[:, j] for g, m in cells.items()} + r = run_between(column, metric, "welch" if len(column) == 2 else "anova") + rows.append({"effect": "group", "within": epoch, "test": r.test, "result": r}) + frame = pd.DataFrame( + { + "effect": [row["effect"] for row in rows], + "within": [row["within"] for row in rows], + "test": [row["test"] for row in rows], + "statistic": [row["result"].statistic_text() for row in rows], + "p": [row["result"].p for row in rows], + } + ) + frame["p_adj"] = adjust(frame["p"], correction) + return frame + + +def compare_mixed( + tidy: pd.DataFrame, metric: str, epochs: list[str], *, correction: str = "holm" +) -> MixedResult: + """Group x epoch on ``metric`` from the tidy frame; incomplete animals are dropped.""" + groups = list(between_cells(tidy[tidy["range"].isin(epochs)], metric)[0]) + cells, dropped = {}, [] + for group in groups: + matrix, _kept, lost = within_matrix(tidy, metric, list(epochs), group=group) + dropped += lost + if len(matrix): + cells[group] = matrix + return run_mixed(cells, list(epochs), metric, correction=correction, dropped=dropped) diff --git a/src/glider/analysis/stats/result.py b/src/glider/analysis/stats/result.py new file mode 100644 index 00000000..abf95cd3 --- /dev/null +++ b/src/glider/analysis/stats/result.py @@ -0,0 +1,87 @@ +"""What a test returns: enough to report it, and enough to trust it.""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field + +import pandas as pd + +from glider.analysis.stats.frames import Dropped + +__all__ = ["TestResult", "format_p"] + + +def format_p(p: float) -> str: + """APA style: exact to three places, ``< .001`` below that.""" + if p is None or math.isnan(p): + return "p = n/a" + if p < 0.001: + return "p < .001" + return f"p = {p:.3f}".replace("0.", ".", 1) + + +def _num(value: float, places: int = 2) -> str: + return "n/a" if value is None or math.isnan(value) else f"{value:.{places}f}" + + +@dataclass +class TestResult: + """One test on one metric. + + ``df`` is ``()`` for a test without degrees of freedom (Mann-Whitney), + one value for t, two for F. ``effect_ci`` is ``None`` where no interval is + computed. ``posthoc`` holds pairwise follow-ups for a several-group test: + one row per pair with ``a, b, diff, statistic, p, p_adj``. + """ + + __test__ = False # not a pytest class, whatever its name + + test: str + metric: str + statistic_name: str + statistic: float + df: tuple[float, ...] + p: float + effect_name: str + effect: float + effect_ci: tuple[float, float] | None + n: dict[str, int] + #: Mean difference (first minus second) and its bootstrap 95% CI; two-cell tests only. + difference: float | None = None + difference_ci: tuple[float, float] | None = None + #: Why this test, when picked automatically; empty when the user chose it. + reason: str = "" + #: Assumption checks: ``"Shapiro "`` / ``"Levene"`` -> p. + checks: dict[str, float] = field(default_factory=dict) + warnings: list[str] = field(default_factory=list) + dropped: list[Dropped] = field(default_factory=list) + posthoc: pd.DataFrame | None = None + posthoc_name: str = "" + #: Further numbers a careful reader may want: the uncorrected and + #: epsilon-corrected p values and epsilons of a repeated-measures test. + details: dict[str, float] = field(default_factory=dict) + + @property + def ok(self) -> bool: + """Whether the test produced a p value at all.""" + return not math.isnan(self.p) + + def statistic_text(self) -> str: + """``t(12.31) = 2.41``, ``F(2, 21) = 5.10``, ``U = 12.0``.""" + if not self.df: + return f"{self.statistic_name} = {_num(self.statistic)}" + dfs = ", ".join(f"{d:g}" if float(d).is_integer() else f"{d:.2f}" for d in self.df) + return f"{self.statistic_name}({dfs}) = {_num(self.statistic)}" + + def methods(self) -> str: + """A sentence fit for a results section.""" + cells = ", ".join(f"{name} (n = {n})" for name, n in self.n.items()) + if not self.ok: + return f"{self.test} on {self.metric} not run: {'; '.join(self.warnings)}" + text = f"{self.test} on {self.metric}, {cells}: {self.statistic_text()}, {format_p(self.p)}" + if self.effect_name: + text += f", {self.effect_name} = {_num(self.effect)}" + if self.effect_ci is not None: + text += f" [95% CI {_num(self.effect_ci[0])}, {_num(self.effect_ci[1])}]" + return text + "." diff --git a/src/glider/analysis/stats/trend.py b/src/glider/analysis/stats/trend.py new file mode 100644 index 00000000..af2d94c3 --- /dev/null +++ b/src/glider/analysis/stats/trend.py @@ -0,0 +1,358 @@ +"""Change across time bins: does a measure drift over the session, and do groups differ? + +The input is the time-bin export's tidy frame (one row per animal per bin, +``bin, bin_start_s, bin_end_s, …``) or any tidy frame whose ``range`` +column is in time order. Three readings of one question, each with its own +strength: + +* **Slope per animal** -- an ordinary least-squares line through each + animal's bins, in units per minute. The slopes are then one number per + animal: tested against zero ("does it habituate?") and between groups + ("does one group habituate faster?"). Uses every bin an animal has. +* **Polynomial contrasts** -- linear and quadratic trend within the + repeated-measures design, each tested on its own error term (t² is the + contrast's F). The quadratic answers "does it level off?". +* **First vs last bin** -- the plain habituation index: change and ratio per + animal, the change tested as a paired comparison. + +Plus the repeated-measures ANOVA across bins for one group, and the group x +bin mixed ANOVA for several, from :mod:`.repeated`. + +A partial last bin (the session ended mid-bin) holds less time, so for a +total such as distance it reads as the animal slowing down -- a habituation +that is not there. Partial bins are left out unless asked for, and named. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field + +import numpy as np +import pandas as pd +from scipy import stats + +from glider.analysis.stats.compare import run_between, run_one_sample, run_paired +from glider.analysis.stats.corrections import METHODS, adjust +from glider.analysis.stats.frames import UNGROUPED, Dropped +from glider.analysis.stats.repeated import MixedResult, run_mixed, run_repeated +from glider.analysis.stats.result import TestResult, format_p + +__all__ = [ + "GroupTrend", + "TrendResult", + "animal_slopes", + "compare_trend", + "polynomial_contrasts", + "run_trend", + "trend_contrasts", +] + +#: A bin shorter than this fraction of the full width is partial. +_PARTIAL = 0.999 + + +# --------------------------------------------------------------------------- +# arithmetic on arrays +# --------------------------------------------------------------------------- + + +def polynomial_contrasts(times, degree: int = 2) -> np.ndarray: + """Orthonormal polynomial contrast weights for (possibly unequal) ``times``. + + ``(k, degree)``: column 0 linear, column 1 quadratic. Each sums to zero, + they are mutually orthogonal, and each is signed to rise with time (its + own power of time), so a positive contrast means "goes up" / "curves up". + """ + t = np.asarray(times, dtype=float) + degree = min(degree, len(t) - 1) + vander = np.vander(t - t.mean(), degree + 1, increasing=True) + q, _ = np.linalg.qr(vander) + weights = q[:, 1 : degree + 1] + for j in range(weights.shape[1]): + power = (t - t.mean()) ** (j + 1) + if weights[:, j] @ (power - power.mean()) < 0: + weights[:, j] = -weights[:, j] + return weights + + +def animal_slopes(matrix, times) -> pd.DataFrame: + """One OLS line per row of an animals x bins matrix, NaN bins skipped. + + Columns ``slope, intercept, r2, n_bins``; NaN where an animal has fewer + than three bins (two points always lie on a line, so their slope would + carry no information about fit). + """ + matrix = np.asarray(matrix, dtype=float) + t = np.asarray(times, dtype=float) + rows = [] + for values in matrix: + ok = np.isfinite(values) + if ok.sum() < 3: + rows.append( + {"slope": np.nan, "intercept": np.nan, "r2": np.nan, "n_bins": int(ok.sum())} + ) + continue + x, y = t[ok], values[ok] + if np.ptp(y) == 0: + slope, intercept, r2 = 0.0, float(y[0]), float("nan") + else: + fit = stats.linregress(x, y) + slope, intercept, r2 = float(fit.slope), float(fit.intercept), float(fit.rvalue**2) + rows.append({"slope": slope, "intercept": intercept, "r2": r2, "n_bins": int(ok.sum())}) + return pd.DataFrame(rows, columns=["slope", "intercept", "r2", "n_bins"]) + + +def trend_contrasts(matrix, times, *, correction: str = "holm") -> pd.DataFrame: + """Linear and quadratic trend over a complete animals x bins matrix. + + Each contrast is the one-sample t-test of the animals' contrast scores + against zero (Maxwell & Delaney's separate error term), so t² is the + contrast's F on (1, n - 1) df. ``partial_eta2 = t² / (t² + df)``. + """ + matrix = np.asarray(matrix, dtype=float) + weights = polynomial_contrasts(times) + names = ["linear", "quadratic"][: weights.shape[1]] + rows = [] + n = matrix.shape[0] + for name, w in zip(names, weights.T, strict=True): + scores = matrix @ w + if n < 2 or np.ptp(scores) == 0: + t = p = float("nan") + else: + res = stats.ttest_1samp(scores, 0.0) + t, p = float(res.statistic), float(res.pvalue) + rows.append( + { + "contrast": name, + "estimate": float(scores.mean()) if n else float("nan"), + "t": t, + "df": float(n - 1), + "p": p, + "partial_eta2": t**2 / (t**2 + n - 1) if not math.isnan(t) else float("nan"), + } + ) + frame = pd.DataFrame(rows, columns=["contrast", "estimate", "t", "df", "p", "partial_eta2"]) + frame["p_adj"] = adjust(frame["p"], correction) + return frame + + +# --------------------------------------------------------------------------- +# results +# --------------------------------------------------------------------------- + + +@dataclass +class GroupTrend: + """One group's trend, every reading of it.""" + + group: str + #: Slopes against zero; ``difference`` is the mean slope, per minute. + slope_test: TestResult + #: First bin against last, paired; ``difference`` is first minus last. + habituation: TestResult + #: Linear and quadratic contrasts over the complete animals; None under two bins. + contrasts: pd.DataFrame | None = None + #: Repeated-measures ANOVA across bins; None under three bins. + anova: TestResult | None = None + + +@dataclass +class TrendResult: + """Trend of one metric across time bins, per group and between groups.""" + + metric: str + bins: list[str] + #: Bin midpoints, in ``time_unit``. + times: list[float] + time_unit: str + #: One row per animal: ``session, group, slope, intercept, r2, n_bins, + #: first, last, change, ratio``. + per_animal: pd.DataFrame + groups: dict[str, GroupTrend] + #: Several groups only: slopes, and first-to-last change, compared between groups. + slope_comparison: TestResult | None = None + change_comparison: TestResult | None = None + #: Several groups only: group x bin, complete animals. + mixed: MixedResult | None = None + excluded_bins: list[str] = field(default_factory=list) + warnings: list[str] = field(default_factory=list) + dropped: list[Dropped] = field(default_factory=list) + + def methods(self) -> str: + """A short paragraph: what was fitted, and the headline numbers.""" + lines = [ + f"Trend of {self.metric} across {len(self.bins)} time bins " + f"({self.bins[0]} to {self.bins[-1]})" + + ( + f", excluding partial bin(s) {', '.join(self.excluded_bins)}" + if self.excluded_bins + else "" + ) + + "." + ] + for name, g in self.groups.items(): + slope = g.slope_test + if slope.ok: + lines.append( + f"{name}: mean slope {slope.difference:.3g} per {self.time_unit} " + f"({slope.test}, {slope.statistic_text()}, " + f"{format_p(slope.p)})." + ) + if self.slope_comparison is not None and self.slope_comparison.ok: + lines.append(f"Slopes between groups: {self.slope_comparison.methods()}") + if self.mixed is not None: + lines.append(self.mixed.methods()) + return " ".join(lines) + + +# --------------------------------------------------------------------------- +# entry points +# --------------------------------------------------------------------------- + + +def run_trend( + cells: dict[str, tuple[np.ndarray, list[str]]], + bins: list[str], + times, + metric: str = "", + *, + time_unit: str = "min", + correction: str = "holm", + dropped: list[Dropped] | None = None, +) -> TrendResult: + """``cells`` maps group -> (animals x bins matrix, NaN where missing; session ids).""" + if correction not in METHODS: + raise ValueError(f"unknown correction {correction!r}; one of {sorted(METHODS)}") + times = [float(t) for t in times] + if len(bins) < 2: + raise ValueError("a trend needs two or more bins") + dropped = list(dropped or []) + warnings: list[str] = [] + per_animal = [] + groups: dict[str, GroupTrend] = {} + complete: dict[str, np.ndarray] = {} + for group, (matrix, sessions) in cells.items(): + matrix = np.asarray(matrix, dtype=float) + slopes = animal_slopes(matrix, times) + first, last = matrix[:, 0], matrix[:, -1] + with np.errstate(divide="ignore", invalid="ignore"): + ratio = np.where(first != 0, last / first, np.nan) + for i, session in enumerate(sessions): + per_animal.append( + { + "session": session, + "group": group, + **slopes.iloc[i].to_dict(), + "first": first[i], + "last": last[i], + "change": last[i] - first[i], + "ratio": ratio[i], + } + ) + if math.isnan(slopes.iloc[i]["slope"]): + dropped.append(Dropped(session, f"fewer than 3 bins with {metric} (no slope)")) + has_slope = slopes["slope"].notna().to_numpy() + slope_test = run_one_sample( + slopes["slope"].to_numpy()[has_slope], 0.0, f"{metric} slope", name=group + ) + ends = np.isfinite(first) & np.isfinite(last) + habituation = run_paired(first[ends], last[ends], (bins[0], bins[-1]), metric, dropped=[]) + full = np.isfinite(matrix).all(axis=1) + trend = GroupTrend(group, slope_test, habituation) + if full.sum() >= 2: + trend.contrasts = trend_contrasts(matrix[full], times, correction=correction) + if len(bins) >= 3: + trend.anova = run_repeated(matrix[full], list(bins), metric, correction=correction) + complete[group] = matrix[full] + groups[group] = trend + frame = pd.DataFrame(per_animal) + result = TrendResult( + metric=metric, + bins=list(bins), + times=times, + time_unit=time_unit, + per_animal=frame, + groups=groups, + warnings=warnings, + dropped=dropped, + ) + if len(groups) >= 2: + slopes_by = { + g: frame.loc[(frame["group"] == g) & frame["slope"].notna(), "slope"].to_numpy() + for g in groups + } + result.slope_comparison = run_between(slopes_by, f"{metric} slope", correction=correction) + change_by = { + g: frame.loc[(frame["group"] == g) & frame["change"].notna(), "change"].to_numpy() + for g in groups + } + result.change_comparison = run_between(change_by, f"{metric} change", correction=correction) + usable = {g: m for g, m in complete.items() if m.shape[0] >= 2} + if len(usable) >= 2: + result.mixed = run_mixed( + usable, list(bins), metric, correction=correction, factor="bin" + ) + else: + warnings.append("too few complete animals per group for the group x bin ANOVA") + return result + + +def _group_of(value) -> str: + return UNGROUPED if value is None or str(value) in ("", "nan") else str(value) + + +def compare_trend( + tidy: pd.DataFrame, + metric: str, + *, + include_partial: bool = False, + correction: str = "holm", +) -> TrendResult: + """Trend of ``metric`` from a tidy time-bin frame. + + Bins are ordered and placed by ``bin_start_s`` / ``bin_end_s`` when the + frame has them (the time-bin export), in minutes; otherwise ``range`` (or + ``bin``) labels are taken in order of first appearance, one bin apart. + """ + label = "bin" if "bin" in tidy.columns else "range" + timed = {"bin_start_s", "bin_end_s"} <= set(tidy.columns) + if timed: + edges = tidy.groupby(label, sort=False)[["bin_start_s", "bin_end_s"]].first() + edges = edges.sort_values("bin_start_s") + widths = edges["bin_end_s"] - edges["bin_start_s"] + partial = widths < _PARTIAL * widths.max() + excluded = [] if include_partial else [str(b) for b in edges.index[partial]] + edges = edges[~edges.index.isin(excluded)] + bins = [str(b) for b in edges.index] + times = ((edges["bin_start_s"] + edges["bin_end_s"]) / 2 / 60.0).tolist() + unit = "min" + else: + bins = [str(b) for b in pd.unique(tidy[label])] + times, unit, excluded = list(range(len(bins))), "bin", [] + rows = tidy[tidy[label].astype(str).isin(bins)] + cells: dict[str, tuple[list, list[str]]] = {} + for session, block in rows.groupby("session", sort=False): + group = _group_of(block["group"].iloc[0] if "group" in block else None) + by_bin = dict(zip(block[label].astype(str), pd.to_numeric(block[metric]), strict=True)) + values = [float(by_bin.get(b, np.nan)) for b in bins] + matrix, sessions = cells.setdefault(group, ([], [])) + matrix.append(values) + sessions.append(str(session)) + order = sorted(cells, key=lambda g: (g == UNGROUPED, g)) + result = run_trend( + {g: (np.asarray(cells[g][0], dtype=float), cells[g][1]) for g in order}, + bins, + times, + metric, + time_unit=unit, + correction=correction, + ) + result.excluded_bins = excluded + if excluded: + result.warnings.insert( + 0, + f"left out partial bin {', '.join(excluded)}: it holds less time, so a total " + "would read as a decline", + ) + return result diff --git a/tests/unit/analysis/test_stats.py b/tests/unit/analysis/test_stats.py new file mode 100644 index 00000000..3f3fa2db --- /dev/null +++ b/tests/unit/analysis/test_stats.py @@ -0,0 +1,353 @@ +"""The statistics engine, pinned to independent reference values. + +The reference numbers below were produced once by pingouin 0.7.0 and +scikit-posthocs 0.17.1 (scipy 1.18.1) on the datasets beside them, in a +throwaway environment. GLIDER does not depend on either library; these +constants are what keeps its own arithmetic honest. +""" + +from __future__ import annotations + +import math + +import numpy as np +import pandas as pd +import pytest + +from glider.analysis.stats import ( + UNGROUPED, + adjust, + between_cells, + compare_epochs, + compare_groups, + dunn, + games_howell, + run_between, + run_paired, + tukey, + within_matrix, +) +from glider.analysis.stats.result import format_p + +# Unequal n, unequal spread, rounded so ranks tie. +A = np.array([9.8, 9.4, 9.8, 10.2, 9.0, 11.1, 10.7, 10.1, 10.6]) +B = np.array([9.5, 10.2, 12.8, 9.4, 13.3, 13.6, 13.1]) +C = np.array([14.5, 21.6, 13.4, 12.7, 13.3, 13.2, 15.0, 10.5, 7.1, 10.3, 16.0]) +# The same ten animals, before and during stimulation. +BASE = np.array([18.0, 22.0, 24.2, 15.9, 23.7, 21.3, 21.1, 19.9, 24.3, 17.2]) +STIM = np.array([21.1, 22.9, 23.6, 17.2, 26.5, 25.8, 26.6, 22.1, 24.5, 17.8]) + +TOL = {"rel": 1e-9, "abs": 1e-9} + + +def _pair(frame: pd.DataFrame, a: str, b: str) -> pd.Series: + return frame[(frame["a"] == a) & (frame["b"] == b)].iloc[0] + + +class TestTwoGroups: + def test_welch(self): + r = run_between({"A": A, "B": B}, "x", "welch") + assert r.statistic == pytest.approx(-2.1568514551531504, **TOL) + assert r.df[0] == pytest.approx(7.132295867688487, **TOL) + assert r.p == pytest.approx(0.06720500593969006, **TOL) + + def test_student(self): + r = run_between({"A": A, "B": B}, "x", "student") + assert r.statistic == pytest.approx(-2.3985588414152024, **TOL) + assert r.df == (14.0,) + + def test_hedges_g_and_its_interval(self): + r = run_between({"A": A, "B": B}, "x", "welch") + assert r.effect_name == "Hedges' g" + assert r.effect == pytest.approx(-1.142827672262066, **TOL) + assert r.effect_ci == pytest.approx((-2.307315454348, 0.021660109824), abs=1e-9) + + def test_mann_whitney_and_rank_biserial(self): + r = run_between({"A": A, "B": B}, "x", "mannwhitney") + assert r.statistic == 18.0 + assert r.p == pytest.approx(0.16786039631458505, **TOL) + assert r.effect == pytest.approx(-0.4285714285714286, **TOL) + assert r.df == () + + def test_the_difference_and_its_bootstrap_interval(self): + r = run_between({"A": A, "B": B}, "x", "welch") + assert r.difference == pytest.approx(A.mean() - B.mean()) + low, high = r.difference_ci + assert low < r.difference < high + + def test_permutation_is_reproducible(self): + first = run_between({"A": A, "B": B}, "x", "permutation") + again = run_between({"A": A, "B": B}, "x", "permutation") + assert first.p == again.p + assert first.difference_ci == again.difference_ci + # Close to Welch's answer on data this well behaved. + assert first.p == pytest.approx(0.067, abs=0.03) + + def test_a_two_group_test_refuses_three_groups(self): + r = run_between({"A": A, "B": B, "C": C}, "x", "welch") + assert not r.ok + assert "two groups" in r.warnings[0] + + +class TestSeveralGroups: + CELLS = {"A": A, "B": B, "C": C} + + def test_one_way_anova(self): + r = run_between(self.CELLS, "x", "anova") + assert r.statistic == pytest.approx(4.12870599851146, **TOL) + assert r.df == (2.0, 24.0) + assert r.p == pytest.approx(0.028772701951883733, **TOL) + assert r.effect == pytest.approx(0.25598494999490373, **TOL) + assert r.posthoc_name == "Tukey HSD" + + def test_welch_anova(self): + r = run_between(self.CELLS, "x", "welch_anova") + assert r.statistic == pytest.approx(5.974963959822446, **TOL) + assert r.df == pytest.approx((2.0, 11.342202939183386), **TOL) + assert r.p == pytest.approx(0.01689327405596058, **TOL) + assert r.posthoc_name == "Games-Howell" + + def test_kruskal_wallis(self): + r = run_between(self.CELLS, "x", "kruskal") + assert r.statistic == pytest.approx(7.613604610548366, **TOL) + assert r.p == pytest.approx(0.022219115428219315, **TOL) + assert r.effect == pytest.approx(7.613604610548366 / 26) + assert r.posthoc_name == "Dunn (Holm)" + + @pytest.mark.parametrize( + "a, b, diff, t, p", + [ + ("A", "B", -1.6222222222222236, -1.2419403543118823, 0.44085695205351405), + ("A", "C", -3.3404040404040405, -2.867356608864776, 0.022278181867078417), + ("B", "C", -1.718181818181817, -1.3710659270338121, 0.37143552962877135), + ], + ) + def test_tukey(self, a, b, diff, t, p): + row = _pair(tukey(self.CELLS), a, b) + assert row["diff"] == pytest.approx(diff, **TOL) + assert row["statistic"] == pytest.approx(t, **TOL) + # pingouin and scipy.stats.tukey_hsd agree on these to 1e-15. + assert row["p_adj"] == pytest.approx(p, rel=1e-7) + + @pytest.mark.parametrize( + "a, b, t, p", + [ + ("A", "B", -2.156851455153151, 0.146319810276771), + ("A", "C", -2.9467072855530065, 0.03341297723416825), + ("B", "C", -1.297469055214847, 0.4173043299973618), + ], + ) + def test_games_howell(self, a, b, t, p): + row = _pair(games_howell(self.CELLS), a, b) + assert row["statistic"] == pytest.approx(t, **TOL) + assert row["p_adj"] == pytest.approx(p, rel=1e-7) + + @pytest.mark.parametrize( + "a, b, p, p_holm", + [ + ("A", "B", 0.20035390989521307, 0.40070781979042613), + ("A", "C", 0.005815153961925041, 0.017445461885775122), + ("B", "C", 0.21899949966394483, 0.40070781979042613), + ], + ) + def test_dunn_is_tie_corrected_and_holm_adjusted(self, a, b, p, p_holm): + row = _pair(dunn(self.CELLS), a, b) + assert row["p"] == pytest.approx(p, **TOL) + assert row["p_adj"] == pytest.approx(p_holm, **TOL) + + def test_two_groups_get_no_posthoc(self): + assert run_between({"A": A, "B": B}, "x", "anova").posthoc is None + + def test_anova_on_two_groups_is_students_t_squared(self): + f = run_between({"A": A, "B": B}, "x", "anova") + t = run_between({"A": A, "B": B}, "x", "student") + assert f.statistic == pytest.approx(t.statistic**2) + assert f.p == pytest.approx(t.p) + + +class TestPaired: + def test_paired_t(self): + r = run_paired(BASE, STIM, ("Base", "Stim"), "x", "paired_t") + assert r.statistic == pytest.approx(-3.3340802494058734, **TOL) + assert r.df == (9.0,) + + def test_cohen_dz_and_its_interval(self): + r = run_paired(BASE, STIM, ("Base", "Stim"), "x", "paired_t") + assert r.effect == pytest.approx(-1.0543287489904811, **TOL) + assert r.effect_ci == pytest.approx((-1.946606971655, -0.162050526326), abs=1e-9) + + def test_wilcoxon_and_matched_rank_biserial(self): + r = run_paired(BASE, STIM, ("Base", "Stim"), "x", "wilcoxon") + assert r.statistic == 2.0 + assert r.p == pytest.approx(0.005859375, **TOL) + assert r.effect == pytest.approx(-0.9272727272727272, **TOL) + + def test_paired_permutation_is_reproducible(self): + first = run_paired(BASE, STIM, test="permutation") + assert first.p == run_paired(BASE, STIM, test="permutation").p + assert first.p < 0.05 + + def test_no_change_at_all_is_reported_not_raised(self): + r = run_paired(BASE, BASE, test="wilcoxon") + assert not r.ok and "zero" in r.warnings[0] + + def test_unequal_lengths_are_an_error(self): + with pytest.raises(ValueError): + run_paired([1, 2, 3], [1, 2]) + + +class TestAuto: + # Symmetric and unimodal: Shapiro-Wilk p = 0.999. (B above is bimodal, + # p = 0.044, so auto rightly sends it to a rank test.) + NORMAL = np.array([9.1, 9.8, 10.0, 10.3, 10.6, 10.9, 11.5, 9.6, 10.2]) + + def test_a_bimodal_group_is_not_called_normal(self): + assert run_between({"A": A, "B": B}, "x").test == "Mann-Whitney U" + + def test_normal_two_groups_get_welch(self): + r = run_between({"A": self.NORMAL, "B": self.NORMAL + 1.5}, "x") + assert r.test == "Welch's t-test" + assert "Shapiro-Wilk p ≥ .05" in r.reason + + def test_a_skewed_group_gets_mann_whitney_and_names_it(self): + skewed = np.array([1.0, 1.1, 1.0, 1.2, 1.1, 1.0, 9.0, 12.0, 30.0]) + r = run_between({"A": A, "Skewed": skewed}, "x") + assert r.test == "Mann-Whitney U" + assert "in Skewed" in r.reason + + def test_unequal_variances_get_welch_anova(self): + wide = np.array([-40.0, 30.0, -25.0, 45.0, 0.0, -10.0, 20.0, 35.0, -30.0, 5.0]) + r = run_between({"A": self.NORMAL, "B": self.NORMAL + 1.5, "Wide": wide}, "x") + assert r.test == "Welch's ANOVA" + assert "Levene" in r.reason + + def test_equal_variances_get_anova_and_tukey(self): + shifted = {"A": A, "A+1": A + 1, "A+2": A + 2} + r = run_between(shifted, "x") + assert r.test == "One-way ANOVA" + assert r.posthoc_name == "Tukey HSD" + + def test_normal_differences_get_the_paired_t(self): + assert run_paired(BASE, STIM).test == "Paired t-test" + + def test_skewed_differences_get_wilcoxon(self): + stim = BASE + np.array([0.1, 0.2, 0.1, 0.3, 0.2, 0.1, 0.2, 0.1, 9.0, 15.0]) + assert run_paired(BASE, stim).test == "Wilcoxon signed-rank" + + def test_small_groups_are_warned_about(self): + r = run_between({"A": A[:4], "B": B[:4]}, "x") + assert any("fewer than 5" in w for w in r.warnings) + + def test_a_chosen_test_carries_no_reason(self): + assert run_between({"A": A, "B": B}, "x", "student").reason == "" + + +class TestRefusals: + def test_one_group(self): + r = run_between({"A": A}, "x") + assert not r.ok and "two groups" in r.warnings[0] + + def test_a_group_of_one(self): + r = run_between({"A": A, "B": B[:1]}, "x") + assert not r.ok and "(B)" in r.warnings[0] + + def test_an_unknown_test(self): + with pytest.raises(ValueError): + run_between({"A": A, "B": B}, "x", "chi-by-eye") + + def test_constant_groups_give_nan_effect_not_an_exception(self): + r = run_between({"A": [1.0, 1.0, 1.0], "B": [1.0, 1.0, 1.0]}, "x", "student") + assert math.isnan(r.effect) + + +class TestCorrections: + P = [0.01, 0.04, 0.03, 0.005, 0.2, 0.04] + + def test_holm(self): + assert adjust(self.P, "holm") == pytest.approx([0.05, 0.12, 0.12, 0.03, 0.2, 0.12]) + + def test_bonferroni(self): + expected = [0.06, 0.24, 0.18, 0.03, 1.0, 0.24] + assert adjust(self.P, "bonferroni") == pytest.approx(expected) + + def test_benjamini_hochberg(self): + expected = [0.03, 0.048, 0.048, 0.03, 0.2, 0.048] + assert adjust(self.P, "fdr_bh") == pytest.approx(expected) + + def test_nan_is_kept_out_of_the_family(self): + out = adjust([0.01, float("nan"), 0.04], "bonferroni") + assert out[0] == pytest.approx(0.02) and math.isnan(out[1]) + + def test_an_unknown_method(self): + with pytest.raises(ValueError): + adjust([0.1], "sidak-ish") + + +def _tidy() -> pd.DataFrame: + """Two groups, two epochs; m4 has no Stim value and m6 no group.""" + rows = [] + values = { + "m1": ("ChR2", 10.0, 14.0), + "m2": ("ChR2", 11.0, 16.0), + "m3": ("ChR2", 9.5, 13.0), + "m4": ("ChR2", 10.5, None), + "m5": ("eYFP", 10.2, 10.1), + "m6": ("", 9.9, 10.4), + } + for session, (group, base, stim) in values.items(): + rows.append({"session": session, "group": group, "range": "Base", "dist": base}) + rows.append({"session": session, "group": group, "range": "Stim", "dist": stim}) + return pd.DataFrame(rows) + + +class TestFromTheEpochTable: + def test_between_cells_split_by_group_with_ungrouped_last(self): + cells, sessions, dropped = between_cells(_tidy(), "dist", "Stim") + assert list(cells) == ["ChR2", "eYFP", UNGROUPED] + assert sessions["ChR2"] == ["m1", "m2", "m3"] + assert [d.session for d in dropped] == ["m4"] + assert "Stim" in dropped[0].reason + + def test_within_drops_incomplete_animals_listwise(self): + matrix, kept, dropped = within_matrix(_tidy(), "dist", ["Base", "Stim"]) + assert kept == ["m1", "m2", "m3", "m5", "m6"] + assert matrix.shape == (5, 2) + assert [d.session for d in dropped] == ["m4"] + + def test_within_one_group(self): + _matrix, kept, _ = within_matrix(_tidy(), "dist", ["Base", "Stim"], group="ChR2") + assert kept == ["m1", "m2", "m3"] + + def test_compare_groups_carries_who_was_dropped(self): + r = compare_groups(_tidy(), "dist", "Base", "welch") + assert r.n == {"ChR2": 4, "eYFP": 1, UNGROUPED: 1} + assert not r.ok # two-group test, three groups + + def test_compare_epochs(self): + r = compare_epochs(_tidy(), "dist", ["Base", "Stim"], "paired_t", group="ChR2") + assert r.n == {"Base": 3, "Stim": 3} + assert [d.session for d in r.dropped] == ["m4"] + assert r.difference < 0 + + def test_compare_epochs_wants_two(self): + with pytest.raises(ValueError): + compare_epochs(_tidy(), "dist", ["Base"]) + + +class TestReporting: + def test_methods_sentence(self): + text = run_between({"A": A, "B": B}, "distance", "welch").methods() + assert text.startswith("Welch's t-test on distance, A (n = 9), B (n = 7): t(7.13) = -2.16") + assert "p = .067" in text and "Hedges' g = -1.14 [95% CI -2.31, 0.02]" in text + + def test_integer_df_print_without_decimals(self): + assert run_between({"A": A, "B": B}, "x", "student").statistic_text().startswith("t(14)") + + @pytest.mark.parametrize( + "p, text", [(0.0004, "p < .001"), (0.04567, "p = .046"), (float("nan"), "p = n/a")] + ) + def test_format_p(self, p, text): + assert format_p(p) == text + + def test_a_failed_test_says_why(self): + assert "not run" in run_between({"A": A}, "x").methods() diff --git a/tests/unit/analysis/test_stats_repeated.py b/tests/unit/analysis/test_stats_repeated.py new file mode 100644 index 00000000..29ed625f --- /dev/null +++ b/tests/unit/analysis/test_stats_repeated.py @@ -0,0 +1,407 @@ +"""Repeated-measures, Friedman and mixed ANOVA, pinned to independent references. + +Reference values: pingouin 0.7.0 (``rm_anova``, ``sphericity``, ``epsilon``, +``friedman``, ``pairwise_tests``, ``mixed_anova``), computed once in a throwaway +environment on the data below. pingouin follows R's ``mauchly.test`` and takes +a mixed design's sphericity from the pooled within-group covariance, as SPSS +and afex do. +""" + +from __future__ import annotations + +import math + +import numpy as np +import pandas as pd +import pytest + +from glider.analysis.stats import ( + compare_epochs, + compare_mixed, + run_mixed, + run_repeated, + sphericity, +) + +# Twelve animals, four epochs; e3 swings with each animal: sphericity fails. +SPH = np.array( + [ + [20.06, 21.08, 22.98, 19.44], + [20.13, 21.5, 24.13, 31.79], + [19.44, 19.19, 22.1, 23.4], + [18.57, 18.95, 21.15, 11.86], + [18.42, 20.07, 22.16, 19.57], + [17.79, 19.07, 21.05, 18.36], + [19.17, 20.36, 22.39, 21.61], + [22.04, 23.44, 25.73, 32.77], + [18.09, 19.53, 22.83, 18.78], + [18.64, 19.35, 20.83, 20.95], + [20.35, 22.51, 24.5, 30.69], + [20.85, 21.31, 23.79, 21.44], + ] +) +# Ten animals, three epochs, an additive model: sphericity holds. +OK = np.array( + [ + [11.3, 10.7, 12.1], + [6.6, 8.8, 10.8], + [9.9, 11.2, 11.5], + [6.3, 6.8, 9.0], + [9.5, 8.6, 9.8], + [10.1, 8.8, 11.5], + [12.4, 13.2, 12.9], + [13.1, 13.5, 15.3], + [5.8, 7.7, 5.8], + [9.1, 8.8, 7.1], + ] +) +# Mixed: 7 controls and 9 ChR2 animals; ChR2 responds in e1. +CTL = np.array( + [ + [9.88, 9.71, 9.39], + [8.9, 9.22, 7.37], + [10.25, 10.86, 9.47], + [14.75, 14.13, 15.66], + [8.53, 8.93, 7.95], + [8.63, 7.25, 7.82], + [10.77, 8.78, 11.46], + ] +) +CHR2 = np.array( + [ + [6.23, 10.03, 6.53], + [12.61, 13.88, 12.46], + [7.52, 10.43, 8.03], + [12.82, 15.1, 13.52], + [10.28, 12.18, 10.94], + [6.37, 10.95, 7.55], + [12.2, 12.36, 12.69], + [14.54, 14.87, 12.83], + [10.21, 18.09, 9.41], + ] +) +# Mixed: three groups (6, 8, 5) over four epochs, nothing going on. +HX = np.array( + [ + [3.93, 5.14, 4.15, 4.54], + [2.94, 2.84, 5.39, 3.54], + [2.65, 2.32, 1.92, 1.85], + [6.73, 5.8, 5.95, 6.12], + [6.14, 5.65, 5.35, 4.4], + [3.73, 6.06, 6.07, 4.97], + ] +) +HY = np.array( + [ + [6.31, 6.69, 7.69, 7.3], + [5.03, 4.48, 6.99, 5.16], + [6.23, 7.09, 4.6, 4.14], + [6.64, 6.43, 6.14, 6.42], + [3.46, 2.8, 4.15, 3.34], + [5.63, 7.78, 7.21, 7.68], + [2.14, 2.47, 2.13, 2.24], + [6.49, 5.51, 6.65, 6.63], + ] +) +HZ = np.array( + [ + [6.59, 7.58, 8.78, 7.96], + [2.99, 2.62, 4.07, 2.89], + [4.13, 5.64, 4.36, 3.08], + [4.34, 6.2, 4.91, 3.69], + [4.04, 6.11, 5.14, 5.03], + ] +) + +RM = { + "SPH": { + "F": 4.569840299068531, + "p_unc": 0.008755902565414194, + "p_gg": 0.054035496642435794, + "np2": 0.2935059198610998, + "eps_gg": 0.34494665198442204, + "eps_hf": 0.3484834304513767, + "W": 0.0016846191600016023, + "p_spher": 7.112047303375329e-12, + "chi2": 62.08820876345587, + "Q": 20.7, + "kendall_w": 0.575, + "friedman_p": 0.00012151040342188821, + "pairs": [ + ("e0", "e1", -5.692842658277161, 0.0005586444587153528), + ("e0", "e2", -15.59999478861659, 3.768342256776994e-08), + ("e0", "e3", -1.9921394569627595, 0.21529057425689085), + ("e1", "e2", -16.005725886292534, 3.4470823213255915e-08), + ("e1", "e3", -1.3738547477116823, 0.39366993783402343), + ("e2", "e3", 0.16965181282896502, 0.8683637375889692), + ], + }, + "OK": { + "F": 3.0979202491807514, + "p_unc": 0.06978886407996222, + "p_gg": 0.07939686506940949, + "np2": 0.25607048032826446, + "eps_gg": 0.8656375137331306, + "eps_hf": 1.0, + "W": 0.8447820431355607, + "p_spher": 0.5093058860319419, + "chi2": 1.3494129761131404, + "Q": 5.282051282051282, + "kendall_w": 0.2641025641025641, + "friedman_p": 0.071288116031728, + "pairs": [ + ("e0", "e1", -1.0715994395536397, 0.3325584400764982), + ("e0", "e2", -2.1945645403483645, 0.16749969478788607), + ("e1", "e2", -1.5062050138151357, 0.3325584400764982), + ], + }, +} + +MIXED_TWO = [ + { + "F": 1.1260673668046084, + "df1": 1, + "df2": 14, + "p_unc": 0.30657974740684835, + "np2": 0.07444548139960393, + "p_gg": None, + "eps": None, + "W_spher": None, + "p_spher": None, + }, + { + "F": 7.078723785962117, + "df1": 2, + "df2": 28, + "p_unc": 0.00325066427677455, + "np2": 0.33582316737203816, + "p_gg": 0.012110729965951417, + "eps": 0.6223684443608658, + "W_spher": 0.393234732479185, + "p_spher": 0.0023186604520896107, + }, + { + "F": 7.452743966506648, + "df1": 2, + "df2": 28, + "p_unc": 0.0025411634797635734, + "np2": 0.3474028300595175, + "p_gg": 0.010286558313837454, + "eps": 0.6223684443608658, + "W_spher": 0.393234732479185, + "p_spher": 0.0023186604520896107, + }, +] +MIXED_THREE = [ + { + "F": 0.5454978458663371, + "df1": 2, + "df2": 16, + "p_unc": 0.5899577948536606, + "np2": 0.06383453084950545, + "p_gg": None, + "eps": None, + "W_spher": None, + "p_spher": None, + }, + { + "F": 2.704679796264302, + "df1": 3, + "df2": 48, + "p_unc": 0.055704993628794615, + "np2": 0.14459909636114063, + "p_gg": 0.07529725340002384, + "eps": 0.7413970001695399, + "W_spher": 0.543073140390134, + "p_spher": 0.11025454219636883, + }, + { + "F": 0.6625728763206146, + "df1": 6, + "df2": 48, + "p_unc": 0.6800145858156408, + "np2": 0.07648684585751377, + "p_gg": 0.6373898140778078, + "eps": 0.7413970001695399, + "W_spher": 0.543073140390134, + "p_spher": 0.11025454219636883, + }, +] + +TOL = {"rel": 1e-9, "abs": 1e-12} + + +def _epochs(k: int) -> list[str]: + return [f"e{j}" for j in range(k)] + + +@pytest.mark.parametrize("key, data", [("SPH", SPH), ("OK", OK)]) +class TestRepeatedMeasuresAnova: + def test_f_p_and_partial_eta(self, key, data): + ref = RM[key] + r = run_repeated(data, _epochs(data.shape[1]), "x", "rm_anova") + assert r.statistic == pytest.approx(ref["F"], **TOL) + assert r.details["p uncorrected"] == pytest.approx(ref["p_unc"], **TOL) + assert r.details["p GG"] == pytest.approx(ref["p_gg"], **TOL) + assert r.effect == pytest.approx(ref["np2"], **TOL) + + def test_sphericity(self, key, data): + ref = RM[key] + s = sphericity(data - data.mean(axis=0), data.shape[0] - 1) + assert s.gg == pytest.approx(ref["eps_gg"], **TOL) + assert s.hf == pytest.approx(ref["eps_hf"], **TOL) + assert s.w == pytest.approx(ref["W"], **TOL) + assert s.chi2 == pytest.approx(ref["chi2"], **TOL) + assert s.p == pytest.approx(ref["p_spher"], rel=1e-9, abs=1e-20) + + def test_pairwise_follow_ups_are_holm_adjusted(self, key, data): + r = run_repeated(data, _epochs(data.shape[1]), "x", "rm_anova") + for (a, b, t, p_holm), (_, row) in zip(RM[key]["pairs"], r.posthoc.iterrows(), strict=True): + assert (row["a"], row["b"]) == (a, b) + assert row["statistic"] == pytest.approx(t, **TOL) + assert row["p_adj"] == pytest.approx(p_holm, **TOL) + + def test_friedman_and_kendalls_w(self, key, data): + ref = RM[key] + r = run_repeated(data, _epochs(data.shape[1]), "x", "friedman") + assert r.statistic == pytest.approx(ref["Q"], **TOL) + assert r.effect == pytest.approx(ref["kendall_w"], **TOL) + assert r.p == pytest.approx(ref["friedman_p"], **TOL) + assert r.posthoc_name.startswith("Wilcoxon") + + +class TestSphericityDecisions: + def test_a_violation_reports_the_corrected_p_and_df(self): + r = run_repeated(SPH, _epochs(4), "x", "rm_anova") + assert r.p == pytest.approx(RM["SPH"]["p_gg"], **TOL) + assert r.df == pytest.approx((3 * RM["SPH"]["eps_gg"], 33 * RM["SPH"]["eps_gg"])) + assert "Greenhouse-Geisser" in r.reason + # The correction is what turns p = .009 into p = .054 here. + assert r.details["p uncorrected"] < 0.01 < 0.05 < r.p + + def test_no_violation_reports_the_plain_test(self): + r = run_repeated(OK, _epochs(3), "x", "rm_anova") + assert r.p == pytest.approx(RM["OK"]["p_unc"], **TOL) + assert r.df == (2.0, 18.0) + assert "Greenhouse" not in r.reason + + def test_two_levels_are_always_spherical(self): + s = sphericity(OK[:, :2] - OK[:, :2].mean(axis=0), 9) + assert (s.w, s.gg, s.hf, s.p) == (1.0, 1.0, 1.0, 1.0) + + def test_fewer_animals_than_contrasts_cannot_test_it(self): + few = SPH[:3] + s = sphericity(few - few.mean(axis=0), 2) + assert math.isnan(s.p) and not s.violated + r = run_repeated(few, _epochs(4), "x", "rm_anova") + assert any("sphericity" in w for w in r.warnings) + + +class TestRepeatedAuto: + def test_normal_residuals_get_the_anova(self): + r = run_repeated(OK, _epochs(3), "x") + assert r.test == "Repeated-measures ANOVA" + assert "residuals" in r.reason + + def test_skewed_residuals_get_friedman(self): + skewed = OK.copy() + skewed[0, 2] += 40.0 + r = run_repeated(skewed, _epochs(3), "x") + assert r.test == "Friedman test" + + def test_two_epochs_belong_to_the_paired_tests(self): + with pytest.raises(ValueError): + run_repeated(OK[:, :2], _epochs(2)) + + +def _check_mixed(result, reference): + for (name, r), ref in zip(result.effects.items(), reference, strict=True): + assert r.statistic == pytest.approx(ref["F"], **TOL), name + assert r.details["p uncorrected"] == pytest.approx(ref["p_unc"], **TOL), name + assert r.effect == pytest.approx(ref["np2"], **TOL), name + if ref["eps"] is not None: + assert r.details["p GG"] == pytest.approx(ref["p_gg"], **TOL), name + assert result.sphericity.gg == pytest.approx(ref["eps"], **TOL) + assert result.sphericity.w == pytest.approx(ref["W_spher"], **TOL) + assert result.sphericity.p == pytest.approx(ref["p_spher"], **TOL) + + +class TestMixedAnova: + def test_two_unbalanced_groups(self): + result = run_mixed({"ctl": CTL, "chr2": CHR2}, _epochs(3), "x") + _check_mixed(result, MIXED_TWO) + assert [r.df[:2] for r in result.effects.values()][0] == (1.0, 14.0) + + def test_three_groups(self): + result = run_mixed({"x": HX, "y": HY, "z": HZ}, _epochs(4), "v") + _check_mixed(result, MIXED_THREE) + + def test_a_violation_corrects_the_within_effects_only(self): + result = run_mixed({"ctl": CTL, "chr2": CHR2}, _epochs(3), "x") + assert result.sphericity.violated + assert result.effects["group"].p == pytest.approx(MIXED_TWO[0]["p_unc"], **TOL) + assert result.effects["epoch"].p == pytest.approx(MIXED_TWO[1]["p_gg"], **TOL) + assert result.effects["group × epoch"].p == pytest.approx(MIXED_TWO[2]["p_gg"], **TOL) + assert "Greenhouse-Geisser corrected" in result.methods() + + def test_a_significant_interaction_brings_simple_effects(self): + result = run_mixed({"ctl": CTL, "chr2": CHR2}, _epochs(3), "x") + simple = result.simple_effects + assert list(simple["effect"]) == ["epoch", "epoch", "group", "group", "group"] + assert list(simple["within"]) == ["ctl", "chr2", "e0", "e1", "e2"] + # Only ChR2 changes across epochs. + by = simple.set_index("within") + assert by.loc["chr2", "p_adj"] < 0.05 < by.loc["ctl", "p_adj"] + assert (simple["p_adj"] >= simple["p"]).all() + + def test_no_interaction_no_simple_effects(self): + result = run_mixed({"x": HX, "y": HY, "z": HZ}, _epochs(4), "v") + assert result.simple_effects is None + + def test_the_table(self): + table = run_mixed({"ctl": CTL, "chr2": CHR2}, _epochs(3), "x").table() + assert list(table["effect"]) == ["group", "epoch", "group × epoch"] + assert table["p_uncorrected"].tolist() == pytest.approx([r["p_unc"] for r in MIXED_TWO]) + + def test_non_normal_residuals_are_warned_about(self): + result = run_mixed({"ctl": CTL, "chr2": CHR2}, _epochs(3), "x") + assert any("not normal" in w for w in result.warnings) + + def test_one_group_is_refused(self): + with pytest.raises(ValueError): + run_mixed({"ctl": CTL}, _epochs(3)) + + +def _tidy() -> pd.DataFrame: + """The two-group mixed data as the epoch table writes it, m03 missing e2.""" + rows = [] + for group, data in (("ctl", CTL), ("chr2", CHR2)): + for i, animal in enumerate(data): + for j, value in enumerate(animal): + session = f"{group}{i:02d}" + missing = session == "ctl03" and j == 2 + rows.append( + { + "session": session, + "group": group, + "range": f"e{j}", + "dist": float("nan") if missing else value, + } + ) + return pd.DataFrame(rows) + + +class TestFromTheEpochTable: + def test_compare_mixed_matches_the_matrices_and_names_who_was_dropped(self): + result = compare_mixed(_tidy(), "dist", _epochs(3)) + assert result.n == {"chr2": 9, "ctl": 6} + assert [d.session for d in result.dropped] == ["ctl03"] + expected = run_mixed({"chr2": CHR2, "ctl": np.delete(CTL, 3, axis=0)}, _epochs(3), "dist") + assert result.effects["group × epoch"].statistic == pytest.approx( + expected.effects["group × epoch"].statistic + ) + + def test_compare_epochs_with_three_is_repeated_measures(self): + r = compare_epochs(_tidy(), "dist", _epochs(3), "rm_anova", group="chr2") + assert r.test == "Repeated-measures ANOVA" + assert r.n == {"e0": 9, "e1": 9, "e2": 9} diff --git a/tests/unit/analysis/test_stats_trend.py b/tests/unit/analysis/test_stats_trend.py new file mode 100644 index 00000000..ea6f7546 --- /dev/null +++ b/tests/unit/analysis/test_stats_trend.py @@ -0,0 +1,227 @@ +"""Trends across time bins: slopes, polynomial contrasts, first vs last. + +No outside library computes repeated-measures polynomial contrasts the way +these do, so they are pinned by identities that must hold exactly: on +equally spaced bins the linear contrast *is* the slope test, and with three +bins the two contrasts together *are* the repeated-measures ANOVA. +""" + +from __future__ import annotations + +import math + +import numpy as np +import pandas as pd +import pytest +from scipy import stats + +from glider.analysis.stats import ( + animal_slopes, + compare_trend, + polynomial_contrasts, + run_one_sample, + run_paired, + run_repeated, + run_trend, + trend_contrasts, +) + +RNG = np.random.default_rng(11) +# Twelve animals over six one-minute bins, declining about 1.5 a minute. +DECLINE = (50 + RNG.normal(0, 4, (12, 1)) - 1.5 * np.arange(6) + RNG.normal(0, 1.5, (12, 6))).round( + 2 +) +MINUTES = np.arange(6) + 0.5 + + +class TestOneSample: + VALUES = np.array([0.4, -0.2, 1.1, 0.8, 0.5, 0.9, 0.3, 1.4, -0.1, 0.6]) + + def test_t_matches_scipy(self): + r = run_one_sample(self.VALUES, 0.0, "x", "t") + ref = stats.ttest_1samp(self.VALUES, 0.0) + assert r.statistic == pytest.approx(ref.statistic) + assert r.p == pytest.approx(ref.pvalue) + assert r.effect == pytest.approx(self.VALUES.mean() / self.VALUES.std(ddof=1)) + assert r.test == "One-sample t-test" and r.effect_name == "Cohen's d" + + def test_mu_shifts_the_comparison(self): + r = run_one_sample(self.VALUES, 1.0, "ratio", "t") + assert r.statistic == pytest.approx(stats.ttest_1samp(self.VALUES, 1.0).statistic) + assert r.difference == pytest.approx(self.VALUES.mean() - 1.0) + + def test_wilcoxon_matches_scipy(self): + r = run_one_sample(self.VALUES, 0.0, "x", "wilcoxon") + assert r.p == pytest.approx(stats.wilcoxon(self.VALUES).pvalue) + + def test_its_mean_ci_brackets_the_mean(self): + low, high = run_one_sample(self.VALUES).difference_ci + assert low < self.VALUES.mean() < high + + def test_paired_is_one_sample_on_the_differences(self): + a, b = self.VALUES + 3, self.VALUES[::-1] + paired = run_paired(a, b, test="paired_t") + single = run_one_sample(a - b, test="t") + assert paired.statistic == pytest.approx(single.statistic) + assert paired.p == pytest.approx(single.p) + + +class TestContrastWeights: + def test_orthonormal_zero_sum_and_signed_to_rise(self): + w = polynomial_contrasts(MINUTES) + assert w.shape == (6, 2) + assert w.T @ w == pytest.approx(np.eye(2), abs=1e-12) + assert w.sum(axis=0) == pytest.approx([0, 0], abs=1e-12) + assert w[-1, 0] > w[0, 0] # linear rises with time + assert w[0, 1] > w[2, 1] # quadratic is U-shaped + + def test_unequal_spacing_is_respected(self): + times = np.array([0.5, 1.5, 2.5, 3.25]) # a short last bin + linear = polynomial_contrasts(times)[:, 0] + centred = times - times.mean() + assert linear == pytest.approx(centred / np.linalg.norm(centred)) + + def test_two_bins_have_only_a_linear_contrast(self): + assert polynomial_contrasts([0.5, 1.5]).shape == (2, 1) + + +class TestIdentities: + def test_slopes_match_linregress(self): + slopes = animal_slopes(DECLINE, MINUTES) + for row, values in zip(slopes.itertuples(), DECLINE, strict=True): + fit = stats.linregress(MINUTES, values) + assert row.slope == pytest.approx(fit.slope) + assert row.r2 == pytest.approx(fit.rvalue**2) + + def test_equal_spacing_makes_the_linear_contrast_the_slope_test(self): + linear = trend_contrasts(DECLINE, MINUTES).iloc[0] + slopes = animal_slopes(DECLINE, MINUTES)["slope"] + slope_t = stats.ttest_1samp(slopes, 0.0) + assert linear["t"] == pytest.approx(slope_t.statistic) + assert linear["p"] == pytest.approx(slope_t.pvalue) + + def test_with_three_bins_the_contrasts_sum_to_the_anova(self): + three = DECLINE[:, :3] + contrasts = trend_contrasts(three, MINUTES[:3]) + n = three.shape[0] + ss_contrasts = n * (contrasts["estimate"] ** 2).sum() + grand = three.mean() + ss_bins = n * ((three.mean(axis=0) - grand) ** 2).sum() + assert ss_contrasts == pytest.approx(ss_bins) + # And partial eta squared of the linear trend sits under 1. + assert 0 < contrasts.iloc[0]["partial_eta2"] < 1 + + def test_contrasts_are_holm_adjusted_together(self): + c = trend_contrasts(DECLINE, MINUTES) + assert (c["p_adj"] >= c["p"]).all() + + def test_a_missing_bin_still_gets_a_slope(self): + holes = DECLINE.copy() + holes[0, 2] = np.nan + slopes = animal_slopes(holes, MINUTES) + fit = stats.linregress(np.delete(MINUTES, 2), np.delete(holes[0], 2)) + assert slopes.iloc[0]["slope"] == pytest.approx(fit.slope) + assert slopes.iloc[0]["n_bins"] == 5 + + def test_under_three_bins_there_is_no_slope(self): + sparse = np.full((1, 6), np.nan) + sparse[0, :2] = [1.0, 2.0] + assert math.isnan(animal_slopes(sparse, MINUTES).iloc[0]["slope"]) + + +def _tidy(groups: dict[str, np.ndarray], *, last_end: float | None = None) -> pd.DataFrame: + """As the time-bin export writes it: one row per animal per one-minute bin.""" + rows = [] + for group, matrix in groups.items(): + for i, animal in enumerate(matrix): + for b, value in enumerate(animal): + start = 60.0 * b + end = 60.0 * (b + 1) + if b == len(animal) - 1 and last_end is not None: + end = last_end + rows.append( + { + "session": f"{group}{i:02d}", + "group": group, + "bin": f"{b}–{b + 1} min", + "bin_start_s": start, + "bin_end_s": end, + "dist": value, + } + ) + return pd.DataFrame(rows) + + +class TestCompareTrend: + def test_one_group(self): + r = compare_trend(_tidy({"all": DECLINE}), "dist") + trend = r.groups["all"] + assert r.time_unit == "min" + assert r.times == pytest.approx(MINUTES.tolist()) + assert trend.slope_test.difference == pytest.approx( + animal_slopes(DECLINE, MINUTES)["slope"].mean() + ) + assert trend.slope_test.p < 0.001 + assert trend.anova.test == "Repeated-measures ANOVA" + assert r.slope_comparison is None and r.mixed is None + + def test_first_vs_last_is_a_paired_test_of_the_end_bins(self): + r = compare_trend(_tidy({"all": DECLINE}), "dist") + expected = run_paired(DECLINE[:, 0], DECLINE[:, -1], test="auto") + assert r.groups["all"].habituation.statistic == pytest.approx(expected.statistic) + animal = r.per_animal.iloc[0] + assert animal["change"] == pytest.approx(DECLINE[0, -1] - DECLINE[0, 0]) + assert animal["ratio"] == pytest.approx(DECLINE[0, -1] / DECLINE[0, 0]) + + def test_groups_are_compared_and_crossed_with_bins(self): + flat = (50 + RNG.normal(0, 4, (10, 1)) + RNG.normal(0, 1.5, (10, 6))).round(2) + r = compare_trend(_tidy({"chr2": DECLINE, "ctl": flat}), "dist") + assert set(r.groups) == {"chr2", "ctl"} + assert r.slope_comparison.p < 0.001 + assert r.mixed is not None + assert list(r.mixed.effects) == ["group", "bin", "group × bin"] + assert r.mixed.effects["group × bin"].p < 0.001 + assert "per min" in r.methods() + + def test_a_partial_last_bin_is_left_out_by_default(self): + """A total over half a minute halves: it must not read as habituation.""" + steady = (50 + RNG.normal(0, 4, (12, 1)) + RNG.normal(0, 1.0, (12, 6))).round(2) + partial = steady.copy() + partial[:, -1] *= 0.5 # the last bin is 30 s of a steady animal + tidy = _tidy({"all": partial}, last_end=330.0) + + kept = compare_trend(tidy, "dist", include_partial=True) + assert kept.groups["all"].slope_test.p < 0.05 # the fake decline + + r = compare_trend(tidy, "dist") + assert r.excluded_bins == ["5–6 min"] + assert len(r.bins) == 5 + assert "partial bin" in r.warnings[0] + assert r.groups["all"].slope_test.p > 0.05 # steady, as it was + + def test_an_animal_with_too_few_bins_is_named(self): + tidy = _tidy({"all": DECLINE}) + short = tidy["session"].eq("all03") & tidy["bin_start_s"].ge(120) + tidy.loc[short, "dist"] = np.nan + r = compare_trend(tidy, "dist") + assert any(d.session == "all03" and "slope" in d.reason for d in r.dropped) + # The contrasts and ANOVA use complete animals only. + assert r.groups["all"].anova.n["0–1 min"] == 11 + + def test_without_bin_times_labels_are_taken_in_order(self): + tidy = _tidy({"all": DECLINE}).drop(columns=["bin_start_s", "bin_end_s"]) + tidy = tidy.rename(columns={"bin": "range"}) + r = compare_trend(tidy, "dist") + assert r.time_unit == "bin" + assert r.times == [0, 1, 2, 3, 4, 5] + assert r.excluded_bins == [] + + def test_one_bin_is_refused(self): + with pytest.raises(ValueError): + run_trend({"all": (DECLINE[:, :1], ["a"] * 12)}, ["0–1 min"], [0.5]) + + def test_the_anova_matches_running_it_directly(self): + r = compare_trend(_tidy({"all": DECLINE}), "dist") + bins = r.bins + direct = run_repeated(DECLINE, bins, "dist") + assert r.groups["all"].anova.statistic == pytest.approx(direct.statistic) From 19e75e278a92e8b4d9f6202d0579c53441df9534 Mon Sep 17 00:00:00 2001 From: gbradham Date: Thu, 29 Oct 2026 13:39:33 -0500 Subject: [PATCH 4/4] fix(review): exported CSVs open cleanly in Excel Bin labels use a plain hyphen (0-1 min), and the window, epoch and time-bin exports are written with a UTF-8 byte-order mark, so Excel no longer shows an en dash or other non-ASCII range names as mojibake. --- src/glider/analysis/epochs.py | 3 ++- src/glider/gui/behavior/analysis_window.py | 11 ++++++++--- tests/unit/analysis/test_epochs.py | 6 +++--- tests/unit/gui/behavior/test_analysis_window.py | 10 ++++++---- 4 files changed, 19 insertions(+), 11 deletions(-) diff --git a/src/glider/analysis/epochs.py b/src/glider/analysis/epochs.py index 284a9983..63cc9ba1 100644 --- a/src/glider/analysis/epochs.py +++ b/src/glider/analysis/epochs.py @@ -255,7 +255,8 @@ def time_bins(start_s: float, end_s: float, bin_s: float) -> list[Epoch]: k = first while k * bin_s < end_s - 1e-9: lo, hi = max(k * bin_s, start_s), min((k + 1) * bin_s, end_s) - name = f"{_bin_edge(lo, unit)}–{_bin_edge(hi, unit)} {suffix}" + # A plain hyphen: an en dash reads as "–" in a CSV opened by Excel. + name = f"{_bin_edge(lo, unit)}-{_bin_edge(hi, unit)} {suffix}" bins.append(Epoch(f"bin_{k - first + 1}", name, lo, hi)) k += 1 return bins diff --git a/src/glider/gui/behavior/analysis_window.py b/src/glider/gui/behavior/analysis_window.py index 8b833f32..3a9bc09d 100644 --- a/src/glider/gui/behavior/analysis_window.py +++ b/src/glider/gui/behavior/analysis_window.py @@ -134,6 +134,11 @@ def _short_path(path: Path, keep: int = 3) -> str: return str(path) if len(parts) <= keep else "…/" + "/".join(parts[-keep:]) +#: Exported tables carry a byte-order mark: without one Excel reads UTF-8 as +#: Windows-1252, and a range named "Stim–Post" opens as "Stim–Post". +_EXCEL_CSV = "utf-8-sig" + + def _unique(ids: list[str]) -> list[str]: """Session ids made unique, later repeats suffixed `` (2)``, `` (3)``. @@ -2908,7 +2913,7 @@ def _export_window(self) -> None: # session's empty cells already say it), are not measurements. frame = frame[[c for c in frame.columns if not str(c).startswith("_") and c != "outside"]] try: - frame.to_csv(path, index=False) + frame.to_csv(path, index=False, encoding=_EXCEL_CSV) except OSError as e: QMessageBox.critical(self, "Export window", f"Could not write {path}: {e}") return @@ -3084,7 +3089,7 @@ def _export_epochs(self, shape: str) -> None: if not path: return try: - frame.to_csv(path, index=False) + frame.to_csv(path, index=False, encoding=_EXCEL_CSV) except OSError as e: QMessageBox.critical(self, "Export epoch table", f"Could not write {path}: {e}") return @@ -3189,7 +3194,7 @@ def _export_bins(self) -> None: QMessageBox.information(self, "Export time bins", "The span holds no whole frame.") return try: - frame.to_csv(path, index=False) + frame.to_csv(path, index=False, encoding=_EXCEL_CSV) except OSError as e: QMessageBox.critical(self, "Export time bins", f"Could not write {path}: {e}") return diff --git a/tests/unit/analysis/test_epochs.py b/tests/unit/analysis/test_epochs.py index 8c6bc909..470bc6c0 100644 --- a/tests/unit/analysis/test_epochs.py +++ b/tests/unit/analysis/test_epochs.py @@ -182,17 +182,17 @@ def test_retrieval_metrics_join_the_catalog_only_when_present(): class TestTimeBins: def test_whole_minutes_are_named_in_minutes(self): bins = time_bins(0.0, 180.0, 60.0) - assert [b.name for b in bins] == ["0–1 min", "1–2 min", "2–3 min"] + assert [b.name for b in bins] == ["0-1 min", "1-2 min", "2-3 min"] assert [(b.start_s, b.end_s) for b in bins] == [(0, 60), (60, 120), (120, 180)] assert len({b.key for b in bins}) == 3 def test_the_last_bin_is_clipped_to_the_end(self): bins = time_bins(0.0, 150.0, 60.0) - assert bins[-1].name == "2–2.5 min" + assert bins[-1].name == "2-2.5 min" assert bins[-1].end_s == 150.0 def test_other_widths_are_named_in_seconds(self): - assert [b.name for b in time_bins(0.0, 60.0, 30.0)] == ["0–30 s", "30–60 s"] + assert [b.name for b in time_bins(0.0, 60.0, 30.0)] == ["0-30 s", "30-60 s"] def test_edges_stay_on_multiples_of_the_width(self): bins = time_bins(90.0, 200.0, 60.0) diff --git a/tests/unit/gui/behavior/test_analysis_window.py b/tests/unit/gui/behavior/test_analysis_window.py index fde970e8..ca2ceaac 100644 --- a/tests/unit/gui/behavior/test_analysis_window.py +++ b/tests/unit/gui/behavior/test_analysis_window.py @@ -2278,8 +2278,8 @@ def test_tidy_is_a_row_per_session_per_bin(self, qtbot, tmp_path): frame = win.bin_frame(0.0, 10.0, 5.0, "tidy") assert len(frame) == 4 assert list(frame.columns[:5]) == ["session", "group", "bin", "bin_start_s", "bin_end_s"] - first = frame[frame["bin"] == "0–5 s"] - second = frame[frame["bin"] == "5–10 s"] + first = frame[frame["bin"] == "0-5 s"] + second = frame[frame["bin"] == "5-10 s"] assert first["groom_s"].tolist() == pytest.approx([5.0, 5.0]) assert first["rear_s"].tolist() == pytest.approx([0.0, 0.0]) assert second["rear_s"].tolist() == pytest.approx([5.0, 5.0]) @@ -2288,7 +2288,7 @@ def test_wide_is_a_row_per_session(self, qtbot, tmp_path): win = self._win(qtbot, tmp_path) frame = win.bin_frame(0.0, 10.0, 5.0, "wide") assert len(frame) == 2 - assert "5–10 s__rear_s" in frame.columns + assert "5-10 s__rear_s" in frame.columns def test_the_whole_session_span_reaches_the_last_frame(self, qtbot, tmp_path): win = self._win(qtbot, tmp_path) @@ -2309,8 +2309,10 @@ def test_export_writes_the_bins(self, qtbot, tmp_path, monkeypatch): win._export_bins() written = pd.read_csv(out) # 10 s in 4 s bins: two whole and one clipped, for each of two sessions. - assert sorted(set(written["bin"])) == ["0–4 s", "4–8 s", "8–10 s"] + assert sorted(set(written["bin"])) == ["0-4 s", "4-8 s", "8-10 s"] assert len(written) == 6 + # Excel reads it as UTF-8 only with the byte-order mark. + assert out.read_bytes().startswith(b"\xef\xbb\xbf") def test_cancelling_the_dialog_writes_nothing(self, qtbot, tmp_path, monkeypatch): win = self._win(qtbot, tmp_path)