<jupyter_start><jupyter_text>Helpful functions<jupyter_code>def test_normality(arr_pd, show_qq_plot=False, normality_test='kolmogorov'):
    """
    Test normality for each group in array
    """
    num_groups = arr_pd.shape[1]
    rejected_normal_groups_l = [False] * num_groups
    normality_results = []
    # Step 1 : Check normality of each interval
    for group_i in range(num_groups):
        sample = arr_pd.iloc[:, group_i].dropna().values

        if show_qq_plot:
            sm.qqplot(sample)
            plt.show()

        if normality_test == 'kolmogorov':
            norm_result = stats.kstest(sample, 'norm')
            norm_pvalue = norm_result.pvalue
        elif normality_test == 'shapiro':
            norm_result = stats.shapiro(sample)
            norm_pvalue = norm_result[1]
        else:
            raise ValueError(f'{normality_test} value error, kolmogorov or shapiro options available')
        # Reject two-sided null hypothesis (identical to norm distribution) for p < 0.05
        if norm_pvalue < 0.05:
            print(f'Interval: {arr_pd.columns[group_i]} != normal')
            rejected_normal_groups_l[group_i] = True
        else:
            print(f'Interval: {arr_pd.columns[group_i]} == normal')

        normality_results.append(norm_pvalue)
    rejected_normal_groups_np = np.array(rejected_normal_groups_l)

    return rejected_normal_groups_np, normality_results

def two_tailed_ttest(arr_pd, col1, col2, is_unpaired=True):
    if is_unpaired:
        return stats.ttest_ind(arr_pd[col1].dropna(), arr_pd[col2].dropna()).pvalue
    else:
        return stats.ttest_rel(arr_pd[col1].dropna(), arr_pd[col2].dropna()).pvalue
    

def get_stats_multiple_groups(arr_pd, show_qq_plot=False, normality_test='kolmogorov', force_parametric=False, within_group_anova=False):
    assert arr_pd.shape[1] > 2

    rejected_normal_groups_np, kstest_result = test_normality(arr_pd, show_qq_plot=show_qq_plot, normality_test=normality_test)

    stat_result = None

    arr_pd_melt = pd.melt(arr_pd, var_name='interval', value_name='value').dropna()
    # Case: Some or all intervals are NOT normal:
    # We may choose to force parametric approach anyways
    if np.sum(rejected_normal_groups_np) > 0 and (not force_parametric):
        print('Non-normal groups -> Kruskal')
        # use kruskal-wallis

        #x = [arr_pd.iloc[:, i].dropna() for i in range(arr_pd.shape[1])]
        #print(stats.kruskal(*x))
        kruskal_result = stats.kruskal(*[group['value'].values for name, group in arr_pd_melt.groupby("interval")])
        print(kruskal_result)

        # On rejection apply post-hoc Dunn w\ bonferroni p adjustment correction
        if kruskal_result.pvalue < 0.05:
            print('Rejection: kruskal')
            post_hoc_result = sp.posthoc_dunn(arr_pd_melt, val_col='value', group_col='interval', p_adjust='bonferroni')
            print('Dunn-bonferroni output', post_hoc_result)
    #Case: All intervals are normal
    else:
        if force_parametric:
            print(f'Force parametric: Anova, num normal groups: {np.sum(~rejected_normal_groups_np)}/{len(rejected_normal_groups_np)}')
        else:
            print('All normal groups -> Anova')
        # Apply anova one-way
        anova_result = f_oneway(*[group['value'].values for name, group in arr_pd_melt.groupby("interval")])
        print(anova_result)

        if anova_result.pvalue < 0.05:
            print('Rejection: anova')
            # On rejection apply tukey's test
            tukey_result = pairwise_tukeyhsd(arr_pd_melt['value'], arr_pd_melt['interval'])

            #display results
            print('Tukey output', tukey_result)
    
base_folder = '/Users/achilleasgeorgiou/Desktop/data_output/'
show_qq_plots =True
show_qq_plots_dense = False<jupyter_output><empty_output><jupyter_text>Figure 1<jupyter_code>plot_folder = 'misc'

data_path = 'capsid_density_data.csv'

full_path = os.path.join(base_folder, plot_folder, data_path)
print(full_path)
arr_pd = pd.read_csv(full_path, index_col=0)
print(arr_pd.head())

rejected_normal_groups_np, normality_results = test_normality(arr_pd, show_qq_plot=show_qq_plots, normality_test='shapiro')

