Same Sentences, Different Grammars, Different Brain Responses?: An MEG Study on Case and Agreement Encoding in Hindi and Nepali Split-Ergative Structures.

Code ↔ Paper

5 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 5 matches

[1] § METHODS › Data Preprocessing ↔ scripts.zip/scripts/1_SgTrialRegressions_withMean_looseOrientation.py, lines 239–298 · score 0.87 · loose orientation, covariance matrices, inverse solution, warped, source space, vertices
[2] § METHODS › Data Analysis ↔ scripts.zip/scripts/1_SgTrialRegressions_withMean_looseOrientation.py, lines 492–549 · score 0.57 · linear regression, verb stem cloze, epochs
[3] § METHODS › Data Analysis ↔ scripts.zip/scripts/1_SgTrialRegressions_withMean_looseOrientation.py, lines 39–61 · score 0.54 · language network, sensor space, permutation, threshold, temporal, brain
[4] § RESULTS › MEG Results › Subject epoch analyses ↔ scripts.zip/scripts/1_SgTrialRegressions_withMean_looseOrientation.py, lines 39–61 · score 0.51 · temporal lobe, Sensor space, post, thresholded, onset, brain
[5] § METHODS › Data Analysis ↔ scripts.zip/scripts/1_SgTrialRegressions_withMean_looseOrientation.py, lines 111–176 · score 0.50 · 0–1000 ms, onset, space, language, Nepali, Hindi

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

Python · 1,336 lines · 47 KB · no license · 5 matches

