Shin et al. (2025) from dandiset 000248: decoding

The decoding claims of the Neuropixels data (Figs. 3a, 3c, 3d; Extended Data Fig. 7), on DANDI:000248 version 0.240502.2344. The header of reproduce_shin_2025.ipynb gives the paper, the text the claims file was written from, the choice of version and the other claims. Loading and the decoder are in illusion.py.

Claim	Tested here
fig3a-npx-v1-decoding	yes
fig3c-npx-v1-inference	yes (further seeds in decoding_checks_shin_2025.ipynb, reading a, and decoding_seed1_shin_2025.ipynb, decoding_seed2_shin_2025.ipynb, reading b)
fig3d-npx-areas-inference	yes (further seeds as for fig3c)
ed7b-npx-xre-inference	yes
ed7a-npx-no-area-above-v1	yes

import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import time
import warnings

import numpy as np
import pandas as pd
from joblib import Parallel, delayed
from scipy import stats

import illusion as il

warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
N_JOBS = 4

started = time.time()
sessions = il.load_all(ccg=False)
print(f"loaded {len(sessions)} sessions of dandiset {il.DANDISET_ID} version {il.DANDISET_VERSION} "
      f"in {time.time() - started:.0f} s")

loaded 12 sessions of dandiset 000248 version 0.240502.2344 in 132 s

Data check

RS units per area and session, and ICwcfg1 trials of the decoded images.

rows = []
for s in sessions:
    u, b = s["units"], s["blocks"]["ICwcfg1"]
    n = {c: int(np.sum(b["codes"] == c)) for c in [il.IC1, il.LC1, il.LC2, il.IC2, il.TRE1, il.TRE2, il.XRE1, il.XRE2]}
    rows.append({"subject": s["subject"], **{f"{a}_RS": int(((u.area == a) & u.rs).sum()) for a in il.AREAS},
                 "trials IC1/LC1/LC2/IC2": "/".join(str(n[c]) for c in [il.IC1, il.LC1, il.LC2, il.IC2]),
                 "TRE1/TRE2": f"{n[il.TRE1]}/{n[il.TRE2]}", "XRE1/XRE2": f"{n[il.XRE1]}/{n[il.XRE2]}"})
units_tab = pd.DataFrame(rows)
print(units_tab.to_string(index=False))
allR = np.concatenate([s["blocks"]["ICwcfg1"]["R"].ravel() for s in sessions])
print(f"IC response (spikes/s): min {allR.min():.1f}, median {np.median(allR):.1f}, max {allR.max():.1f}")
print("Session-areas with >= 50 RS units:", {a: int((units_tab[f"{a}_RS"] >= 50).sum()) for a in il.AREAS},
      "(paper 12, 8, 11, 10, 12, 11)")

subject  V1_RS  LM_RS  RL_RS  AL_RS  PM_RS  AM_RS trials IC1/LC1/LC2/IC2 TRE1/TRE2 XRE1/XRE2
 619293     89      0    156    106    151    107        400/400/400/400   400/400   400/400
 619296    168      4     87     94     95     93        400/400/400/400   400/400   400/400
 620333    114    119    128     91    102    129        400/400/400/400   400/400   400/400
 620334     68    183    100      0     90     82        400/400/400/400   400/400   400/400
 625545    190    181    122    160    173    134        400/400/400/400   400/400   400/400
 625554    179      0    136    139    193    183        400/400/400/400   400/400   400/400
 625555    314      0    208    159    214    198        400/400/400/400   400/400   400/400
 630506     58    131    140    121    208      0        400/400/400/400   400/400   400/400
 631510    129    105     75     87     70     96        400/400/400/400   400/400   400/400
 631570    142     74      0    105    158    109        400/400/400/400   400/400   400/400
 633229    181    172    130    179    210    142        400/400/400/400   400/400   400/400
 637484    172    122    279      0     67     96        400/400/400/400   400/400   400/400

IC response (spikes/s): min 0.0, median 0.0, max 525.0
Session-areas with >= 50 RS units: {'V1': 12, 'LM': 8, 'RL': 11, 'AL': 10, 'PM': 12, 'AM': 11} (paper 12, 8, 11, 10, 12, 11)

The counts of RS units and trials are those of the main notebook's data check (400 trials of each decoded image per session). The clock of spike times and trial table is checked there.

Decoder

The decoder of the claims file: ICwcfg1, RS units of one area, sessions with at least 50 RS units in the area; IC response; classes IC1, LC1, LC2, IC2 balanced to the smallest class by its first trials; stratified tenfold partition; z-scoring by the training trials; linear SVM, one-vs-one. Reading (a): C = 1; reading (b): C chosen by an inner threefold cross-validation on the training trials over 10^-3 to 1, the Python counterpart of the code's hyperparameter search. The search over C above 1, the one-vs-all coding and the kernel scale (equivalent to C for a linear kernel) is examined in decoding_checks_shin_2025.ipynb on V1, PM and AM. TRE trials are decoded by each fold model and the fractions averaged. The XRE decoders are trained on XRE1 against XRE2 and probed with IC1 and IC2. Each (session, area) has its own seed. Both readings are run on every area.

t0 = time.time()
jobs = []
for i, s in enumerate(sessions):
    u = s["units"]
    for j, a in enumerate(il.AREAS):
        m = ((u.area == a) & u.rs).values
        if m.sum() < 50:
            continue
        b = s["blocks"]["ICwcfg1"]
        jobs.append((s["subject"], a, b["R"][:, m], b["codes"], 1000 * i + j))

def run_decoders(subject, a, R, codes, seed):
    b = {"R": R, "codes": codes}
    m = np.ones(R.shape[1], bool)
    out = []
    for rd in ["a", "b"]:
        r = il.decode(b["R"][:, m], b["codes"], [il.IC1, il.LC1, il.LC2, il.IC2], [il.TRE1, il.TRE2], seed=seed, reading=rd)
        score, icv, lcv = il.inference_score(r["probe"])
        x = il.decode(b["R"][:, m], b["codes"], [il.XRE1, il.XRE2], [il.IC1, il.IC2], seed=seed, reading=rd)
        out.append({"subject": subject, "area": a, "n_units": int(m.sum()), "reading": rd,
                    "acc": r["acc"], "tre_ic": icv, "tre_lc": lcv, "inference": score, "n_per_class": r["n_per_class"],
                    "xre_acc": x["acc"], "xre_score": (x["probe"][0, 0] + x["probe"][1, 1]) / 2,
                    "C_tre": r["C_chosen"], "C_xre": x["C_chosen"], "unconverged": r["n_unconverged"] + x["n_unconverged"]})
    return out

dec = pd.DataFrame([r for rs in Parallel(n_jobs=N_JOBS)(delayed(run_decoders)(*j) for j in jobs) for r in rs])
warnings.filterwarnings("ignore")
print(f"{len(jobs)} session-area decoders, both readings, in {time.time() - t0:.0f} s")
print("Sessions per area:", dec[dec.reading == "a"].groupby("area").size().reindex(il.AREAS).to_dict(),
      "(paper 12, 8, 11, 10, 12, 11)")
print("Trials per class:", sorted(dec.n_per_class.unique()))
print("Fits stopped at max_iter without converging:", dec.groupby("reading").unconverged.sum().to_dict(),
      "of", len(jobs) * 2 * 10, "outer fits per reading (reading b: each with 12 inner fits and a refit)")
cb = dec[dec.reading == "b"]
print("Reading (b), C chosen per fold, TRE decoder:", pd.Series(np.concatenate(cb.C_tre.values)).value_counts().sort_index().to_dict(),
      "; XRE decoder:", pd.Series(np.concatenate(cb.C_xre.values)).value_counts().sort_index().to_dict())

64 session-area decoders, both readings, in 1023 s
Sessions per area: {'V1': 12, 'LM': 8, 'RL': 11, 'AL': 10, 'PM': 12, 'AM': 11} (paper 12, 8, 11, 10, 12, 11)
Trials per class: [np.int64(400)]
Fits stopped at max_iter without converging: {'a': 0, 'b': 2} of 1280 outer fits per reading (reading b: each with 12 inner fits and a refit)
Reading (b), C chosen per fold, TRE decoder: {0.001: 41, 0.01: 489, 0.1: 84, 1.0: 26} ; XRE decoder: {0.001: 87, 0.01: 395, 0.1: 89, 1.0: 69}

fig3a-npx-v1-decoding

A linear SVM decodes the four images IC1, LC1, LC2 and IC2 from V1 activity in the Neuropixels data with tenfold cross-validation accuracy above the chance of 0.25.

for rd in ["a", "b"]:
    x = dec[(dec.reading == rd)]
    v = x[x.area == "V1"].acc
    print(f"reading ({rd}) V1: mean accuracy {v.mean():.3f} +/- {v.sem():.3f} (s.e.m.), range {v.min():.3f} to {v.max():.3f}, "
          f"n = {len(v)}, signed-rank against 0.25 p = {il.signrank_p(v - 0.25):.4g}; sessions above 0.25: {(v > 0.25).sum()}")
    if rd == "a":
        print("   other areas:", {a: (round(x[x.area == a].acc.mean(), 3), round(il.signrank_p(x[x.area == a].acc - 0.25), 4))
                                 for a in il.AREAS[1:]})
print(dec[(dec.reading == "a") & (dec.area == "V1")][["subject", "n_units", "acc"]].round(3).to_string(index=False))

reading (a) V1: mean accuracy 0.697 +/- 0.052 (s.e.m.), range 0.268 to 0.887, n = 12, signed-rank against 0.25 p = 0.0004883; sessions above 0.25: 12
   other areas: {'LM': (np.float64(0.709), 0.0078), 'RL': (np.float64(0.543), 0.001), 'AL': (np.float64(0.648), 0.002), 'PM': (np.float64(0.506), 0.0005), 'AM': (np.float64(0.48), 0.001)}
reading (b) V1: mean accuracy 0.719 +/- 0.056 (s.e.m.), range 0.264 to 0.922, n = 12, signed-rank against 0.25 p = 0.0004883; sessions above 0.25: 12
subject  n_units   acc
 619293       89 0.268
 619296      168 0.887
 620333      114 0.682
 620334       68 0.694
 625545      190 0.781
 625554      179 0.781
 625555      314 0.486
 630506       58 0.755
 631510      129 0.821
 631570      142 0.798
 633229      181 0.856
 637484      172 0.553

Result: reproduced. V1 decoding accuracy is 0.697 under reading (a) and 0.719 under reading (b), against the chance of 0.25, with a signed-rank p = 0.00049 under both (12 sessions, every one above 0.25).

fig3c-npx-v1-inference

In Neuropixels V1, TRE trials are decoded as the matching IC image more often than as the matching LC image (0.51 against 0.22), so IC inference is represented in V1.

for rd in ["a", "b"]:
    x = dec[(dec.reading == rd) & (dec.area == "V1")]
    print(f"reading ({rd}): P(TRE -> IC) {x.tre_ic.mean():.3f} +/- {x.tre_ic.sem():.3f}, P(TRE -> LC) {x.tre_lc.mean():.3f} +/- "
          f"{x.tre_lc.sem():.3f}; IC > LC in {(x.tre_ic > x.tre_lc).sum()} of {len(x)} sessions; "
          f"signed-rank IC vs LC p = {il.signrank_p(x.tre_ic, x.tre_lc):.4g}")
print("Paper: 0.51 (tolerance 0.051) and 0.22 (0.022), P < 0.05")
print(dec[(dec.area == "V1")].pivot_table(index="subject", columns="reading", values=["tre_ic", "tre_lc"]).round(3).to_string())

reading (a): P(TRE -> IC) 0.516 +/- 0.062, P(TRE -> LC) 0.197 +/- 0.026; IC > LC in 10 of 12 sessions; signed-rank IC vs LC p = 0.002441
reading (b): P(TRE -> IC) 0.518 +/- 0.067, P(TRE -> LC) 0.210 +/- 0.033; IC > LC in 10 of 12 sessions; signed-rank IC vs LC p = 0.002441
Paper: 0.51 (tolerance 0.051) and 0.22 (0.022), P < 0.05
        tre_ic        tre_lc       
reading      a      b      a      b
subject                            
619293   0.254  0.245  0.257  0.270
619296   0.901  0.935  0.033  0.024
620333   0.408  0.366  0.273  0.318
620334   0.435  0.426  0.139  0.166
625545   0.611  0.672  0.250  0.234
625554   0.791  0.785  0.177  0.183
625555   0.301  0.292  0.217  0.222
630506   0.335  0.322  0.208  0.201
631510   0.610  0.684  0.062  0.023
631570   0.719  0.722  0.143  0.162
633229   0.549  0.482  0.309  0.424
637484   0.280  0.289  0.294  0.293

Result: reproduced. Under reading (b), TRE is decoded as the matching IC with probability 0.518 and as the matching LC with 0.210, both within 0.51 +/- 0.051 and 0.22 +/- 0.022, with IC above LC at signed-rank p = 0.0024. Under reading (a), IC is 0.516 and LC 0.197, just below the lower bound 0.198, with the same p. With two further cross-validation seeds, both values are within tolerance under both readings: reading (a) gives IC 0.513 and 0.511, LC 0.198 and 0.199, p = 0.0024 (decoding_checks_shin_2025.ipynb); reading (b) gives IC 0.519 and 0.514, LC 0.209 and 0.215, p = 0.0015 and 0.0068 (decoding_seed1_shin_2025.ipynb, decoding_seed2_shin_2025.ipynb).

fig3d-npx-areas-inference

In the Neuropixels data almost every visual area decodes TRE as IC more than as LC, significantly in V1, LM, RL, AL and PM and not in AM, with V1, LM and AL the strongest.

paper_p = {"V1": 0.0132, "LM": 0.0156, "RL": 0.0244, "AL": 0.0039, "PM": 0.0342, "AM": 0.1748}
for rd in ["a", "b"]:
    rows = []
    for a in il.AREAS:
        x = dec[(dec.reading == rd) & (dec.area == a)]
        rows.append({"area": a, "sessions": len(x), "IC": x.tre_ic.mean(), "LC": x.tre_lc.mean(), "inference_mean": x.inference.mean(),
                     "inference_sem": x.inference.sem(), "p": il.signrank_p(x.inference), "paper_p": paper_p[a]})
    t = pd.DataFrame(rows)
    t["sig_IC_above_LC"] = (t.p < 0.05) & (t.inference_mean > 0)
    print(f"reading ({rd})")
    print(t.round(4).to_string(index=False))
    print(f"  {int(t.sig_IC_above_LC.sum())} areas significant")
    print("  significant areas:", t[t.sig_IC_above_LC].area.tolist(),
          "; three largest mean scores:", t.nlargest(3, "inference_mean").area.tolist(),
          "; three smallest p:", t.nsmallest(3, "p").area.tolist())
print("Per-session inference scores, reading (a) and (b):")
print(dec.pivot_table(index="subject", columns=["reading", "area"], values="inference").round(3).to_string())

reading (a)
area  sessions     IC     LC  inference_mean  inference_sem      p  paper_p  sig_IC_above_LC
  V1        12 0.5162 0.1966          0.3195         0.0805 0.0024   0.0132             True
  LM         8 0.5371 0.2356          0.3015         0.0798 0.0156   0.0156             True
  RL        11 0.4220 0.3306          0.0914         0.0250 0.0020   0.0244             True
  AL        10 0.5111 0.2526          0.2585         0.0747 0.0020   0.0039             True
  PM        12 0.3690 0.3111          0.0578         0.0255 0.0923   0.0342            False
  AM        11 0.3926 0.3319          0.0607         0.0138 0.0029   0.1748             True
  5 areas significant
  significant areas: ['V1', 'LM', 'RL', 'AL', 'AM'] ; three largest mean scores: ['V1', 'LM', 'AL'] ; three smallest p: ['RL', 'AL', 'V1']
reading (b)
area  sessions     IC     LC  inference_mean  inference_sem      p  paper_p  sig_IC_above_LC
  V1        12 0.5184 0.2100          0.3084         0.0907 0.0024   0.0132             True
  LM         8 0.5654 0.2195          0.3459         0.0859 0.0156   0.0156             True
  RL        11 0.4248 0.3407          0.0841         0.0264 0.0049   0.0244             True
  AL        10 0.5366 0.2513          0.2853         0.0830 0.0020   0.0039             True
  PM        12 0.3810 0.3091          0.0718         0.0347 0.0771   0.0342            False
  AM        11 0.3958 0.3465          0.0494         0.0279 0.0830   0.1748            False
  4 areas significant
  significant areas: ['V1', 'LM', 'RL', 'AL'] ; three largest mean scores: ['LM', 'V1', 'AL'] ; three smallest p: ['AL', 'V1', 'RL']
