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()