Distinct Roles of Deep and Superficial Cortical Layers in Tone Prediction, Comparison, and Adaptation in Human Auditory Cortices.

Code ↔ Paper

13 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 13 matches

[1] § Results › Superficial layers represent shared expectation-stimulus preference dynamics ↔ finalized/PredAdapt_fMRI/ses2_regressions/ROI_and_layer_analysis.ipynb, lines 841–949 · score 0.78 · Planum Polare, Planum Temporale, posterior STG, FDR corrected, repetition suppression, Heschl
[2] § Results › Distinct laminar profiles of expectations trending towards deep layers and layer-uniform repetition suppression ↔ finalized/PredAdapt_fMRI/ses2_regressions/ROI_and_layer_analysis.ipynb, lines 841–949 · score 0.77 · Planum Polare, Planum Temporale, posterior STG, FDR corrected, Heschl, gyrus
[3] § Methods › Calibration of the auditory loudness waveform ↔ projects/PredAdapt_MEG/psychtoolbox_exp/mainpredsound_MEG.m, lines 83–152 · score 0.68 · piecewise cubic interpolation, equal loudness curve, intensity, tones, stimulus
[4] § Methods › Functional data analysis › Preprocessing and alignment of functional data ↔ finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/functional_plot.py, lines 445–570 · score 0.67 · distortion correction, motion correction, preprocessing, Topup, FSL, slice
[5] § Methods › Calibration of the auditory loudness waveform ↔ finalized/PredAdapt_fMRI/psychtoolbox_exp/mainpredsound.m, lines 83–144 · score 0.66 · piecewise cubic interpolation, loudness curve, intensity, tones, stimulus
[6] § Methods › Functional data analysis › Preprocessing and alignment of functional data ↔ finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/functional.py, lines 197–290 · score 0.63 · motion correction, preprocessing, rotations, translations, Topup, BrainVoyager
[7] § Results ↔ finalized/PredAdapt_fMRI/ses2_modelstims/DREX/run_DREX_model.m, lines 1–55 · score 0.60 · context hypotheses, REX model, predictive distribution, mixtures, belief, window
[8] § Results › Distinct laminar profiles of expectations trending towards deep layers and layer-uniform repetition suppression ↔ finalized/AssociativeLearning/analysis/Laminar_UnivariateAndStats.ipynb, lines 344–454 · score 0.54 · aSTG, FDR corrected, pSTG, laminar, PT, HG
[9] § Methods › Stimulus modelling ↔ finalized/PredAdapt_fMRI/psychtoolbox_exp/settings_tonotopy.m, lines 1–46 · score 0.54 · inter stimulus interval, ISI, offsets, onsets, duration, TR
[10] § Results › Superficial layers represent shared expectation-stimulus preference dynamics ↔ finalized/AssociativeLearning/analysis/Laminar_UnivariateAndStats.ipynb, lines 778–901 · score 0.54 · aSTG, FDR corrected, pSTG, interactions, PT, HG
[11] § Results › Superficial layers represent shared expectation-stimulus preference dynamics ↔ finalized/AssociativeLearning/analysis/Laminar_UnivariateAndStats.ipynb, lines 44–61 · score 0.53 · aSTG, pSTG, middle layers, deep layers, surface, PT
[12] § Methods › Stimulus modelling ↔ projects/PredAdapt_MEG/psychtoolbox_exp/settings_localizer.m, lines 1–45 · score 0.52 · inter stimulus interval, ISI, offsets, onsets, duration, repetition
[13] § Methods › Computational modelling › D-REX model ↔ finalized/PredAdapt_fMRI/ses2_modelstims/DREX/run_DREX_model.m, lines 1–55 · score 0.51 · Gaussian Mixture Model, GMM, REX, sequences, probability

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 1,858 lines · 88 KB · no license · 3 matches