Per-session inference scores, reading (a) and (b):
reading      a                                         b                                   
area        AL     AM     LM     PM     RL     V1     AL     AM     LM     PM     RL     V1
subject                                                                                    
619293   0.423  0.029    NaN -0.009  0.118 -0.003  0.491 -0.166    NaN -0.125  0.109 -0.025
619296   0.716  0.071    NaN  0.164 -0.009  0.867  0.755  0.072    NaN  0.179 -0.033  0.911
620333   0.136  0.104  0.542 -0.015  0.251  0.135  0.151  0.087  0.550 -0.033  0.253  0.048
620334     NaN  0.055  0.545  0.195  0.228  0.296    NaN -0.008  0.605  0.207  0.225  0.260
625545   0.066  0.060  0.352  0.053  0.046  0.361  0.116  0.092  0.327  0.046  0.035  0.438
625554   0.565  0.115    NaN  0.141  0.048  0.614  0.667  0.146    NaN  0.259  0.089  0.602
625555   0.228  0.009    NaN  0.150  0.054  0.084  0.271  0.024    NaN  0.149  0.032  0.070
630506   0.258    NaN  0.448 -0.075  0.113  0.128  0.217    NaN  0.613 -0.080  0.101  0.121
631510   0.063  0.117  0.110 -0.032  0.052  0.549  0.046  0.172  0.071 -0.033  0.034  0.661
631570   0.014 -0.023 -0.019 -0.005    NaN  0.576  0.002 -0.006 -0.017  0.111    NaN  0.559
633229   0.117  0.097  0.042  0.051  0.009  0.240  0.137  0.108  0.202  0.108  0.010  0.059
637484     NaN  0.032  0.392  0.076  0.094 -0.013    NaN  0.023  0.416  0.076  0.070 -0.003

Result: reproduced, resting on two assumptions. The readings differ on the seed of this notebook. Under reading (a), five areas decode TRE as IC significantly more than as LC (V1, LM, RL, AL and AM; PM p = 0.092), and V1, LM and AL have the three largest mean scores, which meets the criterion; the two further seeds give the same five areas (decoding_checks_shin_2025.ipynb). Under reading (b), adopted, four areas are significant (V1, LM, RL and AL; PM p = 0.077, AM p = 0.083), and LM, V1 and AL have the three largest mean scores, a partial reproduction. Reading (b) is adopted, as the paper's code chooses the SVM hyperparameters by a search, which reading (a) leaves out. With the two further seeds, reading (b) gives five significant areas each, with LM, V1 and AL the three largest: V1, LM, RL, AL and PM with seed +100000 (PM p = 0.0425, AM p = 0.083), and V1, LM, RL, AL and AM with seed +200000 (AM p = 0.0186, PM p = 0.052) (decoding_seed1_shin_2025.ipynb, decoding_seed2_shin_2025.ipynb). The criterion is therefore met in five of the six runs (two readings, three seeds); only the seed of this notebook under reading (b), where PM and AM fall just above 0.05, does not meet it. The outcome follows the majority of the seeds under reading (b), a rule chosen after the first seed's result and not part of the criterion (assumption); a rule requiring every seed would give a partial reproduction. V1, LM, RL and AL are significant and LM, V1 and AL the three largest in all six runs. Which of PM and AM is significant changes with the seed under reading (b); under reading (a), AM is significant and PM is not with every seed, the reverse of the paper.

ed7b-npx-xre-inference

Decoders trained on the XRE real-bar images classify IC1 and IC2 trials as the XRE image whose real bar matches the illusory bar, above chance in all six Neuropixels areas.

XRE inference score (P(IC1 -> XRE1) + P(IC2 -> XRE2))/2, signed-rank against 0.5, the chance of a two-class decoder.

paper_xp = {"V1": 0.0005, "LM": 0.0078, "RL": 0.0049, "AL": 0.0020, "PM": 0.0005, "AM": 0.0137}
for rd in ["a", "b"]:
    rows = []
    for a in il.AREAS:
        x = dec[(dec.reading == rd) & (dec.area == a)]
        rows.append({"area": a, "sessions": len(x), "xre_acc": x.xre_acc.mean(), "xre_score_mean": x.xre_score.mean(),
                     "xre_score_sem": x.xre_score.sem(), "p_vs_0.5": il.signrank_p(x.xre_score - 0.5), "paper_p": paper_xp[a]})
    t = pd.DataFrame(rows)
    t["sig_above_0.5"] = (t["p_vs_0.5"] < 0.05) & (t.xre_score_mean > 0.5)
    print(f"reading ({rd}): areas significant above 0.5: {int(t['sig_above_0.5'].sum())} of 6")
    print(t.round(4).to_string(index=False))

