gin commit from UX410UAK

Modified files: 7

python_code/analysis_PlosCB/create_figure_1.py

@@ -13,47 +13,25 @@ def cm2inch(value):
     return value/2.54
 
 
-def boxplot_focal_contra_onlypaired(pkl_file_pre, pkl_file_inter,
-                                 resultpath, excludeInconsistencies, **kwargs):
+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_results[(mr_results['fit_method'] == 'offset')]
-    conditions = list()
-
-    for key in kwargs.keys():
-        conditions.append(mr_results[key] == kwargs[key])
-    if len(conditions) > 0:
-        mr_results = mr_results[mr_utilities.conjunction(*conditions)]
-
-    if excludeInconsistencies:
-        rejected_inconsistencies = ['H_nonzero', 'H_linear']
-        for incon in rejected_inconsistencies:
-            mr_results = mr_results[[incon not in mr_results[
-                'inconsistencies'].tolist()[i] for i in range(
-                    len(mr_results['inconsistencies'].tolist()))]]
-
-    mr_results['logtau'] = mr_results.apply(lambda row: np.log(row['tau']),
-                                            axis=1)
-    mr_results['logepsilon'] = mr_results.apply(lambda row: abs(np.log(1-row[
-        'm'])), axis=1)
+    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)
     # ================================================================== #
-    truediff = []
     paired_recordings = []
     for rec in mr_results.recording.unique():
-        recframe = mr_results[mr_results['recording']==rec]
+        recframe = mr_results[mr_results['recording'] == rec]
         if not len(recframe['m'].tolist()) == 2:
             continue
-        focalval = recframe[recframe['binning_focus']=='focal'][quantity].item()
-        contraval = recframe[recframe['binning_focus'] == 'contra-lateral'][
-            quantity].item()
-        truediff.append(focalval - contraval)
         paired_recordings.append(rec)
 
     mr_results_paired = mr_results[mr_results['recording'].isin(paired_recordings)]
@@ -78,14 +56,14 @@ def boxplot_focal_contra_onlypaired(pkl_file_pre, pkl_file_inter,
                   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)])
+    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())
@@ -111,17 +89,13 @@ def boxplot_focal_contra_onlypaired(pkl_file_pre, pkl_file_inter,
     sns.despine(fig)
     fig.subplots_adjust(bottom=None, right=None,
                         wspace=0.8, hspace=0.5, left=0.15, top=1.2)
-    kmin = mr_results["kmin"].unique()[0]
-    kmax = mr_results["kmax"].unique()[0]
-    resultfile = '{}Fig1_boxplot_{}_kmin{}_kmax{}_onlypaired.pdf'.format(
-        resultpath,quantity, kmin, kmax)
+    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, excludeInconsistencies, **kwargs):
+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)
@@ -130,27 +104,11 @@ def boxplot_focal_contra_Engel(pkl_file_pre, pkl_file_inter,
 
     patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
     mr_results = mr_results[(mr_results['patient'].isin(patients1a))]
-    mr_results = mr_results[(mr_results['fit_method'] == 'offset')]
-
-    conditions = list()
-    for key in kwargs.keys():
-        conditions.append(mr_results[key] == kwargs[key])
-    if len(conditions) > 0:
-        mr_results = mr_results[mr_utilities.conjunction(*conditions)]
-
-    if excludeInconsistencies:
-        rejected_inconsistencies = ['H_nonzero', 'H_linear']
-        for incon in rejected_inconsistencies:
-            mr_results = mr_results[[incon not in mr_results[
-                'inconsistencies'].tolist()[i] for i in range(
-                    len(mr_results['inconsistencies'].tolist()))]]
-
-    mr_results['logtau'] = mr_results.apply(lambda row: np.log(row['tau']),
-                                            axis=1)
-    mr_results['logepsilon'] = mr_results.apply(lambda row: abs(np.log(1-row[
-        'm'])), axis=1)
-    mr_results['epsilon'] = mr_results.apply(lambda row: 1-row[
-        'm'], axis=1)
+    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)
@@ -174,8 +132,8 @@ def boxplot_focal_contra_Engel(pkl_file_pre, pkl_file_inter,
 
     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)