print('Normality p-values', normality_results)
# Quick check that all are normal (not rejected)
assert np.sum(rejected_normal_groups_np) == 0


# Use unpaired-t-test between groups [hour_1, day_10] (which is day_12 is paper), [day_10, cell_body]
hour_day_result = two_tailed_ttest(arr_pd, 'hour_1', 'day_10', is_unpaired=True)
day_cell_result = two_tailed_ttest(arr_pd, 'day_10', 'cell_body', is_unpaired=True)

if hour_day_result < 0.05:
    print(f'hour_1 <-> day_10 ttest hypothesis rejected : {hour_day_result:.3}')

if hour_day_result < 0.05:
    print(f'day_10 <-> cell_body ttest hypothesis rejected: {day_cell_result:.3}')<jupyter_output>/Users/achilleasgeorgiou/Desktop/data_output/misc/capsid_density_data.csv
     hour_1     day_10  cell_body
0  9.400000  13.700000   31.19128
1  7.800000   8.747743   41.00000
2  7.645113  23.882290   29.00000
3  7.506638  21.589300   26.90802
4  6.388628  16.303270   36.84701<jupyter_text>Figure 2  rest to run<jupyter_code>plot_folder = 'dataset_1_all/all/plots'
rest_to_run_data_paths = ["rest_to_run_durations/range_30_61-step_10range_30_61-step_10-data.csv",
                            "rest_to_run_proportions/range_30_61-step_10range_30_61-step_10-data.csv",
                            "rest_to_run_sizes/range_30_61-step_10range_30_61-step_10-data.csv", 
                            "rest_to_run_speed/range_30_61-step_5range_30_61-step_5-data.csv"]

for data_path in rest_to_run_data_paths:
    full_path = os.path.join(base_folder, plot_folder, data_path)
    print(full_path)
    arr_pd = pd.read_csv(full_path, index_col=0)
    
    print(f'For plot: {data_path}')
    get_stats_multiple_groups(arr_pd, show_qq_plot=show_qq_plots_dense)

    print('------------\n\n\n\n')<jupyter_output>/Users/achilleasgeorgiou/Desktop/data_output/dataset_1_all/all/plots/rest_to_run_durations/range_30_61-step_10range_30_61-step_10-data.csv
For plot: rest_to_run_durations/range_30_61-step_10range_30_61-step_10-data.csv
Interval: -30 != normal
Interval: -20 != normal
Interval: -10 != normal
Interval: 0 != normal
Interval: 10 != normal
Interval: 20 != normal
Interval: 30 != normal
Interval: 40 != normal
Interval: 50 != normal
Non-normal groups -> Kruskal
KruskalResult(statistic=36.498526592266735, pvalue=1.4228515672765761e-05)
Rejection: kruskal
Dunn-bonferroni output           -10       -20       -30         0        10   20        30  \
-10  1.000000  1.000000  1.000000  1.000000  1.000000  1.0  0.258258   
-20  1.000000  1.000000  1.000000  1.000000  1.000000  1.0  0.457798   
-30  1.000000  1.000000  1.000000  1.000000  1.000000  1.0  0.154906   
0    1.000000  1.000000  1.000000  1.000000  1.000000  1.0  0.513682   
10   1.000000  1.000000  1.000000  1.000000  1.000000  1.0  0.0538[...]<jupyter_text>run stick run<jupyter_code>plot_folder = 'dataset_1_all/all/plots'
run_stick_run_data_paths = ["run_stick_run_amplitudes/range_20_41-step_10range_20_41-step_10-data.csv",
                            "run_stick_run_durations/range_20_41-step_10range_20_41-step_10-data.csv",
                            "run_stick_run_proportions/range_20_41-step_10range_20_41-step_10-data.csv",
                            "run_stick_run_sizes/range_20_41-step_10range_20_41-step_10-data.csv",
                            "run_stick_run_speed/range_20_41-step_10range_20_41-step_10-data.csv"]

for data_path in run_stick_run_data_paths:
    full_path = os.path.join(base_folder, plot_folder, data_path)
    print(full_path)
    arr_pd = pd.read_csv(full_path, index_col=0)
    
    print(f'For plot: {data_path}')
    get_stats_multiple_groups(arr_pd, show_qq_plot=show_qq_plots_dense)

    print('------------\n\n\n\n')
