Schuck, Cai, Wilson & Niv (2016) on CRCNS data set ofc-4¶
Paper: Schuck NW, Cai MB, Wilson RC, Niv Y (2016). Human orbitofrontal cortex represents a cognitive
map of state space. Neuron 91:1402-1412, doi:10.1016/j.neuron.2016.08.019. The claims file
(claims/crcns-ofc-4/10.1016-j.neuron.2016.08.019.yaml) was written from the version of record with
its Supplemental Information.
Data set: CRCNS ofc-4, "fMRI time series in orbitofrontal cortex during face-house state-space task"
(Schuck, Cai, Wilson & Niv 2019), doi:10.6080/K0FN14CH, read from the public S3 bucket
crcnsarchive under ofc-4/live/. Version: unversioned. CRCNS data sets carry no version numbers
(docs/archives.md); the bucket holds one copy, under live/ (no old/ folder), described by version
0.6 of the data descriptor (April 19, 2019), the only version. The paper cites no version of the data
set, which was published three years after it, so this copy is the one to pin. The notebook prints
the data set's checksums.md5 and the MD5 of the zip as streamed; the files read are FHSS_OFC.zip
(MD5 4372953c3900c5a965ad09526a9510ed) and checksums.md5, so a later change of the zip is detectable.
The data set holds, for the 24 participants of the paper's trial-wise analysis (the paper analysed 27), the OFC time series in native space with no voxel coordinates, the design matrices (16 state columns, 'first', 'err') with drift and motion nuisance time courses, and the event files (onset, reaction time, stimulus duration per state, 'first' and 'err'). It holds no whole-brain data, no other ROI, no training session, no response side or response mapping, and no voxel coordinates.
| Claim | Tested here |
|---|---|
| behavior-error-rate | tested |
| supp-errors-by-trial-class | tested |
| supp-errors-stable-over-blocks | tested |
| supp-errors-not-by-age-category | tested |
| supp-training-errors-higher | blocked: data_missing (no training session) |
| supp-training-rt-unchanged | blocked: data_missing (no training session) |
| supp-rt-stable-over-blocks | tested |
| supp-rt-faces-faster | tested |
| supp-rt-by-trial-class | tested |
| fig2a-16way-decoding | tested |
| fig2-two-back-state-space | tested |
| fig2b-hidden-components | tested |
| fig2b-current-age-absent | tested |
| fig2b-two-back-absent | tested |
| figs2c-8way-decoding | tested |
| figs2a-permutation-control | tested |
| synthetic-data-control | blocked: metadata_insufficient (no voxel coordinates) |
| fig2c-enter-cluster | tested |
| fig2d-enter-nonenter-decoding | tested |
| fig3-searchlight-conjunction | blocked: data_missing (no whole-brain data) |
| other-rois-incomplete | blocked: data_missing (no other ROIs) |
| fig4a-decoding-error-correlation | tested |
| fig4b-error-trial-drop | tested |
| fig4b-decline-before-errors | tested |
| fig4b-error-vs-matched | tested |
| fig5-similarity-predicts-errors | tested |
| figs5-euclidean-distance | tested |
| fig5-synthetic-control | blocked: metadata_insufficient (depends on synthetic-data-control) |
| supp-motor-action-not-in-ofc | blocked: data_missing (no response side; no whole-brain data) |
Method details taken from the paper, readings and assumptions are given in each section and in
ofc4.py. One assumption runs through the decoding claims: the beta maps are z-scored per voxel
across the 16 state maps of a run ('voxel'), the operation of the Princeton MVPA toolbox's
zscore_runs, which the Methods name ("Data Preprocessing"); the other reading of O2, each map
across voxels ('map'), is shown everywhere beside it. The notebook was revised after the first
verification round: the readings of O3 (post-error trials), a composition of the error blocks, a
class-matched comparison of error trials, a control for trial counts in the distance analysis and
a by-participant slope variant of the RT model were added. No map or volume is smoothed (O1): the files
hold no voxel coordinates.
import hashlib
import time
import warnings
import fsspec
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy import stats
import ofc4
from statsmodels.tools.sm_exceptions import ConvergenceWarning
# Mixed-model fits print whether they converged; the optimiser's own warnings are silenced.
warnings.filterwarnings("ignore", category=ConvergenceWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning, module="statsmodels")
warnings.filterwarnings("ignore", category=UserWarning, module="statsmodels")
pd.set_option("display.width", 160)
pd.set_option("display.max_columns", 30)
T0 = time.time()
DATASET = "crcns:ofc-4/unversioned" # doi:10.6080/K0FN14CH
# Files and checksums of the data set, and the zip's MD5 computed from the stream.
print(ofc4.read_text(f"{ofc4.BASE}/checksums.md5"))
md5 = hashlib.md5()
with fsspec.open(ofc4.ZIP_URL, "rb", block_size=2**22) as f:
while chunk := f.read(2**22):
md5.update(chunk)
print("MD5 of FHSS_OFC.zip as streamed:", md5.hexdigest())
names = ofc4.members()
print(len(names), "zip members")
SUBJ = ofc4.subjects()
print("participants:", len(SUBJ), SUBJ)
print("runs per participant:", {s: len(ofc4.runs(s)) for s in SUBJ})
33f768f9516d6fc66927c38d50236b2c docs/crcns_ofc-4_data_description.docx aa90a399744ca98a66acce4f6b73a53d docs/crcns_ofc-4_data_description.pdf 4372953c3900c5a965ad09526a9510ed FHSS_OFC.zip 279ea8d1457e645c45b40b9c1c66c263 filelist.txt
MD5 of FHSS_OFC.zip as streamed: 4372953c3900c5a965ad09526a9510ed
4454 zip members
participants: 24 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]
runs per participant: {1: 4, 2: 5, 3: 5, 4: 5, 5: 5, 6: 5, 7: 5, 8: 5, 9: 5, 10: 5, 11: 5, 12: 5, 13: 5, 14: 5, 15: 5, 16: 5, 17: 5, 18: 5, 19: 5, 20: 5, 21: 5, 22: 5, 23: 5, 24: 5}
The zip's MD5 as streamed equals the one in checksums.md5. The data set holds 24 participants, one
with four runs and 23 with five, as the descriptor says.
# Every field of every run file: shapes, labels, TR, and checks of the stated contents.
rows = []
for s in SUBJ:
for r in ofc4.runs(s):
d = ofc4.load_run(s, r)
rows.append({"subject": s, "run": r, "n_TR": d["n_TR"], "Y_rows": d["Y"].shape[0], "n_voxel": d["n_voxel"],
"Y_cols": d["Y"].shape[1], "X_cols": d["X"].shape[1], "TR": d["TR"],
"labels": ",".join(l for l in d["labels"] if l not in ofc4.STATES),
"nuis_cols": d["nuis"].shape[1], "nuis0_const": bool(np.allclose(d["nuis"][:, 0], 1)),
"Y_min": d["Y"].min(), "Y_max": d["Y"].max(), "Y_mean": d["Y"].mean(),
"n_zero_vox": int((d["Y"].std(0) == 0).sum()),
"states_present": sum(l in ofc4.STATES for l in d["labels"]), "keys": ",".join(d["keys"]),
"md5_Y": hashlib.md5(d["Y"].tobytes()).hexdigest(),
"md5_Y_head": hashlib.md5(np.round(d["Y"][:20, :50], 3).tobytes()).hexdigest()})
files = pd.DataFrame(rows)
print("fields of the run files (all variables in each .mat file):", files["keys"].value_counts().to_dict())
print(files.drop(columns=["md5_Y", "md5_Y_head", "keys"]).describe().round(2).to_string())
print("non-state columns:", files.labels.value_counts().to_dict())
print("TR values:", files.TR.unique(), "; n_TR equals rows:", (files.n_TR == files.Y_rows).all(),
"; n_voxel equals columns:", (files.n_voxel == files.Y_cols).all())
print("distinct time series:", files.md5_Y.nunique(), "of", len(files), "; distinct first 20 x 50 blocks:", files.md5_Y_head.nunique())
print("voxels constant across runs within participant:", files.groupby("subject").n_voxel.nunique().eq(1).all())
print("voxels per participant: min", files.n_voxel.min(), "max", files.n_voxel.max())
print("all 16 states in every run:", (files.states_present == 16).all(), "; first nuisance column constant:", files.nuis0_const.all())
print("runs without a 'first' column:", files.loc[~files.labels.str.contains("first"), ["subject", "run"]].values.tolist())
print("runs without an 'err' column:", files.loc[~files.labels.str.contains("err"), ["subject", "run"]].values.tolist())
fields of the run files (all variables in each .mat file): {'OFC_timecourse,TR,condition_label,design_matrix_task,n_TR,n_voxel,nuisance_timecourse': 119}
subject run n_TR Y_rows n_voxel Y_cols X_cols TR nuis_cols Y_min Y_max Y_mean n_zero_vox states_present
count 119.0 119.00 119.00 119.00 119.00 119.00 119.00 119.0 119.0 119.00 119.00 119.00 119.0 119.0
mean 12.6 2.98 191.13 191.13 4107.90 4107.90 17.71 2.4 10.0 381.82 1629.82 793.64 0.0 16.0
std 6.9 1.41 8.73 8.73 331.86 331.86 0.45 0.0 0.0 35.67 132.43 56.53 0.0 0.0
min 1.0 1.00 168.00 168.00 3450.00 3450.00 17.00 2.4 10.0 278.00 1400.00 655.27 0.0 16.0
25% 7.0 2.00 186.00 186.00 3823.00 3823.00 17.00 2.4 10.0 356.50 1519.50 763.17 0.0 16.0
50% 13.0 3.00 189.00 189.00 4167.00 4167.00 18.00 2.4 10.0 391.00 1622.00 798.98 0.0 16.0
75% 18.5 4.00 194.00 194.00 4366.00 4366.00 18.00 2.4 10.0 408.00 1729.00 840.04 0.0 16.0
max 24.0 5.00 222.00 222.00 4635.00 4635.00 18.00 2.4 10.0 452.00 1954.00 903.78 0.0 16.0
non-state columns: {'first,err': 85, 'first': 31, 'err': 3}
TR values: [2.4] ; n_TR equals rows: True ; n_voxel equals columns: True
distinct time series: 119 of 119 ; distinct first 20 x 50 blocks: 119
voxels constant across runs within participant: True
voxels per participant: min 3450 max 4635
all 16 states in every run: True ; first nuisance column constant: True
runs without a 'first' column: [[11, 3], [18, 4], [18, 5]]
runs without an 'err' column: [[3, 1], [3, 3], [3, 4], [3, 5], [4, 1], [4, 3], [4, 4], [4, 5], [5, 3], [5, 5], [6, 5], [7, 3], [10, 1], [10, 5], [12, 5], [15, 1], [16, 1], [16, 4], [16, 5], [17, 5], [19, 1], [19, 2], [19, 4], [20, 5], [21, 5], [23, 2], [23, 4], [23, 5], [24, 3], [24, 4], [24, 5]]
Checks of the descriptor against the files. Every run file holds the same seven variables
(OFC_timecourse, TR, condition_label, design_matrix_task, n_TR, n_voxel,
nuisance_timecourse) and no mask or voxel coordinates. TR is 2.4 s in every run; n_TR and n_voxel match the
array shapes; each participant's voxel count is the same in every run (3450 to 4635 voxels), consistent
with one native-space mask per participant, though the files give no coordinates to show it. The 119
time series are all distinct, so no run was deposited twice. All 16 state columns are present in every
run. The descriptor says the 'err' column is missing where a run had no error (31 runs here, all
without error events); it does not say that the 'first' column is missing in three runs (Subj-11 run 3,
Subj-18 runs 4 and 5). The first nuisance column is constant, so the GLM's run mean is in the nuisance
set. Intensities are raw scanner values (278 to 1954).
# Events: one row per labelled trial, 'first' trial and error block.
EV = pd.concat([ofc4.events(s, r) for s in SUBJ for r in ofc4.runs(s)], ignore_index=True)
EV["is_state"] = EV.cond.isin(ofc4.STATES)
per_run = EV.groupby(["subject", "run"]).agg(labelled=("is_state", "sum"),
first=("cond", lambda c: (c == "first").sum()),
err_blocks=("cond", lambda c: (c == "err").sum()),
last_onset=("onset", "max")).reset_index()
print(per_run.describe().round(2).to_string())
print("RT range (s):", EV.loc[EV.is_state, "rt"].min(), EV.loc[EV.is_state, "rt"].max(),
"; RT values are multiples of 0.1 s:", np.allclose(EV.rt * 10, np.round(EV.rt * 10)))
print("stimulus duration range (s):", EV.loc[EV.is_state, "dur"].min(), EV.loc[EV.is_state, "dur"].max())
print("onsets equal in RT and dur files:", np.allclose(EV.onset, EV.onset_dur))
# Labelled trials per run plus trials taken up by error blocks: a run of the design has 96 trials after the first.
per_run["designed_in_blocks"] = 96 - per_run.labelled
print("runs with no error block and 96 labelled trials:", ((per_run.err_blocks == 0) & (per_run.labelled == 96)).sum(),
"of", (per_run.err_blocks == 0).sum(), "runs with no error block")
# Consecutive labelled trials must chain: the previous components of a trial are the current ones of the trial before.
bad = []
for (s, r), g in EV.groupby(["subject", "run"]):
c = g.cond.tolist()
for a, b in zip(c, c[1:]):
if a in ofc4.STATES and b in ofc4.STATES and b[2:] != a[:2]:
bad.append((s, r, a, b))
print("adjacent labelled trials that do not chain:", bad)
cls = EV[EV.is_state].assign(klass=lambda d: d.cond.map(ofc4.trial_class))
print("labelled trials per class and run (design: 32 Enter, 32 Internal, 32 Exit per run of 96):")
print(cls.groupby(["subject", "run"]).klass.value_counts().unstack().describe().round(1).to_string())
subject run labelled first err_blocks last_onset count 119.0 119.00 119.00 119.00 119.00 119.00 mean 12.6 2.98 93.03 0.97 1.92 450.83 std 6.9 1.41 3.34 0.16 1.95 19.96 min 1.0 1.00 79.00 0.00 0.00 397.10 25% 7.0 2.00 92.00 1.00 0.00 439.55 50% 13.0 3.00 94.00 1.00 1.00 446.50 75% 18.5 4.00 96.00 1.00 3.00 459.50 max 24.0 5.00 96.00 1.00 10.00 520.90 RT range (s): 0.4 2.7 ; RT values are multiples of 0.1 s: True stimulus duration range (s): 2.8 5.0 onsets equal in RT and dur files: True runs with no error block and 96 labelled trials: 31 of 31 runs with no error block adjacent labelled trials that do not chain: [(1, 4, 'HoHy', 'HoHy')] labelled trials per class and run (design: 32 Enter, 32 Internal, 32 Exit per run of 96): klass Enter Exit Internal count 119.0 119.0 119.0 mean 30.4 31.4 31.2 std 1.8 0.9 1.3 min 23.0 29.0 23.0 25% 30.0 31.0 31.0 50% 31.0 32.0 32.0 75% 32.0 32.0 32.0 max 32.0 32.0 32.0
Every run without an error block has 96 labelled trials after the first, the 96 of the design (Methods, Design), and labelled trials per class are close to the 32 Enter, 32 Internal and 32 Exit per run the design implies. Labels read as current then previous trial chain in every pair of adjacent labelled trials but one: in Subj-01 run 4 a 'HoHy' trial follows a 'HoHy' trial, which the task rules do not allow. RTs are given to 0.1 s (0.4 to 2.7 s), and stimulus durations run from 2.8 to 5.0 s, inside the 2.75 to 5 s of the Methods.
# Error blocks: span from onset to the end of the last stimulus, by the class of the trial before.
ERR = pd.concat([ofc4.error_trials(g) for _, g in EV.groupby(["subject", "run"])], ignore_index=True)
ERR["after_exit"] = ERR.P.map(lambda x: x in ofc4.STATES and ofc4.trial_class(x) == "Exit")
print("error blocks:", len(ERR), "; after an Exit trial:", int(ERR.after_exit.sum()), "; on a run's first trial:", int(ERR.P.isna().sum()))
print("span (s) after an Exit trial:", ERR.loc[ERR.after_exit, "span"].describe().round(1).to_dict())
print("span (s) after another trial:", ERR.loc[~ERR.after_exit, "span"].describe().round(1).to_dict())
print("blocks in runs whose labelled trials and blocks add up to 96:", int(ERR.one_per_block.sum()),
"; of these after an Exit trial:", int((ERR.one_per_block & ERR.after_exit).sum()),
"; reconstructed states legal:", ERR.legal.value_counts().to_dict())
print("error trials with a state, strict reading:", int(ERR.state.notna().sum()), "; span reading:", int(ERR.state_span.notna().sum()),
"; the two disagree on:", int((ERR.state.notna() & ERR.state_span.notna() & (ERR.state != ERR.state_span)).sum()))
print("class known (Enter after Exit, or reconstructed):", int(ERR.klass.notna().sum()), "of", len(ERR))
print("runs whose labelled trials and blocks add up to 96, among runs with blocks:",
int((per_run.err_blocks.gt(0) & (per_run.labelled + per_run.err_blocks == 96)).sum()), "of", int(per_run.err_blocks.gt(0).sum()))
# designed trials taken up per block, in runs whose blocks all follow an Exit trial versus others
per_run2 = per_run.merge(ERR.groupby(["subject", "run"]).after_exit.all().rename("all_after_exit").reset_index(), on=["subject", "run"])
per_run2["extra_per_block"] = (96 - per_run2.labelled - per_run2.err_blocks) / per_run2.err_blocks
print("designed trials beyond one per block, per block: runs with only Enter errors",
per_run2.loc[per_run2.all_after_exit, "extra_per_block"].mean().round(2),
"; runs with other errors", per_run2.loc[~per_run2.all_after_exit, "extra_per_block"].mean().round(2))
# Composition of the blocks. Trial slots per block: time from the block's onset to the next event,
# over the run's mean interval between consecutive labelled trials.
slots = []
for (s, r), g in EV.groupby(["subject", "run"]):
g = g.reset_index(drop=True)
iti = np.mean([g.onset[i] - g.onset[i - 1] for i in range(1, len(g)) if g.is_state[i] and (g.is_state[i - 1] or g.cond[i - 1] == "first")])
for i in np.flatnonzero(g.cond.values == "err"):
P = g.cond[i - 1] if i > 0 else None
slots.append({"before": ofc4.trial_class(P) if P in ofc4.STATES else "none", "rt": g.rt[i], "span": g.dur[i],
"slots": (g.onset[i + 1] - g.onset[i]) / iti if i + 1 < len(g) else np.nan})
slots = pd.DataFrame(slots)
print("trial slots per block (onset to next event over the mean trial interval), by the class of the trial before the block:")
print(slots.groupby("before").slots.describe().round(2).to_string())
print("rounded slots per block:", slots.groupby("before").slots.apply(lambda x: x.dropna().round().astype(int).value_counts().sort_index().to_dict()).dropna().to_dict())
print("'err' durations: RT file", slots.rt.min(), "to", slots.rt.max(), "s; dur file (span)", slots.span.min(), "to", slots.span.max(),
"s; span less RT duration, median", round((slots.span - slots.rt).median(), 2), "s")
# Designed trials inside the blocks, per block, by least squares over runs: 96 less the labelled trials
# against the counts of blocks after an Exit trial and of other blocks.
cnt = ERR.assign(k=np.where(ERR.after_exit, "after_exit", "other")).groupby(["subject", "run", "k"]).size().unstack(fill_value=0)
cnt = cnt.join(per_run.set_index(["subject", "run"]).labelled)
coef = np.linalg.lstsq(cnt[["after_exit", "other"]].values, 96 - cnt.labelled.values, rcond=None)[0]
print(f"designed trials per block, least squares over {len(cnt)} runs with blocks: after an Exit trial {coef[0]:.2f}, other blocks {coef[1]:.2f}")
print("share of blocks after an Exit trial:", round(ERR.after_exit.mean(), 3),
"; Enter share of errors implied by the paper's class rates (3.4, 2.3, 1.4% over equal trial counts):", round(3.4 / (3.4 + 2.3 + 1.4), 3))
fig, ax = plt.subplots(figsize=(6, 3))
ax.hist([ERR.loc[ERR.after_exit, "span"], ERR.loc[~ERR.after_exit, "span"]], bins=np.arange(5, 46, 1),
stacked=True, label=["after Exit", "after Enter or Internal"])
ax.set_xlabel("error block span (s)"); ax.set_ylabel("blocks"); ax.legend()
print("span range plotted 5 to 46 s; outside:", int(((ERR.span < 5) | (ERR.span >= 46)).sum()))
plt.show()
error blocks: 228 ; after an Exit trial: 181 ; on a run's first trial: 3
span (s) after an Exit trial: {'count': 181.0, 'mean': 13.3, 'std': 6.1, 'min': 7.0, '25%': 8.5, '50%': 10.1, '75%': 17.4, 'max': 41.0}
span (s) after another trial: {'count': 47.0, 'mean': 18.5, 'std': 5.6, 'min': 6.8, '25%': 15.7, '50%': 16.8, '75%': 19.2, 'max': 36.6}
blocks in runs whose labelled trials and blocks add up to 96: 40 ; of these after an Exit trial: 40 ; reconstructed states legal: {True: 39}
error trials with a state, strict reading: 39 ; span reading: 98 ; the two disagree on: 0
class known (Enter after Exit, or reconstructed): 181 of 228
runs whose labelled trials and blocks add up to 96, among runs with blocks: 27 of 88
designed trials beyond one per block, per block: runs with only Enter errors 0.38 ; runs with other errors 0.69
trial slots per block (onset to next event over the mean trial interval), by the class of the trial before the block:
count mean std min 25% 50% 75% max
before
Enter 22.0 4.66 1.50 0.98 4.00 4.41 4.92 8.18
Exit 172.0 3.29 1.43 1.70 2.16 2.60 4.24 9.33
Internal 21.0 4.23 1.07 3.63 3.75 3.99 4.18 8.71
none 4.0 3.79 1.70 1.90 2.69 3.78 4.89 5.72
rounded slots per block: {('Enter', 1): 1.0, ('Enter', 4): 12.0, ('Enter', 5): 5.0, ('Enter', 7): 2.0, ('Enter', 8): 2.0, ('Exit', 4): 37.0, ('Exit', 5): 21.0, ('Exit', 7): 1.0, ('Exit', 8): 3.0, ('Exit', 2): 82.0, ('Exit', 3): 20.0, ('Exit', 6): 7.0, ('Exit', 9): 1.0, ('Internal', 4): 18.0, ('Internal', 5): 2.0, ('Internal', 9): 1.0, ('none', 5): 1.0, ('none', 2): 1.0, ('none', 3): 1.0, ('none', 6): 1.0}
'err' durations: RT file 4.2 to 37.3 s; dur file (span) 6.8 to 41.0 s; span less RT duration, median 2.1 s
designed trials per block, least squares over 88 runs with blocks: after an Exit trial 1.52, other blocks 1.87
share of blocks after an Exit trial: 0.794 ; Enter share of errors implied by the paper's class rates (3.4, 2.3, 1.4% over equal trial counts): 0.479
span range plotted 5 to 46 s; outside: 0
The 'err' events are blocks, not single trials: their spans (onset to end of the last stimulus) run from 6.8 to 41 s, two to several trials. All 181 blocks that follow an Exit trial hold a trial of the other category, so their first trial is an Enter error, and the short ones (about two trials) hold the error and its repetition. In runs whose labelled trials and blocks add up to 96 (27 of 88 runs with blocks), every block holds one designed trial and the error trial's state follows from its neighbours; all 39 such states are legal, and they are all Enter errors. Elsewhere, blocks take up more designed trials than one (0.38 more per block in runs with only Enter errors, 0.69 in others), so a block can hold more than one error or the trials after it. The number of errors and the state of most error trials therefore cannot be read from the files: the class is known for 181 of 228 blocks, the full state for 39 (strict reading) or 98 (span reading).
Composition of the blocks. Measured in trial slots (onset to the next event over the run's mean trial interval), 82 of the 172 blocks after an Exit trial span about two slots, the error and its repetition, and most others three to five; blocks after Internal or Enter trials span about four (median 4.0 and 4.4), as an Exit error followed by the repeated preceding trial and the repeated Exit trial would. Over runs, each block after an Exit trial takes up 1.52 designed trials and each other block 1.87 (least squares), so a block after an Exit trial holds on average about half a designed trial (0.52) besides its error trial, and another block close to one (0.87). Whether that trial is a second error or a correct trial after the error the files do not say. The 'err' RT-file duration runs to the last response of the block (span less RT duration, median 2.1 s). Blocks after an Exit trial are 79% of all blocks, where the paper's class rates imply that Enter errors are 48% of errors: Internal and Exit errors appear undercounted here, or merged into blocks.
# Timing: the data set's state columns against regressors rebuilt from the event files (SPM HRF),
# for several shifts of the event times; the shift with the best match is the one the helper uses.
chk = []
for s in SUBJ[:6]:
for r in ofc4.runs(s)[:2]:
d = ofc4.load_run(s, r)
ev = EV[(EV.subject == s) & (EV.run == r)]
for shift in (-1.2, -0.6, 0.0, 0.6):
for j, lab in enumerate(d["labels"]):
sel = ev[ev.cond == lab]
if len(sel) == 0:
continue
x = ofc4.regressor(sel.onset, sel.rt, d["n_TR"], shift=shift)
chk.append({"subject": s, "run": r, "shift": shift, "col": lab, "r": np.corrcoef(x, d["X"][:, j])[0, 1]})
chk = pd.DataFrame(chk)
chk["kind"] = np.where(chk.col.isin(ofc4.STATES), "state", chk.col)
print(chk.groupby(["kind", "shift"]).r.agg(["median", "min"]).round(3).to_string())
# The 'err' column against the err durations of the dur files too
errchk = []
for s in SUBJ[:6]:
for r in ofc4.runs(s)[:2]:
d = ofc4.load_run(s, r)
if "err" not in d["labels"]:
continue
ev = EV[(EV.subject == s) & (EV.run == r) & (EV.cond == "err")]
j = d["labels"].index("err")
errchk.append({"subject": s, "run": r,
"r_rt": np.corrcoef(ofc4.regressor(ev.onset, ev.rt, d["n_TR"]), d["X"][:, j])[0, 1],
"r_dur": np.corrcoef(ofc4.regressor(ev.onset, ev.dur, d["n_TR"]), d["X"][:, j])[0, 1],
"r_impulse": np.corrcoef(ofc4.regressor(ev.onset, np.full(len(ev), 0.1), d["n_TR"]), d["X"][:, j])[0, 1],
"err_col_max": d["X"][:, j].max(), "state_col_max": d["X"][:, :16].max()})
print(pd.DataFrame(errchk).round(3).to_string())
median min
kind shift
err -1.2 0.533 0.385
-0.6 0.460 0.322
0.0 0.381 0.257
0.6 0.302 0.190
first -1.2 0.972 0.949
-0.6 0.999 0.994
0.0 0.986 0.953
0.6 0.931 0.874
state -1.2 0.975 0.951
-0.6 0.997 0.988
0.0 0.984 0.946
0.6 0.929 0.871
subject run r_rt r_dur r_impulse err_col_max state_col_max 0 1 1 0.350 0.307 0.958 3.373 0.348 1 1 2 0.368 0.332 0.980 2.737 0.384 2 2 1 0.424 0.349 0.935 2.700 0.377 3 2 2 0.403 0.355 0.942 2.534 0.408 4 3 2 0.728 0.610 0.983 1.153 0.288 5 4 2 0.730 0.629 0.980 1.182 0.326 6 5 1 0.690 0.538 0.977 1.314 0.386 7 5 2 0.322 0.298 0.833 4.338 0.456 8 6 1 0.495 0.428 0.930 2.136 0.277 9 6 2 0.760 0.589 0.979 1.130 0.350
The data set's state and 'first' columns match regressors rebuilt from the RT files with SPM's canonical HRF (median r 0.997 at a shift of -0.6 s, minimum 0.988), so the design matrices and event files share one clock. The 'err' column does not match the error blocks' RT or stimulus durations (r 0.30 to 0.76) but matches a short event at each block onset (r 0.83 to 0.98): it is not the boxcar from stimulus onset to button press that the descriptor describes for the design matrices.
# Behavioural tables. Designed trials: each labelled trial (correct) and each error block's error
# trial (one error per block, a lower bound when a block holds more). The first trial of a run is
# kept apart. Runs per participant from the files.
NRUN = EV.groupby("subject").run.nunique()
lab = EV[EV.is_state].copy()
lab["klass"] = lab.cond.map(ofc4.trial_class)
lab["cur_cat"] = lab.cond.str[0]
lab["cur_age"] = lab.cond.str[1]
lab["error"] = 0
errt = ERR.copy()
errt["error"] = 1
errt["first_trial"] = errt.P.isna()
trials = pd.concat([lab[["subject", "run", "onset", "klass", "cur_cat", "error"]].assign(cur_age=lab.cur_age, state=lab.cond),
errt.loc[~errt.first_trial, ["subject", "run", "onset", "klass", "cur_cat", "error", "state", "state_span"]]
.assign(cur_age=lambda d: d.state_span.str[1])], ignore_index=True)
print("designed trials after the first:", len(trials), "; errors:", trials.error.sum(),
"; error blocks on a run's first trial:", int(errt.first_trial.sum()))
print("errors with class known:", trials.loc[trials.error == 1, "klass"].notna().sum(),
"; with current age known (span reading):", trials.loc[trials.error == 1, "cur_age"].notna().sum())
def wald(fit, names):
"""Wald chi-square of the named coefficients of a fitted statsmodels model."""
b = fit.params[names].values
V = fit.cov_params().loc[names, names].values
chi2 = float(b @ np.linalg.solve(V, b))
return chi2, len(names), float(stats.chi2.sf(chi2, len(names)))
designed trials after the first: 11296 ; errors: 225 ; error blocks on a run's first trial: 3 errors with class known: 181 ; with current age known (span reading): 98
behavior-error-rate¶
Participants performed the task nearly without errors, with a mean error rate of 2.3% in the scanner.
Errors are the error blocks of the event files ('err'); the files do not tell time-outs from wrong responses, so the two readings of the criterion (errors only, errors plus time-outs) cannot be separated, and both are this count. A block that holds more than one error counts once. Denominators: (a) all trials including the repetitions that follow errors (97 per run plus one repetition per error, Methods, Task), (b) all trials but the first of each run (96 per run, with errors on a first trial left out).
nerr_all = ERR.groupby("subject").size().reindex(SUBJ, fill_value=0)
nerr_nofirst = ERR[ERR.P.notna()].groupby("subject").size().reindex(SUBJ, fill_value=0)
rate_a = 100 * nerr_all / (97 * NRUN + nerr_all)
rate_b = 100 * nerr_nofirst / (96 * NRUN)
er = pd.DataFrame({"runs": NRUN, "error_blocks": nerr_all, "rate_a": rate_a, "rate_b": rate_b})
print(er.round(2).to_string())
for k in ("rate_a", "rate_b"):
print(f"{k}: mean {er[k].mean():.2f}% (SEM {er[k].sem():.2f}), median {er[k].median():.2f}, range {er[k].min():.2f} to {er[k].max():.2f}; "
f"expected 2.3 +/- 0.23")
print("bonus threshold 2.5%: participants below it (reading b):", int((er.rate_b < 2.5).sum()), "of", len(er))
# A block can hold more than one error: the designed trials it takes up (96 per run less the labelled
# trials) bound the number of errors from above, if every designed trial inside a block was an error.
up = per_run.assign(designed_in_blocks=96 - per_run.labelled).groupby("subject").designed_in_blocks.sum().reindex(SUBJ)
er["rate_upper"] = 100 * up / (96 * NRUN)
print(f"upper bound, reading b denominator: mean {er.rate_upper.mean():.2f}% (designed trials inside blocks: {int(up.sum())}, blocks: {int(nerr_all.sum())})")
ERRRATE = er.rate_b / 100
runs error_blocks rate_a rate_b subject 1 4 11 2.76 2.86 2 5 18 3.58 3.75 3 5 1 0.21 0.21 4 5 1 0.21 0.21 5 5 8 1.62 1.67 6 5 7 1.42 1.46 7 5 4 0.82 0.83 8 5 10 2.02 2.08 9 5 36 6.91 7.50 10 5 4 0.82 0.83 11 5 15 3.00 2.92 12 5 14 2.81 2.92 13 5 17 3.39 3.54 14 5 24 4.72 5.00 15 5 5 1.02 1.04 16 5 3 0.61 0.62 17 5 6 1.22 1.25 18 5 12 2.41 2.08 19 5 2 0.41 0.42 20 5 6 1.22 1.25 21 5 8 1.62 1.67 22 5 10 2.02 2.08 23 5 2 0.41 0.42 24 5 4 0.82 0.83 rate_a: mean 1.92% (SEM 0.33), median 1.52, range 0.21 to 6.91; expected 2.3 +/- 0.23 rate_b: mean 1.98% (SEM 0.35), median 1.56, range 0.21 to 7.50; expected 2.3 +/- 0.23 bonus threshold 2.5%: participants below it (reading b): 17 of 24 upper bound, reading b denominator: mean 3.10% (designed trials inside blocks: 353, blocks: 228)
Result: blocked, metadata insufficient. The files do not record the number of errors. Counting one error per block gives a mean of 1.92% (reading a) or 1.98% (reading b), below the 2.07 to 2.53% of the tolerance; counting every designed trial inside a block as an error gives an upper bound of 3.10%, above it. The paper's 2.3% lies between these bounds. Time-outs cannot be told from errors, so the two readings of the criterion give the same count.
supp-errors-by-trial-class¶
Error rates differed by trial class (Enter 3.4%, Exit 2.3%, Internal 1.4%).
An error trial's class is Enter when the trial before its block is an Exit trial. Otherwise it is
Internal or Exit, and the data tell which only where the state is reconstructed (ofc4.error_trials).
The errors of unknown class are shown two ways, all Internal and all Exit, which bound each rate. Per participant, the
rate in a class is its errors over its designed trials. Model readings: (a) linear mixed model of
per-participant class rates, random intercept; (b) logistic model of trial-wise errors; statsmodels
has no frequentist logistic mixed model, so (b) uses participant fixed effects in its place.
unk = trials[(trials.error == 1) & trials.klass.isna()]
print("errors of unknown class (Internal or Exit):", len(unk), "of", int(trials.error.sum()))
cls_res = {}
for assign in ("Internal", "Exit"):
t = trials.copy()
t.loc[t.klass.isna(), "klass"] = assign
pr = t.groupby(["subject", "klass"]).error.mean().mul(100).rename("rate").reset_index()
means = pr.groupby("klass").rate.mean()
m = smf.mixedlm("rate ~ C(klass, Treatment('Internal'))", pr, groups=pr["subject"]).fit()
names = [n for n in m.params.index if n.startswith("C(klass")]
chi2, df, p = wald(m, names)
lg = smf.logit("error ~ C(klass, Treatment('Internal')) + C(subject)", t).fit(disp=0)
chi2b, dfb, pb = wald(lg, [n for n in lg.params.index if n.startswith("C(klass")])
cls_res[assign] = means
print(f"unknown errors {assign}: rates Enter {means['Enter']:.2f}%, Exit {means['Exit']:.2f}%, Internal {means['Internal']:.2f}% "
f"| (a) LMM chi2({df}) = {chi2:.1f}, p = {p:.2g} | (b) logit chi2({dfb}) = {chi2b:.1f}, p = {pb:.2g}")
print("expected: Enter 3.4 +/- 0.34, Exit 2.3 +/- 0.23, Internal 1.4 +/- 0.14")
# Enter against the other two classes together, which the data do tell apart
t = trials.assign(enter=np.where(trials.klass.eq("Enter"), "Enter", "other"))
pr = t.groupby(["subject", "enter"]).error.mean().mul(100).rename("rate").reset_index()
m = smf.mixedlm("rate ~ C(enter)", pr, groups=pr["subject"]).fit()
chi2, df, p = wald(m, ["C(enter)[T.other]"])
print(f"Enter {pr[pr.enter=='Enter'].rate.mean():.2f}% against Internal and Exit together {pr[pr.enter=='other'].rate.mean():.2f}%: "
f"chi2(1) = {chi2:.1f}, p = {p:.2g}; participants with a higher Enter rate: "
f"{int((pr.pivot(index='subject', columns='enter', values='rate').Enter > pr.pivot(index='subject', columns='enter', values='rate').other).sum())} of {pr.subject.nunique()}")
print("errors per class, all participants: Enter", int(trials[(trials.error == 1) & trials.klass.eq("Enter")].shape[0]),
"; Internal or Exit", len(unk), "; trials per class:", trials.klass.value_counts().to_dict())
errors of unknown class (Internal or Exit): 44 of 225
unknown errors Internal: rates Enter 4.79%, Exit 0.00%, Internal 1.20% | (a) LMM chi2(2) = 63.8, p = 1.4e-14 | (b) logit chi2(2) = 73.0, p = 1.4e-16
unknown errors Exit: rates Enter 4.79%, Exit 1.18%, Internal 0.00% | (a) LMM chi2(2) = 63.6, p = 1.6e-14 | (b) logit chi2(2) = 74.2, p = 7.9e-17
expected: Enter 3.4 +/- 0.34, Exit 2.3 +/- 0.23, Internal 1.4 +/- 0.14
Enter 4.79% against Internal and Exit together 0.60%: chi2(1) = 33.4, p = 7.5e-09; participants with a higher Enter rate: 24 of 24
errors per class, all participants: Enter 181 ; Internal or Exit 44 ; trials per class: {'Enter': 3801, 'Exit': 3737, 'Internal': 3714}
Result: partially reproduced, resting on one error per block and the class of an error taken from the trial before its block. Error rates differ by class (LMM chi2(2) = 63.6 to 63.8, p < 1e-13, under both assignments of the 44 errors of unknown class), and the Enter rate is higher than the other two classes together in all 24 participants. The rates do not match: Enter is 4.79% (expected 3.4 +/- 0.34), and Exit and Internal together are 0.60%. Given all 44 errors of unknown class, Exit would reach 1.18% or Internal 1.20%, so neither can reach its expected value (2.3% and 1.4%) whichever class the unknown errors belong to, and whether Exit is above Internal cannot be told. Blocks after an Exit trial are 79% of blocks against the 48% the paper's rates imply (data check above), so Internal and Exit errors may be undercounted or merged into blocks.
supp-errors-stable-over-blocks¶
Error rates did not change over the blocks of the main experiment.
Per participant and run, errors over the run's 96 trials after the first; run as a linear term, random intercept and slope by participant (reading a); (b) trial-wise logistic model with participant fixed effects. Judged on (a) at p < 0.05.
rr = trials.groupby(["subject", "run"]).error.mean().mul(100).rename("rate").reset_index()
m = smf.mixedlm("rate ~ run", rr, groups=rr["subject"], re_formula="~run").fit()
chi2, df, p = wald(m, ["run"])
print(f"(a) LMM: slope {m.params['run']:.3f} percentage points per run, chi2({df}) = {chi2:.2f}, p = {p:.3f}; converged {m.converged}")
lg = smf.logit("error ~ run + C(subject)", trials).fit(disp=0)
chi2b, dfb, pb = wald(lg, ["run"])
print(f"(b) logit: log-odds slope {lg.params['run']:.3f} per run, chi2({dfb}) = {chi2b:.2f}, p = {pb:.3f}")
print("mean error rate per run (%):", rr.groupby("run").rate.mean().round(2).to_dict())
print("paper: chi2(1) = 0.08, p = .77")
(a) LMM: slope -0.110 percentage points per run, chi2(1) = 1.44, p = 0.231; converged True
(b) logit: log-odds slope -0.060 per run, chi2(1) = 1.53, p = 0.216
mean error rate per run (%): {1: 1.94, 2: 2.35, 3: 2.27, 4: 2.01, 5: 1.51}
paper: chi2(1) = 0.08, p = .77
Result: reproduced, resting on one error per block. The run effect on per-run error rates is not significant (slope -0.11 percentage points per run, chi2(1) = 1.44, p = 0.23; logistic variant p = 0.22). Errors are counted one per block.
supp-errors-not-by-age-category¶
Error rates were not affected by the age or the category judged, nor by their interaction.
The current category of every error trial follows from the trial before its block; its current age
is known only for reconstructed Enter errors (span reading, state_span). The full model of the
criterion (category, age, interaction; random intercept) therefore runs on Enter trials whose age is
known (labelled Enter trials and the reconstructed Enter errors), which leaves out the other errors;
the category alone is also tested on all trials and errors.
ta = trials[trials.klass.eq("Enter") & trials.cur_age.notna()]
pa = ta.groupby(["subject", "cur_cat", "cur_age"]).error.mean().mul(100).rename("rate").reset_index()
m = smf.mixedlm("rate ~ C(cur_cat) * C(cur_age)", pa, groups=pa["subject"]).fit()
for term in ["C(cur_cat)[T.H]", "C(cur_age)[T.y]", "C(cur_cat)[T.H]:C(cur_age)[T.y]"]:
chi2, df, p = wald(m, [term])
print(f"(a) LMM, Enter trials with known age ({int(ta.error.sum())} errors): {term}: coef {m.params[term]:.3f}, chi2 = {chi2:.2f}, p = {p:.3f}")
lg = smf.logit("error ~ C(cur_cat) * C(cur_age) + C(subject)", ta).fit(disp=0)
for term in ["C(cur_cat)[T.H]", "C(cur_age)[T.y]", "C(cur_cat)[T.H]:C(cur_age)[T.y]"]:
chi2, df, p = wald(lg, [term])
print(f"(b) logit: {term}: chi2 = {chi2:.2f}, p = {p:.3f}")
tc = trials[trials.cur_cat.notna()]
pc = tc.groupby(["subject", "cur_cat"]).error.mean().mul(100).rename("rate").reset_index()
m = smf.mixedlm("rate ~ C(cur_cat)", pc, groups=pc["subject"]).fit()
chi2, df, p = wald(m, ["C(cur_cat)[T.H]"])
print(f"category alone, all {int(tc.error.sum())} errors: faces {pc[pc.cur_cat=='F'].rate.mean():.2f}%, houses {pc[pc.cur_cat=='H'].rate.mean():.2f}%, "
f"chi2 = {chi2:.2f}, p = {p:.3f}")
print("cell rates (errors with known age), mean over participants (%):")
print(pa.groupby(["cur_cat", "cur_age"]).rate.mean().round(2).to_string())
(a) LMM, Enter trials with known age (98 errors): C(cur_cat)[T.H]: coef -0.588, chi2 = 0.67, p = 0.414 (a) LMM, Enter trials with known age (98 errors): C(cur_age)[T.y]: coef -0.352, chi2 = 0.24, p = 0.625 (a) LMM, Enter trials with known age (98 errors): C(cur_cat)[T.H]:C(cur_age)[T.y]: coef -0.069, chi2 = 0.00, p = 0.946
(b) logit: C(cur_cat)[T.H]: chi2 = 0.60, p = 0.437 (b) logit: C(cur_age)[T.y]: chi2 = 0.26, p = 0.610 (b) logit: C(cur_cat)[T.H]:C(cur_age)[T.y]: chi2 = 0.02, p = 0.879 category alone, all 224 errors: faces 1.82%, houses 2.22%, chi2 = 1.45, p = 0.229 cell rates (errors with known age), mean over participants (%):
cur_cat cur_age
F o 3.16
y 2.81
H o 2.58
y 2.15
Result: reproduced, resting on the span reading of error states. On Enter trials whose age is known, category (p = 0.41), age (p = 0.63) and their interaction (p = 0.95) do not affect error rates, and the logistic variant agrees (p = 0.44, 0.61, 0.88). The category alone, over all trials, has p = 0.23. The age test covers only the 98 reconstructed Enter errors.
supp-rt-stable-over-blocks, supp-rt-faces-faster, supp-rt-by-trial-class¶
Reaction times are the durations of the RT event files, given to 0.1 s, for the labelled (correct) trials after the first. Error trials have no RT of their own in the files (an error block's duration is its span), so the reading "all responded trials" cannot be run; only correct trials are used. RTs are medians per participant and factor cell (Methods, Behavioral Analyses).
lab["rt_ms"] = lab.rt * 1000
rb = lab.groupby(["subject", "run"]).rt_ms.median().rename("rt").reset_index()
m = smf.mixedlm("rt ~ run", rb, groups=rb["subject"], re_formula="~run").fit()
chi2, df, p = wald(m, ["run"])
print(f"LMM: slope {m.params['run']:.2f} ms per run, chi2({df}) = {chi2:.2f}, p = {p:.3f}; converged {m.converged}")
print("mean of participant medians per run (ms):", rb.groupby("run").rt.mean().round(1).to_dict())
print("paper: chi2(1) = 1.7, p = .18")
LMM: slope -8.55 ms per run, chi2(1) = 1.74, p = 0.187; converged True
mean of participant medians per run (ms): {1: 991.7, 2: 987.5, 3: 997.9, 4: 995.8, 5: 941.3}
paper: chi2(1) = 1.7, p = .18
Result: reproduced. RT medians do not change over runs (slope -8.6 ms per run, chi2(1) = 1.74, p = 0.19; the paper gives chi2(1) = 1.7, p = .18).
rc = lab.groupby(["subject", "cur_cat", "cur_age"]).rt_ms.median().rename("rt").reset_index()
m = smf.mixedlm("rt ~ C(cur_cat) * C(cur_age)", rc, groups=rc["subject"]).fit()
for term in ["C(cur_cat)[T.H]", "C(cur_age)[T.y]", "C(cur_cat)[T.H]:C(cur_age)[T.y]"]:
chi2, df, p = wald(m, [term])
print(f"LMM {term}: coef {m.params[term]:.1f} ms, chi2(1) = {chi2:.1f}, p = {p:.3g}")
# variant: random slopes of category by participant, which the criterion's random-intercept model leaves out
m2 = smf.mixedlm("rt ~ C(cur_cat) * C(cur_age)", rc, groups=rc["subject"], re_formula="~C(cur_cat)").fit()
chi2, df, p = wald(m2, ["C(cur_cat)[T.H]"])
print(f"variant with by-participant category slopes: category coef {m2.params['C(cur_cat)[T.H]']:.1f} ms, chi2(1) = {chi2:.1f}, p = {p:.3g}; converged {m2.converged}")
cat_med = lab.groupby(["subject", "cur_cat"]).rt_ms.median().unstack()
print("mean over participants of category medians (ms): faces %.1f, houses %.1f" % tuple(cat_med[["F", "H"]].mean()))
print("median of pooled trials (ms): faces %.1f, houses %.1f" % (lab[lab.cur_cat == "F"].rt_ms.median(), lab[lab.cur_cat == "H"].rt_ms.median()))
print("participants faster for faces:", int((cat_med.F < cat_med.H).sum()), "of", len(cat_med),
"; expected faces 944 +/- 94.4, houses 1028 +/- 102.8")
LMM C(cur_cat)[T.H]: coef 87.5 ms, chi2(1) = 58.7, p = 1.8e-14 LMM C(cur_age)[T.y]: coef 4.2 ms, chi2(1) = 0.1, p = 0.715 LMM C(cur_cat)[T.H]:C(cur_age)[T.y]: coef -8.3 ms, chi2(1) = 0.3, p = 0.606
variant with by-participant category slopes: category coef 87.5 ms, chi2(1) = 54.1, p = 1.91e-13; converged True mean over participants of category medians (ms): faces 945.8, houses 1012.5 median of pooled trials (ms): faces 900.0, houses 1000.0 participants faster for faces: 15 of 24 ; expected faces 944 +/- 94.4, houses 1028 +/- 102.8
Result: reproduced. Faces are faster than houses (mean of participant medians 945.8 against 1012.5 ms, both within tolerance; LMM category chi2(1) = 58.7, p = 1.8e-14), with no effect of age (p = 0.72) and no interaction (p = 0.61). With by-participant category slopes added (a variant; the criterion's model has a random intercept only), category chi2(1) = 54.1, p = 1.9e-13. Pooled trial medians are 900 and 1000 ms. RTs are given to 0.1 s, so many participant medians tie: 15 of 24 participants have a lower face median.
rk = lab.groupby(["subject", "klass"]).rt_ms.median().rename("rt").reset_index()
m = smf.mixedlm("rt ~ C(klass, Treatment('Internal'))", rk, groups=rk["subject"]).fit()
chi2, df, p = wald(m, [n for n in m.params.index if n.startswith("C(klass")])
print(f"LMM class effect chi2({df}) = {chi2:.1f}, p = {p:.3g}")
print("mean of participant medians (ms):", rk.groupby("klass").rt.mean().round(1).to_dict())
print("participants with Internal below both Enter and Exit:",
int(((rk.pivot(index="subject", columns="klass", values="rt").Internal < rk.pivot(index="subject", columns="klass", values="rt").Enter) &
(rk.pivot(index="subject", columns="klass", values="rt").Internal < rk.pivot(index="subject", columns="klass", values="rt").Exit)).sum()), "of", rk.subject.nunique())
LMM class effect chi2(2) = 116.6, p = 4.81e-26
mean of participant medians (ms): {'Enter': 1020.8, 'Exit': 1022.9, 'Internal': 883.3}
participants with Internal below both Enter and Exit: 21 of 24
Result: reproduced. Trial class affects RT (chi2(2) = 116.6, p = 4.8e-26); Internal trials are the fastest (883 ms against 1021 for Enter and 1023 for Exit), in 21 of 24 participants.
Decoding: one pass over the participants¶
For each participant (ofc4.analyse_subject): a GLM per run with the data set's 16 state columns,
its 'first' and 'err' columns and the ten nuisance time courses; the beta maps z-scored within run
under both readings of O2 ('voxel', 'map'); a linear SVM (cost 1) trained leave-one-run-out on the
16 states and tested on the held-out run's maps, with 10 permuted-label classifiers per fold; the
same for the 8-state space (pairs of state columns summed, which equals modelling the two states as
one event type) and the two-back space (regressors rebuilt from the events), and the 16-state space
with its regressors rebuilt from the events (a check of the rebuilt regressors); and the 16-state
classifier tested on single volumes of the held-out run. Each participant has its own seed (1000 +
participant number). No beta map or volume is smoothed: the files hold no voxel coordinates (O1).
t = time.time()
R = {s: ofc4.analyse_subject(s, seed=1000 + s) for s in SUBJ}
print(f"decoding done in {time.time() - t:.0f} s")
PRED = pd.concat([r["pred"] for r in R.values()], ignore_index=True)
PERM = pd.concat([r["perm"] for r in R.values()], ignore_index=True)
TRIAL = pd.concat([r["trial"] for r in R.values()], ignore_index=True)
PRED["c"] = PRED.true == PRED.pred
PERM["c"] = PERM.true == PERM.pred
print("test maps per participant, space and reading:", PRED.groupby(["subject", "z", "space"]).size().describe().round(1).to_dict())
miss = {(s, r, k): v for s, res in R.items() for r, d in res["missing_two"].items() for k, v in d.items() if v}
print("two-back states absent from a run (left out of training and test):", miss)
print("rebuilt 16-state regressors against the data set's columns: r from",
round(min(min(v) for res in R.values() for v in res["rebuilt_r"].values()), 3), "to",
round(max(max(v) for res in R.values() for v in res["rebuilt_r"].values()), 3))
decoding done in 196 s
test maps per participant, space and reading: {'count': 336.0, 'mean': 73.7, 'std': 14.2, 'min': 32.0, '25%': 80.0, '50%': 80.0, '75%': 80.0, 'max': 80.0}
two-back states absent from a run (left out of training and test): {}
rebuilt 16-state regressors against the data set's columns: r from 0.983 to 1.0
# Timing check: the 16-state classifier on single volumes of labelled trials, from the volume at the
# stimulus onset (lag 0) to 5 TRs after it. With the event times on the data's clock, state
# information should rise a few seconds after onset (hemodynamic delay).
# Internal trials repeat the state of the trial before, so the volume at onset still carries the
# state of the previous trial; Enter and Exit trials, whose state differs from the previous trial's,
# give the cleaner contrast.
lagt = TRIAL[TRIAL.vol.str.startswith("lag") & TRIAL.cond.isin(ofc4.STATES)].copy()
lagt["c"] = lagt.true == lagt.pred
lagt["switch"] = lagt.cond.map(ofc4.trial_class) != "Internal"
for subset, sel in (("all labelled trials", lagt), ("Enter and Exit trials", lagt[lagt.switch])):
lag_acc = sel.groupby(["z", "vol", "subject"]).c.mean().unstack("vol")
for z in ("voxel", "map"):
m = lag_acc.loc[z]
print(f"{subset}, z by {z}: mean 16-way accuracy (%) by lag:", (100 * m.mean()).round(2).to_dict())
print(" t against 6.25% by lag:", {k: round(float(stats.ttest_1samp(m[k], 1 / 16).statistic), 2) for k in m})
lag_acc = lagt[lagt.switch].groupby(["z", "vol", "subject"]).c.mean().unstack("vol")
fig, ax = plt.subplots(figsize=(5, 3))
for z in ("voxel", "map"):
m = lag_acc.loc[z]
ax.errorbar(range(6), 100 * m.mean(), 100 * m.sem(), label=f"z by {z}", capsize=3)
ax.axhline(100 / 16, ls="--", c="k")
ax.set_xlabel("volume after onset (TRs)"); ax.set_ylabel("16-way accuracy (%), Enter and Exit trials"); ax.legend()
print("plotted means range from", round(100 * lag_acc.mean().min(), 2), "to", round(100 * lag_acc.mean().max(), 2), "%")
plt.show()
va = TRIAL[TRIAL.vol.isin(["after", "nearest"]) & TRIAL.cond.isin(ofc4.STATES)].copy()
va["c"] = va.true == va.pred
va["klass"] = va.cond.map(ofc4.trial_class)
print("16-way accuracy (%) on single volumes 2 TRs after onset, labelled trials by class (mean over participants):")
print((100 * va.groupby(["z", "vol", "klass", "subject"]).c.mean().groupby(["z", "vol", "klass"]).mean()).round(2).unstack("klass").to_string())
SHIFT = pd.concat([r["shift"] for r in R.values()], ignore_index=True)
sh = SHIFT.groupby(["subject", "shift"]).acc.mean().unstack()
print("16-way accuracy (%, z by voxel) with the state regressors shifted by k TRs against the time series:")
print(" mean:", (100 * sh.mean()).round(2).to_dict())
print(" t against 6.25%:", {k: round(float(stats.ttest_1samp(sh[k], 1 / 16).statistic), 2) for k in sh})
all labelled trials, z by voxel: mean 16-way accuracy (%) by lag: {'lag0': 8.49, 'lag1': 9.17, 'lag2': 9.09, 'lag3': 8.74, 'lag4': 8.3, 'lag5': 8.2}
t against 6.25% by lag: {'lag0': 7.03, 'lag1': 9.58, 'lag2': 8.26, 'lag3': 8.94, 'lag4': 6.48, 'lag5': 7.5}
all labelled trials, z by map: mean 16-way accuracy (%) by lag: {'lag0': 6.08, 'lag1': 6.73, 'lag2': 7.04, 'lag3': 7.02, 'lag4': 6.25, 'lag5': 6.1}
t against 6.25% by lag: {'lag0': -0.66, 'lag1': 1.42, 'lag2': 3.26, 'lag3': 2.34, 'lag4': -0.01, 'lag5': -0.55}
Enter and Exit trials, z by voxel: mean 16-way accuracy (%) by lag: {'lag0': 5.68, 'lag1': 6.53, 'lag2': 6.46, 'lag3': 6.66, 'lag4': 5.54, 'lag5': 5.78}
t against 6.25% by lag: {'lag0': -1.33, 'lag1': 0.61, 'lag2': 0.42, 'lag3': 0.94, 'lag4': -1.64, 'lag5': -1.18}
Enter and Exit trials, z by map: mean 16-way accuracy (%) by lag: {'lag0': 5.34, 'lag1': 5.81, 'lag2': 6.54, 'lag3': 6.39, 'lag4': 5.69, 'lag5': 5.65}
t against 6.25% by lag: {'lag0': -2.53, 'lag1': -1.08, 'lag2': 0.82, 'lag3': 0.35, 'lag4': -1.67, 'lag5': -1.85}
plotted means range from 5.51 to 6.52 %
16-way accuracy (%) on single volumes 2 TRs after onset, labelled trials by class (mean over participants):
klass Enter Exit Internal
z vol
map after 8.30 4.85 8.04
nearest 7.87 4.80 8.38
voxel after 0.90 11.79 14.30
nearest 0.84 12.53 14.28
16-way accuracy (%, z by voxel) with the state regressors shifted by k TRs against the time series:
mean: {-6: 6.34, -4: 6.45, -2: 7.32, -1: 7.88, 0: 8.93, 1: 8.42, 2: 9.18, 4: 7.08, 6: 7.79}
t against 6.25%: {-6: 0.24, -4: 0.37, -2: 1.37, -1: 2.49, 0: 3.18, 1: 3.61, 2: 3.97, 4: 1.76, 6: 2.18}
# What single volumes of Enter trials are decoded as. An Enter trial follows an Exit trial, and the
# two states differ, so a volume that still carries the trial before is decoded as the state of the
# Exit trial before. Per participant: share of predictions equal to the state of the trial before,
# against the share of all the participant's predictions (labelled trials, same z-scoring and volume)
# that name that state, the rate expected without lag. Error trials of the span reading as well.
tq = TRIAL[TRIAL.vol.isin(["after", "nearest"] + [f"lag{k}" for k in range(6)])].copy()
prev = tq[["z", "vol", "subject", "run", "pos", "cond"]].assign(pos=tq.pos + 1).rename(columns={"cond": "prev_cond"})
tq = tq.merge(prev, on=["z", "vol", "subject", "run", "pos"], how="left")
tq["klass"] = tq.cond.map(lambda x: ofc4.trial_class(x) if x in ofc4.STATES else ("error" if x == "err" else None))
tq = tq[tq.prev_cond.isin(ofc4.STATES)]
base = tq[tq.cond.isin(ofc4.STATES)].groupby(["z", "vol", "subject"]).pred.value_counts(normalize=True).rename("base").reset_index()
tq = tq.merge(base.rename(columns={"pred": "prev_cond"}), on=["z", "vol", "subject", "prev_cond"], how="left").fillna({"base": 0})
tq["as_prev"] = (tq.pred == tq.prev_cond).astype(float)
tq["as_own"] = (tq.pred == tq.true_span).astype(float)
sel = tq[tq.klass.isin(["Enter", "Exit"]) | (tq.klass.eq("error") & tq.true_span.notna())]
lagtab = sel.groupby(["z", "vol", "klass", "subject"])[["as_own", "as_prev", "base"]].mean().groupby(["z", "vol", "klass"])
out = (100 * lagtab.mean()).round(2)
out["t_prev_vs_base"] = sel.groupby(["z", "vol", "klass", "subject"])[["as_prev", "base"]].mean().assign(d=lambda d: d.as_prev - d.base).groupby(["z", "vol", "klass"]).d.apply(lambda d: stats.ttest_1samp(d, 0).statistic).round(2)
out["n_subj"] = lagtab.size()
print("share (%) of single volumes decoded as the trial's own state, as the state of the trial before, and that state's base rate:")
print(out.to_string())
pc = TRIAL[TRIAL.vol.eq("after") & TRIAL.cond.isin(ofc4.STATES)].assign(pred_class=lambda d: d.pred.map(ofc4.trial_class))
print("share (%) of predictions naming a state of each class, labelled trials, volume 'after' (8 Enter, 4 Exit and 4 Internal states: 50, 25 and 25% without bias):")
print((100 * pc.groupby(["z", "subject"]).pred_class.value_counts(normalize=True).unstack().groupby("z").mean()).round(2).to_string())
share (%) of single volumes decoded as the trial's own state, as the state of the trial before, and that state's base rate:
as_own as_prev base t_prev_vs_base n_subj
z vol klass
map after Enter 8.37 4.10 4.15 -0.16 24
Exit 4.82 7.42 6.95 1.02 24
error 7.06 7.35 3.93 0.70 21
lag0 Enter 7.02 5.22 4.16 3.00 24
Exit 3.80 7.83 6.93 1.63 24
error 3.57 4.74 3.86 0.36 21
lag1 Enter 7.82 4.29 4.02 0.89 24
Exit 3.94 7.69 7.05 1.23 24
error 9.53 5.32 3.69 0.89 21
lag2 Enter 8.37 4.10 4.15 -0.16 24
Exit 4.82 7.42 6.95 1.02 24
error 7.06 7.35 3.93 0.70 21
lag3 Enter 8.20 3.52 4.08 -2.42 24
Exit 4.72 6.53 7.03 -1.24 24
error 6.16 1.19 3.70 -2.87 21
lag4 Enter 7.82 3.86 4.00 -0.70 24
Exit 3.63 6.79 7.03 -0.72 24
error 2.42 3.06 3.71 -0.44 21
lag5 Enter 7.39 3.75 4.07 -1.13 24
Exit 4.02 5.88 7.01 -4.08 24
error 5.63 2.98 3.88 -0.37 21
nearest Enter 7.91 4.26 4.22 0.11 24
Exit 4.78 7.40 6.97 0.71 24
error 8.93 7.31 3.95 0.71 21
voxel after Enter 0.92 10.78 10.90 -0.20 24
Exit 11.72 5.74 6.38 -1.60 24
error 1.59 17.38 11.93 1.01 21
lag0 Enter 0.72 12.63 10.98 3.12 24
Exit 10.58 6.62 6.32 0.93 24
error 0.00 12.19 12.09 0.02 21
lag1 Enter 0.83 11.73 10.82 1.68 24
Exit 11.91 6.13 6.43 -0.88 24
error 1.59 18.21 12.13 1.19 21
lag2 Enter 0.92 10.78 10.90 -0.20 24
Exit 11.72 5.74 6.38 -1.60 24
error 1.59 17.38 11.93 1.01 21
lag3 Enter 0.78 10.22 10.82 -1.51 24
Exit 12.09 5.83 6.45 -1.72 24
error 0.00 11.91 11.98 -0.02 21
lag4 Enter 0.65 11.01 10.80 0.33 24
Exit 10.28 5.94 6.47 -1.58 24
error 1.59 19.85 11.89 1.13 21
lag5 Enter 0.56 10.68 10.84 -0.28 24
Exit 10.73 5.49 6.48 -4.01 24
error 0.00 15.04 12.08 0.52 21
nearest Enter 0.82 11.40 11.05 0.83 24
Exit 12.41 5.74 6.34 -1.86 24
error 1.59 22.62 12.03 1.75 21
share (%) of predictions naming a state of each class, labelled trials, volume 'after' (8 Enter, 4 Exit and 4 Internal states: 50, 25 and 25% without bias):
pred_class Enter Exit Internal
z
map 56.74 16.78 26.48
voxel 5.60 43.69 50.71
Timing of the time series against the events. With the state regressors shifted against the time series, 16-way accuracy is at chance for shifts of -4 and -6 TRs (6.45% and 6.34%), rises to 8.93% at no shift and stays above chance for shifts of +1 and +2 TRs (8.42% and 9.18%). On single volumes, with maps z-scored by map, accuracy on labelled trials is at chance at onset (6.08%) and highest 2 to 3 TRs after it (7.04% and 7.02%). Enter and Exit trials alone, whose state differs from the trial before, decode near chance at every lag (5.34% to 6.66%), so they cannot show the delay. This is consistent with the events and the time series sharing a clock to within about two TRs; it does not show it more precisely. 16-way accuracy also stays above chance with the regressors shifted by +6 TRs (7.79%, t = 2.18) and -1 TR (7.88%, t = 2.49), so part of the decoded signal changes slowly, over several trials, as runs of Internal trials that repeat a state would make it.
On the single volume used for the trial-wise claims (2 TRs after onset), accuracy depends strongly on the trial's class: with maps z-scored by voxel, Enter trials decode at 0.90%, Exit at 11.79% and Internal at 14.30% (volume 'after'); with maps z-scored by map, 8.30%, 4.85% and 8.04%. Any comparison of trials of one class with trials of mixed class is therefore confounded by class, which bears on the error-trial claims (fig4b), where every error trial that can be given a state is an Enter trial.
Why correct Enter trials decode at 0.90% with maps z-scored by voxel, far below chance. If the volume 2 TRs after onset still carried the Exit trial before, Enter trials would be decoded as that Exit state more often than the classifier names it otherwise. They are not: at the volume 'after' the share decoded as the previous trial's state is 10.78%, the same as the base rate at which that state is predicted on all labelled trials (10.90%, t = -0.20); an excess appears only at onset (lag 0, 12.63% against 10.98%, t = 3.12). The low accuracy comes from a bias of the classifier: with maps z-scored by voxel it names an Enter state on 5.60% of labelled volumes, where the 8 Enter states of 16 would take 50% without bias (z by map, 56.74%). On single volumes of this classifier, Enter trials are therefore near 0% whatever the signal, and comparisons of Enter error trials with trials of other classes measure this bias. Error trials of the span reading are decoded as the previous state on 17.38% of volumes against a base rate of 11.93% (t = 1.01, 21 participants), not a significant excess.
def comp_acc(df, comps=ofc4.COMPONENTS, index_map=None):
"""Per participant accuracy on each component; index_map gives the label position of each component."""
index_map = index_map or {"cur_cat": 0, "cur_age": 1, "prev_cat": 2, "prev_age": 3}
out = {}
for comp in comps:
i = index_map[comp]
out[comp] = (df.true.str[i] == df.pred.str[i]).groupby(df.subject).mean()
return pd.DataFrame(out)
def ttest(x, mu, alternative="greater"):
r = stats.ttest_1samp(x, mu, alternative=alternative)
return float(r.statistic), float(r.pvalue)
def boot_p(x, mu, n=10000, seed=0):
"""Bootstrap over participants (Supplemental Results 1.2): share of resampled means at or below mu."""
rng = np.random.default_rng(seed)
x = np.asarray(x)
means = x[rng.integers(0, len(x), (n, len(x)))].mean(1)
return float((means <= mu).mean())
ACC = PRED.groupby(["z", "space", "subject"]).c.mean()
fig2a-16way-decoding¶
The 16 task states can be decoded from multivoxel patterns in anatomically defined OFC above chance (mean 12.2% against 6.25%).
Readings: z-scoring O2 ('voxel', 'map'); the t-test side O4 (one- and two-sided). Smoothing (O1)
can only be run unsmoothed. Post-error trials (O3): the data set's 'err' column follows a short event
at each block onset (data check), so with the data set's columns ('16') the trials inside a block after
the error trial are in no regressor, which is neither reading of O3. Both readings are rebuilt from
the 'err' events (ofc4.subject_designs): (a) '16o3a', the 'err' column with a separate nuisance
boxcar from onset + 4.5 s (the mean trial duration, Methods, Task) to the end of the block; (b)
'16o3b', one 'err' boxcar over the whole block.
for z in ("voxel", "map"):
for space in ("16", "16rebuilt", "16o3a", "16o3b"):
a = ACC.loc[(z, space)]
t1, p1 = ttest(a, 1 / 16)
t2, p2 = ttest(a, 1 / 16, "two-sided")
print(f"z by {z}, {space}: mean {100 * a.mean():.2f}% (SEM {100 * a.sem():.2f}), above chance in {int((a > 1 / 16).sum())} of {len(a)}; "
f"t{len(a) - 1} = {t1:.2f}, one-sided p = {p1:.2g}, two-sided p = {p2:.2g}")
print("spaces: 16 = the data set's columns; 16rebuilt = state regressors rebuilt from the events; 16o3a and 16o3b = the data set's\n"
"state columns with the two readings of O3 for the error blocks (a: separate post-error regressor, b: 'err' over the whole block)")
print("expected 12.2 +/- 1.22; per participant (%):")
print((100 * ACC.loc[(slice(None), "16")].unstack("z")).round(2).T.to_string())
z by voxel, 16: mean 8.93% (SEM 0.84), above chance in 15 of 24; t23 = 3.18, one-sided p = 0.0021, two-sided p = 0.0041 z by voxel, 16rebuilt: mean 9.57% (SEM 0.80), above chance in 18 of 24; t23 = 4.13, one-sided p = 0.00021, two-sided p = 0.00041 z by voxel, 16o3a: mean 8.97% (SEM 0.84), above chance in 17 of 24; t23 = 3.23, one-sided p = 0.0019, two-sided p = 0.0037 z by voxel, 16o3b: mean 9.38% (SEM 0.92), above chance in 14 of 24; t23 = 3.40, one-sided p = 0.0012, two-sided p = 0.0024 z by map, 16: mean 8.92% (SEM 0.69), above chance in 16 of 24; t23 = 3.89, one-sided p = 0.00037, two-sided p = 0.00073 z by map, 16rebuilt: mean 9.19% (SEM 0.73), above chance in 16 of 24; t23 = 4.04, one-sided p = 0.00025, two-sided p = 0.00051 z by map, 16o3a: mean 8.37% (SEM 0.66), above chance in 15 of 24; t23 = 3.24, one-sided p = 0.0018, two-sided p = 0.0036 z by map, 16o3b: mean 8.10% (SEM 0.65), above chance in 15 of 24; t23 = 2.85, one-sided p = 0.0045, two-sided p = 0.0091 spaces: 16 = the data set's columns; 16rebuilt = state regressors rebuilt from the events; 16o3a and 16o3b = the data set's state columns with the two readings of O3 for the error blocks (a: separate post-error regressor, b: 'err' over the whole block) expected 12.2 +/- 1.22; per participant (%): subject 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 z map 7.81 8.75 6.25 12.5 11.25 5.0 12.5 7.5 15.00 6.25 1.25 10.0 3.75 11.25 12.50 12.50 6.25 11.25 10.00 11.25 6.25 11.25 7.5 6.25 voxel 9.38 6.25 3.75 12.5 10.00 5.0 10.0 10.0 8.75 8.75 2.50 15.0 5.00 11.25 11.25 18.75 3.75 12.50 3.75 6.25 11.25 7.50 15.0 6.25
Result: partially reproduced. 16-way decoding is above chance under both z-scorings (8.93% and 8.92%, one-sided p = 0.0021 and 0.00037, two-sided p = 0.0041 and 0.00073), but below the 10.98 to 13.42% of the tolerance. With regressors rebuilt from the events it is 9.57% and 9.19%. Under the two readings of O3 it is 8.97% and 9.38% (z by voxel) and 8.37% and 8.10% (z by map), all above chance (one-sided p <= 0.0045) and all below the tolerance, so the outcome is the same under every reading. The maps are unsmoothed, as the files hold no coordinates; the paper smoothed them.
fig2-two-back-state-space¶
A classifier for an alternative 16-state space built from the task-irrelevant category and age of two trials back decodes worse than the hypothesized state space (8.1% against 12.2%).
Each labelled trial is relabelled by its current category and age and the category and age of the trial two back (the previous components of the trial before it). Readings for a trial whose trial before is the first trial or an error block: (a) left out, as one nuisance regressor; (b) after a block with one reconstructed error, which is an Enter error repeated as the same trial (Methods, Task), the trial two back as presented is the error trial, whose state is that of its repetition, the trial before; others left out. The hypothesized space is compared both as the data set's columns ('16') and as regressors rebuilt from the events the same way as the two-back ones ('16rebuilt').
for z in ("voxel", "map"):
for hyp in ("16", "16rebuilt"):
h = ACC.loc[(z, hyp)]
for reading in ("a", "b"):
a = ACC.loc[(z, "two" + reading)]
t, p = stats.ttest_rel(a, h)
print(f"z by {z}, against {hyp}, reading {reading}: two-back {100 * a.mean():.2f}% against hypothesized {100 * h.mean():.2f}%; "
f"paired t{len(a) - 1} = {t:.2f}, two-sided p = {p:.3g}; two-back lower in {int((a < h).sum())} of {len(a)}")
print("expected two-back 8.1 +/- 0.81")
print("labelled trials left out of the two-back states (nuisance regressor):",
{k: sum(res["two_leftout"][r][k] for res in R.values() for r in res["two_leftout"]) for k in "ab"},
"of", int(EV.is_state.sum()))
z by voxel, against 16, reading a: two-back 5.96% against hypothesized 8.93%; paired t23 = -3.75, two-sided p = 0.00105; two-back lower in 15 of 24
z by voxel, against 16, reading b: two-back 6.17% against hypothesized 8.93%; paired t23 = -3.29, two-sided p = 0.00317; two-back lower in 14 of 24
z by voxel, against 16rebuilt, reading a: two-back 5.96% against hypothesized 9.57%; paired t23 = -4.90, two-sided p = 6.04e-05; two-back lower in 19 of 24
z by voxel, against 16rebuilt, reading b: two-back 6.17% against hypothesized 9.57%; paired t23 = -4.19, two-sided p = 0.000349; two-back lower in 18 of 24
z by map, against 16, reading a: two-back 7.88% against hypothesized 8.92%; paired t23 = -1.58, two-sided p = 0.128; two-back lower in 13 of 24
z by map, against 16, reading b: two-back 7.88% against hypothesized 8.92%; paired t23 = -1.51, two-sided p = 0.144; two-back lower in 16 of 24
z by map, against 16rebuilt, reading a: two-back 7.88% against hypothesized 9.19%; paired t23 = -1.90, two-sided p = 0.0698; two-back lower in 15 of 24
z by map, against 16rebuilt, reading b: two-back 7.88% against hypothesized 9.19%; paired t23 = -1.88, two-sided p = 0.0729; two-back lower in 16 of 24
expected two-back 8.1 +/- 0.81
labelled trials left out of the two-back states (nuisance regressor): {'a': 334, 'b': 295} of 11071
Result: partially reproduced. The two-back space decodes worse than the hypothesized one, but no reading meets both parts of the criterion. With maps z-scored by voxel (the assumption), it is lower (5.96% against 8.93%, t23 = -3.75, p = 0.001; reading b 6.17%, p = 0.003; against the rebuilt hypothesized space 9.57%, p = 6e-5 and 0.0003), but its accuracy is at chance, below the 7.29 to 8.91% of the tolerance. With maps z-scored by map, its accuracy is within tolerance (7.88% in both readings) but not significantly lower (p = 0.13 and 0.14 against '16', 0.070 and 0.073 against '16rebuilt').
fig2b-hidden-components and fig2b-current-age-absent¶
Within the 16-way classification, the three hidden, task-relevant components (previous category, previous age, current category) are each classified above chance.
The observable component, current age, is not classified above chance within the 16-way classification.
Readings: one-sided t-test against 50%, and the bootstrap of participants (10000 resamples, Supplemental Results 1.2), whose p is the share of resampled means at or below 50%.
COMP16 = {}
for z, space in [(z, sp) for sp in ("16", "16o3a", "16o3b") for z in ("voxel", "map")]:
ca = comp_acc(PRED[(PRED.z == z) & (PRED.space == space)])
if space == "16":
COMP16[z] = ca
for comp in ["prev_cat", "prev_age", "cur_cat", "cur_age"]:
t, p = ttest(ca[comp], 0.5)
print(f"{space}, z by {z}, {comp}: mean {100 * ca[comp].mean():.2f}% (SEM {100 * ca[comp].sem():.2f}), t{len(ca) - 1} = {t:.2f}, "
f"one-sided p = {p:.3g}, bootstrap p = {boot_p(ca[comp], 0.5, seed=11):.4f}, above 50% in {int((ca[comp] > 0.5).sum())} of {len(ca)}")
print("paper: previous category p = 0.014, previous age p = 0.003, current category p = 0.007 (t26 > 3); current age p > 0.72")
16, z by voxel, prev_cat: mean 52.03% (SEM 1.17), t23 = 1.74, one-sided p = 0.0476, bootstrap p = 0.0458, above 50% in 15 of 24 16, z by voxel, prev_age: mean 49.08% (SEM 1.26), t23 = -0.74, one-sided p = 0.765, bootstrap p = 0.7682, above 50% in 11 of 24 16, z by voxel, cur_cat: mean 50.26% (SEM 1.19), t23 = 0.22, one-sided p = 0.415, bootstrap p = 0.4137, above 50% in 13 of 24 16, z by voxel, cur_age: mean 50.83% (SEM 1.33), t23 = 0.63, one-sided p = 0.268, bootstrap p = 0.2667, above 50% in 11 of 24 16, z by map, prev_cat: mean 51.78% (SEM 0.96), t23 = 1.86, one-sided p = 0.038, bootstrap p = 0.0301, above 50% in 13 of 24 16, z by map, prev_age: mean 51.16% (SEM 1.15), t23 = 1.01, one-sided p = 0.161, bootstrap p = 0.1557, above 50% in 13 of 24 16, z by map, cur_cat: mean 51.78% (SEM 1.35), t23 = 1.32, one-sided p = 0.0997, bootstrap p = 0.0901, above 50% in 15 of 24 16, z by map, cur_age: mean 47.92% (SEM 1.18), t23 = -1.76, one-sided p = 0.954, bootstrap p = 0.9692, above 50% in 8 of 24 16o3a, z by voxel, prev_cat: mean 51.63% (SEM 1.32), t23 = 1.23, one-sided p = 0.116, bootstrap p = 0.1081, above 50% in 13 of 24 16o3a, z by voxel, prev_age: mean 49.80% (SEM 1.11), t23 = -0.18, one-sided p = 0.569, bootstrap p = 0.5701, above 50% in 11 of 24 16o3a, z by voxel, cur_cat: mean 50.51% (SEM 0.94), t23 = 0.54, one-sided p = 0.296, bootstrap p = 0.2868, above 50% in 15 of 24 16o3a, z by voxel, cur_age: mean 51.72% (SEM 1.27), t23 = 1.35, one-sided p = 0.0945, bootstrap p = 0.0873, above 50% in 16 of 24 16o3a, z by map, prev_cat: mean 52.11% (SEM 0.89), t23 = 2.38, one-sided p = 0.0129, bootstrap p = 0.0093, above 50% in 16 of 24 16o3a, z by map, prev_age: mean 49.71% (SEM 1.02), t23 = -0.28, one-sided p = 0.609, bootstrap p = 0.6232, above 50% in 10 of 24 16o3a, z by map, cur_cat: mean 51.05% (SEM 1.04), t23 = 1.02, one-sided p = 0.16, bootstrap p = 0.1517, above 50% in 12 of 24 16o3a, z by map, cur_age: mean 48.50% (SEM 1.74), t23 = -0.86, one-sided p = 0.801, bootstrap p = 0.8101, above 50% in 10 of 24 16o3b, z by voxel, prev_cat: mean 52.60% (SEM 1.30), t23 = 2.00, one-sided p = 0.0288, bootstrap p = 0.0219, above 50% in 14 of 24 16o3b, z by voxel, prev_age: mean 49.19% (SEM 1.15), t23 = -0.70, one-sided p = 0.756, bootstrap p = 0.7637, above 50% in 12 of 24 16o3b, z by voxel, cur_cat: mean 50.99% (SEM 1.10), t23 = 0.90, one-sided p = 0.188, bootstrap p = 0.1819, above 50% in 15 of 24 16o3b, z by voxel, cur_age: mean 51.47% (SEM 1.31), t23 = 1.12, one-sided p = 0.137, bootstrap p = 0.1285, above 50% in 16 of 24 16o3b, z by map, prev_cat: mean 51.99% (SEM 1.04), t23 = 1.92, one-sided p = 0.0339, bootstrap p = 0.0246, above 50% in 14 of 24
16o3b, z by map, prev_age: mean 49.97% (SEM 1.18), t23 = -0.02, one-sided p = 0.509, bootstrap p = 0.5117, above 50% in 9 of 24 16o3b, z by map, cur_cat: mean 51.51% (SEM 1.16), t23 = 1.31, one-sided p = 0.102, bootstrap p = 0.0926, above 50% in 12 of 24 16o3b, z by map, cur_age: mean 48.15% (SEM 1.66), t23 = -1.12, one-sided p = 0.862, bootstrap p = 0.8741, above 50% in 11 of 24 paper: previous category p = 0.014, previous age p = 0.003, current category p = 0.007 (t26 > 3); current age p > 0.72
Result: partially reproduced. Of the three hidden components, only the previous category is above 50% (52.03%, one-sided p = 0.048; bootstrap p = 0.046; z by map 51.78%, p = 0.038). Previous age (49.08%, p = 0.77) and current category (50.26%, p = 0.42) are not; z-scored by map they are 51.16% (p = 0.16) and 51.78% (p = 0.10). Under the two readings of O3 the pattern is the same: previous category reaches p < 0.05 in three of four (p = 0.013 to 0.034; 0.12 for O3 a, z by voxel), more strongly than in the main reading for O3 (a) with maps z-scored by map (p = 0.013 against 0.038), previous age and current category in none (p >= 0.10).
Result: reproduced. Current age is not above chance (50.83%, one-sided p = 0.27; z by map 47.92%, p = 0.95; bootstrap p = 0.27 and 0.97), nor under the two readings of O3 (p = 0.094 to 0.86).
fig2b-two-back-absent¶
The task-irrelevant category and age from two trials back are not classified above chance by the two-back classifier.
for z in ("voxel", "map"):
for reading in ("a", "b"):
ca = comp_acc(PRED[(PRED.z == z) & (PRED.space == "two" + reading)])
for comp, name in (("prev_cat", "two-back category"), ("prev_age", "two-back age")):
t, p = ttest(ca[comp], 0.5)
print(f"z by {z}, reading {reading}, {name}: mean {100 * ca[comp].mean():.2f}%, t{len(ca) - 1} = {t:.2f}, one-sided p = {p:.3g}, "
f"bootstrap p = {boot_p(ca[comp], 0.5, seed=12):.4f}")
print(f" current category {100 * ca['cur_cat'].mean():.2f}%, current age {100 * ca['cur_age'].mean():.2f}%")
z by voxel, reading a, two-back category: mean 50.03%, t23 = 0.02, one-sided p = 0.491, bootstrap p = 0.5160 z by voxel, reading a, two-back age: mean 49.05%, t23 = -0.67, one-sided p = 0.744, bootstrap p = 0.7589 current category 49.83%, current age 51.13% z by voxel, reading b, two-back category: mean 49.97%, t23 = -0.02, one-sided p = 0.509, bootstrap p = 0.5298 z by voxel, reading b, two-back age: mean 48.74%, t23 = -0.91, one-sided p = 0.814, bootstrap p = 0.8265 current category 49.78%, current age 50.46% z by map, reading a, two-back category: mean 51.37%, t23 = 1.67, one-sided p = 0.0546, bootstrap p = 0.0422 z by map, reading a, two-back age: mean 50.16%, t23 = 0.14, one-sided p = 0.445, bootstrap p = 0.4445 current category 51.78%, current age 50.55%
z by map, reading b, two-back category: mean 51.16%, t23 = 1.38, one-sided p = 0.0899, bootstrap p = 0.0805 z by map, reading b, two-back age: mean 50.99%, t23 = 0.93, one-sided p = 0.18, bootstrap p = 0.1770 current category 51.47%, current age 49.77%
Result: reproduced. Neither two-back component is above chance in any reading by the t-test (smallest one-sided p = 0.055, the two-back category with maps z-scored by map, reading a). The bootstrap variant gives p = 0.042 for that component (0.081 in reading b), and p > 0.5 with maps z-scored by voxel.
figs2c-8way-decoding¶
Decoding an 8-state space without current age is above chance (16.7%), and its three components are decoded above baseline (53.5%, 52.9%, 55.5%).
t-tests against 12.5% and 50%, one- and two-sided (the side is not stated; O4).
for z in ("voxel", "map"):
a = ACC.loc[(z, "8")]
t1, p1 = ttest(a, 1 / 8)
t2, p2 = ttest(a, 1 / 8, "two-sided")
print(f"z by {z}: 8-way mean {100 * a.mean():.2f}% (SEM {100 * a.sem():.2f}), t{len(a) - 1} = {t1:.2f}, one-sided p = {p1:.3g}, two-sided p = {p2:.3g}")
ca = comp_acc(PRED[(PRED.z == z) & (PRED.space == "8")], ["prev_cat", "prev_age", "cur_cat"],
{"cur_cat": 0, "prev_cat": 1, "prev_age": 2})
for comp in ca:
t1, p1 = ttest(ca[comp], 0.5)
t2, p2 = ttest(ca[comp], 0.5, "two-sided")
print(f" {comp}: {100 * ca[comp].mean():.2f}%, t = {t1:.2f}, one-sided p = {p1:.3g}, two-sided p = {p2:.3g}")
print("expected: 8-way 16.7 +/- 1.67; previous category 53.5, previous age 52.9, current category 55.5 (+/- 10%)")
z by voxel: 8-way mean 15.89% (SEM 1.24), t23 = 2.73, one-sided p = 0.00603, two-sided p = 0.0121
prev_cat: 55.36%, t = 3.12, one-sided p = 0.00243, two-sided p = 0.00486 prev_age: 49.09%, t = -0.45, one-sided p = 0.67, two-sided p = 0.66 cur_cat: 54.48%, t = 2.95, one-sided p = 0.00361, two-sided p = 0.00723 z by map: 8-way mean 14.79% (SEM 1.10), t23 = 2.07, one-sided p = 0.0247, two-sided p = 0.0495 prev_cat: 51.22%, t = 0.87, one-sided p = 0.196, two-sided p = 0.391
prev_age: 50.26%, t = 0.14, one-sided p = 0.443, two-sided p = 0.886 cur_cat: 50.78%, t = 0.55, one-sided p = 0.295, two-sided p = 0.589 expected: 8-way 16.7 +/- 1.67; previous category 53.5, previous age 52.9, current category 55.5 (+/- 10%)
Result: reproduced, resting on the voxel z-scoring. With maps z-scored by voxel, 8-way accuracy is 15.89% (within 15.03 to 18.37%), above 12.5% (one-sided p = 0.006, two-sided p = 0.012), and previous category (55.36%, p = 0.0024) and current category (54.48%, p = 0.0036) are above 50%. Previous age is not (49.09%, p = 0.67), which fits neither the main text ("above baseline") nor the supplement's marginal p = .09. With maps z-scored by map, accuracy is 14.79% (outside tolerance, one-sided p = 0.025) and no component is above 50% (p = 0.20 to 0.44).
figs2a-permutation-control¶
Classifiers trained on permuted labels decode at chance, and the real decoding is above the permutation results for the 16-way classification and each binary component.
10 permutations of the training labels per fold and participant, test labels unchanged, accuracies averaged within participant (Supplemental Results 1.2). Readings: (a) 95th percentile over the participants' averaged permutation accuracies, paired t-test of real against permuted; (b) 95th percentile over all participant and permutation values, compared with the real group mean.
for z in ("voxel", "map"):
q = PERM[PERM.z == z]
real16 = ACC.loc[(z, "16")]
perm_sub = q.groupby(["subject", "perm"]).c.mean()
pa = perm_sub.groupby("subject").mean()
t, p = ttest(pa, 1 / 16)
print(f"z by {z}: permuted 16-way mean {100 * pa.mean():.2f}%, t against 6.25% = {t:.2f}, one-sided p = {p:.3f}")
rows = [("16-way", real16, pa, perm_sub)]
for comp in ["prev_cat", "prev_age", "cur_cat", "cur_age"]:
i = {"cur_cat": 0, "cur_age": 1, "prev_cat": 2, "prev_age": 3}[comp]
cc = (q.true.str[i] == q.pred.str[i]).groupby([q.subject, q.perm]).mean()
rows.append((comp, COMP16[z][comp], cc.groupby("subject").mean(), cc))
for name, real, permavg, pooled in rows:
tt, pp = stats.ttest_rel(real, permavg.reindex(real.index), alternative="greater")
print(f" {name}: real {100 * real.mean():.2f}%, permuted {100 * permavg.mean():.2f}%; (a) 95th pct {100 * np.percentile(permavg, 95):.2f}%, "
f"paired t = {tt:.2f}, one-sided p = {pp:.3g}; (b) 95th pct {100 * np.percentile(pooled, 95):.2f}%, real mean above it: {real.mean() > np.percentile(pooled, 95)}")
print("expected 95th percentiles: 16-way 7.04, previous category 51.5, previous age 51.7, current category 51.4, current age 51.7 (+/- 10%)")
z by voxel: permuted 16-way mean 5.98%, t against 6.25% = -2.03, one-sided p = 0.973 16-way: real 8.93%, permuted 5.98%; (a) 95th pct 7.09%, paired t = 3.47, one-sided p = 0.00104; (b) 95th pct 10.00%, real mean above it: False prev_cat: real 52.03%, permuted 49.71%; (a) 95th pct 51.69%, paired t = 1.92, one-sided p = 0.0334; (b) 95th pct 57.50%, real mean above it: False prev_age: real 49.08%, permuted 49.97%; (a) 95th pct 52.09%, paired t = -0.72, one-sided p = 0.76; (b) 95th pct 57.52%, real mean above it: False cur_cat: real 50.26%, permuted 49.88%; (a) 95th pct 51.84%, paired t = 0.32, one-sided p = 0.374; (b) 95th pct 57.50%, real mean above it: False cur_age: real 50.83%, permuted 49.53%; (a) 95th pct 52.00%, paired t = 0.98, one-sided p = 0.168; (b) 95th pct 57.86%, real mean above it: False z by map: permuted 16-way mean 5.94%, t against 6.25% = -2.02, one-sided p = 0.973
16-way: real 8.92%, permuted 5.94%; (a) 95th pct 6.86%, paired t = 3.96, one-sided p = 0.000314; (b) 95th pct 10.00%, real mean above it: False prev_cat: real 51.78%, permuted 49.88%; (a) 95th pct 51.73%, paired t = 1.96, one-sided p = 0.0311; (b) 95th pct 57.50%, real mean above it: False prev_age: real 51.16%, permuted 49.98%; (a) 95th pct 52.21%, paired t = 1.06, one-sided p = 0.151; (b) 95th pct 57.50%, real mean above it: False cur_cat: real 51.78%, permuted 49.50%; (a) 95th pct 51.36%, paired t = 1.72, one-sided p = 0.0491; (b) 95th pct 57.56%, real mean above it: False cur_age: real 47.92%, permuted 49.63%; (a) 95th pct 52.07%, paired t = -1.37, one-sided p = 0.908; (b) 95th pct 57.50%, real mean above it: False expected 95th percentiles: 16-way 7.04, previous category 51.5, previous age 51.7, current category 51.4, current age 51.7 (+/- 10%)
Result: partially reproduced. Permuted classifiers decode at chance (5.98% and 5.94%, not above 6.25%), and real 16-way accuracy exceeds the permuted one (paired one-sided p = 0.001 and 0.0003). No component meets p < 0.025: previous category comes closest (p = 0.033 and 0.031), then current category with maps z-scored by map (p = 0.049). The 95th percentiles of the participants' averaged permutation accuracies (reading a) are 7.09% for the 16-way and 51.69, 52.09, 51.84 and 52.00% for the components, all within tolerance of the paper's 7.04, 51.5, 51.7, 51.4 and 51.7%; over pooled participant and permutation values (reading b) they are 10.0% and about 57.5%.
Pattern similarity (fig2c-enter-cluster, fig5-similarity-predicts-errors, figs5-euclidean-distance)¶
Methods, State Space Similarities: each run's state maps against all maps of the other runs, the matrices averaged over runs, negative correlations set to NaN. Readings for the maps (O1, O2): raw betas, or z-scored within run by voxel or by map as for classification.
ENTER = np.array([ofc4.trial_class(s) == "Enter" for s in ofc4.STATES])
def between_run(betas, rs, prep, metric):
"""16 x 16 between-run similarity (Pearson r) or distance (Euclidean), averaged over ordered run pairs."""
maps = {r: betas[("16", r)] if prep == "raw" else ofc4.zscore_maps(betas[("16", r)], prep) for r in rs}
mats = []
for r in rs:
for r2 in rs:
if r2 == r:
continue
A, B = maps[r], maps[r2]
if metric == "corr":
Az = (A - A.mean(1, keepdims=True)) / A.std(1, keepdims=True)
Bz = (B - B.mean(1, keepdims=True)) / B.std(1, keepdims=True)
mats.append(Az @ Bz.T / A.shape[1])
else:
mats.append(np.sqrt(((A[:, None, :] - B[None, :, :]) ** 2).sum(-1)))
M = np.mean(mats, axis=0)
if metric == "corr":
M[M < 0] = np.nan
return M
SIM = {(prep, metric, s): between_run(R[s]["betas"], sorted({k[1] for k in R[s]["betas"]}), prep, metric)
for prep in ("raw", "voxel", "map") for metric in ("corr", "dist") for s in SUBJ}
for prep in ("raw", "voxel", "map"):
allc = np.array([SIM[(prep, "corr", s)] for s in SUBJ])
print(f"{prep}: share of NaN (negative) correlations {np.isnan(allc).mean():.3f}; correlations from {np.nanmin(allc):.3f} to {np.nanmax(allc):.3f}")
raw: share of NaN (negative) correlations 0.200; correlations from 0.000 to 0.257 voxel: share of NaN (negative) correlations 0.501; correlations from 0.000 to 0.133 map: share of NaN (negative) correlations 0.200; correlations from 0.000 to 0.257
fig2c-enter-cluster¶
In the similarity structure of OFC state patterns, pairs of states are more similar when both are Enter or both non-Enter states than when they differ.
Per participant the mean within-class (Enter with Enter, non-Enter with non-Enter) against the mean between-class correlation, paired t-test. Readings: the diagonal left out or kept; maps raw or z-scored.
same = ENTER[:, None] == ENTER[None, :]
for prep in ("raw", "voxel", "map"):
for diag in (False, True):
w, b = [], []
for s in SUBJ:
M = SIM[(prep, "corr", s)]
mask_w = same.copy()
if not diag:
np.fill_diagonal(mask_w, False)
w.append(np.nanmean(M[mask_w]))
b.append(np.nanmean(M[~same]))
w, b = np.array(w), np.array(b)
t, p = stats.ttest_rel(w, b)
print(f"{prep}, diagonal {'in' if diag else 'out'}: within {w.mean():.4f}, between {b.mean():.4f}, paired t{len(w) - 1} = {t:.2f}, "
f"two-sided p = {p:.3g}; within > between in {int((w > b).sum())} of {len(w)}")
# dendrogram of the group matrix (shown only; the linkage method is not given, average linkage is a standard choice)
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy.spatial.distance import squareform
G = np.nanmean([SIM[("voxel", "corr", s)] for s in SUBJ], axis=0)
G = (G + G.T) / 2
D = 1 - np.nan_to_num(G, nan=0.0)
np.fill_diagonal(D, 0)
fig, ax = plt.subplots(figsize=(7, 3))
labels_paper = [f"({s[2:]}){s[:2]}" + ("*" if e else "") for s, e in zip(ofc4.STATES, ENTER)]
dendrogram(linkage(squareform(D, checks=False), "average"), labels=labels_paper, ax=ax)
ax.set_ylabel("1 - r (z by voxel)"); ax.set_title("* Enter states")
print("group 1 - r ranges from", round(D[~np.eye(16, dtype=bool)].min(), 3), "to", round(D.max(), 3))
plt.show()
raw, diagonal out: within 0.0355, between 0.0320, paired t23 = 2.82, two-sided p = 0.0097; within > between in 20 of 24 raw, diagonal in: within 0.0365, between 0.0320, paired t23 = 3.50, two-sided p = 0.00194; within > between in 22 of 24 voxel, diagonal out: within 0.0119, between 0.0102, paired t23 = 2.05, two-sided p = 0.0521; within > between in 16 of 24 voxel, diagonal in: within 0.0130, between 0.0102, paired t23 = 3.43, two-sided p = 0.00227; within > between in 19 of 24 map, diagonal out: within 0.0355, between 0.0320, paired t23 = 2.82, two-sided p = 0.0097; within > between in 20 of 24 map, diagonal in: within 0.0365, between 0.0320, paired t23 = 3.50, two-sided p = 0.00194; within > between in 22 of 24 group 1 - r ranges from 0.977 to 0.994
Result: reproduced. Within-class correlations exceed between-class ones (raw or map-z-scored maps, diagonal left out: 0.0355 against 0.0320, t23 = 2.82, p = 0.0097, 20 of 24 participants; diagonal in: p = 0.0019). With maps z-scored by voxel, p = 0.052 without the diagonal and 0.0023 with it. Correlations are small (at most 0.257) and 20% of them (50% after voxel z-scoring) are negative and set to NaN.
fig2d-enter-nonenter-decoding¶
All three hidden components are decoded above chance even when the 16-way classifier is tested separately on Enter and on non-Enter trials.
The predictions of the 16-way classifiers, split by the class of the test map; component accuracy in each subset. Readings: (a) per component the mean of the Enter and non-Enter accuracies; (b) each of the six. One-sided t-tests against 50% (two-sided p shown).
for z, space in [(z, sp) for sp in ("16", "16o3a", "16o3b") for z in ("voxel", "map")]:
q = PRED[(PRED.z == z) & (PRED.space == space)]
q = q.assign(enter=q.true.map(lambda s: ofc4.trial_class(s) == "Enter"))
sub = {k: comp_acc(q[q.enter == v], ["prev_cat", "prev_age", "cur_cat"]) for k, v in (("Enter", True), ("non-Enter", False))}
avg = (sub["Enter"] + sub["non-Enter"]) / 2
for comp in avg:
t, p = ttest(avg[comp], 0.5)
print(f"{space}, z by {z}, (a) {comp}: {100 * avg[comp].mean():.2f}%, t = {t:.2f}, one-sided p = {p:.3g}, two-sided p = {min(1, 2 * min(p, 1 - p)):.3g}")
for k, ca in sub.items():
for comp in ca:
t, p = ttest(ca[comp], 0.5)
print(f" (b) {k} {comp}: {100 * ca[comp].mean():.2f}%, t = {t:.2f}, one-sided p = {p:.3g}")
16, z by voxel, (a) prev_cat: 52.03%, t = 1.74, one-sided p = 0.0476, two-sided p = 0.0953 16, z by voxel, (a) prev_age: 49.08%, t = -0.74, one-sided p = 0.765, two-sided p = 0.469 16, z by voxel, (a) cur_cat: 50.26%, t = 0.22, one-sided p = 0.415, two-sided p = 0.829 (b) Enter prev_cat: 52.34%, t = 1.43, one-sided p = 0.0826 (b) Enter prev_age: 48.83%, t = -0.74, one-sided p = 0.766
(b) Enter cur_cat: 48.49%, t = -1.14, one-sided p = 0.867 (b) non-Enter prev_cat: 51.72%, t = 0.99, one-sided p = 0.167 (b) non-Enter prev_age: 49.32%, t = -0.37, one-sided p = 0.642 (b) non-Enter cur_cat: 52.03%, t = 1.17, one-sided p = 0.128 16, z by map, (a) prev_cat: 51.78%, t = 1.86, one-sided p = 0.038, two-sided p = 0.076 16, z by map, (a) prev_age: 51.16%, t = 1.01, one-sided p = 0.161, two-sided p = 0.323 16, z by map, (a) cur_cat: 51.78%, t = 1.32, one-sided p = 0.0997, two-sided p = 0.199 (b) Enter prev_cat: 50.99%, t = 0.74, one-sided p = 0.232 (b) Enter prev_age: 52.14%, t = 1.29, one-sided p = 0.105 (b) Enter cur_cat: 51.33%, t = 1.01, one-sided p = 0.162 (b) non-Enter prev_cat: 52.58%, t = 1.75, one-sided p = 0.0467 (b) non-Enter prev_age: 50.18%, t = 0.14, one-sided p = 0.446 (b) non-Enter cur_cat: 52.24%, t = 1.22, one-sided p = 0.118 16o3a, z by voxel, (a) prev_cat: 51.63%, t = 1.23, one-sided p = 0.116, two-sided p = 0.231 16o3a, z by voxel, (a) prev_age: 49.80%, t = -0.18, one-sided p = 0.569, two-sided p = 0.862 16o3a, z by voxel, (a) cur_cat: 50.51%, t = 0.54, one-sided p = 0.296, two-sided p = 0.592 (b) Enter prev_cat: 51.64%, t = 1.14, one-sided p = 0.132 (b) Enter prev_age: 49.17%, t = -0.59, one-sided p = 0.721 (b) Enter cur_cat: 49.92%, t = -0.06, one-sided p = 0.523 (b) non-Enter prev_cat: 51.61%, t = 0.93, one-sided p = 0.181 (b) non-Enter prev_age: 50.44%, t = 0.25, one-sided p = 0.404 (b) non-Enter cur_cat: 51.09%, t = 0.65, one-sided p = 0.262 16o3a, z by map, (a) prev_cat: 52.11%, t = 2.38, one-sided p = 0.0129, two-sided p = 0.0258 16o3a, z by map, (a) prev_age: 49.71%, t = -0.28, one-sided p = 0.609, two-sided p = 0.781 16o3a, z by map, (a) cur_cat: 51.05%, t = 1.02, one-sided p = 0.16, two-sided p = 0.32 (b) Enter prev_cat: 52.16%, t = 1.48, one-sided p = 0.0761 (b) Enter prev_age: 49.95%, t = -0.05, one-sided p = 0.518 (b) Enter cur_cat: 50.42%, t = 0.29, one-sided p = 0.386 (b) non-Enter prev_cat: 52.06%, t = 1.67, one-sided p = 0.0541 (b) non-Enter prev_age: 49.48%, t = -0.37, one-sided p = 0.642 (b) non-Enter cur_cat: 51.69%, t = 1.07, one-sided p = 0.147 16o3b, z by voxel, (a) prev_cat: 52.60%, t = 2.00, one-sided p = 0.0288, two-sided p = 0.0576 16o3b, z by voxel, (a) prev_age: 49.19%, t = -0.70, one-sided p = 0.756, two-sided p = 0.488 16o3b, z by voxel, (a) cur_cat: 50.99%, t = 0.90, one-sided p = 0.188, two-sided p = 0.376 (b) Enter prev_cat: 52.24%, t = 1.32, one-sided p = 0.0994 (b) Enter prev_age: 48.15%, t = -1.21, one-sided p = 0.881 (b) Enter cur_cat: 49.64%, t = -0.23, one-sided p = 0.591 (b) non-Enter prev_cat: 52.97%, t = 1.77, one-sided p = 0.0447 (b) non-Enter prev_age: 50.23%, t = 0.13, one-sided p = 0.448 (b) non-Enter cur_cat: 52.34%, t = 1.49, one-sided p = 0.0747 16o3b, z by map, (a) prev_cat: 51.99%, t = 1.92, one-sided p = 0.0339, two-sided p = 0.0678 16o3b, z by map, (a) prev_age: 49.97%, t = -0.02, one-sided p = 0.509, two-sided p = 0.983 16o3b, z by map, (a) cur_cat: 51.51%, t = 1.31, one-sided p = 0.102, two-sided p = 0.204 (b) Enter prev_cat: 51.43%, t = 0.94, one-sided p = 0.18 (b) Enter prev_age: 50.78%, t = 0.51, one-sided p = 0.309 (b) Enter cur_cat: 51.43%, t = 0.96, one-sided p = 0.174 (b) non-Enter prev_cat: 52.55%, t = 1.52, one-sided p = 0.0716 (b) non-Enter prev_age: 49.17%, t = -0.62, one-sided p = 0.731 (b) non-Enter cur_cat: 51.59%, t = 1.07, one-sided p = 0.149
Result: not reproduced. Tested separately on Enter and non-Enter maps, only the previous category reaches p < 0.05 in reading (a) (52.03%, one-sided p = 0.048; z by map 51.78%, p = 0.038); previous age and current category do not (p >= 0.10). In reading (b) at most one of six reaches p < 0.05. Under the two readings of O3 the same holds: in reading (a) only previous category reaches p < 0.05 (three of four), and in reading (b) at most one of six.
fig4a-decoding-error-correlation¶
Across participants, mean classification accuracy of the three hidden components in OFC correlates negatively with error rate.
The mean of the three hidden-component accuracies (fig2b-hidden-components) against the error rate (behavior-error-rate, reading b, the error blocks over the 96 trials after the first of each run); Pearson correlation and bootstrap 95% CI over participants (percentile intervals, 10000 and 1000 resamples).
for z in ("voxel", "map"):
hid = COMP16[z][["prev_cat", "prev_age", "cur_cat"]].mean(axis=1)
for name, err in (("reading b", ERRRATE), ("reading a", er.rate_a / 100)):
x, y = hid.reindex(SUBJ).values, err.reindex(SUBJ).values
r, p = stats.pearsonr(x, y)
line = f"z by {z}, error rate {name}: r = {r:.3f}, t{len(x) - 2} = {r * np.sqrt((len(x) - 2) / (1 - r ** 2)):.2f}, p = {p:.3g}"
for nb, seed in ((10000, 21), (1000, 22)):
rng = np.random.default_rng(seed)
bs = []
for _ in range(nb):
i = rng.integers(0, len(x), len(x))
if np.std(x[i]) > 0 and np.std(y[i]) > 0:
bs.append(np.corrcoef(x[i], y[i])[0, 1])
lo, hi = np.percentile(bs, [2.5, 97.5])
line += f"; {nb} resamples CI [{lo:.3f}, {hi:.3f}]"
print(line)
sp = stats.spearmanr(hid.reindex(SUBJ), ERRRATE.reindex(SUBJ))
print(f" Spearman rho = {sp.statistic:.3f}, p = {sp.pvalue:.3g} (variant; participant 9 has the highest error rate, {100 * ERRRATE.max():.2f}%)")
print("expected r = -0.58 +/- 0.2632 (paper CI -0.8432 to -0.0169)")
fig, ax = plt.subplots(figsize=(4, 3))
hid = COMP16["voxel"][["prev_cat", "prev_age", "cur_cat"]].mean(axis=1)
ax.scatter(100 * ERRRATE.reindex(SUBJ), 100 * hid.reindex(SUBJ))
ax.set_xlabel("error rate (%)"); ax.set_ylabel("mean hidden-component accuracy (%), z by voxel")
print("plotted ranges: error rate", round(100 * ERRRATE.min(), 2), "to", round(100 * ERRRATE.max(), 2),
"; accuracy", round(100 * hid.min(), 2), "to", round(100 * hid.max(), 2))
plt.show()
z by voxel, error rate reading b: r = -0.268, t22 = -1.30, p = 0.205; 10000 resamples CI [-0.717, 0.338]; 1000 resamples CI [-0.700, 0.327]
z by voxel, error rate reading a: r = -0.259, t22 = -1.26, p = 0.221; 10000 resamples CI [-0.715, 0.337]; 1000 resamples CI [-0.700, 0.328] Spearman rho = -0.206, p = 0.333 (variant; participant 9 has the highest error rate, 7.50%)
z by map, error rate reading b: r = -0.004, t22 = -0.02, p = 0.984; 10000 resamples CI [-0.453, 0.359]; 1000 resamples CI [-0.446, 0.377]
z by map, error rate reading a: r = -0.028, t22 = -0.13, p = 0.897; 10000 resamples CI [-0.471, 0.342]; 1000 resamples CI [-0.464, 0.348] Spearman rho = -0.125, p = 0.56 (variant; participant 9 has the highest error rate, 7.50%) expected r = -0.58 +/- 0.2632 (paper CI -0.8432 to -0.0169) plotted ranges: error rate 0.21 to 7.5 ; accuracy 44.17 to 55.0
Result: partially reproduced, resting on the voxel z-scoring. The correlation is negative but weak and not significant (r = -0.268, p = 0.21, bootstrap CI -0.717 to 0.338); it lies inside the paper's CI (-0.843 to -0.017) but outside the tolerance (-0.843 to -0.317). With maps z-scored by map, r = -0.004, outside the paper's CI. The error rate counts one error per block, a lower bound.
Trial-wise decoding around errors (fig4b-error-trial-drop, fig4b-decline-before-errors, fig4b-error-vs-matched)¶
The 16-way classifier of each fold tested on single volumes of the held-out run, z-scored per voxel over the run (Methods), unsmoothed (no voxel coordinates). Readings: the volume (first at or after onset + 4.8 s, or nearest), the z-scoring of the training maps, and the error trial's state (strict: runs where every block holds one designed trial; span: Enter errors whose block spans two trials). An error enters when its five preceding trials are labelled trials of the same run (Methods, Participants); trial +1 is the first labelled trial after the block. Matched courses: (a) every non-overlapping run of seven consecutive labelled trials, (b) one such run per error, drawn from the error's run outside its course.
def courses(z, vol, err_reading, seed):
rng = np.random.default_rng(seed)
q = TRIAL[(TRIAL.z == z) & (TRIAL.vol == vol)].copy()
q["truth"] = q["true"] if err_reading == "strict" else q["true_span"]
q["c"] = (q.truth == q.pred).astype(float)
q["cls"] = q.truth.map(lambda x: ofc4.trial_class(x) if isinstance(x, str) else None)
err_rows, mat_a, mat_b = [], [], []
for (s, r), g in q.groupby(["subject", "run"]):
g = g.set_index("pos").sort_index()
lab_pos = set(g.index[g.cond.isin(ofc4.STATES)])
used = set()
for p, e in g[g.cond == "err"].iterrows():
if pd.isna(e.truth) or not all((p - k) in lab_pos for k in range(1, 6)):
continue
eid = f"{s}-{r}-{p}"
for k in range(-5, 2):
pp = p + k
if k == 0:
err_rows.append({"subject": s, "eid": eid, "rel": 0, "c": e.c, "cls": e.cls})
elif pp in lab_pos:
err_rows.append({"subject": s, "eid": eid, "rel": k, "c": g.loc[pp, "c"], "cls": g.loc[pp, "cls"]})
course = set(range(p - 5, p + 2))
starts = [st for st in sorted(lab_pos) if all((st + k) in lab_pos for k in range(7)) and not (set(range(st, st + 7)) & course)]
if starts:
st = starts[rng.integers(len(starts))]
for k in range(7):
mat_b.append({"subject": s, "eid": eid, "rel": k - 5, "c": g.loc[st + k, "c"]})
st = min(lab_pos) if lab_pos else None
positions = sorted(lab_pos)
i = 0
while i < len(positions):
st = positions[i]
if all((st + k) in lab_pos for k in range(7)):
for k in range(7):
mat_a.append({"subject": s, "eid": f"{s}-{r}-m{st}", "rel": k - 5, "c": g.loc[st + k, "c"], "cls": g.loc[st + k, "cls"]})
i += 7
else:
i += 1
return pd.DataFrame(err_rows), pd.DataFrame(mat_a), pd.DataFrame(mat_b)
COURSES = {(z, vol, er_): courses(z, vol, er_, seed=31) for z in ("voxel", "map") for vol in ("after", "nearest") for er_ in ("strict", "span")}
for key, (e, ma, mb) in COURSES.items():
print(key, f"errors {e.eid.nunique()} from {e.subject.nunique()} participants, by class {e[e.rel == 0].cls.value_counts().to_dict()}; "
f"matched (a) {ma.eid.nunique()} courses, (b) {mb.eid.nunique()}; (a) with an Enter trial at position 0 {ma[(ma.rel == 0) & (ma.cls == 'Enter')].eid.nunique()}")
('voxel', 'after', 'strict') errors 37 from 16 participants, by class {'Enter': 37}; matched (a) 1427 courses, (b) 37; (a) with an Enter trial at position 0 509
('voxel', 'after', 'span') errors 81 from 21 participants, by class {'Enter': 81}; matched (a) 1427 courses, (b) 81; (a) with an Enter trial at position 0 509
('voxel', 'nearest', 'strict') errors 37 from 16 participants, by class {'Enter': 37}; matched (a) 1429 courses, (b) 37; (a) with an Enter trial at position 0 509
('voxel', 'nearest', 'span') errors 81 from 21 participants, by class {'Enter': 81}; matched (a) 1429 courses, (b) 81; (a) with an Enter trial at position 0 509
('map', 'after', 'strict') errors 37 from 16 participants, by class {'Enter': 37}; matched (a) 1427 courses, (b) 37; (a) with an Enter trial at position 0 509
('map', 'after', 'span') errors 81 from 21 participants, by class {'Enter': 81}; matched (a) 1427 courses, (b) 81; (a) with an Enter trial at position 0 509
('map', 'nearest', 'strict') errors 37 from 16 participants, by class {'Enter': 37}; matched (a) 1429 courses, (b) 37; (a) with an Enter trial at position 0 509
('map', 'nearest', 'span') errors 81 from 21 participants, by class {'Enter': 81}; matched (a) 1429 courses, (b) 81; (a) with an Enter trial at position 0 509
fig4b-error-trial-drop¶
Trial-wise 16-way decoding is lower on error trials (4.2%) than in the five trials preceding an error (11%).
for key, (e, ma, mb) in COURSES.items():
pre = e[e.rel.between(-5, -1)].groupby("subject").c.mean()
at = e[e.rel == 0].groupby("subject").c.mean()
t, p = stats.ttest_rel(pre, at.reindex(pre.index))
print(f"{key}: n = {len(pre)} participants, {e.eid.nunique()} errors; five preceding {100 * pre.mean():.2f}%, error trials {100 * at.mean():.2f}%; "
f"paired t{len(pre) - 1} = {t:.2f}, two-sided p = {p:.3g}")
print("expected: preceding 11 +/- 1.1, error 4.2 +/- 0.42")
lab_trials = TRIAL[(TRIAL.vol == "after") & TRIAL.cond.isin(ofc4.STATES)]
print("all labelled trials, volume 'after':", (100 * (lab_trials.true == lab_trials.pred).groupby(lab_trials.z).mean()).round(2).to_dict())
('voxel', 'after', 'strict'): n = 16 participants, 37 errors; five preceding 8.21%, error trials 2.08%; paired t15 = 1.41, two-sided p = 0.18
('voxel', 'after', 'span'): n = 21 participants, 81 errors; five preceding 8.02%, error trials 1.59%; paired t20 = 1.99, two-sided p = 0.0605
('voxel', 'nearest', 'strict'): n = 16 participants, 37 errors; five preceding 9.04%, error trials 2.08%; paired t15 = 1.62, two-sided p = 0.126
('voxel', 'nearest', 'span'): n = 21 participants, 81 errors; five preceding 8.08%, error trials 1.59%; paired t20 = 1.98, two-sided p = 0.0613
('map', 'after', 'strict'): n = 16 participants, 37 errors; five preceding 4.42%, error trials 9.38%; paired t15 = -1.17, two-sided p = 0.262
('map', 'after', 'span'): n = 21 participants, 81 errors; five preceding 5.67%, error trials 7.70%; paired t20 = -0.70, two-sided p = 0.495
('map', 'nearest', 'strict'): n = 16 participants, 37 errors; five preceding 3.54%, error trials 10.62%; paired t15 = -1.45, two-sided p = 0.169
('map', 'nearest', 'span'): n = 21 participants, 81 errors; five preceding 4.47%, error trials 10.17%; paired t20 = -1.66, two-sided p = 0.112
expected: preceding 11 +/- 1.1, error 4.2 +/- 0.42
all labelled trials, volume 'after': {'map': 7.04, 'voxel': 9.11}
Result: not reproduced. Error trials decode worse than the five trials before them with maps z-scored by voxel (1.59% against 8.02%, 21 participants, 81 errors), but not significantly (p = 0.061; 0.18 with the 37 strictly reconstructed errors of 16 participants), and both values are outside tolerance (expected 4.2% and 11%). With maps z-scored by map, error trials decode better than the preceding ones (7.70% against 5.67%). Only error trials whose state can be reconstructed enter, and all are Enter trials, which decode at 0.90% (z by voxel) on correct trials too, as that classifier rarely names an Enter state (data check), so the comparison with the mixed-class trials before them is confounded by class.
fig4b-decline-before-errors¶
Decoding accuracy declines linearly over the trials leading up to and including an error, and no such decline is seen in matched time courses without errors.
Linear mixed models of accuracy on trial position (-5 to 0), random intercept and slope by participant, Wald chi-square; accuracy entered per trial (0 or 1) or as the mean per participant and position; the interaction from a joint model of error and matched courses.
def slope_test(df, formula, term, re="~rel"):
"""Random intercept and slope; when that fit does not converge, the random-intercept model
is used and the result is marked."""
try:
m = smf.mixedlm(formula, df, groups=df["subject"], re_formula=re).fit()
slope_ok = bool(m.converged)
except np.linalg.LinAlgError: # singular fit: counted as not converged
slope_ok = False
if not slope_ok:
try:
m = smf.mixedlm(formula, df, groups=df["subject"]).fit()
except np.linalg.LinAlgError:
print(" a random-intercept fit was singular; its result is left out:", formula, f"({len(df)} rows)")
return np.nan, np.nan, np.nan, False
chi2, _, p = wald(m, [term])
return m.params[term], chi2, p, slope_ok
DECL = {}
for key, (e, ma, mb) in COURSES.items():
for mname, mdf in (("a", ma), ("b", mb)):
for unit in ("trial", "mean"):
ee = e[e.rel <= 0].assign(cond="error")
mm = mdf[mdf.rel <= 0].assign(cond="matched")
if unit == "mean":
ee = ee.groupby(["subject", "rel", "cond"]).c.mean().reset_index()
mm = mm.groupby(["subject", "rel", "cond"]).c.mean().reset_index()
b1, c1, p1, v1 = slope_test(ee, "c ~ rel", "rel")
b2, c2, p2, v2 = slope_test(mm, "c ~ rel", "rel")
j = pd.concat([ee, mm])
j["is_err"] = (j.cond == "error").astype(int)
b3, c3, p3, v3 = slope_test(j, "c ~ rel * is_err", "rel:is_err")
DECL[(key, mname, unit)] = (b1, p1, b3, p3, b2, p2)
print(f"{key}, matched {mname}, {unit}{'' if (v1 and v2 and v3) else ' (random intercept only in a model)'}: error slope {100 * b1:.2f} pts/trial chi2 = {c1:.2f} p = {p1:.3g} | "
f"interaction {100 * b3:.2f} chi2 = {c3:.2f} p = {p3:.3g} | matched slope {100 * b2:.2f} chi2 = {c2:.2f} p = {p2:.3g}")
print("paper: error slope chi2(1) = 4.7, p = 0.03; interaction chi2(1) = 4.3, p = 0.04; non-error slope p = 0.33")
e, ma, mb = COURSES[("voxel", "after", "span")]
fig, ax = plt.subplots(figsize=(5, 3))
for lbl, df in (("error courses", e), ("matched (a)", ma)):
g = df.groupby(["subject", "rel"]).c.mean().unstack()
ax.errorbar(g.columns, 100 * g.mean(), 100 * g.sem(), label=lbl, capsize=3, marker="o")
print(lbl, "mean accuracy (%) by position:", (100 * g.mean()).round(2).to_dict())
ax.axhline(100 / 16, ls="--", c="k"); ax.set_xlabel("trial relative to error"); ax.set_ylabel("16-way accuracy (%)"); ax.legend()
plt.show()
('voxel', 'after', 'strict'), matched a, trial (random intercept only in a model): error slope -0.46 pts/trial chi2 = 0.24 p = 0.622 | interaction -0.46 chi2 = 0.16 p = 0.689 | matched slope -0.03 chi2 = 0.02 p = 0.875
('voxel', 'after', 'strict'), matched a, mean (random intercept only in a model): error slope -0.42 pts/trial chi2 = 0.15 p = 0.699 | interaction -0.37 chi2 = 0.14 p = 0.705 | matched slope -0.05 chi2 = 0.07 p = 0.798
('voxel', 'after', 'strict'), matched b, trial (random intercept only in a model): error slope -0.46 pts/trial chi2 = 0.24 p = 0.622 | interaction -0.08 chi2 = 0.00 p = 0.961 | matched slope -0.39 chi2 = 0.10 p = 0.752
('voxel', 'after', 'strict'), matched b, mean (random intercept only in a model): error slope -0.42 pts/trial chi2 = 0.15 p = 0.699 | interaction 0.89 chi2 = 0.23 p = 0.633 | matched slope -1.31 chi2 = 0.78 p = 0.378
('voxel', 'after', 'span'), matched a, trial (random intercept only in a model): error slope -0.63 pts/trial chi2 = 0.98 p = 0.321 | interaction -0.57 chi2 = 0.52 p = 0.473 | matched slope -0.03 chi2 = 0.02 p = 0.875
('voxel', 'after', 'span'), matched a, mean (random intercept only in a model): error slope -0.24 pts/trial chi2 = 0.08 p = 0.776 | interaction -0.19 chi2 = 0.05 p = 0.828 | matched slope -0.05 chi2 = 0.07 p = 0.798
('voxel', 'after', 'span'), matched b, trial (random intercept only in a model): error slope -0.63 pts/trial chi2 = 0.98 p = 0.321 | interaction -0.35 chi2 = 0.15 p = 0.701 | matched slope -0.40 chi2 = 0.21 p = 0.65
('voxel', 'after', 'span'), matched b, mean (random intercept only in a model): error slope -0.24 pts/trial chi2 = 0.08 p = 0.776 | interaction 0.80 chi2 = 0.43 p = 0.51 | matched slope -1.04 chi2 = 1.86 p = 0.173
('voxel', 'nearest', 'strict'), matched a, trial (random intercept only in a model): error slope -0.62 pts/trial chi2 = 0.42 p = 0.519 | interaction -0.65 chi2 = 0.31 p = 0.58 | matched slope 0.00 chi2 = 0.00 p = 0.983
('voxel', 'nearest', 'strict'), matched a, mean (random intercept only in a model): error slope -1.01 pts/trial chi2 = 0.66 p = 0.418 | interaction -0.98 chi2 = 0.81 p = 0.367 | matched slope -0.03 chi2 = 0.02 p = 0.899
('voxel', 'nearest', 'strict'), matched b, trial (random intercept only in a model): error slope -0.62 pts/trial chi2 = 0.42 p = 0.519 | interaction -1.47 chi2 = 0.74 p = 0.391 | matched slope 0.74 chi2 = 0.22 p = 0.642
('voxel', 'nearest', 'strict'), matched b, mean (random intercept only in a model): error slope -1.01 pts/trial chi2 = 0.66 p = 0.418 | interaction -0.95 chi2 = 0.21 p = 0.645 | matched slope -0.05 chi2 = 0.00 p = 0.973
('voxel', 'nearest', 'span'), matched a, trial (random intercept only in a model): error slope -0.42 pts/trial chi2 = 0.50 p = 0.479 | interaction -0.39 chi2 = 0.24 p = 0.624 | matched slope 0.00 chi2 = 0.00 p = 0.983
('voxel', 'nearest', 'span'), matched a, mean: error slope -0.87 pts/trial chi2 = 0.80 p = 0.372 | interaction -0.85 chi2 = 0.81 p = 0.368 | matched slope -0.03 chi2 = 0.02 p = 0.899
/tmp/ipykernel_2502/3304350975.py:11: SingularMatrixWarning: The random effects covariance matrix is singular. m = smf.mixedlm(formula, df, groups=df["subject"]).fit()
('voxel', 'nearest', 'span'), matched b, trial (random intercept only in a model): error slope -0.42 pts/trial chi2 = 0.50 p = 0.479 | interaction 0.04 chi2 = 0.00 p = 0.97 | matched slope -0.46 chi2 = 0.39 p = 0.534
('voxel', 'nearest', 'span'), matched b, mean (random intercept only in a model): error slope -0.87 pts/trial chi2 = 0.80 p = 0.372 | interaction 0.52 chi2 = 0.15 p = 0.699 | matched slope -1.39 chi2 = 2.15 p = 0.143
('map', 'after', 'strict'), matched a, trial (random intercept only in a model): error slope 1.39 pts/trial chi2 = 2.45 p = 0.117 | interaction 1.59 chi2 = 2.41 p = 0.121 | matched slope -0.08 chi2 = 0.27 p = 0.606
('map', 'after', 'strict'), matched a, mean (random intercept only in a model): error slope 1.30 pts/trial chi2 = 1.83 p = 0.176 | interaction 1.39 chi2 = 3.22 p = 0.0728 | matched slope -0.09 chi2 = 0.33 p = 0.566
('map', 'after', 'strict'), matched b, trial (random intercept only in a model): error slope 1.39 pts/trial chi2 = 2.45 p = 0.117 | interaction 1.47 chi2 = 1.22 p = 0.27 | matched slope -0.08 chi2 = 0.01 p = 0.937
('map', 'after', 'strict'), matched b, mean (random intercept only in a model): error slope 1.30 pts/trial chi2 = 1.83 p = 0.176 | interaction 1.21 chi2 = 0.74 p = 0.391 | matched slope 0.09 chi2 = 0.01 p = 0.936
('map', 'after', 'span'), matched a, trial (random intercept only in a model): error slope 1.16 pts/trial chi2 = 3.07 p = 0.08 | interaction 1.25 chi2 = 3.19 p = 0.0739 | matched slope -0.08 chi2 = 0.27 p = 0.606
('map', 'after', 'span'), matched a, mean (random intercept only in a model): error slope 1.21 pts/trial chi2 = 2.79 p = 0.0951 | interaction 1.30 chi2 = 3.53 p = 0.0603 | matched slope -0.09 chi2 = 0.33 p = 0.566
/tmp/ipykernel_2502/3304350975.py:11: SingularMatrixWarning: The random effects covariance matrix is singular. m = smf.mixedlm(formula, df, groups=df["subject"]).fit()
('map', 'after', 'span'), matched b, trial (random intercept only in a model): error slope 1.16 pts/trial chi2 = 3.07 p = 0.08 | interaction 1.69 chi2 = 2.96 p = 0.0852 | matched slope -0.53 chi2 = 0.55 p = 0.457
('map', 'after', 'span'), matched b, mean (random intercept only in a model): error slope 1.21 pts/trial chi2 = 2.79 p = 0.0951 | interaction 1.86 chi2 = 3.16 p = 0.0756 | matched slope -0.65 chi2 = 0.73 p = 0.393
('map', 'nearest', 'strict'), matched a, trial (random intercept only in a model): error slope 2.08 pts/trial chi2 = 6.09 p = 0.0136 | interaction 2.04 chi2 = 4.08 p = 0.0434 | matched slope -0.02 chi2 = 0.01 p = 0.93
('map', 'nearest', 'strict'), matched a, mean (random intercept only in a model): error slope 1.85 pts/trial chi2 = 3.63 p = 0.0567 | interaction 1.88 chi2 = 6.48 p = 0.0109 | matched slope -0.03 chi2 = 0.02 p = 0.879
('map', 'nearest', 'strict'), matched b, trial (random intercept only in a model): error slope 2.08 pts/trial chi2 = 6.09 p = 0.0136 | interaction 3.01 chi2 = 5.23 p = 0.0222 | matched slope -0.83 chi2 = 0.59 p = 0.441
('map', 'nearest', 'strict'), matched b, mean (random intercept only in a model): error slope 1.85 pts/trial chi2 = 3.63 p = 0.0567 | interaction 2.18 chi2 = 2.70 p = 0.101 | matched slope -0.33 chi2 = 0.11 p = 0.736
('map', 'nearest', 'span'), matched a, trial: error slope 1.47 pts/trial chi2 = 4.24 p = 0.0394 | interaction 1.51 chi2 = 4.70 p = 0.0302 | matched slope -0.02 chi2 = 0.01 p = 0.93
a random-intercept fit was singular; its result is left out: c ~ rel * is_err (270 rows)
('map', 'nearest', 'span'), matched a, mean (random intercept only in a model): error slope 1.49 pts/trial chi2 = 3.81 p = 0.0509 | interaction nan chi2 = nan p = nan | matched slope -0.03 chi2 = 0.02 p = 0.879
('map', 'nearest', 'span'), matched b, trial (random intercept only in a model): error slope 1.47 pts/trial chi2 = 4.24 p = 0.0394 | interaction 2.26 chi2 = 5.20 p = 0.0226 | matched slope -0.69 chi2 = 0.64 p = 0.425
('map', 'nearest', 'span'), matched b, mean (random intercept only in a model): error slope 1.49 pts/trial chi2 = 3.81 p = 0.0509 | interaction 2.23 chi2 = 4.46 p = 0.0347 | matched slope -0.73 chi2 = 0.89 p = 0.346
paper: error slope chi2(1) = 4.7, p = 0.03; interaction chi2(1) = 4.3, p = 0.04; non-error slope p = 0.33
error courses mean accuracy (%) by position: {-5: 4.44, -4: 6.43, -3: 13.65, -2: 3.97, -1: 11.59, 0: 1.59, 1: 13.06}
matched (a) mean accuracy (%) by position: {-5: 10.83, -4: 7.23, -3: 10.17, -2: 8.93, -1: 9.37, 0: 9.4, 1: 8.13}
Result: not reproduced. In no reading is the slope over trials -5 to 0 negative with p < 0.05 together with a significant interaction. With maps z-scored by voxel, the error slope is negative but not significant (smallest p = 0.32) and the interaction is not significant. With maps z-scored by map and the nearest volume, the slope and interaction are significant but positive (accuracy rises towards the error; span reading, matched courses (a), slope 1.47 percentage points per trial, p = 0.039; interaction p = 0.030; strict reading, slope 2.08, p = 0.014). Most random-slope fits did not converge, and the random-intercept models are reported in their place, as marked; one singular fit is left out. Position 0 is always an Enter trial and position -1 an Exit trial, so the course also reflects the decoding of each class (data check). See the printed table for every reading.
fig4b-error-vs-matched¶
Decoding on the error trial is lower than on matched correct trials, and marginally lower on the trial before the error.
for key, (e, ma, mb) in COURSES.items():
for mname, mdf in (("a", ma), ("b", mb)):
out = []
for rel in (0, -1):
x = e[e.rel == rel].groupby("subject").c.mean()
y = mdf[mdf.rel == rel].groupby("subject").c.mean().reindex(x.index)
ok = y.notna()
t, p1 = stats.ttest_rel(x[ok], y[ok], alternative="less")
_, p2 = stats.ttest_rel(x[ok], y[ok])
out.append(f"trial {rel}: error {100 * x[ok].mean():.2f}% vs matched {100 * y[ok].mean():.2f}%, t{int(ok.sum()) - 1} = {t:.2f}, one-sided p = {p1:.3g}, two-sided p = {p2:.3g}")
print(key, f"matched {mname}:", " | ".join(out))
# Reading (c): every error trial here is an Enter trial, so position 0 is compared with matched courses (a)
# whose position 0 is an Enter trial, and with all labelled Enter trials of the participant.
for key, (e, ma, mb) in COURSES.items():
x = e[e.rel == 0].groupby("subject").c.mean()
y = ma[(ma.rel == 0) & (ma.cls == "Enter")].groupby("subject").c.mean().reindex(x.index)
ok = y.notna()
t, p1 = stats.ttest_rel(x[ok], y[ok], alternative="less")
q = TRIAL[(TRIAL.z == key[0]) & (TRIAL.vol == key[1]) & TRIAL.cond.isin(ofc4.STATES)]
en = (q.true == q.pred)[q.cond.map(ofc4.trial_class) == "Enter"].groupby(q.subject).mean().reindex(x.index)
t2, p2 = stats.ttest_rel(x, en, alternative="less")
print(key, f"(c) trial 0: error {100 * x[ok].mean():.2f}% vs matched with Enter at 0 {100 * y[ok].mean():.2f}%, t{int(ok.sum()) - 1} = {t:.2f}, one-sided p = {p1:.3g} | "
f"error {100 * x.mean():.2f}% vs all labelled Enter trials {100 * en.mean():.2f}%, t{len(x) - 1} = {t2:.2f}, one-sided p = {p2:.3g}")
print("paper: trial 0 t24 = 2.3, p = 0.01; trial -1 t23 = 1.5, p = 0.07")
('voxel', 'after', 'strict') matched a: trial 0: error 2.08% vs matched 8.93%, t15 = -2.95, one-sided p = 0.005, two-sided p = 0.00999 | trial -1: error 13.54% vs matched 10.37%, t15 = 0.47, one-sided p = 0.677, two-sided p = 0.645
('voxel', 'after', 'strict') matched b: trial 0: error 2.08% vs matched 8.33%, t15 = -1.14, one-sided p = 0.135, two-sided p = 0.27 | trial -1: error 13.54% vs matched 16.67%, t15 = -0.34, one-sided p = 0.369, two-sided p = 0.738
('voxel', 'after', 'span') matched a: trial 0: error 1.59% vs matched 9.21%, t20 = -4.07, one-sided p = 0.000296, two-sided p = 0.000591 | trial -1: error 11.59% vs matched 8.81%, t20 = 0.54, one-sided p = 0.703, two-sided p = 0.593
('voxel', 'after', 'span') matched b: trial 0: error 1.59% vs matched 3.57%, t20 = -1.42, one-sided p = 0.0855, two-sided p = 0.171 | trial -1: error 11.59% vs matched 6.43%, t20 = 0.90, one-sided p = 0.81, two-sided p = 0.379
('voxel', 'nearest', 'strict') matched a: trial 0: error 2.08% vs matched 10.73%, t15 = -3.53, one-sided p = 0.0015, two-sided p = 0.00301 | trial -1: error 11.46% vs matched 9.85%, t15 = 0.24, one-sided p = 0.594, two-sided p = 0.812
('voxel', 'nearest', 'strict') matched b: trial 0: error 2.08% vs matched 8.33%, t15 = -1.14, one-sided p = 0.135, two-sided p = 0.27 | trial -1: error 11.46% vs matched 13.54%, t15 = -0.21, one-sided p = 0.419, two-sided p = 0.837
('voxel', 'nearest', 'span') matched a: trial 0: error 1.59% vs matched 9.69%, t20 = -4.28, one-sided p = 0.000185, two-sided p = 0.000369 | trial -1: error 9.68% vs matched 9.20%, t20 = 0.09, one-sided p = 0.537, two-sided p = 0.926
('voxel', 'nearest', 'span') matched b: trial 0: error 1.59% vs matched 4.76%, t20 = -1.79, one-sided p = 0.044, two-sided p = 0.0881 | trial -1: error 9.68% vs matched 8.02%, t20 = 0.26, one-sided p = 0.6, two-sided p = 0.8
('map', 'after', 'strict') matched a: trial 0: error 9.38% vs matched 6.24%, t15 = 0.72, one-sided p = 0.758, two-sided p = 0.485 | trial -1: error 3.12% vs matched 6.73%, t15 = -1.15, one-sided p = 0.134, two-sided p = 0.268
('map', 'after', 'strict') matched b: trial 0: error 9.38% vs matched 3.12%, t15 = 1.10, one-sided p = 0.856, two-sided p = 0.287 | trial -1: error 3.12% vs matched 12.50%, t15 = -2.18, one-sided p = 0.0227, two-sided p = 0.0454
('map', 'after', 'span') matched a: trial 0: error 7.70% vs matched 6.23%, t20 = 0.52, one-sided p = 0.695, two-sided p = 0.61 | trial -1: error 4.05% vs matched 6.61%, t20 = -1.33, one-sided p = 0.0996, two-sided p = 0.199
('map', 'after', 'span') matched b: trial 0: error 7.70% vs matched 4.96%, t20 = 1.12, one-sided p = 0.863, two-sided p = 0.275 | trial -1: error 4.05% vs matched 5.79%, t20 = -0.53, one-sided p = 0.3, two-sided p = 0.6
('map', 'nearest', 'strict') matched a: trial 0: error 10.62% vs matched 6.77%, t15 = 0.82, one-sided p = 0.788, two-sided p = 0.424 | trial -1: error 3.12% vs matched 7.18%, t15 = -1.22, one-sided p = 0.121, two-sided p = 0.243
('map', 'nearest', 'strict') matched b: trial 0: error 10.62% vs matched 6.25%, t15 = 0.54, one-sided p = 0.7, two-sided p = 0.599 | trial -1: error 3.12% vs matched 5.42%, t15 = -0.50, one-sided p = 0.312, two-sided p = 0.624
('map', 'nearest', 'span') matched a: trial 0: error 10.17% vs matched 6.33%, t20 = 1.14, one-sided p = 0.867, two-sided p = 0.267 | trial -1: error 2.38% vs matched 6.52%, t20 = -2.55, one-sided p = 0.00957, two-sided p = 0.0191
('map', 'nearest', 'span') matched b: trial 0: error 10.17% vs matched 6.03%, t20 = 1.81, one-sided p = 0.957, two-sided p = 0.0852 | trial -1: error 2.38% vs matched 7.78%, t20 = -1.52, one-sided p = 0.072, two-sided p = 0.144
('voxel', 'after', 'strict') (c) trial 0: error 2.08% vs matched with Enter at 0 1.47%, t15 = 0.53, one-sided p = 0.697 | error 2.08% vs all labelled Enter trials 0.81%, t15 = 0.68, one-sided p = 0.745
('voxel', 'after', 'span') (c) trial 0: error 1.59% vs matched with Enter at 0 1.42%, t20 = 0.18, one-sided p = 0.57 | error 1.59% vs all labelled Enter trials 0.88%, t20 = 0.49, one-sided p = 0.684
('voxel', 'nearest', 'strict') (c) trial 0: error 2.08% vs matched with Enter at 0 1.25%, t15 = 0.57, one-sided p = 0.712 | error 2.08% vs all labelled Enter trials 0.89%, t15 = 0.63, one-sided p = 0.731
('voxel', 'nearest', 'span') (c) trial 0: error 1.59% vs matched with Enter at 0 1.25%, t20 = 0.29, one-sided p = 0.614 | error 1.59% vs all labelled Enter trials 0.90%, t20 = 0.47, one-sided p = 0.679
('map', 'after', 'strict') (c) trial 0: error 9.38% vs matched with Enter at 0 7.30%, t15 = 0.56, one-sided p = 0.708 | error 9.38% vs all labelled Enter trials 7.97%, t15 = 0.32, one-sided p = 0.625
('map', 'after', 'span') (c) trial 0: error 7.70% vs matched with Enter at 0 7.42%, t20 = 0.11, one-sided p = 0.542 | error 7.70% vs all labelled Enter trials 8.02%, t20 = -0.12, one-sided p = 0.454
('map', 'nearest', 'strict') (c) trial 0: error 10.62% vs matched with Enter at 0 6.64%, t15 = 0.99, one-sided p = 0.831 | error 10.62% vs all labelled Enter trials 7.56%, t15 = 0.71, one-sided p = 0.755
('map', 'nearest', 'span') (c) trial 0: error 10.17% vs matched with Enter at 0 6.44%, t20 = 1.14, one-sided p = 0.865 | error 10.17% vs all labelled Enter trials 7.93%, t20 = 0.69, one-sided p = 0.751
paper: trial 0 t24 = 2.3, p = 0.01; trial -1 t23 = 1.5, p = 0.07
Result: not reproduced, resting on the class-matched reading (c). Every error trial that can be given a state is an Enter trial, where the paper's errors were of all three classes, and correct Enter trials decode far worse than others (0.90% against 11.79% and 14.30%, z by voxel), as that classifier names an Enter state on 5.60% of volumes (data check). Matched to their class, error trials are not lower: 1.59% against 1.42% at position 0 of matched courses (a) with an Enter trial there (t20 = 0.18, one-sided p = 0.57) and against 0.88% on all labelled Enter trials (p = 0.68), and p >= 0.45 in every reading. The readings of the criterion, which do not match class, give different outcomes: with maps z-scored by voxel and matched courses (a), error trials are lower (1.59% against 9.21%, t20 = -4.07, one-sided p = 0.0003; strict 2.08% against 8.93%, p = 0.005); with one matched course per error (b), p = 0.044 to 0.14; with maps z-scored by map the error trial is not lower (p >= 0.69). Reading (c) is adopted because it compares like with like. That the paper's matched courses had the class mix of its error trials is an inference: the Figure 4B legend says only "seven consecutive trials with no behavioral error". With maps z-scored by voxel, reading (c) is at floor (1.59% against 1.42% and 0.88%, below the 6.25% chance level) and could not show a drop. The reading that could is the map z-scoring, where labelled Enter trials decode at 8.02%: there the error trial is not lower under the criterion's readings (p >= 0.69) or under (c) (p >= 0.45). Trial -1 is not lower with maps z-scored by voxel (11.59% against 8.81%, one-sided p = 0.70; paper p = 0.07). The degrees of freedom are 20 (21 participants), against t24 in the paper.
fig5-similarity-predicts-errors and figs5-euclidean-distance¶
Across the eight Exit to Enter transitions, transitions between states whose OFC patterns are more correlated have lower error rates.
The relation between pattern similarity and error rate holds with between-run Euclidean distances in place of correlations.
For each of the eight Enter states, the Exit state before it (its previous category and age, after an age change) and the error rate on the Enter trial: errors whose state is reconstructed (span or strict reading) over those errors plus the correct trials of that state. Enter errors whose current age is unknown cannot be given a transition and are left out. Readings: (a) regression over the eight transition types of the group-mean error rate on the group-mean similarity; (b) over the eight ordinal positions of Figure 5B (NaN correlations, the negative ones, ranked lowest); Figure 5C, each participant's Pearson correlation over the transitions, t-test against 0; and a mixed model of error rate on similarity, random intercept and slope by participant.
ENTER_STATES = [s for s in ofc4.STATES if ofc4.trial_class(s) == "Enter"]
EXIT_OF = {s: s[2:] + s[2] + ("y" if s[3] == "o" else "o") for s in ENTER_STATES}
assert all(ofc4.trial_class(x) == "Exit" for x in EXIT_OF.values())
def transition_table(err_reading, prep, metric):
col = "state_span" if err_reading == "span" else "state"
rows = []
for s in SUBJ:
M = SIM[(prep, metric, s)]
M = (M + M.T) / 2
for en in ENTER_STATES:
ex = EXIT_OF[en]
n_ok = int(((lab.subject == s) & (lab.cond == en)).sum())
n_err = int(((ERR.subject == s) & (ERR[col] == en)).sum())
rows.append({"subject": s, "enter": en, "exit": ex, "sim": M[ofc4.STATES.index(ex), ofc4.STATES.index(en)],
"rate": n_err / (n_ok + n_err), "n_ok": n_ok})
return pd.DataFrame(rows)
def fig5_tests(tt, metric):
g = tt.groupby("enter")[["sim", "rate"]].mean()
ra = stats.linregress(g.sim, g.rate)
ordd = tt.assign(sim_rank=tt.sim.fillna(-np.inf)).sort_values(["subject", "sim_rank"])
ordd["pos"] = ordd.groupby("subject").cumcount()
go = ordd.groupby("pos")[["sim", "rate"]].mean().dropna()
rb = stats.linregress(go.sim, go.rate) if len(go) >= 3 else None
rs_ = []
for s, d in tt.dropna(subset=["sim"]).groupby("subject"):
if len(d) >= 3 and d.rate.std() > 0 and d.sim.std() > 0:
rs_.append(stats.pearsonr(d.sim, d.rate).statistic)
rs_ = np.array(rs_)
alt = "less" if metric == "corr" else "greater"
tc, pc = stats.ttest_1samp(rs_, 0, alternative=alt)
_, pc2 = stats.ttest_1samp(rs_, 0)
# mixed model on similarity standardised over the table (raw distances are in the thousands)
d = tt.dropna(subset=["sim"]).assign(sim=lambda x: (x.sim - x.sim.mean()) / x.sim.std())
try:
m = smf.mixedlm("rate ~ sim", d, groups=d["subject"], re_formula="~sim").fit()
ok = bool(m.converged)
except np.linalg.LinAlgError:
ok = False
if not ok:
m = smf.mixedlm("rate ~ sim", d, groups=d["subject"]).fit()
chi2, _, pm = wald(m, ["sim"])
return ra, rb, rs_, tc, pc, pc2, chi2, pm, m.params["sim"]
FIG5 = {}
for metric in ("corr", "dist"):
for prep in (("raw", "voxel") if metric == "corr" else ("raw", "voxel", "map")):
for err_reading in ("span", "strict"):
tt = transition_table(err_reading, prep, metric)
ra, rb, rs_, tc, pc, pc2, chi2, pm, bm = fig5_tests(tt, metric)
FIG5[(metric, prep, err_reading)] = (ra, rb, rs_, tc, pc)
print(f"{metric}, maps {prep}, errors {err_reading} ({int(tt.rate.gt(0).sum())} participant-transitions with errors, "
f"{tt.sim.isna().sum()} NaN of {len(tt)}):\n"
f" (a) types slope {ra.slope:.4g}, t6 = {ra.slope / ra.stderr:.2f}, p = {ra.pvalue:.3g} | (b) ordinal "
+ (f"slope {rb.slope:.4g}, t = {rb.slope / rb.stderr:.2f}, p = {rb.pvalue:.3g}" if rb else "not computable (too few positions with a correlation)") + "\n"
f" 5C: mean within-participant r {rs_.mean():.3f} (n = {len(rs_)}), t = {tc:.2f}, one-sided p = {pc:.3g}, two-sided p = {pc2:.3g} | "
f"mixed model slope {bm:.4g} per SD, chi2(1) = {chi2:.2f}, p = {pm:.3g}")
# Control for trial counts (variant): errors take trials out of a state's regressor, so a state with
# more errors has fewer correct trials and noisier beta maps. Within participant, the correlation of
# distance with the Enter state's correct-trial count, and the partial correlation of distance and
# error rate with that count held constant.
for prep in ("raw", "voxel", "map"):
tt = transition_table("span", prep, "dist")
r_n, r_p = [], []
for s, d in tt.groupby("subject"):
r_n.append(stats.pearsonr(d.sim, d.n_ok).statistic)
if d.rate.std() > 0:
X = np.column_stack([np.ones(len(d)), d.n_ok])
res = lambda v: v - X @ np.linalg.lstsq(X, v, rcond=None)[0]
r_p.append(stats.pearsonr(res(d.sim.values), res(d.rate.values)).statistic)
r_n, r_p = np.array(r_n), np.array(r_p)
print(f"dist, maps {prep}, span errors: r(distance, correct trials) mean {r_n.mean():.3f}, two-sided p = {stats.ttest_1samp(r_n, 0).pvalue:.3g} | "
f"partial r(distance, error rate | correct trials) mean {r_p.mean():.3f} (n = {len(r_p)}), one-sided p = {stats.ttest_1samp(r_p, 0, alternative='greater').pvalue:.3g}")
print("paper, correlations: t6 = -3.23, p = 0.02; Figure 5B p = 0.04; Figure 5C p < 0.01; mixed model chi2(1) = 3.0, p = 0.08")
print("paper, Euclidean: Figure S5B p = 0.03, Figure S5C p < 0.01")
corr, maps raw, errors span (69 participant-transitions with errors, 37 NaN of 192): (a) types slope -0.6178, t6 = -1.50, p = 0.185 | (b) ordinal slope 0.1605, t = 0.87, p = 0.418 5C: mean within-participant r -0.010 (n = 21), t = -0.09, one-sided p = 0.465, two-sided p = 0.93 | mixed model slope -0.003521 per SD, chi2(1) = 1.37, p = 0.242
corr, maps raw, errors strict (35 participant-transitions with errors, 37 NaN of 192): (a) types slope -0.6656, t6 = -2.09, p = 0.0822 | (b) ordinal slope 0.2607, t = 2.47, p = 0.0485 5C: mean within-participant r 0.041 (n = 16), t = 0.29, one-sided p = 0.614, two-sided p = 0.773 | mixed model slope -0.002172 per SD, chi2(1) = 0.97, p = 0.326
corr, maps voxel, errors span (69 participant-transitions with errors, 123 NaN of 192): (a) types slope -1.44, t6 = -1.20, p = 0.274 | (b) ordinal slope 0.5345, t = 1.71, p = 0.163 5C: mean within-participant r 0.157 (n = 9), t = 0.80, one-sided p = 0.777, two-sided p = 0.447 | mixed model slope -0.001173 per SD, chi2(1) = 0.01, p = 0.915
corr, maps voxel, errors strict (35 participant-transitions with errors, 123 NaN of 192): (a) types slope -1.031, t6 = -0.96, p = 0.375 | (b) ordinal slope 1.3, t = 3.04, p = 0.0385 5C: mean within-participant r 0.053 (n = 6), t = 0.19, one-sided p = 0.572, two-sided p = 0.855 | mixed model slope -0.006854 per SD, chi2(1) = 1.08, p = 0.299
dist, maps raw, errors span (69 participant-transitions with errors, 0 NaN of 192): (a) types slope 4.074e-05, t6 = 1.03, p = 0.343 | (b) ordinal slope 3.626e-05, t = 6.17, p = 0.000834 5C: mean within-participant r 0.260 (n = 21), t = 3.13, one-sided p = 0.00263, two-sided p = 0.00526 | mixed model slope -0.006561 per SD, chi2(1) = 2.06, p = 0.151
dist, maps raw, errors strict (35 participant-transitions with errors, 0 NaN of 192): (a) types slope 3.092e-05, t6 = 0.88, p = 0.411 | (b) ordinal slope 1.013e-05, t = 1.58, p = 0.166 5C: mean within-participant r 0.147 (n = 16), t = 1.31, one-sided p = 0.105, two-sided p = 0.21 | mixed model slope 0.0009559 per SD, chi2(1) = 0.28, p = 0.598
dist, maps voxel, errors span (69 participant-transitions with errors, 0 NaN of 192): (a) types slope 0.002112, t6 = 1.05, p = 0.332 | (b) ordinal slope 0.002532, t = 4.10, p = 0.00638 5C: mean within-participant r 0.313 (n = 21), t = 4.01, one-sided p = 0.000344, two-sided p = 0.000689 | mixed model slope 0.01074 per SD, chi2(1) = 6.78, p = 0.0092
dist, maps voxel, errors strict (35 participant-transitions with errors, 0 NaN of 192): (a) types slope 0.001322, t6 = 0.73, p = 0.493 | (b) ordinal slope 0.0007569, t = 2.10, p = 0.08 5C: mean within-participant r 0.197 (n = 16), t = 2.03, one-sided p = 0.0302, two-sided p = 0.0604 | mixed model slope 0.003717 per SD, chi2(1) = 3.35, p = 0.0672
dist, maps map, errors span (69 participant-transitions with errors, 0 NaN of 192): (a) types slope -0.0001097, t6 = -0.01, p = 0.994 | (b) ordinal slope -0.003176, t = -1.42, p = 0.206 5C: mean within-participant r -0.088 (n = 21), t = -0.82, one-sided p = 0.79, two-sided p = 0.419 | mixed model slope -0.005183 per SD, chi2(1) = 1.04, p = 0.308
dist, maps map, errors strict (35 participant-transitions with errors, 0 NaN of 192): (a) types slope -0.0005489, t6 = -0.05, p = 0.965 | (b) ordinal slope -0.003723, t = -4.04, p = 0.00678 5C: mean within-participant r -0.142 (n = 16), t = -1.20, one-sided p = 0.875, two-sided p = 0.25 | mixed model slope 0.001628 per SD, chi2(1) = 0.70, p = 0.404 dist, maps raw, span errors: r(distance, correct trials) mean -0.212, two-sided p = 0.0133 | partial r(distance, error rate | correct trials) mean 0.072 (n = 21), one-sided p = 0.215
dist, maps voxel, span errors: r(distance, correct trials) mean -0.258, two-sided p = 0.00298 | partial r(distance, error rate | correct trials) mean 0.076 (n = 21), one-sided p = 0.175 dist, maps map, span errors: r(distance, correct trials) mean 0.093, two-sided p = 0.394 | partial r(distance, error rate | correct trials) mean -0.005 (n = 21), one-sided p = 0.523 paper, correlations: t6 = -3.23, p = 0.02; Figure 5B p = 0.04; Figure 5C p < 0.01; mixed model chi2(1) = 3.0, p = 0.08 paper, Euclidean: Figure S5B p = 0.03, Figure S5C p < 0.01
Result: not reproduced. No reading gives a negative regression of error rate on pattern correlation at p < 0.05: over the eight transition types the slope is negative but not significant (p = 0.08 to 0.37), and over ordinal positions it is positive where it reaches p < 0.05 (strict errors, raw maps p = 0.049, voxel z-scoring p = 0.039). The within-participant correlations are not below 0 (mean r -0.010 to 0.157, one-sided p >= 0.46), and the mixed model finds no effect (p >= 0.24; paper p = 0.08). Error rates rest on the reconstructed Enter errors only.
Result: reproduced, resting on the span reading of error states. With Euclidean distances and the span reading, error rate rises with distance over the ordinal positions (raw maps t = 6.17, p = 0.0008; voxel z-scoring t = 4.10, p = 0.0064) and the within-participant correlations are above 0 (mean r 0.260, p = 0.003; 0.313, p = 0.0003). With the strict reading (39 errors in 35 participant-transitions, 16 participants) neither holds at p < 0.05 for raw maps, and with maps z-scored by map the relation is absent or reversed. The regression over transition types is not significant in any reading (p >= 0.33). The paper reports the same result for correlations, which is not found here (fig5-similarity-predicts-errors).
Errors take trials out of a state's regressor, so a state with more errors has fewer correct trials and noisier maps. Distance falls with the Enter state's correct-trial count (mean within-participant r -0.212, p = 0.013 raw; -0.258, p = 0.003 voxel z-scored), and with that count held constant the distance and the error rate are not related (mean partial r 0.072, one-sided p = 0.22; 0.076, p = 0.18). The paper's analysis has the same coupling, as its error trials were also left out of the state regressors, so the outcome follows the criterion; the result may come from trial counts.
Summary¶
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.
| Claim | Outcome |
|---|---|
| behavior-error-rate | blocked: metadata insufficient (errors not countable) |
| supp-errors-by-trial-class | partially reproduced (assumption: one error per block) |
| supp-errors-stable-over-blocks | reproduced (assumption: one error per block) |
| supp-errors-not-by-age-category | reproduced (assumption: span reading) |
| supp-training-errors-higher | blocked: data missing |
| supp-training-rt-unchanged | blocked: data missing |
| supp-rt-stable-over-blocks | reproduced |
| supp-rt-faces-faster | reproduced |
| supp-rt-by-trial-class | reproduced |
| fig2a-16way-decoding | partially reproduced |
| fig2-two-back-state-space | partially reproduced |
| fig2b-hidden-components | partially reproduced |
| fig2b-current-age-absent | reproduced |
| fig2b-two-back-absent | reproduced |
| figs2c-8way-decoding | reproduced (assumption: voxel z-scoring) |
| figs2a-permutation-control | partially reproduced |
| synthetic-data-control | blocked: metadata insufficient |
| fig2c-enter-cluster | reproduced |
| fig2d-enter-nonenter-decoding | not reproduced |
| fig3-searchlight-conjunction | blocked: data missing |
| other-rois-incomplete | blocked: data missing |
| fig4a-decoding-error-correlation | partially reproduced (assumption: voxel z-scoring) |
| fig4b-error-trial-drop | not reproduced |
| fig4b-decline-before-errors | not reproduced |
| fig4b-error-vs-matched | not reproduced (assumption: class-matched comparison) |
| fig5-similarity-predicts-errors | not reproduced |
| figs5-euclidean-distance | reproduced (assumption: span reading) |
| fig5-synthetic-control | blocked: metadata insufficient |
| supp-motor-action-not-in-ofc | blocked: data missing |
How the analysis differs from the paper's:
- 24 of the 27 participants (those of the paper's trial-wise analysis, as the data set says).
- No spatial smoothing of beta maps or volumes: the files hold no voxel coordinates.
- The GLM's nuisance set is the data set's (four drift terms, six motion parameters), without the 26 cardiac regressors of the paper; it uses the data set's design matrices, built with AFNI.
- Beta maps are z-scored per voxel across the state maps of a run (an assumption; the per-map reading is shown throughout).
- Errors are the error blocks of the event files, one per block, a lower bound; the state of an error trial is reconstructed from its neighbours (strict and span readings), and only Enter errors can be reconstructed.
- The logistic readings of the behavioural models use participant fixed effects, as statsmodels has no frequentist logistic mixed model; the linear mixed models use statsmodels' MixedLM in place of lmer.
- The data set's 'err' column puts the trials inside an error block after the error trial in no regressor; both readings of O3 are rebuilt from the error blocks and shown for the 16-way decoding and its components, with the same outcomes.
- Every error trial that can be given a state is an Enter trial, and with maps z-scored by voxel Enter trials decode worse than others on single volumes, so the error trial is compared with Enter trials of the matched courses (an assumption; the criterion's readings are shown).
- The trial-wise mixed models fall back to a random intercept where the random-slope fit did not converge.
What was found about the data set: the error events merge each error with its repetitions and often the trials after it, so the number of errors, time-outs, error RTs and the state of most error trials cannot be recovered; Enter errors are 181 of 228 error blocks, 79% against the 48% the paper's class rates imply; the 'err' column of the design matrices is not the boxcar the descriptor describes; three runs lack a 'first' column; one pair of adjacent labelled trials in Subj-01 run 4 breaks the task rules; and the files hold no voxel coordinates, no training session, no whole-brain data or other ROIs, and no response side or mapping.
print(f"notebook ran in {(time.time() - T0) / 60:.1f} min")
notebook ran in 5.5 min