+                  '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())
@@ -208,42 +166,22 @@ def boxplot_focal_contra_Engel(pkl_file_pre, pkl_file_inter,
     sns.despine(fig)
     fig.subplots_adjust(bottom=None, right=None,
                         wspace=0.8, hspace=0.5, left=0.15, top=1.2)
-    kmin = mr_results["kmin"].unique()[0]
-    kmax = mr_results["kmax"].unique()[0]
-    resultfile = '{}Fig1_boxplot_{}_kmin{}_kmax{}1aEngel.pdf'.format(
-        resultpath,quantity, kmin, kmax)
+    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,
-                         excludeInconsistencies, **kwargs):
+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_results[(mr_results['fit_method'] == 'offset')]
-    conditions = list()
-    for key in kwargs.keys():
-        conditions.append(mr_results[key] == kwargs[key])
-    if len(conditions) > 0:
-        mr_results = mr_results[mr_utilities.conjunction(*conditions)]
-
-    if excludeInconsistencies:
-        rejected_inconsistencies = ['H_nonzero', 'H_linear']
-        for incon in rejected_inconsistencies:
-            mr_results = mr_results[[incon not in mr_results[
-                'inconsistencies'].tolist()[i] for i in range(
-                    len(mr_results['inconsistencies'].tolist()))]]
-
-    mr_results['logtau'] = mr_results.apply(lambda row: np.log(row['tau']),
-                                            axis=1)
-    mr_results['logepsilon'] = mr_results.apply(lambda row: abs(np.log(1-row[
-                                'm'])), axis=1)
+    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)
@@ -265,8 +203,7 @@ def boxplot_focal_contra(pkl_file_pre, pkl_file_inter, resultpath,
         patch.set_edgecolor('gray')
 
     sns.stripplot(ax=ax, x="binning_focus", y=quantity,
-                  data=mr_results, alpha=0.9, order=['focal',
-                                                     'contra-lateral'],
+                  data=mr_results, alpha=0.9, order=['focal','contra-lateral'],
                   size=2, dodge=True, palette=['darkblue', 'slategray'],
                   jitter=False)
 
@@ -305,17 +242,13 @@ def boxplot_focal_contra(pkl_file_pre, pkl_file_inter, resultpath,
     sns.despine(fig)
     fig.subplots_adjust(bottom=None, right=None,
                         wspace=0.8, hspace=0.5, left=0.15, top=1.2)
-    kmin = mr_results["kmin"].unique()[0]
-    kmax = mr_results["kmax"].unique()[0]
-    resultfile = '{}Fig1_boxplot_{}_kmin{}_kmax{}.pdf'.format(
-        resultpath,quantity, kmin, kmax)
+    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, excludeInconsistencies, **kwargs):
+def patientwise_pointplots(pkl_file_pre, pkl_file_inter, result_path):
     # -----------------------------------------------------------------#
     # Read data
     # -----------------------------------------------------------------#
@@ -327,22 +260,8 @@ def patientwise_pointplots(pkl_file_pre, pkl_file_inter,
                                                           axis=1)
     mr_results = pd.concat([mr_results_pre, mr_results_inter],
                            ignore_index=True)
-    mr_results = mr_results[(mr_results['fit_method'] == 'offset')]
-    conditions = list()
-    for key in kwargs.keys():
-        conditions.append(mr_results[key] == kwargs[key])
-    if len(conditions) > 0:
-        mr_results = mr_results[mr_utilities.conjunction(*conditions)]
-
-    if excludeInconsistencies:
-        rejected_inconsistencies = ['H_nonzero', 'H_linear']
-        for incon in rejected_inconsistencies:
-            mr_results = mr_results[[incon not in mr_results[
-                'inconsistencies'].tolist()[i] for i in range(
-                    len(mr_results['inconsistencies'].tolist()))]]
-
-    mr_results['logepsilon'] = mr_results.apply(lambda row: abs(
-        np.log(1-row['m'])), axis=1)
+    mr_results = mr_utilities.exclude_inconsistencies(mr_results)
+
     mr_results['epsilon'] = mr_results.apply(lambda row: 1-row['m'], axis=1)
     quantity = 'epsilon'
     # -----------------------------------------------------------------#
@@ -363,7 +282,8 @@ def patientwise_pointplots(pkl_file_pre, pkl_file_inter,
                     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')
+        stat, p_unpaired = stats.mannwhitneyu(all_SOZ, all_nSOZ,
+                                              alternative='two-sided')
 
         for rec in recordings:
             mr_rec = mr_patient[mr_patient['recording'] == rec]
@@ -377,27 +297,28 @@ def patientwise_pointplots(pkl_file_pre, pkl_file_inter,
                     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:
+            elif len(mr_rec[mr_rec['binning_focus'] == 'focal'].index) == 1:
 
-                q_SOZ = mr_rec[mr_rec['binning_focus']=='focal'][
+                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'][
+                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)
+                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),
+            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)
@@ -422,39 +343,20 @@ def patientwise_pointplots(pkl_file_pre, pkl_file_inter,
     sns.despine(fig)
     fig.subplots_adjust(bottom=None, right=None,
                         wspace=0.7, hspace=0.7, left=0.15, top=1.2)
-
-    kmin = mr_results_pre['kmin'].unique()[0]
-    kmax = mr_results_pre['kmax'].unique()[0]
-
-    plotfile_mx = '{}S2_{}_patients_kmin{}_kmax{}.pdf'.format(
-        result_path, quantity, kmin, kmax)
+    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,
-                                          excludeInconsistencies, **kwargs):
+                                          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_results[(mr_results['fit_method'] == 'offset')]
-    conditions = list()
-    for key in kwargs.keys():
-        conditions.append(mr_results[key] == kwargs[key])
-    if len(conditions) > 0:
-        mr_results = mr_results[mr_utilities.conjunction(*conditions)]
-
-    if excludeInconsistencies:
-        rejected_inconsistencies = ['H_nonzero', 'H_linear']
-        for incon in rejected_inconsistencies:
-            mr_results = mr_results[[incon not in mr_results[
-                'inconsistencies'].tolist()[i] for i in range(
-                    len(mr_results['inconsistencies'].tolist()))]]
+    mr_results = mr_utilities.exclude_inconsistencies(mr_results)
 
     # -----------------------------------------------------------------#
     # Write the brainregion, independent of hemisphere, into frame
@@ -467,56 +369,36 @@ def boxplot_brainregions_focus_onlypaired(pkl_file_pre, pkl_file_inter,
     # -----------------------------------------------------------------#
     def give_region_and_focus(row):
         regionname = row['binning'][1:]
-        focus = row['binning_focus'][0] # f for focal and c for contralateral
+        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['logtau'] = mr_results.apply(lambda row: np.log(row['tau']),
+    mr_results['region'] = mr_results.apply(lambda row: give_region(row),
                                             axis=1)
-    mr_results['logepsilon'] = mr_results.apply(lambda row: abs(np.log(
-        1-row['m'])), 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)
-    differences = dict({'A':[], 'H':[], 'PHC':[], 'EC':[]})
-
     for rec in mr_results.recording.unique():
-        recframe = mr_results[mr_results['recording']==rec]
+        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
-            focalval = area_frame[area_frame['binning_focus'] == 'focal'][
-                quantity].item()
-            contraval = \
-            area_frame[area_frame['binning_focus'] == 'contra-lateral'][
-                quantity].item()
-            differences[area].append(focalval - contraval)
             mr_results_paired = mr_results_paired.append(area_frame,
-                                                       ignore_index=True)
+                                                         ignore_index=True)
     mr_results = mr_results_paired
 
     # -----------------------------------------------------------------#
-    # Compute p-values of pairwise comparison
+    # Compute p-values
     # -----------------------------------------------------------------#
-    pvals = dict()
-    medians = dict()
-
-    for area in mr_results['region'].unique():
-        region_differences = differences[area]
-        stat, p_paired = stats.wilcoxon(region_differences)  # wilcoxon is by
-        # default two-sided
-        pvals[area] = p_paired
-        medians['f{}'.format(area)] = np.median(focalval)
-        medians['c{}'.format(area)] = np.median(contraval)
-
-    global_pvals, pvals_regions = hierarchical_model.hierarchical_model_regions(
+    global_pvals, pvals = hierarchical_model.hierarchical_model_regions(
         mr_results)
-    pvals = pvals_regions
 
     # -----------------------------------------------------------------#
     # Boxplot of regions
@@ -524,8 +406,8 @@ def boxplot_brainregions_focus_onlypaired(pkl_file_pre, pkl_file_inter,
     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)
+                  '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))
@@ -541,13 +423,13 @@ def boxplot_brainregions_focus_onlypaired(pkl_file_pre, pkl_file_inter,
     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'][
+        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())
@@ -568,7 +450,7 @@ def boxplot_brainregions_focus_onlypaired(pkl_file_pre, pkl_file_inter,
                     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,
+        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
@@ -578,37 +460,20 @@ def boxplot_brainregions_focus_onlypaired(pkl_file_pre, pkl_file_inter,
     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)
-    kmin = mr_results['kmin'].unique()[0]
-    kmax = mr_results['kmax'].unique()[0]
-    plotfile_mx = '{}Fig1_regions_{}_kmin{}_kmax{}_exclude_{}_onlypaired.pdf'.format(
-        result_path, quantity, kmin, kmax, excludeInconsistencies)
+    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,
-                               excludeInconsistencies, **kwargs):
+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_results[(mr_results['fit_method'] == 'offset')]
-    conditions = list()
-    for key in kwargs.keys():
-        conditions.append(mr_results[key] == kwargs[key])
-    if len(conditions) > 0:
-        mr_results = mr_results[mr_utilities.conjunction(*conditions)]
-
-    if excludeInconsistencies:
-        rejected_inconsistencies = ['H_nonzero', 'H_linear']
-        for incon in rejected_inconsistencies:
-            mr_results = mr_results[[incon not in mr_results[
-                'inconsistencies'].tolist()[i] for i in range(
-                len(mr_results['inconsistencies'].tolist()))]]
+    mr_results = mr_utilities.exclude_inconsistencies(mr_results)
 
     # -----------------------------------------------------------------#
     # Write the brainregion, independent of hemisphere, into frame
@@ -624,35 +489,19 @@ def boxplot_brainregions_focus(pkl_file_pre, pkl_file_inter, result_path,
         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['logtau'] = mr_results.apply(lambda row: np.log(row['tau']),axis=1)
-    mr_results['logepsilon'] = mr_results.apply(lambda row: -np.log(1-row['m']), axis=1)
+    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 of pairwise comparison
+    # Compute p-values
     # -----------------------------------------------------------------#
-    pvals = dict()
-    medians = dict()
-
-    for area in mr_results['region'].unique():
-        area_frame = mr_results[mr_results['region'] == area]
-        f_area_frame = area_frame[area_frame['binning_focus'] == 'focal']
-        c_area_frame = area_frame[area_frame['binning_focus'] ==
-                                 'contra-lateral']
-        focalval = f_area_frame[quantity].tolist()
-        contraval = c_area_frame[quantity].tolist()
-
-        stat, p = stats.mannwhitneyu(focalval, contraval,
-                                     alternative='two-sided')
-        pvals[area] = p
-        medians['f{}'.format(area)] = np.median(focalval)
-        medians['c{}'.format(area)] = np.median(contraval)
-
-    global_pvals, pvals_regions = hierarchical_model.hierarchical_model_regions(
-        mr_results)
-    pvals = pvals_regions
+    global_pvals, pvals = hierarchical_model.hierarchical_model_regions(
+                                    mr_results)
     # -----------------------------------------------------------------#
     # Boxplot of regions
     # -----------------------------------------------------------------#
@@ -660,8 +509,8 @@ def boxplot_brainregions_focus(pkl_file_pre, pkl_file_inter, result_path,
 
     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)
+                  '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()
@@ -678,12 +527,11 @@ def boxplot_brainregions_focus(pkl_file_pre, pkl_file_inter, result_path,
 
     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'][
+        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)
 
@@ -703,12 +551,11 @@ def boxplot_brainregions_focus(pkl_file_pre, pkl_file_inter, result_path,
     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,
+        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)
 
     # -----------------------------------------------------------------#
@@ -719,17 +566,13 @@ def boxplot_brainregions_focus(pkl_file_pre, pkl_file_inter, result_path,
     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)
-    kmin = mr_results['kmin'].unique()[0]
-    kmax = mr_results['kmax'].unique()[0]
-    plotfile_mx = '{}Fig1_regions_{}_kmin{}_kmax{}_exclude_{}.pdf'.format(
-        result_path, quantity, kmin, kmax, excludeInconsistencies)
+    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,
-                                     excludeInconsistencies, **kwargs):
+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)
@@ -738,19 +581,7 @@ def boxplot_brainregions_focus_Engel(pkl_file_pre,pkl_file_inter,result_path,
                            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_results[(mr_results['fit_method'] == 'offset')]
-    conditions = list()
-    for key in kwargs.keys():
-        conditions.append(mr_results[key] == kwargs[key])
-    if len(conditions) > 0:
-        mr_results = mr_results[mr_utilities.conjunction(*conditions)]
-
-    if excludeInconsistencies:
-        rejected_inconsistencies = ['H_nonzero', 'H_linear']
-        for incon in rejected_inconsistencies:
-            mr_results = mr_results[[incon not in mr_results[
-                'inconsistencies'].tolist()[i] for i in range(
-                    len(mr_results['inconsistencies'].tolist()))]]
+    mr_results = mr_utilities.exclude_inconsistencies(mr_results)
 
     # -----------------------------------------------------------------#
     # Write the brainregion, independent of hemisphere, into frame
@@ -763,38 +594,23 @@ def boxplot_brainregions_focus_Engel(pkl_file_pre,pkl_file_inter,result_path,
     # -----------------------------------------------------------------#
     def give_region_and_focus(row):
         regionname = row['binning'][1:]
-        focus = row['binning_focus'][0] # f for focal and c for contralateral
+        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['logtau'] = mr_results.apply(lambda row: np.log(row['tau']),axis=1)
-    mr_results['logepsilon'] = mr_results.apply(lambda row: abs(np.log(1-row['m'])), axis=1)
+    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 of pairwise comparison
+    # Compute p-values
     # -----------------------------------------------------------------#
-    pvals = dict()
-    medians = dict()
-
-    for area in mr_results['region'].unique():
-        area_frame = mr_results[mr_results['region'] == area]
-        f_area_frame = area_frame[area_frame['binning_focus'] == 'focal']
-        c_area_frame = area_frame[area_frame['binning_focus'] ==
-                                 'contra-lateral']
-        focalval = f_area_frame[quantity].tolist()
-        contraval = c_area_frame[quantity].tolist()
-
-        stat, p = stats.mannwhitneyu(focalval, contraval,
-                                     alternative='two-sided')
-        pvals[area] = p
-        medians['f{}'.format(area)] = np.median(focalval)
-        medians['c{}'.format(area)] = np.median(contraval)
-
-    global_pvals, pvals_regions = hierarchical_model.hierarchical_model_regions(
+    global_pvals, pvals = hierarchical_model.hierarchical_model_regions(
         mr_results)
-    pvals = pvals_regions
 
     # -----------------------------------------------------------------#
     # Boxplot of regions
@@ -849,7 +665,7 @@ def boxplot_brainregions_focus_Engel(pkl_file_pre,pkl_file_inter,result_path,
                     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,
+        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)
 
     # -----------------------------------------------------------------#
@@ -860,11 +676,8 @@ def boxplot_brainregions_focus_Engel(pkl_file_pre,pkl_file_inter,result_path,
     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)
-    kmin = mr_results['kmin'].unique()[0]
-    kmax = mr_results['kmax'].unique()[0]
 
-    plotfile_mx = '{}Fig1_regions_{}_kmin{}_kmax{}_exclude_{}_Engel1a.pdf'.format(
-        result_path, quantity, kmin, kmax, excludeInconsistencies)
+    plotfile_mx = '{}Fig1_regions_Engel1a.pdf'.format(result_path)
     fig.savefig(plotfile_mx, bbox_inches='tight')
     plt.show()
     plt.close('all')
@@ -873,21 +686,16 @@ def boxplot_brainregions_focus_Engel(pkl_file_pre,pkl_file_inter,result_path,
 def boxplot_hemispheres():
     result_path = '../Results/preictal/singleunits/'
     pkl_file = result_path + 'mr_results.pkl'
-    ksteps = (1, 400)
     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, excludeInconsistencies=True,
-                         kmin=ksteps[0], kmax=ksteps[1])
+                         result_path)
     boxplot_focal_contra_Engel(pkl_file_pre, pkl_file_inter,
-                         result_path, excludeInconsistencies=True,
-                         kmin=ksteps[0], kmax=ksteps[1])
-
+                         result_path)
     boxplot_focal_contra_onlypaired(pkl_file_pre, pkl_file_inter,
-                         result_path, excludeInconsistencies=True,
-                         kmin=ksteps[0], kmax=ksteps[1])
+                         result_path)
 
 
 def pointplot_patients_hemispheres():
@@ -896,10 +704,7 @@ def pointplot_patients_hemispheres():
     result_path_pointplot = '../Results/preictal/singleunits/'
 
     patientwise_pointplots(pkl_file_pre, pkl_file_inter,
-                           result_path_pointplot,
-                           excludeInconsistencies=True,
-                           kmin=1, kmax=400,
-                           fit_method='offset')
+                           result_path_pointplot)
 
 
 def boxplot_regions():
@@ -908,29 +713,11 @@ def boxplot_regions():
 
     result_path_inter = '../Results/interictal/singleunits_regions/'
     pkl_file_inter = result_path_inter + 'mr_results.pkl'
-    ksteps = (1, 400)
-
-    boxplot_brainregions_focus(pkl_file, pkl_file_inter,
-                                            result_path,
-                                            kmin=ksteps[0],
-                                            kmax=ksteps[1],
-                                            excludeInconsistencies=True,
-                                            fit_method='offset')
-
-    boxplot_brainregions_focus_Engel(pkl_file, pkl_file_inter,
-                                            result_path,
-                                            kmin=ksteps[0],
-                                            kmax=ksteps[1],
-                                            excludeInconsistencies=True,
-                                            fit_method='offset')
-
-    boxplot_brainregions_focus_onlypaired(pkl_file,
-                                            pkl_file_inter,
-                                            result_path,
-                                            kmin=ksteps[0],
-                                            kmax=ksteps[1],
-                                            excludeInconsistencies=True,
-                                            fit_method='offset')
+
+    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():
@@ -941,3 +728,4 @@ def main():
 
 if __name__ == '__main__':
     main()
+

python_code/analysis_PlosCB/create_figure_2.py

@@ -14,13 +14,13 @@ def cm2inch(value):
 
 
 def mt_examples(patient, save_path, fit_method, kmin, kmax,
-                windowsize, windowstep, dt, excludeInconsistencies=True):
+                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
+    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{}_' \
@@ -35,20 +35,10 @@ def mt_examples(patient, save_path, fit_method, kmin, kmax,
             # read m(t) data
             # -----------------------------------------------------------------#
             mt_frame = pd.read_pickle(mt_file)
-            if excludeInconsistencies:
-                rejected_inconsistencies = ['H_nonzero', 'H_linear']
-                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()))]]
+            mt_frame = mr_utilities.exclude_inconsistencies(mt_frame)
             if mt_frame.empty:
                 continue
 