plot_folder = 'dataset_1_all/all/plots'
run_to_rest_data_paths = ["run_to_rest_amplitudes/range_30_61-step_10range_30_61-step_10-data.csv",
                            "run_to_rest_durations/range_30_61-step_10range_30_61-step_10-data.csv",
                            "run_to_rest_proportions/range_30_61-step_10range_30_61-step_10-data.csv",
                            "run_to_rest_sizes/range_30_61-step_10range_30_61-step_10-data.csv",
                            "run_to_rest_speed/range_30_61-step_5range_30_61-step_5-data.csv"]

for data_path in run_to_rest_data_paths:
    full_path = os.path.join(base_folder, plot_folder, data_path)
    print(full_path)
    arr_pd = pd.read_csv(full_path, index_col=0)
    
    print(f'For plot: {data_path}')
    get_stats_multiple_groups(arr_pd, show_qq_plot=show_qq_plots_dense)

    print('------------\n\n\n\n')<jupyter_output>/Users/achilleasgeorgiou/Desktop/data_output/dataset_1_all/all/plots/run_to_rest_amplitudes/range_30_61-step_10range_30_61-step_10-data.csv
For plot: run_to_rest_amplitudes/range_30_61-step_10range_30_61-step_10-data.csv
Interval: -30 != normal
Interval: -20 != normal
Interval: -10 != normal
Interval: 0 != normal
Interval: 10 != normal
Interval: 20 != normal
Interval: 30 != normal
Interval: 40 != normal
Interval: 50 != normal
Non-normal groups -> Kruskal
KruskalResult(statistic=12.421308627409417, pvalue=0.13337274433150498)
------------




/Users/achilleasgeorgiou/Desktop/data_output/dataset_1_all/all/plots/run_to_rest_durations/range_30_61-step_10range_30_61-step_10-data.csv
For plot: run_to_rest_durations/range_30_61-step_10range_30_61-step_10-data.csv
Interval: -30 != normal
Interval: -20 != normal
Interval: -10 != normal
Interval: 0 != normal
Interval: 10 != normal
Interval: 20 != normal
Interval: 30 != normal
Interval: 40 != normal
Interval: 50 != normal
Non-normal groups -> Kru[...]<jupyter_text>Figure S3  Section C - Astrocyte density in BX layers<jupyter_code>plot_folder = 'misc'

data_path = 'layer_density_astrocyte_data.csv'
full_path = os.path.join(base_folder, plot_folder, data_path)
print(full_path)

arr_pd = pd.read_csv(full_path, index_col=0)
print(arr_pd.head())

get_stats_multiple_groups(arr_pd, show_qq_plot=show_qq_plots, normality_test='shapiro')<jupyter_output>/Users/achilleasgeorgiou/Desktop/data_output/misc/layer_density_astrocyte_data.csv
   astro_layer_1  astro_layer_2  astro_layer_3  astro_layer_4
0            0.0       3.172428      10.812469        5.39728
1            0.0       4.650038      11.912600        2.95308
2            0.0       2.069057      15.193643        4.20275
3            0.0       6.608708      23.680269        3.80695
4            0.0       4.647070      32.186379        4.86940<jupyter_text>Section F - Neuron density in BX layers<jupyter_code>plot_folder = 'misc'
data_path = os.path.join(base_folder, plot_folder, 'layer_density_neuron_data.csv')

full_path = os.path.join(base_folder, plot_folder, data_path)
arr_pd = pd.read_csv(full_path, index_col=0)

print(arr_pd.head())

get_stats_multiple_groups(arr_pd, show_qq_plot=show_qq_plots, normality_test='shapiro')<jupyter_output>axon_layer_1  axon_layer_2  axon_layer_3  axon_layer_4
0      0.000000      0.000000      5.241551      8.727168
1      0.000000      0.683069      5.181616      6.143110
2      0.000000      0.487914      4.472532      4.601446
3     10.381305      2.476657     16.368975     10.039829
4      0.000000      0.921493     14.150535      7.966867<jupyter_text>Fig S5  Section A - rest to run<jupyter_code>col_filt = ['rest', 'running']

plot_folder = 'dataset_1_all/all/plots'

identifiers = ['Amplitude', 'Duration', 'Event number', 'Size']
rest_run_data_paths = ["transition_dots_rest_run/dots_amplitude-data.csv",
                            "transition_dots_rest_run/dots_duration-data.csv",
                            "transition_dots_rest_run/dots_Event number (per minute)-data.csv",
                            "transition_dots_rest_run/dots_size-data.csv"]