reading (a): areas significant above 0.5: 6 of 6
area  sessions  xre_acc  xre_score_mean  xre_score_sem  p_vs_0.5  paper_p  sig_above_0.5
  V1        12   0.8624          0.7184         0.0328    0.0005   0.0005           True
  LM         8   0.8125          0.7131         0.0364    0.0078   0.0078           True
  RL        11   0.6439          0.6086         0.0268    0.0020   0.0049           True
  AL        10   0.7644          0.6924         0.0397    0.0020   0.0020           True
  PM        12   0.6520          0.5983         0.0253    0.0010   0.0005           True
  AM        11   0.6331          0.5764         0.0251    0.0137   0.0137           True
reading (b): areas significant above 0.5: 6 of 6
area  sessions  xre_acc  xre_score_mean  xre_score_sem  p_vs_0.5  paper_p  sig_above_0.5
  V1        12   0.8811          0.7621         0.0345    0.0005   0.0005           True
  LM         8   0.8270          0.7241         0.0384    0.0078   0.0078           True
  RL        11   0.6530          0.6068         0.0297    0.0049   0.0049           True
  AL        10   0.7749          0.7053         0.0403    0.0020   0.0020           True
  PM        12   0.6632          0.5984         0.0272    0.0010   0.0005           True
  AM        11   0.6274          0.5847         0.0304    0.0137   0.0137           True

Result: reproduced. The XRE inference score is above 0.5 with signed-rank p < 0.05 in all six areas under both readings, with p from 0.0005 to 0.0137, close to the paper's values (0.0005 to 0.0137).

ed7a-npx-no-area-above-v1

No higher visual area shows a substantially larger inference performance than V1, although decoding accuracy and inference scores differ across areas.

Kruskal-Wallis across the six areas of the per-session values, then (a) two-sided signed-rank of V1 against each HVA over the sessions with both, and (b) Tukey-Kramer comparisons after the Kruskal-Wallis. A claim of absence: reproduced when no HVA is significantly larger than V1.

def tukey_vs_v1(vals_by_area):
    # Tukey-Kramer on mean ranks after Kruskal-Wallis (MATLAB multcompare default), HVA minus V1
    allv = np.concatenate([vals_by_area[a] for a in il.AREAS])
    lab = np.concatenate([[a] * len(vals_by_area[a]) for a in il.AREAS])
    r = stats.rankdata(allv)
    N = len(allv)
    _, c = np.unique(allv, return_counts=True)
    var = N * (N + 1) / 12 * (1 - np.sum(c ** 3 - c) / (N ** 3 - N))
    qc = il.tukey_qcrit(6)
    out = {}
    for a in il.AREAS[1:]:
        d = r[lab == a].mean() - r[lab == "V1"].mean()
        se = np.sqrt(var * (1 / np.sum(lab == a) + 1 / np.sum(lab == "V1")) / 2)
        out[a] = (round(d, 1), bool(abs(d) / se > qc))
    return out

for rd in ["a", "b"]:
    x = dec[dec.reading == rd]
    print(f"reading ({rd})")
    for meas, paper in [("acc", 0.0019), ("inference", 0.0249), ("xre_acc", 0.0005), ("xre_score", 0.0022)]:
        vals = {a: x[x.area == a][meas].values for a in il.AREAS}
        kwp = stats.kruskal(*vals.values()).pvalue
        piv = x.pivot(index="subject", columns="area", values=meas)
        pair = {}
        for a in il.AREAS[1:]:
            both = piv[["V1", a]].dropna()
            pair[a] = (round(both[a].mean() - both["V1"].mean(), 3), round(il.signrank_p(both[a], both["V1"]), 4))
        print(f"  {meas}: Kruskal-Wallis p = {kwp:.4f} (paper {paper}); means", {a: round(v.mean(), 3) for a, v in vals.items()})
        print(f"     (a) HVA minus V1 and signed-rank p:", pair)
        print(f"     (b) Tukey-Kramer mean-rank difference HVA minus V1, significant:", tukey_vs_v1(vals))

reading (a)
  acc: Kruskal-Wallis p = 0.0006 (paper 0.0019); means {'V1': np.float64(0.697), 'LM': np.float64(0.709), 'RL': np.float64(0.543), 'AL': np.float64(0.648), 'PM': np.float64(0.506), 'AM': np.float64(0.48)}
     (a) HVA minus V1 and signed-rank p: {'LM': (np.float64(-0.033), 0.8438), 'RL': (np.float64(-0.144), 0.0869), 'AL': (np.float64(-0.063), 0.3223), 'PM': (np.float64(-0.191), 0.0024), 'AM': (np.float64(-0.211), 0.001)}
     (b) Tukey-Kramer mean-rank difference HVA minus V1, significant: {'LM': (np.float64(0.9), False), 'RL': (np.float64(-19.4), False), 'AL': (np.float64(-5.7), False), 'PM': (np.float64(-23.1), True), 'AM': (np.float64(-25.6), True)}

  inference: Kruskal-Wallis p = 0.0071 (paper 0.0249); means {'V1': np.float64(0.32), 'LM': np.float64(0.302), 'RL': np.float64(0.091), 'AL': np.float64(0.258), 'PM': np.float64(0.058), 'AM': np.float64(0.061)}
     (a) HVA minus V1 and signed-rank p: {'LM': (np.float64(0.018), 1.0), 'RL': (np.float64(-0.205), 0.123), 'AL': (np.float64(-0.097), 0.375), 'PM': (np.float64(-0.262), 0.0049), 'AM': (np.float64(-0.276), 0.0098)}
     (b) Tukey-Kramer mean-rank difference HVA minus V1, significant: {'LM': (np.float64(-2.2), False), 'RL': (np.float64(-16.8), False), 'AL': (np.float64(-2.0), False), 'PM': (np.float64(-21.3), False), 'AM': (np.float64(-20.5), False)}
  xre_acc: Kruskal-Wallis p = 0.0005 (paper 0.0005); means {'V1': np.float64(0.862), 'LM': np.float64(0.812), 'RL': np.float64(0.644), 'AL': np.float64(0.764), 'PM': np.float64(0.652), 'AM': np.float64(0.633)}
     (a) HVA minus V1 and signed-rank p: {'LM': (np.float64(-0.102), 0.0781), 'RL': (np.float64(-0.215), 0.0137), 'AL': (np.float64(-0.094), 0.1934), 'PM': (np.float64(-0.21), 0.0015), 'AM': (np.float64(-0.221), 0.0029)}
     (b) Tukey-Kramer mean-rank difference HVA minus V1, significant: {'LM': (np.float64(-5.5), False), 'RL': (np.float64(-26.0), True), 'AL': (np.float64(-9.4), False), 'PM': (np.float64(-25.4), True), 'AM': (np.float64(-26.3), True)}

  xre_score: Kruskal-Wallis p = 0.0077 (paper 0.0022); means {'V1': np.float64(0.718), 'LM': np.float64(0.713), 'RL': np.float64(0.609), 'AL': np.float64(0.692), 'PM': np.float64(0.598), 'AM': np.float64(0.576)}
     (a) HVA minus V1 and signed-rank p: {'LM': (np.float64(-0.023), 0.7422), 'RL': (np.float64(-0.096), 0.1016), 'AL': (np.float64(-0.031), 0.625), 'PM': (np.float64(-0.12), 0.0034), 'AM': (np.float64(-0.148), 0.0049)}
     (b) Tukey-Kramer mean-rank difference HVA minus V1, significant: {'LM': (np.float64(0.5), False), 'RL': (np.float64(-16.1), False), 'AL': (np.float64(-3.7), False), 'PM': (np.float64(-19.3), False), 'AM': (np.float64(-21.8), False)}
reading (b)

  acc: Kruskal-Wallis p = 0.0017 (paper 0.0019); means {'V1': np.float64(0.719), 'LM': np.float64(0.738), 'RL': np.float64(0.556), 'AL': np.float64(0.677), 'PM': np.float64(0.525), 'AM': np.float64(0.498)}
     (a) HVA minus V1 and signed-rank p: {'LM': (np.float64(-0.028), 0.7422), 'RL': (np.float64(-0.152), 0.123), 'AL': (np.float64(-0.059), 0.3223), 'PM': (np.float64(-0.193), 0.0034), 'AM': (np.float64(-0.218), 0.002)}
     (b) Tukey-Kramer mean-rank difference HVA minus V1, significant: {'LM': (np.float64(1.4), False), 'RL': (np.float64(-19.4), False), 'AL': (np.float64(-5.4), False), 'PM': (np.float64(-20.6), False), 'AM': (np.float64(-24.0), True)}
  inference: Kruskal-Wallis p = 0.0156 (paper 0.0249); means {'V1': np.float64(0.308), 'LM': np.float64(0.346), 'RL': np.float64(0.084), 'AL': np.float64(0.285), 'PM': np.float64(0.072), 'AM': np.float64(0.049)}
     (a) HVA minus V1 and signed-rank p: {'LM': (np.float64(0.078), 0.8438), 'RL': (np.float64(-0.201), 0.2061), 'AL': (np.float64(-0.059), 0.7695), 'PM': (np.float64(-0.237), 0.0122), 'AM': (np.float64(-0.276), 0.0186)}
     (b) Tukey-Kramer mean-rank difference HVA minus V1, significant: {'LM': (np.float64(5.8), False), 'RL': (np.float64(-13.2), False), 'AL': (np.float64(3.5), False), 'PM': (np.float64(-13.4), False), 'AM': (np.float64(-16.4), False)}
  xre_acc: Kruskal-Wallis p = 0.0002 (paper 0.0005); means {'V1': np.float64(0.881), 'LM': np.float64(0.827), 'RL': np.float64(0.653), 'AL': np.float64(0.775), 'PM': np.float64(0.663), 'AM': np.float64(0.627)}
     (a) HVA minus V1 and signed-rank p: {'LM': (np.float64(-0.108), 0.0234), 'RL': (np.float64(-0.224), 0.0098), 'AL': (np.float64(-0.1), 0.1934), 'PM': (np.float64(-0.218), 0.001), 'AM': (np.float64(-0.245), 0.002)}
     (b) Tukey-Kramer mean-rank difference HVA minus V1, significant: {'LM': (np.float64(-7.2), False), 'RL': (np.float64(-26.2), True), 'AL': (np.float64(-10.7), False), 'PM': (np.float64(-25.4), True), 'AM': (np.float64(-30.0), True)}
  xre_score: Kruskal-Wallis p = 0.0015 (paper 0.0022); means {'V1': np.float64(0.762), 'LM': np.float64(0.724), 'RL': np.float64(0.607), 'AL': np.float64(0.705), 'PM': np.float64(0.598), 'AM': np.float64(0.585)}
     (a) HVA minus V1 and signed-rank p: {'LM': (np.float64(-0.061), 0.4609), 'RL': (np.float64(-0.143), 0.0244), 'AL': (np.float64(-0.051), 0.2754), 'PM': (np.float64(-0.164), 0.0005), 'AM': (np.float64(-0.178), 0.0029)}
     (b) Tukey-Kramer mean-rank difference HVA minus V1, significant: {'LM': (np.float64(-5.1), False), 'RL': (np.float64(-22.2), True), 'AL': (np.float64(-7.8), False), 'PM': (np.float64(-24.2), True), 'AM': (np.float64(-25.5), True)}

Result: reproduced. No higher visual area has a significantly larger TRE or XRE inference score than V1 under either pairwise test or either decoder reading: every significant difference has the HVA below V1. These are PM and AM in several measures, RL in XRE accuracy (both readings; signed-rank and Tukey-Kramer) and XRE score (reading b; signed-rank p = 0.0244 and Tukey-Kramer), and LM in XRE accuracy under reading (b) (signed-rank p = 0.0234). The Kruskal-Wallis tests across areas are significant, as in the paper (reading a: p = 0.0006, 0.0071, 0.0005 and 0.0077 for accuracy, inference score, XRE accuracy and XRE score; paper 0.0019, 0.0249, 0.0005, 0.0022).

print(f"Notebook ran in {(time.time() - started) / 60:.1f} min")

Notebook ran in 19.3 min

Summary

Claim	Outcome
fig3a-npx-v1-decoding	reproduced
fig3c-npx-v1-inference	reproduced (within tolerance under reading b with three seeds; under reading a with two of three seeds)
fig3d-npx-areas-inference	reproduced (assumptions: decoder reading b, and the majority of three seeds; reading a meets the criterion with every seed, reading b with two of three)
ed7b-npx-xre-inference	reproduced
ed7a-npx-no-area-above-v1	reproduced

A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.

How the analysis differs from the paper's: the decoders use scikit-learn's libsvm SVM in place of MATLAB's fitcecoc. The code's Bayesian hyperparameter search is replaced by a grid search over C from 10^-3 to 1 with one-vs-one coding (reading b); decoding_checks_shin_2025.ipynb shows, on the first outer fold of V1, PM and AM, that C above 1 gave a higher inner accuracy in 1 of 35 sessions (by 0.001, from fits that mostly did not converge) and that the one-vs-all coding changed the best inner accuracy by -0.026 to 0.027. Extended Data Fig. 7c pools Neuropixels and 2p areas; its Kruskal-Wallis values are compared here on the six Neuropixels areas only.