# %% [markdown]
# # <center> $\Large{\text{Univariate Analysis}}$ <br> Laminar Level (1.96 Sound) </center>
# %%
%autoreload 2
# %%
import sys
import itertools
from glob import glob
from tqdm import tqdm
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import bvbabel
import matlab.engine
import scipy.io as sio
from scipy.stats import ttest_rel, ttest_1samp, friedmanchisquare
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.multitest import fdrcorrection
import scikit_posthocs as sp
sys.path.insert(0, '/mnt/hdd2/associative_learning/notebooks/toolbox/')
from fmri_toolbox import FMRIToolbox
# matlab and seaborn config
import matplotlib as mpl
mpl.rcParams['text.usetex'] = True
sns.set(style="whitegrid")
import matplotlib.lines as mlines
import matplotlib.gridspec as gridspec
def add_significance(ax, x1, x2, y, h, text):
if text=='ns': return
ax.plot([x1, x1, x2, x2], [y, y+h, y+h, y], lw=1.5, color='black')
ax.text((x1 + x2) * .5, y + h, text, ha='center', va='bottom', color='black', fontsize=20)
plt.rcParams['text.usetex'] = True
# %% [markdown]
# - Layer 1 (D) Deep
# - Layer 2 (M) Middle
# - Layer 3 (S) Surface
# %% [markdown]
# ### Load data
# %%
data_dir = f'/mnt/hdd2/associative_learning/analysis/estimates/single_trial_gm_mask'
mask_dir = '/mnt/hdd2/associative_learning/masks/drawn_regions/'
tbx = FMRIToolbox(data_dir, mask_dir=mask_dir)
subjects = ['S04', 'S05', 'S06', 'S07', 'S08', 'S09', 'S10', 'S11', 'S12', 'S15', 'S17']
regions = ['PP','HG', 'PT', 'aSTG', 'pSTG', 'pMTG', 'parsOrbitalis', 'parsTriangularis', 'TPOj']
hemispheres = ['LH', 'RH']
# load dataset
dataset = tbx.load_dataset(subjects, hemispheres, regions, extract_layers=True)
# %% [markdown]
# ### Plot voxel counts
# %%
fig = plt.figure(figsize=(14,4))
palette = sns.color_palette()
palette1 = sns.color_palette("muted")
palette2 = sns.color_palette("pastel")
regions_layers_list = np.array([[region + ' (D)', region + ' (M)', region + ' (S)'] for region in regions]).ravel()
ax = plt.subplot(1,2,1)
voxel_counts_lh_l1 = {region:np.mean([dataset[subject]['LH'][region]['betas']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l1_se = {region:np.std([dataset[subject]['LH'][region]['betas']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l1_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_lh_l1_se.values()]
voxel_counts_lh_l2 = {region:np.mean([dataset[subject]['LH'][region]['betas']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l2_se = {region:np.std([dataset[subject]['LH'][region]['betas']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l2_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_lh_l2_se.values()]
voxel_counts_lh_l3 = {region:np.mean([dataset[subject]['LH'][region]['betas']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l3_se = {region:np.std([dataset[subject]['LH'][region]['betas']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l3_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_lh_l3_se.values()]
x_positions1 = np.arange(len(voxel_counts_lh_l1))
ax.bar(x=x_positions1, height=voxel_counts_lh_l1.values(), capsize=2, width=0.3, color=palette, yerr=voxel_counts_lh_l1_ci)
x_positions2 = x_positions1 + 0.3
ax.bar(x=x_positions2, height=voxel_counts_lh_l2.values(), capsize=2, width=0.3, color=palette1, yerr=voxel_counts_lh_l2_ci)
x_positions3 = x_positions1 + 0.6
ax.bar(x=x_positions3, height=voxel_counts_lh_l3.values(), capsize=2, width=0.3, color=palette2, yerr=voxel_counts_lh_l3_ci)
ax.set_xticks(np.sort(np.concatenate([x_positions1, x_positions2, x_positions3])), regions_layers_list)
plt.xticks(rotation=90)
ax.set_title('LH')
##############
ax = plt.subplot(1,2,2)
voxel_counts_lh_l1 = {region:np.mean([dataset[subject]['RH'][region]['betas']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l1_se = {region:np.std([dataset[subject]['RH'][region]['betas']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l1_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_lh_l1_se.values()]
voxel_counts_lh_l2 = {region:np.mean([dataset[subject]['RH'][region]['betas']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l2_se = {region:np.std([dataset[subject]['RH'][region]['betas']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l2_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_lh_l2_se.values()]
voxel_counts_lh_l3 = {region:np.mean([dataset[subject]['RH'][region]['betas']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l3_se = {region:np.std([dataset[subject]['RH'][region]['betas']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l3_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_lh_l3_se.values()]
x_positions1 = np.arange(len(voxel_counts_lh_l1))
ax.bar(x=x_positions1, height=voxel_counts_lh_l1.values(), capsize=2, width=0.3, color=palette, yerr=voxel_counts_lh_l1_ci)
ax.bar(x=x_positions2, height=voxel_counts_lh_l2.values(), capsize=2, width=0.3, color=palette1, yerr=voxel_counts_lh_l2_ci)
ax.bar(x=x_positions3, height=voxel_counts_lh_l3.values(), capsize=2, width=0.3, color=palette2, yerr=voxel_counts_lh_l3_ci)
ax.set_xticks(np.sort(np.concatenate([x_positions1, x_positions2, x_positions3])), regions_layers_list)
plt.xticks(rotation=90)
ax.set_title('RH')
plt.suptitle('Voxel Counts')
plt.show()
# %% [markdown]
# ### Perform T-Tests
# %%
tbx.utils.calculate_tstats(dataset, trials='sound')
# %% [markdown]
# ### Check voxels count using one threshold
# %%
t_threshold = 1.96
tbx.utils.select_voxels(dataset, selection={'method': 'tstat', 'strategy': 't_val', 'threshold':t_threshold, 'drop_negative':True})
# %%
fig = plt.figure(figsize=(14,4))
palette = sns.color_palette()
palette1 = sns.color_palette("muted")
palette2 = sns.color_palette("pastel")
regions_layers_list = np.array([[region + ' (D)', region + ' (M)', region + ' (S)'] for region in regions]).ravel()
ax = plt.subplot(1,2,1)
voxel_counts_lh_l1 = {region:np.nanmean([dataset[subject]['LH'][region]['betas_selected']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l1_min = {region:np.nanmin([dataset[subject]['LH'][region]['betas_selected']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l1_se = {region:np.std([dataset[subject]['LH'][region]['betas_selected']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l1_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_lh_l1_se.values()]
voxel_counts_lh_l2 = {region:np.nanmean([dataset[subject]['LH'][region]['betas_selected']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l2_min = {region:np.nanmin([dataset[subject]['LH'][region]['betas_selected']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l2_se = {region:np.std([dataset[subject]['LH'][region]['betas_selected']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l2_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_lh_l2_se.values()]
voxel_counts_lh_l3 = {region:np.nanmean([dataset[subject]['LH'][region]['betas_selected']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l3_min = {region:np.nanmin([dataset[subject]['LH'][region]['betas_selected']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l3_se = {region:np.std([dataset[subject]['LH'][region]['betas_selected']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l3_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_lh_l3_se.values()]
x_positions1 = np.arange(len(voxel_counts_lh_l1))
ax.bar(x=x_positions1, height=voxel_counts_lh_l1.values(), capsize=2, width=0.3, color=palette[0], yerr=voxel_counts_lh_l1_ci)
x_positions2 = x_positions1 + 0.3
ax.bar(x=x_positions2, height=voxel_counts_lh_l2.values(), capsize=2, width=0.3, color=palette1[0], yerr=voxel_counts_lh_l2_ci)
x_positions3 = x_positions1 + 0.6
ax.bar(x=x_positions3, height=voxel_counts_lh_l3.values(), capsize=2, width=0.3, color=palette2[0], yerr=voxel_counts_lh_l3_ci)
ax.set_xticks(np.sort(np.concatenate([x_positions1, x_positions2, x_positions3])), regions_layers_list)
plt.xticks(rotation=90)
ax.set_title('LH')
##############
regions_layers_list = np.array([[region + ' (D)', region + ' (M)', region + ' (S)'] for region in regions]).ravel()
ax = plt.subplot(1,2,2)
voxel_counts_rh_l1 = {region:np.nanmean([dataset[subject]['RH'][region]['betas_selected']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l1_min = {region:np.nanmin([dataset[subject]['RH'][region]['betas_selected']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l1_se = {region:np.std([dataset[subject]['RH'][region]['betas_selected']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l1_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_rh_l1_se.values()]
voxel_counts_rh_l2 = {region:np.nanmean([dataset[subject]['RH'][region]['betas_selected']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l2_min = {region:np.nanmin([dataset[subject]['RH'][region]['betas_selected']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l2_se = {region:np.std([dataset[subject]['RH'][region]['betas_selected']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l2_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_rh_l2_se.values()]
voxel_counts_rh_l3 = {region:np.nanmean([dataset[subject]['RH'][region]['betas_selected']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l3_min = {region:np.nanmin([dataset[subject]['RH'][region]['betas_selected']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l3_se = {region:np.std([dataset[subject]['RH'][region]['betas_selected']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l3_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_rh_l3_se.values()]
x_positions1 = np.arange(len(voxel_counts_rh_l1))
ax.bar(x=x_positions1, height=voxel_counts_rh_l1.values(), capsize=2, width=0.3, color=palette[0], yerr=voxel_counts_rh_l1_ci)
ax.bar(x=x_positions2, height=voxel_counts_rh_l2.values(), capsize=2, width=0.3, color=palette1[0], yerr=voxel_counts_rh_l2_ci)
ax.bar(x=x_positions3, height=voxel_counts_rh_l3.values(), capsize=2, width=0.3, color=palette2[0], yerr=voxel_counts_rh_l3_ci)
ax.set_xticks(np.sort(np.concatenate([x_positions1, x_positions2, x_positions3])), regions_layers_list)
plt.xticks(rotation=90)
ax.set_title('RH')
plt.suptitle(r'Voxel Counts after Selection ($F \geq$ ' + str(t_threshold) + ')')
plt.show()
# %% [markdown]
# # Display Percentage of Voxels
# %%
fig = plt.figure(figsize=(14,4))
palette = sns.color_palette()
palette1 = sns.color_palette("muted")
palette2 = sns.color_palette("pastel")
regions_layers_list = np.array([[region + ' (D)', region + ' (M)', region + ' (S)'] for region in regions]).ravel()
ax = plt.subplot(1,2,1)
voxel_counts_lh_l1 = {region:np.nanmean([dataset[subject]['LH'][region]['betas_selected']['L1'].shape[0]*100/dataset[subject]['LH'][region]['betas']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l1_se = {region:np.std([dataset[subject]['LH'][region]['betas_selected']['L1'].shape[0]*100/dataset[subject]['LH'][region]['betas']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l1_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_lh_l1_se.values()]
voxel_counts_lh_l2 = {region:np.nanmean([dataset[subject]['LH'][region]['betas_selected']['L2'].shape[0]*100/dataset[subject]['LH'][region]['betas']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l2_se = {region:np.std([dataset[subject]['LH'][region]['betas_selected']['L2'].shape[0]*100/dataset[subject]['LH'][region]['betas']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l2_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_lh_l2_se.values()]
voxel_counts_lh_l3 = {region:np.nanmean([dataset[subject]['LH'][region]['betas_selected']['L3'].shape[0]*100/dataset[subject]['LH'][region]['betas']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l3_se = {region:np.std([dataset[subject]['LH'][region]['betas_selected']['L3'].shape[0]*100/dataset[subject]['LH'][region]['betas']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_lh_l3_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_lh_l3_se.values()]
x_positions1 = np.arange(len(voxel_counts_lh_l1))
ax.bar(x=x_positions1, height=voxel_counts_lh_l1.values(), capsize=2, width=0.3, color=palette[0], yerr=voxel_counts_lh_l1_ci)
x_positions2 = x_positions1 + 0.3
ax.bar(x=x_positions2, height=voxel_counts_lh_l2.values(), capsize=2, width=0.3, color=palette1[0], yerr=voxel_counts_lh_l2_ci)
x_positions3 = x_positions1 + 0.6
ax.bar(x=x_positions3, height=voxel_counts_lh_l3.values(), capsize=2, width=0.3, color=palette2[0], yerr=voxel_counts_lh_l3_ci)
ax.set_xticks(np.sort(np.concatenate([x_positions1, x_positions2, x_positions3])), regions_layers_list)
plt.xticks(rotation=90)
ax.set_title('LH')
##############
regions_layers_list = np.array([[region + ' (D)', region + ' (M)', region + ' (S)'] for region in regions]).ravel()
ax = plt.subplot(1,2,2)
voxel_counts_rh_l1 = {region:np.nanmean([dataset[subject]['RH'][region]['betas_selected']['L1'].shape[0]*100/dataset[subject]['RH'][region]['betas']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l1_se = {region:np.std([dataset[subject]['RH'][region]['betas_selected']['L1'].shape[0]*100/dataset[subject]['RH'][region]['betas']['L1'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l1_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_rh_l1_se.values()]
voxel_counts_rh_l2 = {region:np.nanmean([dataset[subject]['RH'][region]['betas_selected']['L2'].shape[0]*100/dataset[subject]['RH'][region]['betas']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l2_se = {region:np.std([dataset[subject]['RH'][region]['betas_selected']['L2'].shape[0]*100/dataset[subject]['RH'][region]['betas']['L2'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l2_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_rh_l2_se.values()]
voxel_counts_rh_l3 = {region:np.nanmean([dataset[subject]['RH'][region]['betas_selected']['L3'].shape[0]*100/dataset[subject]['RH'][region]['betas']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l3_se = {region:np.std([dataset[subject]['RH'][region]['betas_selected']['L3'].shape[0]*100/dataset[subject]['RH'][region]['betas']['L3'].shape[0] for subject in subjects]) for region in regions}
voxel_counts_rh_l3_ci = [1.96 * (std_dev / np.sqrt(len(subjects))) for std_dev in voxel_counts_rh_l3_se.values()]
x_positions1 = np.arange(len(voxel_counts_rh_l1))
ax.bar(x=x_positions1, height=voxel_counts_rh_l1.values(), capsize=2, width=0.3, color=palette[0], yerr=voxel_counts_rh_l1_ci)
ax.bar(x=x_positions2, height=voxel_counts_rh_l2.values(), capsize=2, width=0.3, color=palette1[0], yerr=voxel_counts_rh_l2_ci)
ax.bar(x=x_positions3, height=voxel_counts_rh_l3.values(), capsize=2, width=0.3, color=palette2[0], yerr=voxel_counts_rh_l3_ci)
ax.set_xticks(np.sort(np.concatenate([x_positions1, x_positions2, x_positions3])), regions_layers_list)
plt.xticks(rotation=90)
ax.set_title('RH')
plt.suptitle(r'Voxel Percentage after Selection ($F \geq$ ' + str(t_threshold) + ')')
plt.show()
# %% [markdown]
# ---
#
# # Plot Mean Activations
# %%
fig, ax = plt.subplots(2, 5, figsize=(12,5))
for region_idx, axis in enumerate(ax.ravel()):
if region_idx==9:continue
activations = {
val: [np.nanmean(dataset[subject]['LH'][regions[region_idx]]['betas_selected'][f'L{idx+1}'].values) for subject in subjects]
for idx, val in enumerate(['D', 'M', 'S'])
}
sns.barplot(activations, ax=axis)
axis.set_title(regions[region_idx])
plt.suptitle('Mean Beta Activations (LH)')
plt.tight_layout()
plt.show()
fig, ax = plt.subplots(2, 5, figsize=(12,5))
for region_idx, axis in enumerate(ax.ravel()):
if region_idx==9:continue
activations = {
val: [np.nanmean(dataset[subject]['RH'][regions[region_idx]]['betas_selected'][f'L{idx+1}'].values) for subject in subjects]
for idx, val in enumerate(['D', 'M', 'S'])
}
sns.barplot(activations, ax=axis)
axis.set_title(regions[region_idx])
plt.suptitle('Mean Beta Activations (RH)')
plt.tight_layout()
plt.show()
# %%
def non_parametric_permutation_ttest(betas):
n_subjects = betas.size
permutation_matrix = np.array(list(itertools.product([1,-1], repeat=n_subjects))).T
permutation_product = betas @ permutation_matrix
p_val = (sum(abs(permutation_product[0]) <= abs(permutation_product[1:-1])) + 1) / (permutation_product.size)
return p_val
# %% [markdown]
# ---
#
# # <center> $\Large\textbf{Congruent vs. Incongruent}$ </center>
#
# **Test the effect of prediction error.**
# %%
regex1 = '_cong$'
regex2 = '_incong$'
layers = ['L1', 'L2', 'L3']
tstats = {}; pvals = {}; set1 = {};set2 = {}
for layer_name, layer_id in zip(['Deep', 'Middle', 'Superficial'], ['L1', 'L2', 'L3']):
tstats[layer_id] = []; pvals[layer_id] = []; df_set1 = {}; df_set2 = {}
for region in regions:
df_set1[region] = []
df_set2[region] = []
for subject in subjects:
stimuli = tbx.utils.extract_stimuli(subject)
betas_l = tbx.utils.extract_condition_betas(stimuli=stimuli,
betas=dataset[subject]['LH'][region]['betas_selected'][layer_id].values.T,
stimulus_type='sound')
betas_r = tbx.utils.extract_condition_betas(stimuli=stimuli,
betas=dataset[subject]['RH'][region]['betas_selected'][layer_id].values.T,
stimulus_type='sound')
df_set1[region].append(np.mean([
np.nanmean(betas_l.filter(regex=regex1).values),
np.nanmean(betas_r.filter(regex=regex1).values)]))
df_set2[region].append(np.mean([
np.nanmean(betas_l.filter(regex=regex2).values),
np.nanmean(betas_r.filter(regex=regex2).values)]))
set1[layer_id] = df_set1
set2[layer_id] = df_set2
all_non_parametric_pvals = np.array([[non_parametric_permutation_ttest(np.subtract(set2[layer][region], set1[layer][region])) for layer in layers] for region in regions])
all_non_parametric_pvals_flat = np.array(all_non_parametric_pvals).ravel()
reject_fdr, p_adjusted = fdrcorrection(all_non_parametric_pvals_flat)
p_adjusted = p_adjusted.reshape(np.array(all_non_parametric_pvals).shape)
# %%
all_layer_pvals = []
for region in regions:
for layerA, layerB in itertools.combinations(layers, 2):
diff_betas_layerA = np.subtract(set2[layerA][region], set1[layerA][region]) # difference within layerA
diff_betas_layerB = np.subtract(set2[layerB][region], set1[layerB][region]) # difference within layerB
all_layer_pvals.append(non_parametric_permutation_ttest(diff_betas_layerA - diff_betas_layerB))
reject_fdr, p_layers_adjusted = fdrcorrection(all_layer_pvals)
all_non_parametric_pvals_for_layers = {}
idx = 0
for region in regions:
all_non_parametric_pvals_for_layers[region] = {}
for layerA, layerB in itertools.combinations(layers, 2):
all_non_parametric_pvals_for_layers[region]['_'.join([layerA, layerB])] = p_layers_adjusted[idx]
idx += 1
fig = plt.figure(figsize=(9, 9))
gs = gridspec.GridSpec(3, 3, height_ratios=[1, 1, 1], hspace=0.6, wspace=0.4)
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[0, 2])
ax4 = fig.add_subplot(gs[1, 0])
ax5 = fig.add_subplot(gs[1, 1])
ax6 = fig.add_subplot(gs[1, 2])
ax7 = fig.add_subplot(gs[2, 0])
ax8 = fig.add_subplot(gs[2, 1])
ax9 = fig.add_subplot(gs[2, 2])
for idx, axis in enumerate(fig.get_axes()): # regions
means1 = np.array([np.nanmean(set1[layer][regions[idx]]) for layer in layers])
stderr1 = np.array([np.nanstd(set1[layer][regions[idx]]) for layer in layers])/np.sqrt(11)
sns.lineplot(means1, ax=axis, c='b', linewidth=3)
axis.errorbar(x=range(len(means1)), y=means1, yerr=stderr1, fmt='o', capsize=2, c='b')
means2 = np.array([np.nanmean(set2[layer][regions[idx]]) for layer in layers])
stderr2 = np.array([np.nanstd(set2[layer][regions[idx]]) for layer in layers])/np.sqrt(11)
sns.lineplot(means2, ax=axis, c='r', linewidth=3)
axis.errorbar(x=range(len(means2)), y=means2, yerr=stderr2, fmt='o', capsize=2, c='r')
# for x, star in enumerate(tbx.viz.pvals2stars(p_adjusted[idx])):
# if star == 'ns': continue
# axis.text(x, 2.5, '*', ha='center', va='bottom', color='k', weight='bold', fontsize=20)
# # layer statistics
# for layer_comb in all_non_parametric_pvals_for_layers[regions[idx]]:
# lA, lB = layer_comb.split('_')
# xA, xB = int(lA[1])-1, int(lB[1])-1
# p_layer = all_non_parametric_pvals_for_layers[regions[idx]][layer_comb]
# star = tbx.viz.pvals2stars(p_layer)
# if star != 'ns':
# if xA==0 and xB==2:
# axis.plot([xA+.1,xB-.1], [2.2,2.2], c='k')
# axis.text((xB+xA)/2, 2.2, '*', ha='center', va='bottom', color='k', weight='bold', fontsize=20)
# axis.plot([xA+.1, xA+.1], [2.1,2.2], c='k')
# axis.plot([xB-.1, xB-.1], [2.1,2.2], c='k')
# else:
# axis.plot([xA+.1,xB-.1], [2,2], c='k')
# axis.text((xB+xA)/2, 2, '*', ha='center', va='bottom', color='k', weight='bold', fontsize=20)
# axis.plot([xA+.1, xA+.1], [2,1.9], c='k')
# axis.plot([xB-.1, xB-.1], [2,1.9], c='k')
axis.set_xticks([0,1,2], [r'$\textbf{D}$', r'$\textbf{M}$', r'$\textbf{S}$'])
axis.set_ylim([0.5, 2.5])
if idx in [0,3,6]:
axis.set_ylabel(r'$\mathbf{\beta\ (\%)}$', fontsize=20)
axis.set_title(r'$\textbf{'+regions[idx]+r'}$')
# make border thicker when there is an interaction
if regions[idx] in ['pSTG', 'pMTG', 'parsOrbitalis', 'parsTriangularis']:
for spine in axis.spines.values():
spine.set_linewidth(2)
spine.set_color('black')
# add layer interactions
if regions[idx] == 'pSTG':
axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5)
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5)
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5)
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5)
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5)
axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5)
if regions[idx] == 'pMTG':
axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5)
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5)
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5)
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5)
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5)
axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5)
if regions[idx] == 'parsOrbitalis':
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5)
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5)
axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5)
if regions[idx] == 'parsTriangularis':
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5)
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5)
axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5)
# if regions[idx] in ['PP', 'HG', 'PT', 'aSTG', 'TPOj']:
# axis.text(1.95, 2.35, '*', color='k', weight='bold', fontsize=20)
blue_line = mlines.Line2D([], [], color='b', linestyle='-', linewidth=2, label='Valid')
red_line = mlines.Line2D([], [], color='r', linestyle='-', linewidth=2, label='Invalid')
fig.legend(handles=[red_line, blue_line], loc='center', bbox_to_anchor=(.5, 0.02))
# Mean Activations of Congruent vs. Invalid Trials (11 Subjects, Averaged Hemisphere, Non-Parametric FDR Corrected
plt.savefig('./figures/univariate_laminar_cong_incong_nonparametric_fdr_friedman.png', dpi=300, bbox_inches='tight')
plt.show()
# %%
from scipy.stats import ttest_rel
ttest_rel(diff_betas_layerA, diff_betas_layerB)
# %%
region = 'parsOrbitalis'
means_diff = np.array([np.mean(np.array(set2[layer][region]) - np.array(set1[layer][region])) for layer in layers])
stderr_diff = np.array([np.std(np.array(set2[layer][region]) - np.array(set1[layer][region])) for layer in layers])/np.sqrt(11)
sns.lineplot(means_diff, ax=axis, c=palette[-2], linewidth=3)
axis.errorbar(x=range(len(means_diff)), y=means_diff, yerr=stderr_diff, fmt='o', capsize=2, c=palette[-2])
# friedman
diffMeas = np.array([set2[item][region] for item in ['L1', 'L2', 'L3']]) - np.array([set1[item][region] for item in ['L1', 'L2', 'L3']]).astype(np.float64)
friedman = friedmanchisquare(*diffMeas)
print(friedman)
# apply multcompare
eng = matlab.engine.start_matlab()
diffMeas_ml = matlab.double(diffMeas.tolist())
eng.workspace['diffMeas'] = diffMeas_ml
eng.eval("[p,tbl,stats] = friedman(diffMeas',1);", nargout=0)
c = eng.eval("multcompare(stats);", nargout=1)
c_py = np.array(c)
eng.quit()
print(pd.DataFrame(c_py[:,2:], index=['D-M', 'D-S', 'M-S'], columns=['Lower CI', 'Mean rank Diff', 'Upper CI', 'p-value']).to_latex())
# %%
sns.heatmap(diffMeas)
# %%
from pprint import pprint
pprint(diffMeas.T)
# %%
# do manual wilcoxon
region = 'parsOrbitalis'
from scipy.stats import wilcoxon
print(wilcoxon(diffMeas[0], diffMeas[1]))
print(wilcoxon(diffMeas[0], diffMeas[2]))
print(wilcoxon(diffMeas[1], diffMeas[2]))
# %%
regions_selected = ['pSTG', 'pMTG', 'parsOrbitalis', 'parsTriangularis']
palette = sns.color_palette('Grays')
all_layer_pvals = []
for region in regions_selected:
for layerA, layerB in itertools.combinations(layers, 2):
diff_betas_layerA = np.subtract(set2[layerA][region], set1[layerA][region]) # difference within layerA
diff_betas_layerB = np.subtract(set2[layerB][region], set1[layerB][region]) # difference within layerB
all_layer_pvals.append(non_parametric_permutation_ttest(diff_betas_layerA - diff_betas_layerB))
reject_fdr, p_layers_adjusted = fdrcorrection(all_layer_pvals)
all_non_parametric_pvals_for_layers = {}
idx = 0
for region in regions_selected:
all_non_parametric_pvals_for_layers[region] = {}
for layerA, layerB in itertools.combinations(layers, 2):
all_non_parametric_pvals_for_layers[region]['_'.join([layerA, layerB])] = p_layers_adjusted[idx]
idx += 1
fig = plt.figure(figsize=(6, 6))
gs = gridspec.GridSpec(2, 2, height_ratios=[1, 1], hspace=0.6, wspace=0.4)
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax4 = fig.add_subplot(gs[1, 0])
ax5 = fig.add_subplot(gs[1, 1])
for idx, axis in enumerate(fig.get_axes()): # regions_selected
means_diff = np.array([np.mean(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])
stderr_diff = np.array([np.std(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])/np.sqrt(11)
sns.lineplot(means_diff, ax=axis, c=palette[-2], linewidth=3)
axis.errorbar(x=range(len(means_diff)), y=means_diff, yerr=stderr_diff, fmt='o', capsize=2, c=palette[-2])
# friedman
diffMeas = np.array([set2[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']]) - np.array([set1[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']]).astype(np.float64)
friedman = friedmanchisquare(*diffMeas)
print(friedman)
# apply multcompare
eng = matlab.engine.start_matlab()
diffMeas_ml = matlab.double(diffMeas.tolist())
eng.workspace['diffMeas'] = diffMeas_ml
eng.eval("[p,tbl,stats] = friedman(diffMeas',1);", nargout=0)
c = eng.eval("multcompare(stats);", nargout=1)
c_py = np.array(c)
eng.quit()
print(pd.DataFrame(c_py[:,2:], index=['D-M', 'D-S', 'M-S'], columns=['Lower CI', 'Mean rank Diff', 'Upper CI', 'p-value']).to_latex())
if regions_selected[idx] == 'pSTG':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(0.85, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
# axis.axhline(0.25, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
if regions_selected[idx] == 'pMTG':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(0.85, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
# axis.axhline(0.25, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
if regions_selected[idx] == 'parsOrbitalis':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
axis.axhline(0.85, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
axis.axhline(0.85, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
if regions_selected[idx] == 'parsTriangularis':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(0.85, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
# axis.axhline(0.85, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
axis.set_xticks([0,1,2], [r'$\textbf{D}$', r'$\textbf{M}$', r'$\textbf{S}$'])
axis.set_ylim([0.0, 1.])
if idx in [0,2]:
axis.set_ylabel(r'$\mathbf{\Delta\beta\ (\%)}$', fontsize=20)
axis.set_title(r'$\textbf{'+regions_selected[idx]+r'}$')
# blue_line = mlines.Line2D([], [], color='b', linestyle='-', linewidth=2, label='Valid')
# red_line = mlines.Line2D([], [], color='r', linestyle='-', linewidth=2, label='Invalid')
# fig.legend(handles=[red_line, blue_line], loc='center', bbox_to_anchor=(.5, 0.02))
# Mean Activations of Congruent vs. Invalid Trials (11 Subjects, Averaged Hemisphere, Non-Parametric FDR Corrected
plt.savefig('./figures/univariate_laminar_cong_incong_nonparametric_fdr_friedman_posthoc.png', dpi=300, bbox_inches='tight')
plt.show()
# %%
diffs = []
for idx, axis in enumerate(fig.get_axes()): # regions_selected
means_diff = np.array([np.mean(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])
stderr_diff = np.array([np.std(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])/np.sqrt(11)
sns.lineplot(means_diff, ax=axis, c=palette[-2], linewidth=3)
axis.errorbar(x=range(len(means_diff)), y=means_diff, yerr=stderr_diff, fmt='o', capsize=2, c=palette[-2])
# friedman
diffMeas = np.array([set2[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']]) - np.array([set1[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']]).astype(np.float64)
print(np.mean(diffMeas, axis=1), '\t', regions_selected[idx])
diffs.append(np.mean(diffMeas, axis=1))
diffs = np.array(diffs)
np.min(diffs), np.max(diffs)
# %%
# (0.16910681941292502, 0.6572449369864031
# (-1.1287472214211116, 0.6597942059690302
# (-0.49109694124622777, -0.08562372760339217
scales = [0.16910681941292502, 0.6572449369864031, -1.1287472214211116, 0.6597942059690302, -0.49109694124622777, -0.08562372760339217]
mincol, maxcol = np.min(scales), np.max(scales)
mincol, maxcol
maxcol = 0.4
mincol = 0.1
# %%
from copy import deepcopy
smp_head, smp_data = bvbabel.smp.read_smp('/mnt/hdd2/associative_learning/paper_draft_v1/paper_draft_ALPHA/templates/group_patches_LH_template.smp')
smp_data[smp_data<6] = 0
smp_data[smp_data>0] = 6
# %%
all_regs
# %%
all_regs = [entry['Name'] for entry in smp_head['Map']]
dep = ['pSTG', 'pMTG', 'parsOrbitalis', 'parsTriangularis']
mid = ['pSTG', 'pMTG', 'parsOrbitalis', 'parsTriangularis']
dep = ['pSTG', 'pMTG', 'parsOrbitalis', 'parsTriangularis']
# Superficial
_, indcs, _ = np.intersect1d(all_regs, dep, return_indices=True)
smp_head_cp = deepcopy(smp_head)
smp_head_cp['Nr maps'] = len(indcs)
smp_head_cp['Map'] = list(np.array(smp_head_cp['Map'])[indcs])
data = deepcopy(smp_data[:, indcs])
for idx, tval in enumerate(diffs[:,0]):
data[:, idx][data[:, idx]==6] = tval
smp_head_cp['Map'][idx]['Threshold min'] = mincol
smp_head_cp['Map'][idx]['Threshold max'] = maxcol
print(smp_head_cp['Map'][idx]['Threshold min'], smp_head_cp['Map'][idx]['Threshold max'])
bvbabel.smp.write_smp('/mnt/hdd2/associative_learning/paper_draft_v1/figures/group_patches_LH_univariate_valid_invalid_dep.smp', smp_head_cp, data)
# Middle
_, indcs, _ = np.intersect1d(all_regs, mid, return_indices=True)
smp_head_cp = deepcopy(smp_head)
smp_head_cp['Nr maps'] = len(indcs)
smp_head_cp['Map'] = list(np.array(smp_head_cp['Map'])[indcs])
data = deepcopy(smp_data[:, indcs])
for idx, tval in enumerate(diffs[:,1]):
data[:, idx][data[:, idx]==6] = tval
smp_head_cp['Map'][idx]['Threshold min'] = mincol
smp_head_cp['Map'][idx]['Threshold max'] = maxcol
print(smp_head_cp['Map'][idx]['Threshold min'], smp_head_cp['Map'][idx]['Threshold max'])
bvbabel.smp.write_smp('/mnt/hdd2/associative_learning/paper_draft_v1/figures/group_patches_LH_univariate_valid_invalid_mid.smp', smp_head_cp, data)
# Deep
_, indcs, _ = np.intersect1d(all_regs, sup, return_indices=True)
smp_head_cp = deepcopy(smp_head)
smp_head_cp['Nr maps'] = len(indcs)
smp_head_cp['Map'] = list(np.array(smp_head_cp['Map'])[indcs])
data = deepcopy(smp_data[:, indcs])
for idx, tval in enumerate(diffs[:,2]):
data[:, idx][data[:, idx]==6] = tval
smp_head_cp['Map'][idx]['Threshold min'] = mincol
smp_head_cp['Map'][idx]['Threshold max'] = maxcol
print(smp_head_cp['Map'][idx]['Threshold min'], smp_head_cp['Map'][idx]['Threshold max'])
bvbabel.smp.write_smp('/mnt/hdd2/associative_learning/paper_draft_v1/figures/group_patches_LH_univariate_valid_invalid_sup.smp', smp_head_cp, data)
# %%
import colorsys
import matplotlib.pyplot as plt
def hex_to_rgb(hex_color):
"""Convert hex color to an RGB tuple."""
return tuple(int(hex_color[i:i+2], 16) for i in (0, 2, 4))
def generate_gradient_hsv(colors, steps=100):
"""
Generate a gradient passing through multiple RGB colors using HSV interpolation.
Parameters:
colors (list of tuples): List of RGB colors [(R, G, B), ...] with values between 0-255.
steps (int): The total number of gradient steps.
Returns:
List of RGB tuples forming the gradient.
"""
num_segments = len(colors) - 1
steps_per_segment = steps // num_segments
gradient = []
for i in range(num_segments):
color1 = colors[i]
color2 = colors[i + 1]
hsv1 = colorsys.rgb_to_hsv(color1[0] / 255.0, color1[1] / 255.0, color1[2] / 255.0)
hsv2 = colorsys.rgb_to_hsv(color2[0] / 255.0, color2[1] / 255.0, color2[2] / 255.0)
if abs(hsv2[0] - hsv1[0]) > 0.5:
if hsv1[0] > hsv2[0]:
hsv2 = (hsv2[0] + 1, hsv2[1], hsv2[2])
else:
hsv1 = (hsv1[0] + 1, hsv1[1], hsv1[2])
for j in range(steps_per_segment):
t = j / steps_per_segment
hue = (hsv1[0] + (hsv2[0] - hsv1[0]) * t) % 1
saturation = hsv1[1] + (hsv2[1] - hsv1[1]) * t
value = hsv1[2] + (hsv2[2] - hsv1[2]) * t
r, g, b = colorsys.hsv_to_rgb(hue, saturation, value)
gradient.append((int(r * 255), int(g * 255), int(b * 255)))
return gradient
def display_gradient(gradient, fn):
"""Display the gradient using Matplotlib."""
fig, ax = plt.subplots(figsize=(10, 1))
ax.imshow([gradient], extent=[0, len(gradient), 0, 1], aspect='auto')
ax.set_xticks([])
ax.set_yticks([])
plt.savefig(fn, dpi=300)
plt.show()
# %%
# hex_colors = np.array(['087bed','086ee0','0861d3','0854c7','0848ba','083aae','082ea1','082194','081588','08087b','b50808','bb1b08','c12e08','c84208','ce5408','d46708','db7b08','e18e08','e7a108','edb508'])#[::-1]
hex_colors = np.array(['b50808','bb1b08','c12e08','c84208','ce5408','d46707','db7b08','e18e08','edb507','edb507'])#[::-1]
colors = [hex_to_rgb(hex_color) for hex_color in hex_colors]
steps = 2000
gradient = generate_gradient_hsv(colors, steps)
display_gradient(gradient, '/mnt/hdd2/associative_learning/paper_draft_v1/figures/gradient_valid_invalid.png')
# %% [markdown]
# ---
#
# # <center> $\Large\textbf{Congruent vs. Omitted}$ </center>
#
# **Test the effect of prediction error.**
# %%
regex1 = '_cong$'
regex2 = '^(?!random).*_omis$'
layers = ['L1', 'L2', 'L3']
tstats = {}; pvals = {}; set1 = {};set2 = {}
for layer_name, layer_id in zip(['Deep', 'Middle', 'Superficial'], ['L1', 'L2', 'L3']):
tstats[layer_id] = []; pvals[layer_id] = []; df_set1 = {}; df_set2 = {}
for region in regions:
df_set1[region] = []
df_set2[region] = []
for subject in subjects:
stimuli = tbx.utils.extract_stimuli(subject)
betas_l = tbx.utils.extract_condition_betas(stimuli=stimuli,
betas=dataset[subject]['LH'][region]['betas_selected'][layer_id].values.T,
stimulus_type='sound')
betas_r = tbx.utils.extract_condition_betas(stimuli=stimuli,
betas=dataset[subject]['RH'][region]['betas_selected'][layer_id].values.T,
stimulus_type='sound')
df_set1[region].append(np.mean([
np.nanmean(betas_l.filter(regex=regex1).values),
np.nanmean(betas_r.filter(regex=regex1).values)]))
df_set2[region].append(np.mean([
np.nanmean(betas_l.filter(regex=regex2).values),
np.nanmean(betas_r.filter(regex=regex2).values)]))
set1[layer_id] = df_set1
set2[layer_id] = df_set2
all_non_parametric_pvals = np.array([[non_parametric_permutation_ttest(np.subtract(set2[layer][region], set1[layer][region])) for layer in layers] for region in regions])
all_non_parametric_pvals_flat = np.array(all_non_parametric_pvals).ravel()
reject_fdr, p_adjusted = fdrcorrection(all_non_parametric_pvals_flat)
p_adjusted = p_adjusted.reshape(np.array(all_non_parametric_pvals).shape)
# %%
all_layer_pvals = []
for region in regions:
for layerA, layerB in itertools.combinations(layers, 2):
diff_betas_layerA = np.subtract(set2[layerA][region], set1[layerA][region]) # difference within layerA
diff_betas_layerB = np.subtract(set2[layerB][region], set1[layerB][region]) # difference within layerB
all_layer_pvals.append(non_parametric_permutation_ttest(diff_betas_layerA - diff_betas_layerB))
reject_fdr, p_layers_adjusted = fdrcorrection(all_layer_pvals)
all_non_parametric_pvals_for_layers = {}
idx = 0
for region in regions:
all_non_parametric_pvals_for_layers[region] = {}
for layerA, layerB in itertools.combinations(layers, 2):
all_non_parametric_pvals_for_layers[region]['_'.join([layerA, layerB])] = p_layers_adjusted[idx]
idx += 1
fig = plt.figure(figsize=(9, 9))
gs = gridspec.GridSpec(3, 3, height_ratios=[1, 1, 1], hspace=0.6, wspace=0.4)
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[0, 2])
ax4 = fig.add_subplot(gs[1, 0])
ax5 = fig.add_subplot(gs[1, 1])
ax6 = fig.add_subplot(gs[1, 2])
ax7 = fig.add_subplot(gs[2, 0])
ax8 = fig.add_subplot(gs[2, 1])
ax9 = fig.add_subplot(gs[2, 2])
for idx, axis in enumerate(fig.get_axes()): # regions
means1 = np.array([np.nanmean(set1[layer][regions[idx]]) for layer in layers])
stderr1 = np.array([np.nanstd(set1[layer][regions[idx]]) for layer in layers])/np.sqrt(11)
sns.lineplot(means1, ax=axis, c='b', linewidth=3)
axis.errorbar(x=range(len(means1)), y=means1, yerr=stderr1, fmt='o', capsize=2, c='b')
means2 = np.array([np.nanmean(set2[layer][regions[idx]]) for layer in layers])
stderr2 = np.array([np.nanstd(set2[layer][regions[idx]]) for layer in layers])/np.sqrt(11)
sns.lineplot(means2, ax=axis, c='r', linewidth=3)
axis.errorbar(x=range(len(means2)), y=means2, yerr=stderr2, fmt='o', capsize=2, c='r')
axis.set_xticks([0,1,2], [r'$\textbf{D}$', r'$\textbf{M}$', r'$\textbf{S}$'])
axis.set_ylim([0.5, 2.5])
if idx in [0,3,6]:
axis.set_ylabel(r'$\mathbf{\beta\ (\%)}$', fontsize=20)
axis.set_title(r'$\textbf{'+regions[idx]+r'}$')
# make border thicker when there is an interaction
if regions[idx] in ['PP', 'HG', 'PT', 'aSTG', 'pSTG', 'parsOrbitalis', 'parsTriangularis', 'TPOj']:
for spine in axis.spines.values():
spine.set_linewidth(2)
spine.set_color('black')
# add layer interactions
if regions[idx] == 'PP':
# axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
# axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'HG':
axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
# axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'PT':
axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
# axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
# axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
# axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'aSTG':
# axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
# axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
# axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
# axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
# axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'pSTG':
axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
# axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'parsOrbitalis':
# axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
# axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
# axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'parsTriangularis':
# axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
# axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'TPOj':
# axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
# axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
# axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
# axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
# axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
# if regions[idx] in ['pMTG']:
# axis.text(1.95, 2.35, '*', color='k', weight='bold', fontsize=20)
blue_line = mlines.Line2D([], [], color='b', linestyle='-', linewidth=2, label='Valid')
red_line = mlines.Line2D([], [], color='r', linestyle='-', linewidth=2, label='Omitted (Predictable)')
fig.legend(handles=[red_line, blue_line], loc='center', bbox_to_anchor=(.5, 0.02))
# Mean Activations of Congruent vs. Invalid Trials (11 Subjects, Averaged Hemisphere, Non-Parametric FDR Corrected
plt.savefig('./figures/univariate_laminar_cong_omis_nonparametric_fdr_friedman.png', dpi=300, bbox_inches='tight')
plt.show()
# %%
regions_selected = ['PP', 'HG', 'PT', 'aSTG', 'pSTG', 'parsOrbitalis', 'parsTriangularis', 'TPOj']
palette = sns.color_palette('Grays')
all_layer_pvals = []
for region in regions_selected:
for layerA, layerB in itertools.combinations(layers, 2):
diff_betas_layerA = np.subtract(set2[layerA][region], set1[layerA][region]) # difference within layerA
diff_betas_layerB = np.subtract(set2[layerB][region], set1[layerB][region]) # difference within layerB
all_layer_pvals.append(non_parametric_permutation_ttest(diff_betas_layerA - diff_betas_layerB))
reject_fdr, p_layers_adjusted = fdrcorrection(all_layer_pvals)
all_non_parametric_pvals_for_layers = {}
idx = 0
for region in regions_selected:
all_non_parametric_pvals_for_layers[region] = {}
for layerA, layerB in itertools.combinations(layers, 2):
all_non_parametric_pvals_for_layers[region]['_'.join([layerA, layerB])] = p_layers_adjusted[idx]
idx += 1
fig = plt.figure(figsize=(12, 6))
gs = gridspec.GridSpec(2, 4, height_ratios=[1, 1], hspace=0.6, wspace=0.4)
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[0, 2])
ax4 = fig.add_subplot(gs[0, 3])
ax5 = fig.add_subplot(gs[1, 0])
ax6 = fig.add_subplot(gs[1, 1])
ax7 = fig.add_subplot(gs[1, 2])
ax8 = fig.add_subplot(gs[1, 3])
for idx, axis in enumerate(fig.get_axes()): # regions_selected
means_diff = np.array([np.nanmean(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])
stderr_diff = np.array([np.nanstd(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])/np.sqrt(11)
sns.lineplot(means_diff, ax=axis, c=palette[-2], linewidth=3)
axis.errorbar(x=range(len(means_diff)), y=means_diff, yerr=stderr_diff, fmt='o', capsize=2, c=palette[-2])
friedman
diffMeas = np.array([set2[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']]) - np.array([set1[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']])
friedman = friedmanchisquare(*diffMeas)
print(friedman)
# apply multcompare
eng = matlab.engine.start_matlab()
diffMeas_ml = matlab.double(diffMeas.tolist())
eng.workspace['diffMeas'] = diffMeas_ml
eng.eval("[p,tbl,stats] = friedman(diffMeas',1);", nargout=0)
c = eng.eval("multcompare(stats);", nargout=1)
c_py = np.array(c)
eng.quit()
print(pd.DataFrame(c_py[:,2:], index=['D-M', 'D-S', 'M-S'], columns=['Lower CI', 'Mean rank Diff', 'Upper CI', 'p-value']).to_latex())
if regions_selected[idx] == 'PP':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(0.8, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
# axis.axhline(0.8, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
if regions_selected[idx] == 'HG':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
axis.axhline(0.8, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
# axis.axhline(0.8, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
if regions_selected[idx] == 'PT':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(0.8, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
axis.axhline(0.8, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
if regions_selected[idx] == 'aSTG':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(0.8, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
axis.axhline(0.8, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
if regions_selected[idx] == 'pSTG':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(0.8, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
axis.axhline(0.8, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
if regions_selected[idx] == 'parsOrbitalis':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(0.8, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
axis.axhline(0.8, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
if regions_selected[idx] == 'parsTriangularis':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(0.8, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
# axis.axhline(0.8, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
if regions_selected[idx] == 'TPOj':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(0.8, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
# axis.axhline(0.8, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
axis.set_xticks([0,1,2], [r'$\textbf{D}$', r'$\textbf{M}$', r'$\textbf{S}$'])
axis.set_ylim([-1.5, 1.0])
if idx in [0,4]:
axis.set_ylabel(r'$\mathbf{\Delta\beta\ (\%)}$', fontsize=20)
axis.set_title(r'$\textbf{'+regions_selected[idx]+r'}$')
# blue_line = mlines.Line2D([], [], color='b', linestyle='-', linewidth=2, label='Valid')
# red_line = mlines.Line2D([], [], color='r', linestyle='-', linewidth=2, label='Invalid')
# fig.legend(handles=[red_line, blue_line], loc='center', bbox_to_anchor=(.5, 0.02))
# Mean Activations of Congruent vs. Invalid Trials (11 Subjects, Averaged Hemisphere, Non-Parametric FDR Corrected
plt.savefig('./figures/univariate_laminar_cong_omis_nonparametric_fdr_friedman_posthoc.png', dpi=300, bbox_inches='tight')
plt.show()
# %%
diffs = []
for idx, axis in enumerate(fig.get_axes()): # regions_selected
means_diff = np.array([np.mean(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])
stderr_diff = np.array([np.std(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])/np.sqrt(11)
sns.lineplot(means_diff, ax=axis, c=palette[-2], linewidth=3)
axis.errorbar(x=range(len(means_diff)), y=means_diff, yerr=stderr_diff, fmt='o', capsize=2, c=palette[-2])
# friedman
diffMeas = np.array([set2[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']]) - np.array([set1[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']]).astype(np.float64)
print(np.mean(diffMeas, axis=1), '\t', regions_selected[idx])
diffs.append(np.mean(diffMeas, axis=1))
diffs = np.array(diffs)
print(diffs.shape)
np.min(diffs), np.max(diffs)
# %%
mincol = 0.1
maxcol = .4
# %%
all_regs = [entry['Name'] for entry in smp_head['Map']]
indcs = np.array([0,1,2,3,4,8,7,6])
selected = ['HG', 'PP', 'PT', 'aSTG', 'pSTG', 'parsOrbitalis', 'parsTriangularis', 'TPOj']
# Superficial
# _, indcs, _ = np.intersect1d(all_regs, selected, return_indices=True)
smp_head_cp = deepcopy(smp_head)
smp_head_cp['Nr maps'] = len(indcs)
smp_head_cp['Map'] = list(np.array(smp_head_cp['Map'])[indcs])
data = deepcopy(smp_data[:, indcs])
for idx, tval in enumerate(diffs[:,0]):
data[:, idx][data[:, idx]==6] = tval
smp_head_cp['Map'][idx]['Threshold min'] = mincol
smp_head_cp['Map'][idx]['Threshold max'] = maxcol
print(smp_head_cp['Map'][idx]['Threshold min'], smp_head_cp['Map'][idx]['Threshold max'])
bvbabel.smp.write_smp('/mnt/hdd2/associative_learning/paper_draft_v1/figures/group_patches_LH_univariate_valid_predomit_dep.smp', smp_head_cp, data)
# Middle
# _, indcs, _ = np.intersect1d(all_regs, selected, return_indices=True)
smp_head_cp = deepcopy(smp_head)
smp_head_cp['Nr maps'] = len(indcs)
smp_head_cp['Map'] = list(np.array(smp_head_cp['Map'])[indcs])
data = deepcopy(smp_data[:, indcs])
for idx, tval in enumerate(diffs[:,1]):
data[:, idx][data[:, idx]==6] = tval
smp_head_cp['Map'][idx]['Threshold min'] = mincol
smp_head_cp['Map'][idx]['Threshold max'] = maxcol
print(smp_head_cp['Map'][idx]['Threshold min'], smp_head_cp['Map'][idx]['Threshold max'])
bvbabel.smp.write_smp('/mnt/hdd2/associative_learning/paper_draft_v1/figures/group_patches_LH_univariate_valid_predomit_mid.smp', smp_head_cp, data)
# Deep
# _, indcs, _ = np.intersect1d(all_regs, selected, return_indices=True)
smp_head_cp = deepcopy(smp_head)
smp_head_cp['Nr maps'] = len(indcs)
smp_head_cp['Map'] = list(np.array(smp_head_cp['Map'])[indcs])
data = deepcopy(smp_data[:, indcs])
for idx, tval in enumerate(diffs[:,2]):
data[:, idx][data[:, idx]==6] = tval
smp_head_cp['Map'][idx]['Threshold min'] = mincol
smp_head_cp['Map'][idx]['Threshold max'] = maxcol
print(smp_head_cp['Map'][idx]['Threshold min'], smp_head_cp['Map'][idx]['Threshold max'])
bvbabel.smp.write_smp('/mnt/hdd2/associative_learning/paper_draft_v1/figures/group_patches_LH_univariate_valid_predomit_sup.smp', smp_head_cp, data)
# %%
def display_gradient(gradient, fn, figsize=(10, 1)):
"""Display the gradient using Matplotlib."""
fig, ax = plt.subplots(figsize=figsize)
ax.imshow([gradient], extent=[0, len(gradient), 0, 1], aspect='auto')
ax.set_xticks([])
ax.set_yticks([])
plt.savefig(fn, dpi=300)
plt.show()
# %%
steps = 2000
fig, ax = plt.subplots(1,2, figsize=(10,1))
plt.subplots_adjust(wspace=0.01)
hex_colors = np.array(['b50808','bb1b08','c12e08','c84208','ce5408','d46708','db7b08','e18e08','e7a108','edb508'])#[::-1]
colors = [hex_to_rgb(hex_color) for hex_color in hex_colors]
gradient = generate_gradient_hsv(colors, steps)
fn = '/mnt/hdd2/associative_learning/paper_draft_v1/figures/gradient_valid_predomit.png'
ax[1].imshow([gradient], extent=[0, len(gradient), 0, 1], aspect='auto')
ax[1].set_xticks([])
ax[1].set_yticks([])
hex_colors = np.array(['087bed','086ee0','0861d3','0854c7','0848ba','083aae','082ea1','082194','081588','08087b'])
colors = [hex_to_rgb(hex_color) for hex_color in hex_colors]
gradient = generate_gradient_hsv(colors, steps)
ax[0].imshow([gradient], extent=[0, len(gradient), 0, 1], aspect='auto')
ax[0].set_xticks([])
ax[0].set_yticks([])
plt.savefig(fn, dpi=300)
plt.show()
# %% [markdown]
# ---
#
# # <center> $\Large\textbf{Congruent vs. Random}$ </center>
#
# **Test the effect of prediction.**
# %%
regex1 = '_cong$'
regex2 = 'MANUAL'
set1 = {};set2 = {}
for layer_name, layer_id in zip(['Deep', 'Middle', 'Superficial'], ['L1', 'L2', 'L3']):
tstats[layer_id] = []; pvals[layer_id] = []; df_set1 = {}; df_set2 = {}
for region in regions:
df_set1[region] = []
df_set2[region] = []
for subject in subjects:
stimuli = tbx.utils.extract_stimuli(subject)
betas_l = tbx.utils.extract_condition_betas(stimuli=stimuli,
betas=dataset[subject]['LH'][region]['betas_selected'][layer_id].values.T,
stimulus_type='sound')
betas_r = tbx.utils.extract_condition_betas(stimuli=stimuli,
betas=dataset[subject]['RH'][region]['betas_selected'][layer_id].values.T,
stimulus_type='sound')
df_set1[region].append(np.mean([
np.nanmean(betas_l.filter(regex=regex1).values),
np.nanmean(betas_r.filter(regex=regex1).values)]))
df_set2[region].append(np.mean([
np.nanmean(betas_l.filter(regex='^random').filter(regex='^(?!.*_omis$)').values),
np.nanmean(betas_r.filter(regex='^random').filter(regex='^(?!.*_omis$)').values)]))
set1[layer_id] = df_set1
set2[layer_id] = df_set2
all_non_parametric_pvals = np.array([[non_parametric_permutation_ttest(np.subtract(set2[layer][region], set1[layer][region])) for layer in layers] for region in regions])
all_non_parametric_pvals_flat = np.array(all_non_parametric_pvals).ravel()
reject_fdr, p_adjusted = fdrcorrection(all_non_parametric_pvals_flat)
p_adjusted = p_adjusted.reshape(np.array(all_non_parametric_pvals).shape)
# %%
all_layer_pvals = []
for region in regions:
for layerA, layerB in itertools.combinations(layers, 2):
diff_betas_layerA = np.subtract(set2[layerA][region], set1[layerA][region]) # difference within layerA
diff_betas_layerB = np.subtract(set2[layerB][region], set1[layerB][region]) # difference within layerB
all_layer_pvals.append(non_parametric_permutation_ttest(diff_betas_layerA - diff_betas_layerB))
reject_fdr, p_layers_adjusted = fdrcorrection(all_layer_pvals)
all_non_parametric_pvals_for_layers = {}
idx = 0
for region in regions:
all_non_parametric_pvals_for_layers[region] = {}
for layerA, layerB in itertools.combinations(layers, 2):
all_non_parametric_pvals_for_layers[region]['_'.join([layerA, layerB])] = p_layers_adjusted[idx]
idx += 1
fig = plt.figure(figsize=(9, 9))
gs = gridspec.GridSpec(3, 3, height_ratios=[1, 1, 1], hspace=0.6, wspace=0.4)
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[0, 2])
ax4 = fig.add_subplot(gs[1, 0])
ax5 = fig.add_subplot(gs[1, 1])
ax6 = fig.add_subplot(gs[1, 2])
ax7 = fig.add_subplot(gs[2, 0])
ax8 = fig.add_subplot(gs[2, 1])
ax9 = fig.add_subplot(gs[2, 2])
for idx, axis in enumerate(fig.get_axes()): # regions
means1 = np.array([np.nanmean(set1[layer][regions[idx]]) for layer in layers])
stderr1 = np.array([np.nanstd(set1[layer][regions[idx]]) for layer in layers])/np.sqrt(11)
sns.lineplot(means1, ax=axis, c='b', linewidth=3)
axis.errorbar(x=range(len(means1)), y=means1, yerr=stderr1, fmt='o', capsize=2, c='b')
means2 = np.array([np.nanmean(set2[layer][regions[idx]]) for layer in layers])
stderr2 = np.array([np.nanstd(set2[layer][regions[idx]]) for layer in layers])/np.sqrt(11)
sns.lineplot(means2, ax=axis, c='r', linewidth=3)
axis.errorbar(x=range(len(means2)), y=means2, yerr=stderr2, fmt='o', capsize=2, c='r')
axis.set_xticks([0,1,2], [r'$\textbf{D}$', r'$\textbf{M}$', r'$\textbf{S}$'])
axis.set_ylim([0.5, 2.5])
if idx in [0,3,6]:
axis.set_ylabel(r'$\mathbf{\beta\ (\%)}$', fontsize=20)
axis.set_title(r'$\textbf{'+regions[idx]+r'}$')
# make border thicker when there is an interaction
# if regions[idx] in ['PP', 'HG', 'PT', 'aSTG', 'pSTG', 'parsOrbitalis', 'parsTriangularis', 'TPOj']:
# for spine in axis.spines.values():
# spine.set_linewidth(2)
# spine.set_color('black')
# add layer interactions
# if regions[idx] == 'PP':
# axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
# axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
# axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
# axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
# axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
# axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] in ['PP']:
axis.text(1.95, 2.35, '*', color='k', weight='bold', fontsize=20)
blue_line = mlines.Line2D([], [], color='b', linestyle='-', linewidth=2, label='Valid')
red_line = mlines.Line2D([], [], color='r', linestyle='-', linewidth=2, label='Unpredictable')
fig.legend(handles=[red_line, blue_line], loc='center', bbox_to_anchor=(.5, 0.02))
# Mean Activations of Congruent vs. Invalid Trials (11 Subjects, Averaged Hemisphere, Non-Parametric FDR Corrected
plt.savefig('./figures/univariate_laminar_cong_rand_nonparametric_fdr_friedman.png', dpi=300, bbox_inches='tight')
plt.show()
# %%
diffs = []
for idx, axis in enumerate(fig.get_axes()): # regions_selected
means_diff = np.array([np.mean(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])
stderr_diff = np.array([np.std(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])/np.sqrt(11)
sns.lineplot(means_diff, ax=axis, c=palette[-2], linewidth=3)
axis.errorbar(x=range(len(means_diff)), y=means_diff, yerr=stderr_diff, fmt='o', capsize=2, c=palette[-2])
# friedman
diffMeas = np.array([set2[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']]) - np.array([set1[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']]).astype(np.float64)
print(np.mean(diffMeas, axis=1), '\t', regions_selected[idx])
diffs.append(np.mean(diffMeas, axis=1))
diffs = np.array(diffs)
np.min(diffs), np.max(diffs)
# %% [markdown]
# ---
#
# # <center> $\Large\textbf{Omitted Predictable vs. Omitted Random}$ </center>
#
# **Test the effect of prediction.**
# %%
regex1 = '^(?!random).*_omis$'
regex2 = '^random.*_omis$'
set1 = {};set2 = {}
for layer_name, layer_id in zip(['Deep', 'Middle', 'Superficial'], ['L1', 'L2', 'L3']):
tstats[layer_id] = []; pvals[layer_id] = []; df_set1 = {}; df_set2 = {}
for region in regions:
df_set1[region] = []
df_set2[region] = []
for subject in subjects:
stimuli = tbx.utils.extract_stimuli(subject)
betas_l = tbx.utils.extract_condition_betas(stimuli=stimuli,
betas=dataset[subject]['LH'][region]['betas_selected'][layer_id].values.T,
stimulus_type='sound')
betas_r = tbx.utils.extract_condition_betas(stimuli=stimuli,
betas=dataset[subject]['RH'][region]['betas_selected'][layer_id].values.T,
stimulus_type='sound')
df_set1[region].append(np.mean([
np.nanmean(betas_l.filter(regex=regex1).values),
np.nanmean(betas_r.filter(regex=regex1).values)]))
df_set2[region].append(np.mean([
np.nanmean(betas_l.filter(regex=regex2).values),
np.nanmean(betas_r.filter(regex=regex2).values)]))
set1[layer_id] = df_set1
set2[layer_id] = df_set2
all_non_parametric_pvals = np.array([[non_parametric_permutation_ttest(np.subtract(set2[layer][region], set1[layer][region])) for layer in layers] for region in regions])
all_non_parametric_pvals_flat = np.array(all_non_parametric_pvals).ravel()
reject_fdr, p_adjusted = fdrcorrection(all_non_parametric_pvals_flat)
p_adjusted = p_adjusted.reshape(np.array(all_non_parametric_pvals).shape)
# %%
all_layer_pvals = []
for region in regions:
for layerA, layerB in itertools.combinations(layers, 2):
diff_betas_layerA = np.subtract(set2[layerA][region], set1[layerA][region]) # difference within layerA
diff_betas_layerB = np.subtract(set2[layerB][region], set1[layerB][region]) # difference within layerB
all_layer_pvals.append(non_parametric_permutation_ttest(diff_betas_layerA - diff_betas_layerB))
reject_fdr, p_layers_adjusted = fdrcorrection(all_layer_pvals)
all_non_parametric_pvals_for_layers = {}
idx = 0
for region in regions:
all_non_parametric_pvals_for_layers[region] = {}
for layerA, layerB in itertools.combinations(layers, 2):
all_non_parametric_pvals_for_layers[region]['_'.join([layerA, layerB])] = p_layers_adjusted[idx]
idx += 1
fig = plt.figure(figsize=(9, 9))
gs = gridspec.GridSpec(3, 3, height_ratios=[1, 1, 1], hspace=0.6, wspace=0.4)
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[0, 2])
ax4 = fig.add_subplot(gs[1, 0])
ax5 = fig.add_subplot(gs[1, 1])
ax6 = fig.add_subplot(gs[1, 2])
ax7 = fig.add_subplot(gs[2, 0])
ax8 = fig.add_subplot(gs[2, 1])
ax9 = fig.add_subplot(gs[2, 2])
for idx, axis in enumerate(fig.get_axes()): # regions
means1 = np.array([np.nanmean(set1[layer][regions[idx]]) for layer in layers])
stderr1 = np.array([np.nanstd(set1[layer][regions[idx]]) for layer in layers])/np.sqrt(11)
sns.lineplot(means1, ax=axis, c='b', linewidth=3)
axis.errorbar(x=range(len(means1)), y=means1, yerr=stderr1, fmt='o', capsize=2, c='b')
means2 = np.array([np.nanmean(set2[layer][regions[idx]]) for layer in layers])
stderr2 = np.array([np.nanstd(set2[layer][regions[idx]]) for layer in layers])/np.sqrt(11)
sns.lineplot(means2, ax=axis, c='r', linewidth=3)
axis.errorbar(x=range(len(means2)), y=means2, yerr=stderr2, fmt='o', capsize=2, c='r')
axis.set_xticks([0,1,2], [r'$\textbf{D}$', r'$\textbf{M}$', r'$\textbf{S}$'])
axis.set_ylim([0., 2.5])
if idx in [0,3,6]:
axis.set_ylabel(r'$\mathbf{\beta\ (\%)}$', fontsize=20)
axis.set_title(r'$\textbf{'+regions[idx]+r'}$')
# make border thicker when there is an interaction
if regions[idx] in ['HG', 'pSTG', 'pMTG']:
for spine in axis.spines.values():
spine.set_linewidth(2)
spine.set_color('black')
# add layer interactions
if regions[idx] == 'HG':
axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
# axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'pSTG':
axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
# axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
# axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
# axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'pMTG':
axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
# axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
# axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
# axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
blue_line = mlines.Line2D([], [], color='b', linestyle='-', linewidth=2, label='Omitted (Predictable)')
red_line = mlines.Line2D([], [], color='r', linestyle='-', linewidth=2, label='Omitted (Unpredictable)')
fig.legend(handles=[red_line, blue_line], loc='center', bbox_to_anchor=(.5, 0.))
# Mean Activations of Congruent vs. Invalid Trials (11 Subjects, Averaged Hemisphere, Non-Parametric FDR Corrected
# plt.savefig('./results/univariate_laminar_congomit_randomit_nonparametric_fdr.png', dpi=300, bbox_inches='tight')
plt.savefig('./figures/univariate_laminar_congomis_randomis_nonparametric_fdr_friedman.png', dpi=300, bbox_inches='tight')
plt.show()
# %%
regions_selected = ['HG', 'pSTG', 'pMTG']
palette = sns.color_palette('Grays')
all_layer_pvals = []
for region in regions_selected:
for layerA, layerB in itertools.combinations(layers, 2):
diff_betas_layerA = np.subtract(set2[layerA][region], set1[layerA][region]) # difference within layerA
diff_betas_layerB = np.subtract(set2[layerB][region], set1[layerB][region]) # difference within layerB
all_layer_pvals.append(non_parametric_permutation_ttest(diff_betas_layerA - diff_betas_layerB))
reject_fdr, p_layers_adjusted = fdrcorrection(all_layer_pvals)
all_non_parametric_pvals_for_layers = {}
idx = 0
for region in regions_selected:
all_non_parametric_pvals_for_layers[region] = {}
for layerA, layerB in itertools.combinations(layers, 2):
all_non_parametric_pvals_for_layers[region]['_'.join([layerA, layerB])] = p_layers_adjusted[idx]
idx += 1
fig = plt.figure(figsize=(9, 6))
gs = gridspec.GridSpec(2, 3, height_ratios=[1, 1], hspace=0.6, wspace=0.4)
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[0, 2])
for idx, axis in enumerate(fig.get_axes()): # regions_selected
# means_diff = np.array([np.nanmean(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])
# stderr_diff = np.array([np.nanstd(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])/np.sqrt(11)
means_diff = np.array([np.nanmean(np.array(set1[layer][regions_selected[idx]]) - np.array(set2[layer][regions_selected[idx]])) for layer in layers])
stderr_diff = np.array([np.nanstd(np.array(set1[layer][regions_selected[idx]]) - np.array(set2[layer][regions_selected[idx]])) for layer in layers])/np.sqrt(11)
sns.lineplot(means_diff, ax=axis, c=palette[-2], linewidth=3)
axis.errorbar(x=range(len(means_diff)), y=means_diff, yerr=stderr_diff, fmt='o', capsize=2, c=palette[-2])
friedman
diffMeas = np.array([set2[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']]) - np.array([set1[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']])
friedman = friedmanchisquare(*diffMeas)
print(friedman)
# apply multcompare
eng = matlab.engine.start_matlab()
diffMeas_ml = matlab.double(diffMeas.tolist())
eng.workspace['diffMeas'] = diffMeas_ml
eng.eval("[p,tbl,stats] = friedman(diffMeas',1);", nargout=0)
c = eng.eval("multcompare(stats);", nargout=1)
c_py = np.array(c)
eng.quit()
print(pd.DataFrame(c_py[:,2:], index=['D-M', 'D-S', 'M-S'], columns=['Lower CI', 'Mean rank Diff', 'Upper CI', 'p-value']).to_latex())
if regions_selected[idx] == 'HG':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(2.3-2.3+0.2, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
# axis.axhline(2.3-2.3+0.2, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
if regions_selected[idx] == 'pSTG':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(2.3-2.3+0.2, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
axis.axhline(0.85, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
if regions_selected[idx] == 'pMTG':
axis.axhline(0.9, 0.05, 0.95, c='gray', linewidth=1.5) # D-S
# axis.axhline(2.3-2.3+0.2, 0.05, 0.45, c='gray', linewidth=1.5) # D-M
# axis.axhline(2.3-2.3+0.2, 0.55, 0.95, c='gray', linewidth=1.5) # M-S
axis.set_xticks([0,1,2], [r'$\textbf{D}$', r'$\textbf{M}$', r'$\textbf{S}$'])
axis.set_ylim([0.0, 1.0])
if idx in [0,4]:
axis.set_ylabel(r'$\mathbf{\Delta\beta\ (\%)}$', fontsize=20)
axis.set_title(r'$\textbf{'+regions_selected[idx]+r'}$')
# blue_line = mlines.Line2D([], [], color='b', linestyle='-', linewidth=2, label='Valid')
# red_line = mlines.Line2D([], [], color='r', linestyle='-', linewidth=2, label='Invalid')
# fig.legend(handles=[red_line, blue_line], loc='center', bbox_to_anchor=(.5, 0.02))
# Mean Activations of Congruent vs. Invalid Trials (11 Subjects, Averaged Hemisphere, Non-Parametric FDR Corrected
plt.savefig('./figures/univariate_laminar_congomis_randomis_nonparametric_fdr_friedman_posthoc.png', dpi=300, bbox_inches='tight')
plt.show()
# %%
diffs = []
for idx, axis in enumerate(fig.get_axes()): # regions_selected
means_diff = np.array([np.mean(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])
stderr_diff = np.array([np.std(np.array(set2[layer][regions_selected[idx]]) - np.array(set1[layer][regions_selected[idx]])) for layer in layers])/np.sqrt(11)
sns.lineplot(means_diff, ax=axis, c=palette[-2], linewidth=3)
axis.errorbar(x=range(len(means_diff)), y=means_diff, yerr=stderr_diff, fmt='o', capsize=2, c=palette[-2])
# friedman
diffMeas = np.array([set2[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']]) - np.array([set1[item][regions_selected[idx]] for item in ['L1', 'L2', 'L3']]).astype(np.float64)
print(np.mean(diffMeas, axis=1), '\t', regions_selected[idx])
diffs.append(np.mean(diffMeas, axis=1))
diffs = np.array(diffs) * -1
np.min(diffs), np.max(diffs)
# %%
mincol, maxcol
# %%
all_regs = [entry['Name'] for entry in smp_head['Map']]
indcs = np.array([0,4,5])
selected = ['HG', 'pSTG', 'pMTG']
# Superficial
# _, indcs, _ = np.intersect1d(all_regs, selected, return_indices=True)
smp_head_cp = deepcopy(smp_head)
smp_head_cp['Nr maps'] = len(indcs)
smp_head_cp['Map'] = list(np.array(smp_head_cp['Map'])[indcs])
data = deepcopy(smp_data[:, indcs])
for idx, tval in enumerate(diffs[:,0]):
data[:, idx][data[:, idx]==6] = tval
smp_head_cp['Map'][idx]['Threshold min'] = mincol
smp_head_cp['Map'][idx]['Threshold max'] = maxcol
print(smp_head_cp['Map'][idx]['Threshold min'], smp_head_cp['Map'][idx]['Threshold max'])
bvbabel.smp.write_smp('/mnt/hdd2/associative_learning/paper_draft_v1/figures/group_patches_LH_univariate_predomit_unpredomit_dep.smp', smp_head_cp, data)
# Middle
# _, indcs, _ = np.intersect1d(all_regs, selected, return_indices=True)
smp_head_cp = deepcopy(smp_head)
smp_head_cp['Nr maps'] = len(indcs)
smp_head_cp['Map'] = list(np.array(smp_head_cp['Map'])[indcs])
data = deepcopy(smp_data[:, indcs])
for idx, tval in enumerate(diffs[:,1]):
data[:, idx][data[:, idx]==6] = tval
smp_head_cp['Map'][idx]['Threshold min'] = mincol
smp_head_cp['Map'][idx]['Threshold max'] = maxcol
print(smp_head_cp['Map'][idx]['Threshold min'], smp_head_cp['Map'][idx]['Threshold max'])
bvbabel.smp.write_smp('/mnt/hdd2/associative_learning/paper_draft_v1/figures/group_patches_LH_univariate_predomit_unpredomit_mid.smp', smp_head_cp, data)
# Deep
# _, indcs, _ = np.intersect1d(all_regs, selected, return_indices=True)
smp_head_cp = deepcopy(smp_head)
smp_head_cp['Nr maps'] = len(indcs)
smp_head_cp['Map'] = list(np.array(smp_head_cp['Map'])[indcs])
data = deepcopy(smp_data[:, indcs])
for idx, tval in enumerate(diffs[:,2]):
data[:, idx][data[:, idx]==6] = tval
smp_head_cp['Map'][idx]['Threshold min'] = mincol
smp_head_cp['Map'][idx]['Threshold max'] = maxcol
print(smp_head_cp['Map'][idx]['Threshold min'], smp_head_cp['Map'][idx]['Threshold max'])
bvbabel.smp.write_smp('/mnt/hdd2/associative_learning/paper_draft_v1/figures/group_patches_LH_univariate_predomit_unpredomit_sup.smp', smp_head_cp, data)
# %%
# %%
# %%
# %% [markdown]
# ---
#
# # <center> $\Large\textbf{Unpredictable vs. Omitted (Unpredictable)}$ </center>
#
# **Test the effect of prediction.**
# %%
regex1 = 'MANUAL' #'^random.*_(?!omis)$'
regex2 = '^random.*_omis$'
set1 = {};set2 = {}
for layer_name, layer_id in zip(['Deep', 'Middle', 'Superficial'], ['L1', 'L2', 'L3']):
tstats[layer_id] = []; pvals[layer_id] = []; df_set1 = {}; df_set2 = {}
for region in regions:
df_set1[region] = []
df_set2[region] = []
for subject in subjects:
stimuli = tbx.utils.extract_stimuli(subject)
betas_l = tbx.utils.extract_condition_betas(stimuli=stimuli,
betas=dataset[subject]['LH'][region]['betas_selected'][layer_id].values.T,
stimulus_type='sound')
betas_r = tbx.utils.extract_condition_betas(stimuli=stimuli,
betas=dataset[subject]['RH'][region]['betas_selected'][layer_id].values.T,
stimulus_type='sound')
df_set1[region].append(np.mean([
np.nanmean(betas_l.filter(regex='^random').filter(regex='^(?!.*_omis$)').values),
np.nanmean(betas_r.filter(regex='^random').filter(regex='^(?!.*_omis$)').values)]))
df_set2[region].append(np.mean([
np.nanmean(betas_l.filter(regex=regex2).values),
np.nanmean(betas_r.filter(regex=regex2).values)]))
set1[layer_id] = df_set1
set2[layer_id] = df_set2
all_non_parametric_pvals = np.array([[non_parametric_permutation_ttest(np.subtract(set2[layer][region], set1[layer][region])) for layer in layers] for region in regions])
all_non_parametric_pvals_flat = np.array(all_non_parametric_pvals).ravel()
reject_fdr, p_adjusted = fdrcorrection(all_non_parametric_pvals_flat)
p_adjusted = p_adjusted.reshape(np.array(all_non_parametric_pvals).shape)
# %%
all_layer_pvals = []
for region in regions:
for layerA, layerB in itertools.combinations(layers, 2):
diff_betas_layerA = np.subtract(set2[layerA][region], set1[layerA][region]) # difference within layerA
diff_betas_layerB = np.subtract(set2[layerB][region], set1[layerB][region]) # difference within layerB
all_layer_pvals.append(non_parametric_permutation_ttest(diff_betas_layerA - diff_betas_layerB))
reject_fdr, p_layers_adjusted = fdrcorrection(all_layer_pvals)
all_non_parametric_pvals_for_layers = {}
idx = 0
for region in regions:
all_non_parametric_pvals_for_layers[region] = {}
for layerA, layerB in itertools.combinations(layers, 2):
all_non_parametric_pvals_for_layers[region]['_'.join([layerA, layerB])] = p_layers_adjusted[idx]
idx += 1
fig = plt.figure(figsize=(9, 9))
gs = gridspec.GridSpec(3, 3, height_ratios=[1, 1, 1], hspace=0.6, wspace=0.4)
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[0, 2])
ax4 = fig.add_subplot(gs[1, 0])
ax5 = fig.add_subplot(gs[1, 1])
ax6 = fig.add_subplot(gs[1, 2])
ax7 = fig.add_subplot(gs[2, 0])
ax8 = fig.add_subplot(gs[2, 1])
ax9 = fig.add_subplot(gs[2, 2])
for idx, axis in enumerate(fig.get_axes()): # regions
means1 = np.array([np.nanmean(set1[layer][regions[idx]]) for layer in layers])
stderr1 = np.array([np.nanstd(set1[layer][regions[idx]]) for layer in layers])/np.sqrt(11)
sns.lineplot(means1, ax=axis, c='b', linewidth=3)
axis.errorbar(x=range(len(means1)), y=means1, yerr=stderr1, fmt='o', capsize=2, c='b')
means2 = np.array([np.nanmean(set2[layer][regions[idx]]) for layer in layers])
stderr2 = np.array([np.nanstd(set2[layer][regions[idx]]) for layer in layers])/np.sqrt(11)
sns.lineplot(means2, ax=axis, c='r', linewidth=3)
axis.errorbar(x=range(len(means2)), y=means2, yerr=stderr2, fmt='o', capsize=2, c='r')
axis.set_xticks([0,1,2], [r'$\textbf{D}$', r'$\textbf{M}$', r'$\textbf{S}$'])
axis.set_ylim([0., 2.5])
if idx in [0,3,6]:
axis.set_ylabel(r'$\mathbf{\beta\ (\%)}$', fontsize=20)
axis.set_title(r'$\textbf{'+regions[idx]+r'}$')
# make border thicker when there is an interaction
if regions[idx] in ['PT', 'HG', 'aSTG', 'pSTG', 'pMTG', 'parsOrbitalis', 'TPOj']:
for spine in axis.spines.values():
spine.set_linewidth(2)
spine.set_color('black')
# add layer interactions
if regions[idx] == 'HG':
axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'PT':
axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
# axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
# if regions[idx] == 'aSTG':
# axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
# axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
# axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
# axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
# axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
# axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'pSTG':
axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
# axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
# axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'pMTG':
axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
# axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
# axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
# axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'parsOrbitalis':
# axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
if regions[idx] == 'TPOj':
# axis.axhline(2.4, 0.05, 0.95, c='b', linewidth=1.5) # D-S
axis.axhline(2.3, 0.05, 0.45, c='b', linewidth=1.5) # D-M
# axis.axhline(2.3, 0.55, 0.95, c='b', linewidth=1.5) # M-S
# axis.axhline(2.2, 0.05, 0.95, c='r', linewidth=1.5) # D-S
# axis.axhline(2.1, 0.05, 0.45, c='r', linewidth=1.5) # D-M
# axis.axhline(2.1, 0.55, 0.95, c='r', linewidth=1.5) # M-S
blue_line = mlines.Line2D([], [], color='b', linestyle='-', linewidth=2, label='Unpredictable')
red_line = mlines.Line2D([], [], color='r', linestyle='-', linewidth=2, label='Omitted (Unpredictable)')
fig.legend(handles=[red_line, blue_line], loc='center', bbox_to_anchor=(.5, 0.))
# Mean Activations of Congruent vs. Invalid Trials (11 Subjects, Averaged Hemisphere, Non-Parametric FDR Corrected
# plt.savefig('./results/univariate_laminar_congomit_randomit_nonparametric_fdr.png', dpi=300, bbox_inches='tight')
plt.savefig('./figures/univariate_laminar_rand_randomis_nonparametric_fdr_friedman.png', dpi=300, bbox_inches='tight')
plt.show()
# %%
# friedman
pvals = []
stats = []
for roi in regions:
cong = np.array([set1[item][roi] for item in ['L1', 'L2', 'L3']])
incong = np.array([set2[item][roi] for item in ['L1', 'L2', 'L3']])
statistic, pvalue = friedmanchisquare(*(incong-cong))
pvals.append(pvalue)
stats.append(statistic)
# print(f'{roi:<16}: statistic={statistic:0.8f}, \t pvalue={pvalue:0.8f}', end='\t')
# if pvalue < 0.05:
# print('*')
# else:
# print()
data = [{'ROI':roi, 'stat':f'{stats[idx]:0.8f}', 'p':f'{pvals[idx]:0.8f}'} for idx, roi in enumerate(regions)]
headers = data[0].keys()
latex_table_uncorrected = (
# "\\begin{table}[h!]\n"
# "\\centering\n"
"\\caption{Friedman}\n"
"\\begin{tabular}{" + " ".join(["l"] * len(headers)) + "}\n"
"\\toprule\n"
)
latex_table_uncorrected += " & ".join(headers) + " \\\\\n\\midrule\n"
for row in data:
# if float(row['p']) < 0.05:
# latex_table_uncorrected += " & ".join('\\textbf{'+str(row[h])+'}' for h in headers) + " \\\\\n"
# else:
latex_table_uncorrected += " & ".join(str(row[h]) for h in headers) + " \\\\\n"
latex_table_uncorrected += (
"\\bottomrule\n"
"\\end{tabular}\n"
# "\\end{table}"
)
#%%%%%%%%%%%%%
rejections, pvals_corrected = fdrcorrection(pvals)
for idx, roi in enumerate(regions):
statistic = stats[idx]
pvalue = pvals_corrected[idx]
# print(f'{roi:<16}: statistic={statistic:0.8f}, \t pvalue={pvalue:0.8f}', end='\t')
# if pvalue < 0.05:
# print('*')
# else:
# print()
data = [{'ROI':roi, 'stat':f'{stats[idx]:0.8f}', 'p':f'{pvals_corrected[idx]:0.8f}'} for idx, roi in enumerate(regions)]
headers = data[0].keys()
latex_table = (
# "\\begin{table}[h!]\n"
# "\\centering\n"
"\\caption{Friedman (FDR)}\n"
"\\begin{tabular}{" + " ".join(["l"] * len(headers)) + "}\n"
"\\toprule\n"
)
latex_table += " & ".join(headers) + " \\\\\n\\midrule\n"
for row in data:
# if float(row['p']) < 0.05:
# latex_table += " & ".join('\\textbf{'+str(row[h])+'}' for h in headers) + " \\\\\n"
# else:
latex_table += " & ".join(str(row[h]) for h in headers) + " \\\\\n"
latex_table += (
"\\bottomrule\n"
"\\end{tabular}\n"
# "\\end{table}"
)
print('''\\begin{table}[htbp]
\\centering''')
print('\\begin{minipage}{0.45\\textwidth}')
print(latex_table_uncorrected)
print('\\end{minipage}')
print('\\hfill\\begin{minipage}{0.45\\textwidth}')
print(latex_table)
print('\\end{minipage}')
print('\\end{table}')
print()
# %%
# interactions
inter = ['PT', 'HG', 'aSTG', 'pSTG', 'pMTG', 'parsOrbitalis', 'TPOj']
# anova
all_results = []
for roi in inter:
cong = np.array([set1[item][roi] for item in ['L1', 'L2', 'L3']])
incong = np.array([set2[item][roi] for item in ['L1', 'L2', 'L3']])
result = {
'D-M (cong)': non_parametric_permutation_ttest(cong[0,:]-cong[1,:]),
'D-S (cong)': non_parametric_permutation_ttest(cong[0,:]-cong[2,:]),
'M-S (cong)': non_parametric_permutation_ttest(cong[1,:]-cong[2,:]),
'D-M (incong)': non_parametric_permutation_ttest(incong[0,:]-incong[1,:]),
'D-S (incong)': non_parametric_permutation_ttest(incong[0,:]-incong[2,:]),
'M-S (incong)': non_parametric_permutation_ttest(incong[1,:]-incong[2,:])
}
all_results.append(list(result.values()))
# correct and print
rej, pcorr = fdrcorrection(np.array(all_results).ravel())
pcorr = pcorr.reshape(len(inter), 6)
for idx, roi in enumerate(inter):
cong = np.array([set1[item][roi] for item in ['L1', 'L2', 'L3']])
incong = np.array([set2[item][roi] for item in ['L1', 'L2', 'L3']])
result = {
'D-M (cong)': non_parametric_permutation_ttest(cong[0,:]-cong[1,:]),
'D-S (cong)': non_parametric_permutation_ttest(cong[0,:]-cong[2,:]),
'M-S (cong)': non_parametric_permutation_ttest(cong[1,:]-cong[2,:]),
'D-M (incong)': non_parametric_permutation_ttest(incong[0,:]-incong[1,:]),
'D-S (incong)': non_parametric_permutation_ttest(incong[0,:]-incong[2,:]),
'M-S (incong)': non_parametric_permutation_ttest(incong[1,:]-incong[2,:])
}
stats = {
'D-M (cong)': np.round(ttest_rel(cong[0,:],cong[1,:]).statistic, 6),
'D-S (cong)': np.round(ttest_rel(cong[0,:],cong[2,:]).statistic, 6),
'M-S (cong)': np.round(ttest_rel(cong[1,:],cong[2,:]).statistic, 6),
'D-M (incong)': np.round(ttest_rel(incong[0,:],incong[1,:]).statistic, 6),
'D-S (incong)': np.round(ttest_rel(incong[0,:],incong[2,:]).statistic, 6),
'M-S (incong)': np.round(ttest_rel(incong[1,:],incong[2,:]).statistic, 6)
}
print('''\\begin{table}[htbp]
\\centering
\\caption{Layer Effects for separate conditions (\\textbf{''' + roi + '''})}
\\label{tab:my_label}''')
print('\\begin{minipage}{0.45\\textwidth}')
print('\\caption{Uncorrected}')
print(pd.DataFrame([list(stats.values()), list(result.values())], columns=result.keys(),index=['t', 'p']).T.to_latex(index=True))
print('\\end{minipage}')
print('\\begin{minipage}{0.45\\textwidth}')
print('\\caption{FDR}')
print(pd.DataFrame([list(stats.values()), pcorr[idx]], columns=result.keys(), index=['t', 'q']).T.to_latex(index=True))
print('\\end{minipage}')
print('\\end{table}')
print()
# %% [markdown]
# ---

Laminar_UnivariateAndStats.ipynb at commit 9b91341, no license · at the source

Overview

Authors: Jorie J G van Haren^1,2, Floris P de Lange^2,3, Sonja A Kotz^1,4, Federico De Martino^1,5

ORCID iDs: Jorie J G van Haren, Floris P de Lange, Sonja A Kotz, Federico De Martino

Institutions: Radboud University Nijmegen (Netherlands); Maastricht University (Netherlands); Donders Institute for Brain, Cognition and Behaviour (Netherlands); University of Bonn (Germany); Max Planck Institute for Human Cognitive and Brain Sciences (Germany); University of Minnesota (United States)

Journal: Nature communications, volume 17, issue 1, article 8748

Dates: received 26 November 2025; accepted 6 July 2026; published online 17 July 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1038/s41467-026-75662-w · PMID 42469220 · PMCID PMC13493366 · OpenAlex W7169569694

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: human (organism), cognitive (subfield)

Methods: Connectivity, Statistics, Machine learning, Preprocessing, fMRI & imaging, Single-unit activity, calcium imaging, Smoothing, state filtering, decompositions

Keywords: Perception, Neural encoding, Cortex

MeSH: Adaptation, Physiological*, Auditory Cortex*, Auditory Perception*, Acoustic Stimulation, Adult, Brain Mapping, Female, Humans, Magnetic Resonance Imaging, Male, Young Adult (* major topic)

Topic: Multisensory perception and integration (Experimental and Cognitive Psychology, Psychology), according to OpenAlex

Funding: EC | Horizon 2020 Framework Programme (101001270)

Citations: cited by 1 paper (Europe PMC); 67 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

Its files are read in the Code ↔ Paper reader above, with 13 matches between paragraphs and lines of code.

mesoscopic-computational-audition-lab/researchprojects

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 9b91341b5f6d97f089846c761d11751b7accf187, 3 July 2026

Languages: MATLAB (251), Python (82), Jupyter (48)

Size: 796 files, 381 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, 48 notebooks

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: NumPy (60 files), Psychtoolbox (46 files), Matplotlib (37 files), SciPy (37 files), pandas (36 files), seaborn (28 files), scikit-learn (14 files), statsmodels (11 files), NiBabel (8 files), pydicom (6 files), Statistics and Machine Learning Toolbox (5 files), Nilearn (3 files), scikit-posthocs (2 files), FSL (1 file), ggplot2 (1 file), imageio (1 file), MNE-Python (1 file), Nighres (1 file), Pillow (1 file), Pingouin (1 file), tidyverse (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

165 files

finalized/AssociativeLearning/analysis/Behavioral_ReactionTimes.ipynb, Jupyter, 448 lines
finalized/AssociativeLearning/analysis/Hippocampus_Decoding.ipynb, Jupyter, 377 lines
finalized/AssociativeLearning/analysis/Hippocampus_Decoding_revision1.ipynb, Jupyter, 705 lines
finalized/AssociativeLearning/analysis/Laminar_UnivariateAndStats.ipynb, Jupyter, 1,858 lines, 3 matches
finalized/AssociativeLearning/analysis/Laminar_UnivariateAndStats_revision1.ipynb, Jupyter, 2,055 lines
finalized/AssociativeLearning/preprocessing/bvpreproc.py, Python, 1,582 lines
finalized/PredAdapt_fMRI/misc/Tools/BIDS_Structure_Generator.py, Python, 177 lines
finalized/PredAdapt_fMRI/misc/Tools/Isovoxel_Nearest.py, Python, 494 lines
finalized/PredAdapt_fMRI/misc/Tools/LabelMap_to_WMGM_LHRH.py, Python, 107 lines
finalized/PredAdapt_fMRI/misc/Tools/Map_POI.py, Python, 172 lines
finalized/PredAdapt_fMRI/misc/Tools/NifTi_Tools.py, Python, 594 lines
finalized/PredAdapt_fMRI/misc/Tools/VMP_Cortical_depth.py, Python, 330 lines
finalized/PredAdapt_fMRI/misc/Tools/VOI_Tools.py, Python, 116 lines
finalized/PredAdapt_fMRI/misc/Tools/VTC_Box.py, Python, 160 lines
finalized/PredAdapt_fMRI/misc/abstract visualisation.ipynb, Jupyter, 629 lines
finalized/PredAdapt_fMRI/misc/mainexp visualisation.ipynb, Jupyter, 597 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/__init__.py, Python, 12 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/anatomical.py, Python, 367 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/anatomical_nighres.py, Python, 126 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/anatomical_plot.py, Python, 115 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/bvbabel_light.py, Python, 332 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/coregister.py, Python, 85 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/coregister_plot.py, Python, 12 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/functional.py, Python, 424 lines, 1 match
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/functional_nordic.py, Python, 217 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/functional_plot.py, Python, 679 lines, 1 match
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/functional_topup.py, Python, 495 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/prepdicoms.py, Python, 307 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/utils.py, Python, 183 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/voxeltimecourse.py, Python, 233 lines
finalized/PredAdapt_fMRI/preproc_fmri/bv_preproc/voxeltimecourse_plot.py, Python, 93 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/STARTEXP.m, MATLAB, 66 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/Bitsi.m, MATLAB, 263 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/ang2pix.m, MATLAB, 38 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/calc_logistic_growth.m, MATLAB, 9 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/clean_timings.m, MATLAB, 42 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/compareloudness_MEG.m, MATLAB, 87 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/compareloudness_MRI.m, MATLAB, 91 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/compareloudness_leftright.m, MATLAB, 95 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/compareloudness_leftright_MEG.m, MATLAB, 90 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/counterbalance.m, MATLAB, 128 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/create_tones_main.m, MATLAB, 47 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/create_tones_tonotopy.m, MATLAB, 68 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/createwaveform.m, MATLAB, 10 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/disptext.m, MATLAB, 20 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/generate_frequencies_main.m, MATLAB, 48 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/get_bull_tex.m, MATLAB, 129 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/get_bull_tex_2.m, MATLAB, 141 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/logistic_func.m, MATLAB, 14 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/mainpredstims.m, MATLAB, 63 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/multilinetext.m, MATLAB, 16 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/quasirandom_sequence.m, MATLAB, 96 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/save_clean_timings.m, MATLAB, 10 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/savetempfile.m, MATLAB, 12 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/soundstim.m, MATLAB, 6 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/startup1.m, MATLAB, 10 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/waitforbitsi.m, MATLAB, 33 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/waitforbitsi_backup.m, MATLAB, 20 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/waitfornokey.m, MATLAB, 7 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/functions/waitforpulse.m, MATLAB, 21 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/leftrightequalization.m, MATLAB, 103 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/mainequalization.m, MATLAB, 120 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/mainpredsound.m, MATLAB, 305 lines, 1 match
finalized/PredAdapt_fMRI/psychtoolbox_exp/maintonotopy.m, MATLAB, 139 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/settings_main.m, MATLAB, 275 lines
finalized/PredAdapt_fMRI/psychtoolbox_exp/settings_tonotopy.m, MATLAB, 232 lines, 1 match
finalized/PredAdapt_fMRI/revision/ARTAnova.ipynb, Jupyter, 340 lines
finalized/PredAdapt_fMRI/revision/regression_blocked_effects.ipynb, Jupyter, 253 lines
finalized/PredAdapt_fMRI/revision/revision_analysis.ipynb, Jupyter, 588 lines
finalized/PredAdapt_fMRI/ses1_prf_estimation/explore_prf_with_permutations_noCV_bulk_f.m, MATLAB, 199 lines
finalized/PredAdapt_fMRI/ses1_prf_estimation/fitGaussianPrfbasedonRanks_MM.m, MATLAB, 93 lines
finalized/PredAdapt_fMRI/ses1_prf_estimation/get_gausssian_weigthsZ_MM_f.m, MATLAB, 30 lines
finalized/PredAdapt_fMRI/ses1_prf_estimation/get_penaltybasedonRanks_MM.m, MATLAB, 24 lines
finalized/PredAdapt_fMRI/ses1_prf_estimation/saveICAMap.m, MATLAB, 40 lines
finalized/PredAdapt_fMRI/ses1_prf_estimation/tonotopy - create prt.ipynb, Jupyter, 452 lines
finalized/PredAdapt_fMRI/ses1_single_trial_betas/BetaEstimation.m, MATLAB, 86 lines
finalized/PredAdapt_fMRI/ses1_single_trial_betas/BulkBetaEstimation.m, MATLAB, 101 lines
finalized/PredAdapt_fMRI/ses1_single_trial_betas/saveICAMap.m, MATLAB, 31 lines
finalized/PredAdapt_fMRI/ses2_modelstims/Adaptation/Long Trace PRF Demo.ipynb, Jupyter, 502 lines
finalized/PredAdapt_fMRI/ses2_modelstims/Adaptation/backup/Long Trace PRF - Forward Model-LAPTOP-U412C4ST.ipynb, Jupyter, 958 lines
finalized/PredAdapt_fMRI/ses2_modelstims/Adaptation/backup/Long Trace PRF - Forward Model.ipynb, Jupyter, 1,132 lines
finalized/PredAdapt_fMRI/ses2_modelstims/Adaptation/backup/Untitled.ipynb, Jupyter, 435 lines
finalized/PredAdapt_fMRI/ses2_modelstims/Adaptation/backup/hiergausfilt.ipynb, Jupyter, 24 lines
finalized/PredAdapt_fMRI/ses2_modelstims/Adaptation/backup/initiate_hgf.m, MATLAB, 36 lines
finalized/PredAdapt_fMRI/ses2_modelstims/Adaptation/backup/longtrace_adaptation_old.py, Python, 146 lines
finalized/PredAdapt_fMRI/ses2_modelstims/Adaptation/backup/tapas_fitModel2.m, MATLAB, 595 lines
finalized/PredAdapt_fMRI/ses2_modelstims/Adaptation/longtrace_adaptation.py, Python, 433 lines
finalized/PredAdapt_fMRI/ses2_modelstims/Adaptation/longtrace_adaptation_pres.py, Python, 129 lines
finalized/PredAdapt_fMRI/ses2_modelstims/Adaptation/longtrace_adaptation_timedomain.py, Python, 151 lines
finalized/PredAdapt_fMRI/ses2_modelstims/DREX/Untitled.ipynb, Jupyter, 13 lines
finalized/PredAdapt_fMRI/ses2_modelstims/DREX/demo.m, MATLAB, 228 lines
finalized/PredAdapt_fMRI/ses2_modelstims/DREX/display_DREX_output.m, MATLAB, 168 lines
finalized/PredAdapt_fMRI/ses2_modelstims/DREX/estimate_suffstat.m, MATLAB, 198 lines
finalized/PredAdapt_fMRI/ses2_modelstims/DREX/examples/extract_envelope.m, MATLAB, 39 lines
finalized/PredAdapt_fMRI/ses2_modelstims/DREX/examples/extract_spectralmoments.m, MATLAB, 71 lines
finalized/PredAdapt_fMRI/ses2_modelstims/DREX/examples/run_examples.m, MATLAB, 361 lines
finalized/PredAdapt_fMRI/ses2_modelstims/DREX/post_DREX_beliefdynamics.m, MATLAB, 59 lines
finalized/PredAdapt_fMRI/ses2_modelstims/DREX/post_DREX_changedecision.m, MATLAB, 58 lines
finalized/PredAdapt_fMRI/ses2_modelstims/DREX/post_DREX_prediction.m, MATLAB, 65 lines
finalized/PredAdapt_fMRI/ses2_modelstims/DREX/run_DREX_model.m, MATLAB, 992 lines, 2 matches
finalized/PredAdapt_fMRI/ses2_modelstims/DREX/rundrex_stims.m, MATLAB, 87 lines
finalized/PredAdapt_fMRI/ses2_modelstims/HGF/HGF Demo.ipynb, Jupyter, 320 lines
finalized/PredAdapt_fMRI/ses2_modelstims/HGF/__init__.py, Python, 11 lines
finalized/PredAdapt_fMRI/ses2_modelstims/HGF/hgf.py, Python, 614 lines
finalized/PredAdapt_fMRI/ses2_modelstims/HGF/hgf_config.py, Python, 392 lines
finalized/PredAdapt_fMRI/ses2_modelstims/HGF/hgf_fit.py, Python, 337 lines
finalized/PredAdapt_fMRI/ses2_modelstims/HGF/hgf_pres.py, Python, 345 lines
finalized/PredAdapt_fMRI/ses2_modelstims/HGF/hgf_sim.py, Python, 165 lines
finalized/PredAdapt_fMRI/ses2_modelstims/stim_io.py, Python, 1,326 lines
finalized/PredAdapt_fMRI/ses2_modelstims/stim_io_plotting.py, Python, 269 lines
finalized/PredAdapt_fMRI/ses2_modelstims/stim_moddeling.ipynb, Jupyter, 354 lines
finalized/PredAdapt_fMRI/ses2_modelstims/vtc.py, Python, 263 lines
finalized/PredAdapt_fMRI/ses2_modelstims/vtc_masked.py, Python, 263 lines
finalized/PredAdapt_fMRI/ses2_regressions/ROI_and_layer_analysis.ipynb, Jupyter, 1,987 lines, 2 matches
finalized/PredAdapt_fMRI/ses2_regressions/regression.py, Python, 705 lines
finalized/PredAdapt_fMRI/ses2_regressions/regression_IdealObserver.ipynb, Jupyter, 247 lines
finalized/PredAdapt_fMRI/ses2_regressions/regression_adaptation_grid.ipynb, Jupyter, 718 lines
finalized/PredAdapt_fMRI/ses2_regressions/regression_main.ipynb, Jupyter, 318 lines
finalized/PredAdapt_fMRI/ses2_regressions/regression_split_pred.ipynb, Jupyter, 373 lines
finalized/PredAdapt_fMRI/ses2_regressions/save_maps.ipynb, Jupyter, 752 lines
finalized/PredAdapt_fMRI/ses2_regressions/stats.py, Python, 483 lines
finalized/PredAdapt_fMRI/ses2_regressions/varpar.py, Python, 579 lines
projects/PredAdapt_MEG/psychtoolbox_exp/STARTEXP.m, MATLAB, 69 lines
projects/PredAdapt_MEG/psychtoolbox_exp/Z_BITSITEST.m, MATLAB, 40 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/Bitsi.m, MATLAB, 263 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/ang2pix.m, MATLAB, 38 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/calc_logistic_growth.m, MATLAB, 9 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/clean_timings.m, MATLAB, 42 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/compareloudness_MEG.m, MATLAB, 87 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/compareloudness_MRI.m, MATLAB, 91 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/compareloudness_leftright.m, MATLAB, 95 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/compareloudness_leftright_MEG.m, MATLAB, 90 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/counterbalance.m, MATLAB, 130 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/create_tones_main.m, MATLAB, 47 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/create_tones_tonotopy.m, MATLAB, 68 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/createwaveform.m, MATLAB, 10 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/disptext.m, MATLAB, 20 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/generate_frequencies_main.m, MATLAB, 48 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/get_bull_tex.m, MATLAB, 129 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/get_bull_tex_2.m, MATLAB, 141 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/logistic_func.m, MATLAB, 14 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/mainpredstims.m, MATLAB, 63 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/multilinetext.m, MATLAB, 16 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/quasirandom_sequence.m, MATLAB, 96 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/save_clean_timings.m, MATLAB, 10 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/savetempfile.m, MATLAB, 12 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/soundstim.m, MATLAB, 6 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/startup1.m, MATLAB, 10 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/waitforbitsi.m, MATLAB, 33 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/waitforbitsi_backup.m, MATLAB, 20 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/waitfornokey.m, MATLAB, 7 lines
projects/PredAdapt_MEG/psychtoolbox_exp/functions/waitforpulse.m, MATLAB, 21 lines
projects/PredAdapt_MEG/psychtoolbox_exp/leftrightequalization.m, MATLAB, 103 lines
projects/PredAdapt_MEG/psychtoolbox_exp/mainequalization.m, MATLAB, 120 lines
projects/PredAdapt_MEG/psychtoolbox_exp/mainlocalizer.m, MATLAB, 160 lines
projects/PredAdapt_MEG/psychtoolbox_exp/mainpredsound.m, MATLAB, 309 lines
projects/PredAdapt_MEG/psychtoolbox_exp/mainpredsound_MEG.m, MATLAB, 339 lines, 1 match
projects/PredAdapt_MEG/psychtoolbox_exp/maintonotopy.m, MATLAB, 139 lines
projects/PredAdapt_MEG/psychtoolbox_exp/settings_localizer.m, MATLAB, 249 lines, 1 match
projects/PredAdapt_MEG/psychtoolbox_exp/settings_main.m, MATLAB, 275 lines
projects/PredAdapt_MEG/psychtoolbox_exp/settings_tonotopy.m, MATLAB, 232 lines
projects/PredAdapt_MEG/scripts_proc/IdealObserver_Demo.ipynb, Jupyter, 171 lines
projects/PredAdapt_MEG/scripts_proc/MNE and timingdata Demo.ipynb, Jupyter, 369 lines
projects/PredAdapt_MEG/scripts_proc/Regression.ipynb, Jupyter, 2,453 lines
repository limit reached (2,000 files or 30 MB): the rest is at the source (217 files)
README.md, Text, 7 lines

At the source: github.com/mesoScopic-Computational-AuditioN-lab/ResearchProjects/tree/main/finalized/PredAdapt_fMRI

mesoScopic-Computational-AuditioN-lab

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: the link answers

Software Heritage: not checked

Found in: “Code availability”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)

27 September 2026: the link answers (HTTP 200)

At the source: github.com/mesoScopic-Computational-AuditioN-lab

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: mesoScopic-Computational-AuditioN-lab, mesoscopic-computational-audition-lab/researchprojects

Read it in the paper: doi.org/10.1038/s41467-026-75662-w.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
164 scripts, each with its path and the digest of its content;
13 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

doi:10.18112/openneuro.ds006928.v1.0.0, at OpenNeuro; found in the references
openneuro:ds006928, at OpenNeuro; found in “Data availability”

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to a dataset: OpenNeuro ds006928
it points to the authors' code: mesoScopic-Computational-AuditioN-lab, mesoscopic-computational-audition-lab/researchprojects

Read it in the paper: doi.org/10.1038/s41467-026-75662-w.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 3 keywords, 11 MeSH terms, 1 funder, 63 references.

Cite

This paper

van Haren, J. J. G., de Lange, F. P., Kotz, S. A., & De Martino, F. (2026). Distinct Roles of Deep and Superficial Cortical Layers in Tone Prediction, Comparison, and Adaptation in Human Auditory Cortices. Nature communications, 17(1), 8748. https://doi.org/10.1038/s41467-026-75662-w

BibTeX

@article{vanharen2026distinct, author = {van Haren, Jorie J G and de Lange, Floris P and Kotz, Sonja A and De Martino, Federico}, title = {{Distinct Roles of Deep and Superficial Cortical Layers in Tone Prediction, Comparison, and Adaptation in Human Auditory Cortices}}, journal = {Nature communications}, year = {2026}, month = jul, volume = {17}, number = {1}, pages = {8748}, publisher = {Nature Publishing Group}, issn = {2041-1723}, doi = {10.1038/s41467-026-75662-w}, url = {https://doi.org/10.1038/s41467-026-75662-w}, pmid = {42469220}, pmcid = {PMC13493366} }

RIS

TY - JOUR AU - van Haren, Jorie J G AU - de Lange, Floris P AU - Kotz, Sonja A AU - De Martino, Federico TI - Distinct Roles of Deep and Superficial Cortical Layers in Tone Prediction, Comparison, and Adaptation in Human Auditory Cortices T2 - Nature communications J2 - Nat Commun PY - 2026 DA - 2026/07/17 VL - 17 IS - 1 SP - 8748 SN - 2041-1723 PB - Nature Publishing Group DO - 10.1038/s41467-026-75662-w UR - https://doi.org/10.1038/s41467-026-75662-w LA - en ER -

CSL-JSON

{ "id": "10.1038/s41467-026-75662-w", "type": "article-journal", "title": "Distinct Roles of Deep and Superficial Cortical Layers in Tone Prediction, Comparison, and Adaptation in Human Auditory Cortices", "container-title": "Nature communications", "author": [ { "family": "van Haren", "given": "Jorie J G" }, { "family": "de Lange", "given": "Floris P" }, { "family": "Kotz", "given": "Sonja A" }, { "family": "De Martino", "given": "Federico" } ], "container-title-short": "Nat Commun", "volume": "17", "issue": "1", "page": "8748", "DOI": "10.1038/s41467-026-75662-w", "PMID": "42469220", "PMCID": "PMC13493366", "ISSN": "2041-1723", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41467-026-75662-w", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 17 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.