for i, data_path in enumerate(rest_run_data_paths):
    full_path = os.path.join(base_folder, plot_folder, data_path)
    print(full_path)
    arr_pd = pd.read_csv(full_path, index_col=0)

    arr_pd = arr_pd[col_filt]
    print(arr_pd.head())

    rejected_normal_groups_np, normality_results = test_normality(arr_pd, show_qq_plot=True, normality_test='shapiro')

    print('Normality p-values', normality_results)

    #Paired two tailed ttest
    ttest_pvalue = two_tailed_ttest(arr_pd, col_filt[0], col_filt[1], is_unpaired=False)

    print('--------------------')
    print(f'Ttest pvalue {col_filt[0]} {col_filt[1]} - {identifiers[i]}: ttest_pvalue {ttest_pvalue} -> {"rejected" if ttest_pvalue < 0.05 else "not-rejected"}')
    print('--------------------')<jupyter_output>/Users/achilleasgeorgiou/Desktop/data_output/dataset_1_all/all/plots/transition_dots_rest_run/dots_amplitude-data.csv
                                 rest   running
m181129_d190111_c001_day_0   1.794823  1.848458
m2000219_d190411_c003_day_0  1.343808  1.635869
m2000219_d190411_c003_day_1  1.259689  1.288583
m2000219_d190411_c004_day_0  1.393615  1.493766
m2000219_d190411_c004_day_1  1.307464  1.397168<jupyter_text>Section B - Rest to vibrisa exploration<jupyter_code>col_filt = ['rest', 'stick_rest']

plot_folder = 'dataset_1_all/all/plots'

identifiers = ['Amplitude', 'Duration', 'Event number', 'Size']
rest_rest_stick_data_paths = ["transition_dots_rest_rest_stick/dots_amplitude-data.csv",
                            "transition_dots_rest_rest_stick/dots_duration-data.csv",
                            "transition_dots_rest_rest_stick/dots_Event number (per minute)-data.csv",
                            "transition_dots_rest_rest_stick/dots_size-data.csv"]

for i, data_path in enumerate(rest_rest_stick_data_paths):
    full_path = os.path.join(base_folder, plot_folder, data_path)
    print(full_path)
    arr_pd = pd.read_csv(full_path, index_col=0)

    arr_pd = arr_pd[col_filt]
    print(arr_pd.head())

    rejected_normal_groups_np, normality_results = test_normality(arr_pd, show_qq_plot=True, normality_test='shapiro')

    print('Normality p-values', normality_results)

    #Paired two tailed ttest
    ttest_pvalue = two_tailed_ttest(arr_pd, col_filt[0], col_filt[1], is_unpaired=False)

    print('--------------------')
    print(f'Ttest pvalue {col_filt[0]} {col_filt[1]} - {identifiers[i]}: ttest_pvalue {ttest_pvalue} -> {"rejected" if ttest_pvalue < 0.05 else "not-rejected"}')
    print('--------------------')<jupyter_output>/Users/achilleasgeorgiou/Desktop/data_output/dataset_1_all/all/plots/transition_dots_rest_rest_stick/dots_amplitude-data.csv
                                 rest  stick_rest
m181129_d190111_c001_day_0   1.794823    1.816852
m2000219_d190411_c003_day_0  1.343808    1.374178
m2000219_d190411_c003_day_1  1.259689    1.212675
m2000219_d190411_c004_day_0  1.393615    1.381925
m2000219_d190411_c004_day_1  1.307464    1.330125<jupyter_text>Section C - Run to vibrisa stimulation<jupyter_code>col_filt = ['running', 'stick_run_ind_15']

plot_folder = 'dataset_1_all/all/plots'

identifiers = ['Amplitude', 'Duration', 'Event number', 'Size']
run_run_stick_data_paths = ["transition_dots_run_run_stick/dots_amplitude-data.csv",
                            "transition_dots_run_run_stick/dots_duration-data.csv",
                            "transition_dots_run_run_stick/dots_Event number (per minute)-data.csv",
                            "transition_dots_run_run_stick/dots_size-data.csv"]

for i, data_path in enumerate(run_run_stick_data_paths):
    full_path = os.path.join(base_folder, plot_folder, data_path)
    print(full_path)
    arr_pd = pd.read_csv(full_path, index_col=0)

    arr_pd = arr_pd[col_filt]
    print(arr_pd.head())

    rejected_normal_groups_np, normality_results = test_normality(arr_pd, show_qq_plot=True, normality_test='shapiro')

    print('Normality p-values', normality_results)

    #Paired two tailed ttest
    ttest_pvalue = two_tailed_ttest(arr_pd, col_filt[0], col_filt[1], is_unpaired=False)

    print('--------------------')
    print(f'Ttest pvalue {col_filt[0]} {col_filt[1]} - {identifiers[i]}: ttest_pvalue {ttest_pvalue} -> {"rejected" if ttest_pvalue < 0.05 else "not-rejected"}')
    print('--------------------')<jupyter_output>/Users/achilleasgeorgiou/Desktop/data_output/dataset_1_all/all/plots/transition_dots_run_run_stick/dots_amplitude-data.csv
                              running  stick_run_ind_15
m181129_d190111_c001_day_0   1.848458          1.891501
m2000219_d190411_c003_day_0  1.635869          1.729852
m2000219_d190411_c003_day_1  1.288583          1.262413
m2000219_d190411_c004_day_0  1.493766          1.503124
m2000219_d190411_c004_day_1  1.397168          1.392630<jupyter_text>Figure S7  (Section B)<jupyter_code>def repeated_measures_anova_S7(arr_pd, show_qq_plot=False, normality_test='kolmogorov'):
    rejected_normal_groups_np, kstest_result = test_normality(arr_pd, show_qq_plot=show_qq_plot, normality_test=normality_test)

    arr_pd["astro_id"] = arr_pd.index + 1
    arr_pd = pd.melt(arr_pd, var_name='interval', value_name='value', id_vars=['astro_id'])
    
    print(pg.rm_anova(data=arr_pd, dv='value', within='interval', subject='astro_id', detailed=True))

    post_hocs = pg.pairwise_ttests(dv='value', within='interval', subject='astro_id', data=arr_pd)
    print('Post hocs \n', post_hocs)

plot_folder = 'single_multi_stick'

data_paths = ['plots/bar_comparison/area-data.csv', 
             'plots/bar_comparison/count-data.csv', 
             'plots/bar_comparison/dffMax2-data.csv', 
             'plots/bar_comparison/time_s-data.csv'
             ]

for i, data_path in enumerate(data_paths):
    full_path = os.path.join(base_folder, plot_folder, data_path)
    print(full_path)
    arr_pd = pd.read_csv(full_path, index_col=0)

    repeated_measures_anova_S7(arr_pd, show_qq_plot=False, normality_test='shapiro')<jupyter_output>/Users/achilleasgeorgiou/Desktop/data_output/single_multi_stick/plots/bar_comparison/area-data.csv
Interval: running == normal
Interval: single == normal
Interval: multi != normal
     Source        SS  DF        MS         F     p-unc     np2       eps
0  interval  0.250781   2  0.125390  8.646306  0.006604  0.6336  0.838032
1     Error  0.145022  10  0.014502       NaN       NaN     NaN       NaN
Post hocs 
    Contrast        A        B  Paired  Parametric         T  dof alternative  \
0  interval    multi  running    True        True  3.463137  5.0   two-sided   
1  interval    multi   single    True        True  3.491265  5.0   two-sided   
2  interval  running   single    True        True  1.369704  5.0   two-sided   

      p-unc   BF10    hedges  
0  0.017982  4.444  0.823452  
1  0.017447  4.545  1.183787  
2  0.229089  0.724  0.449870  
/Users/achilleasgeorgiou/Desktop/data_output/single_multi_stick/plots/bar_comparison/count-data.csv
Interval: running == normal
Interval: sin[...]<jupyter_text>Figure S9  (Section B)<jupyter_code>plot_folder = 'dataset_1_all/all/plots'
#Using the exact .csv file used in figure paper
data_path = 'split_correlation_all/splits-data-old-random.csv'

full_path = os.path.join(base_folder, plot_folder, data_path)
print(full_path)

arr_pd = pd.read_csv(full_path, index_col=0)
arr_pd.head()

get_stats_multiple_groups(arr_pd, show_qq_plot=True, normality_test='shapiro')<jupyter_output>/Users/achilleasgeorgiou/Desktop/data_output/dataset_1_all/all/plots/split_correlation_all/splits-data-old-random.csv
{'fraction_non_alphanumeric': 0.11023459654550141, 'fraction_numerical': 0.08837329208558907, 'mean_word_length': 3.4619277662755925, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '299952', 'license': 'MIT', 'n_tokens_mistral': 8201, 'n_tokens_neox': 7112, 'n_words': 1120, 'path': 'statistical_analysis.ipynb', 'repo_name': 'Achilleas/aqua-py-analysis', 'size': '19395'}