import os import pandas as pd import numpy as np import matplotlib.pylab as plt import seaborn as sns import mr_utilities from scipy import stats mycol = ['darkgreen', 'mediumaquamarine', 'silver', 'darkblue','steelblue'] def cm2inch(value): return value/2.54 def mt_examples(patient, save_path, fit_method, kmin, kmax, windowsize, windowstep, dt): recordings = mr_utilities.list_preictrecordings(patient)[0] fig, ax = plt.subplots(2, 1, figsize=(cm2inch(13.5), cm2inch(5))) binnings = ['right', 'left'] quantity = 'm' t_onset = 600 # seizure onset at the end of each 10 min recording for irec, rec in enumerate(recordings): for ibin, binning in enumerate(binnings): mt_file = '{}{}/mt_results_{}_{}_kmin{}_kmax{}_winsize{}_' \ 'winstep{}_dt{}.pkl'.format(save_path, rec, binning, fit_method, kmin, kmax, windowsize, windowstep, dt) if not os.path.isfile(mt_file): print('{} not existing'.format(mt_file)) continue # -----------------------------------------------------------------# # read m(t) data # -----------------------------------------------------------------# mt_frame = pd.read_pickle(mt_file) mt_frame = mr_utilities.exclude_inconsistencies(mt_frame) if mt_frame.empty: continue kmin = mt_frame.kmin.unique()[0] kmax = mt_frame.kmax.unique()[0] fit_method = mt_frame.fit_method.unique()[0] winsize = mt_frame.windowsize.unique()[0] t0_list = mt_frame.t0.tolist() qt = mt_frame[quantity].tolist() t0_list = [t0_ms * 0.001 for t0_ms in t0_list] half_win_sec = int((winsize / 2) / 1000 * dt) t_middle = [t0 + half_win_sec for t0 in t0_list] # -----------------------------------------------------------------# # Plotting # -----------------------------------------------------------------# ax[ibin].plot(t_middle, qt, '-', linewidth=1, label='Rec {}'.format(irec+1), color=mycol[irec]) ax[ibin].axhline(y=1, color='black', linestyle='--', alpha=0.3) ax[ibin].set_ylabel(r"$\hat{m}$") ax[ibin].set_xlim( (t0_list[0] - 10, t0_list[-1] + 2 * half_win_sec + 10)) binning_focus = mr_utilities.binning_focus(binning, int(patient)) if binning_focus == 'focal': focus_label = 'ipsi' else: focus_label = 'contra' ax[ibin].set_title("{}".format(focus_label), loc='left') ax[0].legend(loc='center left', bbox_to_anchor=(1.04, -0.2), fontsize='small') ax[0].xaxis.grid(False) ax[1].xaxis.grid(False) ax[1].set_xlabel("time (s)") ax[0].annotate("", xy=(t_onset, 1), xytext=(t_onset, 1.1), arrowprops=dict(arrowstyle="->", color='black')) ax[1].annotate("", xy=(t_onset, 1), xytext=(t_onset, 1.04), arrowprops=dict(arrowstyle="->", color='black')) ax[0].annotate("seizure\n onset", xy=(t_onset, 1), xytext=(t_onset, 1.03)) sns.despine() plt.subplots_adjust(bottom=None, right=None, wspace=0.3, hspace=1.3, left=None, top=None) plt.savefig('{}/Fig2_mt_patient{}.pdf'.format(save_path, patient), bbox_inches='tight') plt.close() def boxplot_mx(result_path, x, ksteps, fit_method, Engel=False): # -----------------------------------------------------------------# # Read and filter data # -----------------------------------------------------------------# mx_results = '{}mt_{}_results_allrecords_kmin{}_kmax{}_{}.pkl'.format( result_path, x, ksteps[0], ksteps[1], fit_method) mx_allresults = pd.read_pickle(mx_results) mx_allresults = mr_utilities.exclude_inconsistencies(mx_allresults) x = mx_allresults['partno'].max() mx_allresults['logepsilon'] = mx_allresults.apply(lambda row: np.log( 1-row['m']), axis=1) mx_allresults['epsilon'] = mx_allresults.apply(lambda row: 1-row['m'], axis=1) measure = 'epsilon' patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19] if Engel: mx_allresults = mx_allresults[(mx_allresults['patient'].isin(patients1a))] # -----------------------------------------------------------------# # write estimates to lists # -----------------------------------------------------------------# focal_diff_list = [] contralateral_diff_list = [] all_last_contra = [] all_first_contra = [] all_last_focal = [] all_first_focal = [] for ir, recording in enumerate(mx_allresults['recording'].unique()): mx_recording = mx_allresults[mx_allresults['recording'] == recording] mx_focal = mx_recording[mx_recording['binning_focus'] == 'focal'] mx_contralateral = mx_recording[ mx_recording['binning_focus'] == 'contra-lateral'] if not mx_focal[mx_focal['partno'] == 1].empty and not mx_focal[ mx_focal['partno'] == x].empty: m_focal_first = mx_focal[mx_focal['partno'] == 1][ measure].item() m_focal_last = mx_focal[mx_focal['partno'] == x][measure].item() focal_diff_log = np.log(m_focal_last) - np.log(m_focal_first) focal_diff_list.append(focal_diff_log) all_first_focal.append(m_focal_first) all_last_focal.append(m_focal_last) if not mx_contralateral[mx_contralateral['partno'] == 1].empty \ and not mx_contralateral[ mx_contralateral['partno'] == x].empty: m_contralateral_first = \ mx_contralateral[mx_contralateral['partno'] == 1][measure].item() m_contralateral_last = \ mx_contralateral[mx_contralateral['partno'] == x][measure].item() contralateral_diff_log = np.log(m_contralateral_last) - \ np.log(m_contralateral_first) contralateral_diff_list.append(contralateral_diff_log) all_first_contra.append(m_contralateral_first) all_last_contra.append(m_contralateral_last) # -----------------------------------------------------------------# # Wilcoxon paired test # -----------------------------------------------------------------# stat, p_contra_w = stats.wilcoxon(contralateral_diff_list) stat, p_focal_w = stats.wilcoxon(focal_diff_list) # -----------------------------------------------------------------# # Create dataframe for plotting # -----------------------------------------------------------------# quantity = measure data1 = { quantity: all_first_focal, 'rec_type': ['10 - 5 min' for _ in range(len(all_first_focal))], 'binning_focus': ['SOZ' for _ in range(len(all_first_focal))] } df = pd.DataFrame(data1) data2 = pd.DataFrame({ quantity: all_last_focal, 'rec_type': ['5 - 0 min' for _ in range(len(all_last_focal))], 'binning_focus': ['SOZ' for _ in range(len(all_last_focal))] }) df = df.append(data2, ignore_index=True) data3 = pd.DataFrame({ quantity: all_first_contra, 'rec_type': ['10 - 5 min' for _ in range(len(all_first_contra))], 'binning_focus': ['nSOZ' for _ in range(len(all_first_contra))] }) df = df.append(data3, ignore_index=True) data4 = pd.DataFrame({ quantity: all_last_contra, 'rec_type': ['5 - 0 min' for _ in range(len(all_last_contra))], 'binning_focus': ['nSOZ' for _ in range(len(all_last_contra))] }) df = df.append(data4, ignore_index=True) # -----------------------------------------------------------------# # Plotting # -----------------------------------------------------------------# fig, ax = plt.subplots(1, 1, figsize=(cm2inch(5.5), cm2inch(4.5))) g = sns.boxplot(ax=ax, x='binning_focus', y=quantity, hue='rec_type', data=df, hue_order = ['10 - 5 min', '5 - 0 min'], order=[ 'SOZ', 'nSOZ'], dodge=True, color='slategray', palette=[ 'lightgray', 'slategray'], linewidth=1.1, width=0.6) g.set_yscale("log") plt.grid(axis='x', b=False) plt.rcParams['legend.title_fontsize'] = 9 plt.legend(loc="upper left", ncol=1, fontsize='x-small', bbox_to_anchor=(0, 1.7), title='time to ' 'seizure onset') n_f_pre = len(all_first_focal) n_f_inter = len(all_last_focal) n_c_pre = len(all_first_contra) n_c_inter = len(all_last_contra) ax.set_xticklabels(['ipsi\n({}, {})'.format(n_f_pre, n_f_inter), 'contra\n({}, {})'.format(n_c_pre, n_c_inter)]) # -----------------------------------------------------------------# # quantity labels # -----------------------------------------------------------------# quantity_label = r'$\hat{m}$' ymax = 1e-3 ymin = 1 ax.set_ylim(ymin, ymax) T = [1, 1e-1, 1e-2, 1e-3] ax.set_yticks(T) ax.set_yticklabels([0, 0.9, 0.99, 0.999]) # -----------------------------------------------------------------# # annotate p-values # -----------------------------------------------------------------# y_max = 2.4*ymax ax.annotate("", xy=(-0.2, y_max), xycoords='data', xytext=(0.2, y_max), textcoords='data', arrowprops=dict(arrowstyle="-", ec='black', connectionstyle="bar,fraction=0.2")) ax.text(-0.25, y_max*0.6, 'p={:.3f}'.format(p_focal_w), fontsize=8) ax.annotate("", xy=(0.8, y_max), xycoords='data', xytext=(1.2, y_max), textcoords='data', arrowprops=dict(arrowstyle="-", ec='black', connectionstyle="bar,fraction=0.2")) ax.text(0.75, y_max*0.6, 'p={:.3f}'.format(p_contra_w), fontsize=8) ax.set_ylim(ymin, ymax) plt.ylabel(quantity_label) plt.xlabel('') sns.despine() # -----------------------------------------------------------------# # Saving # -----------------------------------------------------------------# plotfile_mx = '{}Fig2_2parts_Engel{}.pdf'.format(result_path, Engel) plt.savefig(plotfile_mx, bbox_inches='tight') def pointplot_mx_regions_Engel(result_path, x, ksteps, fit_method): # -----------------------------------------------------------------# # Read and filter data # -----------------------------------------------------------------# mx_results = '{}mt_{}_results_allrecords_kmin{}_kmax{}_{}.pkl'.format( result_path, x, ksteps[0], ksteps[1], fit_method) mx_allresults = pd.read_pickle(mx_results) mx_allresults = mr_utilities.exclude_inconsistencies(mx_allresults) x = mx_allresults['partno'].max() mx_allresults['epsilon'] = mx_allresults.apply(lambda row: 1-row['m'], axis=1) measure = 'epsilon' def give_region(row): return row['binning'][1:] mx_allresults['region'] = mx_allresults.apply(lambda row: give_region(row), axis=1) # -----------------------------------------------------------------# # write paired estimates to lists # -----------------------------------------------------------------# fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(4*1.2))) df = None patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19] mx_allresults = mx_allresults[(mx_allresults['patient'].isin(patients1a))] regions = ['H', 'A', 'PHC', 'EC'] labels = [] for i, reg in enumerate(regions): mx_region = mx_allresults[mx_allresults['region'] == reg] recordings = mx_region.recording.unique() all_first_focal = [] all_last_focal = [] all_first_contra = [] all_last_contra = [] recID = [] for recording in recordings: mx_recording = mx_region[mx_region['recording'] == recording] mx_focal = mx_recording[mx_recording['binning_focus'] == 'focal'] mx_contralateral = mx_recording[ mx_recording['binning_focus'] == 'contra-lateral'] if not mx_focal[mx_focal['partno'] == 1].empty and not mx_focal[ mx_focal['partno'] == x].empty: m_focal_first = mx_focal[mx_focal['partno'] == 1][ measure].item() m_focal_last = mx_focal[mx_focal['partno'] == x][measure].item() all_first_focal.append(m_focal_first) all_last_focal.append(m_focal_last) recID.append(recording) if not mx_contralateral[mx_contralateral['partno'] == 1].empty \ and not \ mx_contralateral[ mx_contralateral['partno'] == x].empty: m_contralateral_first = mx_contralateral[mx_contralateral[ 'partno'] == 1][measure].item() m_contralateral_last = mx_contralateral[mx_contralateral[ 'partno'] == x][measure].item() all_first_contra.append(m_contralateral_first) all_last_contra.append(m_contralateral_last) # -----------------------------------------------------------------# # Create dataframe for plotting # -----------------------------------------------------------------# quantity = measure if df is None: data1 = { 'region': reg, quantity: all_first_focal, 'recording': recID, 'rec_type': ['10 - 5 min' for _ in range(len(all_first_focal))], 'binning_focus': ['SOZ' for _ in range(len(all_first_focal))] } df = pd.DataFrame(data1) else: data1 = pd.DataFrame({ 'region': reg, quantity: all_first_focal, 'recording': recID, 'rec_type': ['10 - 5 min' for _ in range(len(all_last_focal))], 'binning_focus': ['SOZ' for _ in range(len(all_last_focal))] }) df = df.append(data1, ignore_index=True) data2 = pd.DataFrame({ 'region': reg, quantity: all_last_focal, 'recording': recID, 'rec_type': ['5 - 0 min' for _ in range(len(all_last_focal))], 'binning_focus': ['SOZ' for _ in range(len(all_last_focal))] }) df = df.append(data2, ignore_index=True) n_f_first = len(all_first_focal) n_f_last = len(all_last_focal) labels.append('{}\n({}, {})'.format(reg, n_f_first, n_f_last)) x0 = (-0.2+i) * np.ones(len(all_first_focal)) x1 = (0.2+i) * np.ones(len(all_last_focal)) plt.plot([x0, x1], [all_first_focal, all_last_focal], color='grey', linewidth=0.5, linestyle='--') # -----------------------------------------------------------------# # Plotting # -----------------------------------------------------------------# df = df[df['binning_focus'] == 'SOZ'] g = sns.stripplot(ax=ax, x='region', y=quantity, hue='rec_type', data=df, hue_order=['10 - 5 min', '5 - 0 min'], dodge=True, order=regions, jitter=False, color='slategray', palette=['slategray', 'darkblue'], size=2.5, alpha=0.7) ax.set_yscale("log") ax.grid(axis='x', b=False) plt.rcParams['legend.title_fontsize'] = 9 plt.legend(loc='center left', bbox_to_anchor=(1, 0.5), fontsize='x-small', title='time to seizure onset') ax.set_xticklabels(labels) # -----------------------------------------------------------------# # quantity labels # -----------------------------------------------------------------# quantity_label = r'$\hat{m}$' ymax = 1e-3 ymin = 1 ax.set_ylim(ymin, ymax) T = [1, 1e-1, 1e-2, 1e-3] ax.set_yticks(T) ax.set_yticklabels([0, 0.9, 0.99, 0.999]) # -----------------------------------------------------------------# # Annotate p-values of pairwise comparison # -----------------------------------------------------------------# ax.set_ylim(ymin, ymax) plt.ylabel(quantity_label) plt.xlabel('brain area') sns.despine() # -----------------------------------------------------------------# # Saving # -----------------------------------------------------------------# plotfile_mx = '{}Fig2_2parts_regions_Engel1a.pdf'.format(result_path) plt.savefig(plotfile_mx, bbox_inches='tight') def boxplot_mx_regions(result_path, x, ksteps, fit_method, Engel=False): # -----------------------------------------------------------------# # Read and filter data # -----------------------------------------------------------------# mx_results = '{}mt_{}_results_allrecords_kmin{}_kmax{}_{}.pkl'.format( result_path, x, ksteps[0], ksteps[1], fit_method) mx_allresults = pd.read_pickle(mx_results) mx_allresults = mr_utilities.exclude_inconsistencies(mx_allresults) x = mx_allresults['partno'].max() mx_allresults['epsilon'] = mx_allresults.apply(lambda row: 1-row['m'], axis=1) measure = 'epsilon' def give_region(row): return row['binning'][1:] mx_allresults['region'] = mx_allresults.apply(lambda row: give_region(row), axis=1) if Engel: patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19] mx_allresults = mx_allresults[ (mx_allresults['patient'].isin(patients1a))] # -----------------------------------------------------------------# # write estimates to lists # -----------------------------------------------------------------# df = None regions = ['H', 'A', 'PHC', 'EC'] labels = [] pvals = dict() for reg in regions: mx_region = mx_allresults[mx_allresults['region'] == reg] recordings = mx_region.recording.unique() focal_diff_list = [] contralateral_diff_list = [] all_first_focal = [] all_last_focal = [] all_first_contra = [] all_last_contra = [] recID = [] for recording in recordings: mx_recording = mx_region[mx_region['recording'] == recording] mx_focal = mx_recording[mx_recording['binning_focus'] == 'focal'] mx_contralateral = mx_recording[ mx_recording['binning_focus'] == 'contra-lateral'] if not mx_focal[mx_focal['partno'] == 1].empty and not mx_focal[ mx_focal['partno'] == x].empty: m_focal_first = mx_focal[mx_focal['partno'] == 1][ measure].item() m_focal_last = mx_focal[mx_focal['partno'] == x][measure].item() focal_diff_log = np.log(m_focal_last) - np.log(m_focal_first) all_first_focal.append(m_focal_first) all_last_focal.append(m_focal_last) focal_diff_list.append(focal_diff_log) recID.append(recording) if not mx_contralateral[mx_contralateral['partno'] == 1].empty \ and not mx_contralateral[mx_contralateral['partno'] == x].empty: m_contralateral_first = mx_contralateral[mx_contralateral['partno'] == 1][measure].item() m_contralateral_last = mx_contralateral[mx_contralateral['partno'] == x][measure].item() contralateral_diff_log = np.log(m_contralateral_last) - \ np.log(m_contralateral_first) all_first_contra.append(m_contralateral_first) all_last_contra.append(m_contralateral_last) contralateral_diff_list.append(contralateral_diff_log) # -----------------------------------------------------------------# # Wilcoxon paired test # -----------------------------------------------------------------# stat, p_contra_w = stats.wilcoxon(contralateral_diff_list) stat, p_focal_w = stats.wilcoxon(focal_diff_list) pvals[reg] = p_focal_w # -----------------------------------------------------------------# # Create dataframe for plotting # -----------------------------------------------------------------# quantity = measure if df is None: data1 = { 'region': reg, quantity: all_first_focal, 'recording': recID, 'rec_type': ['10 - 5 min' for _ in range(len(all_first_focal))], 'binning_focus': ['SOZ' for _ in range(len(all_first_focal))] } df = pd.DataFrame(data1) else: data1 = pd.DataFrame({ 'region': reg, quantity: all_first_focal, 'recording': recID, 'rec_type': ['10 - 5 min' for _ in range(len(all_last_focal))], 'binning_focus': ['SOZ' for _ in range(len(all_last_focal))] }) df = df.append(data1, ignore_index=True) data2 = pd.DataFrame({ 'region': reg, quantity: all_last_focal, 'recording': recID, 'rec_type': ['5 - 0 min' for _ in range(len(all_last_focal))], 'binning_focus': ['SOZ' for _ in range(len(all_last_focal))] }) df = df.append(data2, ignore_index=True) data3 = pd.DataFrame({ 'region': reg, quantity: all_first_contra, 'recording': ['-' for _ in range(len(all_first_contra))], 'rec_type': ['10 - 5 min' for _ in range(len(all_first_contra))], 'binning_focus': ['nSOZ' for _ in range(len(all_first_contra))] }) df = df.append(data3, ignore_index=True) data4 = pd.DataFrame({ 'region': reg, quantity: all_last_contra, 'recording': ['-' for _ in range(len(all_first_contra))], 'rec_type': ['5 - 0 min' for _ in range(len(all_last_contra))], 'binning_focus': ['nSOZ' for _ in range(len(all_last_contra))] }) df = df.append(data4, ignore_index=True) n_f_first = len(all_first_focal) n_f_last = len(all_last_focal) labels.append('{}\n({}, {})'.format(reg, n_f_first, n_f_last)) # -----------------------------------------------------------------# # Plotting # -----------------------------------------------------------------# fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(4*1.2))) df = df[df['binning_focus'] == 'SOZ'] g = sns.boxplot(ax=ax, x='region', y=quantity, hue='rec_type', data=df, hue_order=['10 - 5 min', '5 - 0 min'], order=regions, dodge=True, color='slategray', palette=[ 'lightgray','slategray'], linewidth=1.1, width=0.8) g.set_yscale("log") ax.grid(axis='x', b=False) plt.rcParams['legend.title_fontsize'] = 9 plt.legend(loc='center left', bbox_to_anchor=(1, 0.5), fontsize='x-small', title='time to seizure onset') ax.set_xticklabels(labels) # -----------------------------------------------------------------# # quantity labels # -----------------------------------------------------------------# quantity_label = r'$\hat{m}$' ymax = 1e-3 ymin = 1 ax.set_ylim(ymin, ymax) T = [1, 1e-1, 1e-2, 1e-3] ax.set_yticks(T) ax.set_yticklabels([0, 0.9, 0.99, 0.999]) # -----------------------------------------------------------------# # Annotate p-values of pairwise comparison # -----------------------------------------------------------------# y_max = 1.3 * ymax y_min = ymin for ir, region in enumerate(regions): if not np.isnan(pvals[region]): ax.annotate("", xy=(-0.2 + ir * 1.0, 2.0 * y_max), xycoords='data', xytext=(0.2 + ir * 1.0, 2.0 * y_max), textcoords='data', arrowprops=dict(arrowstyle="-", ec='black', connectionstyle="bar,fraction=0.3")) ax.text(-0.35 + ir * 1.0, 0.9*y_max * 1. - abs(y_max - y_min) * 0.0, 'p' + '={:.3f}'.format(pvals[region]), fontsize=8.5) ax.set_ylim(ymin, ymax) plt.ylabel(quantity_label) plt.xlabel('brain area') sns.despine() # -----------------------------------------------------------------# # Saving # -----------------------------------------------------------------# plotfile_mx = '{}Fig2_2parts_regions_Engel{}.pdf'.format(result_path, Engel) plt.savefig(plotfile_mx, bbox_inches='tight') def binning_label(patient, soz): focifile = '../Data/patients.txt' f = open(focifile, 'r') foci = [line.rstrip('\n').split('\t') for line in f.readlines()] focus = 'NA' for i, idf in enumerate(foci[1:]): if int(idf[0]) == int(patient): focus = idf[1] if soz == "SOZ": if focus == "L": binning = "left" elif focus == "R": binning = "right" else: raise Exception("No clear focus found.") elif soz == "nSOZ": if focus == "L": binning = "right" elif focus == "R": binning = "left" else: raise Exception("No clear focus found.") return binning def get_mt(patient, recording, soz, save_path, fit_method, kmin, kmax, windowsize, windowstep, dt): binning = binning_label(patient, soz) mt_file = '{}{}/mt_results_{}_{}_kmin{}_kmax{}_winsize{}_winstep{}_dt{}.pkl'.format( save_path, recording, binning, fit_method, kmin, kmax, windowsize, windowstep, dt) if not os.path.isfile(mt_file): raise Exception('{} not existing'.format(mt_file)) # -----------------------------------------------------------------# # read m(t) data # -----------------------------------------------------------------# mt_frame = pd.read_pickle(mt_file) original_length = len(mt_frame.m.tolist()) rejected_inconsistencies = ['H_nonzero', 'H_linear', 'not_converged'] for incon in rejected_inconsistencies: if not mt_frame.empty: mt_frame = mt_frame[[incon not in mt_frame[ 'inconsistencies'].tolist()[i] for i in range( len(mt_frame['inconsistencies'].tolist()))]] if len(mt_frame.m.tolist()) < 0.95 * original_length or len( mt_frame.m.tolist()) < 200: # in general, the recording should be # 260 windows long. if it is much shorter, it is unclear where to # divide it. raise Exception('Too many inconsistencies') mt_list = mt_frame.m.tolist() return mt_list def mad(data): return np.median(np.abs(data-np.median(data))) def variance_mt(save_path, data_path, fit_method, kmin, kmax, windowsize, windowstep, dt, logepsilon=True, Engel=False): patients = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20] if Engel: patients = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19] laterality = ["SOZ", "nSOZ"] # -----------------------------------------------------------------# # Extract m(t) values an compute variances # -----------------------------------------------------------------# df = None for j, soz in enumerate(laterality): variance_first = [] variance_last = [] diff_list = [] for ir, patient in enumerate(patients): recordings = mr_utilities.list_preictrecordings(patient)[0] for irec, rec in enumerate(recordings): try: m_list = get_mt(patient, rec, soz, data_path, fit_method, kmin, kmax, windowsize, windowstep, dt) m_first = m_list[:int(len(m_list)/2)] m_last = m_list[-int(len(m_list) / 2):] if logepsilon: m_first = [np.log(1 - m) for m in m_first] m_last = [np.log(1 - m) for m in m_last] variance_first.append(mad(m_first)) variance_last.append(mad(m_last)) diff_list.append(mad(m_first)-mad(m_last)) except Exception: # If too many window estimates were excluded in the # function get_mt(), an exception will be risen. # We will simply continue with the remaining recordings. continue # -----------------------------------------------------------------# # create boxplot of variance # -----------------------------------------------------------------# stat, p = stats.wilcoxon(diff_list) print('p-value:'+str(p)) quantity = 'variance' data1 = pd.DataFrame({ quantity: variance_first, 'rec_type': ['10 - 5 min' for _ in range(len(variance_first))], 'binning_focus': [soz for _ in range(len(variance_first))], 'pval': [p for _ in range(len(variance_last))] }) if df is None: df = data1 else: df = df.append(data1, ignore_index=True) data2 = pd.DataFrame({ quantity: variance_last, 'rec_type': ['5 - 0 min' for _ in range(len(variance_last))], 'binning_focus': [soz for _ in range(len(variance_last))], 'pval': [p for _ in range(len(variance_last))] }) df = df.append(data2, ignore_index=True) # -----------------------------------------------------------------# # Plotting # -----------------------------------------------------------------# fig, ax = plt.subplots(1, 1, figsize=(cm2inch(5.5), cm2inch(4.5))) g = sns.boxplot(ax=ax, x='binning_focus', y='variance', hue='rec_type', data=df, hue_order=['10 - 5 min', '5 - 0 min'], order=[ 'SOZ', 'nSOZ'], dodge=True, color='slategray', palette=[ 'lightgray', 'slategray'], linewidth=1.1, width=0.6) g.set_yscale("log") plt.grid(axis='x', b=False) plt.rcParams['legend.title_fontsize'] = 9 plt.legend(loc="upper left", ncol=1, fontsize='x-small', bbox_to_anchor=(0, 1.7), title='time to seizure onset') n_f_pre = len(df[(df['binning_focus'] == 'SOZ') & (df['rec_type'] == '5 - 0 min')].index) n_f_inter = len(df[(df['binning_focus'] == 'SOZ') & (df['rec_type'] == '10 - 5 min')].index) n_c_pre = len(df[(df['binning_focus'] == 'nSOZ') & (df['rec_type'] == '5 - 0 min')].index) n_c_inter = len(df[(df['binning_focus'] == 'nSOZ') & (df['rec_type'] == '10 - 5 min')].index) ax.set_xticklabels(['ipsi\n({}, {})'.format(n_f_pre, n_f_inter), 'contra\n({}, {})'.format(n_c_pre, n_c_inter)]) # -----------------------------------------------------------------# # annotate p-values # -----------------------------------------------------------------# y_max = 2.4 * np.max(variance_first) y_min = 0.4 * np.min(variance_first) ax.annotate("", xy=(-0.2, y_max), xycoords='data', xytext=(0.2, y_max), textcoords='data', arrowprops=dict(arrowstyle="-", ec='black', connectionstyle="bar,fraction=0.2")) ax.text(-0.2, y_max + abs(y_max - y_min) * 0.1, 'p={:.3f}'.format(df[df['binning_focus'] == 'SOZ'].pval.tolist()[ 0]), fontsize=8) ax.annotate("", xy=(0.8, y_max), xycoords='data', xytext=(1.2, y_max), textcoords='data', arrowprops=dict(arrowstyle="-", ec='black', connectionstyle="bar,fraction=0.2")) ax.text(0.8, y_max + abs(y_max - y_min) * 0.1, 'p={:.3f}'.format(df[df['binning_focus'] == 'nSOZ'].pval.tolist( )[0]), fontsize=8) plt.xlabel('') sns.despine() plt.ylabel(r'Variability in $\hat{m}$') plt.yscale('log') plt.savefig('{}Fig2_variance_mt.pdf'.format(save_path), bbox_inches='tight') plt.show() def timeresolved_example(): result_path = '../Results/preictal/singleunits/' ksteps = (1, 400) kmin = ksteps[0] kmax = ksteps[1] windowsize = 20000 windowstep = 500 dt = 4 patients = [13] for patient in patients: mt_examples(patient, result_path, 'offset', kmin, kmax, windowsize, windowstep, dt) def boxplot_2parts_hemispheres(): x = 2 ksteps = (1, 400) fit_method = 'offset' result_path = '../Results/preictal/singleunits/' boxplot_mx(result_path, x, ksteps, fit_method, Engel=False) boxplot_mx(result_path, x, ksteps, fit_method, Engel=True) def boxplot_2parts_regions(): x = 2 ksteps = (1, 400) fit_method = 'offset' result_path = '../Results/preictal/singleunits_regions/' boxplot_mx_regions(result_path, x, ksteps, fit_method, Engel=False) pointplot_mx_regions_Engel(result_path, x, ksteps, fit_method) def visualize_mt_variance(): save_path = '../Results/preictal/singleunits/' data_path = '../Results/preictal/singleunits/' variance_mt(save_path, data_path, fit_method='offset', kmin=1, kmax=400, windowsize=20000, windowstep=500, dt=4, Engel=False, logepsilon=True) def main(): timeresolved_example() boxplot_2parts_hemispheres() boxplot_2parts_regions() visualize_mt_variance() if __name__ == '__main__': main()