Window size benchmarking

Below is a comparison of 4 machine learning models (LDA, SVC, kNN, ANN) on the 5 datasets with a subject-independent approach (testing with unseen subjects), with a range of different durations per trial (2 to 10 seconds) to study the influence of this parameter on the classification accuracy [1].

import datetime
import matplotlib.pyplot as plt
import os
import pandas as pd
import seaborn as sns
import torch

from scipy import stats

from benchnirs.load import load_dataset
from benchnirs.process import process_epochs, extract_features
from benchnirs.learn import machine_learn, deep_learn


ALL_DATA_PATH = '../../data/dataset_'  # path to the datasets
DATASETS = {'herff_2014_nb': ['1-back', '2-back', '3-back'],
            'shin_2018_nb': ['0-back', '2-back', '3-back'],
            'shin_2018_wg': ['baseline', 'word generation'],
            'shin_2016_ma': ['baseline', 'mental arithmetic'],
            'bak_2019_me': ['right', 'left', 'foot']}
WINDOW_SIZES = [1.9, 2.9, 3.9, 4.9, 5.9, 6.9, 7.8, 8.9, 9.9]
CONFIDENCE = 0.05  # stat confidence at 95 %


start_time = datetime.datetime.now()
date = start_time.strftime('%Y_%m_%d_%H%M')
out_folder = f'./results/window_size_{date}'
os.makedirs(out_folder)
print(f'Main output folder: {out_folder}/')

print(f'Number of GPUs: {torch.cuda.device_count()}')

with open(f'{out_folder}/results.csv', 'w') as w:
    w.write('dataset;window_size;model;fold;accuracy;hyperparameters\n')

for dataset in DATASETS.keys():
    print(f'=====\n{dataset}\n=====')
    data_path = f'{ALL_DATA_PATH}{dataset[:-3]}/'
    out_path = f'{out_folder}/{dataset}_'

    # Load and preprocess data
    epochs = load_dataset(dataset, data_path, bandpass=[0.01, 0.5],
                          baseline=(-2, 0), roi_sides=True, tddr=True)
    classes = DATASETS[dataset]
    epochs_lab = epochs[classes]

    dict_window_size = {'Window size': [], 'Chance': [],
                        'LDA': [], 'SVC': [], 'kNN': [], 'ANN': []}
    for ws in WINDOW_SIZES:
        print(f'-----\nWindow size {ws}\n-----')
        nirs, labels, groups = process_epochs(epochs_lab, tmax=ws)
        nirs_features = extract_features(nirs, ['mean', 'std', 'slope'])

        # Run models
        lda, hps_lda, _ = machine_learn(
            'lda', nirs_features, labels, groups,
            output_folder=f'{out_path}{ws}_lda')
        svc, hps_svc, _ = machine_learn(
            'svc', nirs_features, labels, groups,
            output_folder=f'{out_path}{ws}_svc')
        knn, hps_knn, _ = machine_learn(
            'knn', nirs_features, labels, groups,
            output_folder=f'{out_path}{ws}_knn')
        ann, hps_ann, _ = deep_learn(
            'ann', nirs_features, labels, groups,
            output_folder=f'{out_path}{ws}_ann')
        dict_window_size['Chance'] += [1/len(classes) for _ in lda]
        dict_window_size['LDA'] += list(lda)
        dict_window_size['SVC'] += list(svc)
        dict_window_size['kNN'] += list(knn)
        dict_window_size['ANN'] += list(ann)
        dict_window_size['Window size'] += [ws for _ in lda]

        # Write results
        results = {'LDA': [lda, hps_lda], 'SVC': [svc, hps_svc],
                   'kNN': [knn, hps_knn], 'ANN': [ann, hps_ann]}
        with open(f'{out_folder}/results.csv', 'a') as w:
            for model in results.keys():
                for fold, accuracy in enumerate(results[model][0]):
                    w.write(f'{dataset};{ws};{model};{fold+1};{accuracy};'
                            f'"{results[model][1][fold]}"\n')

    df_window_size = pd.DataFrame(dict_window_size)
    df_window_size = df_window_size.melt(
        id_vars=['Window size'], value_vars=list(results.keys())+['Chance'],
        var_name='Model', value_name='Accuracy')
    lp = sns.lineplot(data=df_window_size, y='Accuracy', x='Window size',
                      hue='Model', palette='colorblind')
    if len(classes) == 2:
        lp.set(ylim=(0.45, 0.65))
    elif len(classes) == 3:
        lp.set(ylim=(0.25, 0.6))
    plt.legend(bbox_to_anchor=(1.01, 0.5), loc=6)
    plt.savefig(f'{out_path}summary.png', bbox_inches='tight')
    plt.close()


# Stats
print('Stats...')
with open(f'{out_folder}/stats.md', 'w') as w:
    df = pd.read_csv(f'{out_folder}/results.csv', delimiter=';')
    w.write('## Correlation of model accuracies to window size\n\n')
    w.write('|Dataset|Model|Shapiro p-value|Test|p-value|rho|\n')
    w.write('|:---:|:---:|:---:|:---:|:---:|:---:|\n')
    for dataset in DATASETS.keys():
        for model in results.keys():
            w.write(f'|{dataset}|{model}|')
            sub_df = df[(df['dataset'] == dataset) & (df['model'] == model)]
            accuracies = sub_df['accuracy'].to_numpy()
            window_sizes = sub_df['window_size'].to_numpy()
            # Check normality of the distribution
            _, p_shap = stats.shapiro(accuracies)
            w.write(f'{p_shap}|')
            if p_shap > CONFIDENCE:
                # Pearson's correlation test
                rho, p_pears = stats.pearsonr(accuracies, window_sizes)
                w.write(f'Pearson|{p_pears}|{rho}|\n')
            else:
                # Spearman's correlation test
                rho, p_spear = stats.spearmanr(accuracies, window_sizes)
                w.write(f'Spearman|{p_spear}|{rho}|\n')


end_time = datetime.datetime.now()
elapsed_time = end_time - start_time
print(f'===\nElapsed time: {elapsed_time}')

References

[1]

Benerradi, J., Clos, J., Landowska, A., Valstar, M. F., & Wilson, M. L. (2023). Benchmarking framework for machine learning classification from fNIRS data. Frontiers in Neuroergonomics, 4, 994969.