123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731 |
- import mr_utilities
- import pandas as pd
- import numpy as np
- import matplotlib.pylab as plt
- import seaborn as sns
- from scipy import stats
- import hierarchical_model
- colors = ['steelblue', 'slategray', 'firebrick', 'darkblue', 'darkgreen']
- def cm2inch(value):
- return value/2.54
- def boxplot_focal_contra_onlypaired(pkl_file_pre, pkl_file_inter, resultpath):
- mr_results_pre = pd.read_pickle(pkl_file_pre)
- mr_results_inter = pd.read_pickle(pkl_file_inter)
- mr_results = pd.concat([mr_results_pre, mr_results_inter],
- ignore_index=True)
- mr_results = mr_utilities.exclude_inconsistencies(mr_results)
- mr_results['logepsilon'] = mr_results.apply(lambda row: np.log(1-row['m']), axis=1)
- mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
- quantity = 'epsilon'
- # ================================================================== #
- # compute pvalue (paired, just uses data with pair)
- # ================================================================== #
- paired_recordings = []
- for rec in mr_results.recording.unique():
- recframe = mr_results[mr_results['recording'] == rec]
- if not len(recframe['m'].tolist()) == 2:
- continue
- paired_recordings.append(rec)
- mr_results_paired = mr_results[mr_results['recording'].isin(paired_recordings)]
- p_hierarchical = hierarchical_model.hierarchical_model_hemispheres(
- mr_results_paired)
- # ================================================================== #
- # Plotting
- # ================================================================== #
- fig, ax = plt.subplots(1, 1, figsize=(cm2inch(4), cm2inch(3)))
- gb = sns.boxplot(ax=ax, x="binning_focus", y=quantity,
- data=mr_results_paired, order=['focal', 'contra-lateral'],
- dodge=True, palette=['darkblue', 'slategray'])
- for patch in ax.artists:
- r, g, b, a = patch.get_facecolor()
- patch.set_facecolor((r, g, b, .4))
- patch.set_edgecolor('gray')
- sns.stripplot(ax=ax, x="binning_focus", y=quantity,
- data=mr_results_paired, alpha=0.9, order=['focal',
- 'contra-lateral'], size=2, dodge=True, palette=[
- 'darkblue', 'slategray'], jitter=False)
- gb.set_yscale("log")
- plt.grid(axis='x', b=False)
- n_SOZ = len(mr_results_paired[mr_results_paired['binning_focus']
- == 'focal']['m'].tolist())
- n_nSOZ = len(mr_results_paired[mr_results_paired['binning_focus'] ==
- 'contra-lateral']['m'].tolist())
- ax.set_xticklabels(['ipsi\n({})'.format(n_SOZ), 'contra\n({})'.format(n_nSOZ)])
- ax.set_xlabel('')
- quantity_label = r'$\hat{m}$'
- ymin = 1.4 * np.max(mr_results[quantity].tolist())
- ymax = 0.6*np.min(mr_results[quantity].tolist())
- 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])
- plt.ylabel(quantity_label)
- sns.despine()
- # ================================================================== #
- # Annotate pvalues
- # ================================================================== #
- x1, x2 = 0, 1
- y, h = ymax * (1.1), -0.1 * ymax
- ax.plot([x1, x1, x2, x2], [y, y + h, y + h, y], lw=1.5, c='black')
- ax.text((x1 + x2) * .5, y + 4*h, r"$p$"+"={:.3f}".format(
- p_hierarchical), ha='center', va='bottom', color='black', fontsize=10)
- # ================================================================== #
- # Final figure adjustments + saving
- # ================================================================== #
- sns.despine(fig)
- fig.subplots_adjust(bottom=None, right=None,
- wspace=0.8, hspace=0.5, left=0.15, top=1.2)
- resultfile = '{}Fig1_boxplot_onlypaired.pdf'.format(resultpath)
- plt.savefig(resultfile, bbox_inches='tight')
- plt.show()
- plt.close()
- def boxplot_focal_contra_Engel(pkl_file_pre, pkl_file_inter, resultpath):
- mr_results_pre = pd.read_pickle(pkl_file_pre)
- mr_results_inter = pd.read_pickle(pkl_file_inter)
- mr_results = pd.concat([mr_results_pre, mr_results_inter],
- ignore_index=True)
- patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
- mr_results = mr_results[(mr_results['patient'].isin(patients1a))]
- mr_results = mr_utilities.exclude_inconsistencies(mr_results)
- mr_results['logepsilon'] = mr_results.apply(lambda row:
- np.log(1-row['m']), axis=1)
- mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
- quantity = 'epsilon'
- # ================================================================== #
- # compute pvalue (unpaired, all data)
- # ================================================================== #
- p_hierarchical = hierarchical_model.hierarchical_model_hemispheres(
- mr_results)
- # ================================================================== #
- # Plotting
- # ================================================================== #
- fig, ax = plt.subplots(1, 1, figsize=(cm2inch(4), cm2inch(3)))
- gb = sns.boxplot(ax=ax, x="binning_focus", y=quantity,
- data=mr_results, order=['focal', 'contra-lateral'],
- dodge=True, palette=['darkblue', 'slategray'])
- for patch in ax.artists:
- r, g, b, a = patch.get_facecolor()
- patch.set_facecolor((r, g, b, .4))
- patch.set_edgecolor('gray')
- sns.stripplot(ax=ax, x="binning_focus", y=quantity,
- data=mr_results, alpha=0.9, order=['focal',
- 'contra-lateral'], size=2, dodge=True, palette=[
- 'darkblue', 'slategray'], jitter=False)
- gb.set_yscale("log")
- plt.grid(axis='x', b=False)
- n_SOZ = len(mr_results[mr_results['binning_focus']=='focal']['m'].tolist())
- n_nSOZ = len(mr_results[mr_results['binning_focus'] == 'contra-lateral'][
- 'm'].tolist())
- ax.set_xticklabels(['ipsi\n({})'.format(n_SOZ), 'contra\n({})'.format(
- n_nSOZ)])
- ax.set_xlabel('')
- quantity_label = r'$\hat{m}$'
- ymin = 1.4 * np.max(mr_results[quantity].tolist())
- ymax = 0.6*np.min(mr_results[quantity].tolist())
- 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])
- plt.ylabel(quantity_label)
- sns.despine()
- # ================================================================== #
- # Annotate pvalues
- # ================================================================== #
- x1, x2 = 0, 1
- y, h = ymax * (1.1), -0.1 * ymax
- ax.plot([x1, x1, x2, x2], [y, y + h, y + h, y], lw=1.5, c='black')
- ax.text((x1 + x2) * .5, y + 4*h, r"$p$"+"={:.3f}".format(p_hierarchical),
- ha='center', va='bottom', color='black', fontsize=10)
- # ================================================================== #
- # Final figure adjustments + saving
- # ================================================================== #
- sns.despine(fig)
- fig.subplots_adjust(bottom=None, right=None,
- wspace=0.8, hspace=0.5, left=0.15, top=1.2)
- resultfile = '{}Fig1_boxplot_1aEngel.pdf'.format(resultpath)
- plt.savefig(resultfile, bbox_inches='tight')
- plt.show()
- plt.close()
- def boxplot_focal_contra(pkl_file_pre, pkl_file_inter, resultpath):
- mr_results_pre = pd.read_pickle(pkl_file_pre)
- mr_results_inter = pd.read_pickle(pkl_file_inter)
- mr_results = pd.concat([mr_results_pre, mr_results_inter],
- ignore_index=True)
- mr_results = mr_utilities.exclude_inconsistencies(mr_results)
- mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
- mr_results['logepsilon'] = mr_results.apply(lambda row: np.log(1-row['m']),
- axis=1)
- quantity = 'epsilon'
- # ================================================================== #
- # compute p-value (unpaired, all data)
- # ================================================================== #
- p_hierarchical = hierarchical_model.hierarchical_model_hemispheres(mr_results)
- # ================================================================== #
- # Plotting
- # ================================================================== #
- fig, ax = plt.subplots(1, 1, figsize=(cm2inch(4), cm2inch(3)))
- gb = sns.boxplot(ax=ax, x="binning_focus", y=quantity,
- data=mr_results, order=['focal', 'contra-lateral'],
- dodge=True, palette=['darkblue', 'slategray'])
- for patch in ax.artists:
- r, g, b, a = patch.get_facecolor()
- patch.set_facecolor((r, g, b, .4))
- patch.set_edgecolor('gray')
- sns.stripplot(ax=ax, x="binning_focus", y=quantity,
- data=mr_results, alpha=0.9, order=['focal','contra-lateral'],
- size=2, dodge=True, palette=['darkblue', 'slategray'],
- jitter=False)
- gb.set_yscale("log")
- plt.grid(axis='x', b=False)
- n_SOZ = len(mr_results[mr_results['binning_focus'] == 'focal'][
- 'm'].tolist())
- n_nSOZ = len(mr_results[mr_results['binning_focus'] == 'contra-lateral'][
- 'm'].tolist())
- ax.set_xticklabels(['ipsi\n({})'.format(n_SOZ), 'contra\n({})'.format(
- n_nSOZ)])
- ax.set_xlabel('')
- quantity_label = r'$\hat{m}$'
- ymin = 1.4 * np.max(mr_results[quantity].tolist())
- ymax = 0.6*np.min(mr_results[quantity].tolist())
- 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])
- plt.ylabel(quantity_label)
- sns.despine()
- # ================================================================== #
- # Annotate pvalues
- # ================================================================== #
- x1, x2 = 0, 1
- y, h = ymax * (1.1), -0.1 * ymax
- ax.plot([x1, x1, x2, x2], [y, y + h, y + h, y], lw=1.5, c='black')
- ax.text((x1 + x2) * .5, y + 4*h, r"$p$"+"={:.3f}".format(p_hierarchical),
- ha='center', va='bottom', color='black', fontsize=10)
- # ================================================================== #
- # Final figure adjustments + saving
- # ================================================================== #
- sns.despine(fig)
- fig.subplots_adjust(bottom=None, right=None,
- wspace=0.8, hspace=0.5, left=0.15, top=1.2)
- resultfile = '{}Fig1_boxplot.pdf'.format(resultpath)
- plt.savefig(resultfile, bbox_inches='tight')
- plt.show()
- plt.close()
- def patientwise_pointplots(pkl_file_pre, pkl_file_inter, result_path):
- # -----------------------------------------------------------------#
- # Read data
- # -----------------------------------------------------------------#
- mr_results_inter = pd.read_pickle(pkl_file_inter)
- mr_results_inter['rec_type'] = mr_results_inter.apply(lambda row: 'inter',
- axis=1)
- mr_results_pre = pd.read_pickle(pkl_file_pre)
- mr_results_pre['rec_type'] = mr_results_pre.apply(lambda row: 'pre',
- axis=1)
- mr_results = pd.concat([mr_results_pre, mr_results_inter],
- ignore_index=True)
- mr_results = mr_utilities.exclude_inconsistencies(mr_results)
- mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
- quantity = 'epsilon'
- # -----------------------------------------------------------------#
- # Run through all patients that also have preictal data
- # get their interictal and preictal results plot them as a pointplot
- # -----------------------------------------------------------------#
- patients = mr_results['patient'].unique()
- patients_sorted = np.sort(patients)
- fig, ax = plt.subplots(5, 4, figsize=(12, 1.5 * 5))
- v = [0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3]
- h = [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4]
- for ip, patient in enumerate(patients_sorted):
- mr_patient = mr_results[mr_results['patient'] == patient]
- recordings = mr_patient.recording.unique()
- all_SOZ = mr_patient[mr_patient['binning_focus'] == 'focal'][
- quantity].tolist()
- all_nSOZ = mr_patient[mr_patient['binning_focus'] == 'contra-lateral'][
- quantity].tolist()
- stat, p_unpaired = stats.mannwhitneyu(all_SOZ, all_nSOZ,
- alternative='two-sided')
- for rec in recordings:
- mr_rec = mr_patient[mr_patient['recording'] == rec]
- if len(mr_rec[quantity].tolist()) == 2:
- q_SOZ = mr_rec[mr_rec['binning_focus'] == 'focal'][
- quantity].item()
- q_nSOZ = mr_rec[mr_rec['binning_focus'] == 'contra-lateral'][
- quantity].item()
- color = 'darkorange'
- if mr_rec['rec_type'].tolist()[0] == 'pre':
- color = 'indianred'
- ax[h[ip], v[ip]].plot([1, 2], [q_SOZ, q_nSOZ], c=color,
- marker='.', markersize=6)
- elif len(mr_rec[mr_rec['binning_focus'] == 'focal'].index) == 1:
- q_SOZ = mr_rec[mr_rec['binning_focus'] == 'focal'][
- quantity].item()
- color = 'darkorange'
- if mr_rec['rec_type'].tolist()[0] == 'pre':
- color = 'indianred'
- ax[h[ip], v[ip]].plot([1], [q_SOZ], c=color, marker='.',
- markersize=6)
- elif len(mr_rec[mr_rec['binning_focus'] ==
- 'contra-lateral'].index) == 1:
- q_nSOZ = mr_rec[mr_rec['binning_focus'] == 'contra-lateral'][
- quantity].item()
- color = 'darkorange'
- if mr_rec['rec_type'].tolist()[0] == 'pre':
- color = 'indianred'
- ax[h[ip], v[ip]].plot([2], [q_nSOZ], c=color, marker='.',
- markersize=6)
- # if the number of recordings is sufficiently large, output p-value
- if len(all_SOZ)+len(all_nSOZ) >= 8:
- ax[h[ip], v[ip]].text(0.5, 0.9, 'p={:.3f}'.format(p_unpaired),
- horizontalalignment='center', verticalalignment='center',
- transform=ax[h[ip], v[ip]].transAxes, color='black',
- fontsize=10)
- ax[h[ip], v[ip]].set_title('Patient {}'.format(patient))
- ax[h[ip], v[ip]].set_yscale("log")
- quantity_label = r'$\hat{m}$'
- ax[h[ip], v[ip]].set_ylabel(quantity_label)
- ymin = 1.1
- ymax = 1e-3
- ax[h[ip], v[ip]].set_ylim(ymin, ymax)
- ax[h[ip], v[ip]].set_xlim(0.7, 2.3)
- T = [1, 2]
- ax[h[ip], v[ip]].set_xticks(T)
- ax[h[ip], v[ip]].set_xticklabels(['ipsi', 'contra'])
- T = [1, 1e-1, 1e-2, 1e-3]
- ax[h[ip], v[ip]].set_yticks(T)
- ax[h[ip], v[ip]].set_yticklabels([0, 0.9, 0.99, 0.999])
- ax[h[ip], v[ip]].set_xlabel(' ')
- ax[h[ip], v[ip]].set_xticklabels(['ipsi', 'contra'])
- sns.despine(fig)
- fig.subplots_adjust(bottom=None, right=None,
- wspace=0.7, hspace=0.7, left=0.15, top=1.2)
- plotfile_mx = '{}S2_patients.pdf'.format(result_path)
- fig.savefig(plotfile_mx, bbox_inches='tight')
- plt.show()
- def boxplot_brainregions_focus_onlypaired(pkl_file_pre, pkl_file_inter,
- result_path):
- regions = ['H', 'A', 'PHC', 'EC']
- mr_results_pre = pd.read_pickle(pkl_file_pre)
- mr_results_inter = pd.read_pickle(pkl_file_inter)
- mr_results = pd.concat([mr_results_pre, mr_results_inter],
- ignore_index=True)
- mr_results = mr_utilities.exclude_inconsistencies(mr_results)
- # -----------------------------------------------------------------#
- # Write the brainregion, independent of hemisphere, into frame
- # -----------------------------------------------------------------#
- def give_region(row):
- return row['binning'][1:]
- # -----------------------------------------------------------------#
- # Write the brainregion and focus into frame
- # -----------------------------------------------------------------#
- def give_region_and_focus(row):
- regionname = row['binning'][1:]
- focus = row['binning_focus'][0] # f for focal and c for contralateral
- return '{}{}'.format(focus, regionname)
- mr_results['region'] = mr_results.apply(lambda row: give_region(row),
- axis=1)
- mr_results['regionfocus'] = mr_results.apply(lambda row:
- give_region_and_focus(row), axis=1)
- mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
- mr_results['logepsilon'] = mr_results.apply(
- lambda row: np.log(1 - row['m']), axis=1)
- quantity = 'epsilon'
- # -----------------------------------------------------------------#
- # Filter to recordings and regions for which paired data exists
- # -----------------------------------------------------------------#
- mr_results_paired = pd.DataFrame(columns=mr_results.columns)
- for rec in mr_results.recording.unique():
- recframe = mr_results[mr_results['recording'] == rec]
- for area in recframe['region'].unique():
- area_frame = recframe[recframe['region'] == area]
- if not len(area_frame['m'].tolist()) == 2:
- continue
- mr_results_paired = mr_results_paired.append(area_frame,
- ignore_index=True)
- mr_results = mr_results_paired
- # -----------------------------------------------------------------#
- # Compute p-values
- # -----------------------------------------------------------------#
- global_pvals, pvals = hierarchical_model.hierarchical_model_regions(
- mr_results)
- # -----------------------------------------------------------------#
- # Boxplot of regions
- # -----------------------------------------------------------------#
- fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(3*1.2)))
- gb = sns.boxplot(ax=ax, x="region", y=quantity, hue='binning_focus',
- data=mr_results, order=regions, hue_order=['focal',
- 'contra-lateral'], dodge=True, palette=['darkblue',
- 'slategray'], linewidth=1)
- for patch in ax.artists:
- r, g, b, a = patch.get_facecolor()
- patch.set_facecolor((r, g, b, .4))
- patch.set_edgecolor('gray')
- sns.stripplot(ax=ax, x="region", y=quantity, hue='binning_focus',
- data=mr_results, alpha=0.9, order=regions,
- hue_order=['focal', 'contra-lateral'], jitter=False,
- size=2, palette=['darkblue', 'slategray'], dodge=True)
- ax.get_legend().set_visible(False)
- plt.grid(axis='x', b=False)
- gb.set_yscale("log")
- labels = []
- for ireg, reg in enumerate(regions):
- mr_results_r = mr_results[mr_results['region'] == reg]
- n_reg_f = len(mr_results_r[mr_results_r['binning_focus'] == 'focal'][
- 'm'].tolist())
- n_reg_c = len(mr_results_r[mr_results_r['binning_focus'] ==
- 'contra-lateral']['m'].tolist())
- labels.append('{}\n({}, {})'.format(reg, n_reg_f, n_reg_c))
- ax.set_xticklabels(labels)
- quantity_label = r'$\hat{m}$'
- ymin = 1.4 * np.max(mr_results[quantity].tolist())
- ymax = 0.4 * np.min(mr_results[quantity].tolist())
- 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])
- plt.ylabel(quantity_label)
- sns.despine()
- # -----------------------------------------------------------------#
- # Annotate p-values of pairwise comparison
- # -----------------------------------------------------------------#
- y_max = 1.3*ymax
- y_min = ymin
- for ir, region in enumerate(regions):
- ax.annotate("", xy=(-0.2+ir*1.0, 1.5*y_max), xycoords='data',
- xytext=(0.2+ir*1.0, 1.5*y_max), textcoords='data',
- arrowprops=dict(arrowstyle="-", ec='black',
- connectionstyle="bar,fraction=0.3"))
- ax.text(-0.35+ir*1.0, y_max * 0.9 - abs(y_max - y_min) * 0.0,
- 'p'+'={:.3f}'.format(pvals[region]), fontsize=8.5)
- # -----------------------------------------------------------------#
- # Adjust some plot settings and save figures
- # -----------------------------------------------------------------#
- ax.set_xlabel('brain area')
- ax.set_ylim((ymin, ymax))
- ax.set_title('overall model: ' + str(global_pvals), fontsize=4)
- fig.subplots_adjust(bottom=None, right=None,
- wspace=0.2, hspace=0.5, left=0.15, top=1.2)
- plotfile_mx = '{}Fig1_regions_onlypaired.pdf'.format(result_path)
- fig.savefig(plotfile_mx, bbox_inches='tight')
- plt.show()
- plt.close('all')
- def boxplot_brainregions_focus(pkl_file_pre, pkl_file_inter, result_path):
- regions = ['H', 'A', 'PHC', 'EC']
- mr_results_pre = pd.read_pickle(pkl_file_pre)
- mr_results_inter = pd.read_pickle(pkl_file_inter)
- mr_results = pd.concat([mr_results_pre, mr_results_inter],
- ignore_index=True)
- mr_results = mr_utilities.exclude_inconsistencies(mr_results)
- # -----------------------------------------------------------------#
- # Write the brainregion, independent of hemisphere, into frame
- # -----------------------------------------------------------------#
- def give_region(row):
- return row['binning'][1:]
- # -----------------------------------------------------------------#
- # Write the brainregion and focus into frame
- # -----------------------------------------------------------------#
- def give_region_and_focus(row):
- regionname = row['binning'][1:]
- focus = row['binning_focus'][0] # f for focal and c for contralateral
- return '{}{}'.format(focus, regionname)
- mr_results['region'] = mr_results.apply(lambda row: give_region(row),
- axis=1)
- mr_results['regionfocus'] = mr_results.apply(lambda row:
- give_region_and_focus(row), axis=1)
- mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
- mr_results['logepsilon'] = mr_results.apply(
- lambda row: np.log(1 - row['m']), axis=1)
- quantity = 'epsilon'
- # -----------------------------------------------------------------#
- # Compute p-values
- # -----------------------------------------------------------------#
- global_pvals, pvals = hierarchical_model.hierarchical_model_regions(
- mr_results)
- # -----------------------------------------------------------------#
- # Boxplot of regions
- # -----------------------------------------------------------------#
- fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(3*1.2)))
- gb = sns.boxplot(ax=ax, x="region", y=quantity, hue='binning_focus',
- data=mr_results, order=regions, hue_order=['focal',
- 'contra-lateral'], dodge=True, palette=['darkblue',
- 'slategray'], linewidth=1)
- gb.set_yscale("log")
- for patch in ax.artists:
- r, g, b, a = patch.get_facecolor()
- patch.set_facecolor((r, g, b, .4))
- patch.set_edgecolor('gray')
- sns.stripplot(ax=ax, x="region", y=quantity, hue='binning_focus',
- data=mr_results, alpha=0.9, order=regions,
- hue_order=['focal', 'contra-lateral'], jitter=False,
- size=2, palette=['darkblue', 'slategray'], dodge=True)
- ax.get_legend().set_visible(False)
- plt.grid(axis='x', b=False)
- labels = []
- for ireg, reg in enumerate(regions):
- mr_results_r = mr_results[mr_results['region'] == reg]
- n_reg_f = len(mr_results_r[mr_results_r['binning_focus'] == 'focal'][
- 'm'].tolist())
- n_reg_c = len(mr_results_r[mr_results_r['binning_focus']
- == 'contra-lateral']['m'].tolist())
- labels.append('{}\n({}, {})'.format(reg, n_reg_f, n_reg_c))
- ax.set_xticklabels(labels)
- quantity_label = r'$\hat{m}$'
- ymin = 1.4 * np.max(mr_results[quantity].tolist())
- ymax = 0.4 * np.min(mr_results[quantity].tolist())
- 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])
- plt.ylabel(quantity_label)
- sns.despine()
- # -----------------------------------------------------------------#
- # Annotate p-values of pairwise comparison
- # -----------------------------------------------------------------#
- y_max = 1.3*ymax
- y_min = ymin
- for ir, region in enumerate(regions):
- ax.annotate("", xy=(-0.2+ir*1.0, 1.5*y_max), xycoords='data',
- xytext=(0.2+ir*1.0, 1.5*y_max), textcoords='data',
- arrowprops=dict(arrowstyle="-", ec='black',
- connectionstyle="bar,fraction=0.3"))
- ax.text(-0.35+ir*1.0, y_max * 0.9 - abs(y_max - y_min) * 0.0,
- 'p'+'={:.3f}'.format(pvals[region]), fontsize=8.5)
- # -----------------------------------------------------------------#
- # Adjust some plot settings and save figures
- # -----------------------------------------------------------------#
- ax.set_xlabel('brain area')
- ax.set_ylim((ymin, ymax))
- ax.set_title('overall model: ' + str(global_pvals), fontsize=4)
- fig.subplots_adjust(bottom=None, right=None,
- wspace=0.2, hspace=0.5, left=0.15, top=1.2)
- plotfile_mx = '{}Fig1_regions.pdf'.format(result_path)
- fig.savefig(plotfile_mx, bbox_inches='tight')
- plt.show()
- plt.close('all')
- def boxplot_brainregions_focus_Engel(pkl_file_pre,pkl_file_inter,result_path):
- regions = ['H', 'A', 'PHC', 'EC']
- mr_results_pre = pd.read_pickle(pkl_file_pre)
- mr_results_inter = pd.read_pickle(pkl_file_inter)
- mr_results = pd.concat([mr_results_pre, mr_results_inter],
- ignore_index=True)
- patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
- mr_results = mr_results[(mr_results['patient'].isin(patients1a))]
- mr_results = mr_utilities.exclude_inconsistencies(mr_results)
- # -----------------------------------------------------------------#
- # Write the brainregion, independent of hemisphere, into frame
- # -----------------------------------------------------------------#
- def give_region(row):
- return row['binning'][1:]
- # -----------------------------------------------------------------#
- # Write the brainregion and focus into frame
- # -----------------------------------------------------------------#
- def give_region_and_focus(row):
- regionname = row['binning'][1:]
- focus = row['binning_focus'][0] # f for focal and c for contralateral
- return '{}{}'.format(focus, regionname)
- mr_results['region'] = mr_results.apply(lambda row: give_region(row),
- axis=1)
- mr_results['regionfocus'] = mr_results.apply(lambda row:
- give_region_and_focus(row), axis=1)
- mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
- mr_results['logepsilon'] = mr_results.apply(
- lambda row: np.log(1 - row['m']), axis=1)
- quantity = 'epsilon'
- # -----------------------------------------------------------------#
- # Compute p-values
- # -----------------------------------------------------------------#
- global_pvals, pvals = hierarchical_model.hierarchical_model_regions(
- mr_results)
- # -----------------------------------------------------------------#
- # Boxplot of regions
- # -----------------------------------------------------------------#
- fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(3*1.2)))
- gb = sns.boxplot(ax=ax, x="region", y=quantity, hue='binning_focus',
- data=mr_results, order=regions, hue_order=['focal',
- 'contra-lateral'],
- dodge=True, palette=['darkblue', 'slategray'], linewidth=1)
- for patch in ax.artists:
- r, g, b, a = patch.get_facecolor()
- patch.set_facecolor((r, g, b, .4))
- patch.set_edgecolor('gray')
- sns.stripplot(ax=ax, x="region", y=quantity, hue='binning_focus',
- data=mr_results, alpha=0.9, order=regions,
- hue_order=['focal', 'contra-lateral'], jitter=False,
- size=2, palette=['darkblue', 'slategray'], dodge=True)
- ax.get_legend().set_visible(False)
- plt.grid(axis='x', b=False)
- gb.set_yscale("log")
- labels = []
- for ireg, reg in enumerate(regions):
- mr_results_r = mr_results[mr_results['region']==reg]
- n_reg_f = len(mr_results_r[mr_results_r['binning_focus']=='focal'][
- 'm'].tolist())
- n_reg_c = len(mr_results_r[mr_results_r['binning_focus']==
- 'contra-lateral']['m'].tolist())
- labels.append('{}\n({}, {})'.format(reg, n_reg_f, n_reg_c))
- ax.set_xticklabels(labels)
- quantity_label = r'$\hat{m}$'
- ymin = 1.4 * np.max(mr_results[quantity].tolist())
- ymax = 0.4 * np.min(mr_results[quantity].tolist())
- 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])
- plt.ylabel(quantity_label)
- sns.despine()
- # -----------------------------------------------------------------#
- # Annotate p-values of pairwise comparison
- # -----------------------------------------------------------------#
- y_max = 0.85*ymax
- if quantity == 'epsilon':
- y_max = 1.3*ymax
- y_min = ymin
- for ir, region in enumerate(regions):
- ax.annotate("", xy=(-0.2+ir*1.0, 1.5*y_max), xycoords='data',
- xytext=(0.2+ir*1.0, 1.5*y_max), textcoords='data',
- arrowprops=dict(arrowstyle="-", ec='black',
- connectionstyle="bar,fraction=0.3"))
- ax.text(-0.35+ir*1.0, y_max * 0.9 - abs(y_max - y_min) * 0.0,
- 'p'+'={:.3f}'.format(pvals[region]), fontsize=8.5)
- # -----------------------------------------------------------------#
- # Adjust some plot settings and save figures
- # -----------------------------------------------------------------#
- ax.set_xlabel('brain area')
- ax.set_ylim((ymin, ymax))
- ax.set_title('overall model: ' + str(global_pvals), fontsize=4)
- fig.subplots_adjust(bottom=None, right=None,
- wspace=0.2, hspace=0.5, left=0.15, top=1.2)
- plotfile_mx = '{}Fig1_regions_Engel1a.pdf'.format(result_path)
- fig.savefig(plotfile_mx, bbox_inches='tight')
- plt.show()
- plt.close('all')
- def boxplot_hemispheres():
- result_path = '../Results/preictal/singleunits/'
- pkl_file = result_path + 'mr_results.pkl'
- pkl_file_pre = pkl_file
- result_path_inter = '../Results/interictal/singleunits/'
- pkl_file_inter = result_path_inter + 'mr_results.pkl'
- boxplot_focal_contra(pkl_file_pre, pkl_file_inter,
- result_path)
- boxplot_focal_contra_Engel(pkl_file_pre, pkl_file_inter,
- result_path)
- boxplot_focal_contra_onlypaired(pkl_file_pre, pkl_file_inter,
- result_path)
- def pointplot_patients_hemispheres():
- pkl_file_pre = '../Results/preictal/singleunits/mr_results.pkl'
- pkl_file_inter = '../Results/interictal/singleunits/mr_results.pkl'
- result_path_pointplot = '../Results/preictal/singleunits/'
- patientwise_pointplots(pkl_file_pre, pkl_file_inter,
- result_path_pointplot)
- def boxplot_regions():
- result_path = '../Results/preictal/singleunits_regions/'
- pkl_file = result_path + 'mr_results.pkl'
- result_path_inter = '../Results/interictal/singleunits_regions/'
- pkl_file_inter = result_path_inter + 'mr_results.pkl'
- boxplot_brainregions_focus(pkl_file, pkl_file_inter, result_path)
- boxplot_brainregions_focus_Engel(pkl_file, pkl_file_inter, result_path)
- boxplot_brainregions_focus_onlypaired(pkl_file, pkl_file_inter,
- result_path)
- def main():
- boxplot_hemispheres()
- pointplot_patients_hemispheres()
- boxplot_regions()
- if __name__ == '__main__':
- main()
|