In [1]:
from __future__ import division
import os.path as op
import itertools
import numpy as np
import scipy as sp
import pandas as pd
import nibabel as nib
from scipy import stats
import statsmodels.api as sm
import matplotlib as mpl
import matplotlib.pyplot as plt
import nitime as nit
In [2]:
import lyman
from lyman import mvpa, evoked
import seaborn as sns
import moss
%matplotlib inline
sns.set(context="paper", style="ticks", font="Arial")
mpl.rcParams.update({"xtick.major.width": 1, "ytick.major.width": 1, "savefig.dpi": 150})
pd.set_option('display.precision', 3)
In [3]:
from IPython.parallel import Client, TimeoutError
try:
dv = Client()[:]
dv4 = Client()[:4]
except (IOError, TimeoutError):
dv = None
dv4 = None
import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning)
In [69]:
def calculate_evoked_ts(data, n_bins, problem=None, events=None, tr=2,
calc_method="FIR", offset=0, upsample=1,
percent_change=True, correct_baseline=True,
event_names=None):
project = lyman.gather_project_info()
design_template = op.join(project["data_dir"], "%s",
"design/%s.csv" % problem)
col_names = ['subj', 'event', 'tr', 'response']
evoked_ds = pd.DataFrame(columns = col_names)
for i, data_i in enumerate(data):
# Get event information
subj = data_i["subj"]
events_i = pd.read_csv(design_template % subj)
# Map from event names to integer index values
event_names = sorted(events_i.condition.unique())
event_map = pd.Series(range(1, len(event_names) + 1),
index=event_names)
# Create the timeseries of event occurances
calc_tr = float(tr) / upsample
event_list = []
data_list = []
for run, run_data in enumerate(data_i["data"], 1):
# Possibly upsample the data
if upsample != 1:
time_points = len(run_data)
x = np.linspace(0, time_points - 1, time_points)
xx = np.linspace(0, time_points,
time_points * upsample + 1)[:-upsample]
interpolator = interp1d(x, run_data, "cubic", axis=0)
run_data = interpolator(xx)
run_events = events_i[events_i.run == run]
run_events.onset += offset
event_id = np.zeros(len(run_data), int)
event_index = np.array(run_events.onset / calc_tr).astype(int)
event_id[event_index] = run_events.condition.map(event_map)
event_list.append(event_id)
if percent_change:
run_data = nit.utils.percent_change(run_data, ax=0)
data_list.append(run_data)
# Set up the Nitime objects
event_info = np.concatenate(event_list)
data = np.concatenate(data_list, axis=0)
# Do the calculations
calc_bins = n_bins * upsample
if data.ndim == 1:
evoked_data = evoked._evoked_1d(data, event_info, calc_bins, calc_tr,
calc_method, correct_baseline)
elif data.ndim == 2:
evoked_data = evoked._evoked_2d(data, event_info, n_bins, calc_tr,
calc_method, correct_baseline)
# evoked_ds.append(evoked_data)
for event_num, d in enumerate(evoked_data):
evoked_ds = evoked_ds.append(pd.DataFrame({'subj': subj,
'event': event_names[event_num],
'tr': np.arange(1, n_bins+1),
'response': pd.Series(evoked_data[event_num])}))
# return np.array(evoked_ds).squeeze()
return evoked_ds
In [4]:
class Args(object):
def __init__(self, **kwargs):
self.__dict__.update(kwargs)
args = Args(subjects=['/Users/steph-backup/Experiments/ObjFam/data/subids_subset_no23or19.txt'],
experiment = 'objfam',
level = 'group',
altmodel = 'submem')
In [5]:
project = lyman.gather_project_info()
exp = lyman.gather_experiment_info(args.experiment, args.altmodel)
anal_dir = project["analysis_dir"]
data_dir = project["data_dir"]
subjects = pd.Series(lyman.determine_subjects(args.subjects), name="subj")
all_rois = ['lh.lateraloccipital', 'rh.lateraloccipital']
In [ ]:
exp
In [6]:
def objfam_roi_sizes():
roi_sizes = pd.DataFrame(columns=all_rois, index=subjects, dtype=float)
mask_template = op.join(data_dir, "%s/masks/%s.nii.gz")
for roi in roi_sizes:
for subj in subjects:
vox_count = nib.load(mask_template % (subj, roi.lower())).get_data().sum()
roi_sizes.loc[subj, roi] = vox_count
return roi_sizes.describe().loc[["mean", "std"]].astype(int)
objfam_roi_sizes()
Out[6]:
In [7]:
def objfam_roi_signal():
roi_signal = pd.DataFrame(columns=all_rois, index=subjects, dtype=float)
for roi in all_rois:
print roi.lower()
signal = evoked.extract_group(roi.lower(), dv=dv4)
# signal = [np.concatenate(d["data"]).mean() for d in signal]
# roi_signal[roi] = signal
# return (roi_signal / 100).describe().loc[["mean", "std"]]
objfam_roi_signal()
In [208]:
def extract_subject(subj, mask_name, experiment, altmodel):
exp_name = experiment+'-'+altmodel
d = evoked.extract_subject(subj, mask_name, exp_name=exp_name)
return d
In [166]:
signal = extract_subject('subj12', 'lh.lateraloccipital', args.experiment, args.altmodel)
sns.tsplot(signal['data'], color='dodgerblue')
Out[166]:
In [8]:
def extract_group(subjects, mask_name):
exp_name = args.experiment+'-'+args.altmodel
d = evoked.extract_group(mask_name, subjects=subjects, exp_name=exp_name, dv=dv4)
return d
In [227]:
data = extract_group(subjects, 'Left-Cerebral-Cortex-hit-fa')
In [230]:
n_bins = 12
signal = calculate_evoked_ts(data, n_bins,
problem=args.altmodel,
tr=int(exp['TR']))
In [231]:
ev_cols = pd.DataFrame(signal.event.str.split('_').tolist(), columns = ['morph', 'resp'])
signal = signal.reset_index(); signal = signal.join(ev_cols)
data_subset = signal.query('resp != "guess" & morph != "Test"')
palette_resp = dict(new='dodgerblue', old='mediumseagreen')
f, axes = plt.subplots(1, len(data_subset.morph.unique()),
figsize=(10, 5),
sharex=True, sharey=True)
for ax, morph in zip(axes, signal.morph.unique()):
sns.tsplot(data=data_subset[data_subset.morph==morph], time='tr',
unit='subj', ci=68,
condition='resp', value='response',
estimator=np.nanmean, ax=ax, color=palette_resp)
ax.set_title('morph ' + morph)
# filename = op.join(project['analysis_dir'],
# args.experiment+'-'+args.altmodel, 'plots',
# 'peristim_mask_rh_hit-miss.png')
# plt.savefig(filename, dpi=300)
In [232]:
data_avg = data_subset.query('2 < tr < 6').groupby(['subj', 'morph', 'resp']).mean().reset_index()
data_avg.head()
Out[232]:
In [234]:
sns.factorplot('morph', y='response',
hue='resp', units='subj',
kind='point',
ci=68, data=data_avg, dodge=0.1)
Out[234]:
In [100]:
data = extract_group(subjects, 'lateraloccipital-linear-study')
n_bins = 12
signal_linstudy = calculate_evoked_ts(data, n_bins,
problem=args.altmodel,
tr=int(exp['TR']))
In [101]:
signal_subset = signal_linstudy[signal_linstudy.event.isin(['1_old', '1_new'])]
# signal_subset = signal
sns.tsplot(signal_subset, time='tr', unit='subj',
condition='event', value='response')
Out[101]:
In [ ]: