diff --git a/examples/catfr/cohort.py b/examples/catfr/cohort.py index 630838c..2c6259c 100644 --- a/examples/catfr/cohort.py +++ b/examples/catfr/cohort.py @@ -105,8 +105,11 @@ def main(): parser.add_argument("--source-lag", type=int, default=2) parser.add_argument("--n-perm", type=int, default=5000) parser.add_argument("--threads", type=int, default=4) + parser.add_argument("--estimator", default="miller_madow", help="layer estimator of the comparisons") + parser.add_argument("--tag", default="", help="write results to cohort//") args = parser.parse_args() - COHORT.mkdir(parents=True, exist_ok=True) + out_dir = COHORT / args.tag if args.tag else COHORT + out_dir.mkdir(parents=True, exist_ok=True) rows = choose(args.min_trials, args.per_roi) print(f"{len(rows)} eligible subjects", flush=True) @@ -142,6 +145,7 @@ def main(): n_perm=500, curtail=10, map_fn=thread_map(args.threads), + estimator=args.estimator, prng=0, ) per_subject.append( @@ -161,22 +165,35 @@ def main(): fixed = fixed_network(args.source_lag) between = compare_networks( - fixed, fixed, data_rec, data_forg, design="between", dependent=True, n_perm=args.n_perm, prng=0 + fixed, + fixed, + data_rec, + data_forg, + design="between", + dependent=True, + n_perm=args.n_perm, + estimator=args.estimator, + prng=0, ) - between.dataset.to_netcdf(COHORT / "between.nc") + between.dataset.to_netcdf(out_dir / "between.nc") from infoflow.compare import _edge_values per_unit = {} for (s, t), (x_vars, w_vars) in between.edges.items(): per_unit[f"{ROIS[s]}->{ROIS[t]}"] = { - "recalled": [[_edge_values(d, t, x_vars, w_vars)[k] for k in LAYERS] for d in data_rec], - "forgotten": [[_edge_values(d, t, x_vars, w_vars)[k] for k in LAYERS] for d in data_forg], + "recalled": [ + [_edge_values(d, t, x_vars, w_vars, estimator=args.estimator)[k] for k in LAYERS] for d in data_rec + ], + "forgotten": [ + [_edge_values(d, t, x_vars, w_vars, estimator=args.estimator)[k] for k in LAYERS] for d in data_forg + ], } summary = { "subjects": per_subject, "skipped": skipped, "source_lag": args.source_lag, "n_perm": args.n_perm, + "estimator": args.estimator, "between": { f"{ROIS[s]}->{ROIS[t]}": { k: { @@ -189,7 +206,7 @@ def main(): }, "per_unit": per_unit, } - (COHORT / "summary.json").write_text(json.dumps(summary, indent=2, default=str)) + (out_dir / "summary.json").write_text(json.dumps(summary, indent=2, default=str)) print(json.dumps(summary["between"], indent=1), flush=True) diff --git a/examples/catfr/common.py b/examples/catfr/common.py index c0ee3de..909bf85 100644 --- a/examples/catfr/common.py +++ b/examples/catfr/common.py @@ -13,6 +13,8 @@ import re from pathlib import Path +import numpy as np + DATASET = "ds004809" BUCKET = f"s3://openneuro.org/{DATASET}/" ROOT = Path(os.environ.get("CATFR_ROOT", Path.home() / "data" / "openneuro" / DATASET)) @@ -114,6 +116,58 @@ def word_events(events_tsv, sfreq): return words +def retrieval_windows(events_tsv, sfreq, length=1.6, clear=1.8, margin=1.0, step=1.6, after=1.6): + """ + Start times (seconds from the recording start) of retrieval windows. + + Returns ``(recall, deliberation)``. A recall window is the `length` seconds before a + correct, first recall of a list item, inside its recall period, with no other + vocalization starting in the `clear` seconds before the recall. Deliberation windows + tile the recall periods every `step` seconds with no vocalization from `margin` + seconds before to `after` seconds after them. + + With the defaults no two windows share a sample: deliberation windows do not + overlap each other (``step >= length``) or any recall window (``after >= length``). + Overlapping windows duplicate samples, which bias-corrected entropy estimators read + as better-sampled distributions, inflating transfer entropy and synergy. + """ + rows = [r for r in read_tsv(events_tsv) if int(r["sample"]) >= 0] + words = {} + for r in rows: + if r["trial_type"] == "WORD" and int(r["list"]) >= 1: + words.setdefault(r["list"], set()).add(r["item_name"].upper()) + recall, deliberation = [], [] + starts = {r["list"]: int(r["sample"]) / sfreq for r in rows if r["trial_type"] == "REC_START"} + ends = {r["list"]: int(r["sample"]) / sfreq for r in rows if r["trial_type"] == "REC_END"} + for lst, start in starts.items(): + if lst not in ends or lst not in words: + continue + end = ends[lst] + vocal = sorted( + int(r["sample"]) / sfreq + for r in rows + if r["list"] == lst + and r["trial_type"] in ("REC_WORD", "REC_WORD_VV") + and start <= int(r["sample"]) / sfreq <= end + ) + said = set() + for r in rows: + if r["list"] != lst or r["trial_type"] != "REC_WORD": + continue + t, item = int(r["sample"]) / sfreq, r["item_name"].upper() + first = item in words[lst] and item not in said + said.add(item) + earlier = [v for v in vocal if v < t] + if first and t - length >= start and (not earlier or t - earlier[-1] >= clear): + recall.append(t - length) + s = start + margin + while s + length + after <= end: + if not any(s - margin <= v <= s + length + after for v in vocal): + deliberation.append(s) + s += step + return np.array(recall), np.array(deliberation) + + def session_files(subject, session, root=ROOT): base = root / f"sub-{subject}" / f"ses-{session}" / "ieeg" stem = f"sub-{subject}_ses-{session}_task-catFR1" diff --git a/examples/catfr/flow.py b/examples/catfr/flow.py new file mode 100644 index 0000000..989493e --- /dev/null +++ b/examples/catfr/flow.py @@ -0,0 +1,224 @@ +""" +Pre-specified directed flow between MTL and PFC, across the catFR1 cohort. + +Hypothesis (fixed before looking at the data): information flows between the medial +temporal lobe and prefrontal cortex in either direction, and that flow differs with +memory. For every eligible subject (``cohort.choose``), each session is downloaded, +epoched by ``prepare.prepare_phases`` and deleted again, and the regions' HFA is +averaged into three nodes (``cohort.region_trials``) and cut into 4 equal-frequency bins. + +Tested edges: MTL -> PFC and PFC -> MTL with the source at lag 1, 2, 5 or 10 samples +(20, 40, 100, 200 ms), conditioned on the target's past (lags 1-2) and on LTC at the +source lag. Layers are estimated with Chao--Wang--Jost. Two questions per phase: + +* *existence*: real data against the same data with the source's trials shuffled + (within-trial coupling broken, each region's stimulus-locked time course kept); +* *contrast*: ``encoding``, recalled vs forgotten words; ``retrieval``, the 1.6 s before + correct recalls vs deliberation windows. + +Both are paired across subjects (``compare_networks(design="between", dependent=True)``), +with Benjamini--Hochberg over the eight direction x lag tests of each layer. The +primary layer is intrinsic flow. Subjects need ``--min-epochs`` epochs per condition. + + uv run --project ../.. --with mne --with netcdf4 --with "decorator<5.3" python flow.py prepare + uv run --project ../.. --with netcdf4 python flow.py analyze +""" + +import argparse +import json +import subprocess +from types import SimpleNamespace + +import numpy as np +from common import BUCKET, OUT, ROIS, session_files + +from infoflow.compare import _edge_values, compare_networks +from infoflow.data import DiscreteData +from infoflow.layers import LAYERS +from infoflow.preprocess import EqualFrequency, preprocess +from infoflow.selection import TargetSkeleton +from infoflow.stats import benjamini_hochberg + +FLOW = OUT / "flow" +LAGS = (1, 2, 5, 10) +DIRECTIONS = ((0, 2), (2, 0)) # (source, target) indices into ROIS: MTL -> PFC, PFC -> MTL +ESTIMATOR = "chao_wang_jost" +PHASES = {"encoding": ("recalled", "forgotten"), "retrieval": ("pre-recall", "deliberation")} + + +def prepare(n_subjects, per_roi): + from cohort import choose + from prepare import prepare_phases + + FLOW.mkdir(parents=True, exist_ok=True) + rows = choose(40, per_roi)[:n_subjects] + for row in rows: + stem = f"sub-{row['subject']}_ses-{row['session']}" + if (FLOW / f"{stem}_retrieval.nc").exists(): + continue + files = session_files(row["subject"], row["session"]) + downloaded = not files["edf"].exists() + try: + if downloaded: + subprocess.run(["aws", "s3", "cp", "--no-sign-request", "--only-show-errors", BUCKET + files["key"], + str(files["edf"])], check=True) # fmt: skip + e, r = prepare_phases(row["subject"], row["session"], per_roi, FLOW) + print(f"{stem}: encoding {int(e.recalled.sum())}/{int((~e.recalled).sum())}, " + f"retrieval {int(r.recalled.sum())}/{int((~r.recalled).sum())}", flush=True) # fmt: skip + except (ValueError, subprocess.CalledProcessError) as error: + print(f"{stem}: skipped ({error})", flush=True) + finally: + if downloaded: + files["edf"].unlink(missing_ok=True) + + +def edge_network(source, target, lag): + skeleton = { + target: TargetSkeleton( + target=target, target_past=[(target, 1), (target, 2)], sources=[(source, lag)], + conditionals=[(1, lag)], significant=True, + ) + } # fmt: skip + return SimpleNamespace(names=list(ROIS), skeleton=skeleton) + + +def shuffle_source(data, source, rng): + """`data` with the source process's trials permuted (other processes untouched).""" + order = rng.permutation(data.n_trials) + + def swap(arrays): + out = [a.copy() for a in arrays] + for i, j in enumerate(order): + out[i][source] = arrays[j][source] + return out + + return DiscreteData( + past=swap(data.past), present=swap(data.present), past_offset=data.past_offset, + present_offset=data.present_offset, past_alphabet=data.past_alphabet, + present_alphabet=data.present_alphabet, lag_step=data.lag_step, names=data.names, + discretizers=data.discretizers, + ) # fmt: skip + + +def retrieval_mask(path): + """ + The retrieval epochs that :func:`common.retrieval_windows` keeps. + + Epochs written before their start times were stored came from overlapping + deliberation windows (every 0.8 s, 1 s from vocalizations); their starts are + rebuilt from the events and only the current, non-overlapping windows are kept. + """ + import xarray as xr + from common import retrieval_windows, sampling_rate + + ds = xr.open_dataset(path) + files = session_files(ds.attrs["subject"], ds.attrs["session"]) + recall, deliberation = retrieval_windows(files["events"], sampling_rate(files)) + if "start" in ds.coords: + starts = ds.start.values + else: + old = np.concatenate(retrieval_windows(files["events"], sampling_rate(files), step=0.8, after=1.0)) + if len(old) != ds.sizes["trial"]: + raise ValueError(f"{path.name}: cannot match {ds.sizes['trial']} epochs to {len(old)} windows") + starts = old + keep = np.isin(np.round(starts, 6), np.round(np.concatenate([recall, deliberation]), 6)) + return keep + + +def subject_data(path, min_epochs, demean=False): + from cohort import region_trials + + trials, label = region_trials(path) + if path.name.endswith("_retrieval.nc"): + keep = retrieval_mask(path) + trials, label = [t for t, k in zip(trials, keep, strict=True) if k], label[keep] + a, b = np.flatnonzero(label), np.flatnonzero(~label) + n = min(len(a), len(b)) + if n < min_epochs: + return None + rng = np.random.default_rng(0) + a, b = np.sort(rng.choice(a, n, replace=False)), np.sort(rng.choice(b, n, replace=False)) + amplitude = [ + float(np.stack(trials)[a][:, :, i].mean() - np.stack(trials)[b][:, :, i].mean()) for i in range(len(ROIS)) + ] + if demean: + # Remove each condition's mean time course per region: what is left cannot come + # from a shared stimulus- or response-locked trajectory. + x = np.stack(trials) + for idx in (a, b): + x[idx] -= x[idx].mean(axis=0) + trials = list(x) + data = preprocess(trials, names=list(ROIS), max_lag=12, discretizer=EqualFrequency(4), prng=0).data + both = data.subset_trials(list(np.concatenate([a, b]))) + return {"a": data.subset_trials(list(a)), "b": data.subset_trials(list(b)), "both": both, "n": int(n), + "amplitude": amplitude} # fmt: skip + + +def analyze(min_epochs, n_perm, demean=False): + results = {} + for phase, (cond_a, cond_b) in PHASES.items(): + subjects = {} + for path in sorted(FLOW.glob(f"sub-*_{phase}.nc")): + d = subject_data(path, min_epochs, demean) + if d is not None: + subjects[path.name.split(f"_{phase}")[0]] = d + names = sorted(subjects) + units = [subjects[s] for s in names] + amplitude = np.array([u["amplitude"] for u in units]) + out = {"subjects": names, "n_per_condition": [u["n"] for u in units], + "amplitude": {r: {"mean": float(amplitude[:, i].mean()), "positive": int((amplitude[:, i] > 0).sum())} + for i, r in enumerate(ROIS)}, "tests": []} # fmt: skip + rng = np.random.default_rng(1) + for s, t in DIRECTIONS: + for lag in LAGS: + net = edge_network(s, t, lag) + x, w = [(s, lag)], [(1, lag), (t, 1), (t, 2)] + real = [u["both"] for u in units] + surrogate = [shuffle_source(u["both"], s, rng) for u in units] + for question, (da, db) in { + "existence": (real, surrogate), + "contrast": ([u["a"] for u in units], [u["b"] for u in units]), + }.items(): + cmp = compare_networks(net, net, da, db, design="between", dependent=True, n_perm=n_perm, + estimator=ESTIMATOR, prng=0) # fmt: skip + values = np.array([[[_edge_values(d, t, x, w, estimator=ESTIMATOR)[k] for k in LAYERS] + for d in pair] for pair in (da, db)]) # fmt: skip + out["tests"].append({ + "question": question, "edge": f"{ROIS[s]}->{ROIS[t]}", "lag": lag, + **{k: {"a": float(cmp.dataset["a"].sel(layer=k).values[s, t]), + "b": float(cmp.dataset["b"].sel(layer=k).values[s, t]), + "pvalue": float(cmp.dataset["pvalue"].sel(layer=k).values[s, t]), + "positive": int((values[0, :, i] > values[1, :, i]).sum())} + for i, k in enumerate(LAYERS)}, + }) # fmt: skip + print(phase, question, ROIS[s], "->", ROIS[t], lag, + {k: round(out["tests"][-1][k]["pvalue"], 3) for k in LAYERS}, flush=True) # fmt: skip + for question in ("existence", "contrast"): + tests = [x for x in out["tests"] if x["question"] == question] + for k in LAYERS: + _, q = benjamini_hochberg([x[k]["pvalue"] for x in tests], 0.05) + for x, qv in zip(tests, q, strict=True): + x[k]["qvalue"] = float(qv) + out["conditions"] = [cond_a, cond_b] + results[phase] = out + (FLOW / f"results{'_demeaned' if demean else ''}.json").write_text(json.dumps(results, indent=2)) + return results + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("step", choices=("prepare", "analyze")) + parser.add_argument("--subjects", type=int, default=200) + parser.add_argument("--per-roi", type=int, default=10) + parser.add_argument("--min-epochs", type=int, default=30) + parser.add_argument("--n-perm", type=int, default=10000) + parser.add_argument("--demean", action="store_true", help="remove each condition's mean time course first") + args = parser.parse_args() + if args.step == "prepare": + prepare(args.subjects, args.per_roi) + else: + analyze(args.min_epochs, args.n_perm, args.demean) + + +if __name__ == "__main__": + main() diff --git a/examples/catfr/prepare.py b/examples/catfr/prepare.py index 5394f6e..fd3955f 100644 --- a/examples/catfr/prepare.py +++ b/examples/catfr/prepare.py @@ -20,11 +20,12 @@ import json from collections import defaultdict from itertools import zip_longest +from pathlib import Path import mne import numpy as np import xarray as xr -from common import OUT, ROIS, pair_regions, session_files, word_events +from common import OUT, ROIS, pair_regions, retrieval_windows, session_files, word_events EPOCH = 1.6 PAD = 1.0 @@ -64,6 +65,57 @@ def epochs_of(raw, onsets, tmin, tmax): return np.stack([data[:, s : s + n] for s in starts[keep]]), keep +def hfa_epochs(hf, starts): + """ + HFA of ``EPOCH``-second windows starting at `starts`, from a 500 Hz recording. + + Morlet power is computed on windows padded by ``PAD`` seconds each side and cropped, + so the window has no edge artifacts. + """ + padded, keep = epochs_of(hf, np.asarray(starts), -PAD, EPOCH + PAD) + power = mne.time_frequency.tfr_array_morlet(padded, 500.0, HFA_FREQS, n_cycles=7, output="power", verbose="error") + lo = int(PAD * 500) + env = np.log10(power.mean(axis=2))[:, :, lo : lo + int(EPOCH * 500)] + step = int(500 / HFA_SFREQ) + return env.reshape(*env.shape[:2], -1, step).mean(axis=-1), keep + + +def prepare_phases(subject, session, per_roi, out_dir): + """ + Encoding and retrieval HFA epochs of one session, for the directed-flow analysis. + + ``encoding``: every non-practice WORD, 0-1.6 s after onset, ``recalled`` as usual. + ``retrieval``: the 1.6 s before correct recalls (``recalled`` True) and deliberation + windows of the recall periods (False), from :func:`common.retrieval_windows`. + """ + files = session_files(subject, session) + regions = pair_regions(files["channels"], files["electrodes"]) + raw = mne.io.read_raw_edf(files["edf"], preload=False, verbose="error") + sfreq = raw.info["sfreq"] + regions = {k: v for k, v in regions.items() if k in raw.ch_names} + picks = pick_pairs(regions, per_roi) + raw.pick(picks).load_data(verbose="error") + raw.notch_filter(np.arange(60, sfreq / 2 - 1, 60), verbose="error") + hf = raw.resample(500.0, verbose="error") + attrs = {"subject": subject, "session": session, "source_sfreq": sfreq, "signal": "hfa"} + + words = word_events(files["events"], sfreq) + env, keep = hfa_epochs(hf, [w["onset"] for w in words]) + encoding = to_dataset(env, HFA_SFREQ, words, keep, picks, regions, {**attrs, "phase": "encoding"}) + + recall, deliberation = retrieval_windows(files["events"], sfreq) + starts = np.concatenate([recall, deliberation]) + labels = [{"recalled": i < len(recall), "list": -1, "serialpos": -1, "item": ""} for i in range(len(starts))] + env, keep = hfa_epochs(hf, starts) + retrieval = to_dataset(env, HFA_SFREQ, labels, keep, picks, regions, {**attrs, "phase": "retrieval"}) + retrieval = retrieval.assign_coords(start=("trial", starts[keep])) + + stem = Path(out_dir) / f"sub-{subject}_ses-{session}" + encoding.to_netcdf(f"{stem}_encoding.nc") + retrieval.to_netcdf(f"{stem}_retrieval.nc") + return encoding, retrieval + + def to_dataset(values, sfreq, words, keep, picks, regions, attrs): kept = [w for w, k in zip(words, keep) if k] return xr.Dataset( @@ -101,12 +153,7 @@ def prepare(subject, session, per_roi): broadband = to_dataset(values, BROADBAND_SFREQ, words, keep, picks, regions, {**attrs, "signal": "broadband"}) hf = raw.copy().resample(500.0, verbose="error") - padded, keep = epochs_of(hf, onsets, -PAD, EPOCH + PAD) - power = mne.time_frequency.tfr_array_morlet(padded, 500.0, HFA_FREQS, n_cycles=7, output="power", verbose="error") - lo = int(PAD * 500) - env = np.log10(power.mean(axis=2))[:, :, lo : lo + int(EPOCH * 500)] - step = int(500 / HFA_SFREQ) - env = env.reshape(*env.shape[:2], -1, step).mean(axis=-1) + env, keep = hfa_epochs(hf, onsets) hfa = to_dataset(env, HFA_SFREQ, words, keep, picks, regions, {**attrs, "signal": "hfa"}) stem = OUT / f"sub-{subject}_ses-{session}" diff --git a/infoflow/compare.py b/infoflow/compare.py index aaea8ea..3094677 100644 --- a/infoflow/compare.py +++ b/infoflow/compare.py @@ -43,13 +43,13 @@ class NetworkComparison: edges: dict -def _edge_values(data, target, x_vars, w_vars, restarts=1): +def _edge_values(data, target, x_vars, w_vars, restarts=1, estimator="miller_madow"): variables = list(dict.fromkeys(list(x_vars) + list(w_vars))) r = realizations(data, target, variables, max_lag=max(v[1] for v in variables)) if len(r) < 10: return dict.fromkeys(LAYERS, np.nan) w = r.columns(list(w_vars)) if w_vars else None - flow = edge_flows(r.present, r.columns(list(x_vars)), w, estimator="miller_madow", restarts=restarts, prng=0) + flow = edge_flows(r.present, r.columns(list(x_vars)), w, estimator=estimator, restarts=restarts, prng=0) return {k: flow.as_dict()[k] for k in LAYERS} @@ -95,6 +95,7 @@ def compare_networks( alpha=0.05, curtail=None, map_fn=map, + estimator="miller_madow", prng=None, ): """ @@ -125,6 +126,10 @@ def compare_networks( map_fn : callable ``map``-like function over edges (e.g. from :mod:`infoflow.parallel`). Each edge's permutations use a seed drawn up front, so results do not depend on it. + estimator : str + The layers' estimator (see :func:`~infoflow.measures.cmi_from_joint`). + Miller--Madow by default; :data:`~infoflow.measures.DEFAULT_LAYER_ESTIMATOR` + has less bias in large contexts but costs about 13 times as much per entropy. prng : None, int, Generator Returns @@ -149,7 +154,10 @@ def compare_networks( if design == "between": # Per-subject values are fixed; permutations reassign subjects to conditions. - per_unit = {e: np.array([[_edge_values(u, e[1], *edges[e])[k] for k in LAYERS] for u in units]) for e in edges} + per_unit = { + e: np.array([[_edge_values(u, e[1], *edges[e], estimator=estimator)[k] for k in LAYERS] for u in units]) + for e in edges + } def statistic(lab, e): v = per_unit[e] @@ -158,8 +166,12 @@ def statistic(lab, e): else: def statistic(lab, e): - a = _edge_values(_pool([u for u, g in zip(units, lab, strict=True) if g == 0]), e[1], *edges[e]) - b = _edge_values(_pool([u for u, g in zip(units, lab, strict=True) if g == 1]), e[1], *edges[e]) + a = _edge_values( + _pool([u for u, g in zip(units, lab, strict=True) if g == 0]), e[1], *edges[e], estimator=estimator + ) + b = _edge_values( + _pool([u for u, g in zip(units, lab, strict=True) if g == 1]), e[1], *edges[e], estimator=estimator + ) return np.array([a[k] for k in LAYERS]), np.array([b[k] for k in LAYERS]) def permuted(edge_rng): @@ -206,6 +218,12 @@ def test(job): "significant": (("layer", "source", "target"), np.stack([significant[k] for k in LAYERS])), }, coords={"layer": list(LAYERS), "source": names, "target": names}, - attrs={"design": design, "dependent": int(dependent), "n_perm": n_perm, "curtail": curtail or 0}, + attrs={ + "design": design, + "dependent": int(dependent), + "n_perm": n_perm, + "curtail": curtail or 0, + "estimator": estimator, + }, # fmt: skip ) return NetworkComparison(ds, edges) diff --git a/reports/catfr_directed_flow.md b/reports/catfr_directed_flow.md new file mode 100644 index 0000000..7012698 --- /dev/null +++ b/reports/catfr_directed_flow.md @@ -0,0 +1,141 @@ +# Directed flow between MTL and PFC in catFR1 + +*Branch `catfr-directed-flow`, November 2026. Scripts: `examples/catfr/flow.py` (prepare, +analyze), `validation/synergy_null.py`; follows `reports/catfr_pilot.md` and +`reports/dit_estimators.md`.* + +## Summary + +- **No directed MTL-PFC flow that differs with memory.** Across 88 subjects (encoding) and + 56 (retrieval), the pre-specified test of intrinsic flow MTL -> PFC and PFC -> MTL at + 20-200 ms found no difference between recalled and forgotten words or between the + 1.6 s before correct recalls and deliberation, in any layer. Nor did intrinsic flow + exceed a shuffled-source surrogate in the pre-specified analysis. +- **Shared flow is there throughout**, in both directions, at every lag and in both phases + (real > shuffled source, q < 0.01): MTL and PFC co-fluctuate within trials, beyond what + each region's time course explains, but not in a way the target's past and LTC do not + already carry or that is directed. +- **Two artifacts found and removed along the way.** Overlapping retrieval windows made + bias-corrected estimates look like strong synergy (and a strong recall effect); a shared-flow + recall effect at 200 ms came from both regions following the same pre-recall trajectory. +- **Exploratory:** with each condition's mean time course removed, intrinsic flow exceeded + the surrogate in 3 of 8 encoding tests (MTL -> PFC at 20 and 40 ms, PFC -> MTL at 200 ms, + q = 0.012-0.031), by about 0.0001 bits. This was not pre-specified and is not a finding. +- **Open issues closed:** `compare_networks(estimator=...)` (the 25-subject memory contrast + rerun with Chao-Wang-Jost: still null); the synergistic null excess traced to chance; + infoflow passes dit's new `rng` keyword everywhere (#14). + +## 1. Design (fixed before the analysis) + +- **Subjects:** every catFR1 session that `cohort.choose` accepts (one per subject, two or + more non-overlapping bipolar pairs in each of MTL, LTC and PFC, 40 or more trials per + encoding condition): 91 subjects, 88 prepared (3 have no usable epochs). Each EDF is + downloaded, epoched and deleted (`flow.py prepare`). A subject enters a phase with 30 or + more epochs per condition after balancing: 88 for encoding, 56 for retrieval (median 71 and + 43 epochs per condition). +- **Nodes:** each region's HFA (70-150 Hz log power, 50 Hz), every pair z-scored and then + averaged, cut into 4 equal-frequency bins per subject. +- **Edges:** MTL -> PFC and PFC -> MTL, source at lag 1, 2, 5, 10 samples (20, 40, 100, + 200 ms), conditioned on the target's past (lags 1-2) and LTC at the source lag; layers + estimated with Chao-Wang-Jost. +- **Phases:** *encoding*, 0-1.6 s after each word (recalled vs forgotten); *retrieval*, the + 1.6 s before correct first recalls with no vocalization in the preceding 1.8 s, against + non-overlapping deliberation windows of the recall periods with no vocalization from 1 s + before to 1.6 s after (`common.retrieval_windows`). +- **Questions:** *existence*, real data against the source's trials shuffled (within-trial + coupling broken, stimulus-locked time courses kept); *contrast*, condition A against B. + Both paired across subjects (`compare_networks(design="between", dependent=True)`, + 10 000 sign flips), Benjamini-Hochberg over the eight direction x lag tests per layer. + Primary layer: intrinsic. + +## 2. Results + +Means over subjects, bits; q per layer over the 8 tests. + +**Encoding (88 subjects).** + +| Question | Intrinsic | Synergistic | Shared | +|---|---|---|---| +| Existence (real vs shuffled) | 0.0001-0.0002 vs 0.0001; smallest q 0.23 | 0.0012-0.0021 vs 0.0014-0.0023; smallest q 0.16 | 0.0005-0.0017 vs 0.0001; q <= 0.002 in all 8 | +| Recalled vs forgotten | smallest q 0.93 | smallest q 0.48 | smallest q 0.92 | + +**Retrieval (56 subjects), non-overlapping windows.** + +| Question | Intrinsic | Synergistic | Shared | +|---|---|---|---| +| Existence | smallest q 0.81 | smallest q 0.77 | 0.0008-0.0016 vs 0.0001-0.0003; q <= 0.006 in all 8 | +| Pre-recall vs deliberation | smallest q 0.52 | smallest q 0.55 | PFC -> MTL at 200 ms: 0.0019 vs 0.0008, q = 0.019; others q >= 0.54 | + +The single surviving retrieval contrast disappears once each condition's mean time course +is removed from each region (section 3), so it reflects the shared pre-recall rise (PFC HFA +is higher before recalls in 46 of 56 subjects) rather than coupling. + +**Amplitude, for comparison.** Encoding: PFC HFA higher for recalled words in 67 of 88 +subjects, MTL in 54. Retrieval: PFC higher before recalls in 46 of 56. + +## 3. Checks and artifacts + +- **Overlapping deliberation windows.** The first retrieval run tiled deliberation windows + every 0.8 s (1.6 s long), so neighbours shared half their samples. Duplicated samples + turn singletons into doubletons, which bias-corrected estimators such as Chao-Wang-Jost + read as a better-sampled distribution, inflating transfer entropy and synergy. That run + reported synergy existing at q < 0.001 in 6 of 8 tests and deliberation synergy twice + pre-recall synergy (0.022-0.026 against 0.010-0.015 bits, q < 0.001 in 7 of 8 and 0.021 + in the eighth). With + non-overlapping windows (the same epochs, subselected) both vanish: deliberation synergy + drops to 0.009-0.012 bits and every synergistic q is >= 0.55. `retrieval_windows` now + produces non-overlapping windows and stores each epoch's start. +- **Condition-locked trajectories.** Removing each condition's mean time course per region + before discretization (`flow.py analyze --demean`) leaves no contrast in any layer in + either phase, and keeps shared-flow existence (encoding 8 of 8, retrieval 6 of 8). +- **Exploratory intrinsic flow.** In the de-meaned encoding data, intrinsic flow exceeded the + surrogate for MTL -> PFC at 20 ms (0.00018 vs 0.00005, q = 0.012) and 40 ms (q = 0.031), + and PFC -> MTL at 200 ms (0.00021 vs 0.00008, q = 0.012), with real above shuffled in only + 25-32 of 88 subjects (most others tie at zero, where the estimator clips). Not + pre-specified, inconsistent across lags, and about a ten-thousandth of a bit. + +## 4. Open issues closed + +- **`compare_networks(estimator=...)`.** The layers' estimator is now an option (default + Miller-Madow, unchanged). Rerunning the 25-subject region-level memory contrast with + Chao-Wang-Jost lowers the baseline synergy per condition about fourfold (0.0013-0.0025 + bits, against 0.0056-0.0075 with Miller-Madow) and still finds nothing: every + between-subject q >= 0.31, and no within-subject edge (Miller-Madow had one). +- **Synergistic null excess.** `validation/synergy_null.py` split R1332M's 170 recalled + trials 1000 times, with the network selected once on all of them. Miller-Madow synergy + had 7.0% of p-values below 0.05 (intrinsic 4.6%, shared 4.7%). An independent exact test + on 5000 fresh random splits reproduced the p-values (correlation 0.998, 6.5% below 0.05), + the statistic is deterministic and order-independent, and the observed splits' statistics + follow the permutation distribution (KS p = 0.58): the test is exact and the excess is + chance (binomial p of about 0.02 across six rates examined). With Chao-Wang-Jost, + synergy is nominal (4.6%) but 56% of its p-values are exactly 1: at about 85 trials per + group its synergy is usually clipped to zero, so it has little power for synergy at this + sample size. +- **dit's `prng` keyword.** Rather than restoring the alias in dit, infoflow uses `rng=` + for every dit call: `total_correlation_ksg` since #12, and the surrogates and + `select_markov_order` in #14. + +## 5. Where this leaves directed flow + +With region-averaged HFA, 50 Hz, four bins and 43-71 epochs per condition, any memory-related +change in directed MTL-PFC information flow is below what this design can see; the shared +co-fluctuation is robust but undirected. Directions that might still change that: + +- finer spatial resolution (hippocampal and prefrontal pairs with the strongest amplitude + effects, chosen on a held-out session), since region averaging mixes signals; +- broadband phase-based measures at theta/gamma, where MTL-PFC coupling is usually reported, + instead of HFA envelopes; +- pooling sessions within subjects to raise epochs per condition. + +## Reproducing + +```bash +cd examples/catfr +export CATFR_OUT=~/data/catfr_cohort +uv run --project ../.. --with mne --with netcdf4 --with "decorator<5.3" python flow.py prepare +uv run --project ../.. --with mne --with netcdf4 --with "decorator<5.3" python flow.py analyze +uv run --project ../.. --with mne --with netcdf4 --with "decorator<5.3" python flow.py analyze --demean +uv run --project ../.. --with mne --with netcdf4 --with "decorator<5.3" python cohort.py --estimator chao_wang_jost --tag cwj +uv run --project ../.. --with netcdf4 python ../../validation/synergy_null.py \ + --epochs ~/data/catfr_pilot_v3/sub-R1332M_ses-0_hfa.nc --splits 1000 --estimator miller_madow --out mm.json +``` diff --git a/tests/test_engineering.py b/tests/test_engineering.py index 87cb14a..4ac5e2d 100644 --- a/tests/test_engineering.py +++ b/tests/test_engineering.py @@ -63,6 +63,18 @@ def test_compare_curtailed_and_parallel(): assert np.allclose(5 / p[tested], np.round(5 / p[tested])) +def test_compare_estimator_option(): + a, b = _coupled_trials(0.9, 10, 300, seed=0), _coupled_trials(0.5, 10, 300, seed=1) + net_a, net_b = _net(a), _net(b) + mm = compare_networks(net_a, net_b, a, b, n_perm=39, prng=0) + cwj = compare_networks(net_a, net_b, a, b, n_perm=39, estimator="chao_wang_jost", prng=0) + assert mm.dataset.attrs["estimator"] == "miller_madow" and cwj.dataset.attrs["estimator"] == "chao_wang_jost" + diff = cwj.dataset["difference"].sel(layer="intrinsic").values[0, 1] + assert diff > 0.2 and bool(cwj.dataset["significant"].sel(layer="intrinsic").values[0, 1]) + layer = dict(layer="intrinsic", source="x0", target="x1") + assert float(mm.dataset["a"].sel(**layer)) != float(cwj.dataset["a"].sel(**layer)) + + def test_compare_within_null_is_calibrated_with_pooled_selection(): pooled = _coupled_trials(0.8, 40, 200, seed=7) halves = [] diff --git a/validation/synergy_null.py b/validation/synergy_null.py new file mode 100644 index 0000000..e7c38d2 --- /dev/null +++ b/validation/synergy_null.py @@ -0,0 +1,85 @@ +""" +Calibration of ``compare_networks`` on random splits of one condition (catFR1). + +All of one subject's recalled trials are split into two halves at random many times. +The tested network is selected once on all of them, so it is the same for every split, +and the split is one draw from the permutation distribution that the test compares it +with: for each edge and layer, ``P(p <= alpha) <= alpha`` exactly. The script records +every raw p-value, so their distribution can be checked directly. + + python validation/synergy_null.py --epochs ~/data/catfr_pilot_v3/sub-R1332M_ses-0_hfa.nc \\ + --splits 1000 --estimator miller_madow --out synergy_null_mm.json +""" + +import argparse +import json +import sys +from concurrent.futures import ProcessPoolExecutor +from pathlib import Path + +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "examples" / "catfr")) + +STATE = {} + + +def _init(data, net): + STATE["data"], STATE["net"] = data, net + + +def _split(args): + from infoflow.compare import compare_networks + from infoflow.layers import LAYERS + + seed, pool, half, n_perm, estimator = args + draw = np.random.default_rng(seed).permutation(pool) + data, net = STATE["data"], STATE["net"] + cmp = compare_networks( + net, net, data.subset_trials(list(np.sort(draw[:half]))), data.subset_trials(list(np.sort(draw[half:]))), + n_perm=n_perm, estimator=estimator, prng=seed, + ) # fmt: skip + out = {} + for k in LAYERS: + p = cmp.dataset["pvalue"].sel(layer=k).values + out[k] = [float(p[s, t]) for s, t in zip(*np.nonzero(np.isfinite(p)), strict=True)] + return seed, out + + +def main(): + import xarray as xr + from run import condition_network + + from infoflow.parallel import thread_map + from infoflow.preprocess import preprocess + + parser = argparse.ArgumentParser() + parser.add_argument("--epochs", required=True) + parser.add_argument("--splits", type=int, default=1000) + parser.add_argument("--n-perm", type=int, default=500) + parser.add_argument("--estimator", default="miller_madow") + parser.add_argument("--workers", type=int, default=12) + parser.add_argument("--out", required=True) + args = parser.parse_args() + + ds = xr.open_dataset(args.epochs).load() + names = [f"{r}:{c}" for r, c in zip(ds.region.values, ds.channel.values, strict=True)] + trials = [ds.signal.values[i].astype(float) for i in range(ds.sizes["trial"])] + pool = np.flatnonzero(ds.recalled.values) + pre = preprocess(trials, names=names, max_lag=15, network_resolution="auto", map_fn=thread_map(8), prng=0) + net = condition_network(pre, pool, names, 8) + edges = [(names[p], names[t]) for t, sk in net.skeleton.items() if sk.significant for p in sk.parents()] + print(f"{len(pool)} trials, edges {edges}", flush=True) + jobs = [(1000 + r, pool, len(pool) // 2, args.n_perm, args.estimator) for r in range(args.splits)] + rows = [] + with ProcessPoolExecutor(args.workers, initializer=_init, initargs=(pre.data, net)) as pool_: + for seed, out in pool_.map(_split, jobs, chunksize=4): + rows.append({"seed": seed, "pvalues": out}) + Path(args.out).write_text(json.dumps({"edges": edges, "estimator": args.estimator, "splits": rows})) + for k in ("intrinsic", "synergistic", "shared"): + p = np.array([v for r in rows for v in r["pvalues"][k]]) + print(f"{k:12s} n={len(p)} P(p<0.05)={np.mean(p < 0.05):.3f} P(p<0.01)={np.mean(p < 0.01):.3f}", flush=True) + + +if __name__ == "__main__": + main()