diff --git a/actions/stimulus-lfp-response-mec-no-zscore/attributes.yaml b/actions/stimulus-lfp-response-mec-no-zscore/attributes.yaml new file mode 100644 index 000000000..ae0308d0b --- /dev/null +++ b/actions/stimulus-lfp-response-mec-no-zscore/attributes.yaml @@ -0,0 +1,4 @@ +registered: '2020-11-28T15:05:27' +data: + notebook: 20_stimulus-lfp-response.ipynb + html: 20_stimulus-lfp-response.html diff --git a/actions/stimulus-lfp-response-mec-no-zscore/data/20_stimulus-lfp-response-mec.html b/actions/stimulus-lfp-response-mec-no-zscore/data/20_stimulus-lfp-response-mec.html new file mode 100644 index 000000000..4262e4830 --- /dev/null +++ b/actions/stimulus-lfp-response-mec-no-zscore/data/20_stimulus-lfp-response-mec.html @@ -0,0 +1,14624 @@ + + +
+ +%load_ext autoreload
+%autoreload 2
+
import os
+import expipe
+import pathlib
+import numpy as np
+import spatial_maps.stats as stats
+import septum_mec
+import septum_mec.analysis.data_processing as dp
+import septum_mec.analysis.registration
+import head_direction.head as head
+import spatial_maps as sp
+import speed_cells.speed as spd
+import re
+import joblib
+import multiprocessing
+import shutil
+import psutil
+import pandas as pd
+import matplotlib.pyplot as plt
+import matplotlib
+from distutils.dir_util import copy_tree
+from neo import SpikeTrain
+import scipy
+import seaborn as sns
+from tqdm.notebook import tqdm_notebook as tqdm
+tqdm.pandas()
+
+from spike_statistics.core import permutation_resampling_test
+
+from spikewaveform.core import calculate_waveform_features_from_template, cluster_waveform_features
+
+from septum_mec.analysis.plotting import violinplot, despine
+
#############################
+
+perform_zscore = False
+
+if not perform_zscore:
+ zscore_str = "-no-zscore"
+else:
+ zscore_str = ""
+
+#################################
+
%matplotlib inline
+plt.rc('axes', titlesize=12)
+plt.rcParams.update({
+ 'font.size': 12,
+ 'figure.figsize': (6, 4),
+ 'figure.dpi': 150
+})
+
+output_path = pathlib.Path("output") / ("stimulus-lfp-response-mec" + zscore_str)
+(output_path / "statistics").mkdir(exist_ok=True, parents=True)
+(output_path / "figures").mkdir(exist_ok=True, parents=True)
+output_path.mkdir(exist_ok=True)
+
data_loader = dp.Data()
+actions = data_loader.actions
+project = data_loader.project
+
identify_neurons = actions['identify-neurons']
+sessions = pd.read_csv(identify_neurons.data_path('sessions'))
+
lfp_action = actions['stimulus-lfp-response' + zscore_str]
+lfp_results = pd.read_csv(lfp_action.data_path('results'))
+
lfp_results = pd.merge(sessions, lfp_results, how='left')
+
lfp_results = lfp_results.query('stim_location!="ms"')
+
def action_group(row):
+ a = int(row.channel_group in [0,1,2,3])
+ return f'{row.action}-{a}'
+lfp_results['action_side_a'] = lfp_results.apply(action_group, axis=1)
+
lfp_results['stim_strength'] = lfp_results['stim_p_max'] / lfp_results['theta_energy']
+
# lfp_results_hemisphere = lfp_results.sort_values(
+# by=['action_side_a', 'stim_strength', 'signal_to_noise'], ascending=[True, False, False]
+lfp_results_hemisphere = lfp_results.sort_values(
+ by=['action_side_a', 'channel_group'], ascending=[True, False]
+).drop_duplicates(subset='action_side_a', keep='first')
+lfp_results_hemisphere.loc[:,['action_side_a','channel_group', 'signal_to_noise', 'stim_strength']].head()
+
colors = ['#1b9e77','#d95f02','#7570b3','#e7298a']
+labels = ['Baseline I', '11 Hz', 'Baseline II', '30 Hz']
+# Hz11 means that the baseline session was indeed before an 11 Hz session
+queries = ['baseline and i and Hz11', 'frequency==11', 'baseline and ii and Hz30', 'frequency==30']
+
# prepare pairwise comparison: same animal same side same date different sessions
+
def make_entity_date_side(row):
+ s = row.action_side_a.split('-')
+ del s[2]
+ return '-'.join(s)
+
lfp_results_hemisphere['entity_date_side'] = lfp_results_hemisphere.apply(make_entity_date_side, axis=1)
+
density = True
+cumulative = True
+histtype = 'step'
+lw = 2
+if perform_zscore:
+ bins = {
+ 'theta_energy': np.arange(0, .7, .03),
+ 'theta_peak': np.arange(0, .7, .03),
+ 'theta_freq': np.arange(4, 10, .5),
+ 'theta_half_width': np.arange(0, 15, .5)
+ }
+else:
+ bins = {
+ 'theta_energy': np.arange(0, .008, .0003),
+ 'theta_peak': np.arange(0, .007, .0003),
+ 'theta_freq': np.arange(4, 12, .5),
+ 'theta_half_width': np.arange(0, 15, .5)
+ }
+xlabel = {
+ 'theta_energy': 'Theta energy (dB)',
+ 'theta_peak': 'Peak PSD (dB/Hz)',
+ 'theta_freq': '(Hz)',
+ 'theta_half_width': '(Hz)',
+}
+# key = 'theta_energy'
+# key = 'theta_peak'
+results = {}
+for key in bins:
+ results[key] = list()
+ fig = plt.figure(figsize=(3.5,2))
+ plt.suptitle(key)
+ legend_lines = []
+ for color, query, label in zip(colors, queries, labels):
+ values = lfp_results_hemisphere.query(query).loc[:,['entity_date_side', key]]
+ results[key].append(values.rename({key: label}, axis=1))
+ values[key].hist(
+ bins=bins[key], density=density, cumulative=cumulative, lw=lw,
+ histtype=histtype, color=color)
+ legend_lines.append(matplotlib.lines.Line2D([0], [0], color=color, lw=lw, label=label))
+
+ plt.legend(
+ handles=legend_lines,
+ bbox_to_anchor=(1.04,1), borderaxespad=0, frameon=False)
+ plt.tight_layout()
+ plt.grid(False)
+ plt.xlim(right=bins[key].max() - bins[key].max()*0.025)
+ despine()
+ plt.xlabel(xlabel[key])
+ figname = f'lfp-psd-histogram-{key}'
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.png',
+ bbox_inches='tight', transparent=True)
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.svg',
+ bbox_inches='tight', transparent=True)
+
density = True
+cumulative = True
+histtype = 'step'
+lw = 2
+if perform_zscore:
+ bins = {
+ 'stim_energy': np.arange(0, .7, .01),
+ 'stim_half_width': np.arange(0, 10, .5),
+ 'stim_p_max': np.arange(0, 4, .01),
+ 'stim_strength': np.arange(0, 160, 1)
+ }
+else:
+ bins = {
+ 'stim_energy': np.arange(0, .008, .0001),
+ 'stim_half_width': np.arange(0, 0.5, .001),
+ 'stim_p_max': np.arange(0, .06, .0001),
+ 'stim_strength': np.arange(0, 160, 1)
+ }
+xlabel = {
+ 'stim_energy': 'Energy (dB)',
+ 'stim_half_width': '(Hz)',
+ 'stim_p_max': 'Peak PSD (dB/Hz)',
+ 'stim_strength': 'Ratio',
+}
+# key = 'theta_energy'
+# key = 'theta_peak'
+for key in bins:
+ results[key] = list()
+ fig = plt.figure(figsize=(3.2,2))
+ plt.suptitle(key)
+ legend_lines = []
+ for color, query, label in zip(colors[1::2], queries[1::2], labels[1::2]):
+ values = lfp_results_hemisphere.query(query).loc[:,['entity_date_side', key]]
+ results[key].append(values.rename({key: label}, axis=1))
+ values[key].hist(
+ bins=bins[key], density=density, cumulative=cumulative, lw=lw,
+ histtype=histtype, color=color)
+ legend_lines.append(matplotlib.lines.Line2D([0], [0], color=color, lw=lw, label=label))
+
+ plt.legend(
+ handles=legend_lines,
+ bbox_to_anchor=(1.04,1), borderaxespad=0, frameon=False)
+ plt.tight_layout()
+ plt.grid(False)
+ plt.xlim(right=bins[key].max() - bins[key].max()*0.025)
+ despine()
+ plt.xlabel(xlabel[key])
+ figname = f'lfp-psd-histogram-{key}'
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.png',
+ bbox_inches='tight', transparent=True)
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.svg',
+ bbox_inches='tight', transparent=True)
+
from functools import reduce
+
for key, val in results.items():
+ df = reduce(lambda left,right: pd.merge(left, right, on='entity_date_side', how='outer'), val)
+ results[key] = df.drop('entity_date_side',axis=1)
+
def summarize(data):
+ return "{:.1e} ± {:.1e} ({})".format(data.mean(), data.sem(), sum(~np.isnan(data)))
+
+
+def MWU(df, keys):
+ '''
+ Mann Whitney U
+ '''
+ Uvalue, pvalue = scipy.stats.mannwhitneyu(
+ df[keys[0]].dropna(),
+ df[keys[1]].dropna(),
+ alternative='two-sided')
+
+ return "{:.2f}, {:.3f}".format(Uvalue, pvalue)
+
+
+def PRS(df, keys):
+ '''
+ Permutation ReSampling
+ '''
+ pvalue, observed_diff, diffs = permutation_resampling(
+ df[keys[0]].dropna(),
+ df[keys[1]].dropna(), statistic=np.median)
+
+ return "{:.2f}, {:.3f}".format(observed_diff, pvalue)
+
+
+def wilcoxon(df, keys):
+ dff = df.loc[:,keys].dropna()
+ if len(dff[keys].dropna()) == 0:
+ statistic, pvalue = np.nan, np.nan
+ else:
+ statistic, pvalue = scipy.stats.wilcoxon(
+ dff[keys[0]],
+ dff[keys[1]],
+ alternative='two-sided')
+
+ return "{:.2f}, {:.3f}, ({})".format(statistic, pvalue, len(dff))
+
+
+def paired_t(df, keys):
+ dff = df.loc[:,[keys[0], keys[1]]].dropna()
+ statistic, pvalue = scipy.stats.ttest_rel(
+ dff[keys[0]],
+ dff[keys[1]])
+
+ return "{:.2f}, {:.3f}".format(statistic, pvalue)
+
+
+def normality(df, key):
+ if len(df[key].dropna()) < 8:
+ statistic, pvalue = np.nan, np.nan
+ else:
+ statistic, pvalue = scipy.stats.normaltest(
+ df[key].dropna())
+
+ return "{:.2f}, {:.3f}".format(statistic, pvalue)
+
+
+def shapiro(df, key):
+ if len(df[key].dropna()) < 8:
+ statistic, pvalue = np.nan, np.nan
+ else:
+ statistic, pvalue = scipy.stats.shapiro(
+ df[key].dropna())
+
+ return "{:.2f}, {:.3f}".format(statistic, pvalue)
+
+def rename(name):
+ return name.replace("_field", "-field").replace("_", " ").capitalize()
+
results
+
stat = pd.DataFrame()
+
+for key, df in results.items():
+ Key = rename(key)
+ stat[Key] = df.agg(summarize)
+ stat[Key] = df.agg(summarize)
+
+ for i, c1 in enumerate(df.columns):
+ stat.loc[f'Normality {c1}', Key] = normality(df, c1)
+# stat.loc[f'Shapiro {c1}', Key] = shapiro(df, c1)
+ for c2 in df.columns[i+1:]:
+# stat.loc[f'MWU {c1} {c2}', Key] = MWU(df, [c1, c2])
+# stat.loc[f'PRS {c1} {c2}', Key] = PRS(df, [c1, c2])
+ stat.loc[f'Wilcoxon {c1} {c2}', Key] = wilcoxon(df, [c1, c2])
+# stat.loc[f'Paired T {c1} {c2}', Key] = paired_t(df, [c1, c2])
+
+stat.sort_index()
+
stat.to_latex(output_path / "statistics" / "statistics.tex")
+stat.to_csv(output_path / "statistics" / "statistics.csv")
+
psd = pd.read_feather(pathlib.Path("output") / ("stimulus-lfp-response" + zscore_str) / 'data' / 'psd.feather')
+freqs = pd.read_feather(pathlib.Path("output") / ("stimulus-lfp-response" + zscore_str) / 'data' / 'freqs.feather')
+
from septum_mec.analysis.plotting import plot_bootstrap_timeseries
+
freq = freqs.T.iloc[0].values
+
+mask = (freq < 49)
+
fig, axs = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(5,2))
+axs = axs.repeat(2)
+for i, (ax, query) in enumerate(zip(axs.ravel(), queries)):
+ selection = [
+ f'{r.action}_{r.channel_group}'
+ for i, r in lfp_results_hemisphere.query(query).iterrows()]
+ values = psd.loc[mask, selection].to_numpy()
+ values = 10 * np.log10(values)
+ plot_bootstrap_timeseries(freq[mask], values, ax=ax, lw=1, label=labels[i], color=colors[i])
+# ax.set_title(titles[i])
+ ax.set_xlabel('Frequency Hz')
+ ax.legend(frameon=False)
+axs[0].set_ylabel('PSD (dB/Hz)')
+# axs[0].set_ylim(-31, 1)
+despine()
+
+figname = 'lfp-psd'
+fig.savefig(
+ output_path / 'figures' / f'{figname}.png',
+ bbox_inches='tight', transparent=True)
+fig.savefig(
+ output_path / 'figures' / f'{figname}.svg',
+ bbox_inches='tight', transparent=True)
+
action = project.require_action("stimulus-lfp-response-mec" + zscore_str)
+
copy_tree(output_path, str(action.data_path()))
+
septum_mec.analysis.registration.store_notebook(action, "20_stimulus-lfp-response-mec.ipynb")
+
+
\n", + " | action_side_a | \n", + "channel_group | \n", + "signal_to_noise | \n", + "stim_strength | \n", + "
---|---|---|---|---|
71 | \n", + "1833-010719-1-0 | \n", + "7 | \n", + "0.001902 | \n", + "NaN | \n", + "
67 | \n", + "1833-010719-1-1 | \n", + "3 | \n", + "0.003522 | \n", + "NaN | \n", + "
583 | \n", + "1833-020719-1-0 | \n", + "7 | \n", + "-0.002942 | \n", + "NaN | \n", + "
579 | \n", + "1833-020719-1-1 | \n", + "3 | \n", + "0.012323 | \n", + "NaN | \n", + "
375 | \n", + "1833-020719-3-0 | \n", + "7 | \n", + "-0.002042 | \n", + "NaN | \n", + "
\n", + " | Theta energy | \n", + "Theta peak | \n", + "Theta freq | \n", + "Theta half width | \n", + "Stim energy | \n", + "Stim half width | \n", + "Stim p max | \n", + "Stim strength | \n", + "
---|---|---|---|---|---|---|---|---|
11 Hz | \n", + "1.8e-03 ± 2.8e-04 (8) | \n", + "1.4e-03 ± 2.6e-04 (8) | \n", + "8.1e+00 ± 5.3e-02 (8) | \n", + "8.9e-01 ± 7.4e-02 (8) | \n", + "8.4e-04 ± 3.6e-04 (8) | \n", + "2.2e+00 ± 1.4e+00 (8) | \n", + "1.9e-03 ± 7.0e-04 (8) | \n", + "1.2e+00 ± 5.6e-01 (8) | \n", + "
30 Hz | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "
Baseline I | \n", + "2.3e-03 ± 2.1e-04 (50) | \n", + "1.8e-03 ± 1.8e-04 (50) | \n", + "7.8e+00 ± 5.9e-02 (50) | \n", + "1.1e+00 ± 1.8e-01 (49) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Baseline II | \n", + "2.3e-03 ± 2.4e-04 (32) | \n", + "1.8e-03 ± 2.3e-04 (32) | \n", + "8.1e+00 ± 4.7e-02 (32) | \n", + "9.1e-01 ± 3.9e-02 (31) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Normality 11 Hz | \n", + "8.13, 0.017 | \n", + "6.61, 0.037 | \n", + "6.30, 0.043 | \n", + "0.23, 0.890 | \n", + "3.45, 0.178 | \n", + "7.42, 0.024 | \n", + "1.81, 0.405 | \n", + "6.43, 0.040 | \n", + "
Normality 30 Hz | \n", + "nan, nan | \n", + "nan, nan | \n", + "nan, nan | \n", + "nan, nan | \n", + "nan, nan | \n", + "nan, nan | \n", + "nan, nan | \n", + "nan, nan | \n", + "
Normality Baseline I | \n", + "41.72, 0.000 | \n", + "31.09, 0.000 | \n", + "29.47, 0.000 | \n", + "81.74, 0.000 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Normality Baseline II | \n", + "13.17, 0.001 | \n", + "20.78, 0.000 | \n", + "0.96, 0.618 | \n", + "13.33, 0.001 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Wilcoxon 11 Hz 30 Hz | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "
Wilcoxon 11 Hz Baseline II | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Wilcoxon Baseline I 11 Hz | \n", + "5.00, 0.069, (8) | \n", + "2.00, 0.025, (8) | \n", + "0.00, 0.034, (8) | \n", + "6.00, 0.093, (8) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Wilcoxon Baseline I 30 Hz | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Wilcoxon Baseline I Baseline II | \n", + "264.00, 1.000, (32) | \n", + "256.00, 0.881, (32) | \n", + "0.00, 0.000, (32) | \n", + "203.00, 0.544, (30) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Wilcoxon Baseline II 30 Hz | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "nan, nan, (0) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
%load_ext autoreload
+%autoreload 2
+
import os
+import expipe
+import pathlib
+import numpy as np
+import spatial_maps.stats as stats
+import septum_mec
+import septum_mec.analysis.data_processing as dp
+import septum_mec.analysis.registration
+import head_direction.head as head
+import spatial_maps as sp
+import speed_cells.speed as spd
+import re
+import joblib
+import multiprocessing
+import shutil
+import psutil
+import pandas as pd
+import matplotlib.pyplot as plt
+import matplotlib
+from distutils.dir_util import copy_tree
+from neo import SpikeTrain
+import scipy
+import statsmodels
+import seaborn as sns
+from tqdm.notebook import tqdm_notebook as tqdm
+tqdm.pandas()
+
+from spike_statistics.core import permutation_resampling_test
+
+from spikewaveform.core import calculate_waveform_features_from_template, cluster_waveform_features
+
+from septum_mec.analysis.plotting import violinplot, despine
+
#############################
+
+zscore_str = "-no-zscore"
+# zscore_str = ""
+
+stim_loc = 'mec'
+# stim_loc = 'ms'
+
+#################################
+
%matplotlib inline
+plt.rc('axes', titlesize=12)
+plt.rcParams.update({
+ 'font.size': 12,
+ 'figure.figsize': (6, 4),
+ 'figure.dpi': 150
+})
+
+output_path = pathlib.Path("output") / ("stimulus-lfp-response" + '-' + stim_loc + zscore_str)
+(output_path / "statistics").mkdir(exist_ok=True, parents=True)
+(output_path / "figures").mkdir(exist_ok=True, parents=True)
+output_path.mkdir(exist_ok=True)
+
data_loader = dp.Data()
+actions = data_loader.actions
+project = data_loader.project
+
identify_neurons = actions['identify-neurons']
+sessions = pd.read_csv(identify_neurons.data_path('sessions'))
+
lfp_action = actions['stimulus-lfp-response' + zscore_str]
+lfp_results = pd.read_csv(lfp_action.data_path('results'))
+
lfp_results = pd.merge(sessions, lfp_results, how='left')
+
if stim_loc == 'ms':
+ lfp_results = lfp_results.query('stim_location!="mecl" and stim_location!="mecr"')
+elif stim_loc == 'mec':
+ lfp_results = lfp_results.query('stim_location!="ms"')
+else:
+ raise AssertionError('')
+
def action_group(row):
+ a = int(row.channel_group in [0,1,2,3])
+ return f'{row.action}-{a}'
+lfp_results['action_side_a'] = lfp_results.apply(action_group, axis=1)
+
lfp_results['stim_strength'] = lfp_results['stim_p_max'] / lfp_results['theta_bandpower']
+
# lfp_results_hemisphere = lfp_results.sort_values(
+# by=['action_side_a', 'stim_strength', 'signal_to_noise'], ascending=[True, False, False]
+lfp_results_hemisphere = lfp_results.sort_values(
+ by=['action_side_a', 'channel_group'], ascending=[True, False]
+).drop_duplicates(subset='action_side_a', keep='first')
+lfp_results_hemisphere.loc[:,['action_side_a','channel_group', 'signal_to_noise', 'stim_strength']].head()
+
colors = ['#1b9e77','#d95f02','#7570b3','#e7298a']
+labels = ['Baseline I', '11 Hz', 'Baseline II', '30 Hz']
+# Hz11 means that the baseline session was indeed before an 11 Hz session
+queries = ['baseline and i and Hz11', 'frequency==11', 'baseline and ii and Hz30', 'frequency==30']
+
# prepare pairwise comparison: same animal same side same date different sessions
+
def make_entity_date_side(row):
+ s = row.action_side_a.split('-')
+ del s[2]
+ return '-'.join(s)
+
lfp_results_hemisphere['entity_date_side'] = lfp_results_hemisphere.apply(make_entity_date_side, axis=1)
+
from functools import reduce
+
keys = [
+ 'theta_bandpower',
+ 'theta_relpower',
+ 'theta_relpeak',
+ 'theta_peak',
+ 'theta_freq',
+ 'theta_half_width',
+ 'stim_bandpower',
+ 'stim_relpower',
+ 'stim_relpeak',
+ 'stim_half_width',
+ 'stim_p_max',
+ 'stim_strength',
+]
+
+results = {}
+for key in keys:
+ results[key] = list()
+ for query, label in zip(queries, labels):
+ values = lfp_results_hemisphere.query(query).loc[:,['entity_date_side', key]]
+ results[key].append(values.rename({key: label}, axis=1))
+
+for key, val in results.items():
+ df = reduce(lambda left,right: pd.merge(left, right, on='entity_date_side', how='outer'), val)
+ results[key] = df.drop('entity_date_side', axis=1)
+
xlabel = {
+ 'theta_bandpower': 'Theta bandpower (V$^2$/Hz)',
+ 'theta_relpower': 'Theta relative power',
+ 'theta_relpeak': 'Theta relative power',
+ 'theta_peak': 'Peak PSD (V$^2$/Hz)',
+ 'theta_freq': '(Hz)',
+ 'theta_half_width': '(Hz)',
+}
+
+for key in xlabel:
+ fig = plt.figure(figsize=(3.5,2))
+ plt.suptitle(key)
+ legend_lines = []
+ for color, label in zip(colors, labels):
+ legend_lines.append(matplotlib.lines.Line2D([0], [0], color=color, label=label))
+ sns.kdeplot(data=results[key].loc[:,labels], cumulative=True, legend=False, palette=colors, common_norm=False)
+ plt.legend(
+ handles=legend_lines,
+ bbox_to_anchor=(1.04,1), borderaxespad=0, frameon=False)
+ plt.tight_layout()
+ plt.grid(False)
+ despine()
+ plt.xlabel(xlabel[key])
+ figname = f'lfp-psd-histogram-{key}'
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.png',
+ bbox_inches='tight', transparent=True)
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.svg',
+ bbox_inches='tight', transparent=True)
+
xlabel = {
+ 'stim_bandpower': 'Energy (V$^2$/Hz)',
+ 'stim_relpower': 'Relative power',
+ 'stim_relpeak': 'Relative power',
+ 'stim_half_width': '(Hz)',
+ 'stim_p_max': 'Peak PSD (V$^2$/Hz)',
+ 'stim_strength': 'Ratio',
+}
+for key in xlabel:
+ fig = plt.figure(figsize=(3.2,2))
+ plt.suptitle(key)
+ legend_lines = []
+ for color, label in zip(colors[1::2], labels[1::2]):
+ legend_lines.append(matplotlib.lines.Line2D([0], [0], color=color, label=label))
+ sns.kdeplot(data=results[key].loc[:, labels[1::2]], cumulative=True, legend=False, palette=colors[1::2], common_norm=False)
+ plt.legend(
+ handles=legend_lines,
+ bbox_to_anchor=(1.04,1), borderaxespad=0, frameon=False)
+ plt.tight_layout()
+ plt.grid(False)
+ despine()
+ plt.xlabel(xlabel[key])
+ figname = f'lfp-psd-histogram-{key}'
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.png',
+ bbox_inches='tight', transparent=True)
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.svg',
+ bbox_inches='tight', transparent=True)
+
def summarize(data):
+ return "{:.1e} ± {:.1e} ({})".format(data.mean(), data.sem(), sum(~np.isnan(data)))
+
+
+def MWU(df, keys):
+ '''
+ Mann Whitney U
+ '''
+ Uvalue, pvalue = scipy.stats.mannwhitneyu(
+ df[keys[0]].dropna(),
+ df[keys[1]].dropna(),
+ alternative='two-sided')
+
+ return "{:.2f}, {:.3f}".format(Uvalue, pvalue)
+
+
+def PRS(df, keys):
+ '''
+ Permutation ReSampling
+ '''
+ pvalue, observed_diff, diffs = permutation_resampling(
+ df[keys[0]].dropna(),
+ df[keys[1]].dropna(), statistic=np.median)
+
+ return "{:.2f}, {:.3f}".format(observed_diff, pvalue)
+
+
+def wilcoxon(df, keys):
+ dff = df.loc[:,[keys[0], keys[1]]].dropna()
+ statistic, pvalue = scipy.stats.wilcoxon(
+ dff[keys[0]],
+ dff[keys[1]],
+ alternative='two-sided')
+
+ return "{:.1e}, {:.1e}, ({})".format(statistic, pvalue, len(dff))
+
+
+def summarize_wilcoxon(df, keys):
+ dff = df.loc[:,[keys[0], keys[1]]].dropna()
+
+
+ return"{:.1e} ± {:.1e}, {:.1e} ± {:.1e} ({})".format(dff[keys[0]].mean(), dff[keys[0]].sem(), dff[keys[1]].mean(), dff[keys[1]].sem(), len(dff))
+
+
+def paired_t(df, keys):
+ dff = df.loc[:,[keys[0], keys[1]]].dropna()
+ statistic, pvalue = scipy.stats.ttest_rel(
+ dff[keys[0]],
+ dff[keys[1]])
+
+ return "{:.2f}, {:.3f}".format(statistic, pvalue)
+
+
+def normality(df, key):
+ statistic, pvalue = scipy.stats.normaltest(
+ df[key].dropna())
+
+ return "{:.1e}, {:.1e}".format(statistic, pvalue)
+
+
+def shapiro(df, key):
+ statistic, pvalue = scipy.stats.shapiro(
+ df[key].dropna())
+
+ return "{:.2f}, {:.3f}".format(statistic, pvalue)
+
+def rename(name):
+ return name.replace("_field", "-field").replace("_", " ").capitalize()
+
stat = pd.DataFrame()
+
+for key, df in results.items():
+ Key = rename(key)
+ stat[Key] = df.agg(summarize)
+ stat[Key] = df.agg(summarize)
+
+ for i, c1 in enumerate(df.columns):
+ try:
+ stat.loc[f'Normality {c1}', Key] = normality(df, c1)
+ except:
+ stat.loc[f'Normality {c1}', Key] = np.nan
+# stat.loc[f'Shapiro {c1}', Key] = shapiro(df, c1)
+ for c2 in df.columns[i+1:]:
+# stat.loc[f'MWU {c1} {c2}', Key] = MWU(df, [c1, c2])
+# stat.loc[f'PRS {c1} {c2}', Key] = PRS(df, [c1, c2])
+ try:
+ stat.loc[f'Wilcoxon {c1} {c2}', Key] = wilcoxon(df, [c1, c2])
+ except:
+ stat.loc[f'Wilcoxon {c1} {c2}', Key] = np.nan
+# stat.loc[f'Paired T {c1} {c2}', Key] = paired_t(df, [c1, c2])
+ stat.loc[f'Paired summary {c1} {c2}', Key] = summarize_wilcoxon(df, [c1, c2])
+
+stat.sort_index()
+
stat.to_latex(output_path / "statistics" / f"statistics.tex")
+stat.to_csv(output_path / "statistics" / f"statistics.csv")
+
for key, result in results.items():
+ result.to_latex(output_path / "statistics" / f"values_{key}.tex")
+ result.to_csv(output_path / "statistics" / f"values_{key}.csv")
+
psd = pd.read_feather(pathlib.Path("output") / ("stimulus-lfp-response" + zscore_str) / 'data' / 'psd.feather')
+freqs = pd.read_feather(pathlib.Path("output") / ("stimulus-lfp-response" + zscore_str) / 'data' / 'freqs.feather')
+
from septum_mec.analysis.plotting import plot_bootstrap_timeseries
+
freq = freqs.T.iloc[0].values
+
+mask = (freq < 49)
+
fig, axs = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(5,2))
+axs = axs.repeat(2)
+for i, (ax, query) in enumerate(zip(axs.ravel(), queries)):
+ selection = [
+ f'{r.action}_{r.channel_group}'
+ for i, r in lfp_results_hemisphere.query(query).iterrows()]
+ values = psd.loc[mask, selection].to_numpy()
+ values = 10 * np.log10(values)
+ plot_bootstrap_timeseries(freq[mask], values, ax=ax, lw=1, label=labels[i], color=colors[i])
+# ax.set_title(titles[i])
+ ax.set_xlabel('Frequency Hz')
+ ax.legend(frameon=False)
+axs[0].set_ylabel('PSD (dB/Hz)')
+# axs[0].set_ylim(-31, 1)
+despine()
+
+figname = 'lfp-psd'
+fig.savefig(
+ output_path / 'figures' / f'{figname}.png',
+ bbox_inches='tight', transparent=True)
+fig.savefig(
+ output_path / 'figures' / f'{figname}.svg',
+ bbox_inches='tight', transparent=True)
+
action = project.require_action("stimulus-lfp-response" + '-' + stim_loc + zscore_str)
+
copy_tree(output_path, str(action.data_path()))
+
septum_mec.analysis.registration.store_notebook(action, "20_stimulus-lfp-response.ipynb")
+
+
\n", + " | action_side_a | \n", + "channel_group | \n", + "signal_to_noise | \n", + "stim_strength | \n", + "
---|---|---|---|---|
71 | \n", + "1833-010719-1-0 | \n", + "7 | \n", + "0.001902 | \n", + "NaN | \n", + "
67 | \n", + "1833-010719-1-1 | \n", + "3 | \n", + "0.003522 | \n", + "NaN | \n", + "
695 | \n", + "1833-010719-2-0 | \n", + "7 | \n", + "0.004280 | \n", + "1.401239 | \n", + "
691 | \n", + "1833-010719-2-1 | \n", + "3 | \n", + "0.003974 | \n", + "3.920680 | \n", + "
583 | \n", + "1833-020719-1-0 | \n", + "7 | \n", + "-0.002942 | \n", + "NaN | \n", + "
\n", + " | Theta bandpower | \n", + "Theta relpower | \n", + "Theta relpeak | \n", + "Theta peak | \n", + "Theta freq | \n", + "Theta half width | \n", + "Stim bandpower | \n", + "Stim relpower | \n", + "Stim relpeak | \n", + "Stim half width | \n", + "Stim p max | \n", + "Stim strength | \n", + "
---|---|---|---|---|---|---|---|---|---|---|---|---|
11 Hz | \n", + "8.5e-04 ± 8.0e-05 (44) | \n", + "8.6e-02 ± 4.6e-03 (44) | \n", + "2.3e-01 ± 5.2e-02 (44) | \n", + "3.2e-04 ± 3.1e-05 (44) | \n", + "7.7e+00 ± 1.3e-01 (44) | \n", + "1.7e+00 ± 3.0e-01 (43) | \n", + "9.6e-04 ± 8.0e-05 (44) | \n", + "1.1e-01 ± 7.4e-03 (44) | \n", + "1.6e+01 ± 1.9e+00 (44) | \n", + "3.4e-01 ± 2.5e-03 (44) | \n", + "2.1e-03 ± 2.3e-04 (44) | \n", + "3.3e+00 ± 3.8e-01 (44) | \n", + "
30 Hz | \n", + "5.3e-04 ± 6.3e-05 (34) | \n", + "5.3e-02 ± 4.2e-03 (34) | \n", + "3.0e-02 ± 5.3e-02 (34) | \n", + "1.9e-04 ± 2.7e-05 (34) | \n", + "7.4e+00 ± 1.4e-01 (34) | \n", + "1.2e+01 ± 1.9e+00 (34) | \n", + "1.8e-03 ± 3.3e-04 (34) | \n", + "1.6e-01 ± 2.3e-02 (34) | \n", + "8.6e+01 ± 1.4e+01 (34) | \n", + "3.0e-01 ± 2.3e-03 (23) | \n", + "5.3e-03 ± 1.1e-03 (34) | \n", + "1.3e+01 ± 2.8e+00 (34) | \n", + "
Baseline I | \n", + "2.2e-03 ± 2.2e-04 (46) | \n", + "2.6e-01 ± 1.6e-02 (46) | \n", + "6.5e+00 ± 6.7e-01 (46) | \n", + "1.7e-03 ± 1.8e-04 (46) | \n", + "7.8e+00 ± 4.0e-02 (46) | \n", + "7.7e-01 ± 1.8e-02 (46) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "
Baseline II | \n", + "2.2e-03 ± 2.3e-04 (32) | \n", + "2.7e-01 ± 1.7e-02 (32) | \n", + "6.3e+00 ± 7.4e-01 (32) | \n", + "1.6e-03 ± 2.1e-04 (32) | \n", + "8.1e+00 ± 4.7e-02 (32) | \n", + "8.4e-01 ± 3.2e-02 (32) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "nan ± nan (0) | \n", + "
Normality 11 Hz | \n", + "2.1e+01, 2.6e-05 | \n", + "5.6e+00, 6.2e-02 | \n", + "1.2e+01, 2.3e-03 | \n", + "2.4e+01, 5.3e-06 | \n", + "1.9e+00, 3.9e-01 | \n", + "1.0e+01, 5.5e-03 | \n", + "1.7e+01, 1.8e-04 | \n", + "3.9e+00, 1.4e-01 | \n", + "5.8e+00, 5.6e-02 | \n", + "1.6e+01, 4.2e-04 | \n", + "1.5e+01, 5.8e-04 | \n", + "1.2e+01, 2.3e-03 | \n", + "
Normality 30 Hz | \n", + "3.1e+01, 1.9e-07 | \n", + "7.9e+00, 2.0e-02 | \n", + "1.2e+01, 2.3e-03 | \n", + "3.8e+01, 5.3e-09 | \n", + "7.2e+00, 2.8e-02 | \n", + "1.9e+02, 2.2e-41 | \n", + "1.7e+01, 1.7e-04 | \n", + "1.5e+01, 4.8e-04 | \n", + "6.2e+00, 4.5e-02 | \n", + "6.1e-01, 7.4e-01 | \n", + "1.9e+01, 6.3e-05 | \n", + "4.3e+01, 5.1e-10 | \n", + "
Normality Baseline I | \n", + "4.3e+01, 5.8e-10 | \n", + "1.6e+00, 4.6e-01 | \n", + "2.1e+00, 3.4e-01 | \n", + "3.2e+01, 1.3e-07 | \n", + "5.9e+00, 5.3e-02 | \n", + "1.9e+00, 3.8e-01 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Normality Baseline II | \n", + "1.4e+01, 1.1e-03 | \n", + "3.6e-01, 8.3e-01 | \n", + "4.9e+00, 8.8e-02 | \n", + "2.5e+01, 3.4e-06 | \n", + "4.7e+00, 9.7e-02 | \n", + "1.6e+01, 2.8e-04 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Paired summary 11 Hz 30 Hz | \n", + "8.2e-04 ± 8.8e-05, 5.4e-04 ± 6.6e-05 (32) | \n", + "8.4e-02 ± 4.5e-03, 5.5e-02 ± 4.2e-03 (32) | \n", + "2.0e-01 ± 5.0e-02, 4.9e-02 ± 5.4e-02 (32) | \n", + "3.0e-04 ± 3.1e-05, 1.9e-04 ± 2.8e-05 (32) | \n", + "7.6e+00 ± 1.5e-01, 7.3e+00 ± 1.4e-01 (32) | \n", + "1.5e+00 ± 3.5e-01, 1.1e+01 ± 2.0e+00 (31) | \n", + "9.7e-04 ± 9.3e-05, 1.8e-03 ± 3.4e-04 (32) | \n", + "1.1e-01 ± 8.6e-03, 1.6e-01 ± 2.4e-02 (32) | \n", + "1.8e+01 ± 2.4e+00, 8.5e+01 ± 1.4e+01 (32) | \n", + "3.4e-01 ± 3.6e-03, 3.0e-01 ± 2.3e-03 (21) | \n", + "2.1e-03 ± 2.7e-04, 5.4e-03 ± 1.1e-03 (32) | \n", + "3.5e+00 ± 4.9e-01, 1.3e+01 ± 3.0e+00 (32) | \n", + "
Paired summary 11 Hz Baseline II | \n", + "8.2e-04 ± 8.8e-05, 2.2e-03 ± 2.3e-04 (32) | \n", + "8.4e-02 ± 4.5e-03, 2.7e-01 ± 1.7e-02 (32) | \n", + "2.0e-01 ± 5.0e-02, 6.3e+00 ± 7.4e-01 (32) | \n", + "3.0e-04 ± 3.1e-05, 1.6e-03 ± 2.1e-04 (32) | \n", + "7.6e+00 ± 1.5e-01, 8.1e+00 ± 4.7e-02 (32) | \n", + "1.5e+00 ± 3.5e-01, 8.6e-01 ± 2.5e-02 (31) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "
Paired summary Baseline I 11 Hz | \n", + "2.2e-03 ± 2.3e-04, 8.5e-04 ± 8.0e-05 (44) | \n", + "2.6e-01 ± 1.6e-02, 8.6e-02 ± 4.6e-03 (44) | \n", + "6.3e+00 ± 6.7e-01, 2.3e-01 ± 5.2e-02 (44) | \n", + "1.6e-03 ± 1.9e-04, 3.2e-04 ± 3.1e-05 (44) | \n", + "7.8e+00 ± 4.2e-02, 7.7e+00 ± 1.3e-01 (44) | \n", + "7.8e-01 ± 1.8e-02, 1.7e+00 ± 3.0e-01 (43) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "
Paired summary Baseline I 30 Hz | \n", + "2.3e-03 ± 2.8e-04, 5.4e-04 ± 6.6e-05 (32) | \n", + "2.7e-01 ± 1.8e-02, 5.5e-02 ± 4.2e-03 (32) | \n", + "6.6e+00 ± 8.1e-01, 4.9e-02 ± 5.4e-02 (32) | \n", + "1.7e-03 ± 2.3e-04, 1.9e-04 ± 2.8e-05 (32) | \n", + "7.8e+00 ± 4.4e-02, 7.3e+00 ± 1.4e-01 (32) | \n", + "7.6e-01 ± 2.0e-02, 1.2e+01 ± 2.0e+00 (32) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "
Paired summary Baseline I Baseline II | \n", + "2.3e-03 ± 2.8e-04, 2.2e-03 ± 2.3e-04 (32) | \n", + "2.7e-01 ± 1.8e-02, 2.7e-01 ± 1.7e-02 (32) | \n", + "6.6e+00 ± 8.1e-01, 6.3e+00 ± 7.4e-01 (32) | \n", + "1.7e-03 ± 2.3e-04, 1.6e-03 ± 2.1e-04 (32) | \n", + "7.8e+00 ± 4.4e-02, 8.1e+00 ± 4.7e-02 (32) | \n", + "7.6e-01 ± 2.0e-02, 8.4e-01 ± 3.2e-02 (32) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "
Paired summary Baseline II 30 Hz | \n", + "2.2e-03 ± 2.3e-04, 5.4e-04 ± 6.6e-05 (32) | \n", + "2.7e-01 ± 1.7e-02, 5.5e-02 ± 4.2e-03 (32) | \n", + "6.3e+00 ± 7.4e-01, 4.9e-02 ± 5.4e-02 (32) | \n", + "1.6e-03 ± 2.1e-04, 1.9e-04 ± 2.8e-05 (32) | \n", + "8.1e+00 ± 4.7e-02, 7.3e+00 ± 1.4e-01 (32) | \n", + "8.4e-01 ± 3.2e-02, 1.2e+01 ± 2.0e+00 (32) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "nan ± nan, nan ± nan (0) | \n", + "
Wilcoxon 11 Hz 30 Hz | \n", + "1.1e+02, 3.3e-03, (32) | \n", + "4.1e+01, 3.0e-05, (32) | \n", + "1.3e+02, 1.4e-02, (32) | \n", + "1.2e+02, 5.6e-03, (32) | \n", + "1.2e+02, 4.5e-02, (32) | \n", + "6.7e+01, 3.9e-04, (31) | \n", + "2.1e+02, 2.9e-01, (32) | \n", + "2.0e+02, 2.2e-01, (32) | \n", + "3.0e+01, 1.2e-05, (32) | \n", + "0.0e+00, 9.5e-07, (21) | \n", + "1.6e+02, 5.0e-02, (32) | \n", + "8.8e+01, 1.0e-03, (32) | \n", + "
Wilcoxon 11 Hz Baseline II | \n", + "1.2e+01, 2.5e-06, (32) | \n", + "2.0e+00, 9.6e-07, (32) | \n", + "0.0e+00, 8.0e-07, (32) | \n", + "3.0e+00, 1.1e-06, (32) | \n", + "1.2e+02, 9.3e-03, (32) | \n", + "2.3e+02, 6.7e-01, (31) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Wilcoxon Baseline I 11 Hz | \n", + "3.5e+01, 7.9e-08, (44) | \n", + "1.0e+00, 8.2e-09, (44) | \n", + "2.0e+00, 8.7e-09, (44) | \n", + "3.0e+00, 9.4e-09, (44) | \n", + "3.6e+02, 4.7e-01, (44) | \n", + "4.3e+02, 6.2e-01, (43) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Wilcoxon Baseline I 30 Hz | \n", + "6.0e+00, 1.4e-06, (32) | \n", + "0.0e+00, 8.0e-07, (32) | \n", + "0.0e+00, 8.0e-07, (32) | \n", + "0.0e+00, 8.0e-07, (32) | \n", + "9.5e+01, 1.6e-03, (32) | \n", + "7.1e+01, 3.1e-04, (32) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Wilcoxon Baseline I Baseline II | \n", + "2.4e+02, 7.1e-01, (32) | \n", + "2.4e+02, 7.1e-01, (32) | \n", + "2.4e+02, 5.9e-01, (32) | \n", + "2.3e+02, 5.5e-01, (32) | \n", + "6.0e+00, 9.0e-06, (32) | \n", + "1.4e+02, 2.3e-02, (32) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
Wilcoxon Baseline II 30 Hz | \n", + "1.6e+01, 3.5e-06, (32) | \n", + "0.0e+00, 8.0e-07, (32) | \n", + "0.0e+00, 8.0e-07, (32) | \n", + "3.0e+00, 1.1e-06, (32) | \n", + "5.0e+01, 9.9e-05, (32) | \n", + "7.5e+01, 4.1e-04, (32) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "