Responsive units by the ZETA test of the dataset's code: the character, person and location labels¶
Paper: Gerken, Darcher et al. (2025), eLife reviewed preprint 106758 v1, doi:10.7554/eLife.106758.1, full text and Table S1. Dandiset: DANDI:001616, version 0.260702.0824.
The main notebook (reproduce_gerken_2025.ipynb) tests the responsive-unit claims with the test the
Methods describe (Wilcoxon signed-rank on 19 bins, Simes), finds it far from calibrated on this
data, and blocks those claims. The dataset's code (mormannlab/SUMMER
at commit 59fc242, calc_responsive_units.py) uses the ZETA test (zetapy) in its place. This notebook
runs that test for every unit and the character, person and location labels, to show whether the counts of Table S1 could come from it;
responsive_zeta_cuts.ipynb does the same for camera and scene cuts. It is a variant and decides no outcome.
As in the code: spike times and onsets on the recording clock (patient-aligned annotations), onsets preceded by at least 1 s without the label, a window of 1 s after the onset. The code uses 1000 resamplings; here 100, for time. Chance is checked with the same test at random event times, as many as the label's onsets, for 800 random units. The notebook takes about 10 minutes.
import os
for var in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
os.environ.setdefault(var, "1")
import logging
import time
import warnings
import numpy as np
import pandas as pd
from joblib import Parallel, delayed
from zetapy import zetatest
import summer
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
warnings.filterwarnings("ignore")
started_all = started = time.time()
sessions = summer.load_all()
units = summer.unit_table(sessions)
print(f"loaded {len(sessions)} files, {len(units)} units, from dandiset {summer.DANDISET_ID} "
f"version {summer.DANDISET_VERSION} in {time.time() - started:.0f} s")
loaded 29 files, 2286 units, from dandiset 001616 version 0.260702.0824 in 115 s
# Onsets per patient, on the recording clock: segments of the label preceded by at least 1000 ms
# without it (the code's rule, no minimum duration), with the window inside the logged playback.
LABEL_ORDER = ["Tom", "Summer", "McKenzie", "Persons", "Summer Presence", "Indoor", "Outdoor"]
def code_onsets(s, name):
label, value = summer.LABELS[name]
a = s["annotations_patient_aligned"]
a = a[(a.label == label) & (a.value == value)].sort_values("start")
start, stop = a.start.to_numpy(), a.stop.to_numpy()
w = s["cleaned_watchlogs"]
on = start[start - 1000 > np.r_[-np.inf, stop[:-1]]]
return on[(on - 1000 >= w.rec_ms.iloc[0]) & (on + 1000 <= w.rec_ms.iloc[-1])]
onsets = {(int(s["subject_id"]), name): code_onsets(s, name) for s in sessions for name in LABEL_ORDER}
print("units per region:", units.region.value_counts().to_dict(), "; patients:", units.patient.nunique())
print("onset times of `annotations_patient_aligned` lie between",
f"{sessions[0]['annotations_patient_aligned'].start.min():.0f} and {sessions[0]['annotations_patient_aligned'].start.max():.0f}",
"for patient 1 (milliseconds of the recording clock; its watchlog ends at",
f"{sessions[0]['cleaned_watchlogs'].rec_ms.iloc[-1]:.0f} ms)")
print("\nonsets per patient, range over the 29 patients:")
print(pd.DataFrame({name: {"min": min(len(onsets[(p, name)]) for p in units.patient.unique()),
"max": max(len(onsets[(p, name)]) for p in units.patient.unique())}
for name in LABEL_ORDER}).T.to_string())
units per region: {'H': 846, 'A': 559, 'EC': 493, 'PHC': 293, 'PIC': 68, 'other': 27} ; patients: 29
onset times of `annotations_patient_aligned` lie between 0 and 8221608 for patient 1 (milliseconds of the recording clock; its watchlog ends at 8244131 ms)
onsets per patient, range over the 29 patients:
min max
Tom 158 186
Summer 115 140
McKenzie 52 60
Persons 71 84
Summer Presence 40 51
Indoor 43 54
Outdoor 47 62
ZETA for every unit¶
def zeta_p(spikes_ms, events_ms, seed, resamplings=100):
logging.disable(logging.WARNING) # zetapy logs every unit with too few spikes; those get p = 1
np.random.seed(seed)
if len(events_ms) < 2:
return 1.0
return float(zetatest(spikes_ms / 1000, events_ms / 1000, 1, resampling_number=resamplings)[0])
file_of = {int(s["subject_id"]): s for s in sessions}
def spikes_of(i):
return file_of[units.patient[i]]["spikes_ms"][units["index"][i]]
started = time.time()
tasks = [(i, name) for name in LABEL_ORDER for i in units.index]
ps = Parallel(n_jobs=-1, batch_size=16)(delayed(zeta_p)(spikes_of(i), onsets[(units.patient[i], name)], 7 + j)
for j, (i, name) in enumerate(tasks))
Z = pd.DataFrame({"unit": [i for i, _ in tasks], "label": [n for _, n in tasks], "p": ps}).pivot(index="unit", columns="label", values="p")
print(f"{len(tasks)} ZETA tests in {time.time() - started:.0f} s")
16002 ZETA tests in 264 s
# Chance: random event times on the recording clock within each patient's playback, as many as
# the label's onsets for that patient, for 800 random units.
rng_null = np.random.default_rng(51)
null_units = np.sort(rng_null.choice(units.index, 800, replace=False))
started = time.time()
tasks = []
for name in LABEL_ORDER:
for s in sessions:
p = int(s["subject_id"])
w = s["cleaned_watchlogs"]
ev = np.sort(rng_null.uniform(w.rec_ms.iloc[0] + 1000, w.rec_ms.iloc[-1] - 1000, len(onsets[(p, name)])))
tasks += [(i, name, ev) for i in null_units if units.patient[i] == p]
ps = Parallel(n_jobs=-1, batch_size=16)(delayed(zeta_p)(spikes_of(i), ev, 50000 + j) for j, (i, name, ev) in enumerate(tasks))
Zn = pd.DataFrame({"unit": [i for i, _, _ in tasks], "label": [n for _, n, _ in tasks], "p": ps}).pivot(index="unit", columns="label", values="p")
print(f"{len(tasks)} ZETA tests at random times in {time.time() - started:.0f} s")
print(f"\nfraction of the 800 units 'responsive' at random times, per label:")
print(pd.DataFrame({f"alpha_{a}": (Zn[LABEL_ORDER] <= a).mean() for a in summer.ALPHAS}).round(4).to_string())
print("all labels together:", {a: round(float((Zn.to_numpy() <= a).mean()), 4) for a in summer.ALPHAS})
5600 ZETA tests at random times in 93 s
fraction of the 800 units 'responsive' at random times, per label:
alpha_0.001 alpha_0.01
label
Tom 0.0 0.0038
Summer 0.0 0.0012
McKenzie 0.0 0.0038
Persons 0.0 0.0012
Summer Presence 0.0 0.0012
Indoor 0.0 0.0025
Outdoor 0.0 0.0012
all labels together: {0.001: 0.0, 0.01: 0.0021}
Counts against Table S1¶
Table S1 gives, at alpha 0.001, the totals over all 2286 units and the percentages per region for camera cuts.
table_s1 = {"Tom": 216, "Summer": 46, "McKenzie": 0, "Persons": 6, "Tom Faces": 843, "Summer Faces": 634,
"Summer Presence": 0, "Camera Cuts": 1025, "Scene Cuts": 106, "Indoor": 0, "Outdoor": 4}
rows = []
for name in LABEL_ORDER:
row = {"label": name, "paper_total_0.001": table_s1[name]}
for a in summer.ALPHAS:
r = Z[name] <= a
row[f"total_{a}"] = int(r.sum())
row[f"pct_{a}"] = round(100 * float(r.mean()), 2)
for region in ("A", "EC", "H", "PHC"):
row[f"{region}_{a}"] = int(r[units.region == region].sum())
rows.append(row)
print(pd.DataFrame(rows).to_string(index=False))
label paper_total_0.001 total_0.001 pct_0.001 A_0.001 EC_0.001 H_0.001 PHC_0.001 total_0.01 pct_0.01 A_0.01 EC_0.01 H_0.01 PHC_0.01
Tom 216 13 0.57 6 0 0 4 47 2.06 20 2 7 14
Summer 46 11 0.48 6 0 1 3 44 1.92 17 4 11 8
McKenzie 0 1 0.04 0 0 0 1 8 0.35 2 1 2 3
Persons 6 1 0.04 1 0 0 0 34 1.49 12 5 6 8
Summer Presence 0 0 0.00 0 0 0 0 15 0.66 5 1 6 1
Indoor 0 2 0.09 0 0 1 1 22 0.96 10 1 5 6
Outdoor 4 8 0.35 1 1 1 3 35 1.53 5 4 8 14
Result: a variant; the outcomes are those of the main notebook. The ZETA test is conservative here: none of 800 units responds at random times at alpha 0.001 for any label, and 0.12 to 0.38% at 0.01, below the nominal 1%. Its counts are therefore lower bounds. At alpha 0.001 it counts 11 units responsive to Summer, 13 to Tom, 1 to Persons and 8 to Outdoor (Table S1: 46, 216, 6 and 4), and 2 to Indoor, 1 to McKenzie and 0 to Summer Presence (Table S1: 0 for each).
print(f"total run time {(time.time() - started_all) / 60:.1f} min")
total run time 7.9 min