""" Plotting of behavioral metrics during the full task (biased blocks) per lab Guido Meijer 6 May 2020 """ import seaborn as sns import numpy as np from os.path import join import matplotlib.pyplot as plt from scipy import stats import scikit_posthocs as sp from paper_behavior_functions import (figpath, seaborn_style, group_colors, institution_map, FIGURE_WIDTH, FIGURE_HEIGHT, QUERY, fit_psychfunc, dj2pandas, load_csv) import pandas as pd from statsmodels.stats.multitest import multipletests # Initialize seaborn_style() figpath = figpath() pal = group_colors() institution_map, col_names = institution_map() col_names = col_names[:-1] # %% Process data if QUERY is True: # query sessions from paper_behavior_functions import query_sessions_around_criterion from ibl_pipeline import reference, subject, behavior use_sessions, _ = query_sessions_around_criterion(criterion='ephys', days_from_criterion=[2, 0], force_cutoff=True) session_keys = (use_sessions & 'task_protocol LIKE "%biased%"').fetch('KEY') ses = ((use_sessions & 'task_protocol LIKE "%biased%"') * subject.Subject * subject.SubjectLab * reference.Lab * (behavior.TrialSet.Trial & session_keys)) ses = ses.proj('institution_short', 'subject_nickname', 'task_protocol', 'session_uuid', 'trial_stim_contrast_left', 'trial_stim_contrast_right', 'trial_response_choice', 'task_protocol', 'trial_stim_prob_left', 'trial_feedback_type', 'trial_response_time', 'trial_stim_on_time', 'session_end_time').fetch( order_by='institution_short, subject_nickname,session_start_time, trial_id', format='frame').reset_index() behav = dj2pandas(ses) behav['institution_code'] = behav.institution_short.map(institution_map) else: behav = load_csv('Fig4.csv') biased_fits = pd.DataFrame() for i, nickname in enumerate(behav['subject_nickname'].unique()): if np.mod(i+1, 10) == 0: print('Processing data of subject %d of %d' % (i+1, len(behav['subject_nickname'].unique()))) # Get lab lab = behav.loc[behav['subject_nickname'] == nickname, 'institution_code'].unique()[0] # Fit psychometric curve left_fit = fit_psychfunc(behav[(behav['subject_nickname'] == nickname) & (behav['probabilityLeft'] == 80)]) right_fit = fit_psychfunc(behav[(behav['subject_nickname'] == nickname) & (behav['probabilityLeft'] == 20)]) fits = pd.DataFrame(data={'threshold_l': left_fit['threshold'], 'threshold_r': right_fit['threshold'], 'bias_l': left_fit['bias'], 'bias_r': right_fit['bias'], 'lapselow_l': left_fit['lapselow'], 'lapselow_r': right_fit['lapselow'], 'lapsehigh_l': left_fit['lapsehigh'], 'lapsehigh_r': right_fit['lapsehigh'], 'nickname': nickname, 'lab': lab}) biased_fits = biased_fits.append(fits, sort=False) # %% Statistics stats_tests = pd.DataFrame(columns=['variable', 'test_type', 'p_value']) posthoc_tests = {} for i, var in enumerate(['threshold_l', 'threshold_r', 'lapselow_l', 'lapselow_r', 'lapsehigh_l', 'lapsehigh_r', 'bias_l', 'bias_r']): _, normal = stats.normaltest(biased_fits[var]) if normal < 0.05: test_type = 'kruskal' test = stats.kruskal(*[group[var].values for name, group in biased_fits.groupby('lab')]) if test[1] < 0.05: # Proceed to posthocs posthoc = sp.posthoc_dunn(biased_fits, val_col=var, group_col='lab') else: posthoc = np.nan else: test_type = 'anova' test = stats.f_oneway(*[group[var].values for name, group in biased_fits.groupby('lab')]) if test[1] < 0.05: posthoc = sp.posthoc_tukey(biased_fits, val_col=var, group_col='lab') else: posthoc = np.nan posthoc_tests['posthoc_'+str(var)] = posthoc stats_tests.loc[i, 'variable'] = var stats_tests.loc[i, 'test_type'] = test_type stats_tests.loc[i, 'p_value'] = test[1] # Correct for multiple tests stats_tests['p_value'] = multipletests(stats_tests['p_value'], method='fdr_bh')[1] # Test between left/right blocks for i, var in enumerate(['threshold', 'lapselow', 'lapsehigh', 'bias']): stats_tests.loc[stats_tests.shape[0] + 1, 'variable'] = '%s_blocks' % var stats_tests.loc[stats_tests.shape[0], 'test_type'] = 'wilcoxon' _, stats_tests.loc[stats_tests.shape[0], 'p_value'] = stats.wilcoxon( biased_fits['%s_l' % var], biased_fits['%s_r' % var]) print(stats_tests) # Print the results # %% Plot metrics f, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(FIGURE_WIDTH*0.8, FIGURE_HEIGHT)) lab_colors = group_colors() ax1.plot([10, 20], [10, 20], linestyle='dashed', color=[0.6, 0.6, 0.6]) for i, lab in enumerate(biased_fits['lab'].unique()): ax1.errorbar(biased_fits.loc[biased_fits['lab'] == lab, 'threshold_l'].mean(), biased_fits.loc[biased_fits['lab'] == lab, 'threshold_r'].mean(), xerr=biased_fits.loc[biased_fits['lab'] == lab, 'threshold_l'].sem(), yerr=biased_fits.loc[biased_fits['lab'] == lab, 'threshold_l'].sem(), fmt='.', color=lab_colors[i]) ax1.set(xlabel='80:20 block', ylabel='20:80 block', title='Threshold', yticks=ax1.get_xticks(), ylim=ax1.get_xlim()) ax2.plot([0, 0.1], [0, 0.1], linestyle='dashed', color=[0.6, 0.6, 0.6]) for i, lab in enumerate(biased_fits['lab'].unique()): ax2.errorbar(biased_fits.loc[biased_fits['lab'] == lab, 'lapselow_l'].mean(), biased_fits.loc[biased_fits['lab'] == lab, 'lapselow_r'].mean(), xerr=biased_fits.loc[biased_fits['lab'] == lab, 'lapselow_l'].sem(), yerr=biased_fits.loc[biased_fits['lab'] == lab, 'lapselow_r'].sem(), fmt='.', color=lab_colors[i]) ax2.set(xlabel='80:20 block', ylabel='', title='Lapse left', yticks=ax2.get_xticks(), ylim=ax2.get_xlim()) ax3.plot([0, 0.1], [0, 0.1], linestyle='dashed', color=[0.6, 0.6, 0.6]) for i, lab in enumerate(biased_fits['lab'].unique()): ax3.errorbar(biased_fits.loc[biased_fits['lab'] == lab, 'lapsehigh_l'].mean(), biased_fits.loc[biased_fits['lab'] == lab, 'lapsehigh_r'].mean(), xerr=biased_fits.loc[biased_fits['lab'] == lab, 'lapsehigh_l'].sem(), yerr=biased_fits.loc[biased_fits['lab'] == lab, 'lapsehigh_l'].sem(), fmt='.', color=lab_colors[i]) ax3.set(xlabel='80:20 block', ylabel='', title='Lapse right', yticks=ax3.get_xticks(), ylim=ax3.get_xlim()) ax4.plot([-10, 10], [-10, 10], linestyle='dashed', color=[0.6, 0.6, 0.6]) for i, lab in enumerate(biased_fits['lab'].unique()): ax4.errorbar(biased_fits.loc[biased_fits['lab'] == lab, 'bias_l'].mean(), biased_fits.loc[biased_fits['lab'] == lab, 'bias_r'].mean(), xerr=biased_fits.loc[biased_fits['lab'] == lab, 'bias_l'].sem(), yerr=biased_fits.loc[biased_fits['lab'] == lab, 'bias_l'].sem(), fmt='.', color=lab_colors[i]) ax4.set(xlabel='80:20 block', ylabel='', title='Bias', yticks=ax4.get_xticks(), ylim=ax4.get_xlim()) plt.tight_layout(w_pad=-0.1) sns.despine(trim=True) plt.savefig(join(figpath, 'figure4f-i_metrics_per_lab_full.pdf')) plt.savefig(join(figpath, 'figure4f-i_metrics_per_lab_full.png'), dpi=300)