-            mt_frame['logtau'] = mt_frame.apply(lambda row: np.log(row['tau']),
-                                                axis=1)
-            mt_frame['epsilon'] = mt_frame.apply(lambda row: 1-row['m'],
-                                                 axis=1)
             kmin = mt_frame.kmin.unique()[0]
             kmax = mt_frame.kmax.unique()[0]
             fit_method = mt_frame.fit_method.unique()[0]
@@ -89,111 +79,68 @@ def mt_examples(patient, save_path, fit_method, kmin, kmax,
     sns.despine()
     plt.subplots_adjust(bottom=None, right=None,
                         wspace=0.3, hspace=1.3, left=None, top=None)
-    plt.savefig(
-        '{}/Fig2_mt_{}_{}_winsize{}_winstep{}_kmin{}_kmax{}.pdf'.format(
-            save_path, patient, fit_method,
-            windowsize, windowstep, kmin, kmax), bbox_inches='tight')
+    plt.savefig('{}/Fig2_mt_patient{}.pdf'.format(save_path, patient),
+                bbox_inches='tight')
     plt.close()
 
 
-def boxplot_mx(result_path, x, ksteps, fit_method,
-                     excludeInconsistencies=True, Engel=False):
+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)
-
-    if excludeInconsistencies:
-        rejected_inconsistencies = ['H_nonzero', 'H_linear']
-        for incon in rejected_inconsistencies:
-            mx_allresults = mx_allresults[[incon not in mx_allresults[
-                'inconsistencies'].tolist()[i] for i in range(
-                    len(mx_allresults['inconsistencies'].tolist()))]]
+    mx_allresults = mr_utilities.exclude_inconsistencies(mx_allresults)
 
     x = mx_allresults['partno'].max()
-    mx_allresults['logtau'] = mx_allresults.apply(lambda row: np.log(row['tau']),
-                                                  axis=1)
-    mx_allresults['logepsilon'] = mx_allresults.apply(lambda row: abs(np.log(
-        1-row['m'])), axis=1)
-    mx_allresults['epsilon'] = mx_allresults.apply(lambda row:
-        1-row['m'], axis=1)
+    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_m_first_focal = []
-    all_m_last_focal = []
-    all_m_first_contra = []
-    all_m_last_contra = []
-
-    for ip, patient in enumerate(mx_allresults['patient'].unique()):
-        mx_patient = mx_allresults[mx_allresults['patient'] == patient]
-        recordings = mx_patient.recording.unique()
-        focal_diff_patient = []
-        contralateral_diff_patient = []
-        m_contra_first_patient = []
-        m_contra_last_patient = []
-        m_focal_first_patient = []
-        m_focal_last_patient = []
-
-        for recording in recordings:
-
-            mx_recording = mx_patient[mx_patient['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 = m_focal_last - m_focal_first
-                focal_diff_log = np.log(m_focal_last) - np.log(m_focal_first)
-                focal_diff_patient.append(focal_diff_log)
-                m_focal_first_patient.append(m_focal_first)
-                m_focal_last_patient.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 = m_contralateral_last - m_contralateral_first
-                contralateral_diff_log = np.log(m_contralateral_last) - \
+    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_patient.append(contralateral_diff_log)
-                m_contra_first_patient.append(m_contralateral_first)
-                m_contra_last_patient.append(m_contralateral_last)
-
-        if not len(focal_diff_patient) == 0:
-            focal_diff_list.append(focal_diff_patient)
-            all_m_first_focal.append(m_focal_first_patient)
-            all_m_last_focal.append(m_focal_last_patient)
-        if not len(contralateral_diff_patient) == 0:
-            contralateral_diff_list.append(contralateral_diff_patient)
-            all_m_first_contra.append(m_contra_first_patient)
-            all_m_last_contra.append(m_contra_last_patient)
-
-    # -----------------------------------------------------------------#
-    # Permutation test
-    # -----------------------------------------------------------------#
-    focal_diff_list = [item for sublist in focal_diff_list for item in sublist]
-    contralateral_diff_list = [item for sublist in contralateral_diff_list for item in sublist]
-    all_last_contra = [item for sublist in all_m_last_contra for item in sublist]
-    all_first_contra = [item for sublist in all_m_first_contra for item in
-                         sublist]
-    all_last_focal = [item for sublist in all_m_last_focal for item in
-                         sublist]
-    all_first_focal = [item for sublist in all_m_first_focal for item in
-                         sublist]
+            contralateral_diff_list.append(contralateral_diff_log)
+            all_first_contra.append(m_contralateral_first)
+            all_last_contra.append(m_contralateral_last)
 
     # -----------------------------------------------------------------#
     # Wilcoxon paired test
@@ -271,7 +218,6 @@ def boxplot_mx(result_path, x, ksteps, fit_method,
     # annotate p-values
     # -----------------------------------------------------------------#
     y_max = 2.4*ymax
-    y_min = ymin
     ax.annotate("", xy=(-0.2, y_max), xycoords='data',
                 xytext=(0.2, y_max), textcoords='data',
                 arrowprops=dict(arrowstyle="-", ec='black',
@@ -292,36 +238,22 @@ def boxplot_mx(result_path, x, ksteps, fit_method,
     # -----------------------------------------------------------------#
     # Saving
     # -----------------------------------------------------------------#
-    kmin = mx_allresults['kmin'].unique()[0]
-    kmax = mx_allresults['kmax'].unique()[0]
-    plotfile_mx = '{}Fig2_mx{}_{}_kmin{}_kmax{}_Engel{}.pdf'.format(
-                   result_path, x, measure, kmin, kmax, Engel)
+    plotfile_mx = '{}Fig2_2parts_Engel{}.pdf'.format(result_path, Engel)
     plt.savefig(plotfile_mx, bbox_inches='tight')
 
 
-def pointplot_mx_regions_Engel(region, result_path, x, ksteps,
-                               fit_method, excludeInconsistencies=True):
+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)
-
-    if excludeInconsistencies:
-        rejected_inconsistencies = ['H_nonzero', 'H_linear']
-        for incon in rejected_inconsistencies:
-            mx_allresults = mx_allresults[[incon not in mx_allresults[
-                'inconsistencies'].tolist()[i] for i in range(
-                    len(mx_allresults['inconsistencies'].tolist()))]]
+    mx_allresults = mr_utilities.exclude_inconsistencies(mx_allresults)
 
     x = mx_allresults['partno'].max()
-    mx_allresults['logtau'] = mx_allresults.apply(lambda row: np.log(row['tau']),
-                                                  axis=1)
-    mx_allresults['logepsilon'] = mx_allresults.apply(lambda row: abs(np.log(
-        1-row['m'])), axis=1)
-    mx_allresults['epsilon'] = mx_allresults.apply(lambda row:
-        1-row['m'], axis=1)
+    mx_allresults['epsilon'] = mx_allresults.apply(lambda row: 1-row['m'],
+                                                   axis=1)
     measure = 'epsilon'
 
     def give_region(row):
@@ -329,10 +261,9 @@ def pointplot_mx_regions_Engel(region, result_path, x, ksteps,
 
     mx_allresults['region'] = mx_allresults.apply(lambda row: give_region(row),
                                                   axis=1)
-    if not region == 'all':
-        mx_allresults = mx_allresults[mx_allresults['region'] == region]
+
     # -----------------------------------------------------------------#
-    # write estimates to lists
+    # write paired estimates to lists
     # -----------------------------------------------------------------#
     fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(4*1.2)))
     df = None
@@ -340,12 +271,9 @@ def pointplot_mx_regions_Engel(region, result_path, x, ksteps,
     mx_allresults = mx_allresults[(mx_allresults['patient'].isin(patients1a))]
     regions = ['H', 'A', 'PHC', 'EC']
     labels = []
-    pvals = dict()
     for i, reg in enumerate(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 = []
@@ -362,30 +290,20 @@ def pointplot_mx_regions_Engel(region, result_path, x, ksteps,
                 m_focal_first = mx_focal[mx_focal['partno'] == 1][
                     measure].item()
                 m_focal_last = mx_focal[mx_focal['partno'] == x][measure].item()
-                focal_diff = m_focal_last - m_focal_first
                 all_first_focal.append(m_focal_first)
                 all_last_focal.append(m_focal_last)
-                focal_diff_list.append(focal_diff)
                 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 = m_contralateral_last - m_contralateral_first
-
+                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)
-                contralateral_diff_list.append(contralateral_diff)
-
-        # -----------------------------------------------------------------#
-        # 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
@@ -432,7 +350,7 @@ def pointplot_mx_regions_Engel(region, result_path, x, ksteps,
     # -----------------------------------------------------------------#
     # Plotting
     # -----------------------------------------------------------------#
-    df = df[df['binning_focus']=='SOZ']
+    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,
@@ -467,34 +385,20 @@ def pointplot_mx_regions_Engel(region, result_path, x, ksteps,
     # -----------------------------------------------------------------#
     # Saving
     # -----------------------------------------------------------------#
-    kmin = mx_allresults['kmin'].unique()[0]
-    kmax = mx_allresults['kmax'].unique()[0]
-    plotfile_mx = '{}Fig2_pointmxallregions{}_{}_kmin{}_kmax{}_Engel1a.pdf'.format(
-                   result_path, x, measure, kmin, kmax, region)
+    plotfile_mx = '{}Fig2_2parts_regions_Engel1a.pdf'.format(result_path)
     plt.savefig(plotfile_mx, bbox_inches='tight')
 
 
-def boxplot_mx_regions(region, result_path, x, ksteps, fit_method,
-                       excludeInconsistencies=True):
+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)
-
-    if excludeInconsistencies:
-        rejected_inconsistencies = ['H_nonzero', 'H_linear']
-        for incon in rejected_inconsistencies:
-            mx_allresults = mx_allresults[[incon not in mx_allresults[
-                'inconsistencies'].tolist()[i] for i in range(
-                len(mx_allresults['inconsistencies'].tolist()))]]
+    mx_allresults = mr_utilities.exclude_inconsistencies(mx_allresults)
 
     x = mx_allresults['partno'].max()
-    mx_allresults['logtau'] = mx_allresults.apply(lambda row: np.log(row['tau']),
-                                                  axis=1)
-    mx_allresults['logepsilon'] = mx_allresults.apply(lambda row: abs(np.log(
-        1-row['m'])), axis=1)
     mx_allresults['epsilon'] = mx_allresults.apply(lambda row:
         1-row['m'], axis=1)
     measure = 'epsilon'
@@ -504,9 +408,10 @@ def boxplot_mx_regions(region, result_path, x, ksteps, fit_method,
 
     mx_allresults['region'] = mx_allresults.apply(lambda row: give_region(row),
                                                   axis=1)
-
-    if not region == 'all':
-        mx_allresults = mx_allresults[mx_allresults['region'] == region]
+    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
@@ -536,7 +441,6 @@ def boxplot_mx_regions(region, result_path, x, ksteps, fit_method,
                 m_focal_first = mx_focal[mx_focal['partno'] == 1][
                     measure].item()
                 m_focal_last = mx_focal[mx_focal['partno'] == x][measure].item()
-                focal_diff = m_focal_last - m_focal_first
                 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)
@@ -547,7 +451,6 @@ def boxplot_mx_regions(region, result_path, x, ksteps, fit_method,
                     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 = m_contralateral_last - m_contralateral_first
                 contralateral_diff_log = np.log(m_contralateral_last) - \
                                      np.log(m_contralateral_first)
                 all_first_contra.append(m_contralateral_first)
@@ -620,7 +523,7 @@ def boxplot_mx_regions(region, result_path, x, ksteps, fit_method,
     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'],
+                    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)
 
@@ -655,8 +558,7 @@ def boxplot_mx_regions(region, result_path, x, ksteps, fit_method,
                         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)
+                    0.0, 'p' + '={:.3f}'.format(pvals[region]), fontsize=8.5)
 
     ax.set_ylim(ymin, ymax)
     plt.ylabel(quantity_label)
@@ -665,10 +567,8 @@ def boxplot_mx_regions(region, result_path, x, ksteps, fit_method,
     # -----------------------------------------------------------------#
     # Saving
     # -----------------------------------------------------------------#
-    kmin = mx_allresults['kmin'].unique()[0]
-    kmax = mx_allresults['kmax'].unique()[0]
-    plotfile_mx = '{}Fig2_mxallregions{}_{}_kmin{}_kmax{}.pdf'.format(
-        result_path, x, measure, kmin, kmax, region)
+    plotfile_mx = '{}Fig2_2parts_regions_Engel{}.pdf'.format(result_path,
+                                                             Engel)
     plt.savefig(plotfile_mx, bbox_inches='tight')
 
 
@@ -701,7 +601,7 @@ def binning_label(patient, soz):
 
 
 def get_mt(patient, recording, soz, save_path, fit_method, kmin, kmax,
-            windowsize, windowstep, dt, excludeInconsistencies=True):
+            windowsize, windowstep, dt):
 
     binning = binning_label(patient, soz)
     mt_file = '{}{}/mt_results_{}_{}_kmin{}_kmax{}_winsize{}_winstep{}_dt{}.pkl'.format(
@@ -715,19 +615,20 @@ def get_mt(patient, recording, soz, save_path, fit_method, kmin, kmax,
     # read m(t) data
     # -----------------------------------------------------------------#
     mt_frame = pd.read_pickle(mt_file)
-    if excludeInconsistencies:
-        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, it should be 260
-            # windows long. if it is much shorter, it is unclear where to
-            # divide it.
-            raise Exception('Too many inconsistencies')
+    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
 
@@ -740,7 +641,7 @@ 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]
+                15, 16, 17, 18, 19, 20]
     if Engel:
         patients = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
     laterality = ["SOZ", "nSOZ"]
@@ -757,19 +658,22 @@ def variance_mt(save_path, data_path, fit_method, kmin, kmax,
             for irec, rec in enumerate(recordings):
                 try:
                     m_list = get_mt(patient, rec, soz, data_path, fit_method,
-                                    kmin, kmax, windowsize, windowstep, dt,
-                                    excludeInconsistencies=True)
+                                    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_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
 
         # -----------------------------------------------------------------#
@@ -834,32 +738,23 @@ def variance_mt(save_path, data_path, fit_method, kmin, kmax,
                 arrowprops=dict(arrowstyle="-", ec='black',
                                 connectionstyle="bar,fraction=0.2"))
 
-    if quantity == 'epsilon':
-        ax.text(-0.25, y_max * 0.6,
-                'p={:.3f}'.format(df[df['binning_focus']=='SOZ'].pval.tolist()[0]), fontsize=8)
-    else:
-        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.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"))
 
-    if quantity == 'epsilon':
-        ax.text(0.75, y_max * 0.6,
-                'p={:.3f}'.format(df[df['binning_focus']=='nSOZ'].pval.tolist()[0]), fontsize=8)
-    else:
-        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)
+    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_Engel{}_{}_winsize{}_winstep{}_kmin{}_kmax{}_epsilon{}.pdf'.format(
-            save_path, Engel, fit_method, windowsize, windowstep, kmin,
-            kmax, logepsilon), bbox_inches='tight')
+    plt.savefig('{}Fig2_variance_mt.pdf'.format(save_path), bbox_inches='tight')
     plt.show()
 
 
@@ -874,7 +769,7 @@ def timeresolved_example():
     patients = [13]
     for patient in patients:
         mt_examples(patient, result_path, 'offset', kmin, kmax,
-                    windowsize, windowstep, dt, excludeInconsistencies=True)
+                    windowsize, windowstep, dt)
 
 
 def boxplot_2parts_hemispheres():
@@ -882,10 +777,8 @@ def boxplot_2parts_hemispheres():
     ksteps = (1, 400)
     fit_method = 'offset'
     result_path = '../Results/preictal/singleunits/'
-    boxplot_mx(result_path, x, ksteps, fit_method,
-                excludeInconsistencies=True, Engel=False)
-    boxplot_mx(result_path, x, ksteps, fit_method,
-                excludeInconsistencies=True, Engel=True)
+    boxplot_mx(result_path, x, ksteps, fit_method, Engel=False)
+    boxplot_mx(result_path, x, ksteps, fit_method, Engel=True)
 
 
 def boxplot_2parts_regions():
@@ -893,12 +786,8 @@ def boxplot_2parts_regions():
     ksteps = (1, 400)
     fit_method = 'offset'
     result_path = '../Results/preictal/singleunits_regions/'
-    region = 'all'
-    boxplot_mx_regions(region, result_path, x, ksteps, fit_method,
-                              excludeInconsistencies=True)
-
-    pointplot_mx_regions_Engel(region, result_path, x, ksteps,
-                                      fit_method, excludeInconsistencies=True)
+    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():
@@ -911,10 +800,11 @@ def visualize_mt_variance():
 
 def main():
     timeresolved_example()
-    #boxplot_2parts_hemispheres()
-    #boxplot_2parts_regions()
+    boxplot_2parts_hemispheres()
+    boxplot_2parts_regions()
     visualize_mt_variance()
 
 
 if __name__ == '__main__':
     main()
+

python_code/analysis_PlosCB/create_supp_figures.py

@@ -30,8 +30,8 @@ def distribution_A_recordings(resultpath):
     A_nonzero_frac = []
 
     for recording in recordings:
-
-        outputfile = '{}{}/activity_SU_4ms.pkl'.format(data_path, recording)
+        outputfile = '{}{}/pkl/activity_SU_4ms.pkl'.format(data_path,
+                                                           recording)
 
         if not os.path.isfile(outputfile):
             print("No binning data for recording ", recording)
@@ -39,7 +39,6 @@ def distribution_A_recordings(resultpath):
 
         data = pd.read_pickle(outputfile)
         binnings = data['binning'].tolist()
-
         for binning in binnings:
             activity = data[data['binning'] == binning]['activity'].tolist()[0]
             A_mean.append(np.mean(activity))
@@ -124,8 +123,7 @@ def mt_all_supp(save_path,fit_method, kmin, kmax, windowsize, windowstep, dt):
         for j, soz in enumerate(laterality):
             ax = plt.Subplot(fig, inner[j])
             ax = mt_plot(ax, patient, soz, save_path,fit_method, kmin, kmax,
-                          windowsize, windowstep, dt,
-                          excludeInconsistencies=True)
+                          windowsize, windowstep, dt)
 
             if j == 0:
                 ax.set_title('Patient {}'.format(int(patient)), fontsize=12)
@@ -136,19 +134,7 @@ def mt_all_supp(save_path,fit_method, kmin, kmax, windowsize, windowstep, dt):
 
             fig.add_subplot(ax)
     sns.despine(fig)
-    plt.savefig(
-        '{}/S6_allmt_{}_{}_winsize{}_winstep{}_kmin{}_kmax{}.pdf'.format(
-            save_path, patient, fit_method,
-            windowsize, windowstep, kmin, kmax), bbox_inches='tight')
-
-    plt.savefig(
-        '{}/S6_allmt_{}_{}_winsize{}_winstep{}_kmin{}_kmax{}.png'.format(
-            save_path, patient, fit_method,
-            windowsize, windowstep, kmin, kmax), bbox_inches='tight', dpi=350)
-    plt.savefig(
-        '{}/S6_allmt_{}_{}_winsize{}_winstep{}_kmin{}_kmax{}.svg'.format(
-            save_path, patient, fit_method,
-            windowsize, windowstep, kmin, kmax), bbox_inches='tight')
+    plt.savefig('{}/S6_allmt.pdf'.format(save_path), bbox_inches='tight')
     plt.show()
     plt.close()
 
@@ -182,7 +168,7 @@ def binning_label(patient, soz):
 
 
 def mt_plot(ax, patient, soz, save_path, fit_method, kmin, kmax,
-            windowsize, windowstep, dt, excludeInconsistencies=True):
+            windowsize, windowstep, dt):
 
     recordings = mr_utilities.list_preictrecordings(patient)[0]
     binning = binning_label(patient, soz)
@@ -201,17 +187,17 @@ def mt_plot(ax, patient, soz, save_path, fit_method, kmin, kmax,
         # read m(t) data
         # -----------------------------------------------------------------#
         mt_frame = pd.read_pickle(mt_file)
-        if excludeInconsistencies:
-            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:
-                continue
+
+        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:
+            continue
 
         kmin = mt_frame.kmin.unique()[0]
         kmax = mt_frame.kmax.unique()[0]
@@ -253,7 +239,6 @@ def mt_plot(ax, patient, soz, save_path, fit_method, kmin, kmax,
 
     ax.xaxis.grid(False)
     ax.set_xlabel("time (s)")
-
     ax.tick_params(top='off', bottom='off', left='off', right='off',
                    labelleft='on', labelbottom='off')
     ax.tick_params(top='off', bottom='on', left='off', right='off',
@@ -265,7 +250,7 @@ def mt_plot(ax, patient, soz, save_path, fit_method, kmin, kmax,
 def windowsize_analysis_recs(data_path, recording, binning, outputfilename,
                              resultpath):
 
-    outputfile = '{}{}/{}'.format(data_path, recording, outputfilename)
+    outputfile = '{}{}/pkl/{}'.format(data_path, recording, outputfilename)
 
     if not os.path.isfile(outputfile):
         print("No binning data for recording ", recording)
@@ -358,15 +343,18 @@ def windowsize_analysis():
     outputfilename = 'activity_SU_4ms.pkl'
     resultpath = '../Results/preictal/singleunits/'
 
-    data_path = '../Data/interictal/'
-    recording = '13ref'
-    windowsize_analysis_recs(data_path, recording, binning, outputfilename,
-                             resultpath)
+    distribution_A_recordings(resultpath)
 
-    data_path = '../Data/preictal/'
-    recording = '13_02'
-    windowsize_analysis_recs(data_path, recording, binning, outputfilename,
-                             resultpath)
+    # warning: analysis below has long computation time
+    # data_path = '../Data/interictal/'
+    # recording = '13ref'
+    # windowsize_analysis_recs(data_path, recording, binning, outputfilename,
+    #                          resultpath)
+    #
+    # data_path = '../Data/preictal/'
+    # recording = '13_02'
+    # windowsize_analysis_recs(data_path, recording, binning, outputfilename,
+    #                          resultpath)
 
 
 def all_mt():
@@ -381,23 +369,16 @@ def all_mt():
 
 
 def example_autocorrelations(patient, pkl_file_inter, pkl_file_pre,
-                             resultpath, **kwargs):
+                             resultpath):
     mr_results_pre = pd.read_pickle(pkl_file_pre)
     mr_results_inter = pd.read_pickle(pkl_file_inter)
     mr_results_pre['rectype'] = 'preictal'
     mr_results_inter['rectype'] = 'interictal'
     mr_results = pd.concat([mr_results_pre, mr_results_inter],
                            ignore_index=True)
-    conditions = list()
-    for key in kwargs.keys():
-        conditions.append(mr_results[key] == kwargs[key])
-    if len(conditions) > 0:
-        mr_results = mr_results[mr_utilities.conjunction(*conditions)]
 
     mr_results = mr_results[mr_results['patient'] == patient]
     binning_foci = ['focal', 'contra-lateral']
-    kmin = mr_results.kmin.unique()[0]
-    kmax = mr_results.kmax.unique()[0]
 
     sns.set(style="white")
     fig, ax = plt.subplots(2, 4, figsize=(cm2inch(24), cm2inch(6)))
@@ -460,11 +441,7 @@ def example_autocorrelations(patient, pkl_file_inter, pkl_file_pre,
     plt.subplots_adjust(bottom=None, right=None,
                         wspace=0.5, hspace=1, left=0.15, top=1.2)
     sns.despine(fig)
-    resultpath_rec = '{}S1_autocorrelations_kmin{}_kmax{}_patient{}.pdf'.format(resultpath, kmin, kmax, patient)
-    plt.savefig(resultpath_rec, bbox_inches='tight')
-    resultpath_rec = '{}S1_autocorrelations_kmin{}_kmax{}_patient{}.png'.format(resultpath, kmin, kmax, patient)
-    plt.savefig(resultpath_rec, bbox_inches='tight', dpi=300)
-    resultpath_rec = '{}S1_autocorrelations_kmin{}_kmax{}_patient{}.svg'.format(resultpath, kmin, kmax, patient)
+    resultpath_rec = '{}S1_autocorrelations_patient{}.pdf'.format(resultpath, patient)
     plt.savefig(resultpath_rec, bbox_inches='tight')
     plt.close()
 
@@ -473,20 +450,18 @@ def create_example_ACF_plot():
     pkl_file_pre = '../Results/preictal/singleunits/mr_results.pkl'
     pkl_file_inter = '../Results/interictal/singleunits/mr_results.pkl'
     resultpath = '../Results/preictal/singleunits/'
-    ksteps = (1, 400)
     patient = 12
-    example_autocorrelations(patient, pkl_file_inter, pkl_file_pre,
-         resultpath, kmin=ksteps[0], kmax=ksteps[1])
+    example_autocorrelations(patient, pkl_file_inter, pkl_file_pre, resultpath)
     patient = 13
-    example_autocorrelations(patient, pkl_file_inter, pkl_file_pre,
-         resultpath, kmin=ksteps[0], kmax=ksteps[1])
+    example_autocorrelations(patient, pkl_file_inter, pkl_file_pre, resultpath)
 
 
 def main():
-    #create_example_ACF_plot()
-    #windowsize_analysis()
+    create_example_ACF_plot()
+    windowsize_analysis()
     all_mt()
 
 
 if __name__ == '__main__':
     main()
+

python_code/analysis_PlosCB/mr_epilepsy_timeresolved.py

@@ -13,7 +13,7 @@ def analyse_timeresolved_mr(data_path, result_path, recordings,
                             fit_method="offset"):
 
     for recording in recordings:
-        outputfile = '{}{}/activity_SU_4ms.pkl'.format(data_path, recording)
+        outputfile = '{}{}/pkl/activity_SU_4ms.pkl'.format(data_path, recording)
 
         if os.path.isfile(outputfile):
             data = pd.read_pickle(outputfile)
@@ -222,4 +222,4 @@ def mt_x_parts(x, outputfilename, data_path, result_path, recordings,
         result_path, x, ksteps[0], ksteps[1], fit_method)
 
     if mx_allresults is not None:
-        mx_allresults.to_pickle(datafile_mx)
+        mx_allresults.to_pickle(datafile_mx)

python_code/analysis_PlosCB/mr_utilities.py

@@ -1,11 +1,18 @@
-import sys
-sys.path.insert(0, '../mrestimator/')
 import numpy as np
 import functools
 import mrestimator as mre
 from scipy import stats
 
 
+def exclude_inconsistencies(mr_results):
+    rejected_inconsistencies = ['H_nonzero', 'H_linear']
+    for incon in rejected_inconsistencies:
+        mr_results = mr_results[[incon not in mr_results[
+            'inconsistencies'].tolist()[i] for i in range(
+            len(mr_results['inconsistencies'].tolist()))]]
+    return mr_results
+
+
 def consistency_checks(activity, rk, fitres):
     """ Checks the consistency of MR estimates
     We only use H_nonzero and H_linear as exclusion criteria.
@@ -235,3 +242,4 @@ def binning_focus(binning, patient):
         binning_focus = 'contra-lateral'
 
     return binning_focus
+

python_code/analysis_PlosCB/run_analysis_fullrecs.py

@@ -24,7 +24,7 @@ def run_mr_regions():
         # -----------------------------------------------------------------------#
         # run mr analysis
         # -----------------------------------------------------------------------#
-        outputfilename = 'activity_brainregions_SU_4ms.pkl'
+        outputfilename = 'pkl/activity_brainregions_SU_4ms.pkl'
 
         mr_epilepsy_fullrec.mr_analysis_confint(data_path, result_path,
                     outputfilename=outputfilename, recordings='all',
@@ -52,7 +52,7 @@ def run_mr_hemispheres():
         # -----------------------------------------------------------------------#
         # run mr analysis
         # -----------------------------------------------------------------------#
-        outputfilename = 'activity_SU_4ms.pkl'
+        outputfilename = 'pkl/activity_SU_4ms.pkl'
 
         mr_epilepsy_fullrec.mr_analysis_confint(data_path, result_path,
                             outputfilename=outputfilename, recordings='all',

python_code/analysis_PlosCB/run_analysis_timeresolved.py

@@ -24,7 +24,7 @@ def run_2parts_regions():
     # -----------------------------------------------------------------------#
     # run mr analysis
     # -----------------------------------------------------------------------#
-    outputfilename = 'activity_brainregions_SU_4ms.pkl'
+    outputfilename = 'pkl/activity_brainregions_SU_4ms.pkl'
     x = 2
     mr_epilepsy_timeresolved.mt_x_parts(x, outputfilename, data_path,
        result_path, recordings=recordings, binnings='all', ksteps=(1, 400),
@@ -54,7 +54,7 @@ def run_2parts_hemispheres():
     # run timeresolved mr analysis
     # -----------------------------------------------------------------------#
     x = 2
-    outputfilename = "activity_SU_4ms.pkl"
+    outputfilename = "pkl/activity_SU_4ms.pkl"
     mr_epilepsy_timeresolved.mt_x_parts(x, outputfilename, data_path,
        result_path, recordings=recordings, binnings=['left', 'right'],
        ksteps=(1, 400), fit_method="offset")
@@ -103,9 +103,9 @@ def run_timeresolved_hemispheres():
 
 
 def main():
-    run_timeresolved_hemispheres()
-    #run_2parts_hemispheres()
-    #run_2parts_regions()
+    run_2parts_hemispheres()
+    run_2parts_regions()
+    run_timeresolved_hemispheres()  # warning: timeresolved analysis is computationally costly
 
 
 if __name__ == '__main__':