Responsive units by the ZETA test of the dataset's code: camera cuts and scene cuts¶

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 camera cuts and scene cuts, to show whether the counts of Table S1 could come from it; responsive_zeta_labels.ipynb does the same for the character, person and location labels. 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 400 random units. The notebook takes about 15 minutes.

In [1]:
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 111 s
In [2]:
# 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 = ["Camera Cuts", "Scene Cuts"]


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
Camera Cuts  876  988
Scene Cuts   132  158

ZETA for every unit¶

In [3]:
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")

4572 ZETA tests in 515 s
In [4]:
# Chance: random event times on the recording clock within each patient's playback, as many as
# the label's onsets for that patient, for 400 random units.
rng_null = np.random.default_rng(51)
null_units = np.sort(rng_null.choice(units.index, 400, 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 400 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})

800 ZETA tests at random times in 90 s

fraction of the 400 units 'responsive' at random times, per label:
             alpha_0.001  alpha_0.01
label                               
Camera Cuts          0.0      0.0100
Scene Cuts           0.0      0.0025
all labels together: {0.001: 0.0, 0.01: 0.0063}

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.

In [5]:
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))

paper_cc = {"A": 46.90, "EC": 45.68, "H": 38.92, "PHC": 51.74}
rows = []
for region in ("A", "EC", "H", "PHC", "PIC", "other"):
    idx = units.index[units.region == region]
    row = {"region": region, "units": len(idx), "paper_camera_cuts_pct_0.001": paper_cc.get(region)}
    for name in LABEL_ORDER:
        for a in summer.ALPHAS:
            row[f"{name}_{a}_pct"] = round(100 * float((Z.loc[idx, name] <= a).mean()), 2)
    rows.append(row)
print()
print(pd.DataFrame(rows).to_string(index=False))
print("paper, PHC: camera cuts 193 of 373 (51.74%), scene cuts 25 of 373 (6.70%)")
phc = units.index[units.region == "PHC"]
for a in summer.ALPHAS:
    print(f"PHC at alpha {a}: camera cuts {int((Z.loc[phc, 'Camera Cuts'] <= a).sum())} of {len(phc)}, "
          f"scene cuts {int((Z.loc[phc, 'Scene Cuts'] <= a).sum())} of {len(phc)}")

      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
Camera Cuts               1025           79       3.46       18         3       11         35         189      8.27      47       21      36        67
 Scene Cuts                106           19       0.83        3         2        0         12          95      4.16      24       13      20        29

region  units  paper_camera_cuts_pct_0.001  Camera Cuts_0.001_pct  Camera Cuts_0.01_pct  Scene Cuts_0.001_pct  Scene Cuts_0.01_pct
     A    559                        46.90                   3.22                  8.41                  0.54                 4.29
    EC    493                        45.68                   0.61                  4.26                  0.41                 2.64
     H    846                        38.92                   1.30                  4.26                  0.00                 2.36
   PHC    293                        51.74                  11.95                 22.87                  4.10                 9.90
   PIC     68                          NaN                   5.88                  8.82                  0.00                 2.94
 other     27                          NaN                  29.63                 44.44                  7.41                25.93
paper, PHC: camera cuts 193 of 373 (51.74%), scene cuts 25 of 373 (6.70%)
PHC at alpha 0.001: camera cuts 35 of 293, scene cuts 12 of 293
PHC at alpha 0.01: camera cuts 67 of 293, scene cuts 29 of 293

100 against 1000 resamplings¶

The code uses 1000 resamplings. For the PHC units and scene cuts, the test is repeated with 1000.

In [6]:
started = time.time()
p1000 = np.array(Parallel(n_jobs=-1)(delayed(zeta_p)(spikes_of(i), onsets[(units.patient[i], "Scene Cuts")], 90000 + j, 1000)
                                     for j, i in enumerate(phc)))
p100 = Z.loc[phc, "Scene Cuts"].to_numpy()
print(f"{len(phc)} PHC units, scene cuts, 1000 resamplings in {time.time() - started:.0f} s")
for a in summer.ALPHAS:
    print(f"alpha {a}: responsive with 100 resamplings {int((p100 <= a).sum())}, with 1000 {int((p1000 <= a).sum())}; "
          f"units classed differently {int(((p100 <= a) != (p1000 <= a)).sum())}")
lp = np.log10(np.clip(np.c_[p100, p1000], 1e-12, 1))
print(f"correlation of log10 p, 100 against 1000 resamplings: {np.corrcoef(lp.T)[0, 1]:.3f}")

293 PHC units, scene cuts, 1000 resamplings in 77 s
alpha 0.001: responsive with 100 resamplings 12, with 1000 11; units classed differently 1
alpha 0.01: responsive with 100 resamplings 29, with 1000 31; units classed differently 6
correlation of log10 p, 100 against 1000 resamplings: 0.994

Result: a variant; the outcomes are those of the main notebook. The ZETA test is not anticonservative here: none of 400 units responds at random times at alpha 0.001, and 1.0% (camera cuts) and 0.25% (scene cuts) at 0.01. With 400 units, zero hits at 0.001 cannot show calibration at that level, so its counts are best read as lower bounds. It counts 79 units (3.46%) responsive to camera cuts at alpha 0.001 and 189 at 0.01, against 1025 in Table S1; by region A 3.22%, EC 0.61%, H 1.30% and PHC 11.95% (35 of 293) at 0.001, against 46.90, 45.68, 38.92 and 51.74% in the paper. For scene cuts it counts 19 units at 0.001 (paper 106), 12 of them in PHC (4.10%; paper 6.70%), and 29 PHC units (9.90%) at 0.01. With 1000 resamplings the PHC scene cut counts are 11 and 31 (correlation of log p 0.994). The counts of the paper's Table S1 come from neither this test nor the Methods' test as run in the main notebook.

In [7]:
print(f"total run time {(time.time() - started_all) / 60:.1f} min")

total run time 13.2 min