### RUN THIS IN EELBRAIN
import csv
import mne
import numpy as np
import pandas as pd
import os
import csv
from mne.stats import fdr_correction, linear_regression, spatio_temporal_cluster_1samp_test
import pickle
import scipy
from matplotlib import rc
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import seaborn as sn
from helpers import *
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica Neue']})
rc('font',**{'family':'serif','serif':['Helvetica Neue']})
rc('text', usetex=False)
rePlot = False
decim = 4 # 1000 Hz -> 250 Hz
makeSTCs = False # load in fwd/cov/bem for each participant to make inv?
SNR = 2.0 # for regressions; 3.0 for ANOVAs
fixed = False # False for orientation free (=unsigned), True for fixed orientation (=signed)
# Parameters for spatio-temporal cluster-based permutation test
n_permutations = 10000
pThresh = 0.05
plotThresh = 0.10
tail = 0
searchTMin = int(000) # 200ms post-onset; units are MILISECONDS; 0 = minimal time
searchTMax = int(1000) # 800ms post-onset; units are MILISECONDS; 1000 = maximal time
title = date + ''
formulas = [
'Case_SubjNP',
#'Cloze_V_Morph',
#'Case_ObjNP+VGen',
]
spaces = [
# 'sensor', # sensor-space analysis
'source_wholebrain', # whole-brain analysis
#'source_bilat', # bilateral language network analyses
#'source_temp' # bilateral temporal lobe analyses
]
#'Cloze_V' <- just cloze manipulation
#formulas = ['Cloze_V_Morph']
# Construct the left and right hemisphere 'language network' search spaces for
# exploratory analyses for case manipulations, and left/right temporal lobe
# for more targeted N400 analyses of the cloze manipulation
labels = mne.read_labels_from_annot('fsaverage', 'aparc')
regions = ['bankssts', 'inferiorparietal', 'insula', 'lateralorbitofrontal', #'medialorbitofrontal',
'middletemporal', 'parsopercularis', 'parsorbitalis', 'parstriangularis', 'superiortemporal', 'supramarginal', 'temporalpole', 'transversetemporal']# 'inferiortemporal', 'fusiform', 'rostralmiddlefrontal', 'caudalmiddlefrontal']
tempRegions = ['bankssts', 'superiortemporal', 'middletemporal', 'transversetemporal', 'inferiorparietal', 'supramarginal', 'temporalpole', 'parsopercularis', 'parsorbitalis', 'parstriangularis', 'insula']
leftLngNet = mne.Label(hemi='lh')
leftTemp = mne.Label(hemi='lh')
leftLngNet.subject = 'fsaverage'
leftTemp.subject = 'fsaverage'
rightLngNet = mne.Label(hemi='rh')
rightLngNet.subject = 'fsaverage'
rightTemp = mne.Label(hemi='rh')
rightTemp.subject = 'fsaverage'
for region in regions:
for label in labels:
if region+'-lh' == label.name:
leftLngNet += label
elif region+'-rh' == label.name:
rightLngNet += label
for region in tempRegions:
for label in labels:
if region+'-lh' == label.name:
leftTemp += label
elif region+'-rh' == label.name:
rightTemp += label
bilatLngNet = mne.BiHemiLabel(lh=leftLngNet, rh=rightLngNet)
bilatTemp = mne.BiHemiLabel(lh=leftTemp, rh=rightTemp)
# These are the variables that are categorical (e.g. 1/0s) or
# that are factorial as part of the initial experiment design
catVars = ['UpcomingVGenF', 'GenF', 'ClozeF']
factVars = ['ObjCF', 'SubjCF', 'InteractionF']
# Each analysis is done on a 1000ms epoch
#times = np.arange(0,1001,1) # 0-1000ms in 1ms increments; time scale for all plots
times = np.arange(0,1000,4) # 0-1000ms in 1ms increments at 250Hz (4ms increments)
print('Reading in fsaverage source space...')
srcFName = os.path.join(subjects_dir, 'fsaverage', 'bem', 'fsaverage-ico-4-src.fif')
fsave_src = mne.read_source_spaces(fname=srcFName)
for language in ['Hindi', 'Nepali']:
if language == 'Hindi':
subjects = h_subjects
regressorFname = 'Hindi_regressors_withCorpusFreq_z.csv'
subExpName = 'HObjAgr'
else:
subjects = n_subjects
regressorFname = 'Nepali_regressors_withCorpusFreq.csv'
subExpName = 'NObjAgr'
# Get the regressors;
# these have been z-scored.
# Cloze values are drawn from the norming study
with open(regressorFname, 'r') as f:
regressors = [i for i in csv.DictReader(f)]
for formula in formulas:
if formula == 'Case_SubjNP':
model = ['Intercept', 'SubjCF', 'TrialNum', 'W1Length']#, 'W1Freq']
tmin = -1.0 # Units are in SECONDS; 0 = object onset
tmax = 0.0 # Units are in SECONDS; 0 = object onset
elif formula == 'Case_ObjNP':
model = ['Intercept', 'SubjCF', 'ObjCF', 'InteractionF', 'TrialNum', 'W2Length']#, 'W2Freq']
tmin = 0.0 # Units are in SECONDS; 0 = object onset
tmax = 1.0 # Units are in SECONDS; 0 = object onset
elif formula == 'Case_ObjNP+VGen':
if language == 'Hindi':
model = ['Intercept', 'SubjCF', 'ObjCF', 'InteractionF', 'TrialNum', 'W2Length', 'UpcomingVGenF']#, 'W2Freq']
else:
model = ['Intercept', 'SubjCF', 'ObjCF', 'InteractionF', 'TrialNum', 'W2Length', 'GenF']#, 'W2Freq']
tmin = 0.0 # Units are in SECONDS; 0 = object onset
tmax = 1.0 # Units are in SECONDS; 0 = object onset
elif formula == 'VGen_ObjNP':
if language == 'Hindi':
model = ['Intercept', 'UpcomingVGenF', 'TrialNum', 'W2Length']#, 'W2Freq']
else:
model = ['Intercept', 'GenF', 'TrialNum', 'W2Length']#, 'W2Freq']
tmin = 0.0 # Units are in SECONDS; 0 = object onset
tmax = 1.0 # Units are in SECONDS; 0 = object onset
elif formula == 'Cloze_V':
model = ['Intercept', 'ClozeF', 'TrialNum', 'W3Length']#, 'W3Freq']
tmin = 1.0 # Units are in SECONDS; 0 = object onset
tmax = 2.0 # Units are in SECONDS; 0 = object onset
elif formula == 'Cloze_V_Morph':
if language == 'Hindi':
model = ['Intercept', 'ClozeF', 'TrialNum', 'W3Length', 'UpcomingVGenF', 'SubjCF']#, 'W3Freq']
elif language == 'Nepali':
model = ['Intercept', 'ClozeF', 'TrialNum', 'W3Length', 'GenF', 'SubjCF']#, 'W3Freq']
tmin = 1.0 # Units are in SECONDS; 0 = object onset
tmax = 2.0 # Units are in SECONDS; 0 = object onset
# Check to see if we've already computed this participant's
# betas or not; if we have, then we load those
betas = dict()
betas_src = dict()
raw_conds = dict()
raw_conds_src = dict()
expConds = [
# SubjC:
'ErgPerf', 'NomImpf',
# ObjC:
'Acc', 'Bare',
# SubjC x ObjC:
'ErgAcc', 'ErgBare', 'NomAcc', 'NomBare',
# Cloze:
'LowCloze', 'HighCloze',
# VGen:
'MVerb_Hindi', 'FVerb_Hindi',
# VGen
'MVerb_Nepali', 'FVerb_Nepali']
expCondFactors = {
'ErgPerf':['MHiErgBare', 'MHiErgAcc', 'MLoErgBare', 'MLoErgAcc', 'FHiErgBare', 'FHiErgAcc', 'FLoErgBare', 'FLoErgAcc'],
'NomImpf': ['MHiBareBare', 'MHiBareAcc', 'MLoBareBare', 'MLoBareAcc', 'FHiBareBare', 'FHiBareAcc', 'FLoBareBare', 'FLoBareAcc'],
'Acc' : ['MHiErgAcc', 'MHiBareAcc', 'MLoBareAcc', 'MLoErgAcc', 'FHiBareAcc', 'FHiErgAcc', 'FLoBareAcc', 'FLoErgAcc'],
'Bare': ['MHiErgBare', 'MHiBareBare', 'MLoBareBare', 'MLoErgBare', 'FHiBareBare', 'FHiErgBare', 'FLoBareBare', 'FLoErgBare'],
'NomAcc': ['MHiBareAcc', 'MLoBareAcc', 'FHiBareAcc', 'FLoBareAcc'],
'NomBare': ['MHiBareBare', 'MLoBareBare','FHiBareBare','FLoBareBare'],
'ErgAcc': ['MHiErgAcc', 'MLoErgAcc', 'FHiErgAcc', 'FLoErgAcc'],
'ErgBare': ['MHiErgBare', 'MLoErgBare', 'FHiErgBare', 'FLoErgBare'],
'LowCloze': ['MLoBareAcc', 'FLoBareAcc', 'MLoBareBare', 'FLoBareBare', 'MLoErgAcc', 'FLoErgAcc', 'MLoErgBare', 'FLoErgBare'],
'HighCloze': ['MHiBareAcc', 'FHiBareAcc', 'MHiBareBare', 'FHiBareBare', 'MHiErgAcc', 'FHiErgAcc', 'MHiErgBare', 'FHiErgBare'],
'MVerb_Nepali': ['MLoBareAcc', 'MLoBareBare', 'MLoErgAcc', 'MLoErgBare'],
'FVerb_Nepali': ['FLoBareAcc', 'FLoBareBare', 'FLoErgAcc', 'FLoErgBare'],
'MVerb_Hindi': ['MLoBareAcc', 'MLoBareBare', 'MLoErgAcc', 'FLoErgAcc', 'FLoErgBare'],
'FVerb_Hindi': ['FLoBareAcc', 'FLoBareBare', 'MLoErgBare']
}
for cond in expConds:
raw_conds[cond] = []
raw_conds_src[cond] = []
for cond in model:
betas[cond] = []
betas_src[cond] = []
for subject in subjects:
indPlotFolder = os.path.join('/Volumes', 'falconlab', 'experiments', expName, 'plots', 'individuals', language, subject)
if not os.path.exists(indPlotFolder):
os.makedirs(indPlotFolder)
MEGDir = os.path.join('/Volumes', 'falconlab', 'data', 'meg', subject, subExpName)
for fName in os.listdir(MEGDir):
if fName.endswith('retainedEpochs.csv'):
retEposFname = fName
STCDir = os.path.join('/Volumes', 'falconlab', 'data', 'stc', subject, expName)
if not os.path.exists(STCDir):
os.makedirs(STCDir)
subjBetas = dict()
subjBetas_src = dict()
print('Reading in raw data for subject %s' %subject)
epos = mne.read_epochs(os.path.join(MEGDir, '%s-ica-good-epo.fif' %(subject)))
epos = epos.resample(sfreq=250)
if makeSTCs:
print('Calculating covariance matrix for subject %s' %subject)
cov = mne.compute_covariance(epos,tmin=baseline[0],tmax=baseline[1], method=['shrunk', 'diagonal_fixed', 'empirical'])
print('Reading in forward solution, trans file, source space, and morphing to fsaverage...')
fwd_fname = os.path.join(MEGDir, subject+'-fwd.fif')
trans = os.path.join(MEGDir, subject+'-trans.fif')
# cov_fname = os.path.join(language, 'MEG', subject, subject+'-GrAvgFixed-cov.fif')
print('Reading in source space...')
src = mne.read_source_spaces(os.path.join(subjects_dir, subject, 'bem', subject+'-ico-4-src.fif'))
print('Warping to fsaverage vertices...')
morph = mne.compute_source_morph(src, subject_from=subject, subject_to='fsaverage', spacing=4, subjects_dir=subjects_dir)
try:
print('Reading in BEM solution...')
# bemFName = os.path.join(subjects_dir, subject, 'bem', subject+'-5120-5120-5120-bem-sol.fif')
bemFName = os.path.join(subjects_dir, subject, 'bem', subject+'-inner_skull-bem-sol.fif')
bem = mne.read_bem_solution(bemFName)
except:
print("Can't read in BEM solution for subject %s, calculating it again..." %subject)
conductivity = (0.3, #0.006, 0.3
) # for three layers layer
bemModel = mne.make_bem_model(subject=subject, ico=4, conductivity=conductivity, subjects_dir=subjects_dir)
bem = mne.make_bem_solution(bemModel)
mne.write_bem_solution(bemFName, bem, overwrite=True)
print('Reading in forward solution...')
# if os.path.exists(fwd_fname):
# fwd = mne.read_forward_solution(fwd_fname)
# else:
fwd = mne.make_forward_solution(info=epos.info,
trans=trans, src=src, bem=bem,
eeg=False,
meg=True,
# mag=True,
ignore_ref=True)
print('Calculating inverse solution...')
inv = mne.minimum_norm.make_inverse_operator(epos.info, fwd, cov, depth=None, fixed=fixed, loose=0.20)
print('Checking on STCs...')
if rePlot:
retEposFile = open(os.path.join(MEGDir, retEposFname),'r')
retEpos = pd.DataFrame([i for i in csv.DictReader(retEposFile)])
retEposFile.seek(0)
retEposFile.close()
totalObs = len(retEpos)
totalAcc = len(retEpos.query("`Hit?` == 'Hit'")) / len(retEpos.query("`Hit?` != 'NaN'"))
with open(os.path.join(indPlotFolder, 'behavioral.txt'), 'w') as f:
f.write('Total epochs retained: ' + str(totalObs) + '\n')
f.write('Total correct: ' + str(totalAcc))
print('Replotting by-subject data for subject %s....' %subject)
print('Plotting butterfly plot')
fig = epos.average().plot(gfp=True, show=False)
fig = tidyFig(fig)
fig.savefig(os.path.join(indPlotFolder, 'butterflyPlot_REST.png'))
evokedPlot = epos.average()
whitePlot = mne.viz.plot_evoked_white(evokedPlot, cov, show=False)
whitePlot.savefig(os.path.join(indPlotFolder, 'whitenedCov.png'))
try:
srcAlign = mne.viz.plot_alignment(epos.info, subject=subject, src=src, trans=trans,
#trans = 'fsaverage',
surfaces=dict(white=0.4, outer_skull=0.2, head=0.5),
# ['head-dense', 'outer_skull', 'white'],
fwd = fwd,
bem = bem,
# eeg = 'projected')
# eeg=dict(original=0.0, projected = 0.5))
)
screenshot = srcAlign.plotter.screenshot()
fig, ax = plt.subplots(figsize=(12,12))
ax.imshow(screenshot, origin='upper')
ax.set_axis_off() # Disable axis labels and ticks
fig.tight_layout()
fig.savefig(os.path.join(indPlotFolder, 'srcAlign.png'), dpi=300)
except:
srcAlign = mne.viz.plot_alignment(epos.info, subject=subject, src=src, trans=trans,
#trans = 'fsaverage',
surfaces=dict(white=0.4, head=0.5),
# ['head-dense', 'outer_skull', 'white'],
fwd = fwd,
bem = bem,
# eeg = 'projected')
# eeg=dict(original=0.0, projected = 0.5))
)
screenshot = srcAlign.plotter.screenshot()
fig, ax = plt.subplots(figsize=(12,12))
ax.imshow(screenshot, origin='upper')
ax.set_axis_off() # Disable axis labels and ticks
fig.tight_layout()
fig.savefig(os.path.join(indPlotFolder, 'srcAlign.png'), dpi=300)
cov = mne.compute_covariance(epos,tmin=baseline[0],tmax=baseline[1], method=['shrunk', 'diagonal_fixed', 'empirical'])
epos.crop(tmin=tmin, tmax=tmax)
epos.pick(picks=['meg'])
for cond in expConds:
if language == 'Hindi' and 'Nepali' in cond:
pass
elif language == 'Nepali' and 'Hindi' in cond:
pass
else:
stcFName = os.path.join(STCDir, language + '_' + cond + '_' + str(int(tmin)) + '-' + str(int(tmax)) + '_src_raw_loose_250Hz')
if os.path.exists(stcFName+'-lh.stc'):
print('STC for %s trials for subject %s in time window %s - %s already exists!' %(cond, subject, str(int(tmin)), str(int(tmax))))
stc = mne.read_source_estimate(stcFName)
elif makeSTCs:
print('Calculating STC for %s trials for subject %s in time window %s - %s' %(cond, subject, str(int(tmin)), str(int(tmax))))
stc = morph.apply(mne.minimum_norm.apply_inverse(epos[expCondFactors[cond]].average(), inv, lambda2 = 1.0 / 3.0 ** 2.0, method='dSPM'))
stc.save(stcFName)
raw_conds_src[cond].append(stc)
print('Averaging raw sensor + src space data by key comparison in experiment design')
for cond in expConds:
if language == 'Hindi' and 'Nepali' in cond:
pass
elif language == 'Nepali' and 'Hindi' in cond:
pass
else:
raw_conds[cond].append(epos[expCondFactors[cond]].average())
if os.path.exists(os.path.join(STCDir, language+'_'+formula+'_sns_250Hz_signed.pickle')) and os.path.exists(os.path.join(STCDir, language+'_'+formula+'_src_250Hz_loose.pickle')):
print('**********Loading betas for subject %s' %subject)
with open(os.path.join(STCDir, language+'_'+formula+'_sns_250Hz_signed.pickle'),'rb') as f:
subjBetas = pickle.load(f)
with open(os.path.join(STCDir, language+'_'+formula+'_src_250Hz_loose.pickle'),'rb') as f:
subjBetas_src = pickle.load(f)
else:
MEGDir = os.path.join('/Volumes', 'falconlab', 'data', 'meg', subject, subExpName)
retEposFile = open(os.path.join(MEGDir, retEposFname),'r')
retEpos = pd.DataFrame([i for i in csv.DictReader(retEposFile)])
retEposFile.seek(0)
retEposFile.close()
totalObs = len(retEpos)
totalAcc = len(retEpos.query("`Hit?` == 'Hit'")) / len(retEpos.query("`Hit?` != 'NaN'"))
with open(os.path.join(indPlotFolder, 'behavioral.txt'), 'w') as f:
f.write('Total epochs retained: ' + str(totalObs) + '\n')
f.write('Total correct: ' + str(totalAcc))
evokedPlot = epos.average()
whitePlot = mne.viz.plot_evoked_white(evokedPlot, cov, show=False)
whitePlot.savefig(os.path.join(indPlotFolder, 'whitenedCov.png'))
print('Plotting butterfly plot')
fig = epos.average().plot(gfp=True, show=False)
fig = tidyFig(fig)
fig.savefig(os.path.join(indPlotFolder, 'butterflyPlot_REST.png'))
# Initialize the parameters that we're including as regressors
epos.metadata['W1Freq'] = 0
epos.metadata['W2Freq'] = 0
epos.metadata['W3Freq'] = 0
epos.metadata['W1Length'] = 0
epos.metadata['W2Length'] = 0
epos.metadata['W3Length'] = 0
epos.metadata['VStemCloze'] = 0
epos.metadata['GenF'] = 0 # 0 is for MF
epos.metadata['SubjCF'] = 0 # Nom; 1 = Erg
epos.metadata['ObjCF'] = 0 # Acc; 1 = Bare
epos.metadata['ClozeF'] = 0 # Hi; 1 = Lo
epos.metadata['InteractionF'] = 0 # For all conditions except for Erg:Bare; 1 = Erg:Bare
epos.metadata['UpcomingVGenF'] = 0 # For masculine
# Assign an intercept
epos.metadata['Intercept'] = 1
for i, row in epos.metadata.iterrows():
# I didn't include the itemset number in
# the log files, so we need to reconstruct this
# by looking at the individual words
w1 = row['W1']
w2 = row['W2']
w3 = row['W3']
epos.metadata.at[i, 'W1Length'] = len(w1)
epos.metadata.at[i, 'W2Length'] = len(w2)
epos.metadata.at[i, 'W3Length'] = len(w3)
if row['Gen'] == 'FM' or row['Gen'] == 'F':
epos.metadata.at[i, 'GenF'] = 1
if row['SubjC'] == 'Erg':
epos.metadata.at[i, 'SubjCF'] = 1
if row['ObjC'] == 'Bare':
epos.metadata.at[i, 'ObjCF'] = 1
if row['SubjC'] == 'Erg' and row['ObjC'] == 'Bare':
epos.metadata.at[i, 'InteractionF'] = 1
if row['SubjC'] == 'Bare' and row['Gen'] == 'FM':
epos.metadata.at[i, 'UpcomingVGenF'] = 1
elif row['SubjC'] == 'Erg' and row['ObjC'] == 'Bare' and row['Gen'] == 'MF':
epos.metadata.at[i, 'UpcomingVGenF'] = 1
if row['Cloze'] == 'Lo':
epos.metadata.at[i, 'ClozeF'] = 1
# The regressors file stores the low vs. high cloze
# verb stems in different columns, so this helps
# pick the right value
clozeCond = row['Cloze']
foundRegressors = False
for x in regressors:
if not foundRegressors:
if w1 == x['W1'] and w2 == x['W2'] and w3 == x['W3'] and clozeCond == x['Cloze']:
epos.metadata.at[i, 'W1Freq'] = float(x['W1Freq'])
epos.metadata.at[i, 'W2Freq'] = float(x['W2Freq'])
epos.metadata.at[i, 'W3Freq'] = float(x['W3Freq'])
if clozeCond == 'Lo':
epos.metadata.at[i, 'VStemCloze'] = float(x['LoP'])
elif clozeCond == 'Hi':
epos.metadata.at[i, 'VStemCloze'] = float(x['HiP'])
else:
print("Couldn't find regressor")
epos.metadata['TrialNum'] = epos.metadata['TrialNum'].astype(float)
try:
srcAlign = mne.viz.plot_alignment(epos.info, subject=subject, src=src, trans=trans,
#trans = 'fsaverage',
surfaces=dict(white=0.4, outer_skull=0.2, head=0.5),
# ['head-dense', 'outer_skull', 'white'],
fwd = fwd,
bem = bem,
# eeg = 'projected')
# eeg=dict(original=0.0, projected = 0.5))
)
screenshot = srcAlign.plotter.screenshot()
fig, ax = plt.subplots(figsize=(12,12))
ax.imshow(screenshot, origin='upper')
ax.set_axis_off() # Disable axis labels and ticks
fig.tight_layout()
fig.savefig(os.path.join(indPlotFolder, 'srcAlign.png'), dpi=300)
except:
print("Can't plot source alignment")
res = linear_regression(epos.apply_function(abs), epos.metadata[model], names=model)
for cond in model:
if cond not in subjBetas:
subjBetas[cond] = list()
subjBetas_src[cond] = list()
subjBetas[cond].append(res[cond].beta)
print('Applying inverse solution...')
if makeSTCs:
stcs = mne.minimum_norm.apply_inverse_epochs(epos, inv, lambda2 = 1.0 / SNR ** 2.0, method = 'dSPM', return_generator=True)
res = linear_regression(stcs, epos.metadata[model], names=model)
for cond in model:
newStc = morph.apply(res[cond].beta)
subjBetas_src[cond].append(newStc)
print('Saving betas...')
with open(os.path.join(STCDir, language+'_'+formula+'_sns_250Hz_signed.pickle'), 'wb') as f:
pickle.dump(subjBetas, f)
if makeSTCs:
with open(os.path.join(STCDir, language+'_'+formula+'_src_250Hz_loose.pickle'), 'wb') as f:
pickle.dump(subjBetas_src, f)
print('Assigning subject betas to group-level betas...')
for cond in model:
betas[cond] += subjBetas[cond]
betas_src[cond] += subjBetas_src[cond]
print('Setting up for spatio-temporal analyses...')
try:
src_adjacency = mne.spatial_src_adjacency(src[:1])
src_adjacency_whole = mne.spatial_src_adjacency(src)
except:
srcFName = os.path.join(subjects_dir, 'fsaverage', 'bem', 'fsaverage-ico-4-src.fif')
src = mne.read_source_spaces(fname=srcFName)
src_adjacency = mne.spatial_src_adjacency(fsave_src[:1])
src_adjacency_whole = mne.spatial_src_adjacency(fsave_src)
sns_adjacency, ch_names = mne.channels.find_ch_adjacency(betas['Intercept'][0].info, 'mag')
print('Plotting data...')
for space in spaces:
note = '-' + str(searchTMin) + '-' + str(searchTMax) + 'ms' + '_' + space + '_p' + str(pThresh)
groupPlotFolder = os.path.join('/Volumes', 'falconlab', 'experiments', expName, 'stats', language, formula + '_N=' + str(len(subjects)) + note)
# groupPlotFolder = os.path.join(title, language, formula)
if not os.path.exists(groupPlotFolder):
os.makedirs(groupPlotFolder)
for cond in model:
if cond != 'Intercept' and cond != 'TrialNum' and 'Length' not in cond:
# plot key ROIs...
if space == 'sensor':
adjacency = sns_adjacency
searchSpace = 'allSensors'
searchSpace_Ixs = np.arange(0,208,1)
combined = mne.combine_evoked(betas[cond], weights='equal')
tmp = combined.plot(show=False)
tmp.axes[0].set_title(cond)
tmp.figure.set_size_inches(6,3)
tmp.savefig(os.path.join(groupPlotFolder, '%s_%s_Plot_%s' %(language, formula, cond)), dpi=300)
elif space == 'source_lh' or space == 'source_rh' or space == 'source_wholebrain':
# searchSpace_Ixs = [np.arange(0,2562,1), np.arange(0,2562,1)]
searchSpace_Ixs = np.arange(0,2562,1)
cortex = [(0.8, 0.8, 0.8), (0.55, 0.55, 0.55)]
adjacency = src_adjacency
searchSpace = 'wholebrain'
elif space == 'source_ltemp':
searchSpace = leftTemp
searchSpace_Ixs = leftTemp.get_vertices_used(np.arange(0,2562,1))
cortex = [(0.3, 0.3, 0.3), (0.0, 0.0, 0.0)]
elif space == 'source_rtemp':
searchSpace = rightTemp
searchSpace_Ixs = rightTemp.get_vertices_used(np.arange(0,2562,1))
cortex = [(0.3, 0.3, 0.3), (0.0, 0.0, 0.0)]
elif space == 'source_lhlgnet':
searchSpace = leftLngNet
searchSpace_Ixs = leftLngNet.get_vertices_used(np.arange(0,2562,1))
cortex = [(0.3, 0.3, 0.3), (0.0, 0.0, 0.0)]
elif space == 'source_rhlgnet':
searchSpace = rightLngNet
searchSpace_Ixs = rightLngNet.get_vertices_used(np.arange(0,2562,1))
cortex = [(0.3, 0.3, 0.3), (0.0, 0.0, 0.0)]
elif space == 'source_bilat':
searchSpace = bilatLngNet
searchSpace_Ixs_Lh = bilatLngNet.lh.get_vertices_used(np.arange(0,2562,1))
searchSpace_Ixs_Rh = bilatLngNet.rh.get_vertices_used(np.arange(0,2562,1))
cortex = [(0.3, 0.3, 0.3), (0.0, 0.0, 0.0)]
searchSpace_Ixs = np.concatenate([searchSpace_Ixs_Lh, searchSpace_Ixs_Rh+2562])
elif space == 'source_temp':
searchSpace = bilatTemp
searchSpace_Ixs_Lh = bilatTemp.lh.get_vertices_used(np.arange(0,2562,1))
searchSpace_Ixs_Rh = bilatTemp.rh.get_vertices_used(np.arange(0,2562,1))
cortex = [(0.3, 0.3, 0.3), (0.0, 0.0, 0.0)]
searchSpace_Ixs = np.concatenate([searchSpace_Ixs_Lh, searchSpace_Ixs_Rh+2562])
if space == 'source_lh' or space == 'source_lhlgnet' or space == 'source_ltemp':
hemi = 'lh'
searchSpace_Ixs_List = [searchSpace_Ixs, []]
elif space == 'source_rh' or space == 'source_rhlgnet' or space == 'source_rtemp':
hemi = 'rh'
searchSpace_Ixs_List = [[], searchSpace_Ixs]
elif space == 'source_wholebrain':
searchSpace_Ixs_Lh = np.arange(0,2562,1)
searchSpace_Ixs_Rh = np.arange(0,2562,1)
elif space == 'source_bilat' or space == 'source_wholebrain' or space == 'source_temp':
searchSpace_Ixs_List = [searchSpace_Ixs_Lh, searchSpace_Ixs_Rh]
if space == 'source_bilat' or space == 'source_temp':
tmp_adj = src_adjacency_whole.toarray()
allIxs = np.unique(np.concatenate([searchSpace_Ixs_Lh, 2562+searchSpace_Ixs_Rh]))
tmp_adj = tmp_adj[allIxs, :]
tmp_adj = tmp_adj[:, allIxs]
adjacency = scipy.sparse.coo_matrix(tmp_adj)
elif space != 'sensor':
tmp_adj = src_adjacency.toarray()
tmp_adj = tmp_adj[searchSpace_Ixs, :]
tmp_adj = tmp_adj[:, searchSpace_Ixs]
adjacency = scipy.sparse.coo_matrix(tmp_adj)
elif space == 'source_wholebrain':
allIxs = np.arange(0, 5124, 1)
seaerchSpace_Ixs_List = [np.arange(0,2562,1), np.arange(0,2562,1)]
adjacency = src_adjacency_whole
searchSpace_Ixs = allIxs
# Time indices that we are searching through
timeWindow_Ixs = times[int(searchTMin/decim):int(searchTMax/decim)]
print('Sensor space permutation tests...')
if space == 'sensor':
rawData = np.array([x.get_data() for x in betas[cond]])
elif space == 'source_lh' or space == 'source_lhlgnet' or space == 'source_ltemp':
rawData = np.array([x.lh_data for x in betas_src[cond]])
elif space == 'source_rh' or space == 'source_rhlgnet' or space == 'source_rtemp':
rawData = np.array([x.rh_data for x in betas_src[cond]])
elif space == 'source_bilat' or space == 'source_wholebrain' or space == 'source_temp':
rawData = np.array([x.data for x in betas_src[cond]])
# Original epoching including 0 and 1000ms,
# so decimation creates epochs of different sizes;
# we're reducing these all to 250 time points
if rawData.shape[2] > int(1000/decim):
rawData = rawData[:,:,:int(1000/decim)]
# subset to search times
data = rawData[:, :, int(searchTMin/decim):int(searchTMax/decim)]
# subset to search space
data = data[:, searchSpace_Ixs, :]
# rearrange data to subjects x times x spaces
data = np.transpose(data, [0, 2, 1])
# degrees of freedom and calculating t-threshold
df = len(subjects) - 1
tThresh = scipy.stats.distributions.t.ppf(1 - pThresh / 2, df = df)
print('One sample t-test over beta values for %s in %s' %(cond, space))
T0, clusters, clusterPs, H0 = clu = spatio_temporal_cluster_1samp_test(
data,
threshold = tThresh,
n_permutations = n_permutations,
tail = tail,
adjacency = adjacency,
verbose = True,
n_jobs = 30)
good_cluster_inds = np.where(clusterPs < plotThresh)[0]
for i_clu, clu_idx in enumerate(good_cluster_inds):
clusterInfo = ''
clusterTimes, clusterSpaces = clusters[clu_idx]
clusterPVal = round(clusterPs[clu_idx], 2)
clusterSize = round(sum(T0[clusters[clu_idx][0], clusters[clu_idx][1]]), 2)
clusterTime_Ixs = np.unique(clusterTimes)
clusterSpace_Ixs = np.unique(clusterSpaces)
# Averaging over time dimension for T values for spatial plot
#t_map = np.ma.masked_invalid(T0[clusterTime_Ixs, ...]).mean(0)
t_map = T0[clusterTime_Ixs, ...].mean(0)
if True:
# threshhold the t-values to t > 1.6
subTMap = np.where(abs(t_map) > 1.2)[0] # indices of
clusterSpace_Ixs = clusterSpace_Ixs[np.in1d(clusterSpace_Ixs, subTMap)]
# Subset to just the spatial points that are significant
sig_t_map = t_map[clusterSpace_Ixs]
# Test to see if you have positive and negative values and exclude
# those that aren't relevant
negPts = np.where(sig_t_map < 0)
posPts = np.where(sig_t_map > 0)
if negPts[0].shape[0] != 0 and posPts[0].shape[0] != 0:
print('Cluster has positive and negative t-values; removing the smallest')
if negPts[0].shape[0] > posPts[0].shape[0]:
sig_t_map = sig_t_map[negPts]
clusterSpace_Ixs = clusterSpace_Ixs[negPts]
else:
sig_t_map = sig_t_map[posPts]
clusterSpace_Ixs = clusterSpace_Ixs[posPts]
tValmin = np.percentile(sig_t_map, 2.5)
tValmean = np.percentile(sig_t_map, 50)
tValmax = np.percentile(sig_t_map, 97.5)
# For plotting purposes -- what's the minimum, mean, and maximum t-values that we'll plot?
if tValmin > 0 and tValmax > 0:
cmap = 'Reds'
clim = dict(kind="value", lims=[tValmin, tValmean, tValmax])
elif tValmin < 0 and tValmax < 0:
cmap = 'Blues'
sig_t_map = -1 * sig_t_map
clim = dict(kind="value", lims=[abs(tValmax), abs(tValmean), abs(tValmin)])
# elif 'sensor' in space:
# cmap = 'coolwarm'
# clim = dict(kind="value", lims=[-1*abs(tValmax), 0, abs(tValmax)])
else:
cmap = 'coolwarm'
clim = dict(kind="value", lims=[-1*max(abs(tValmin), abs(tValmax)), abs(tValmean), max(abs(tValmin), abs(tValmax))])
print(i_clu, clu_idx, cmap)
# Now we'll fetch the significant times:
sig_Times = timeWindow_Ixs[clusterTime_Ixs]
# For reporting, we'll want to know the MEG channel names
# / ico-4 src points
sig_Spaces = searchSpace_Ixs[clusterSpace_Ixs]
if space == 'sensor':
views = ['sensor']
elif 'wholebrain' in space:
views = ['lateral', 'ventral', 'medial']
else:
views = ['lateral']
for view in views:
print(view)
# fig = plt.subplots(figsize=(14, 3))
fig = plt.figure(figsize=(20, 3))
gs = plt.GridSpec(nrows=2, ncols=2, width_ratios=[1,3])
ax_topo = fig.add_subplot(gs[0,0])
ax_signals = fig.add_subplot(gs[0,1])
ax_bar = fig.add_subplot(gs[1,0])
ax_rawsigs = fig.add_subplot(gs[1,1])
plt.tight_layout()
if space == 'sensor':
mask = np.zeros((t_map.shape[0], 1), dtype=bool)
mask[sig_Spaces, :] = True
if cmap == 'Blues':
t_evoked = mne.EvokedArray(-1*t_map[:, np.newaxis], betas['Intercept'][0].info, tmin=0)
else:
t_evoked = mne.EvokedArray(t_map[:, np.newaxis], betas['Intercept'][0].info, tmin=0)
t_evoked.plot_topomap(times=0,
mask=mask,
axes=ax_topo,
cmap=cmap,
ch_type = 'mag',
scalings=dict(mag=1),
vlim=(clim['lims'][0], clim['lims'][2]),
#show=False,
#sensors=False,
extrapolate='local',
colorbar=False, mask_params=dict(markersize=9, markerfacecolor='#CCFFFF', markeredgecolor='#00CCCC'))
ax_topo.set_title("")
else:
if min(sig_Spaces) > 2562 and ('bilat' in space or '_temp' in space or 'wholebrain' in space):
hemis = ['rh']
sig_Spaces = sig_Spaces - 2562
if type(searchSpace) is mne.label.BiHemiLabel:
searchSpace = searchSpace.rh
elif min(sig_Spaces) > 2562:
print('hi!!')
hemis = ['rh']
else:
hemis = ['lh']
if type(searchSpace) is mne.label.BiHemiLabel:
searchSpace = searchSpace.lh
for hemi in hemis:
if hemi == 'lh':
t_stc = mne.SourceEstimate(sig_t_map, [sig_Spaces, []], 0, 0.001)
else:
t_stc = mne.SourceEstimate(sig_t_map, [[], sig_Spaces], 0, 0.001)
try:
# save label to file for cross-language comparisons
newLabel = mne.Label(vertices=sig_Spaces, hemi=hemi)
newLabel.fill(src=fsave_src)
#fName = '%s-%s-%s-%s-clNo%s-p=%s-%s_%s-label' %(language, formula, space, cond, str(i_clu), str(clusterPVal), cond, hemi)
fName = '%s-%s-%s-%s-clNo%s-p=%s-%s' %(language, formula, space, cond, str(i_clu), str(clusterPVal), cond)
newLabel.save(os.path.join(groupPlotFolder, fName))
brain = Brain(subject = 'fsaverage',
subjects_dir = subjects_dir,
hemi=hemi,
size = 1000,
surf = 'inflated',
views = view,
cortex = cortex,
background = 'white',
)
if 'wholebrain' not in space:
brain.add_label(searchSpace, color='white', alpha=0.8)
if hemi == 'lh':
tData = t_stc.lh_data
vertNo = t_stc.lh_vertno
else:
tData = t_stc.rh_data
vertNo = t_stc.rh_vertno
# if cmap != 'coolwarm':
brain.add_data(tData,
hemi=hemi,
vertices=vertNo,
colorbar = False,
clim = clim,
colormap = cmap,
smoothing_steps=8,
# fmin = clim['lims'][0],
# fmid = clim['lims'][1],
# fmax = clim['lims'][2],
transparent = True,
alpha=1.0
# # vertices = vertices[0]
)
# else:
# tData_neg = tData[np.where(tData < 0)[0]]
# tData_pos = tData[np.where(tData > 0)[0]]
# tValmin_pos = np.percentile(tData_pos, 2.5)
# tValmean_pos = np.percentile(tData_pos, 50)
# tValmax_pos = np.percentile(tData_pos, 97.5)
# tValmin_neg = np.percentile(tData_pos, 2.5)
# tValmean_neg = np.percentile(tData_pos, 50)
# tValmax_neg = np.percentile(tData_pos, 97.5)
# clim_pos = dict(kind="value", lims=[tValmin_pos, tValmean_pos, tValmax_pos])
# clim_neg = dict(kind="value", lims=[abs(tValmax_neg), abs(tValmean_neg), abs(tValmin_neg)])
# tData_neg = -1 * tData_neg
# brain.add_data(tData_pos,
# hemi=hemi,
# vertices=vertNo[np.where(tData > 0)[0]],
# colorbar = False,
# clim = clim_pos,
# colormap = 'Reds',
# smoothing_steps=8,
# # fmin = clim['lims'][0],
# # fmid = clim['lims'][1],
# # fmax = clim['lims'][2],
# transparent = True,
# alpha=1.0
# # # vertices = vertices[0]
# )
# brain.add_data(tData_neg,
# hemi=hemi,
# vertices=vertNo[np.where(tData < 0)[0]],
# colorbar = False,
# clim = clim_neg,
# colormap = 'Blues',
# smoothing_steps=8,
# # fmin = clim['lims'][0],
# # fmid = clim['lims'][1],
# # fmax = clim['lims'][2],
# transparent = True,
# alpha=1.0
# # # vertices = vertices[0]
# )
img = brain.screenshot()
nonwhite_pix = (img != 255).any(-1)
nonwhite_row = nonwhite_pix.any(1)
nonwhite_col = nonwhite_pix.any(0)
cropped_screenshot = img[nonwhite_row][:, nonwhite_col]
# ax_topo.imshow(cropped_screenshot
# ax_topo.axis('off')
ax_topo.imshow(cropped_screenshot)
ax_topo.set_yticklabels([])
ax_topo.set_xticklabels([])
ax_topo.set_yticks([])
ax_topo.set_xticks([])
ax_topo.spines[:].set_visible(False)
except:
print("Can't plot data in this hemisphere")
divider = make_axes_locatable(ax_topo)
ax_colorbar = divider.append_axes('left', size='3%', pad=0.2)
mne.viz.plot_brain_colorbar(ax_colorbar, clim=clim,
transparent=True, colormap=cmap,
label='', bgcolor='white')
if cmap == 'Blues':
ylabs = ax_colorbar.get_yticklabels()
for x in ylabs:
x.set_text('-' + x.get_text())
ax_colorbar.set_yticklabels(ylabs)
ax_topo.set_xlabel('Averaged t-map %s - %s ms' %(str(round(sig_Times[0])), str(round(sig_Times[-1]))))
#ax_signals = divider.append_axes('right', size='300%', pad=1.2)
title = 'v = {0}, p = {1}, {2} source'.format(clusterSize, clusterPVal, len(sig_Spaces))
if len(sig_Spaces) > 1:
title += "s (mean)"
# Take the raw data and subset to significant spatial indices
# And average over those spatial coodinates
# Then average over the participants and calculate SE
t_Timecourse = rawData[:, sig_Spaces, :]
t_Timecourse = np.average(t_Timecourse,axis=1)
t_tc_means = np.average(t_Timecourse,axis=0)
t_tc_upper = t_tc_means + 1.96*scipy.stats.sem(t_Timecourse, axis=0)
t_tc_lower = t_tc_means - 1.96*scipy.stats.sem(t_Timecourse, axis=0)
# Do this fore the baseline/intercept
#baseline_Timecourse = baselineData[:, sig_Spaces, :]
#baseline_Timecourse = np.average(baseline_Timecourse, axis=1)
#baseline_tc_means = np.average(baseline_Timecourse, axis=0)
#baseline_tc_upper = baseline_tc_means + 1.96*scipy.stats.sem(baseline_Timecourse, axis=0)
#baseline_tc_lower = baseline_tc_means - 1.96*scipy.stats.sem(baseline_Timecourse, axis=0)
#baseline_Avgs = baseline_Timecourse[:, clusterTime_Ixs]
#baseline_Avgs = baseline_Avgs.flatten()
#baseline_Avgs = np.average(baseline_Avgs, axis = 0)
# And other factorial variables
# Baseline/Intercept is 'NomAcc'
color='blue'
label=cond
ax_signals.plot(times, t_tc_means, color, label = label)
ax_signals.fill_between(times, y1=t_tc_upper, y2=t_tc_lower, color=color, alpha=0.1, interpolate=True)
ymin, ymax = ax_signals.get_ylim()
if space == 'sensor':
fig = tidyFig(fig, ax = 1, xmin=times[0], xmax=times[-1], size=(15,5))
else:
fig = tidyFig(fig, ax = 1, xmin=times[0], xmax=times[-1], ylim=ymin, size=(15,5))
ax_signals.fill_betweenx((-1*max(abs(ymin), abs(ymax)), max(abs(ymin), abs(ymax))), sig_Times[0], sig_Times[-1], color='lightblue', alpha=0.2, zorder=500)
ax_signals.fill_betweenx((-1*max(abs(ymin), abs(ymax)), max(abs(ymin), abs(ymax))), times[0], timeWindow_Ixs[0], color='grey', alpha=0.3)
ax_signals.fill_betweenx((-1*max(abs(ymin), abs(ymax)), max(abs(ymin), abs(ymax))), timeWindow_Ixs[-1], times[-1], color='grey', alpha=0.3)
ax_signals.hlines(y=0, color='black', xmin=times[0], xmax=times[-1])
ax_signals.set_xlabel('Times (ms)')
if space == 'sensor':
ax_signals.set_ylabel('Activity (β, fT)')
else:
ax_signals.set_ylabel('Activity (β, dSPM)')
ax_signals.set_title(title)
# fig.subplots_adjust(bottom=.05)
if cond == 'ObjCF':
comparisons = ['Acc', 'Bare']
colors = ['#128BFC', '#9C7535']
elif cond == 'SubjCF':
comparisons = ['NomImpf', 'ErgPerf']
colors = ['#128BFC', '#9C7535']
elif cond == 'InteractionF':
comparisons = ['ErgAcc', 'ErgBare', 'NomAcc', 'NomBare',]
colors = ['#128BFC', '#D89014', '#3B488A', '#9C7535']
elif cond == 'UpcomingVGenF':
comparisons = ['MVerb_Hindi', 'FVerb_Hindi']
colors = ['#D83F23', '#428669']
elif cond == 'GenF':
comparisons = ['MVerb_Nepali', 'FVerb_Nepali']
colors = ['#D83F23', '#428669']
elif cond == 'ClozeF':
comparisons = ['HighCloze', 'LowCloze']
colors = ['#4A5C3D', '#CF23D8']
# Palette
palette = dict()
# For ptcpt x space x time data; constrained to cluster spatial coordinates
data_3d = dict()
# For ptcpt x time data; averaged over space
data_tc = dict()
# SEs for plotting
data_tc_se = dict()
# For ptpct data; averaged over space and time
data_1d = dict()
for i in range(0,len(comparisons)):
color = colors[i]
comparison = comparisons[i]
palette[comparison] = color
for comparison in comparisons:
if space == 'sensor':
# participants x space x time
data_3d[comparison] = np.array([x.get_data() for x in raw_conds[comparison]])
# Subset to significant spaces in the cluster
data_3d[comparison] = data_3d[comparison][:, sig_Spaces, :]
else:
data_3d[comparison] = np.array([x.extract_label_time_course(newLabel, src=fsave_src, mode='mean') for x in raw_conds_src[comparison]])
# Clip off the extra ms if needed
if data_3d[comparison].shape[2] > int(1000/decim):
data_3d[comparison] = data_3d[comparison][:,:,:int(1000/decim)]
# Average over space and participants for the tiemecourse data
data_tc[comparison] = np.average(data_3d[comparison],axis=1)
# Calculate the SE before averaging over participants
data_tc_se[comparison] = scipy.stats.sem(data_tc[comparison], axis=0)
# Average across participants for the mean timecourse
data_tc[comparison] = np.average(data_tc[comparison], axis=0)
upper_tc = data_tc[comparison] + 1.96*data_tc_se[comparison]
lower_tc = data_tc[comparison] - 1.96*data_tc_se[comparison]
# bare_tc_lower = bareTC_means - 1.96*scipy.stats.sem(bareTC, axis=0)
ax_rawsigs.plot(times, data_tc[comparison], color = palette[comparison], label=comparison)
ax_rawsigs.fill_between(times, y1=upper_tc, y2=lower_tc, color=palette[comparison], alpha=0.1, interpolate=True)
ymin, ymax = ax_rawsigs.get_ylim()
ax_rawsigs.fill_betweenx((-1*max(abs(ymin), abs(ymax)), max(abs(ymin), abs(ymax))), sig_Times[0], sig_Times[-1], color='lightblue', alpha=0.2, zorder=500)
ax_rawsigs.fill_betweenx((-1*max(abs(ymin), abs(ymax)), max(abs(ymin), abs(ymax))), times[0], timeWindow_Ixs[0], color='grey', alpha=0.3)
ax_rawsigs.fill_betweenx((-1*max(abs(ymin), abs(ymax)), max(abs(ymin), abs(ymax))), timeWindow_Ixs[-1], times[-1], color='grey', alpha=0.3)
if space == 'sensor':
fig = tidyFig(fig, ax = 3, xmin=times[0], xmax=times[-1], size=(15,5))
else:
fig = tidyFig(fig, ax = 3, xmin=times[0], xmax=times[-1], ylim=ymin, size=(15,5))
ax_rawsigs.hlines(y=0, color='black', xmin=times[0], xmax=times[-1])
ax_rawsigs.set_xlabel('Times (ms)')
if space == 'sensor':
ax_rawsigs.set_ylabel('Activity (fT)')
else:
ax_rawsigs.set_ylabel('Activity (dSPM)')
ax_rawsigs.legend(loc='lower right')
if space == 'sensor' and len(comparisons) == 2:
std = comparisons[0]
contrast = comparisons[1]
stdRaw = mne.combine_evoked(raw_conds[std], weights=[1]*len(subjects))
contrastRaw = mne.combine_evoked(raw_conds[contrast], weights=[1]*len(subjects))
stdRaw = stdRaw.crop(tmin=tmin+sig_Times[0]/1000, tmax=tmin+sig_Times[-1]/1000)
contrastRaw = contrastRaw.crop(tmin=tmin+sig_Times[0]/1000, tmax=tmin+sig_Times[-1]/1000)
midPt = (sig_Times[0]/1000 + sig_Times[-1]/1000)/2
maxMidDiff = (sig_Times[-1]/1000) - midPt
diff = mne.combine_evoked([stdRaw, contrastRaw], weights=[1, -1])
diffData = diff.get_data().mean(axis=1)
diffEvoked = mne.EvokedArray(diffData[:, np.newaxis], betas['Intercept'][0].info, tmin=0)
diffMin = min(diffData)
diffMax = max(diffData)
diffMean = diffData.mean()
extremeVal = max(abs(diffMin), abs(diffMax))
clim = dict(kind="value", lims=[-extremeVal, diffMean, extremeVal])
diffEvoked.plot_topomap(times=0,
mask=mask,
axes=ax_bar,
cmap='coolwarm',
ch_type = 'mag',
scalings=dict(mag=1),
vlim=(-1*extremeVal, 1*extremeVal),
#show=False,
#sensors=False,
extrapolate='local',
colorbar=False, mask_params=dict(markersize=9, markerfacecolor='#CCFFFF', markeredgecolor='#00CCCC'))
ax_bar.set_title("")
ax_bar.set_xlabel("Average signal during cluster: %s - %s" %(std, contrast))
divider = make_axes_locatable(ax_bar)
ax_colorbar2 = divider.append_axes('left', size='3%', pad=0.2)
mne.viz.plot_brain_colorbar(ax_colorbar2, clim=clim,
transparent=True, colormap='coolwarm',
label='', bgcolor='white')
# We'll plot the difference in the MEG topography
else:
# We'll plot the bar graph in the cluster
# The dataframe for plotting in seaborn
avgMeans = dict()
for comparison in comparisons:
# Data has already been constrainted to
# spatial coordinates, so now we will constrain to
# time indices
tmp = data_3d[comparison][:, :, clusterTime_Ixs]
# Then we will average over space...
tmp = np.average(tmp, axis=1)
# ... and time
tmp = np.average(tmp, axis=1)
data_1d[comparison] = tmp
avgMeans[comparison] = data_1d[comparison]
sn.barplot(avgMeans, ax=ax_bar, palette=palette, width=0.8)
# ax_bar.set_
# sn.barplot(avgMeans, x='variable', y='value', ax=ax_bar, palette=colors, hue='variable', width=1, errorbar='se')
avgs = []
for x in list(avgMeans.keys()):
avgs.append(avgMeans[x].mean())
biggestVal = max(avgs)
scale = (max(0,abs(ymin))+min(abs(ymax), biggestVal))/100
yIncrement = 1.1
for i_fact in range(0,len(list(avgMeans.keys()))):
for j_fact in range(0,len(list(avgMeans.keys()))):
if i_fact < j_fact:
ilab = list(avgMeans.keys())[i_fact]
jlab = list(avgMeans.keys())[j_fact]
tval, pval = scipy.stats.ttest_ind(a = avgMeans[ilab], b = avgMeans[jlab])
tTestRes = '%s vs. %s: t = %s, p = %s\n' %(ilab, jlab, str(tval), str(pval))
clusterInfo += tTestRes
iavg = avgMeans[ilab].mean()
javg = avgMeans[jlab].mean()
biggestVal = max(abs(biggestVal), abs(iavg), abs(javg))
if iavg <= 0 and javg <= 0:
yval = min(iavg, javg, -1*biggestVal) - 1.0*((yIncrement-1)*scale)
else:
yval = max(iavg, javg, biggestVal) + 1.0*((yIncrement-1)*scale)
if pval < 0.05:
yIncrement += 1.1
ax_bar.plot([i_fact, i_fact, j_fact, j_fact], [yval, yval+0.1*yval, yval+0.1*yval, yval], lw=1.5, c='black')
ax_bar.text((i_fact+j_fact) * 0.5, yval, '*', ha='center', va='center', color='black', fontsize=20)
avgs.append(yval)
avgs.append(biggestVal)
ax_bar.spines['top'].set_visible(False)
ax_bar.spines['right'].set_visible(False)
ax_bar.spines['bottom'].set_visible(False)
ax_bar.title.set_visible(False)
# ax_bar.set_xticks([])
# ax_bar.xaxis.set_major_formatter(plt.NullFormatter())
ax_bar.set_xticklabels(comparisons)
# newMin, newMax = ax_bar.get_ylim()
# if sum(1 for number in avgs if number < 0) == 0:
# ax_bar.set_ylim(ymin = max(0, ymin, newMin), ymax = min(max(avgs), newMax)+scale)
# else:
# ax_bar.set_ylim(ymin = max(ymin, newMin), ymax = min(max(avgs), newMax)+scale)
fig.subplots_adjust(bottom=.05)
# label_params = ax_signals.get_legend_handles_labels()
# figl, axl = plt.subplots()
# axl.axis(False)
# axl.legend(*label_params, ncols = 2, loc="center", bbox_to_anchor=(0.5, 0.5))
# fName = '%s-%s-%s-%s-clNo%s-p=%s-%s_LABEL_ONLY.png' %(language, formula, space, cond, str(i_clu), str(clusterPVal), cond)
# figl.savefig(os.path.join(groupPlotFolder, fName + "_LABEL_ONLY.png"))
fName = '%s-%s-%s-%s-clNo%s-p=%s-%s_%s.png' %(language, formula, space, cond, str(i_clu), str(clusterPVal), cond, view)
fig.savefig(os.path.join(groupPlotFolder, fName))
clusterInfo += 'formula: ' + formula + '\n'
clusterInfo += 'effect: ' + cond + '\n'
clusterInfo += 'n_permutations: ' + str(n_permutations) + '\n'
clusterInfo += 'pThresh: ' + str(pThresh) + '\n'
clusterInfo += 'plotThresh: ' + str(plotThresh) + '\n'
clusterInfo += 'timeWindow_Ixs: ' + str(timeWindow_Ixs) + '\n'
clusterInfo += 'space: ' + str(space) + '\n'
clusterInfo += 'searchSpace_Ixs: ' + str(searchSpace_Ixs) + '\n'
clusterInfo += 'i_clu: ' + str(i_clu) + '\n'
clusterInfo += 'clusterPVal: ' + str(clusterPVal) + '\n'
clusterInfo += 'sig_Times: ' + str(sig_Times) + '\n'
clusterInfo += 'sig_Spaces: ' + str(sig_Spaces) + '\n'
fName = '%s-%s-%s-%s-clNo%s-p=%s-%s_info.txt' %(language, formula, space, cond, str(i_clu), str(clusterPVal), cond)
with open(os.path.join(groupPlotFolder, fName), 'w') as f:
f.write(clusterInfo)

1_SgTrialRegressions_withMean_looseOrientation.py, no license · at the source

Overview

Authors: Dustin A. Chacón^1,2, Subhekshya Shrestha^1,2, Brian W. Dillon³, Rajesh Bhatt³, Diogo Almeida², Alec Marantz^2,4

ORCID iDs: Dustin A. Chacón

Institutions: New York University Abu Dhabi (United Arab Emirates); University of California, Santa Cruz (United States); University of Massachusetts Amherst (United States); New York University (United States)

Journal: Neurobiology of language (Cambridge, Mass.), volume 7, article NOL.a.247

Dates: received 12 February 2024; accepted 19 February 2026; published online 4 June 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1162/nol.a.247 · PMID 42272623 · PMCID PMC13249467 · OpenAlex W7133501173

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: MEG (modality)

Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Preprocessing, fMRI & imaging, Physiology & signal measures

Keywords: agreement, case, electrophysiology, ergativity, Hindi, magnetoencephalography (MEG), Nepali, neurolinguistics, psycholinguistics, sentence processing

Topic: Neurobiology of Language and Bilingualism (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: Research Institute Centers, New York University Abu Dhabi (G1001)

Citations: not cited yet (Europe PMC); 125 references in the paper

Abstract

At first glance, the brain’s language network appears to be universal, but languages clearly differ. Does the brain adapt to the specific details of individual grammatical systems? Here, we present a magnetoencephalography (MEG) study on case and agreement in Hindi and Nepali. Both languages use split-ergative case systems. However, these systems interact with verb agreement differently—in Hindi, case features conspire to determine which noun phrase (NP) the verb agrees with (subject, object, or neither), but in Nepali the verb always agrees with the subject NP. We found that NPs with different case values elicit different MEG signals around 200–500 and 600–900 ms. In subsequent exploratory analyses, we failed to find a reliable difference in this brain activity between the two languages corresponding to the different relations between case and agreement. However, we identified a portion of the left temporoparietal junction as exhibiting a statistically nonsignificant effect that may warrant further investigation.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

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

OSF qrzpm

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Size: 13 files

Software Heritage: not checked

Found in: “DATA AND CODE AVAILABILITY STATEMENTS”

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

Tools: MNE-Python (3 files), NumPy (3 files), pandas (3 files), Matplotlib (2 files), SciPy (2 files), seaborn (2 files), statsmodels (1 file)

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

27 September 2026: the link answers (HTTP 200)

4 files

scripts.zip/scripts/0_Preprocessing.py, Python, 276 lines
scripts.zip/scripts/1_SgTrialRegressions_withMean_looseOrientation.py, Python, 1,336 lines, 5 matches
scripts.zip/scripts/2_ByLgTimecourses.py, Python, 319 lines
scripts.zip/scripts/params.py, Python, 142 lines

At the source: osf.io/qrzpm

The paper's code and data availability statement is in the Data section.

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:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
4 scripts, each with its path and the digest of its content;
5 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

No dataset and no data link were found in the paper.

Data and code availability statements

Preprocessed MEG data and processing code have been made publicly available on Open Science Framework at https://osf.io/qrzpm.

Reproduced under the paper's license (CC BY), from the paper cited above.

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, pages, dates, 6 authors, 10 keywords, 1 funder, 116 references.

Cite

This paper

Chacón, D. A., Shrestha, S., Dillon, B. W., Bhatt, R., Almeida, D., & Marantz, A. (2026). Same Sentences, Different Grammars, Different Brain Responses?: An MEG Study on Case and Agreement Encoding in Hindi and Nepali Split-Ergative Structures. Neurobiology of language (Cambridge, Mass.), 7, NOL.a.247. https://doi.org/10.1162/nol.a.247

BibTeX

@article{chacon2026same, author = {Chacón, Dustin A. and Shrestha, Subhekshya and Dillon, Brian W. and Bhatt, Rajesh and Almeida, Diogo and Marantz, Alec}, title = {{Same Sentences, Different Grammars, Different Brain Responses?: An MEG Study on Case and Agreement Encoding in Hindi and Nepali Split-Ergative Structures}}, journal = {Neurobiology of language (Cambridge, Mass.)}, year = {2026}, month = jun, volume = {7}, pages = {NOL.a.247}, publisher = {MIT Press}, issn = {2641-4368}, doi = {10.1162/nol.a.247}, url = {https://doi.org/10.1162/nol.a.247}, pmid = {42272623}, pmcid = {PMC13249467} }

RIS

TY - JOUR AU - Chacón, Dustin A. AU - Shrestha, Subhekshya AU - Dillon, Brian W. AU - Bhatt, Rajesh AU - Almeida, Diogo AU - Marantz, Alec TI - Same Sentences, Different Grammars, Different Brain Responses?: An MEG Study on Case and Agreement Encoding in Hindi and Nepali Split-Ergative Structures T2 - Neurobiology of language (Cambridge, Mass.) J2 - Neurobiol Lang (Camb) PY - 2026 DA - 2026/06/04 VL - 7 SP - NOL.a.247 SN - 2641-4368 PB - MIT Press DO - 10.1162/nol.a.247 UR - https://doi.org/10.1162/nol.a.247 LA - en ER -

CSL-JSON

{ "id": "10.1162/nol.a.247", "type": "article-journal", "title": "Same Sentences, Different Grammars, Different Brain Responses?: An MEG Study on Case and Agreement Encoding in Hindi and Nepali Split-Ergative Structures", "container-title": "Neurobiology of language (Cambridge, Mass.)", "author": [ { "family": "Chacón", "given": "Dustin A." }, { "family": "Shrestha", "given": "Subhekshya" }, { "family": "Dillon", "given": "Brian W." }, { "family": "Bhatt", "given": "Rajesh" }, { "family": "Almeida", "given": "Diogo" }, { "family": "Marantz", "given": "Alec" } ], "container-title-short": "Neurobiol Lang (Camb)", "volume": "7", "page": "NOL.a.247", "DOI": "10.1162/nol.a.247", "PMID": "42272623", "PMCID": "PMC13249467", "ISSN": "2641-4368", "publisher": "MIT Press", "URL": "https://doi.org/10.1162/nol.a.247", "language": "en", "issued": { "date-parts": [ [ 2026, 6, 4 ] ] } }

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.