Richardson et al. (2018) on OpenNeuro ds000228, snapshot 1.1.1¶
Paper: Richardson, Lisandrelli, Riobueno-Naylor & Saxe (2018), "Development of the social brain from age three to twelve years", Nature Communications 9, 1027, doi:10.1038/s41467-018-03399-2.
Dataset: OpenNeuro ds000228, snapshot
1.1.1, DOI 10.18112/openneuro.ds000228.v1.1.1. The paper cites the dataset on OpenfMRI
without a version; the claims file names snapshot 1.1.1, which is also the latest published
snapshot. Its CHANGES file says 1.1.1 only updated the citation to the paper (2023), so the
pinned snapshot holds the same data as 1.1.0. It is pinned because the claims file gives it and
the protocol pins the latest published version when the paper cites none.
Claims: claims/openneuro-ds000228/10.1038-s41467-018-03399-2.yaml, written from the version of
record and its Supplementary Information (Supplementary Tables 1 to 4 read from the PDF). Its
header labels the shared pipeline P1 to P9 and the open details O1 to O9; this notebook uses the
same labels.
What the snapshot holds (listed below from the archive): for all 155 participants the raw BOLD and T1 images and the paper's own preprocessed data (SPM8, normalized, smoothed; ART artifact indicators and five CompCor regressors per run; white-matter and brain masks), the group ROIs of every network the paper uses (ToM, Pain, face and scene parcels, M1, V1), and the behavioural scores in participants.tsv. Every file is read into memory from its URL pinned to the snapshot.
| Claim | Tested here |
|---|---|
| participants, supp-table1-behavior-by-age, beh-tom-age, beh-tom-dccs, beh-fb-groups, beh-fb-tom, beh-fb-age, beh-fb-dccs | tested (participants.tsv) |
| motion-child-adult, motion-not-age-tom, motion-not-fb | tested (deposited ART indicators) |
| irc-adult, irc-age-within, irc-age-across, irc-adult-vs-child, irc-3yo-distinct, irc-3yo-not-anticorrelated, irc-3yo-homologous, irc-3yo-nonbilateral, irc-tom-behavior, irc-tom-behavior-age, irc-fb-within, irc-fb-within-dccs, irc-fb-across | tested (deposited preprocessed BOLD and ROIs) |
| irc-control-networks, irc-williams, supp-table4-other-networks, face-rc, face-overlap | tested; the face and scene parcels the claims file lists as outside input are deposited in the snapshot |
| rc-events, rc-overlap, rc-child-adult-timecourse, rc-age-events, rc-t04-tom, rc-fb-no-difference, rc-3yo-events | tested |
| maturity-across, maturity-within, maturity-tom-behavior, maturity-fb | tested |
| replication-sample | blocked: the independent sample of 20 adults is not deposited anywhere |
| hand-coded-events | blocked: the hand coding of Jacoby et al. 2016 is not in the snapshot or the paper, and the lab's script that holds it is not pinned |
In [1]:
import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import matplotlib.pyplot as plt
import social_brain as sb
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
print("dataset", sb.DATASET_ID, "snapshot", sb.SNAPSHOT, "doi", sb.SNAPSHOT_DOI)
files = sb.file_table()
print(len(files), "files,", round(files["size"].sum() / 1e9, 2), "GB")
kinds = files.index.str.replace(r"sub-pixar\d+", "sub-*", regex=True)
print(files.groupby(kinds)["size"].agg(["count", "sum"]).rename(columns={"sum": "bytes"}).to_string())
dataset ds000228 snapshot 1.1.1 doi 10.18112/openneuro.ds000228.v1.1.1
2060 files, 22.22 GB
count bytes
filename
.datalad/.gitattributes 1 132
.datalad/config 1 63
.gitattributes 1 284
CHANGES 1 453
README 1 1593
dataset_description.json 1 1412
derivatives/README.txt 1 5293
derivatives/ROIs/AMCC_9mm_sphere_xyz.mat 1 34256
derivatives/ROIs/DMPFC_9mm_sphere_xyz.mat 1 31696
derivatives/ROIs/LInsula_9mm_sphere_xyz.mat 1 33536
derivatives/ROIs/LMFG_9mm_sphere_xyz.mat 1 34816
derivatives/ROIs/LS2_9mm_sphere_xyz.mat 1 35856
derivatives/ROIs/LTPJ_9mm_sphere_xyz.mat 1 43776
derivatives/ROIs/MMPFC_9mm_sphere_xyz.mat 1 36336
derivatives/ROIs/PC_9mm_sphere_xyz.mat 1 44896
derivatives/ROIs/RInsula_9mm_sphere_xyz.mat 1 39056
derivatives/ROIs/RMFG_9mm_sphere_xyz.mat 1 25696
derivatives/ROIs/RS2_9mm_sphere_xyz.mat 1 43776
derivatives/ROIs/RTPJ_9mm_sphere_xyz.mat 1 44416
derivatives/ROIs/VMPFC_9mm_sphere_xyz.mat 1 30176
derivatives/ROIs/lFFA.mat 1 5096
derivatives/ROIs/lOFA.mat 1 6976
derivatives/ROIs/lPPA.mat 1 5197
derivatives/ROIs/lRSC.mat 1 5419
derivatives/ROIs/lSTS.mat 1 5385
derivatives/ROIs/lTOS.mat 1 5104
derivatives/ROIs/left_CS.mat 1 5013
derivatives/ROIs/left_primotor.mat 1 4976
derivatives/ROIs/rFFA.mat 1 5462
derivatives/ROIs/rOFA.mat 1 5277
derivatives/ROIs/rPPA.mat 1 5096
derivatives/ROIs/rRSC.mat 1 5932
derivatives/ROIs/rSTS_masked.mat 1 6483
derivatives/ROIs/rTOS.mat 1 7936
derivatives/ROIs/right_CS.mat 1 4993
derivatives/ROIs/right_primotor.mat 1 4989
derivatives/mriqc/T1w.csv 1 205539
derivatives/mriqc/T1w_group.html 1 609338
derivatives/mriqc/bold.csv 1 125387
derivatives/mriqc/bold_group.html 1 466253
derivatives/mriqc/derivatives/sub-*_T1w.json 155 21421711
derivatives/mriqc/derivatives/sub-*_task-pixar_bold.json 155 41145626
derivatives/mriqc/reports/T1w_group.html 1 609338
derivatives/mriqc/reports/bold_group.html 1 466253
derivatives/mriqc/reports/sub-*_T1w.html 155 2335064149
derivatives/mriqc/reports/sub-*_task-pixar_bold.html 155 351818802
derivatives/preprocessed_data/sub-*/sub-*_CompCor_mask.nii.gz 155 1104530
derivatives/preprocessed_data/sub-*/sub-*_analysis_mask.nii.gz 155 2588959
derivatives/preprocessed_data/sub-*/sub-*_normed_anat.nii.gz 155 74963282
derivatives/preprocessed_data/sub-*/sub-*_task-pixar_run-001_ART_and_CompCor_nuisance_regressors.mat 155 1177018
derivatives/preprocessed_data/sub-*/sub-*_task-pixar_run-001_swrf_bold.nii.gz 155 14049601947
participants.json 1 4971
participants.tsv 1 16041
sub-*/anat/sub-*_T1w.json 155 22924224
sub-*/anat/sub-*_T1w.nii.gz 155 1028398574
sub-*/func/sub-*_task-pixar_bold.json 155 51641821
sub-*/func/sub-*_task-pixar_bold.nii.gz 155 4230197234
task-pixar_bold.json 1 887
In [2]:
# The files read whole, with their checksums (OpenNeuro `id`); the per-participant files read are
# the swrf BOLD image and the nuisance-regressor file of each participant, listed after extraction.
for f in ["dataset_description.json", "README", "CHANGES", "participants.tsv", "participants.json",
"task-pixar_bold.json", "derivatives/README.txt"]:
print(f, files.loc[f, "id"])
for r in sb.ALL_ROIS:
f = [k for k in files.index if k.startswith("derivatives/ROIs/" + r)]
print(f[0], files.loc[f[0], "id"])
print()
print(sb.text("CHANGES"))
print(sb.text("README"))
dataset_description.json 940a08748c32d0bc604d9b6174a56b675b27a8b8 README 6d37c17444778fd32d906a840b2eb2ec0c36bf70 CHANGES b83ad35939cb9083805d1451e2e6fe5affe0835b participants.tsv 4c6a181900081abf6f5bf6cb49b34e3a18e5798e participants.json 0a010b8ad59370db587afec9d0df4ad9a7bb999f task-pixar_bold.json a7997d30dece9e24eed6e0f4e1172e29cb1b1c3e derivatives/README.txt MD5E-s5293--955f1d5b9a262d6402bde466879e6fc4.txt derivatives/ROIs/RTPJ_9mm_sphere_xyz.mat MD5E-s44416--2c059836fe8ab19783f679d2d542a600.mat derivatives/ROIs/LTPJ_9mm_sphere_xyz.mat MD5E-s43776--a913be491eb40c1119160ff2c1b76a84.mat derivatives/ROIs/PC_9mm_sphere_xyz.mat MD5E-s44896--bc64e5c7bc401f2568d1b9eda5ee24dc.mat derivatives/ROIs/DMPFC_9mm_sphere_xyz.mat MD5E-s31696--92846d239f831eec8b83a8190f4c0cc7.mat derivatives/ROIs/MMPFC_9mm_sphere_xyz.mat MD5E-s36336--e13ab79f6f0b65e0019a9de1985351df.mat derivatives/ROIs/VMPFC_9mm_sphere_xyz.mat MD5E-s30176--ad46c858f495ebceb04cee51261c6c6a.mat derivatives/ROIs/RS2_9mm_sphere_xyz.mat MD5E-s43776--b703459602859b3612bd8e70e2bc6bc0.mat derivatives/ROIs/LS2_9mm_sphere_xyz.mat MD5E-s35856--fa96c0ee2c6526f9925f1eb67ddbe7b4.mat derivatives/ROIs/RInsula_9mm_sphere_xyz.mat MD5E-s39056--c975ac2e0b7e07e4e3f426b8ff62190d.mat derivatives/ROIs/LInsula_9mm_sphere_xyz.mat MD5E-s33536--f7c1270d9e9d4bb532f11f003df92890.mat derivatives/ROIs/RMFG_9mm_sphere_xyz.mat MD5E-s25696--12cc5c49fc71cee03911f64dad7307a7.mat derivatives/ROIs/LMFG_9mm_sphere_xyz.mat MD5E-s34816--b90b4b03a79c5db3bd777dc0ab543298.mat derivatives/ROIs/AMCC_9mm_sphere_xyz.mat MD5E-s34256--aa44694e164ba1930b6426e414c80242.mat derivatives/ROIs/rSTS_masked.mat MD5E-s6483--99d06d9447cd1cabe654253914921d65.mat derivatives/ROIs/lSTS.mat MD5E-s5385--fb2e11d3b3f183cd278de13274c4d5ba.mat derivatives/ROIs/rOFA.mat MD5E-s5277--3dba9ef02b05bbf872f77e553cc380d9.mat derivatives/ROIs/lOFA.mat MD5E-s6976--589e9d9354d3388633e0336143bfd44c.mat derivatives/ROIs/rFFA.mat MD5E-s5462--aff09920c32deab5670b68e0488c887a.mat derivatives/ROIs/lFFA.mat MD5E-s5096--41c1ab9b27c2984d02ff80075976cd96.mat derivatives/ROIs/rRSC.mat MD5E-s5932--99661673fdda7b849e0b54df127005b9.mat derivatives/ROIs/lRSC.mat MD5E-s5419--76bca374f91d72dbb916654f5efdc5ac.mat derivatives/ROIs/rTOS.mat MD5E-s7936--3df2beee3f5579daddb8af5cef097c00.mat derivatives/ROIs/lTOS.mat MD5E-s5104--a9b0b1a54edda5cf09bc2d9f979362a1.mat derivatives/ROIs/rPPA.mat MD5E-s5096--0e6b185df5cfcfd28c4c1ea9406362d0.mat derivatives/ROIs/lPPA.mat MD5E-s5197--abc1885231e976c756a1dd71cc97ed1b.mat derivatives/ROIs/right_primotor.mat MD5E-s4989--5a8c4f59d19a5fd604a6b7ee1efc9b78.mat derivatives/ROIs/left_primotor.mat MD5E-s4976--d0d7a280af97002e990558e326bc519d.mat derivatives/ROIs/right_CS.mat MD5E-s4993--8db7a46c296674dcf1449884466462cf.mat derivatives/ROIs/left_CS.mat MD5E-s5013--60aed22c52f48a6ef07ff07e7b499011.mat
1.1.1 2023-09-27 - Updated citation for paper to acknowledge. 1.1.0 2019-05-15 - Renamed preprocessed data folder to make clear that data were not preprocessed through fmriprep, but rather were preprocessed in SPM8. Clarified timing of reverse correlation events, and provided timings for modeling conditions based on these events, accounting for the hemodynamic lag. Updated paper reference/acknowledgements. 1.0.0 2019-04-03 - Initial snapshot Experiment Details Participants watched Disney Pixar’s “Partly Cloudy” while lying in the scanner. There was no task; participants were simply instructed to lie still and watch the movie. The movie began after 10s of rest (black screen; TRs 1-5). The first 10s of the movie are the opening credits (disney castle, pixar logo; TRs 6-10). IPS = 168 TRs TR = 2s Experiment length: 5.6 minutes Pixar Movie Reverse Correlation Events Events defined by conducting reverse correlation analysis in two separate adult samples, using the average response in ToM brain regions (ToM events) and in the pain matrix (Pain events). Events listed are those that replicated across the two samples. Onsets and Durations are noted in TRs (1 TR = 2s); scanner trigger = TR 1. Event types: Theory of Mind (ToM), Physical Sensation/Pain (Pain) TRs identified by RC analysis, as reported in Richardson et al. 2018 ToM Event Onsets; Durations 46; 2 52; 3 63; 2 91; 8 122; 3 129; 4 153; 3 Pain Event Onsets; Durations 38; 2 49; 1 56; 2 71; 5 100; 2 108; 6 117; 3 134; 3 159; 2 TIMING OF EVENTS FOR MODELING (taking into account hemodynamic lag (4s), scanner trigger = 0): All timings assume 10s from trigger until movie begins to play. In seconds: Mental Event Onsets; Durations 86; 4 98; 6 120; 4 176; 16 238; 6 252; 8 300; 6 Pain Event Onsets; Durations 70; 4 92; 2 106; 4 136; 10 194; 4 210; 12 228; 6 262; 6 312; 4 In TRs when TR=2: Mental Event Onsets; Durations 43; 2 49; 3 60; 2 88; 8 119; 3 126; 4 150; 3 Pain Event Onsets; Durations 35; 2 46; 1 53; 2 68; 5 97; 2 105; 6 114; 3 131; 3 156; 2
In [3]:
print(sb.text("derivatives/README.txt"))
print(sb.text("participants.json"))
print(sb.text("task-pixar_bold.json"))
Derivatives fMRI Prep sub-[]_task-pixar_run-001_swrf_bold.nii.gz BOLD images collected during movie viewing. FMRI data were analyzed using SPM8 (http://www.fil.ion.ucl.ac.uk/spm) and custom software written in Matlab. Functional images were registered to the first image of the run; that image was registered to each participant’s anatomical image, and each participant’s anatomical image was normalized to the Montreal Neurological Institute (MNI) template. Registration of each individual’s brain to the MNI template was visually inspected, including checking the match of the cortical envelope and internal features like the AC-PC and major sulci. All data were smoothed using a Gaussian filter (5mm kernel). As participants were initially recruited for different studies, there are small differences in voxel size and slice gaps across participants (3.13 mm isotropic with no slice gap (n=5 adults, n=3 7yos, n=20 8-12yo); 3.13 mm isotropic with 10% slice gap (n=28 adults), 3 mm isotropic with 20% slice gap (n=1 3yo, n=3 4yo, n=2 7yo, n=1 9yo); 3 mm isotropic with 10% slice gap (all remaining participants)); all functional data were subsequently upsampled in normalized space to 2mm isotropic voxels. IPS = 168 TRs; 1 TR = 2s TR 1-5: Rest (black screen) TR 6-10: Pixar and Disney logos TR 11-end: Movie sub-[]_normed_anat.nii.gz Individual subject anatomical image, normalized to the Montreal Neurological Institute (MNI) template. sub-[]_analysis_mask.nii.gz Individual subject skull-stripped whole-brain mask. sub-[]_CompCor_mask.nii.gz Individual subject white matter mask, normalized in MNI space, and eroded by two voxels in every direction (in order to avoid partial voluming with grey matter). This white matter mask was used for creation of CompCor regressors. sub-[]_task-pixar_run-001_ART_and_CompCor_nuisance_regressors.mat Individual subject artifact timepoint regressors and CompCor regressors. Artifact timepoint regressors were created using the Artifact Detection Tool (ART; https://www.nitrc.org/projects/artifact_detect/) as timepoints for which there was 1) more than 2mm composite motion relative to the previous timepoint or 2) a fluctuation in global signal that exceeded a threshold of three standard deviations from the mean global signal. Compcor regressors are 5 principle component analysis (PCA)-based noise regressors generated using CompCor within individual subject eroded white matter masks. CompCor regressors were defined using scrubbed data (e.g. artifact timepoints were identified and interpolated over prior to running CompCor). ROIs (Regions of Interest) Center coordinates are in MNI space; size is in voxels. ToM Network RTPJ_9mm_sphere_xyz.mat (center: 48,-60,30; size: 376) LTPJ_9mm_sphere_xyz.mat (center: -48,-62,30; size: 368) PC_9mm_sphere_xyz.mat (center: 0,-54,34; size: 382) DMPFC_9mm_sphere_xyz.mat (center: -6,54,36; size: 217) MMPFC_9mm_sphere_xyz.mat (center: -4,58,16; size: 275) VMPFC_9mm_sphere_xyz.mat (center: -4,56,-16; size: 198) Pain Matrix RS2_9mm_sphere_xyz.mat (center: 60,-28,38; size: 368) LS2_9mm_sphere_xyz.mat (center: -62,-32,34; size: 269) RInsula_9mm_sphere_xyz.mat (center: 42,6,-6; size: 309) LInsula_9mm_sphere_xyz.mat (center: -42,-2,-4; size: 240) RMFG_9mm_sphere_xyz.mat (center: 50,42,12; size: 142) LMFG_9mm_sphere_xyz.mat (center: -46,36,14; size: 256) AMCC_9mm_sphere_xyz.mat (center: 0,2,42; size: 249) Face Network rSTS_masked.mat lSTS.mat rOFA.mat lOFA.mat rFFA.mat lFFA.mat Place Network rRSC.mat lRSC.mat rTOS.mat lTOS.mat rPPA.mat lPPA.mat Primary motor cortex right_primotor.mat left_primotor.mat Primary visual cortex right_CS.mat left_CS.mat ToM and pain matrix group ROIs were created in an independent group of adults (n=20), scanned by Evelina Fedorenko and colleagues. These data were preprocessed and analyzed with procedures identical to those used for participants in the current study. Reverse correlation analyses were conducted in this separate group of adults, using 10mm group ROIs surrounding peaks reported in previous publications (ToM regions: Dufour et al., 2013; Pain matrix: Bruneau et al., 2015). Seven ToM and nine pain events were identified (ToM: 60s total, M(SD) length: 8.6(4.6)s, Pain: 66s total, M(SD) length: 7.3(4.4)s). We subsequently used a general-linear model to analyze BOLD activity of these participants as a function of condition, using these events. Second-level random effects analyses were used to examine the group-level response to Mental > Pain and Pain> Mental (p<.001, k=10, uncorrected). We then drew 9mm spheres surrounding the peak activation in each region. Face and scene parcels are publicly available (http://web.mit.edu/bcs/nklab/GSS.shtml), and described in: Julian, J. B., Fedorenko, E., Webster, J. & Kanwisher, N. An algorithmic method for functionally defining regions of interest in the ventral visual pathway. NeuroImage 60, 2357–2364 (2012). Primary motor and visual cortex ROIs are 10mm spheres drawn around peak coordinates generated with Neurosynth (http://neurosynth.org/; M1: search term: “primary motor,” coordinates: [38,-24,58], [-38,-20,58]; V1: search term: “primary visual,” coordinates: [-10 -86 2], [10 -86 2].
{
"Participant Demographics, Behavioral Data, and Scan Info":{
"Theory of Mind (ToM) Measures Metadata": {
"Description": "ToM Measures are summary scores describing performance on a ToM behavioral battery. The ToM battery is publicly available at the link below.",
"TermURL":"osf.io/g5zpv"
},
"ToM Booklet-Matched":{
"Description":"Proportion of questions answered correctly on n=24 matched items in behavioral battery. These items include the six explicit false-belief questions used to group children based on false-belief task performance."
},
"ToM Booklet-Matched-NOFB":{
"Description":"Proportion of questions answered correctly on n=18 matched items in behavioral battery. These items do not include the six explicit false-belief questions used to group children based on false-belief task performance."
},
"FB_Composite":{
"Description":"Number of explicit false-belief questions answered correctly, out of six questions. This number was used to group children as false belief task passers (5-6 out of 6 correct), inconsistent performers (3-4 out of 6 correct), and failers (0-2 out of 6 correct)."
},
"FB_Group":{
"Description":"Explicit false-belief task groups. Groups were assigned based on the number of explicit FB questions answered correctly (see FB_Composite)",
"Levels":{
"Pass":"5-6 out of 6 FB questions correct",
"Inc":"3-4 out of 6 FB questions correct",
"Fail":"0-2 out of 6 FB questions correct"
}
},
"Nonverbal IQ Measures Metadata":{
"Description":"Nonverbal IQ was measured with the Wechsler Preschool and Primary Scale of Intelligence (WPPSI) Block Design task for all children younger than age five years (n=65). Children five years and older completed the Matrices task of the Kaufman Brief Intelligence Test (KBIT-2)"
},
"WPPSI BD raw":{
"Description": "Raw score on WPPSI block design task"
},
"WPPSI BD scaled":{
"Description": "Scaled score on WPPSI block design task"
},
"KBIT_raw":{
"Description": "Raw score on KBIT nonverbal matrices task"
},
"KBIT_standard":{
"Description": "Standardized score on KBIT nonverbal matrices task"
}
},
"Dimensional Change Card Sort Metadata":{
"Description":"Dimensional Change Card Sort (DCCS) is a task that measures the response inhibition component of executive functions, designd for young children (Zelazo, 2006). Children ages 3-5 years old completed this task. The task is comprised of two versions: the standard version (pre-switch task and post-switch task), where children have to sort cards according to one rule (pre-switch), and then according to a new rule (post-switch), and the border version, where children must sort cards according to two rules, depending on whether the card has a border. The shared data is the summary score (0-3), described below.",
"Levels":{
"0": "Child failed the pre-switch phase of the task (fewer than 6/6 pre-switch trials correct",
"1": "Child passed the pre-switch phase of the task (6/6 correct) but fail the post-switch phase (fewer than 5/6 post-switch trials correct",
"2": "Child passed the pre- and post-switch phase of the task (6/6 pre-switch trials correct, 5-6/6 post-switch trials correct), but fail the border version of the task (fewer than 9/12 border trials correct",
"3": "Child passed the pre- and post-switch phase of the standard task, and passed the border version of the task"
}
},
"Scanlog Metadata":{
"Description":"Whole-brain structural and functional MRI data were acquired on one of two 3-Tesla Siemens Tim Trio scanners located at the Athinoula A. Martinos Imaging Center at MIT, using one of three 32-channel head coils (see Coil column). As participants were initially recruited for different studies, there are small differences in voxel size and slice gaps across participants (voxel size: 3 mm vs. 3.13 mm isotropic; slice gap: 0, 10% and 20%)."
},
"Scanlog: Scanner":{
"Levels": {
"3T1": "Scanner 1",
"3T2": "Scanner 2"
}
},
"Scanlog: Coil":{
"Description":"Children under age five years used one of two custom 32-channel phased-array head coils made for younger (n=3, M(SD)=3.91(.42) years) or older (n=28, M(SD)=4.07(.42) years) children (Keil et al., 2011); all other participants used the standard Siemens 32-channel head coil.",
"Levels":{
"5-6yo 32ch":"Custom 32ch coil made for young children",
"7-8yo 32ch":"Custom 32ch coil made for older children",
"32ch adult": "Standard Siemens 32-channel head coil"
}
},
"Scanlog: Voxel Size":{
"Levels":{
"3mm iso":"3mm isotropic",
"3.13 mm iso": "3.13 mm isotropic"
}
},
"Scanlog: Slice Gap":{
"Levels":{
"None":"No slice gap",
".1":"10% slice gap",
".2":"20% slice gap"
}
}
}
{
"EchoTime": 0.03,
"RepetitionTime": 2,
"MagneticFieldStrength": 3,
"dcmmeta_slice_dim": 2,
"PhaseEncodingDirection": "j-",
"EffectiveEchoSpacing": 0.0005,
"dcmmeta_reorient_transform": [
[
0.0,
-1.0,
0.0,
63.0
],
[
1.0,
0.0,
0.0,
0.0
],
[
0.0,
0.0,
1.0,
0.0
],
[
0.0,
0.0,
0.0,
1.0
]
],
"ManufacturersModelName": "TrioTim",
"dcmmeta_version": 0.6,
"TaskName": "pixar",
"dcmmeta_shape": [
64,
64,
32,
168
],
"Manufacturer": "Siemens",
"CogAtlasID": "http://www.cognitiveatlas.org/task/id/trm_4c898da401420"
}
Loading¶
Each participant's preprocessed image (*_swrf_bold.nii.gz, about 90 MB) is read into memory,
and the voxels of the 29 deposited ROIs are kept; the nuisance file gives the artifact
indicators and the CompCor regressors. This takes about two minutes.
In [4]:
p = sb.participants()
subs = list(p.participant_id)
D = sb.extract_all(subs, n_threads=8)
p["n_art"] = [len(sb.artifact_index(D[s]["artifacts"])) for s in subs]
p = p.set_index("participant_id")
TOM, NOFB, DCCS = "ToM Booklet-Matched", "ToM Booklet-Matched-NOFB", "DCCS Summary"
children = p[p.Child_Adult == "child"].copy()
adults = p[p.Child_Adult == "adult"].copy()
young = children[children.Age < 6].copy()
print(len(subs), "participants loaded")
print(pd.Series({s: files.loc[sb.derivative(s, "task-pixar_run-001_swrf_bold.nii.gz"), "id"]
for s in subs[:3]}).to_string(), "... (checksums of all 155 in `files`)")
155 participants loaded sub-pixar001 MD5E-s91565296--35d849c8a633a2ca31555c7275fcc2... sub-pixar002 MD5E-s102442997--a18fbb4967e6528ea8a1c2352ee78... sub-pixar003 MD5E-s103690359--a4c0b27e85df7ea44c57a90931297... ... (checksums of all 155 in `files`)
Data check¶
In [5]:
print("participants.tsv columns:", list(p.reset_index().columns))
print(p.reset_index().describe(include="all").T.to_string())
print()
print("children", len(children), "adults", len(adults))
print("children: female", (children.Gender == "F").sum(), "mean age %.2f (s.d. %.2f)" % (children.Age.mean(), children.Age.std()),
"range %.2f to %.2f" % (children.Age.min(), children.Age.max()))
print("adults: female", (adults.Gender == "F").sum(), "mean age %.2f (s.d. %.2f)" % (adults.Age.mean(), adults.Age.std()),
"range %.0f to %.0f" % (adults.Age.min(), adults.Age.max()))
print(p.groupby("AgeGroup").Age.agg(["count", "min", "max", "mean", "std"]).round(2).to_string())
print("handedness:", p.groupby("Child_Adult").Handedness.value_counts().to_dict())
print("under 6:", len(young), "FB groups", young.FB_Group.value_counts().to_dict(),
"; FB_Group in children aged 6 and over:", children[children.Age >= 6].FB_Group.value_counts().to_dict())
print("DCCS present:", children[DCCS].notna().sum(), "children, all aged under 6:", bool((children.loc[children[DCCS].notna(), "Age"] < 6).all()))
print("FB_Composite against FB_Group (cut-offs of the Methods):")
print(pd.crosstab(young.FB_Composite, young.FB_Group).to_string())
print("ToM score range %.2f to %.2f; no-FB score range %.2f to %.2f" % (children[TOM].min(), children[TOM].max(), children[NOFB].min(), children[NOFB].max()))
print()
print("voxel size by age group (Methods: 3.13 mm for 5 + 28 adults, 3 7yo, 20 8-12yo):")
print(pd.crosstab(p.AgeGroup, [p["Scanlog: Voxel slize"], p["Scanlog: Slice Gap"].fillna(0)]).to_string())
print("coil:", pd.crosstab(p.AgeGroup, p["Scanlog: Coil"]).to_dict())
participants.tsv columns: ['participant_id', 'Age', 'AgeGroup', 'Child_Adult', 'Gender', 'Handedness', 'ToM Booklet-Matched', 'ToM Booklet-Matched-NOFB', 'FB_Composite', 'FB_Group', 'WPPSI BD raw', 'WPPSI BD scaled', 'KBIT_raw', 'KBIT_standard', 'DCCS Summary', 'Scanlog: Scanner', 'Scanlog: Coil', 'Scanlog: Voxel slize', 'Scanlog: Slice Gap', 'n_art']
count unique top freq mean std min 25% 50% 75% max
participant_id 155 155 sub-pixar001 1 NaN NaN NaN NaN NaN NaN NaN
Age 155.0 NaN NaN NaN 10.555189 8.071957 3.518138 5.3 7.68 10.975 39.0
AgeGroup 155 6 5yo 34 NaN NaN NaN NaN NaN NaN NaN
Child_Adult 155 2 child 122 NaN NaN NaN NaN NaN NaN NaN
Gender 155 2 F 84 NaN NaN NaN NaN NaN NaN NaN
Handedness 155 3 R 142 NaN NaN NaN NaN NaN NaN NaN
ToM Booklet-Matched 122.0 NaN NaN NaN 0.775266 0.199293 0.24 0.65 0.818182 0.919167 1.0
ToM Booklet-Matched-NOFB 122.0 NaN NaN NaN 0.771929 0.189931 0.25 0.631579 0.8125 0.935268 1.0
FB_Composite 122.0 NaN NaN NaN 4.639344 1.676604 0.0 4.0 5.0 6.0 6.0
FB_Group 122 3 pass 84 NaN NaN NaN NaN NaN NaN NaN
WPPSI BD raw 30.0 NaN NaN NaN 16.166667 3.63966 8.0 15.0 16.0 18.0 23.0
WPPSI BD scaled 30.0 NaN NaN NaN 10.033333 2.953471 3.0 8.0 10.0 11.75 16.0
KBIT_raw 91.0 NaN NaN NaN 28.32967 8.463063 14.0 20.0 30.0 35.5 44.0
KBIT_standard 91.0 NaN NaN NaN 116.10989 14.056592 79.0 108.0 116.0 123.5 149.0
DCCS Summary 64.0 NaN NaN NaN 2.171875 0.679745 0.0 2.0 2.0 3.0 3.0
Scanlog: Scanner 155 2 3T1 101 NaN NaN NaN NaN NaN NaN NaN
Scanlog: Coil 155 3 32ch adult 124 NaN NaN NaN NaN NaN NaN NaN
Scanlog: Voxel slize 155 2 3mm iso 99 NaN NaN NaN NaN NaN NaN NaN
Scanlog: Slice Gap 127.0 NaN NaN NaN 0.105512 0.022911 0.1 0.1 0.1 0.1 0.2
n_art 155.0 NaN NaN NaN 8.870968 10.076145 0.0 1.0 5.0 12.0 48.0
children 122 adults 33
children: female 64 mean age 6.71 (s.d. 2.33) range 3.52 to 12.30
adults: female 20 mean age 24.77 (s.d. 5.31) range 18 to 39
count min max mean std
AgeGroup
3yo 17 3.52 3.99 3.75 0.18
4yo 14 4.06 4.86 4.43 0.29
5yo 34 5.01 5.99 5.51 0.29
7yo 23 7.00 7.96 7.54 0.37
8-12yo 34 8.02 12.30 9.77 1.18
Adult 33 18.00 39.00 24.77 5.31
handedness: {('adult', 'R'): 32, ('adult', 'L'): 1, ('child', 'R'): 110, ('child', 'L'): 9, ('child', 'Ambi'): 3}
under 6: 65 FB groups {'pass': 30, 'inc': 20, 'fail': 15} ; FB_Group in children aged 6 and over: {'pass': 54, 'inc': 3}
DCCS present: 64 children, all aged under 6: True
FB_Composite against FB_Group (cut-offs of the Methods):
FB_Group fail inc pass
FB_Composite
0.0 3 0 0
1.0 9 0 0
2.0 3 0 0
3.0 0 12 0
4.0 0 8 0
5.0 0 0 18
6.0 0 0 12
ToM score range 0.24 to 1.00; no-FB score range 0.25 to 1.00
voxel size by age group (Methods: 3.13 mm for 5 + 28 adults, 3 7yo, 20 8-12yo):
Scanlog: Voxel slize 3.13 mm iso 3mm iso
Scanlog: Slice Gap 0.0 0.1 0.1 0.2
AgeGroup
3yo 0 0 16 1
4yo 0 0 11 3
5yo 0 0 34 0
7yo 2 0 19 2
8-12yo 21 0 12 1
Adult 5 28 0 0
coil: {'32ch adult': {'3yo': 0, '4yo': 0, '5yo': 34, '7yo': 23, '8-12yo': 34, 'Adult': 33}, '5-6yo 32ch': {'3yo': 2, '4yo': 1, '5yo': 0, '7yo': 0, '8-12yo': 0, 'Adult': 0}, '7-8yo 32ch': {'3yo': 15, '4yo': 13, '5yo': 0, '7yo': 0, '8-12yo': 0, 'Adult': 0}}
In [6]:
shapes = {tuple(D[s]["shape"]) for s in subs}
zooms = {tuple(float(z) for z in D[s]["zooms"][:3]) for s in subs}
print("preprocessed image shapes:", shapes, "voxel sizes (mm):", zooms, "dtype:", {D[s]["dtype"] for s in subs})
print("affines identical:", all(np.allclose(D[s]["affine"], D[subs[0]]["affine"]) for s in subs))
print("168 volumes in every image (TR 2 s from task-pixar_bold.json):", all(D[s]["shape"][3] == sb.N_TR for s in subs))
print("artifact indicator rows each flag one timepoint:", all((D[s]["artifacts"].sum(axis=1) == 1).all() for s in subs))
print("CompCor rows per participant:", {D[s]["compcor"].shape for s in subs})
print("artifact timepoints per participant: max %d of 168 (%.2f; Methods drop at one third)" % (p.n_art.max(), p.n_art.max() / 168))
print(p.groupby("AgeGroup").n_art.agg(["mean", "std", "min", "max"]).round(2).to_string())
roi_size = pd.Series({r: D[subs[0]]["rois"][r]["n_vox"] for r in sb.ALL_ROIS})
print("ROI sizes (voxels):", roi_size.to_dict())
print("Supplementary Table 2 sizes: RTPJ 376 LTPJ 368 PC 382 DMPFC 217 MMPFC 275 VMPFC 198; RS2 368 LS2 269 RInsula 309 LInsula 240 RMFG 142 LMFG 256 AMCC 249; RFFA 1019 LFFA 531")
cover = pd.DataFrame({s: {r: D[s]["rois"][r]["n_with_data"] for r in sb.ALL_ROIS} for s in subs}).T
print("participants with any ROI voxel without data (zero or non-finite at some timepoint), per ROI:")
print((cover < roi_size).sum().to_string())
print("fewest voxels with data in a face ROI:", int(cover[sb.FACE_ROIS].min().min()),
"(Supplementary Fig. 4: every participant had at least 100 voxels in every face ROI)")
face_incomplete = (cover[sb.FACE_ROIS] < roi_size[sb.FACE_ROIS]).any(axis=1)
print("participants with incomplete face-ROI coverage:", int(face_incomplete.sum()),
p.loc[face_incomplete[face_incomplete].index, "AgeGroup"].value_counts().to_dict(),
"(Supplementary Fig. 4: 17, of them 2 adults, 12 8-12yo, one each 7yo, 4yo, 3yo)")
print("NIfTI scl_slope, scl_inter:", {(str(D[s]["scl_slope"]), str(D[s]["scl_inter"])) for s in subs})
rng_img = pd.DataFrame({s: D[s]["range"] for s in subs}, index=["min", "max"]).T
print("image value range: minimum %.1f to %.1f, maximum %.1f to %.1f over participants; %s %s" % (
rng_img["min"].min(), rng_img["min"].max(), rng_img["max"].min(), rng_img["max"].max(), subs[0], D[subs[0]]["range"]))
# ROI coordinates: centroid and largest distance from the stated center, to check that the masks
# are in MNI mm (derivatives/README.txt and Supplementary Table 2 centers; M1 and V1 from Methods).
STATED = {"RTPJ": (48, -60, 30), "LTPJ": (-48, -62, 30), "PC": (0, -54, 34), "DMPFC": (-6, 54, 36), "MMPFC": (-4, 58, 16),
"VMPFC": (-4, 56, -16), "RS2": (60, -28, 38), "LS2": (-62, -32, 34), "RInsula": (42, 6, -6), "LInsula": (-42, -2, -4),
"RMFG": (50, 42, 12), "LMFG": (-46, 36, 14), "AMCC": (0, 2, 42), "rFFA": (38, -42, -22), "lFFA": (-40, -52, -18),
"right_primotor": (38, -24, 58), "left_primotor": (-38, -20, 58), "right_CS": (10, -86, 2), "left_CS": (-10, -86, 2)}
cen = []
for r in sb.ALL_ROIS:
xyz = sb.all_roi_mm()[r]
row = {"roi": r, "n": xyz.shape[1], "centroid": tuple(np.round(xyz.mean(axis=1), 1))}
if r in STATED:
dist = np.linalg.norm(xyz - np.array(STATED[r])[:, None], axis=0)
row.update({"stated center": STATED[r], "centroid to center mm": round(float(np.linalg.norm(xyz.mean(axis=1) - STATED[r])), 1),
"max distance mm": round(float(dist.max()), 1)})
cen.append(row)
print(pd.DataFrame(cen).to_string(index=False))
preprocessed image shapes: {(79, 95, 68, 168)} voxel sizes (mm): {(2.0, 2.0, 2.0)} dtype: {'int16'}
affines identical: True
168 volumes in every image (TR 2 s from task-pixar_bold.json): True
artifact indicator rows each flag one timepoint: True
CompCor rows per participant: {(5, 168)}
artifact timepoints per participant: max 48 of 168 (0.29; Methods drop at one third)
mean std min max
AgeGroup
3yo 9.82 12.59 0 39
4yo 9.36 8.52 0 26
5yo 13.06 12.62 0 48
7yo 8.91 8.68 0 31
8-12yo 9.85 9.35 0 35
Adult 2.82 4.00 0 15
ROI sizes (voxels): {'RTPJ': 376, 'LTPJ': 368, 'PC': 382, 'DMPFC': 217, 'MMPFC': 275, 'VMPFC': 198, 'RS2': 368, 'LS2': 269, 'RInsula': 309, 'LInsula': 240, 'RMFG': 142, 'LMFG': 256, 'AMCC': 249, 'rSTS_masked': 2503, 'lSTS': 844, 'rOFA': 790, 'lOFA': 211, 'rFFA': 1019, 'lFFA': 531, 'rRSC': 1741, 'lRSC': 1063, 'rTOS': 251, 'lTOS': 133, 'rPPA': 553, 'lPPA': 732, 'right_primotor': 515, 'left_primotor': 515, 'right_CS': 515, 'left_CS': 515}
Supplementary Table 2 sizes: RTPJ 376 LTPJ 368 PC 382 DMPFC 217 MMPFC 275 VMPFC 198; RS2 368 LS2 269 RInsula 309 LInsula 240 RMFG 142 LMFG 256 AMCC 249; RFFA 1019 LFFA 531
participants with any ROI voxel without data (zero or non-finite at some timepoint), per ROI:
RTPJ 0
LTPJ 0
PC 0
DMPFC 0
MMPFC 0
VMPFC 5
RS2 0
LS2 0
RInsula 0
LInsula 0
RMFG 0
LMFG 0
AMCC 0
rSTS_masked 0
lSTS 0
rOFA 16
lOFA 7
rFFA 19
lFFA 10
rRSC 0
lRSC 0
rTOS 0
lTOS 0
rPPA 2
lPPA 1
right_primotor 3
left_primotor 5
right_CS 0
left_CS 0
fewest voxels with data in a face ROI: 101 (Supplementary Fig. 4: every participant had at least 100 voxels in every face ROI)
participants with incomplete face-ROI coverage: 19 {'8-12yo': 14, 'Adult': 2, '4yo': 1, '5yo': 1, '7yo': 1} (Supplementary Fig. 4: 17, of them 2 adults, 12 8-12yo, one each 7yo, 4yo, 3yo)
NIfTI scl_slope, scl_inter: {('None', 'None')}
image value range: minimum -21.0 to 0.0, maximum 931.0 to 4081.0 over participants; sub-pixar001 (-6.0, 4056.0)
roi n centroid stated center centroid to center mm max distance mm
RTPJ 376 (48.0, -60.2, 29.8) (48, -60, 30) 0.3 8.9
LTPJ 368 (-47.7, -61.9, 30.0) (-48, -62, 30) 0.3 8.9
PC 382 (-0.0, -54.1, 34.1) (0, -54, 34) 0.1 8.9
DMPFC 217 (-6.6, 52.3, 35.4) (-6, 54, 36) 1.9 8.9
MMPFC 275 (-3.5, 57.6, 16.3) (-4, 58, 16) 0.7 8.9
VMPFC 198 (-2.9, 53.5, -14.7) (-4, 56, -16) 3.0 8.9
RS2 368 (59.6, -28.0, 37.9) (60, -28, 38) 0.4 8.9
LS2 269 (-60.1, -32.3, 33.1) (-62, -32, 34) 2.1 8.9
RInsula 309 (41.8, 5.5, -5.4) (42, 6, -6) 0.8 8.9
LInsula 240 (-40.8, -1.7, -3.6) (-42, -2, -4) 1.3 8.9
RMFG 142 (46.5, 42.1, 10.7) (50, 42, 12) 3.7 8.9
LMFG 256 (-44.6, 35.6, 14.0) (-46, 36, 14) 1.4 8.9
AMCC 249 (0.4, 2.1, 42.0) (0, 2, 42) 0.4 8.9
rSTS_masked 2503 (52.7, -41.0, 5.6) NaN NaN NaN
lSTS 844 (-56.2, -41.6, 7.6) NaN NaN NaN
rOFA 790 (40.1, -73.3, -14.4) NaN NaN NaN
lOFA 211 (-36.4, -75.0, -11.9) NaN NaN NaN
rFFA 1019 (39.7, -48.1, -22.1) (38, -42, -22) 6.3 22.7
lFFA 531 (-37.4, -57.2, -17.7) (-40, -52, -18) 5.9 19.4
rRSC 1741 (16.3, -55.7, 13.7) NaN NaN NaN
lRSC 1063 (-12.6, -57.1, 14.3) NaN NaN NaN
rTOS 251 (37.4, -77.8, 25.2) NaN NaN NaN
lTOS 133 (-31.7, -80.2, 26.3) NaN NaN NaN
rPPA 553 (23.2, -44.1, -10.7) NaN NaN NaN
lPPA 732 (-20.4, -45.2, -9.3) NaN NaN NaN
right_primotor 515 (38.0, -24.0, 58.0) (38, -24, 58) 0.0 10.0
left_primotor 515 (-38.0, -20.0, 58.0) (-38, -20, 58) 0.0 10.0
right_CS 515 (10.0, -86.0, 2.0) (10, -86, 2) 0.0 10.0
left_CS 515 (-10.0, -86.0, 2.0) (-10, -86, 2) 0.0 10.0
Timing check. There are no events files; the movie timing is fixed (README: TRs 1 to 5 rest, 6 to 10 logos, movie from TR 11). Two checks that the image volumes and the regressors share the same clock: (1) the ROI signal should deviate most at isolated artifact timepoints, compared with the same indicators shifted in time; (2) the adult ToM-network signal should be higher during the paper's ToM events than during its pain events (Supplementary Table 3), and the contrast should drop when the event labels are shifted in time; this is shown in the reverse correlation section, after the network timecourses are built.
In [7]:
def spike_at(s, shift):
# isolated artifact timepoints only; deviation of the mean ToM+Pain ROI signal from the median
# of the 7 timepoints around it, in units of that participant's median deviation
x = np.mean([D[s]["rois"][r]["all"] for r in sb.TOM_ROIS + sb.PAIN_ROIS], axis=0)
med = pd.Series(x).rolling(7, center=True, min_periods=4).median().values
dev = np.abs(x - med); dev = dev / np.median(dev)
a = sb.artifact_index(D[s]["artifacts"])
iso = [t for t in a if t - 1 not in a and t + 1 not in a and 3 <= t + shift < sb.N_TR - 3]
return [dev[t + shift] for t in iso]
for shift in [-2, -1, 0, 1, 2, 10]:
vals = np.concatenate([spike_at(s, shift) for s in subs])
print(f"isolated artifact timepoints shifted by {shift:3d} TRs: n = {len(vals)}, median deviation = {np.median(vals):.2f}")
isolated artifact timepoints shifted by -2 TRs: n = 99, median deviation = 1.44 isolated artifact timepoints shifted by -1 TRs: n = 99, median deviation = 2.14
isolated artifact timepoints shifted by 0 TRs: n = 98, median deviation = 6.45
isolated artifact timepoints shifted by 1 TRs: n = 96, median deviation = 2.58 isolated artifact timepoints shifted by 2 TRs: n = 94, median deviation = 0.91
isolated artifact timepoints shifted by 10 TRs: n = 90, median deviation = 1.00
In [8]:
counts = {"children": len(children), "adults": len(adults), "female children": int((children.Gender == "F").sum()),
"female adults": int((adults.Gender == "F").sum())}
counts.update({f"children {g}": int((children.AgeGroup == g).sum()) for g in sb.AGE_GROUPS})
counts["mean child age"] = round(children.Age.mean(), 2)
expected = [122, 33, 64, 20, 17, 14, 34, 23, 34, 6.7]
print(pd.DataFrame({"measured": pd.Series(counts), "expected": expected}).to_string())
print("Age-group ranges:", children.groupby("AgeGroup").Age.agg(["min", "max"]).round(2).to_dict("index"))
print("The 5 children excluded for motion are not in the dataset (155 participants listed).")
measured expected
children 122.00 122.0
adults 33.00 33.0
female children 64.00 64.0
female adults 20.00 20.0
children 3yo 17.00 17.0
children 4yo 14.00 14.0
children 5yo 34.00 34.0
children 7yo 23.00 23.0
children 8-12yo 34.00 34.0
mean child age 6.71 6.7
Age-group ranges: {'3yo': {'min': 3.52, 'max': 3.99}, '4yo': {'min': 4.06, 'max': 4.86}, '5yo': {'min': 5.01, 'max': 5.99}, '7yo': {'min': 7.0, 'max': 7.96}, '8-12yo': {'min': 8.02, 'max': 12.3}}
The 5 children excluded for motion are not in the dataset (155 participants listed).
Result: reproduced. All nine counts match exactly (122 children, 33 adults, 64 and 20 female; 17, 14, 34, 23 and 34 children by age group) and the mean child age is 6.71 years against 6.7. The age-group ranges match Supplementary Table 1. The five children excluded for motion are not in the dataset.
supp-table1-behavior-by-age¶
Mean ToM score rises across age groups (3yo .54, 4yo .63, 5yo .73, 7yo .88, 8-12yo .96), and the false-belief groups (pass/inconsistent/fail) are 4/4/9 at 3 years, 3/7/4 at 4 years and 23/9/2 at 5 years.
In [9]:
m = children.groupby("AgeGroup")[TOM].agg(["mean", "std", "count"]).loc[sb.AGE_GROUPS]
m["expected"] = [0.54, 0.63, 0.73, 0.88, 0.96]
m["within 10%"] = (m["mean"] - m.expected).abs() <= 0.1 * m.expected
print(m.round(3).to_string())
fb = pd.crosstab(young.AgeGroup, young.FB_Group)[["pass", "inc", "fail"]]
print(fb.to_string())
print("expected 3yo 4/4/9, 4yo 3/7/4, 5yo 23/9/2; totals", fb.sum().to_dict(), "(expected 30, 20, 15)")
mean std count expected within 10%
AgeGroup
3yo 0.517 0.178 17 0.54 True
4yo 0.610 0.152 14 0.63 True
5yo 0.724 0.154 34 0.73 True
7yo 0.875 0.080 23 0.88 True
8-12yo 0.956 0.063 34 0.96 True
FB_Group pass inc fail
AgeGroup
3yo 4 4 9
4yo 3 7 4
5yo 23 9 2
expected 3yo 4/4/9, 4yo 3/7/4, 5yo 23/9/2; totals {'pass': 30, 'inc': 20, 'fail': 15} (expected 30, 20, 15)
Result: reproduced. The five mean ToM scores (.517, .610, .724, .875, .956) are within 10% of .54, .63, .73, .88 and .96, and the nine false-belief group counts by age (4/4/9, 3/7/4, 23/9/2) match exactly.
beh-tom-age¶
Among children, ToM score (proportion correct) increases with age (Kendall tau = .66), and among the 3- to 5-year-olds DCCS summary score also increases with age (Kendall tau = .20).
In [10]:
t1 = sb.kendall(children[TOM], children.Age)
yd = young.dropna(subset=[DCCS])
t2 = sb.kendall(yd[DCCS], yd.Age)
print(f"ToM score ~ age, n = {len(children)}: tau-b = {t1[0]:.3f}, p = {t1[1]:.2g} (paper .66, p < .00001)")
print(f"DCCS ~ age, n = {len(yd)}: tau-b = {t2[0]:.3f}, p = {t2[1]:.3f} (paper .20, p = .049)")
for v in ["c"]:
print(f" variant tau-{v}:", round(stats.kendalltau(children[TOM], children.Age, variant=v).statistic, 3),
round(stats.kendalltau(yd[DCCS], yd.Age, variant=v).statistic, 3))
ToM score ~ age, n = 122: tau-b = 0.665, p = 5.4e-26 (paper .66, p < .00001) DCCS ~ age, n = 64: tau-b = 0.198, p = 0.048 (paper .20, p = .049) variant tau-c: 0.654 0.191
Result: reproduced. ToM score against age: tau-b = .665 (paper .66), p = 5.4e-26. DCCS against age in the 64 children with a DCCS score: tau-b = .198 (paper .20), p = .048 (paper .049, same side of .05).
beh-tom-dccs¶
In 3- to 5-year-olds, ToM and DCCS scores are positively correlated controlling for age (partial Kendall tau = .19).
In [11]:
r, pv = sb.partial_kendall(yd[TOM], yd[DCCS], yd.Age)
print(f"partial tau(ToM, DCCS | age), n = {len(yd)}: {r:.3f}, p = {pv:.3f} (paper .19, p = .03)")
partial tau(ToM, DCCS | age), n = 64: 0.186, p = 0.031 (paper .19, p = .03)
Result: reproduced. Partial Kendall tau of ToM and DCCS given age = .186 (paper .19, tolerance .019), p = .031 (paper .03).
beh-fb-groups¶
Among the 65 children aged 3 to 5, the false-belief composite gives 30 passers (15 female), 20 inconsistent performers (13 female) and 15 failers (6 female).
In [12]:
g = young.groupby("FB_Group").agg(n=("Age", "size"), female=("Gender", lambda x: int((x == "F").sum())),
age_mean=("Age", "mean"), age_sd=("Age", "std"))
print(g.round(2).to_string())
recut = pd.cut(young.FB_Composite, [-1, 2, 4, 6], labels=["fail", "inc", "pass"]).astype(str)
print("groups recomputed from FB_Composite agree with FB_Group:", bool((recut == young.FB_Group).all()))
print("handedness by group:", young.groupby("FB_Group").Handedness.value_counts().to_dict())
n female age_mean age_sd
FB_Group
fail 15 6 4.08 0.56
inc 20 13 4.75 0.73
pass 30 15 5.23 0.70
groups recomputed from FB_Composite agree with FB_Group: True
handedness by group: {('fail', 'R'): 11, ('fail', 'L'): 4, ('inc', 'R'): 17, ('inc', 'L'): 2, ('inc', 'Ambi'): 1, ('pass', 'R'): 26, ('pass', 'Ambi'): 2, ('pass', 'L'): 2}
Result: reproduced. 30 passers (15 female), 20 inconsistent (13 female), 15 failers (6 female), as deposited in FB_Group; the groups recomputed from FB_Composite with the Methods' cut-offs are the same.
beh-fb-tom, beh-fb-age, beh-fb-dccs¶
False-belief failers and inconsistent performers score lower than passers on the remaining ToM items (means .55, .57 and .75). False-belief failers are younger than inconsistent performers and passers (mean ages 4.1, 4.8 and 5.2 years). False-belief failers have lower DCCS scores than the other two groups (means 1.73, 2.26, 2.33).
Tukey HSD over the three groups (statsmodels pairwise_tukeyhsd), diffs under both scales of O9;
Kruskal-Wallis from scipy.
In [13]:
def fb_tests(var):
sub = young.dropna(subset=[var])
out = sub.groupby("FB_Group")[var].agg(["mean", "sem", "count"]).loc[["fail", "inc", "pass"]]
print(out.round(3).to_string())
kw = stats.kruskal(*[sub.loc[sub.FB_Group == g, var] for g in ["pass", "inc", "fail"]])
print(f"Kruskal-Wallis H = {kw.statistic:.2f}, p = {kw.pvalue:.2g}")
for scale in ["z-scored (O9a)", "raw (O9b)"]:
v = (sub[var] - sub[var].mean()) / sub[var].std() if scale.startswith("z") else sub[var]
t = pairwise_tukeyhsd(v.values, sub.FB_Group.values)
df = pd.DataFrame(t.summary().data[1:], columns=t.summary().data[0])
df["p-adj (exact)"] = [f"{x:.2g}" for x in t.pvalues]
print(" Tukey,", scale); print(df.to_string(index=False))
return out, kw
print("== beh-fb-tom, reading (a): 18 non-FB items (ToM Booklet-Matched-NOFB)")
fbtom_a = fb_tests(NOFB)
print("\n== beh-fb-tom, reading (b): full ToM score")
fbtom_b = fb_tests(TOM)
== beh-fb-tom, reading (a): 18 non-FB items (ToM Booklet-Matched-NOFB)
mean sem count
FB_Group
fail 0.546 0.043 15
inc 0.572 0.034 20
pass 0.753 0.021 30
Kruskal-Wallis H = 22.95, p = 1e-05
Tukey,
z-scored (O9a) group1 group2 meandiff p-adj lower upper reject p-adj (exact) fail inc 0.1561 0.8455 -0.5222 0.8344 False 0.85 fail pass 1.2399 0.0000 0.6120 1.8679 True 3.8e-05 inc pass 1.0838 0.0001 0.5106 1.6571 True 7.8e-05
Tukey, raw (O9b)
group1 group2 meandiff p-adj lower upper reject p-adj (exact)
fail inc 0.0262 0.8455 -0.0875 0.1398 False 0.85
fail pass 0.2077 0.0000 0.1025 0.3129 True 3.8e-05
inc pass 0.1816 0.0001 0.0855 0.2776 True 7.8e-05
== beh-fb-tom, reading (b): full ToM score
mean sem count
FB_Group
fail 0.451 0.037 15
inc 0.571 0.026 20
pass 0.792 0.017 30
Kruskal-Wallis H = 41.51, p = 9.7e-10
Tukey, z-scored (O9a) group1 group2 meandiff p-adj lower upper reject p-adj (exact) fail inc 0.6606 0.0073 0.1543 1.1670 True 0.0073 fail pass 1.8833 0.0000 1.4145 2.3521 True 1.6e-13 inc pass 1.2227 0.0000 0.7947 1.6506 True 1.1e-08
Tukey, raw (O9b) group1 group2 meandiff p-adj lower upper reject p-adj (exact) fail inc 0.1197 0.0073 0.0280 0.2114 True 0.0073 fail pass 0.3412 0.0000 0.2563 0.4261 True 1.6e-13 inc pass 0.2215 0.0000 0.1440 0.2990 True 1.1e-08
Result: reproduced. Under reading (a), the 18 non-false-belief items deposited as ToM Booklet-Matched-NOFB, the means are .546, .572 and .753 (paper .55, .57, .75) and Kruskal-Wallis H = 22.95 (paper 22.96), p = 1e-05. The Tukey tests match: Pass-Fail diff 1.24 on the z scale (paper 1.2), p = 3.8e-05; Pass-Inc 1.08 (paper 1.08), p = 7.8e-05; Inc-Fail .16 (paper .16), p = .85. The z-scaled diffs (O9a) are the ones that match the paper's. Reading (b), the full score, gives means .451, .571, .792, out of tolerance for failers and passers.
In [14]:
fbage = fb_tests("Age")
mean sem count FB_Group fail 4.082 0.145 15 inc 4.751 0.164 20 pass 5.228 0.127 30 Kruskal-Wallis H = 20.35, p = 3.8e-05
Tukey, z-scored (O9a) group1 group2 meandiff p-adj lower upper reject p-adj (exact) fail inc 0.8263 0.0149 0.1372 1.5154 True 0.015 fail pass 1.4157 0.0000 0.7777 2.0537 True 4.3e-06 inc pass 0.5894 0.0467 0.0070 1.1718 True 0.047
Tukey, raw (O9b) group1 group2 meandiff p-adj lower upper reject p-adj (exact) fail inc 0.6688 0.0149 0.1110 1.2266 True 0.015 fail pass 1.1459 0.0000 0.6295 1.6623 True 4.3e-06 inc pass 0.4771 0.0467 0.0057 0.9485 True 0.047
Result: reproduced. Mean ages 4.08, 4.75 and 5.23 (paper 4.1, 4.8, 5.2). Tukey on the z scale: Pass-Fail diff 1.42 (paper 1.4), p = 4.3e-06; Inc-Fail .83 (paper .83), p = .015 (paper .01); Pass-Inc .59 (paper .59), p = .047 (paper .047).
In [15]:
fbdccs = fb_tests(DCCS)
mean sem count FB_Group fail 1.733 0.206 15 inc 2.263 0.168 19 pass 2.333 0.088 30 Kruskal-Wallis H = 7.56, p = 0.023
Tukey, z-scored (O9a) group1 group2 meandiff p-adj lower upper reject p-adj (exact) fail inc 0.7794 0.0524 -0.0064 1.5653 False 0.052 fail pass 0.8827 0.0124 0.1632 1.6022 True 0.012 inc pass 0.1032 0.9267 -0.5639 0.7703 False 0.93
Tukey, raw (O9b) group1 group2 meandiff p-adj lower upper reject p-adj (exact) fail inc 0.5298 0.0524 -0.0044 1.0640 False 0.052 fail pass 0.6000 0.0124 0.1109 1.0891 True 0.012 inc pass 0.0702 0.9267 -0.3833 0.5236 False 0.93
Result: reproduced. Mean DCCS 1.733, 2.263 and 2.333 (paper 1.73, 2.26, 2.33), Kruskal-Wallis H = 7.56, p = .023 (paper 7.56, .02); Tukey Pass-Fail diff .88 on the z scale (paper .88), p = .012; Inc-Fail .78 (paper .78), p = .052 (paper .052, reported, not deciding).
In [16]:
w = stats.ttest_ind(children.n_art, adults.n_art, equal_var=False)
print(p.groupby("Child_Adult").n_art.agg(["mean", "std", "count"]).round(2).to_string())
print(f"Welch t = {w.statistic:.2f}, df = {w.df:.1f}, p = {w.pvalue:.2g} (paper t(137.7) = 6.49, p < .000001)")
mean std count Child_Adult adult 2.82 4.0 33 child 10.51 10.6 122 Welch t = 6.49, df = 137.7, p = 1.4e-09 (paper t(137.7) = 6.49, p < .000001)
Result: reproduced. Mean artifact timepoints 10.51 (s.d. 10.6) in children and 2.82 (s.d. 4.0) in adults (paper 10.5 and 2.8); Welch t = 6.49 on 137.7 df, p = 1.4e-09, as in the paper.
motion-not-age-tom¶
Among children, the number of artifact timepoints is not correlated with age or with ToM score.
In [17]:
rs = stats.spearmanr(children.n_art, children.Age)
rk = sb.kendall(children.n_art, children[TOM])
print(f"Spearman(n_art, age) = {rs.statistic:.3f}, p = {rs.pvalue:.2f} (paper .02, p = .86)")
print(f"Kendall(n_art, ToM score) = {rk[0]:.3f}, p = {rk[1]:.2f} (paper -.005, p = .94)")
Spearman(n_art, age) = 0.017, p = 0.86 (paper .02, p = .86) Kendall(n_art, ToM score) = -0.005, p = 0.94 (paper -.005, p = .94)
Result: reproduced. Neither correlation is significant: Spearman with age .017, p = .86; Kendall with ToM score -.005, p = .94, the paper's values.
motion-not-fb¶
In 3- to 5-year-olds, the number of artifact timepoints does not differ between false-belief passers and failers, nor vary with DCCS score.
In [18]:
pf = young[young.FB_Group.isin(["pass", "fail"])].copy()
pf["fail"] = (pf.FB_Group == "fail").astype(float)
for z in [False, True]:
r1 = sb.ols(pf, "n_art", ["fail"], zscore=z)
r2 = sb.ols(yd, "n_art", [DCCS], zscore=z)
print(f"{'z-scored (O4a)' if z else 'raw (O4b)'}: FB group (fail = 1, n = {int(r1.nobs)}) b = {r1.params['fail']:.3f}, t = {r1.tvalues['fail']:.2f}, p = {r1.pvalues['fail']:.2f};"
f" DCCS (n = {int(r2.nobs)}) b = {r2.params[DCCS]:.3f}, t = {r2.tvalues[DCCS]:.2f}, p = {r2.pvalues[DCCS]:.2f}")
print("paper: FB b = -.04, t = -.12, p = .9; DCCS b = .16, t = 1.18, p = .25")
raw (O4b): FB group (fail = 1, n = 45) b = 0.433, t = 0.12, p = 0.90; DCCS (n = 64) b = 2.570, t = 1.17, p = 0.24 z-scored (O4a): FB group (fail = 1, n = 45) b = 0.039, t = 0.12, p = 0.90; DCCS (n = 64) b = 0.148, t = 1.17, p = 0.24 paper: FB b = -.04, t = -.12, p = .9; DCCS b = .16, t = 1.18, p = .25
Result: reproduced. Neither effect is significant: FB group t = .12, p = .90 (z-scaled b = .039; the paper's sign is the opposite coding), DCCS t = 1.17, p = .24 (z-scaled b = .148, paper .16).
Inter-region correlation¶
Pipeline P4 and P5 on the deposited preprocessed images. For each participant, every ROI
timecourse is cleaned (sb.clean) and correlated with every other. Readings run: O1 (all voxels
of the deposited masks, all, or only voxels with data, with_data; the deposited masks match
Supplementary Table 2 so reading O1b, full spheres, is not needed), O2 (one indicator per
artifact timepoint, a, or one regressor, b) and O3 (all 168 TRs, a, or TRs 11 to 168, b). The
main reading is O1 all, O2a (the deposited regressors), O3a.
Where the readings of a pipeline detail give different outcomes, the outcome is taken under the reading the paper best supports, adopted as an assumption (protocol, author step 1.6), and every reading is shown. This applies to the claims of absence too (irc-tom-behavior-age, irc-fb-across, maturity-within, rc-fb-no-difference): the claims file names one test for each, but not the pipeline readings, so the assumed readings decide (O2a and O3a for the inter-region measures, O7a for event magnitudes), and the outcome is marked as resting on them.
In [19]:
READINGS = [(v, s, m) for v in ["all", "with_data"] for s in [False, True] for m in [False, True]]
def reading_name(r):
return f"O1 {r[0]}, O2{'b' if r[1] else 'a'}, O3{'b' if r[2] else 'a'}"
MAIN = ("all", False, False)
TP = sb.TOM_ROIS + sb.PAIN_ROIS
BILATERAL = {"a": [("RTPJ", "LTPJ"), ("RS2", "LS2"), ("RInsula", "LInsula"), ("RMFG", "LMFG")],
"b": [("RTPJ", "LTPJ"), ("RInsula", "LInsula")]}
def irc_measures(reading, m1):
voxels, single, movie = reading
trs = np.arange(sb.FIRST_MOVIE_TR - 1, sb.N_TR) if movie else None
rois = TP if m1 else sb.ALL_ROIS
rows, mats = {}, {}
for s in subs:
tc = sb.subject_timecourses(D[s], rois, voxels, m1, single)
m = sb.irc_matrix(tc, rois, trs)
mats[s] = m
row = {"wToM": sb.within(m, sb.TOM_ROIS), "wPain": sb.within(m, sb.PAIN_ROIS),
"across": sb.across(m, sb.TOM_ROIS, sb.PAIN_ROIS)}
for k, pairs in BILATERAL.items():
row[f"wToM_nonbi_{k}"] = sb.within(m, sb.TOM_ROIS, pairs)
row[f"wPain_nonbi_{k}"] = sb.within(m, sb.PAIN_ROIS, pairs)
if not m1:
for net in ["Face", "Scene", "M1", "V1"]:
row["w" + net] = sb.within(m, sb.NETWORKS[net])
row["aFaceScene"] = sb.across(m, sb.FACE_ROIS, sb.SCENE_ROIS)
row["aM1V1"] = sb.across(m, sb.M1_ROIS, sb.V1_ROIS)
rows[s] = row
return pd.DataFrame(rows).T.join(p), mats
IRC = {r: irc_measures(r, True) for r in READINGS} # M1 regressed: main-text analyses
IRC_NOM1 = {r: irc_measures(r, False) for r in READINGS} # no M1 regression: Supplementary Fig. 3
irc = IRC[MAIN][0]
print(irc.groupby("AgeGroup")[["wToM", "wPain", "across"]].agg(["mean", "sem"]).round(3).to_string())
print("ranges:", irc[["wToM", "wPain", "across"]].agg(["min", "max"]).round(3).to_dict())
wToM wPain across
mean sem mean sem mean sem
AgeGroup
3yo 0.326 0.027 0.179 0.021 0.074 0.020
4yo 0.361 0.042 0.217 0.030 0.064 0.019
5yo 0.417 0.021 0.246 0.011 0.083 0.015
7yo 0.494 0.023 0.284 0.016 0.025 0.014
8-12yo 0.435 0.021 0.316 0.017 0.028 0.017
Adult 0.604 0.021 0.398 0.017 -0.051 0.021
ranges: {'wToM': {'min': 0.106, 'max': 0.85}, 'wPain': {'min': 0.042, 'max': 0.669}, 'across': {'min': -0.293, 'max': 0.281}}
In [20]:
fig, axes = plt.subplots(1, 3, figsize=(12, 3.2), sharey=True)
order = sb.AGE_GROUPS + ["Adult"]
for ax, v in zip(axes, ["wToM", "wPain", "across"]):
ax.boxplot([irc.loc[irc.AgeGroup == g, v] for g in order], tick_labels=order)
ax.axhline(0, color="grey", lw=0.5); ax.set_title(v)
axes[0].set_ylabel("Fisher z (main reading)")
plt.tight_layout(); plt.show()
print("all values shown; y range", round(irc[["wToM", "wPain", "across"]].min().min(), 3), "to", round(irc[["wToM", "wPain", "across"]].max().max(), 3))
all values shown; y range -0.293 to 0.85
irc-adult¶
In adults, within-ToM (mean z .48) and within-Pain (.35) correlations are strongly positive and the across-network correlation is negative (-.17), each within-network value exceeding the across-network value.
In [21]:
def tol(v, e, t):
return abs(v - e) <= t + 1e-12
rows = []
for r in READINGS:
a = IRC[r][0].loc[adults.index]
t1 = stats.ttest_rel(a.wToM, a.across); t2 = stats.ttest_rel(a.wPain, a.across)
rows.append({"reading": reading_name(r), "wToM": a.wToM.mean(), "wPain": a.wPain.mean(), "across": a.across.mean(),
"t ToM-across": t1.statistic, "p1": t1.pvalue, "t Pain-across": t2.statistic, "p2": t2.pvalue})
res = pd.DataFrame(rows)
res["values in tolerance"] = [sum([tol(x.wToM, .48, .02), tol(x.wPain, .35, .02), tol(x.across, -.17, .02),
tol(x["t ToM-across"], 19.1, 1.91), tol(x["t Pain-across"], 23.2, 2.32)]) for _, x in res.iterrows()]
print(res.to_string(float_format=lambda x: f"{x:.3g}"))
print("expected: .48, .35, -.17, t 19.1, t 23.2 (5 values)")
reading wToM wPain across t ToM-across p1 t Pain-across p2 values in tolerance 0 O1 all, O2a, O3a 0.604 0.398 -0.051 19.1 4.97e-19 21.8 9.3e-21 2 1 O1 all, O2a, O3b 0.594 0.397 -0.0525 18.9 6.64e-19 21 2.91e-20 2 2 O1 all, O2b, O3a 0.606 0.402 -0.0475 18.9 5.98e-19 21.6 1.21e-20 2 3 O1 all, O2b, O3b 0.597 0.401 -0.0483 18.7 8.47e-19 20.7 4.38e-20 1 4 O1 with_data, O2a, O3a 0.604 0.398 -0.051 19.1 4.97e-19 21.8 9.3e-21 2 5 O1 with_data, O2a, O3b 0.594 0.397 -0.0525 18.9 6.64e-19 21 2.91e-20 2 6 O1 with_data, O2b, O3a 0.606 0.402 -0.0475 18.9 5.98e-19 21.6 1.21e-20 2 7 O1 with_data, O2b, O3b 0.597 0.401 -0.0483 18.7 8.47e-19 20.7 4.38e-20 1 expected: .48, .35, -.17, t 19.1, t 23.2 (5 values)
In [22]:
# Variants that are not the paper's method, to show how far the adult and 3yo means move with the
# filter and the CompCor regressors (main reading otherwise).
def means_with(clean_fn, idx):
vals = []
for s in idx:
tc = {r: clean_fn(s, r) for r in TP}
m = sb.irc_matrix(tc, TP)
vals.append([sb.within(m, sb.TOM_ROIS), sb.within(m, sb.PAIN_ROIS), sb.across(m, sb.TOM_ROIS, sb.PAIN_ROIS)])
return np.round(np.mean(vals, axis=0), 3)
def m1_of(s):
r, l = D[s]["rois"]["right_primotor"], D[s]["rois"]["left_primotor"]
return (r["all"] * r["n_vox"] + l["all"] * l["n_vox"]) / (r["n_vox"] + l["n_vox"])
def variant(hp_cutoff=sb.HIGHPASS_S, compcor=True):
def f(s, r):
d = D[s]; bad = sb.artifact_index(d["artifacts"])
y = sb.nearest_neighbour_fill(d["rois"][r]["all"], bad) if len(bad) else d["rois"][r]["all"].copy()
parts = ([d["compcor"].T] if compcor else []) + ([d["artifacts"].T] if len(bad) else []) + [m1_of(s)[:, None]]
y = sb.residualize(y, np.column_stack(parts))
if hp_cutoff:
hp = sb.dct_highpass_basis(sb.N_TR, sb.TR_S, hp_cutoff); y = y - hp @ np.linalg.lstsq(hp, y, rcond=None)[0]
y[bad] = np.nan
return y
return f
y3_idx = children[children.AgeGroup == "3yo"].index
for name, fn in {"paper (100 s, CompCor)": variant(), "no high-pass": variant(None), "50 s high-pass": variant(50.0), "no CompCor": variant(compcor=False)}.items():
print(f"{name:24s} adults (wToM, wPain, across) {means_with(fn, adults.index)} 3yo {means_with(fn, y3_idx)}")
print("paper: adults .48, .35, -.17; 3yo .21, .23, .05")
paper (100 s, CompCor) adults (wToM, wPain, across) [ 0.604 0.398 -0.051] 3yo [0.326 0.179 0.074]
no high-pass adults (wToM, wPain, across) [0.626 0.397 0.032] 3yo [0.334 0.19 0.097]
50 s high-pass adults (wToM, wPain, across) [ 0.556 0.38 -0.017] 3yo [0.333 0.17 0.072]
no CompCor adults (wToM, wPain, across) [0.645 0.465 0.059] 3yo [0.401 0.286 0.187] paper: adults .48, .35, -.17; 3yo .21, .23, .05
In [23]:
# Where the offset from the paper lies. The mean 13 x 13 matrices (main reading) of adults and
# 3-year-olds; Fig. 2 shows these matrices as colours only, with no numbers, so a pair-by-pair
# comparison with the paper is not possible.
for name, idx in [("adults", adults.index), ("3yo", y3_idx)]:
M = pd.DataFrame(np.mean([IRC[MAIN][1][s].values for s in idx], axis=0), index=TP, columns=TP)
print(f"mean Fisher z, {name} (n = {len(idx)}), main reading:")
print(M.round(2).to_string())
print()
# Diagnostic variants, none of them the paper's stated method: averaging r in place of Fisher z;
# the high-pass basis in the same regression as the nuisance covariates (as in an SPM GLM); left
# and right M1 as two regressors; and, to test whether a shared signal accounts for the offset,
# the mean ToM and Pain signal added as a regressor.
def diag_variant(kind):
def f(s, r):
d = D[s]; bad = sb.artifact_index(d["artifacts"])
y = sb.nearest_neighbour_fill(d["rois"][r]["all"], bad) if len(bad) else d["rois"][r]["all"].copy()
art = [d["artifacts"].T] if len(bad) else []
m1 = np.column_stack([d["rois"][k]["all"] for k in sb.M1_ROIS]) if kind == "M1 left and right" else m1_of(s)[:, None]
extra = [m1]
if kind == "plus mean ToM and Pain signal":
extra.append(np.mean([d["rois"][k]["all"] for k in TP], axis=0)[:, None])
hp = sb.dct_highpass_basis()
if kind == "high-pass in the nuisance regression":
y = sb.residualize(y, np.column_stack([d["compcor"].T] + art + extra + [hp]))
else:
y = sb.residualize(y, np.column_stack([d["compcor"].T] + art + extra))
y = y - hp @ np.linalg.lstsq(hp, y, rcond=None)[0]
y[bad] = np.nan
return y
return f
def mean_r(idx):
v = [[sb.within(np.tanh(IRC[MAIN][1][s]), sb.TOM_ROIS), sb.within(np.tanh(IRC[MAIN][1][s]), sb.PAIN_ROIS),
sb.across(np.tanh(IRC[MAIN][1][s]), sb.TOM_ROIS, sb.PAIN_ROIS)] for s in idx]
return np.round(np.mean(v, axis=0), 3)
print(f"{'mean of r, not z':40s} adults (wToM, wPain, across) {mean_r(adults.index)} 3yo {mean_r(y3_idx)}")
for kind in ["high-pass in the nuisance regression", "M1 left and right", "plus mean ToM and Pain signal"]:
fn = diag_variant(kind)
print(f"{kind:40s} adults (wToM, wPain, across) {means_with(fn, adults.index)} 3yo {means_with(fn, y3_idx)}")
print("paper: adults .48, .35, -.17; 3yo .21, .23, .05")
mean Fisher z, adults (n = 33), main reading:
RTPJ LTPJ PC DMPFC MMPFC VMPFC RS2 LS2 RInsula LInsula RMFG LMFG AMCC
RTPJ NaN 0.78 0.83 0.38 0.44 0.47 -0.16 -0.16 -0.10 -0.12 0.10 0.05 -0.09
LTPJ 0.78 NaN 0.81 0.68 0.63 0.53 -0.23 -0.13 -0.11 -0.07 -0.01 0.10 -0.07
PC 0.83 0.81 NaN 0.46 0.53 0.57 -0.22 -0.17 -0.14 -0.13 -0.03 -0.02 -0.06
DMPFC 0.38 0.68 0.46 NaN 0.82 0.53 -0.16 -0.09 -0.05 -0.05 0.01 0.10 0.00
MMPFC 0.44 0.63 0.53 0.82 NaN 0.60 -0.10 -0.03 0.03 0.01 0.04 0.06 0.02
VMPFC 0.47 0.53 0.57 0.53 0.60 NaN -0.16 -0.11 0.03 0.04 -0.01 0.03 0.03
RS2 -0.16 -0.23 -0.22 -0.16 -0.10 -0.16 NaN 0.94 0.38 0.37 0.44 0.36 0.36
LS2 -0.16 -0.13 -0.17 -0.09 -0.03 -0.11 0.94 NaN 0.29 0.38 0.33 0.46 0.32
RInsula -0.10 -0.11 -0.14 -0.05 0.03 0.03 0.38 0.29 NaN 0.92 0.31 0.19 0.48
LInsula -0.12 -0.07 -0.13 -0.05 0.01 0.04 0.37 0.38 0.92 NaN 0.24 0.24 0.45
RMFG 0.10 -0.01 -0.03 0.01 0.04 -0.01 0.44 0.33 0.31 0.24 NaN 0.60 0.14
LMFG 0.05 0.10 -0.02 0.10 0.06 0.03 0.36 0.46 0.19 0.24 0.60 NaN 0.16
AMCC -0.09 -0.07 -0.06 0.00 0.02 0.03 0.36 0.32 0.48 0.45 0.14 0.16 NaN
mean Fisher z, 3yo (n = 17), main reading:
RTPJ LTPJ PC DMPFC MMPFC VMPFC RS2 LS2 RInsula LInsula RMFG LMFG AMCC
RTPJ NaN 0.56 0.64 0.09 0.17 0.28 0.14 0.03 -0.04 -0.07 0.18 0.13 -0.02
LTPJ 0.56 NaN 0.64 0.16 0.18 0.18 0.08 0.22 -0.03 -0.10 0.14 0.14 -0.02
PC 0.64 0.64 NaN 0.15 0.26 0.29 0.11 0.11 -0.05 -0.09 0.10 0.13 -0.07
DMPFC 0.09 0.16 0.15 NaN 0.61 0.23 0.09 0.06 0.05 0.04 0.14 0.24 0.24
MMPFC 0.17 0.18 0.26 0.61 NaN 0.47 0.21 0.19 -0.02 -0.05 0.21 0.35 -0.01
VMPFC 0.28 0.18 0.29 0.23 0.47 NaN 0.13 0.08 0.00 0.03 0.05 0.10 -0.04
RS2 0.14 0.08 0.11 0.09 0.21 0.13 NaN 0.66 0.21 0.13 0.19 0.14 0.17
LS2 0.03 0.22 0.11 0.06 0.19 0.08 0.66 NaN 0.18 0.16 0.04 0.17 0.14
RInsula -0.04 -0.03 -0.05 0.05 -0.02 0.00 0.21 0.18 NaN 0.86 0.02 -0.05 0.27
LInsula -0.07 -0.10 -0.09 0.04 -0.05 0.03 0.13 0.16 0.86 NaN -0.10 -0.10 0.31
RMFG 0.18 0.14 0.10 0.14 0.21 0.05 0.19 0.04 0.02 -0.10 NaN 0.34 0.04
LMFG 0.13 0.14 0.13 0.24 0.35 0.10 0.14 0.17 -0.05 -0.10 0.34 NaN -0.00
AMCC -0.02 -0.02 -0.07 0.24 -0.01 -0.04 0.17 0.14 0.27 0.31 0.04 -0.00 NaN
mean of r, not z adults (wToM, wPain, across) [ 0.519 0.354 -0.049] 3yo [0.291 0.159 0.068]
high-pass in the nuisance regression adults (wToM, wPain, across) [ 0.593 0.393 -0.086] 3yo [0.317 0.195 0.07 ]
M1 left and right adults (wToM, wPain, across) [ 0.599 0.39 -0.037] 3yo [0.323 0.177 0.073]
plus mean ToM and Pain signal adults (wToM, wPain, across) [ 0.286 0.258 -0.382] 3yo [ 0.029 0.035 -0.168] paper: adults .48, .35, -.17; 3yo .21, .23, .05
Result: partially reproduced. Within-ToM and within-Pain correlations are strongly positive in adults and above the across-network value (t = 19.1 and 21.8 under the main reading, both p < 1e-18), but the means are .604, .398 and -.051 against .48, .35 and -.17, all three out of tolerance under every reading; the two t values are within tolerance under three of the four O2 and O3 readings. The across-network mean is negative (-.051) but much weaker than the paper's. None of the variants outside the paper's method (no high-pass, 50 s high-pass, no CompCor) brings the means to the paper's. The adult and 3yo mean matrices are printed above; Fig. 2 gives them only as colours, so the offset cannot be traced pair by pair. The masks are in MNI mm (centroids of the 9 mm spheres within 3.7 mm of the stated centers and every voxel within 8.9 mm; M1 and V1 centered exactly, radius 10 mm), and the images are int16 with no scaling (maximum value 931 to 4081 across participants). Of the diagnostic variants, averaging r in place of z gives adults .519, .354, -.049, the high-pass in the nuisance regression .593, .393, -.086, and left and right M1 as two regressors .599, .390, -.037; only adding the mean ToM and Pain signal as a regressor moves the across-network value past the paper's (.286, .258, -.382). A shared signal that the paper's processing removed and this one does not would give an offset in this direction, but no step of the Methods' pipeline accounts for it.
irc-age-within and irc-age-across¶
Among children, within-ToM and within-Pain correlations increase with age, controlling for motion. Among children, the across-network correlation decreases with age.
Partial Spearman correlation with the number of artifact timepoints as covariate (ppcor form).
In [24]:
rows = []
for r in READINGS:
c = IRC[r][0].loc[children.index]
row = {"reading": reading_name(r)}
for v in ["wToM", "wPain", "across"]:
rr, pv, df = sb.partial_spearman(c[v], c.Age, c.n_art)
row[v + " rs"], row[v + " p"] = rr, pv
rows.append(row)
age_irc = pd.DataFrame(rows)
print(age_irc.to_string(float_format=lambda x: f"{x:.3g}"))
print("df =", df, "; expected wToM .37 (p < .00005), wPain .28 (p = .002), across -.35 (p < .0001), tolerance 10%")
reading wToM rs wToM p wPain rs wPain p across rs across p 0 O1 all, O2a, O3a 0.306 0.000639 0.461 1.05e-07 -0.259 0.00407 1 O1 all, O2a, O3b 0.306 0.000639 0.47 5.39e-08 -0.271 0.00263 2 O1 all, O2b, O3a 0.324 0.000294 0.458 1.24e-07 -0.261 0.00379 3 O1 all, O2b, O3b 0.324 0.000285 0.465 7.63e-08 -0.271 0.0026 4 O1 with_data, O2a, O3a 0.306 0.000639 0.461 1.05e-07 -0.259 0.00407 5 O1 with_data, O2a, O3b 0.306 0.000639 0.47 5.39e-08 -0.271 0.00263 6 O1 with_data, O2b, O3a 0.324 0.000294 0.458 1.24e-07 -0.261 0.00379 7 O1 with_data, O2b, O3b 0.324 0.000285 0.465 7.63e-08 -0.271 0.0026 df = 119 ; expected wToM .37 (p < .00005), wPain .28 (p = .002), across -.35 (p < .0001), tolerance 10%
Result: not reproduced (irc-age-within), not reproduced (irc-age-across), by Rule S. Within-ToM rises with age (partial rs = .306, p = .00064) and within-Pain too (.461, p = 1.1e-07), and the across-network correlation falls (-.259, p = .0041), so the directions hold and all three are significant at .05. But rs within-ToM and across are out of tolerance (.37 and -.35, 10%) and their p-values do not reach the paper's bounds (p < .00005 and p < .0001) under any reading; within-Pain is much stronger than the paper's .28. Rule S counts a test that does not meet the paper's bound as failing.
irc-adult-vs-child¶
Adults have higher within-ToM and within-Pain correlations and a more negative across-network correlation than children, controlling for motion.
Regression of each measure on group (child = 1, adult = 0) and number of artifact timepoints.
In [25]:
rows = []
for r in READINGS:
t = IRC[r][0].copy(); t["child"] = (t.Child_Adult == "child").astype(float)
for v in ["wToM", "wPain", "across"]:
f = sb.ols(t, v, ["child", "n_art"]); fz = sb.ols(t, v, ["child", "n_art"], zscore=True)
rows.append({"reading": reading_name(r), "measure": v, "b raw": f.params["child"], "b z (O4a)": fz.params["child"],
"t": f.tvalues["child"], "p": f.pvalues["child"], "t motion": f.tvalues["n_art"], "p motion": f.pvalues["n_art"]})
print(pd.DataFrame(rows).to_string(float_format=lambda x: f"{x:.3g}"))
print("expected t: wToM -5.7, wPain -3.8, across 7.2 (10%); paper motion effect: wToM t = -4.3, others n.s.")
reading measure b raw b z (O4a) t p t motion p motion 0 O1 all, O2a, O3a wToM -0.173 -1.16 -6.54 8.66e-10 -1.64 0.103 1 O1 all, O2a, O3a wPain -0.135 -1.2 -6.72 3.47e-10 -0.512 0.609 2 O1 all, O2a, O3a across 0.0868 0.83 4.49 1.38e-05 2.87 0.00475 3 O1 all, O2a, O3b wToM -0.167 -1.11 -6.23 4.38e-09 -1.77 0.0793 4 O1 all, O2a, O3b wPain -0.132 -1.17 -6.48 1.2e-09 -0.495 0.621 5 O1 all, O2a, O3b across 0.0836 0.798 4.28 3.28e-05 2.77 0.0064 6 O1 all, O2b, O3a wToM -0.18 -1.22 -7.05 5.94e-11 -1.89 0.0608 7 O1 all, O2b, O3a wPain -0.139 -1.21 -6.72 3.42e-10 0.103 0.918 8 O1 all, O2b, O3a across 0.0814 0.769 4.17 5.19e-05 3.3 0.00122 9 O1 all, O2b, O3b wToM -0.174 -1.17 -6.71 3.69e-10 -1.85 0.0665 10 O1 all, O2b, O3b wPain -0.137 -1.19 -6.56 7.96e-10 0.0573 0.954 11 O1 all, O2b, O3b across 0.0774 0.728 3.92 0.000135 3.35 0.00102 12 O1 with_data, O2a, O3a wToM -0.173 -1.16 -6.54 8.66e-10 -1.64 0.103 13 O1 with_data, O2a, O3a wPain -0.135 -1.2 -6.72 3.47e-10 -0.512 0.609 14 O1 with_data, O2a, O3a across 0.0868 0.83 4.49 1.38e-05 2.87 0.00475 15 O1 with_data, O2a, O3b wToM -0.167 -1.11 -6.23 4.38e-09 -1.77 0.0794 16 O1 with_data, O2a, O3b wPain -0.132 -1.17 -6.48 1.2e-09 -0.495 0.621 17 O1 with_data, O2a, O3b across 0.0836 0.798 4.28 3.28e-05 2.77 0.0064 18 O1 with_data, O2b, O3a wToM -0.18 -1.22 -7.05 5.93e-11 -1.89 0.0609 19 O1 with_data, O2b, O3a wPain -0.139 -1.21 -6.72 3.42e-10 0.103 0.918 20 O1 with_data, O2b, O3a across 0.0814 0.769 4.17 5.19e-05 3.3 0.00122 21 O1 with_data, O2b, O3b wToM -0.174 -1.17 -6.71 3.68e-10 -1.85 0.0666 22 O1 with_data, O2b, O3b wPain -0.137 -1.19 -6.56 7.96e-10 0.0572 0.954 23 O1 with_data, O2b, O3b across 0.0774 0.728 3.92 0.000135 3.35 0.00102 expected t: wToM -5.7, wPain -3.8, across 7.2 (10%); paper motion effect: wToM t = -4.3, others n.s.
Result: partially reproduced. Adults have higher within-ToM (t = -6.54 for the child effect, p = 8.7e-10) and within-Pain (t = -6.72, p = 3.5e-10) correlations and a more negative across-network correlation (t = 4.49, p = 1.4e-05) than children, so the tests hold with the paper's signs, but the t values are out of tolerance (10% of -5.7, -3.8 and 7.2): within-Pain (-6.48 to -6.72) and across (3.92 to 4.49) under every reading, and within-ToM under three of the four O2 and O3 readings (-6.54, -6.71, -7.05; -6.23 under O2a with O3b is inside). z-scaled b values are -1.16, -1.20 and .83 (paper -.97, -.75, 1.26). The motion effects differ from the paper's: on within-ToM t = -1.64 (paper -4.3), on across t = 2.87, p = .005 (paper not significant).
irc-control-networks and irc-williams¶
Without M1 regression, within-ToM (rs .39), within-Pain (.38), within-Scene (.23) and within-M1 (.30) correlations rise with age in children, within-Face and within-V1 do not, and the across ToM-Pain correlation falls with age (-.26) while across Face-Scene and M1-V1 do not change.
The face and scene parcels are those deposited in derivatives/ROIs (the same Julian et al. 2012 parcels the claims file names as an outside input).
In [26]:
CTRL = ["wToM", "wPain", "wScene", "wM1", "across", "wFace", "wV1", "aFaceScene", "aM1V1"]
EXP = {"wToM": .39, "wPain": .38, "wScene": .23, "wM1": .30, "across": -.26}
rows = []
for r in READINGS:
c = IRC_NOM1[r][0].loc[children.index]
for v in CTRL:
rr, pv, _ = sb.partial_spearman(c[v], c.Age, c.n_art)
rows.append({"reading": reading_name(r), "measure": v, "rs": rr, "p": pv,
"expected": EXP.get(v, np.nan), "in tolerance": tol(rr, EXP[v], abs(EXP[v]) * .1) if v in EXP else None})
ctrl = pd.DataFrame(rows)
print(ctrl.pivot(index="reading", columns="measure", values="rs")[CTRL].round(3).to_string())
print(ctrl.pivot(index="reading", columns="measure", values="p")[CTRL].map(lambda x: f"{x:.2g}").to_string())
print("paper: wToM .39 (<.0001), wPain .38 (<.0001), wScene .23 (.01), wM1 .30 (.001), across -.26 (.005); wFace .06 (.53), wV1 -.16 (.09), aFaceScene .02 (.8), aM1V1 .18 (.05)")
print("adult means (main reading):", IRC_NOM1[MAIN][0].loc[adults.index, CTRL].mean().round(3).to_dict())
measure wToM wPain wScene wM1 across wFace wV1 aFaceScene aM1V1
reading
O1 all, O2a, O3a 0.367 0.525 0.288 0.227 -0.247 0.106 -0.069 0.065 0.200
O1 all, O2a, O3b 0.369 0.547 0.294 0.242 -0.253 0.119 -0.077 0.071 0.217
O1 all, O2b, O3a 0.371 0.526 0.291 0.235 -0.251 0.118 -0.061 0.071 0.193
O1 all, O2b, O3b 0.379 0.538 0.300 0.246 -0.251 0.125 -0.083 0.076 0.216
O1 with_data, O2a, O3a 0.367 0.525 0.288 0.227 -0.247 0.106 -0.069 0.065 0.200
O1 with_data, O2a, O3b 0.369 0.547 0.294 0.242 -0.253 0.119 -0.077 0.071 0.217
O1 with_data, O2b, O3a 0.371 0.526 0.291 0.235 -0.251 0.118 -0.061 0.071 0.193
O1 with_data, O2b, O3b 0.379 0.538 0.300 0.246 -0.251 0.125 -0.083 0.076 0.216
measure wToM wPain wScene wM1 across wFace wV1 aFaceScene aM1V1
reading
O1 all, O2a, O3a 3.5e-05 6.4e-10 0.0014 0.012 0.0064 0.25 0.45 0.48 0.028
O1 all, O2a, O3b 3.2e-05 8.9e-11 0.001 0.0074 0.0051 0.2 0.4 0.44 0.017
O1 all, O2b, O3a 2.8e-05 5.7e-10 0.0012 0.0096 0.0056 0.2 0.51 0.44 0.034
O1 all, O2b, O3b 1.8e-05 2e-10 0.00083 0.0066 0.0054 0.17 0.37 0.41 0.017
O1 with_data, O2a, O3a 3.5e-05 6.4e-10 0.0014 0.012 0.0064 0.25 0.45 0.48 0.028
O1 with_data, O2a, O3b 3.2e-05 8.9e-11 0.001 0.0074 0.0051 0.2 0.4 0.44 0.017
O1 with_data, O2b, O3a 2.8e-05 5.7e-10 0.0012 0.0096 0.0056 0.2 0.51 0.44 0.034
O1 with_data, O2b, O3b 1.8e-05 2e-10 0.00083 0.0066 0.0054 0.17 0.37 0.41 0.017
paper: wToM .39 (<.0001), wPain .38 (<.0001), wScene .23 (.01), wM1 .30 (.001), across -.26 (.005); wFace .06 (.53), wV1 -.16 (.09), aFaceScene .02 (.8), aM1V1 .18 (.05)
adult means (main reading): {'wToM': 0.603, 'wPain': 0.399, 'wScene': 0.657, 'wM1': 0.736, 'across': -0.07, 'wFace': 0.675, 'wV1': 0.687, 'aFaceScene': 0.291, 'aM1V1': 0.048}
Result: partially reproduced. The ToM and Pain age effects hold (within-ToM rs = .367, p = 3.5e-05; within-Pain .525, p = 6.4e-10; across -.247, p = .0064), with within-ToM and across in tolerance and within-Pain out (.38 expected). Within-Scene (.288, p = .0014) and within-M1 (.227, p = .012) rise with age as in the paper, both out of tolerance. Within-Face (.106, p = .25), within-V1 (-.069, p = .45) and across Face-Scene (.065, p = .48) show no age effect, as in the paper, but across M1-V1 does (.200, p = .028; paper .18, p = .05) under every reading.
In [27]:
PAIRS = [("wToM", "wFace"), ("wPain", "wFace"), ("wToM", "wScene"), ("wPain", "wScene"), ("wToM", "wM1"), ("wPain", "wM1"),
("wToM", "wV1"), ("wPain", "wV1"), ("across", "aFaceScene"), ("across", "aM1V1")]
SHOULD = [True, True, False, False, False, False, True, True, True, True]
rows = []
for r in READINGS:
c = IRC_NOM1[r][0].loc[children.index]
for reading_w in ["a (partial, motion)", "b (no covariate)"]:
for (x, y), should in zip(PAIRS, SHOULD):
if reading_w.startswith("a"):
rx = sb.partial_spearman(c[x], c.Age, c.n_art)[0]; ry = sb.partial_spearman(c[y], c.Age, c.n_art)[0]
rxy = sb.partial_spearman(c[x], c[y], c.n_art)[0]
else:
rx = stats.spearmanr(c[x], c.Age).statistic; ry = stats.spearmanr(c[y], c.Age).statistic
rxy = stats.spearmanr(c[x], c[y]).statistic
w = sb.williams(rx, ry, rxy, len(c))
sig = w["p_t"] < .05
same_dir = (abs(rx) > abs(ry))
rows.append({"reading": reading_name(r), "williams": reading_w, "pair": f"{x} vs {y}", "r age x": rx, "r age y": ry,
"r xy": rxy, "t": w["t"], "p_t": w["p_t"], "z": w["z"], "p_z": w["p_z"], "paper significant": should,
"as paper": (sig and same_dir) if should else (not sig)})
wil = pd.DataFrame(rows)
print(wil[(wil.reading == reading_name(MAIN))].drop(columns="reading").to_string(float_format=lambda x: f"{x:.3g}"))
print(wil.groupby(["reading", "williams"])["as paper"].sum().rename("comparisons as in the paper, of 10").to_string())
williams pair r age x r age y r xy t p_t z p_z paper significant as paper
0 a (partial, motion) wToM vs wFace 0.367 0.106 0.265 2.51 0.0134 2.47 0.0137 True True
1 a (partial, motion) wPain vs wFace 0.525 0.106 0.167 4.08 8.11e-05 3.92 8.97e-05 True True
2 a (partial, motion) wToM vs wScene 0.367 0.288 0.222 0.752 0.454 0.75 0.453 False True
3 a (partial, motion) wPain vs wScene 0.525 0.288 0.342 2.62 0.00982 2.58 0.00996 False False
4 a (partial, motion) wToM vs wM1 0.367 0.227 0.222 1.32 0.191 1.31 0.19 False True
5 a (partial, motion) wPain vs wM1 0.525 0.227 0.167 2.92 0.00425 2.86 0.00428 False False
6 a (partial, motion) wToM vs wV1 0.367 -0.0686 -0.134 3.36 0.00104 3.27 0.00107 True True
7 a (partial, motion) wPain vs wV1 0.525 -0.0686 -0.092 5.03 1.73e-06 4.76 1.97e-06 True True
8 a (partial, motion) across vs aFaceScene -0.247 0.0654 0.277 -2.95 0.00384 -2.87 0.00409 True True
9 a (partial, motion) across vs aM1V1 -0.247 0.2 0.158 -4 0.000109 -3.82 0.000132 True True
10 b (no covariate) wToM vs wFace 0.366 0.103 0.269 2.53 0.0127 2.48 0.013 True True
11 b (no covariate) wPain vs wFace 0.524 0.103 0.154 4.07 8.59e-05 3.9 9.44e-05 True True
12 b (no covariate) wToM vs wScene 0.366 0.276 0.227 0.849 0.398 0.847 0.397 False True
13 b (no covariate) wPain vs wScene 0.524 0.276 0.311 2.68 0.00847 2.63 0.00859 False False
14 b (no covariate) wToM vs wM1 0.366 0.227 0.223 1.31 0.193 1.3 0.193 False True
15 b (no covariate) wPain vs wM1 0.524 0.227 0.164 2.91 0.00436 2.85 0.00439 False False
16 b (no covariate) wToM vs wV1 0.366 -0.0657 -0.139 3.32 0.0012 3.23 0.00124 True True
17 b (no covariate) wPain vs wV1 0.524 -0.0657 -0.0783 5.04 1.71e-06 4.76 1.97e-06 True True
18 b (no covariate) across vs aFaceScene -0.23 0.0665 0.286 -2.81 0.00587 -2.74 0.00618 True True
19 b (no covariate) across vs aM1V1 -0.23 0.195 0.101 -3.65 0.000396 -3.51 0.000445 True True
reading williams
O1 all, O2a, O3a a (partial, motion) 8
b (no covariate) 8
O1 all, O2a, O3b a (partial, motion) 8
b (no covariate) 8
O1 all, O2b, O3a a (partial, motion) 8
b (no covariate) 8
O1 all, O2b, O3b a (partial, motion) 8
b (no covariate) 8
O1 with_data, O2a, O3a a (partial, motion) 8
b (no covariate) 8
O1 with_data, O2a, O3b a (partial, motion) 8
b (no covariate) 8
O1 with_data, O2b, O3a a (partial, motion) 8
b (no covariate) 8
O1 with_data, O2b, O3b a (partial, motion) 8
b (no covariate) 8
Result: not reproduced. The six differences the paper finds significant are significant in the same direction under both readings (Face and V1 against ToM and Pain, t = 2.51 to 5.03; across ToM-Pain against Face-Scene and M1-V1, t = -2.95 and -4.0, main reading, reading a). Of the four the paper finds not significant, ToM against Scene (t = .75) and ToM against M1 (t = 1.32) are not, but Pain against Scene (t = 2.62, p = .0098) and Pain against M1 (t = 2.92, p = .0043) are, so the within-Pain age effect is stronger than the Scene and M1 effects, unlike the paper. 8 of 10 comparisons follow the paper under every reading. The criterion requires all ten and states no partial outcome, and two of the paper's tests (p = .19 and .47) come out on the other side of .05, which Rule S counts as failed tests.
irc-3yo-distinct, irc-3yo-not-anticorrelated, irc-3yo-homologous, irc-3yo-nonbilateral¶
The 17 children of the 3yo group.
In [28]:
y3 = children[children.AgeGroup == "3yo"].index
rows = []
for r in READINGS:
t = IRC[r][0].loc[y3]
t1 = stats.ttest_rel(t.wToM, t.across); t2 = stats.ttest_rel(t.wPain, t.across)
rows.append({"reading": reading_name(r), "wToM": t.wToM.mean(), "sem ToM": t.wToM.sem(), "wPain": t.wPain.mean(),
"sem Pain": t.wPain.sem(), "t ToM-across": t1.statistic, "p1": t1.pvalue, "t Pain-across": t2.statistic, "p2": t2.pvalue,
"in tolerance (of 4)": sum([tol(t.wToM.mean(), .21, .02), tol(t.wPain.mean(), .23, .02),
tol(t1.statistic, 6.2, .62), tol(t2.statistic, 6.9, .69)])})
print(pd.DataFrame(rows).to_string(float_format=lambda x: f"{x:.3g}"))
print("expected wToM .21, wPain .23 (tolerance .02), t 6.2 and 6.9 (10%), n = 17")
reading wToM sem ToM wPain sem Pain t ToM-across p1 t Pain-across p2 in tolerance (of 4) 0 O1 all, O2a, O3a 0.326 0.0273 0.179 0.0212 8.43 2.77e-07 4.54 0.000338 0 1 O1 all, O2a, O3b 0.325 0.0264 0.182 0.0218 8.71 1.82e-07 4.43 0.000418 0 2 O1 all, O2b, O3a 0.321 0.0247 0.185 0.0204 8.37 3.06e-07 5 0.00013 0 3 O1 all, O2b, O3b 0.322 0.0244 0.187 0.0212 8.52 2.42e-07 4.79 0.000201 0 4 O1 with_data, O2a, O3a 0.326 0.0273 0.179 0.0212 8.43 2.77e-07 4.54 0.000338 0 5 O1 with_data, O2a, O3b 0.325 0.0264 0.182 0.0218 8.71 1.82e-07 4.43 0.000418 0 6 O1 with_data, O2b, O3a 0.321 0.0247 0.185 0.0204 8.37 3.06e-07 5 0.00013 0 7 O1 with_data, O2b, O3b 0.322 0.0244 0.187 0.0212 8.52 2.42e-07 4.79 0.000201 0 expected wToM .21, wPain .23 (tolerance .02), t 6.2 and 6.9 (10%), n = 17
Result: not reproduced, by Rule S. In the 17 three-year-olds within-ToM (.326) and within-Pain (.179) exceed the across-network value (t = 8.43, p = 2.8e-07; t = 4.54, p = .00034), so the networks are distinct, but none of the four values is within tolerance (.21, .23, 6.2, 6.9), and the Pain test does not reach the paper's bound of p < .00001 under any reading.
Result: partially reproduced. The 3yo across-network mean is positive, .074 (s.e. .020), not anti-correlated, against adults' -.051; it is just outside the tolerance (.05 ± .02) under every reading (.074 to .075). A one-sample t-test shows it above 0 (t = 3.75, p = .0017).
In [30]:
for r in READINGS:
M = pd.DataFrame(np.mean([IRC[r][1][s].values for s in y3], axis=0), index=TP, columns=TP)
print(reading_name(r))
for net, rois, homol in [("ToM", sb.TOM_ROIS, ("RTPJ", "LTPJ")), ("Pain", sb.PAIN_ROIS, ("RInsula", "LInsula"))]:
pairs = sorted([(M.loc[a, b], a, b) for i, a in enumerate(rois) for b in rois[i + 1:]], reverse=True)
rank = [i for i, (v, a, b) in enumerate(pairs) if {a, b} == set(homol)][0] + 1
print(f" {net} pairs by 3yo mean z:", ", ".join(f"{a}-{b} {v:.3f}" for v, a, b in pairs[:5]),
f"| {homol[0]}-{homol[1]} ranks {rank} of {len(pairs)}")
M = pd.DataFrame(np.mean([IRC[MAIN][1][s].values for s in y3], axis=0), index=TP, columns=TP)
fig, ax = plt.subplots(figsize=(5, 4.5))
im = ax.imshow(M.values, vmin=-0.6, vmax=0.6, cmap="RdBu_r"); ax.set_xticks(range(13), TP, rotation=90); ax.set_yticks(range(13), TP)
plt.colorbar(im, label="mean Fisher z, 3yo"); plt.tight_layout(); plt.show()
print("matrix range (off-diagonal):", round(np.nanmin(M.values), 3), "to", round(np.nanmax(M.values), 3), "; colour scale -0.6 to 0.6")
O1 all, O2a, O3a ToM pairs by 3yo mean z: RTPJ-PC 0.637, LTPJ-PC 0.635, DMPFC-MMPFC 0.609, RTPJ-LTPJ 0.558, MMPFC-VMPFC 0.466 | RTPJ-LTPJ ranks 4 of 15 Pain pairs by 3yo mean z: RInsula-LInsula 0.856, RS2-LS2 0.656, RMFG-LMFG 0.339, LInsula-AMCC 0.313, RInsula-AMCC 0.267 | RInsula-LInsula ranks 1 of 21 O1 all, O2a, O3b ToM pairs by 3yo mean z: RTPJ-PC 0.634, LTPJ-PC 0.630, DMPFC-MMPFC 0.609, RTPJ-LTPJ 0.555, MMPFC-VMPFC 0.473 | RTPJ-LTPJ ranks 4 of 15 Pain pairs by 3yo mean z: RInsula-LInsula 0.855, RS2-LS2 0.662, RMFG-LMFG 0.345, LInsula-AMCC 0.303, RInsula-AMCC 0.260 | RInsula-LInsula ranks 1 of 21 O1 all, O2b, O3a ToM pairs by 3yo mean z: LTPJ-PC 0.629, RTPJ-PC 0.627, DMPFC-MMPFC 0.596, RTPJ-LTPJ 0.553, MMPFC-VMPFC 0.446 | RTPJ-LTPJ ranks 4 of 15 Pain pairs by 3yo mean z: RInsula-LInsula 0.861, RS2-LS2 0.659, RMFG-LMFG 0.333, LInsula-AMCC 0.309, RInsula-AMCC 0.259 | RInsula-LInsula ranks 1 of 21 O1 all, O2b, O3b ToM pairs by 3yo mean z: LTPJ-PC 0.626, RTPJ-PC 0.625, DMPFC-MMPFC 0.595, RTPJ-LTPJ 0.551, MMPFC-VMPFC 0.456 | RTPJ-LTPJ ranks 4 of 15 Pain pairs by 3yo mean z: RInsula-LInsula 0.858, RS2-LS2 0.666, RMFG-LMFG 0.337, LInsula-AMCC 0.303, RInsula-AMCC 0.251 | RInsula-LInsula ranks 1 of 21 O1 with_data, O2a, O3a ToM pairs by 3yo mean z: RTPJ-PC 0.637, LTPJ-PC 0.635, DMPFC-MMPFC 0.609, RTPJ-LTPJ 0.558, MMPFC-VMPFC 0.466 | RTPJ-LTPJ ranks 4 of 15 Pain pairs by 3yo mean z: RInsula-LInsula 0.856, RS2-LS2 0.656, RMFG-LMFG 0.339, LInsula-AMCC 0.313, RInsula-AMCC 0.267 | RInsula-LInsula ranks 1 of 21 O1 with_data, O2a, O3b ToM pairs by 3yo mean z: RTPJ-PC 0.634, LTPJ-PC 0.630, DMPFC-MMPFC 0.609, RTPJ-LTPJ 0.555, MMPFC-VMPFC 0.473 | RTPJ-LTPJ ranks 4 of 15 Pain pairs by 3yo mean z: RInsula-LInsula 0.855, RS2-LS2 0.662, RMFG-LMFG 0.345, LInsula-AMCC 0.303, RInsula-AMCC 0.260 | RInsula-LInsula ranks 1 of 21 O1 with_data, O2b, O3a ToM pairs by 3yo mean z: LTPJ-PC 0.629, RTPJ-PC 0.627, DMPFC-MMPFC 0.596, RTPJ-LTPJ 0.553, MMPFC-VMPFC 0.446 | RTPJ-LTPJ ranks 4 of 15 Pain pairs by 3yo mean z: RInsula-LInsula 0.861, RS2-LS2 0.659, RMFG-LMFG 0.333, LInsula-AMCC 0.309, RInsula-AMCC 0.259 | RInsula-LInsula ranks 1 of 21 O1 with_data, O2b, O3b ToM pairs by 3yo mean z: LTPJ-PC 0.626, RTPJ-PC 0.625, DMPFC-MMPFC 0.595, RTPJ-LTPJ 0.551, MMPFC-VMPFC 0.456 | RTPJ-LTPJ ranks 4 of 15 Pain pairs by 3yo mean z: RInsula-LInsula 0.858, RS2-LS2 0.666, RMFG-LMFG 0.337, LInsula-AMCC 0.303, RInsula-AMCC 0.251 | RInsula-LInsula ranks 1 of 21
matrix range (off-diagonal): -0.103 to 0.856 ; colour scale -0.6 to 0.6
Result: not reproduced. In the 3yo mean matrix the strongest Pain pair is RInsula-LInsula (.856), and the three homologous Pain pairs are the top three, but the strongest ToM pair is RTPJ-PC (.637); RTPJ-LTPJ (.558) ranks fourth of 15. Neither reading of the criterion holds for ToM.
In [31]:
rows = []
for r in READINGS:
t = IRC[r][0].loc[y3]
for k in ["a", "b"]:
t1 = stats.ttest_rel(t[f"wToM_nonbi_{k}"], t.across); t2 = stats.ttest_rel(t[f"wPain_nonbi_{k}"], t.across)
rows.append({"reading": reading_name(r), "pairs dropped": k, "ToM": t[f"wToM_nonbi_{k}"].mean(), "Pain": t[f"wPain_nonbi_{k}"].mean(),
"t ToM": t1.statistic, "p ToM": t1.pvalue, "t Pain": t2.statistic, "p Pain": t2.pvalue,
"in tolerance (of 4)": sum([tol(t[f"wToM_nonbi_{k}"].mean(), .20, .02), tol(t[f"wPain_nonbi_{k}"].mean(), .17, .02),
tol(t1.statistic, 5.1, .51), tol(t2.statistic, 4.4, .44)])})
print(pd.DataFrame(rows).to_string(float_format=lambda x: f"{x:.3g}"))
print("expected ToM .20, Pain .17 (tolerance .02); t 5.1 (p = .0001), 4.4 (p = .0005)")
reading pairs dropped ToM Pain t ToM p ToM t Pain p Pain in tolerance (of 4) 0 O1 all, O2a, O3a a 0.309 0.106 7.64 1.01e-06 1.34 0.198 0 1 O1 all, O2a, O3a b 0.309 0.146 7.64 1.01e-06 3.03 0.00802 0 2 O1 all, O2a, O3b a 0.309 0.109 7.84 7.21e-07 1.34 0.2 0 3 O1 all, O2a, O3b b 0.309 0.149 7.84 7.21e-07 2.99 0.00865 0 4 O1 all, O2b, O3a a 0.304 0.113 7.64 1.01e-06 1.7 0.108 0 5 O1 all, O2b, O3a b 0.304 0.152 7.64 1.01e-06 3.43 0.00346 1 6 O1 all, O2b, O3b a 0.305 0.115 7.73 8.56e-07 1.62 0.125 0 7 O1 all, O2b, O3b b 0.305 0.153 7.73 8.56e-07 3.3 0.00453 1 8 O1 with_data, O2a, O3a a 0.309 0.106 7.64 1.01e-06 1.34 0.198 0 9 O1 with_data, O2a, O3a b 0.309 0.146 7.64 1.01e-06 3.03 0.00802 0 10 O1 with_data, O2a, O3b a 0.309 0.109 7.84 7.21e-07 1.34 0.2 0 11 O1 with_data, O2a, O3b b 0.309 0.149 7.84 7.21e-07 2.99 0.00865 0 12 O1 with_data, O2b, O3a a 0.304 0.113 7.64 1.01e-06 1.7 0.108 0 13 O1 with_data, O2b, O3a b 0.304 0.152 7.64 1.01e-06 3.43 0.00346 1 14 O1 with_data, O2b, O3b a 0.305 0.115 7.73 8.56e-07 1.62 0.125 0 15 O1 with_data, O2b, O3b b 0.305 0.153 7.73 8.56e-07 3.3 0.00453 1 expected ToM .20, Pain .17 (tolerance .02); t 5.1 (p = .0001), 4.4 (p = .0005)
Result: partially reproduced. Under reading (b), dropping only the two pairs the paper names, non-bilateral within-ToM (.309, t = 7.64, p = 1e-06) and within-Pain (.146, t = 3.03, p = .008) exceed the across-network value, so the tests hold, but the means and t values are out of tolerance (except the Pain mean under O2b, .152 and .153). Under reading (a), dropping all three Pain pairs, non-bilateral Pain (.106) does not differ from across (t = 1.34, p = .20).
irc-tom-behavior and irc-tom-behavior-age¶
In [32]:
rows = []
for r in READINGS:
c = IRC[r][0].loc[children.index]
a = sb.partial_kendall(c.wToM, c[TOM], c.n_art); b = sb.partial_kendall(c.across, c[TOM], c.n_art)
rows.append({"reading": reading_name(r), "tau wToM": a[0], "p wToM": a[1], "tau across": b[0], "p across": b[1]})
print(pd.DataFrame(rows).to_string(float_format=lambda x: f"{x:.3g}"))
print("expected wToM .23 (p = .0002), across -.20 (p = .001), tolerance 10%")
reading tau wToM p wToM tau across p across 0 O1 all, O2a, O3a 0.246 6.29e-05 -0.125 0.0416 1 O1 all, O2a, O3b 0.237 0.000113 -0.132 0.0313 2 O1 all, O2b, O3a 0.259 2.61e-05 -0.121 0.0498 3 O1 all, O2b, O3b 0.252 4.02e-05 -0.13 0.0347 4 O1 with_data, O2a, O3a 0.246 6.29e-05 -0.125 0.0416 5 O1 with_data, O2a, O3b 0.237 0.000113 -0.132 0.0313 6 O1 with_data, O2b, O3a 0.259 2.61e-05 -0.121 0.0498 7 O1 with_data, O2b, O3b 0.252 4.02e-05 -0.13 0.0347 expected wToM .23 (p = .0002), across -.20 (p = .001), tolerance 10%
Result: partially reproduced. Within-ToM correlation rises with ToM score controlling for motion (partial tau = .246, p = 6.3e-05; paper .23, within tolerance). The across-network correlation falls with ToM score (tau = -.125, p = .042) but less than the paper's -.20, out of tolerance; both tests are significant with the paper's signs.
In [33]:
rows = []
for r in READINGS:
c = IRC[r][0].loc[children.index]
for v in ["wToM", "across"]:
f = sb.ols(c, v, [TOM, "Age", "n_art"]); fz = sb.ols(c, v, [TOM, "Age", "n_art"], zscore=True)
rows.append({"reading": reading_name(r), "measure": v, "b ToM raw": f.params[TOM], "b ToM z": fz.params[TOM],
"t ToM": f.tvalues[TOM], "p ToM": f.pvalues[TOM], "t age": f.tvalues["Age"], "p age": f.pvalues["Age"]})
print(pd.DataFrame(rows).to_string(float_format=lambda x: f"{x:.3g}"))
print("paper: ToM score ts < 1, ps > .3")
reading measure b ToM raw b ToM z t ToM p ToM t age p age 0 O1 all, O2a, O3a wToM 0.181 0.275 2.04 0.0439 0.621 0.536 1 O1 all, O2a, O3a across 0.0761 0.173 1.29 0.2 -2.98 0.00353 2 O1 all, O2a, O3b wToM 0.167 0.247 1.83 0.0693 0.846 0.399 3 O1 all, O2a, O3b across 0.066 0.149 1.11 0.269 -2.94 0.00393 4 O1 all, O2b, O3a wToM 0.169 0.267 2 0.0477 0.826 0.41 5 O1 all, O2b, O3a across 0.0795 0.178 1.33 0.185 -2.95 0.00383 6 O1 all, O2b, O3b wToM 0.149 0.23 1.72 0.088 1.12 0.266 7 O1 all, O2b, O3b across 0.0711 0.158 1.19 0.235 -2.97 0.00362 8 O1 with_data, O2a, O3a wToM 0.181 0.275 2.04 0.0439 0.621 0.536 9 O1 with_data, O2a, O3a across 0.0761 0.173 1.29 0.2 -2.98 0.00353 10 O1 with_data, O2a, O3b wToM 0.167 0.247 1.83 0.0693 0.846 0.399 11 O1 with_data, O2a, O3b across 0.066 0.149 1.11 0.269 -2.94 0.00393 12 O1 with_data, O2b, O3a wToM 0.169 0.267 2 0.0477 0.827 0.41 13 O1 with_data, O2b, O3a across 0.0795 0.178 1.33 0.185 -2.95 0.00383 14 O1 with_data, O2b, O3b wToM 0.149 0.23 1.72 0.088 1.12 0.265 15 O1 with_data, O2b, O3b across 0.0711 0.158 1.19 0.235 -2.97 0.00362 paper: ToM score ts < 1, ps > .3
Result: not reproduced, resting on an assumption. Under the main reading (O2a, O3a) the ToM score effect on within-ToM is significant with age and motion in the model (t = 2.04, p = .044); on across-network it is not (t = 1.29, p = .20). The readings disagree: O3b gives p = .069, O2b .048, O2b with O3b .088. The main reading is adopted (see deviations).
irc-fb-within, irc-fb-within-dccs, irc-fb-across¶
Passers (30) against failers (15), fail coded 1, with age and motion (and DCCS).
In [34]:
pf_idx = young[young.FB_Group.isin(["pass", "fail"])].index
rows = []
for r in READINGS:
t = IRC[r][0].loc[pf_idx].copy(); t["fail"] = (t.FB_Group == "fail").astype(float)
for v in ["wToM", "across"]:
f = sb.ols(t, v, ["fail", "Age", "n_art"]); fz = sb.ols(t, v, ["fail", "Age", "n_art"], zscore=True)
g = sb.ols(t, v, ["fail", "Age", "n_art", DCCS])
rows.append({"reading": reading_name(r), "measure": v, "mean pass": t.loc[t.fail == 0, v].mean(), "mean fail": t.loc[t.fail == 1, v].mean(),
"b raw": f.params["fail"], "b z": fz.params["fail"], "t": f.tvalues["fail"], "p": f.pvalues["fail"],
"t with DCCS": g.tvalues["fail"], "p with DCCS": g.pvalues["fail"], "t DCCS": g.tvalues[DCCS], "p DCCS": g.pvalues[DCCS]})
fbirc = pd.DataFrame(rows)
print(fbirc.to_string(float_format=lambda x: f"{x:.3g}"))
print("paper: wToM pass .29, fail .25, b = -.70, t = -2.06, p = .046 (pass higher); with DCCS t = -1.80, p = .079; across .04 against .03, t = 1.2, p = .23")
reading measure mean pass mean fail b raw b z t p t with DCCS p with DCCS t DCCS p DCCS 0 O1 all, O2a, O3a wToM 0.398 0.345 0.0345 0.259 0.661 0.512 0.398 0.693 -0.863 0.393 1 O1 all, O2a, O3a across 0.0872 0.0565 -0.0679 -0.87 -2.18 0.035 -1.73 0.0905 1.42 0.164 2 O1 all, O2a, O3b wToM 0.387 0.338 0.0319 0.24 0.603 0.55 0.402 0.69 -0.635 0.529 3 O1 all, O2a, O3b across 0.0794 0.0565 -0.0643 -0.809 -2.01 0.0511 -1.54 0.131 1.55 0.128 4 O1 all, O2b, O3a wToM 0.387 0.339 0.0341 0.277 0.706 0.484 0.513 0.611 -0.585 0.562 5 O1 all, O2b, O3a across 0.0876 0.056 -0.0667 -0.861 -2.19 0.0342 -1.71 0.0955 1.61 0.115 6 O1 all, O2b, O3b wToM 0.377 0.335 0.0333 0.268 0.668 0.508 0.541 0.592 -0.347 0.73 7 O1 all, O2b, O3b across 0.0819 0.056 -0.065 -0.812 -2.05 0.047 -1.53 0.133 1.78 0.0824 8 O1 with_data, O2a, O3a wToM 0.398 0.345 0.0345 0.259 0.661 0.512 0.398 0.693 -0.863 0.393 9 O1 with_data, O2a, O3a across 0.0872 0.0565 -0.0679 -0.87 -2.18 0.035 -1.73 0.0905 1.42 0.164 10 O1 with_data, O2a, O3b wToM 0.387 0.338 0.0319 0.24 0.603 0.55 0.402 0.69 -0.635 0.529 11 O1 with_data, O2a, O3b across 0.0794 0.0565 -0.0643 -0.809 -2.01 0.0511 -1.54 0.131 1.55 0.128 12 O1 with_data, O2b, O3a wToM 0.387 0.339 0.0341 0.277 0.706 0.484 0.513 0.611 -0.585 0.562 13 O1 with_data, O2b, O3a across 0.0876 0.056 -0.0667 -0.861 -2.19 0.0342 -1.71 0.0955 1.61 0.115 14 O1 with_data, O2b, O3b wToM 0.377 0.335 0.0333 0.268 0.668 0.508 0.541 0.592 -0.347 0.73 15 O1 with_data, O2b, O3b across 0.0819 0.056 -0.065 -0.812 -2.05 0.047 -1.53 0.133 1.78 0.0824 paper: wToM pass .29, fail .25, b = -.70, t = -2.06, p = .046 (pass higher); with DCCS t = -1.80, p = .079; across .04 against .03, t = 1.2, p = .23
Result: not reproduced (irc-fb-within, irc-fb-within-dccs, irc-fb-across; irc-fb-across rests on an assumption). Passers have a slightly higher within-ToM correlation than failers (.398 against .345), but with age and motion the group effect is not significant (t = .66, p = .51; paper p = .046) under any reading, and with DCCS added p = .69, outside the .05 to .10 band. Across-network, passers (.087) are above failers (.057) and the group effect is significant under the main reading (t = -2.18, p = .035), where the paper finds none; O3b gives p = .051, O2b .034, O2b with O3b .047.
Reverse correlation¶
Network timecourses (P6) from ROI timecourses cleaned without the M1 regressor. O6: z over TRs 11 to 168 (a) or all 168 (b). O5: p < .05 (a), p < .01 (b), FDR p < .05 over the 158 timepoints (c). O1 is fixed at the main reading (z-normalization makes O1's voxel choice irrelevant up to the few partial voxels; shown in the IRC section).
In [35]:
def network_tcs(net, z=True, z_span="movie", voxels="all"):
rois = sb.FFA_ROIS if net == "FFA" else sb.NETWORKS[net]
out = {}
for s in subs:
tc = sb.subject_timecourses(D[s], rois, voxels, False)
out[s] = sb.network_timecourse(tc, rois, z_span, z)
return pd.DataFrame(out).T
NETS = ["ToM", "Pain", "Face", "Scene", "M1", "V1", "FFA"]
NZ = {(n, span): network_tcs(n, True, span) for n in ["ToM", "Pain", "FFA"] for span in ["movie", "all"]}
RAW = {n: network_tcs(n, False) for n in NETS}
O5 = {"a p<.05": (0.05, False), "b p<.01": (0.01, False), "c FDR<.05": (0.05, True)}
RC = {}
for span in ["movie", "all"]:
for k, (alpha, fdr) in O5.items():
for n in ["ToM", "Pain", "FFA"]:
tab, ev = sb.reverse_correlation(NZ[(n, span)].loc[adults.index].values, alpha, fdr)
RC[(n, span, k)] = (tab, ev, sb.rank_events(ev, tab))
rows = []
for (n, span, k), (tab, ev, rk) in RC.items():
if n == "FFA":
continue
m = sb.match_table_events(ev, n)
rows.append({"network": n, "O6 z span": span, "O5": k, "events": len(ev), "seconds": 2 * sum(len(e) for e in ev),
"table events found": int(m.sum()), "of": len(m), "missed": ",".join(m.index[~m])})
rcsum = pd.DataFrame(rows)
print(rcsum.to_string())
network O6 z span O5 events seconds table events found of missed 0 ToM movie a p<.05 9 80 7 7 1 Pain movie a p<.05 10 84 10 12 P08,P10 2 ToM movie b p<.01 7 56 7 7 3 Pain movie b p<.01 8 62 8 12 P08,P09,P10,P12 4 ToM movie c FDR<.05 7 56 7 7 5 Pain movie c FDR<.05 8 62 8 12 P08,P09,P10,P12 6 ToM all a p<.05 9 80 7 7 7 Pain all a p<.05 10 84 10 12 P08,P10 8 ToM all b p<.01 7 56 7 7 9 Pain all b p<.01 8 62 8 12 P08,P09,P10,P12 10 ToM all c FDR<.05 7 56 7 7 11 Pain all c FDR<.05 8 62 8 12 P08,P09,P10,P12
In [36]:
# Timing check (2): adult network contrast during the paper's ToM against pain events, with the
# event labels shifted in time.
tom_trs = np.array(sorted(tr_set_ := set(t for v in sb.SUPP_TABLE3[sb.SUPP_TABLE3.network == "ToM"].trs for t in v)))
pain_trs = np.array(sorted(set(t for v in sb.SUPP_TABLE3[sb.SUPP_TABLE3.network == "Pain"].trs for t in v)))
for shift in [-10, -4, -2, 0, 2, 4, 10]:
out = []
for n in ["ToM", "Pain"]:
m = np.nanmean(NZ[(n, "movie")].loc[adults.index].values, axis=0)
tt, pp = np.clip(tom_trs + shift, 11, 168) - 1, np.clip(pain_trs + shift, 11, 168) - 1
out.append(m[tt].mean() - m[pp].mean())
print(f"shift {shift:3d} TRs: ToM network, ToM minus pain events = {out[0]:.2f}; Pain network, ToM minus pain events = {out[1]:.2f}")
shift -10 TRs: ToM network, ToM minus pain events = -0.14; Pain network, ToM minus pain events = 0.20 shift -4 TRs: ToM network, ToM minus pain events = 0.11; Pain network, ToM minus pain events = 0.07 shift -2 TRs: ToM network, ToM minus pain events = 0.38; Pain network, ToM minus pain events = -0.32 shift 0 TRs: ToM network, ToM minus pain events = 0.71; Pain network, ToM minus pain events = -0.70 shift 2 TRs: ToM network, ToM minus pain events = 0.44; Pain network, ToM minus pain events = -0.36 shift 4 TRs: ToM network, ToM minus pain events = 0.11; Pain network, ToM minus pain events = -0.17 shift 10 TRs: ToM network, ToM minus pain events = -0.09; Pain network, ToM minus pain events = 0.01
rc-events¶
Reverse correlation on the adult network timecourses finds seven ToM events (68 s, mean length 9.7 s) and twelve pain events (86 s, mean length 7.2 s), at the times of Supplementary Table 3.
In [37]:
for (n, span, k), (tab, ev, rk) in RC.items():
if n != "FFA" and span == "movie":
print(n, k, "events (onset TR, TRs, peak TR):", [(int(a), int(b), int(c)) for a, b, c in zip(rk.onset, rk.length_tr, rk.peak_tr)])
print(" mean length %.1f s (s.d. %.1f)" % (2 * rk.length_tr.mean(), 2 * rk.length_tr.std()) if len(rk) else "")
print(sb.SUPP_TABLE3[["event", "time", "duration_s", "peak_tr", "onset_tr"]].to_string(index=False))
score = rcsum.groupby(["O6 z span", "O5"]).apply(lambda g: pd.Series({
"table events found": g["table events found"].sum(), "count error": (g.events - g["of"]).abs().sum(),
"ToM events": int(g.loc[g.network == "ToM", "events"].iloc[0]), "Pain events": int(g.loc[g.network == "Pain", "events"].iloc[0]),
"ToM s": int(g.loc[g.network == "ToM", "seconds"].iloc[0]), "Pain s": int(g.loc[g.network == "Pain", "seconds"].iloc[0])}), include_groups=False)
print(score.to_string())
# Most table events found, then smallest count error; ties between the two O6 spans go to O6b
# (all 168 TRs): the Methods (paragraph 9) z-normalize each ROI timecourse and only then leave out
# the first 10 timepoints from the t-tests.
score["O6 preference"] = [0 if span == "all" else 1 for span, _ in score.index]
best = score.sort_values(["table events found", "count error", "O6 preference"], ascending=[False, True, True]).index[0]
BEST_SPAN, BEST_O5 = best
print("best-matching reading (most table events found, then smallest count error, then O6b):", best)
ToM a p<.05 events (onset TR, TRs, peak TR): [(129, 4, 131), (90, 9, 95), (63, 4, 64), (149, 7, 150), (45, 3, 46), (122, 5, 125), (52, 4, 54), (163, 2, 163), (24, 2, 25)]
mean length 8.9 s (s.d. 4.6)
Pain a p<.05 events (onset TR, TRs, peak TR): [(107, 7, 112), (116, 5, 118), (71, 8, 77), (38, 4, 39), (158, 4, 160), (134, 3, 135), (83, 2, 84), (56, 5, 56), (142, 2, 142), (100, 2, 100)]
mean length 8.4 s (s.d. 4.2)
ToM b p<.01 events (onset TR, TRs, peak TR): [(130, 3, 131), (91, 7, 95), (63, 4, 64), (150, 6, 150), (46, 2, 46), (122, 4, 125), (53, 2, 54)]
mean length 8.0 s (s.d. 3.8)
Pain b p<.01 events (onset TR, TRs, peak TR): [(107, 7, 112), (116, 5, 118), (71, 7, 77), (38, 3, 39), (158, 3, 160), (134, 2, 135), (83, 2, 84), (56, 2, 56)]
mean length 7.8 s (s.d. 4.3)
ToM c FDR<.05 events (onset TR, TRs, peak TR): [(130, 3, 131), (91, 7, 95), (63, 4, 64), (150, 6, 150), (46, 2, 46), (122, 4, 125), (53, 2, 54)]
mean length 8.0 s (s.d. 3.8)
Pain c FDR<.05 events (onset TR, TRs, peak TR): [(107, 7, 112), (116, 5, 118), (71, 7, 77), (38, 3, 39), (158, 3, 160), (134, 2, 135), (83, 2, 84), (56, 2, 56)]
mean length 7.8 s (s.d. 4.3)
event time duration_s peak_tr onset_tr
T01 4:00 10 131 128
T02 2:46 16 96 91
T03 1:14 6 46 45
T04 4:42 14 150 149
T05 1:28 8 54 52
T06 3:48 10 124 122
T07 1:50 4 64 63
P01 3:36 10 118 116
P02 2:06 18 73 71
P03 3:20 12 111 108
P04 1:00 6 39 38
P05 1:34 6 56 55
P06 4:10 10 135 133
P07 5:02 4 160 159
P08 0:42 4 29 29
P09 4:28 4 143 142
P10 1:20 4 49 48
P11 2:30 4 84 83
P12 3:04 4 100 100
table events found count error ToM events Pain events ToM s Pain s
O6 z span O5
all a p<.05 17 4 9 10 80 84
b p<.01 15 4 7 8 56 62
c FDR<.05 15 4 7 8 56 62
movie a p<.05 17 4 9 10 80 84
b p<.01 15 4 7 8 56 62
c FDR<.05 15 4 7 8 56 62
best-matching reading (most table events found, then smallest count error, then O6b): ('all', 'a p<.05')
In [38]:
fig, axes = plt.subplots(2, 1, figsize=(12, 5), sharex=True)
for ax, n, col in zip(axes, ["ToM", "Pain"], ["tab:red", "tab:green"]):
tab, ev, rk = RC[(n, BEST_SPAN, BEST_O5)]
ax.plot(tab.TR, tab["mean"], color=col)
for e in ev:
ax.axvspan(e[0] - 0.5, e[-1] + 0.5, color=col, alpha=0.2)
for _, row in sb.SUPP_TABLE3[sb.SUPP_TABLE3.network == n].iterrows():
ax.axvspan(row.trs[0] - 0.5, row.trs[-1] + 0.5, fill=False, edgecolor="k", lw=0.8)
ax.set_ylabel(f"adult mean, {n} (z)"); ax.axhline(0, color="grey", lw=0.5)
axes[1].set_xlabel("TR (shaded: found events; outlined: Supplementary Table 3)")
plt.tight_layout(); plt.show()
print("y ranges:", {n: (round(RC[(n, BEST_SPAN, BEST_O5)][0]["mean"].min(), 2), round(RC[(n, BEST_SPAN, BEST_O5)][0]["mean"].max(), 2)) for n in ["ToM", "Pain"]})
y ranges: {'ToM': (np.float64(-0.83), np.float64(0.79)), 'Pain': (np.float64(-0.88), np.float64(0.81))}
Result: partially reproduced. At p < .05 uncorrected (reading O5a, the best match to Supplementary Table 3; O6 makes no difference for these networks), all seven ToM events and 10 of 12 pain events of Supplementary Table 3 are found (P08 and P10 missed), but the counts are 9 ToM events (80 s) and 10 pain events (84 s), against 7 (68 s) and 12 (86 s). The ToM count matches under O5b and O5c (7 events, 56 s), where only 8 of 12 pain events are found. The pain time (84 s) is within tolerance; the ToM time is not under any reading.
rc-overlap¶
ToM and pain events overlap at a single timepoint (the last of T05 and first of P05), less overlap than chance in a permutation test.
Permutations as Methods paragraph 10, 1000 per reading, seed fixed for this section.
In [39]:
rng = np.random.default_rng(20180301)
# Reading (b): 159-TR timecourses, the length implied by the paper's zeros (125 and 116) and
# event totals (34 and 43 TRs); with the Supplementary Table 3 events it gives 125 and 116 zeros.
def overlap_test(ev_a, ev_b, zeros_a, zeros_b, n_perm=1000):
obs = len(set(t for e in ev_a for t in e) & set(t for e in ev_b for t in e))
perm = sb.overlap_permutation([len(e) for e in ev_a], [len(e) for e in ev_b], zeros_a, zeros_b, n_perm, rng)
return obs, (perm <= obs).sum(), perm.mean()
rows = []
evsets = {"found (" + BEST_SPAN + ", " + BEST_O5 + ")": (RC[("ToM", BEST_SPAN, BEST_O5)][1], RC[("Pain", BEST_SPAN, BEST_O5)][1])}
for k in O5:
evsets[f"found (movie, {k})"] = (RC[("ToM", "movie", k)][1], RC[("Pain", "movie", k)][1])
evsets["Supplementary Table 3 events (variant)"] = (list(sb.SUPP_TABLE3[sb.SUPP_TABLE3.network == "ToM"].trs), list(sb.SUPP_TABLE3[sb.SUPP_TABLE3.network == "Pain"].trs))
for name, (et, ep) in evsets.items():
nt, np_ = sum(map(len, et)), sum(map(len, ep))
for reading, (za, zb) in {"(a) 158 TRs": (158 - nt, 158 - np_), "(b) 159 TRs": (159 - nt, 159 - np_)}.items():
obs, k, mean = overlap_test(et, ep, za, zb)
rows.append({"events": name, "reading": reading, "observed overlap TRs": obs, "perms <= observed": k, "p": k / 1000, "mean perm overlap": mean})
ovl = pd.DataFrame(rows)
print(ovl.to_string())
both = set(t for e in RC[("ToM", BEST_SPAN, BEST_O5)][1] for t in e) & set(t for e in RC[("Pain", BEST_SPAN, BEST_O5)][1] for t in e)
print("overlapping TRs of found events:", sorted(both))
events reading observed overlap TRs perms <= observed p mean perm overlap 0 found (all, a p<.05) (a) 158 TRs 0 2 0.002 11.914 1 found (all, a p<.05) (b) 159 TRs 0 0 0.000 11.945 2 found (movie, a p<.05) (a) 158 TRs 0 3 0.003 11.601 3 found (movie, a p<.05) (b) 159 TRs 0 0 0.000 11.854 4 found (movie, b p<.01) (a) 158 TRs 0 19 0.019 6.470 5 found (movie, b p<.01) (b) 159 TRs 0 45 0.045 6.403 6 found (movie, c FDR<.05) (a) 158 TRs 0 36 0.036 6.578 7 found (movie, c FDR<.05) (b) 159 TRs 0 49 0.049 6.360 8 Supplementary Table 3 events (variant) (a) 158 TRs 1 2 0.002 9.515 9 Supplementary Table 3 events (variant) (b) 159 TRs 1 8 0.008 9.692 overlapping TRs of found events: []
Result: partially reproduced. The ToM and pain events found (O5a) do not overlap at all (0 TRs, where the paper has 1), which is less than chance: 2 of 1000 permutations have as little overlap under reading (a), 158 TRs (p = .002, mean permuted overlap 11.9 TRs), and none under (b), 159 TRs (p < .001, mean 11.9). Under O5b and O5c the overlap is also 0, p = .019 to .049. The overlap count differs from the paper's by one timepoint. The criterion states only when the claim is reproduced. The outcome is partial because the claim's test (less overlap than chance, p <= .01) holds and only its expected count is out of tolerance, which is what the claims file's Rule S and the protocol call a partial result; under irc-williams, by contrast, two of the paper's tests fail. With the paper's own events from Supplementary Table 3 the overlap is 1 TR (the last of T05 and first of P05), p = .002 (a) and .008 (b).
rc-child-adult-timecourse and supp-table4-other-networks¶
Network timecourses before z-normalization, averaged over the participants of each age group (artifact timepoints left out), correlated with the adult average over TRs 11 to 168.
In [40]:
MOVIE = np.arange(sb.FIRST_MOVIE_TR - 1, sb.N_TR)
def group_r(net, idx):
a = np.nanmean(RAW[net].loc[adults.index].values, axis=0)[MOVIE]
c = np.nanmean(RAW[net].loc[idx].values, axis=0)[MOVIE]
res = stats.pearsonr(a, c)
return res.statistic, res.pvalue
rows = []
for net in NETS:
for g in sb.AGE_GROUPS:
r, pv = group_r(net, children[children.AgeGroup == g].index)
rows.append({"network": net, "group": g, "r": r, "p": pv})
r, pv = group_r(net, children.index)
rows.append({"network": net, "group": "all children", "r": r, "p": pv})
gr = pd.DataFrame(rows)
print(gr.pivot(index="network", columns="group", values="r").round(3).to_string())
print(gr.pivot(index="network", columns="group", values="p").map(lambda x: f"{x:.1e}").to_string())
group 3yo 4yo 5yo 7yo 8-12yo all children network FFA 0.589 0.621 0.638 0.812 0.835 0.818 Face 0.637 0.634 0.660 0.825 0.848 0.828 M1 0.265 0.186 0.434 0.350 0.315 0.428 Pain 0.531 0.516 0.701 0.797 0.862 0.803 Scene 0.636 0.684 0.781 0.772 0.847 0.814 ToM 0.212 0.329 0.564 0.703 0.785 0.696 V1 0.442 0.654 0.742 0.713 0.773 0.753 group 3yo 4yo 5yo 7yo 8-12yo all children network FFA 3.8e-16 3.4e-18 1.9e-19 2.6e-38 2.2e-42 2.8e-39 Face 2.3e-19 3.7e-19 4.4e-21 1.9e-40 9.2e-45 4.8e-41 M1 7.8e-04 1.9e-02 1.2e-08 6.6e-06 5.5e-05 1.9e-08 Pain 7.1e-13 4.0e-12 1.3e-24 6.4e-36 7.4e-48 6.0e-37 Scene 2.6e-19 4.3e-23 1.1e-33 1.6e-32 1.3e-44 1.0e-38 ToM 7.6e-03 2.5e-05 1.3e-14 7.2e-25 2.5e-34 3.1e-24 V1 6.0e-09 1.2e-20 6.5e-29 8.8e-26 1.3e-32 4.0e-30
In [41]:
EXP_TP = {"ToM": [.28, .31, .60, .72, .82], "Pain": [.60, .56, .73, .83, .89]}
# The paper's bounds (claims file statistics): ToM every bin p < .0005 (3yo exact p = .00046);
# Pain every bin p < 1e-13, 3yo p < 1e-15. Bonferroni alpha .01 for five bins.
BOUND = {"ToM": [.0005] * 5, "Pain": [1e-15, 1e-13, 1e-13, 1e-13, 1e-13]}
for net, e in EXP_TP.items():
g = gr[(gr.network == net) & gr.group.isin(sb.AGE_GROUPS)].copy(); g["expected"] = e
g["in tolerance"] = [tol(a, b, .1 * b) for a, b in zip(g.r, g.expected)]
g["paper bound"] = BOUND[net]
g["meets bound"] = g.p < g["paper bound"]
g["p < .01"] = g.p < .01
print(g.to_string(index=False, float_format=lambda x: f"{x:.3g}"))
network group r p expected in tolerance paper bound meets bound p < .01
ToM 3yo 0.212 0.00761 0.28 False 0.0005 False True
ToM 4yo 0.329 2.5e-05 0.31 True 0.0005 True True
ToM 5yo 0.564 1.27e-14 0.6 True 0.0005 True True
ToM 7yo 0.703 7.23e-25 0.72 True 0.0005 True True
ToM 8-12yo 0.785 2.51e-34 0.82 True 0.0005 True True
network group r p expected in tolerance paper bound meets bound p < .01
Pain 3yo 0.531 7.08e-13 0.6 False 1e-15 False True
Pain 4yo 0.516 4.01e-12 0.56 True 1e-13 False True
Pain 5yo 0.701 1.26e-24 0.73 True 1e-13 True True
Pain 7yo 0.797 6.43e-36 0.83 True 1e-13 True True
Pain 8-12yo 0.862 7.41e-48 0.89 True 1e-13 True True
Result: not reproduced. Every age group's mean timecourse correlates with the adult mean at p < .01 (the paper's Bonferroni alpha): ToM r = .212, .329, .564, .703, .785 and Pain .531, .516, .701, .797, .862. 8 of 10 values are within 10%; ToM 3yo (.212 against .28) and Pain 3yo (.531 against .60) are below. Three correlations miss the paper's bounds: ToM 3yo p = .0076 against p < .0005 (the paper's exact p is .00046), Pain 3yo p = 7.1e-13 against p < 1e-15, and Pain 4yo p = 4.0e-12 against p < 1e-13. The criterion's partial outcome needs every r significant at the paper's bounds, so it does not apply. Judged only at the Bonferroni alpha of .01, every r is significant and the claim would be partially reproduced.
In [42]:
EXP_S4 = {"Face": [.75, .59, .71, .84, .86], "Scene": [.61, .67, .77, .80, .85], "M1": [.11, .06, .35, .44, .53],
"V1": [.53, .61, .78, .74, .76], "FFA": [.72, .60, .68, .82, .83]}
rows = []
for net, e in EXP_S4.items():
g = gr[(gr.network == net) & gr.group.isin(sb.AGE_GROUPS)]
for (_, x), ex in zip(g.iterrows(), e):
rows.append({"network": net, "group": x.group, "r": x.r, "expected": ex, "in tolerance": tol(x.r, ex, .1 * ex),
"p": x.p, "significant (p < .0005)": x.p < .0005})
s4 = pd.DataFrame(rows)
print(s4.to_string(index=False, float_format=lambda x: f"{x:.3g}"))
print("in tolerance:", int(s4["in tolerance"].sum()), "of 25; not significant:", list(s4.loc[~s4["significant (p < .0005)"], ["network", "group"]].itertuples(index=False, name=None)))
network group r expected in tolerance p significant (p < .0005)
Face 3yo 0.637 0.75 False 2.33e-19 True
Face 4yo 0.634 0.59 True 3.68e-19 True
Face 5yo 0.66 0.71 True 4.35e-21 True
Face 7yo 0.825 0.84 True 1.88e-40 True
Face 8-12yo 0.848 0.86 True 9.18e-45 True
Scene 3yo 0.636 0.61 True 2.57e-19 True
Scene 4yo 0.684 0.67 True 4.35e-23 True
Scene 5yo 0.781 0.77 True 1.07e-33 True
Scene 7yo 0.772 0.8 True 1.57e-32 True
Scene 8-12yo 0.847 0.85 True 1.28e-44 True
M1 3yo 0.265 0.11 False 0.000781 False
M1 4yo 0.186 0.06 False 0.019 False
M1 5yo 0.434 0.35 False 1.19e-08 True
M1 7yo 0.35 0.44 False 6.56e-06 True
M1 8-12yo 0.315 0.53 False 5.54e-05 True
V1 3yo 0.442 0.53 False 5.98e-09 True
V1 4yo 0.654 0.61 True 1.18e-20 True
V1 5yo 0.742 0.78 True 6.48e-29 True
V1 7yo 0.713 0.74 True 8.82e-26 True
V1 8-12yo 0.773 0.76 True 1.29e-32 True
FFA 3yo 0.589 0.72 False 3.85e-16 True
FFA 4yo 0.621 0.6 True 3.38e-18 True
FFA 5yo 0.638 0.68 True 1.9e-19 True
FFA 7yo 0.812 0.82 True 2.58e-38 True
FFA 8-12yo 0.835 0.83 True 2.22e-42 True
in tolerance: 17 of 25; not significant: [('M1', '3yo'), ('M1', '4yo')]
Result: not reproduced. 17 of the 25 values are within 10% (the criterion's partial outcome needs 20). The significance pattern holds at p < .0005 (only 3yo M1, p = .00078, and 4yo M1, p = .019, are not), but the M1 values differ throughout (.265, .186, .434, .350, .315 against .11, .06, .35, .44, .53), and the 3yo Face (.637 against .75), V1 (.442 against .53) and FFA (.589 against .72) values are low.
rc-age-events, rc-t04-tom, rc-fb-no-difference¶
Response magnitude at the peak TR of each Supplementary Table 3 event, ToM events from the ToM network and pain events from the Pain network. O7a: z-normalized timecourse (O6 a or b); O7b: before z-normalization. Children whose peak TR is an artifact timepoint (NaN) are left out of that event's test; their number is printed.
In [43]:
MAG = {}
for name, src in {"O7a z movie": lambda n: NZ[(n, "movie")], "O7a z all": lambda n: NZ[(n, "all")], "O7b raw": lambda n: RAW[n]}.items():
MAG[name] = pd.DataFrame({row.event: src(row.network)[row.peak_tr - 1] for _, row in sb.SUPP_TABLE3.iterrows()})
rows = []
for name, mag in MAG.items():
for ev in mag.columns:
c = children.join(mag[[ev]]).dropna(subset=[ev])
rr, pv, _ = sb.partial_spearman(c[ev], c.Age, c.n_art)
rows.append({"reading": name, "event": ev, "n": len(c), "rs": rr, "p": pv, "passes p < .0026": pv < .0026 and rr > 0})
age_ev = pd.DataFrame(rows)
print(age_ev.pivot(index="event", columns="reading", values="rs").round(3).join(
age_ev.pivot(index="event", columns="reading", values="p").map(lambda x: f"{x:.1e}"), rsuffix=" p").to_string())
print("children left out (artifact at peak):", age_ev.groupby("event").n.first().rsub(len(children)).to_dict())
for name in MAG:
passed = list(age_ev[(age_ev.reading == name) & age_ev["passes p < .0026"]].event)
print(name, "events passing:", passed, "; with rs > .28:", list(age_ev[(age_ev.reading == name) & age_ev["passes p < .0026"] & (age_ev.rs > .28)].event))
print("paper: P01, P04, P08, T01, T02")
reading O7a z all O7a z movie O7b raw O7a z all p O7a z movie p O7b raw p
event
P01 0.198 0.193 0.067 3.3e-02 3.8e-02 4.8e-01
P02 0.104 0.101 -0.060 2.6e-01 2.7e-01 5.1e-01
P03 0.186 0.180 0.067 4.5e-02 5.3e-02 4.7e-01
P04 0.325 0.324 0.187 2.9e-04 3.1e-04 4.1e-02
P05 0.157 0.152 0.079 8.7e-02 9.8e-02 3.9e-01
P06 0.251 0.249 0.183 6.3e-03 6.7e-03 4.8e-02
P07 0.171 0.168 0.196 7.0e-02 7.6e-02 3.8e-02
P08 0.021 0.019 -0.108 8.3e-01 8.4e-01 2.5e-01
P09 0.212 0.209 0.177 2.3e-02 2.5e-02 6.0e-02
P10 0.061 0.060 -0.016 5.1e-01 5.2e-01 8.6e-01
P11 0.072 0.073 -0.023 4.5e-01 4.4e-01 8.1e-01
P12 0.009 0.002 -0.072 9.2e-01 9.8e-01 4.4e-01
T01 0.271 0.270 0.173 3.5e-03 3.6e-03 6.5e-02
T02 0.325 0.325 0.251 4.5e-04 4.5e-04 7.4e-03
T03 0.275 0.274 0.249 2.5e-03 2.6e-03 6.4e-03
T04 0.171 0.175 0.135 6.9e-02 6.3e-02 1.6e-01
T05 0.184 0.188 0.171 4.6e-02 4.1e-02 6.4e-02
T06 0.055 0.061 -0.035 5.6e-01 5.2e-01 7.1e-01
T07 0.356 0.352 0.200 7.2e-05 8.5e-05 2.9e-02
children left out (artifact at peak): {'P01': 5, 'P02': 1, 'P03': 4, 'P04': 1, 'P05': 1, 'P06': 4, 'P07': 8, 'P08': 5, 'P09': 7, 'P10': 2, 'P11': 7, 'P12': 6, 'T01': 7, 'T02': 8, 'T03': 2, 'T04': 8, 'T05': 3, 'T06': 7, 'T07': 2}
O7a z movie events passing: ['T02', 'T03', 'T07', 'P04'] ; with rs > .28: ['T02', 'T07', 'P04']
O7a z all events passing: ['T02', 'T03', 'T07', 'P04'] ; with rs > .28: ['T02', 'T07', 'P04']
O7b raw events passing: [] ; with rs > .28: []
paper: P01, P04, P08, T01, T02
/tmp/ipykernel_2510/1048283978.py:11: Pandas4Warning: The 'lsuffix' and 'rsuffix' keywords in DataFrame.join are deprecated and will be removed in a future version. Use the 'suffixes' keyword instead. print(age_ev.pivot(index="event", columns="reading", values="rs").round(3).join(
Result: not reproduced. Under the z-normalized magnitudes (O7a, either O6 span), four events pass p < .0026 with positive rs: T02 (.325), T07 (.352), P04 (.324) and T03 (.274). Of the paper's five (P01, P04, P08, T01, T02), only P04 and T02 pass; P01 (rs = .193, p = .038), T01 (.270, p = .0036) and P08 (.019, p = .84) do not, and T03 and T07 pass where the paper has them not. Under O7b no event passes.
In [44]:
rows = []
for name, mag in MAG.items():
for ev in [e for e in mag.columns if e.startswith("T")]:
c = children.join(mag[[ev]]).dropna(subset=[ev])
f = sb.ols(c, ev, [TOM, "Age", "n_art"]); fz = sb.ols(c, ev, [TOM, "Age", "n_art"], zscore=True)
rows.append({"reading": name, "event": ev, "n": int(f.nobs), "b raw": f.params[TOM], "b z": fz.params[TOM], "t": f.tvalues[TOM], "p": f.pvalues[TOM], "p < .007": f.pvalues[TOM] < .007})
t04 = pd.DataFrame(rows)
print(t04.to_string(index=False, float_format=lambda x: f"{x:.3g}"))
print("paper T04: b = .4, t = 2.98, p = .0035; only T04 marked")
reading event n b raw b z t p p < .007
O7a z movie T01 115 -0.687 -0.224 -1.6 0.113 False
O7a z movie T02 114 0.216 0.0603 0.425 0.672 False
O7a z movie T03 120 0.253 0.0843 0.603 0.547 False
O7a z movie T04 114 0.528 0.171 1.18 0.24 False
O7a z movie T05 119 0.252 0.0848 0.593 0.555 False
O7a z movie T06 115 -0.413 -0.127 -0.851 0.396 False
O7a z movie T07 120 -0.158 -0.0463 -0.335 0.738 False
O7a z all T01 115 -0.719 -0.238 -1.7 0.0928 False
O7a z all T02 114 0.173 0.049 0.345 0.73 False
O7a z all T03 120 0.26 0.0873 0.624 0.534 False
O7a z all T04 114 0.502 0.165 1.14 0.258 False
O7a z all T05 119 0.231 0.0782 0.546 0.586 False
O7a z all T06 115 -0.409 -0.128 -0.86 0.392 False
O7a z all T07 120 -0.149 -0.0443 -0.32 0.749 False
O7b raw T01 115 -4.81 -0.275 -1.89 0.0614 False
O7b raw T02 114 0.809 0.0469 0.323 0.748 False
O7b raw T03 120 0.722 0.0497 0.353 0.725 False
O7b raw T04 114 1.6 0.103 0.705 0.482 False
O7b raw T05 119 -0.749 -0.0456 -0.319 0.75 False
O7b raw T06 115 -2.15 -0.117 -0.791 0.431 False
O7b raw T07 120 1.37 0.0751 0.521 0.604 False
paper T04: b = .4, t = 2.98, p = .0035; only T04 marked
Result: not reproduced. The T04 magnitude does not rise with ToM score with age and motion in the model: t = 1.18, p = .24 (O7a, TRs 11 to 168; t = 1.14 and .71 under the other readings), against t = 2.98, p = .0035. No ToM event passes p < .007 under any reading.
In [45]:
rng_fb = np.random.default_rng(1027)
rows = []
for name, mag in MAG.items():
for ev in [e for e in mag.columns if e.startswith("T")]:
c = young.loc[pf_idx].join(mag[[ev]]).dropna(subset=[ev]).copy(); c["fail"] = (c.FB_Group == "fail").astype(float)
f = sb.ols(c, ev, ["fail", "Age", "n_art"])
row = {"reading": name, "event": ev, "n": int(f.nobs), "mean pass": c.loc[c.fail == 0, ev].mean(), "mean fail": c.loc[c.fail == 1, ev].mean(),
"t": f.tvalues["fail"], "p": f.pvalues["fail"]}
if ev == "T03":
t_obs = f.tvalues["fail"]; perm = []
for _ in range(5000):
cc = c.copy(); cc["fail"] = rng_fb.permutation(cc.fail.values)
perm.append(sb.ols(cc, ev, ["fail", "Age", "n_art"]).tvalues["fail"])
row["p permutation (T03)"] = float(np.mean(np.abs(perm) >= abs(t_obs)))
rows.append(row)
fbev = pd.DataFrame(rows)
print(fbev.to_string(index=False, float_format=lambda x: f"{x:.3g}"))
print("smallest p per reading:", fbev.groupby("reading").p.min().round(4).to_dict(), "(paper: all p > .08)")
reading event n mean pass mean fail t p p permutation (T03)
O7a z movie T01 42 0.101 0.0533 -0.23 0.819 NaN
O7a z movie T02 40 0.619 0.352 -0.0392 0.969 NaN
O7a z movie T03 45 0.00888 0.139 0.366 0.716 0.721
O7a z movie T04 43 0.203 0.0236 0.102 0.919 NaN
O7a z movie T05 44 0.151 0.205 1.73 0.0911 NaN
O7a z movie T06 43 0.167 0.457 1.53 0.135 NaN
O7a z movie T07 44 0.291 0.289 0.418 0.678 NaN
O7a z all T01 42 0.102 0.0657 -0.14 0.889 NaN
O7a z all T02 40 0.611 0.36 0.0286 0.977 NaN
O7a z all T03 45 0.0231 0.151 0.444 0.659 0.663
O7a z all T04 43 0.206 0.0362 0.178 0.86 NaN
O7a z all T05 44 0.157 0.218 1.85 0.0722 NaN
O7a z all T06 43 0.172 0.461 1.59 0.12 NaN
O7a z all T07 44 0.305 0.299 0.478 0.635 NaN
O7b raw T01 42 0.244 0.134 -0.247 0.806 NaN
O7b raw T02 40 2.05 1.58 -0.0235 0.981 NaN
O7b raw T03 45 -0.148 0.368 0.435 0.666 0.669
O7b raw T04 43 0.249 -0.0357 0.783 0.439 NaN
O7b raw T05 44 -0.176 0.49 2.02 0.0497 NaN
O7b raw T06 43 0.823 3.5 1.76 0.0866 NaN
O7b raw T07 44 1.81 1.66 0.017 0.987 NaN
smallest p per reading: {'O7a z all': 0.0722, 'O7a z movie': 0.0911, 'O7b raw': 0.0497} (paper: all p > .08)
Result: reproduced, resting on an assumption. Under the z-normalized magnitudes (O7a), no ToM event differs between passers and failers with age and motion in the model (smallest p = .091 for T05 with O6a, .072 with O6b; T03 permutation p = .72). Under O7b, magnitudes before z-normalization, T05 has p = .0497, just below .05. O7a is adopted (see deviations).
rc-3yo-events¶
Reverse correlation on the 3-year-olds alone finds 4 of the 12 adult pain events and 1 of the 7 adult ToM events; of its 32 s of pain events 14 s fall in adult pain events and 8 s in adult ToM event T04, and of its 8 s of ToM events 4 s fall in an adult ToM event and 4 s in adult pain event P01.
Adult events are those of Supplementary Table 3 (Figure 4a); the events found here in adults are shown as a variant.
In [46]:
def tr_set(evs):
return set(t for e in evs for t in e)
TAB = {n: dict(zip(sb.SUPP_TABLE3[sb.SUPP_TABLE3.network == n].event, sb.SUPP_TABLE3[sb.SUPP_TABLE3.network == n].trs)) for n in ["ToM", "Pain"]}
rows = []
for span in ["movie", "all"]:
for k, (alpha, fdr) in O5.items():
e3 = {n: sb.reverse_correlation(NZ[(n, span)].loc[y3].values, alpha, fdr)[1] for n in ["ToM", "Pain"]}
row = {"O6": span, "O5": k}
for n, other in [("Pain", "ToM"), ("ToM", "Pain")]:
s3 = tr_set(e3[n])
row[f"3yo {n} events"] = len(e3[n]); row[f"3yo {n} s"] = 2 * len(s3)
row[f"adult {n} events found"] = sum(bool(s3 & set(v)) for v in TAB[n].values())
row[f"{n} s in adult {n} events"] = 2 * len(s3 & tr_set(TAB[n].values()))
row[f"{n} s in adult {other} events"] = 2 * len(s3 & tr_set(TAB[other].values()))
row["Pain s in T04"] = 2 * len(tr_set(e3["Pain"]) & set(TAB["ToM"]["T04"]))
row["ToM s in P01"] = 2 * len(tr_set(e3["ToM"]) & set(TAB["Pain"]["P01"]))
row["3yo events (TR lists)"] = {n: [(e[0], len(e)) for e in e3[n]] for n in e3}
row["vs found adult events: Pain s"] = 2 * len(tr_set(e3["Pain"]) & tr_set(RC[("Pain", span, k)][1]))
row["vs found adult events: ToM s"] = 2 * len(tr_set(e3["ToM"]) & tr_set(RC[("ToM", span, k)][1]))
rows.append(row)
r3 = pd.DataFrame(rows)
print(r3.drop(columns="3yo events (TR lists)").T.to_string())
for _, x in r3.iterrows():
print(x.O6, x.O5, x["3yo events (TR lists)"])
print("paper: 3yo pain 32 s, 4 of 12 adult pain events, 14 s in adult pain events, 8 s in T04; 3yo ToM 8 s, 1 of 7 adult ToM events, 4 s in a ToM event, 4 s in P01")
print("reading of rc-events used by the criterion:", BEST_SPAN, BEST_O5)
0 1 2 3 4 5
O6 movie movie movie all all all
O5 a p<.05 b p<.01 c FDR<.05 a p<.05 b p<.01 c FDR<.05
3yo Pain events 7 4 2 7 4 2
3yo Pain s 40 20 8 40 20 8
adult Pain events found 6 4 2 6 4 2
Pain s in adult Pain events 22 12 6 22 12 6
Pain s in adult ToM events 6 0 0 6 0 0
3yo ToM events 4 2 0 5 3 0
3yo ToM s 20 8 0 24 12 0
adult ToM events found 1 0 0 1 1 0
ToM s in adult ToM events 4 0 0 4 2 0
ToM s in adult Pain events 4 4 0 8 4 0
Pain s in T04 6 0 0 6 0 0
ToM s in P01 0 0 0 4 0 0
vs found adult events: Pain s 24 10 6 24 10 6
vs found adult events: ToM s 8 0 0 8 0 0
movie a p<.05 {'ToM': [(np.int64(29), 2), (np.int64(125), 3), (np.int64(139), 2), (np.int64(162), 3)], 'Pain': [(np.int64(40), 3), (np.int64(49), 2), (np.int64(57), 2), (np.int64(73), 3), (np.int64(112), 3), (np.int64(118), 4), (np.int64(152), 3)]}
movie b p<.01 {'ToM': [(np.int64(29), 2), (np.int64(163), 2)], 'Pain': [(np.int64(40), 2), (np.int64(49), 2), (np.int64(57), 2), (np.int64(118), 4)]}
movie c FDR<.05 {'ToM': [], 'Pain': [(np.int64(57), 2), (np.int64(119), 2)]}
all a p<.05 {'ToM': [(np.int64(29), 2), (np.int64(116), 2), (np.int64(125), 3), (np.int64(139), 2), (np.int64(162), 3)], 'Pain': [(np.int64(40), 3), (np.int64(49), 2), (np.int64(57), 2), (np.int64(73), 3), (np.int64(112), 3), (np.int64(118), 4), (np.int64(152), 3)]}
all b p<.01 {'ToM': [(np.int64(29), 2), (np.int64(126), 2), (np.int64(163), 2)], 'Pain': [(np.int64(40), 2), (np.int64(49), 2), (np.int64(57), 2), (np.int64(118), 4)]}
all c FDR<.05 {'ToM': [], 'Pain': [(np.int64(57), 2), (np.int64(119), 2)]}
paper: 3yo pain 32 s, 4 of 12 adult pain events, 14 s in adult pain events, 8 s in T04; 3yo ToM 8 s, 1 of 7 adult ToM events, 4 s in a ToM event, 4 s in P01
reading of rc-events used by the criterion: all a p<.05
Result: partially reproduced, resting on an assumption. No reading reproduced rc-events, so the criterion's "reading that reproduced rc-events" is extended to the reading that matched it best (O5a). Under O5a with O6b, the 3yo pain events (40 s) overlap 6 of the 12 adult pain events (22 s in adult pain events) and adult ToM event T04 (6 s), and the 3yo ToM events (24 s) overlap 1 of 7 adult ToM events (4 s) and adult pain event P01 (4 s). The two cross-network overlaps of the claim hold; the times and the pain-event count do not (paper 32 s, 4 of 12, 14 s, 8 s in T04; 8 s of ToM events). Under O6a the 3yo ToM events do not overlap P01 (0 s), so the claim would not be reproduced; O6b is adopted (see deviations).
Functional maturity¶
Each child's ToM (Pain) network timecourse before z-normalization, TRs 11 to 168, artifact timepoints left out, correlated (Pearson) with the adult mean timecourse (Methods paragraph 13). IRC measures are those of irc-age-within (M1 regressed), under each IRC reading.
In [47]:
def maturity(net):
adult_mean = np.nanmean(RAW[net].loc[adults.index].values, axis=0)[MOVIE]
out = {}
for s in children.index:
x = RAW[net].loc[s].values[MOVIE]; ok = np.isfinite(x)
out[s] = stats.pearsonr(x[ok], adult_mean[ok]).statistic
return pd.Series(out)
MAT = pd.DataFrame({"mat_ToM": maturity("ToM"), "mat_Pain": maturity("Pain")})
print(MAT.describe().round(3).to_string())
mat_ToM mat_Pain count 122.000 122.000 mean 0.241 0.332 std 0.177 0.165 min -0.211 -0.118 25% 0.106 0.231 50% 0.263 0.338 75% 0.369 0.437 max 0.539 0.666
In [48]:
rows = []
for r in READINGS:
c = IRC[r][0].loc[children.index].join(MAT)
for net, wv in [("ToM", "wToM"), ("Pain", "wPain")]:
f = sb.ols(c, f"mat_{net}", ["across", wv, "Age", "n_art"]); fz = sb.ols(c, f"mat_{net}", ["across", wv, "Age", "n_art"], zscore=True)
rows.append({"reading": reading_name(r), "network": net, "b across raw": f.params["across"], "b across z": fz.params["across"],
"t across": f.tvalues["across"], "p across": f.pvalues["across"], "b within z": fz.params[wv], "t within": f.tvalues[wv], "p within": f.pvalues[wv]})
print(pd.DataFrame(rows).to_string(float_format=lambda x: f"{x:.3g}"))
print("paper: ToM across b = -.4, t = -5.5, within t = 1.5 (p = .14); Pain across b = -.51, t = -7.4, within t = 1.9 (p = .06)")
reading network b across raw b across z t across p across b within z t within p within 0 O1 all, O2a, O3a ToM -0.679 -0.336 -4.52 1.47e-05 0.112 1.53 0.129 1 O1 all, O2a, O3a Pain -0.534 -0.284 -3.55 0.000549 0.204 2.4 0.0178 2 O1 all, O2a, O3b ToM -0.655 -0.326 -4.37 2.71e-05 0.119 1.61 0.11 3 O1 all, O2a, O3b Pain -0.577 -0.308 -3.9 0.000162 0.203 2.4 0.0181 4 O1 all, O2b, O3a ToM -0.715 -0.36 -4.85 3.76e-06 0.143 1.95 0.0533 5 O1 all, O2b, O3a Pain -0.575 -0.311 -3.9 0.000159 0.226 2.68 0.00834 6 O1 all, O2b, O3b ToM -0.687 -0.347 -4.62 1e-05 0.147 1.99 0.0492 7 O1 all, O2b, O3b Pain -0.598 -0.324 -4.1 7.69e-05 0.237 2.84 0.00537 8 O1 with_data, O2a, O3a ToM -0.679 -0.336 -4.52 1.47e-05 0.112 1.53 0.129 9 O1 with_data, O2a, O3a Pain -0.534 -0.284 -3.55 0.000549 0.204 2.4 0.0178 10 O1 with_data, O2a, O3b ToM -0.655 -0.326 -4.37 2.71e-05 0.119 1.61 0.11 11 O1 with_data, O2a, O3b Pain -0.576 -0.308 -3.9 0.000162 0.203 2.4 0.0181 12 O1 with_data, O2b, O3a ToM -0.715 -0.36 -4.85 3.76e-06 0.143 1.95 0.0533 13 O1 with_data, O2b, O3a Pain -0.575 -0.311 -3.9 0.000159 0.226 2.68 0.00834 14 O1 with_data, O2b, O3b ToM -0.687 -0.347 -4.62 1e-05 0.147 1.99 0.0492 15 O1 with_data, O2b, O3b Pain -0.598 -0.324 -4.1 7.69e-05 0.237 2.84 0.00537 paper: ToM across b = -.4, t = -5.5, within t = 1.5 (p = .14); Pain across b = -.51, t = -7.4, within t = 1.9 (p = .06)
Result: not reproduced. With age and motion in the model, a stronger across-network anti-correlation predicts a more adult-like timecourse in both networks (ToM t = -4.52, p = 1.5e-05; Pain t = -3.55, p = .00055), with the paper's sign, but both t values are out of tolerance (-5.5 and -7.4; -4.37 to -4.85 and -3.55 to -4.10 under the readings). Within-ToM does not predict ToM maturity under the main reading (t = 1.53, p = .13), as in the paper; under O2b with O3b it does (t = 1.99, p = .049). Within-Pain does predict Pain maturity (t = 2.40, p = .018; p = .0054 to .018 under every reading), where the paper's effect was not significant (p = .06). The criterion requires the within-network effects to be not significant, so the claim is not reproduced.
In [49]:
rows = []
for r in READINGS:
c = IRC[r][0].loc[children.index].join(MAT)
for net, wv in [("ToM", "wToM"), ("Pain", "wPain")]:
f = sb.ols(c, wv, [f"mat_{net}", "Age", "n_art"]); fz = sb.ols(c, wv, [f"mat_{net}", "Age", "n_art"], zscore=True)
rows.append({"reading": reading_name(r), "measure": wv, "b maturity z": fz.params[f"mat_{net}"], "t maturity": f.tvalues[f"mat_{net}"],
"p maturity": f.pvalues[f"mat_{net}"], "b age z": fz.params["Age"], "t age": f.tvalues["Age"], "p age": f.pvalues["Age"]})
print(pd.DataFrame(rows).to_string(float_format=lambda x: f"{x:.3g}"))
print("paper: maturity t = .48 (ToM), .61 (Pain), n.s.; age t = 2.5 (p = .01), 1.96 (p = .05)")
reading measure b maturity z t maturity p maturity b age z t age p age 0 O1 all, O2a, O3a wToM 0.0896 0.831 0.407 0.248 2.39 0.0186 1 O1 all, O2a, O3a wPain 0.168 1.8 0.0745 0.395 4.22 4.87e-05 2 O1 all, O2a, O3b wToM 0.1 0.936 0.351 0.251 2.43 0.0166 3 O1 all, O2a, O3b wPain 0.163 1.75 0.0832 0.403 4.32 3.24e-05 4 O1 all, O2b, O3a wToM 0.112 1.05 0.295 0.257 2.5 0.0137 5 O1 all, O2b, O3a wPain 0.181 1.93 0.0555 0.39 4.18 5.68e-05 6 O1 all, O2b, O3b wToM 0.118 1.11 0.268 0.264 2.59 0.0109 7 O1 all, O2b, O3b wPain 0.194 2.08 0.0395 0.383 4.1 7.53e-05 8 O1 with_data, O2a, O3a wToM 0.0896 0.831 0.407 0.248 2.39 0.0186 9 O1 with_data, O2a, O3a wPain 0.168 1.8 0.0745 0.395 4.22 4.87e-05 10 O1 with_data, O2a, O3b wToM 0.1 0.936 0.351 0.251 2.43 0.0166 11 O1 with_data, O2a, O3b wPain 0.163 1.75 0.0832 0.403 4.32 3.24e-05 12 O1 with_data, O2b, O3a wToM 0.112 1.05 0.295 0.257 2.5 0.0136 13 O1 with_data, O2b, O3a wPain 0.181 1.93 0.0555 0.39 4.18 5.68e-05 14 O1 with_data, O2b, O3b wToM 0.118 1.11 0.268 0.264 2.59 0.0109 15 O1 with_data, O2b, O3b wPain 0.194 2.08 0.0395 0.383 4.1 7.53e-05 paper: maturity t = .48 (ToM), .61 (Pain), n.s.; age t = 2.5 (p = .01), 1.96 (p = .05)
Result: reproduced, resting on an assumption. Under the main reading neither maturity effect is significant (within-ToM t = .83, p = .41; within-Pain t = 1.80, p = .075), while age is (t = 2.39, p = .019; t = 4.22, p = 4.9e-05). Under O2b with O3b the Pain maturity effect has p = .040; the main reading is adopted (see deviations).
In [50]:
c = children.join(MAT)
rr, pv, df = sb.partial_spearman(c.mat_ToM, c[TOM], c.n_art)
print(f"partial Spearman(ToM maturity, ToM score | motion) = {rr:.3f}, p = {pv:.2g}, df = {df} (paper .54, p = 1.3e-10)")
f = sb.ols(c, "mat_ToM", ["Age", TOM, "n_art"]); fz = sb.ols(c, "mat_ToM", ["Age", TOM, "n_art"], zscore=True)
print(f"regression: ToM score b z = {fz.params[TOM]:.2f}, t = {f.tvalues[TOM]:.2f}, p = {f.pvalues[TOM]:.3f}; age b z = {fz.params['Age']:.2f}, t = {f.tvalues['Age']:.2f}, p = {f.pvalues['Age']:.2g}")
print("paper: ToM score b = .1, t = 1.2, p = .3; age b = .5, t = 4.2, p < .00005")
partial Spearman(ToM maturity, ToM score | motion) = 0.498, p = 6.1e-09, df = 119 (paper .54, p = 1.3e-10) regression: ToM score b z = 0.00, t = 0.04, p = 0.966; age b z = 0.52, t = 4.46, p = 1.8e-05 paper: ToM score b = .1, t = 1.2, p = .3; age b = .5, t = 4.2, p < .00005
Result: reproduced. ToM maturity correlates with ToM score controlling for motion (partial rs = .498, p = 6.1e-09; paper .54, tolerance .054). With age in the regression the ToM score effect vanishes (t = .04, p = .97) while age remains (t = 4.46, p = 1.8e-05).
In [51]:
c = young.loc[pf_idx].join(MAT).copy(); c["fail"] = (c.FB_Group == "fail").astype(float)
f = sb.ols(c, "mat_ToM", ["fail", "Age", "n_art"]); fz = sb.ols(c, "mat_ToM", ["fail", "Age", "n_art"], zscore=True)
print(c.groupby("FB_Group").mat_ToM.agg(["mean", "sem", "count"]).round(3).to_string())
print(f"FB group (fail = 1): b raw = {f.params['fail']:.3f}, b z = {fz.params['fail']:.2f}, t = {f.tvalues['fail']:.2f}, p = {f.pvalues['fail']:.2f} (paper b = -.12, t = -.35, p = .73)")
mean sem count FB_Group fail 0.095 0.038 15 pass 0.168 0.032 30 FB group (fail = 1): b raw = 0.082, b z = 0.49, t = 1.53, p = 0.13 (paper b = -.12, t = -.35, p = .73)
Result: reproduced. ToM maturity does not differ between passers and failers with age and motion in the model (t = 1.53, p = .13; failers .095, passers .168).
In [52]:
r_all = group_r("FFA", children.index); r_3 = group_r("FFA", y3)
print(f"children (n = {len(children)}) against adults: r = {r_all[0]:.3f}, p = {r_all[1]:.1e} (paper .86, p < 1e-47)")
print(f"3yo (n = {len(y3)}) against adults: r = {r_3[0]:.3f}, p = {r_3[1]:.1e} (paper .72, p < 1e-26)")
rows = []
for span in ["movie", "all"]:
for k, (alpha, fdr) in O5.items():
tab, ev, rk = RC[("FFA", span, k)]
e3 = sb.reverse_correlation(NZ[("FFA", span)].loc[y3].values, alpha, fdr)[1]
rows.append({"O6": span, "O5": k, "adult events": len(ev), "adult s": 2 * sum(map(len, ev)), "adult lengths (TRs)": [len(e) for e in rk.trs],
"3yo events": len(e3), "3yo s": 2 * sum(map(len, e3)), "3yo TRs in adult events": len(tr_set(e3) & tr_set(ev)), "3yo TR count": len(tr_set(e3))})
face = pd.DataFrame(rows)
print(face.to_string())
print("paper: adults 7 events, 88 s, lengths 16, 4, 5, 7, 4, 4, 4 TRs; 3yo 4 events, 18 s, 8 of 9 TRs in adult events")
tab, ev, rk = RC[("FFA", BEST_SPAN, BEST_O5)]
rows = []
for name, src in {"O7a z": NZ[("FFA", BEST_SPAN)], "O7b raw": RAW["FFA"]}.items():
for i, row in rk.iterrows():
c = children.join(src[row.peak_tr - 1].rename("mag")).dropna(subset=["mag"])
rr, pv, _ = sb.partial_spearman(c.mag, c.Age, c.n_art)
rows.append({"reading": name, "event": f"F{i + 1:02d}", "peak TR": row.peak_tr, "n": len(c), "rs": rr, "p": pv,
"p < .0071": pv < .0071, f"p < .05/{len(rk)}": pv < .05 / len(rk)})
print(f"face-event magnitudes against age (events of the {BEST_SPAN}, {BEST_O5} reading):")
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda x: f"{x:.3g}"))
print("paper: |rs| < .24, ps > .01, none significant at .0071")
children (n = 122) against adults: r = 0.818, p = 2.8e-39 (paper .86, p < 1e-47) 3yo (n = 17) against adults: r = 0.589, p = 3.8e-16 (paper .72, p < 1e-26)
O6 O5 adult events adult s adult lengths (TRs) 3yo events 3yo s 3yo TRs in adult events 3yo TR count 0 movie a p<.05 8 92 [16, 4, 6, 3, 7, 4, 3, 3] 5 24 6 12 1 movie b p<.01 8 72 [4, 3, 8, 5, 3, 7, 3, 3] 1 4 2 2 2 movie c FDR<.05 8 76 [4, 3, 9, 6, 3, 7, 3, 3] 0 0 0 0 3 all a p<.05 10 108 [16, 4, 6, 3, 7, 4, 5, 2, 5, 2] 6 32 10 16 4 all b p<.01 9 88 [6, 4, 9, 6, 3, 7, 3, 3, 3] 1 4 2 2 5 all c FDR<.05 9 90 [6, 4, 9, 6, 3, 7, 4, 3, 3] 0 0 0 0 paper: adults 7 events, 88 s, lengths 16, 4, 5, 7, 4, 4, 4 TRs; 3yo 4 events, 18 s, 8 of 9 TRs in adult events
face-event magnitudes against age (events of the all, a p<.05 reading): reading event peak TR n rs p p < .0071 p < .05/10 O7a z F01 11 97 0.0571 0.581 False False O7a z F02 152 115 0.000794 0.993 False False O7a z F03 160 114 0.203 0.0312 False False O7a z F04 117 117 0.123 0.19 False False O7a z F05 107 112 0.24 0.0111 False False O7a z F06 59 119 0.215 0.0192 False False O7a z F07 71 118 0.0746 0.424 False False O7a z F08 89 115 0.108 0.251 False False O7a z F09 41 119 0.234 0.0106 False False O7a z F10 83 116 -0.0833 0.376 False False O7b raw F01 11 97 -0.17 0.0979 False False O7b raw F02 152 115 -0.191 0.0422 False False O7b raw F03 160 114 -0.00813 0.932 False False O7b raw F04 117 117 -5.34e-05 1 False False O7b raw F05 107 112 0.17 0.0747 False False O7b raw F06 59 119 0.236 0.01 False False O7b raw F07 71 118 -0.0104 0.911 False False O7b raw F08 89 115 0.0249 0.793 False False O7b raw F09 41 119 0.0648 0.485 False False O7b raw F10 83 116 -0.119 0.204 False False paper: |rs| < .24, ps > .01, none significant at .0071
Result: not reproduced. The children's mean fusiform timecourse correlates with the adults' (r = .818, within 10% of .86, p = 2.8e-39), but the 3yo correlation is .589 against .72, outside tolerance and above the paper's bound (p = 3.8e-16 against p < 1e-26). Reverse correlation finds 10 adult face events (108 s) and 6 three-year-old events (32 s) under the rc-events reading (O5a, O6b), against 7 (88 s) and 4 (18 s); no reading gives 7 adult events. No face-event magnitude changes with age at alpha .0071 (smallest p = .010), as in the paper.
In [53]:
rng_face = np.random.default_rng(88)
rows = []
for span, k in [(BEST_SPAN, BEST_O5)] + [("movie", k) for k in O5 if k != BEST_O5 or BEST_SPAN != "movie"]:
ef = RC[("FFA", span, k)][1]
for other in ["ToM", "Pain"]:
eo = RC[(other, span, k)][1]
nf, no = sum(map(len, ef)), sum(map(len, eo))
obs = len(tr_set(ef) & tr_set(eo))
perm = sb.overlap_permutation([len(e) for e in ef], [len(e) for e in eo], 158 - nf, 158 - no, 1000, rng_face)
rows.append({"reading": f"{span}, {k}", "face vs": other, "observed overlap TRs": obs, "perms <= observed": int((perm <= obs).sum()),
"p (fewer than chance)": float((perm <= obs).mean()), "mean perm overlap": perm.mean()})
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda x: f"{x:.3g}"))
print("paper: face-ToM 4 TRs, p = .11; face-pain 15 TRs, p = .93")
reading face vs observed overlap TRs perms <= observed p (fewer than chance) mean perm overlap
all, a p<.05 ToM 7 52 0.052 15
all, a p<.05 Pain 25 983 0.983 15.7
movie, a p<.05 ToM 7 98 0.098 13.1
movie, a p<.05 Pain 20 945 0.945 13.5
movie, b p<.01 ToM 3 166 0.166 7.47
movie, b p<.01 Pain 15 960 0.96 8.39
movie, c FDR<.05 ToM 3 136 0.136 7.85
movie, c FDR<.05 Pain 15 940 0.94 8.8
paper: face-ToM 4 TRs, p = .11; face-pain 15 TRs, p = .93
Result: reproduced. The face events overlap the ToM events by 7 TRs and the pain events by 25 TRs (paper 4 and 15), and neither overlap is smaller than chance (p = .052 and .983 under the rc-events reading; .098 to .166 and .94 to .96 under the others).
Blocked claims¶
- replication-sample: the independent sample of 20 adults (Fedorenko and colleagues) is not in the snapshot, and the paper names no archive for it.
- hand-coded-events: the "mental", "social" and "pain" coding of Jacoby et al. (2016) is not in the snapshot (its README lists only the reverse-correlation events), and the paper does not give the coded times. The Saxe lab's web page offers the experiment script (movloc.zip, version of 29 March 2019) that holds them, with no DOI or repository commit to pin it, so it cannot be used as a pinned input.
In [54]:
# replication-sample: the dataset's README lists the events "that replicated across the two
# samples" (onset TR; duration in TRs). Which Supplementary Table 3 events each overlaps. This
# restates the authors' result; it is not the second sample's data.
import re
readme = sb.text("README")
blk = readme.split("TRs identified by RC analysis")[1].split("TIMING OF EVENTS")[0]
tom_txt, pain_txt = blk.split("Pain Event Onsets")
def parse(s):
return [tuple(int(v) for v in line.split(";")) for line in s.splitlines() if re.fullmatch(r"\s*\d+;\s*\d+\s*", line)]
for net, part in [("ToM", tom_txt), ("Pain", pain_txt)]:
listed = parse(part)
tab = sb.SUPP_TABLE3[sb.SUPP_TABLE3.network == net]
hit = {ev: [(o, d) for o, d in listed if set(range(o, o + d)) & set(trs)] for ev, trs in zip(tab.event, tab.trs)}
print(f"{net}: {len(listed)} replicated events listed in the README:", listed)
print(" table events overlapped:", [ev for ev, h in hit.items() if h], "; not overlapped:", [ev for ev, h in hit.items() if not h])
ToM: 7 replicated events listed in the README: [(46, 2), (52, 3), (63, 2), (91, 8), (122, 3), (129, 4), (153, 3)] table events overlapped: ['T01', 'T02', 'T03', 'T04', 'T05', 'T06', 'T07'] ; not overlapped: [] Pain: 9 replicated events listed in the README: [(38, 2), (49, 1), (56, 2), (71, 5), (100, 2), (108, 6), (117, 3), (134, 3), (159, 2)] table events overlapped: ['P01', 'P02', 'P03', 'P04', 'P05', 'P06', 'P07', 'P10', 'P12'] ; not overlapped: ['P08', 'P09', 'P11']
The README's replicated events overlap all seven ToM events and nine of the twelve pain events of Supplementary Table 3, all but P08, P09 and P11, as the replication-sample claim states, and they are consistent with the onset mapping (stimulus seconds / 2 + 8) used for the table. This restates the authors' result and does not test it on the second sample's data, so the claim stays blocked.
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 |
|---|---|
| participants | reproduced |
| supp-table1-behavior-by-age | reproduced |
| beh-tom-age | reproduced |
| beh-tom-dccs | reproduced |
| beh-fb-groups | reproduced |
| beh-fb-tom | reproduced |
| beh-fb-age | reproduced |
| beh-fb-dccs | reproduced |
| motion-child-adult | reproduced |
| motion-not-age-tom | reproduced |
| motion-not-fb | reproduced |
| irc-adult | partially reproduced |
| irc-age-within | not reproduced |
| irc-age-across | not reproduced |
| irc-adult-vs-child | partially reproduced |
| irc-control-networks | partially reproduced |
| irc-williams | not reproduced |
| irc-3yo-distinct | not reproduced |
| irc-3yo-not-anticorrelated | partially reproduced |
| irc-3yo-homologous | not reproduced |
| irc-3yo-nonbilateral | partially reproduced |
| irc-tom-behavior | partially reproduced |
| irc-tom-behavior-age | not reproduced (assumed reading) |
| irc-fb-within | not reproduced |
| irc-fb-within-dccs | not reproduced |
| irc-fb-across | not reproduced (assumed reading) |
| rc-events | partially reproduced |
| rc-overlap | partially reproduced |
| rc-child-adult-timecourse | not reproduced |
| supp-table4-other-networks | not reproduced |
| rc-age-events | not reproduced |
| rc-t04-tom | not reproduced |
| rc-fb-no-difference | reproduced (assumed reading) |
| rc-3yo-events | partially reproduced (assumed reading) |
| maturity-across | not reproduced |
| maturity-within | reproduced (assumed reading) |
| maturity-tom-behavior | reproduced |
| maturity-fb | reproduced |
| face-rc | not reproduced |
| face-overlap | reproduced |
| replication-sample | blocked: the independent adult sample is not public |
| hand-coded-events | blocked: the hand coding of Jacoby et al. 2016 is in no pinned public source |
How the analysis differs from the paper's. The analysis starts from the paper's own preprocessed images, artifact indicators, CompCor regressors and ROI masks, all deposited in the snapshot, so the differences are in the steps after them: the ROI-mean shortcut (exact for the linear steps the Methods describe), SPM's discrete cosine high-pass filter for "a cutoff of 100 s", an intercept in the nuisance regression, the ppcor forms of the partial correlations, and the readings O1 to O9 of the claims file, each run. The paper's analysis code is not public.
What was found about the dataset. Every behavioural and demographic number of the paper is recovered from participants.tsv: the counts exactly, the other values within tolerance (for example mean ToM scores .517 and .610 against .54 and .63, DCCS p .048 against .049, Inc-Fail p .015 against .01, DCCS b .148 against .16), beh-fb-tom under reading (a) only, and the artifact counts match the Methods to the decimal. The deposited ROI masks have the sizes of Supplementary Table 2 and are in MNI mm, centered on the stated coordinates, and the artifact indicators line up with spikes in the image data (median deviation 6.45 at the flagged timepoints against at most 2.6 one or two TRs away). With these data and the pipeline of the Methods, the adult inter-region correlations are higher within networks and less negative across networks than the paper reports (.604, .398, -.051 against .48, .35, -.17); in the 3-year-olds within-ToM is higher (.326 against .21) but within-Pain is lower (.179 against .23, and non-bilateral Pain .146 against .17), and the event-level tests of age and ToM score (rc-age-events, rc-t04-tom) and the false-belief group effect on within-ToM correlation do not come out as in the paper. The voxel-size column (named "Scanlog: Voxel slize") counts 2 seven-year-olds and 21 8-12-year-olds at 3.13 mm with no slice gap, where the Methods say 3 and 20, and the no-gap level is stored as an empty value where participants.json defines "None". The face-ROI coverage count differs from Supplementary Fig. 4 (19 participants with any face-ROI voxel without data, against 17).