{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# TP2: First Machine Learning Models\nNow that we have visualised the data using Time-frequency and time-series format and understood the dynamic of sleep a bit more in details, we will try to train and interpret the output of some machine learning models to predict sleep stages.\n\nAt the end of this tutorial you will be able to make your first submission to the contest.\n","metadata":{"execution":{"iopub.status.busy":"2022-10-04T11:40:31.534123Z","iopub.execute_input":"2022-10-04T11:40:31.535570Z","iopub.status.idle":"2022-10-04T11:40:44.345563Z","shell.execute_reply.started":"2022-10-04T11:40:31.535524Z","shell.execute_reply":"2022-10-04T11:40:44.343814Z"}}},{"cell_type":"code","source":"pip install dtreeviz","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:25:25.571754Z","iopub.execute_input":"2022-11-27T16:25:25.572148Z","iopub.status.idle":"2022-11-27T16:25:36.610799Z","shell.execute_reply.started":"2022-11-27T16:25:25.572116Z","shell.execute_reply":"2022-11-27T16:25:36.609385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cd \"/kaggle/input/dreem-automated-sleep-staging/\"","metadata":{"id":"XpSZkaqFDrRN","executionInfo":{"status":"ok","timestamp":1603111214265,"user_tz":-120,"elapsed":11781,"user":{"displayName":"Antoine Guillot","photoUrl":"","userId":"13673594054995143914"}},"outputId":"f5fdea53-1c98-41b4-fe67-480d62a48578","execution":{"iopub.status.busy":"2022-11-27T16:25:36.613573Z","iopub.execute_input":"2022-11-27T16:25:36.614555Z","iopub.status.idle":"2022-11-27T16:25:36.622100Z","shell.execute_reply.started":"2022-11-27T16:25:36.614513Z","shell.execute_reply":"2022-11-27T16:25:36.620680Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nfrom datetime import datetime\nfrom scipy.signal import spectrogram\nfrom sklearn.metrics import balanced_accuracy_score, cohen_kappa_score, confusion_matrix\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.preprocessing import StandardScaler\nimport os\nimport random as rd\nimport pandas as pd","metadata":{"id":"ZONWWYI-hKl6","executionInfo":{"status":"ok","timestamp":1603111222409,"user_tz":-120,"elapsed":3666,"user":{"displayName":"Antoine Guillot","photoUrl":"","userId":"13673594054995143914"}},"execution":{"iopub.status.busy":"2022-11-27T16:25:36.623950Z","iopub.execute_input":"2022-11-27T16:25:36.624374Z","iopub.status.idle":"2022-11-27T16:25:36.639203Z","shell.execute_reply.started":"2022-11-27T16:25:36.624336Z","shell.execute_reply":"2022-11-27T16:25:36.638265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Building spectral and statistical features from EEG data\nIn this simple example, we will use features based on the power located in each frequency band, the relative power in the delta, theta, alpha, sigma and beta range.\nWe also provide examples on how to compute several statistic of the EEG signals (stdev, min, max, ...)\n\nWe will then feed each of these feature to several classifier to assess the performances of a simple model.\n","metadata":{}},{"cell_type":"code","source":"frequency_bands = {\n        \"delta\": [0.5, 4],\n        \"theta\": [4, 8],\n        \"alpha\": [8, 12],\n       \"sigma\": [12, 16],\n       \"beta\": [16, 30]\n    }\n\nstatistics = {\n    \"stdev\":lambda x:np.std(x,1),\n    \"mean\":lambda x:np.mean(np.abs(x),1)\n}\nvariable_lists = list(frequency_bands) + list(statistics)\n\ndef get_relative_spectral_power_for_epoch(data):\n    \"\"\"\n    Compute the relative spectral power for each of the frequency bands defined above\n    \"\"\"\n    EEG_FS = 250\n    psds = []\n    sfreqs,t,psd = spectrogram(data,250,nperseg = 1000,noverlap = 750)\n    psd = np.mean(np.abs(psd),-1)\n    spectral_power_band = {}\n    for name, freqband in frequency_bands.items():\n        spec_power = psd[:,(sfreqs >= freqband[0]) & (sfreqs < freqband[1])]\n        spec_power = np.sum(spec_power, 1)\n        spectral_power_band[name] = spec_power / np.sum(psd,1)\n    return spectral_power_band\n\ndef compute_stats(data):\n    \"\"\"\n    Compute the statistics of a signal\n    \"\"\"\n    EEG_FS = 250\n    result = {k:f(data) for k,f in statistics.items()}\n    return result","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:25:36.642385Z","iopub.execute_input":"2022-11-27T16:25:36.642771Z","iopub.status.idle":"2022-11-27T16:25:36.655214Z","shell.execute_reply.started":"2022-11-27T16:25:36.642736Z","shell.execute_reply":"2022-11-27T16:25:36.654008Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In this example, we only use the data from the first EEG channel (Frontal-occipital). You can modify the script to add more EEG channels or take into account accelerometer data.","metadata":{}},{"cell_type":"code","source":"\n\ndata_for_records = {}\nhypnogram_for_records = {}\nhypnograms = pd.read_csv('targets_train.csv')\nfor record in os.listdir(\"training_records\"):\n    record_number = int(record[-5])\n    x = np.load(f'training_records/{record}')\n    data_for_record = x[:,1:250 * 30 + 1]\n    data_for_records[record] = get_relative_spectral_power_for_epoch(data_for_record)\n    data_for_records[record].update(compute_stats(data_for_record))\n    hypnogram_for_records[record] = list(hypnograms[hypnograms['record'] == record_number]['target'])\n    \n\n\nvariable_list = list(data_for_records[record].keys())","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:25:36.656647Z","iopub.execute_input":"2022-11-27T16:25:36.657064Z","iopub.status.idle":"2022-11-27T16:25:41.920269Z","shell.execute_reply.started":"2022-11-27T16:25:36.657028Z","shell.execute_reply":"2022-11-27T16:25:41.919118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x.shape","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:25:41.921769Z","iopub.execute_input":"2022-11-27T16:25:41.922122Z","iopub.status.idle":"2022-11-27T16:25:41.928627Z","shell.execute_reply.started":"2022-11-27T16:25:41.922089Z","shell.execute_reply":"2022-11-27T16:25:41.927728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We split the data between training and test records. ","metadata":{}},{"cell_type":"code","source":"rd.seed(1234)\nrecords_list = list(data_for_records)\nrd.shuffle(records_list)\ntraining_record,test_records = records_list[:4],records_list[4:]\n\nprint('Training records: ',training_record)\nprint('Test records: ', test_records)","metadata":{"id":"KfzUJqPPboAi","executionInfo":{"status":"ok","timestamp":1603111226683,"user_tz":-120,"elapsed":7935,"user":{"displayName":"Antoine Guillot","photoUrl":"","userId":"13673594054995143914"}},"execution":{"iopub.status.busy":"2022-11-27T16:25:41.929937Z","iopub.execute_input":"2022-11-27T16:25:41.931050Z","iopub.status.idle":"2022-11-27T16:25:41.948213Z","shell.execute_reply.started":"2022-11-27T16:25:41.931014Z","shell.execute_reply":"2022-11-27T16:25:41.946863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And build the data under tabular format, each epoch has seven features (one for each frequency band).","metadata":{"id":"68gD5UbGEDod"}},{"cell_type":"code","source":"def build_dataset(records, data_for_records,hypnogram_for_records):\n    X,y = [],[]\n    for record in records:\n        X_for_record = np.expand_dims(np.array([value for value in data_for_records[record].values()]).T,0)\n        y.extend(hypnogram_for_records[record])\n        X.extend(X_for_record)\n    return np.concatenate(X),y\n\n\nX_train,y_train = build_dataset(training_record,data_for_records,hypnogram_for_records)\nX_test,y_test = build_dataset(test_records,data_for_records,hypnogram_for_records)\n\n","metadata":{"id":"uZcvEkUvDjWX","executionInfo":{"status":"ok","timestamp":1603111226686,"user_tz":-120,"elapsed":7934,"user":{"displayName":"Antoine Guillot","photoUrl":"","userId":"13673594054995143914"}},"execution":{"iopub.status.busy":"2022-11-27T16:25:41.950198Z","iopub.execute_input":"2022-11-27T16:25:41.950996Z","iopub.status.idle":"2022-11-27T16:25:41.961671Z","shell.execute_reply.started":"2022-11-27T16:25:41.950950Z","shell.execute_reply":"2022-11-27T16:25:41.960612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(variable_lists)\nprint(X_train[0])\nprint(y_train[0])","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:25:41.963216Z","iopub.execute_input":"2022-11-27T16:25:41.963589Z","iopub.status.idle":"2022-11-27T16:25:41.990658Z","shell.execute_reply.started":"2022-11-27T16:25:41.963557Z","shell.execute_reply":"2022-11-27T16:25:41.989391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## First model, Logistic regression\nFirst, we try to build a Linear classifier to classify the data using the spectral bands.\n\nhttps://scikit-learn.org/stable/modules/linear_model.html#logistic-regression","metadata":{"id":"Jv-Rft5PEcHk"}},{"cell_type":"code","source":"from sklearn.decomposition import PCA\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.metrics import plot_confusion_matrix, f1_score\n# scale input data and reduce dimension\npipeline = make_pipeline(StandardScaler(),LogisticRegression(max_iter = 200))\npipeline.fit(X_train, y_train)\n# linear classifier\n# training: fit the model to the data\n\n# test it\npredictions = pipeline.predict(X_test)\nscores = {'balanced_accuracy': balanced_accuracy_score(y_test, predictions),\n            'cohen_kappa': cohen_kappa_score(y_test, predictions),\n            'macro_f1': f1_score(y_test, predictions,average ='macro')}\n\nprint(scores)\nplot_confusion_matrix(pipeline, X_test,y_test,display_labels = ['Wake','N1','N2','N3','REM'])","metadata":{"id":"Ekji3UiPEO5_","executionInfo":{"status":"ok","timestamp":1603111257565,"user_tz":-120,"elapsed":38808,"user":{"displayName":"Antoine Guillot","photoUrl":"","userId":"13673594054995143914"}},"outputId":"fa0d1599-d5d7-433d-ab61-0fefd4ab82d1","execution":{"iopub.status.busy":"2022-11-27T16:25:41.995746Z","iopub.execute_input":"2022-11-27T16:25:41.996812Z","iopub.status.idle":"2022-11-27T16:25:42.553667Z","shell.execute_reply.started":"2022-11-27T16:25:41.996764Z","shell.execute_reply":"2022-11-27T16:25:42.552555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" **Question: What do you think of the performances ? How can you explain the confusion between the sleep stages ?**\n","metadata":{}},{"cell_type":"markdown","source":"Let's have a look at the LR weights to see how the variable are used to perform prediction for each of the sleep stages.\n- Delta has a positive weight for N3 prediction and negative for N1, REM and WAKE\n- Alpha is correlated with WAKE as well as large stdev\n- Sigma range is correlated with N2 sleep and spindles.","metadata":{}},{"cell_type":"code","source":"df = pd.DataFrame(pipeline[1].coef_)\ndf.index = ['Wake','N1',\"N2\",\"N3\",\"REM\"]\ndf.columns = variable_list\ndf","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:25:42.556412Z","iopub.execute_input":"2022-11-27T16:25:42.557571Z","iopub.status.idle":"2022-11-27T16:25:42.573775Z","shell.execute_reply.started":"2022-11-27T16:25:42.557531Z","shell.execute_reply":"2022-11-27T16:25:42.572780Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Question: Plot some sleep epochs where the model confused some sleep stages, why did these confusions happen ?**","metadata":{}},{"cell_type":"markdown","source":"## Submitting to the contest\n","metadata":{}},{"cell_type":"code","source":"predictions = []\nfor record in os.listdir(\"test_records\"):\n    record_number = int(record[-5])\n    x = np.load(f'test_records/{record}')\n    raw_data_for_record = x[:,1:250 * 30 + 1]\n    data_for_record = get_relative_spectral_power_for_epoch(raw_data_for_record)\n    data_for_record.update(compute_stats(raw_data_for_record))\n    data_for_record = np.array([value for value in data_for_record.values()]).T\n    preds = pipeline.predict(data_for_record)\n    for i, pred in enumerate(preds):\n        predictions.append({\"identifier\":record_number * 10000 + i,'target':pred})\n\npredictions = pd.DataFrame(predictions)\nprint(predictions)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:25:42.575277Z","iopub.execute_input":"2022-11-27T16:25:42.576148Z","iopub.status.idle":"2022-11-27T16:25:44.973486Z","shell.execute_reply.started":"2022-11-27T16:25:42.576114Z","shell.execute_reply":"2022-11-27T16:25:44.972084Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions.to_csv('/kaggle/working/submission.csv',index = None)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:25:44.974676Z","iopub.execute_input":"2022-11-27T16:25:44.974982Z","iopub.status.idle":"2022-11-27T16:25:44.984531Z","shell.execute_reply.started":"2022-11-27T16:25:44.974954Z","shell.execute_reply":"2022-11-27T16:25:44.983372Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Decision Tree\n\n","metadata":{"id":"7---cUpLE-lN"}},{"cell_type":"markdown","source":"Decision Trees (DTs) are a non-parametric supervised learning method used for classification and regression. The goal is to create a model that predicts the value of a target variable by learning simple decision rules inferred from the data features. A tree can be seen as a piecewise constant approximation.\n\nhttps://scikit-learn.org/stable/modules/tree.html#tree","metadata":{}},{"cell_type":"code","source":"from sklearn import tree\nclf = tree.DecisionTreeClassifier(random_state=1234,max_depth = 4, class_weight = 'balanced')\nclf.fit(X_train[:,:len(frequency_bands)],y_train)\n\n\nfrom dtreeviz.trees import dtreeviz # remember to load the package\n\nviz = dtreeviz(clf, X_test[:,:len(frequency_bands)], np.array(y_test),\n                target_name=\"target\",feature_names = list(frequency_bands), \n           class_names=['Wake','N1','N2','N3','REM'],orientation = 'TD')\nviz\n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:25:44.986039Z","iopub.execute_input":"2022-11-27T16:25:44.986449Z","iopub.status.idle":"2022-11-27T16:25:51.759764Z","shell.execute_reply.started":"2022-11-27T16:25:44.986414Z","shell.execute_reply":"2022-11-27T16:25:51.758588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# test it\npredictions = clf.predict(X_test[:,:len(frequency_bands)])\nscores = {'balanced_accuracy': balanced_accuracy_score(y_test, predictions),\n            'cohen_kappa': cohen_kappa_score(y_test, predictions),\n            'macro_f1': f1_score(y_test, predictions,average ='macro')}\n\nprint(scores)\nplot_confusion_matrix(clf, X_test[:,:len(frequency_bands)],y_test,display_labels = ['Wake','N1','N2','N3','REM'])","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:25:51.761565Z","iopub.execute_input":"2022-11-27T16:25:51.762096Z","iopub.status.idle":"2022-11-27T16:25:52.107222Z","shell.execute_reply.started":"2022-11-27T16:25:51.762061Z","shell.execute_reply":"2022-11-27T16:25:52.106372Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Question: How do you interpret the split made by the Decision Tree ? Which additional features could improve the performance ?**","metadata":{}},{"cell_type":"markdown","source":"## Random forest and ensembling","metadata":{}},{"cell_type":"markdown","source":"https://scikit-learn.org/stable/modules/ensemble.html\n\nHere we are going to use the Random Forest method, try to use other ensemble learning functions of the scikit-learn library.\nGo to the online documentation of these functions to set the parameters\n","metadata":{}},{"cell_type":"code","source":"\n# select a classifier and train it\nfrom sklearn.ensemble import RandomForestClassifier\nclf_rf = RandomForestClassifier(random_state=42)\nprint('training...')\nclf_rf.fit(X_train, y_train)\n\n# test it\npredictions = clf_rf.predict(X_test)\nscores = {'balanced_accuracy': balanced_accuracy_score(y_test, predictions),\n            'cohen_kappa': cohen_kappa_score(y_test, predictions),\n            'macro_f1': f1_score(y_test, predictions,average ='macro')}\n\nprint(scores)\nplot_confusion_matrix(clf_rf, X_test,y_test,display_labels = ['Wake','N1','N2','N3','REM'])\nplt.show()","metadata":{"id":"Fp7B5w_GHhYH","executionInfo":{"status":"ok","timestamp":1603111265032,"user_tz":-120,"elapsed":46272,"user":{"displayName":"Antoine Guillot","photoUrl":"","userId":"13673594054995143914"}},"outputId":"6941096b-9b7f-49d4-be93-afb231ec1a78","execution":{"iopub.status.busy":"2022-11-27T16:25:52.108410Z","iopub.execute_input":"2022-11-27T16:25:52.109383Z","iopub.status.idle":"2022-11-27T16:25:53.263691Z","shell.execute_reply.started":"2022-11-27T16:25:52.109339Z","shell.execute_reply":"2022-11-27T16:25:53.262633Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Random forests provides several ways to assess the feature importance.\n- The first one (Mean Decrease impurity) is based on how much a variable can decrease the impurity in the tree.\n- The second one (Permutation importance) is based on how much the permutation between two features decreases the performance\n\nhttps://scikit-learn.org/stable/auto_examples/ensemble/plot_forest_importances.html","metadata":{}},{"cell_type":"code","source":"mdi_importances = pd.Series(\n    clf_rf[-1].feature_importances_, index=variable_lists\n).sort_values(ascending=True)\nfig, ax = plt.subplots(figsize=(25, 10))\nmdi_importances.plot.barh()\nax.set_title(\"Random Forest Feature Importances (MDI)\")\nax.figure.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:25:53.265412Z","iopub.execute_input":"2022-11-27T16:25:53.265862Z","iopub.status.idle":"2022-11-27T16:25:53.630144Z","shell.execute_reply.started":"2022-11-27T16:25:53.265816Z","shell.execute_reply":"2022-11-27T16:25:53.628717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('permutations...')\nfrom sklearn.inspection import permutation_importance\n# We use the build-in sklearn fontion to compute the permutation importance\nresult = permutation_importance(clf_rf, X_test, y_test, n_repeats=50, random_state=0,scoring = 'f1_macro')\nsorted_idx = result.importances_mean.argsort()\n\n# And plot the importance of each variable\nfig, ax = plt.subplots(figsize=(25, 10))\nax.boxplot(result.importances[sorted_idx].T,\n           vert=False, labels=variable_list)\nax.set_title(\"Permutation Importances (Test set)\")\nplt.xlabel('Decrease in Macro-F1')\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:25:53.631700Z","iopub.execute_input":"2022-11-27T16:25:53.632040Z","iopub.status.idle":"2022-11-27T16:26:15.034988Z","shell.execute_reply.started":"2022-11-27T16:25:53.632008Z","shell.execute_reply":"2022-11-27T16:26:15.033703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Question: What do you think of the importance of each feature ?**","metadata":{}},{"cell_type":"markdown","source":"## To begin","metadata":{}},{"cell_type":"markdown","source":"### 1. Loading and plotting the data from a record\nWe seek to plot EEG from the i-ith channel for a specific sleep epoch and the corresponding sleep stage stage","metadata":{}},{"cell_type":"code","source":"def plot_epoch_from_record(record_name,epoch_index,channel):\n    pass","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:26:15.036777Z","iopub.execute_input":"2022-11-27T16:26:15.037501Z","iopub.status.idle":"2022-11-27T16:26:15.042484Z","shell.execute_reply.started":"2022-11-27T16:26:15.037454Z","shell.execute_reply":"2022-11-27T16:26:15.041489Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2. Computing the Fourier transform of an epoch\n\nUsing the previous function, you can now extract the data of a specific sleep epoch and compute the STFT using scipy.","metadata":{}},{"cell_type":"code","source":"def compute_stft_from_record(record_name,epoch_index,channel):\n    pass","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:26:15.044173Z","iopub.execute_input":"2022-11-27T16:26:15.045080Z","iopub.status.idle":"2022-11-27T16:26:15.054957Z","shell.execute_reply.started":"2022-11-27T16:26:15.045016Z","shell.execute_reply":"2022-11-27T16:26:15.053516Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Additional exercises","metadata":{}},{"cell_type":"markdown","source":"**Exercise 1:** \n\na) By having a look at the litterature, what are the features used the most frequently for EEG classification ? Implement these features to see if you manage to improve performances of the model.\n\nb) We only use the data from a single channel of EEG. Can you use the data from the other EEG channels and acceleromter channels to improve the performances ?\n\n**Exercise 2:** We only used a training and a validation dataset and did not optimize for the model hyperparameters. Which method can be used to select the best model hyperparameters ? Implement it for a Random Forest classifier.\n\n**Excercise 3:** A single type of model was used for each of the model above, ensembling, bagging and boosting often improves the performances of the classifier. Implement an ensemble of different types of model to see if performances improves.","metadata":{}},{"cell_type":"markdown","source":"### Appendix 1: Plotting the model mistakes\nThe previous model does mistakes on some sleep epoch. We seek to plot these specific epochs to understand these errors.","metadata":{}},{"cell_type":"markdown","source":"We install yasa in order to be able to plot the spectrogram for the epoch.","metadata":{}},{"cell_type":"code","source":"pip install yasa","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:26:15.056357Z","iopub.execute_input":"2022-11-27T16:26:15.056695Z","iopub.status.idle":"2022-11-27T16:26:26.618713Z","shell.execute_reply.started":"2022-11-27T16:26:15.056659Z","shell.execute_reply":"2022-11-27T16:26:26.617009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In order to be able to plot the raw data on the test set, we iterate over the test records, compute the hypnogram and store the raw data. This way we can match the epoch in the raw data with the one in the predicted and ground truth hypnogram.","metadata":{}},{"cell_type":"code","source":"raw_data_for_records = {}\nfeatures_for_records = {}\nhypnogram_for_records = {}\nhypnograms = pd.read_csv('targets_train.csv')\npredicted_hypnograms_for_record = {}\nfor record in os.listdir(\"training_records\"):\n    if record in test_records:\n        record_number = int(record[-5])\n        x = np.load(f'training_records/{record}')\n        raw_data_for_records[record] = x[:,1 + 250 *30:250 * 30 + 1 + 250 *30]\n        features_for_record = get_relative_spectral_power_for_epoch(raw_data_for_records[record])\n        features_for_record.update(compute_stats(raw_data_for_records[record]))\n        features_for_record = np.array([value for value in features_for_record.values()]).T\n        hypnogram_for_records[record] = list(hypnograms[hypnograms['record'] == record_number]['target'])\n        predicted_hypnograms_for_record[record] = clf_rf.predict(features_for_record)\n\n\nvariable_list = list(data_for_records[record].keys())","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:34:08.707245Z","iopub.execute_input":"2022-11-27T16:34:08.707647Z","iopub.status.idle":"2022-11-27T16:34:11.102323Z","shell.execute_reply.started":"2022-11-27T16:34:08.707615Z","shell.execute_reply":"2022-11-27T16:34:11.101496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"First let's just plot the predicted and true hypnograms for all the records.\nThe model seems to be struggling on record_0.","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import f1_score\nfor record in predicted_hypnograms_for_record:\n    ground_truth = hypnogram_for_records[record]\n    prediction = predicted_hypnograms_for_record[record]\n    f1 = f1_score(ground_truth,prediction,average = 'macro')\n    plt.figure(figsize = (10,5))\n    plt.title(f\"Prediction and ground truth hypnogram for {record.replace('.npy','')} \\n F1: {f1}\")\n    plt.plot(predicted_hypnograms_for_record[record])\n    plt.plot(hypnogram_for_records[record])\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-27T16:50:19.212409Z","iopub.execute_input":"2022-11-27T16:50:19.212801Z","iopub.status.idle":"2022-11-27T16:50:19.878837Z","shell.execute_reply.started":"2022-11-27T16:50:19.212770Z","shell.execute_reply":"2022-11-27T16:50:19.877714Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"These functions are used to plot a specific epoch of data under various format:\n- Raw EEG signals\n- Spectrogram\n- Power frequency plot\n\nWe also display the predicted and ground truth sleep stage.","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport yasa\nfrom lspopt import spectrogram_lspopt\n\ndef plot_epoch(epoch_index,raw_data_for_record,ground_truth_hypnogram,predicted_hypnogram):\n    \"\"\"\n    Plot the raw signal for a given epoch\n    \"\"\"\n    plt.figure(figsize = (20,10))\n    plt.title(f'EEG signal for record \\n Predicted Sleep Stage {predicted_hypnogram[epoch_index]} \\n True Stage: {ground_truth_hypnogram[epoch_index]}')\n    plt.plot(np.arange(0,250 * 30) / 250,raw_data_for_record[epoch_index])\n    plt.xlim(0,30)\n    plt.ylim(-250,250)\n    plt.show()\n    \n\ndef plot_spectrogram(epoch_index,raw_data_for_record,ground_truth_hypnogram,predicted_hypnogram):\n    \"\"\"\n    Plot the spectrogram for a given sleep epoch\n    \"\"\"\n    start = epoch_index\n    end = epoch_index + 1\n    title = f'Spectrogram for record \\n Predicted Sleep Stage {predicted_hypnogram[epoch_index]} | True Stage: {ground_truth_hypnogram[epoch_index]}'\n    f = yasa.plot_spectrogram(np.clip(raw_data_for_record[epoch_index],-200,200), \n                          250,\n                          cmap='Spectral_r', win_sec = 0.5,trimperc = 2)\n    f.set_figheight(7)\n    f.set_figwidth(14)\n    f.suptitle(title)\n    plt.show()\n\n    \n\ndef get_average_spectrum(eeg):\n    \"\"\"\n    Return the average power in each of the fourier bin for several epochs.\n    \"\"\"\n    EEG_FS = 250\n    psds = []\n    freqs,t,psd = spectrogram_lspopt(np.clip(eeg,-200,200),250,nperseg = 1000)\n    return freqs,psd.mean(1)\n\ndef plot_average_spectrum_for_epoch(epoch_index,raw_data_for_record,ground_truth_hypnogram,predicted_hypnogram):\n    \"\"\"\n    Compute and plot the power spectrum of a list of epochs\n    \"\"\"\n    plt.figure(figsize = (12,5))\n    plt.title(f'Spectral power for record \\n Predicted Sleep Stage {predicted_hypnogram[epoch_index]} \\n True Stage: {ground_truth_hypnogram[epoch_index]}')\n    freqs,psd = get_average_spectrum(raw_data_for_record[epoch_index])\n    psd = psd/np.sum(psd)\n    plt.plot(freqs, psd)\n    plt.xlim(0, 18)\n    plt.yscale('log')\n    plt.ylim(1e-6,1)\n    plt.xlabel('Frequency [Hz]')\n    plt.ylabel('Relative power')\n    plt.show()\n    \n","metadata":{"execution":{"iopub.status.busy":"2022-11-27T17:03:36.033702Z","iopub.execute_input":"2022-11-27T17:03:36.034103Z","iopub.status.idle":"2022-11-27T17:03:36.049730Z","shell.execute_reply.started":"2022-11-27T17:03:36.034069Z","shell.execute_reply":"2022-11-27T17:03:36.048528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"With this function, we can plot a few epochs for which the model was right and some other epoch for which the model is wrong.\n\nFirst, we select the data from record 5 - arbitrary decision - and the list of epoch for which the model is right and wrong.","metadata":{}},{"cell_type":"code","source":"record = 'record_3.npy'\nhypnogram = hypnogram_for_records[record]\npredicted_hypnogram = predicted_hypnograms_for_record[record]\nraw_data = raw_data_for_records[record]\ncorrect_indexes = np.where(hypnogram == predicted_hypnogram)[0]\nwrong_indexes = np.where(hypnogram != predicted_hypnogram)[0]","metadata":{"execution":{"iopub.status.busy":"2022-11-27T17:04:44.579112Z","iopub.execute_input":"2022-11-27T17:04:44.579988Z","iopub.status.idle":"2022-11-27T17:04:44.587169Z","shell.execute_reply.started":"2022-11-27T17:04:44.579951Z","shell.execute_reply":"2022-11-27T17:04:44.585844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We randomly select three epochs for which the model is wright and three epoch for which the model is right.","metadata":{}},{"cell_type":"code","source":"correct_indexes_to_plot = np.random.choice(correct_indexes, 3,replace = False)\nwrong_indexes_to_plot = np.random.choice(wrong_indexes, 3,replace = False)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T17:04:46.077819Z","iopub.execute_input":"2022-11-27T17:04:46.078221Z","iopub.status.idle":"2022-11-27T17:04:46.084343Z","shell.execute_reply.started":"2022-11-27T17:04:46.078185Z","shell.execute_reply":"2022-11-27T17:04:46.083149Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For all the true and wrong epochs we will now plot the data.","metadata":{}},{"cell_type":"code","source":"for epoch in correct_indexes_to_plot:\n    print(\"########\")\n    print(epoch)\n    \n    plot_average_spectrum_for_epoch(epoch,raw_data,hypnogram, predicted_hypnogram)\n\n    plot_spectrogram(epoch,raw_data,hypnogram, predicted_hypnogram)\n    plot_epoch(epoch,raw_data,hypnogram, predicted_hypnogram)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T17:04:46.853425Z","iopub.execute_input":"2022-11-27T17:04:46.853919Z","iopub.status.idle":"2022-11-27T17:04:50.903332Z","shell.execute_reply.started":"2022-11-27T17:04:46.853883Z","shell.execute_reply":"2022-11-27T17:04:50.902157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for epoch in wrong_indexes_to_plot:\n    print(\"########\")\n    print(epoch)\n    \n    plot_average_spectrum_for_epoch(epoch,raw_data,hypnogram, predicted_hypnogram)\n\n    plot_spectrogram(epoch,raw_data,hypnogram, predicted_hypnogram)\n    plot_epoch(epoch,raw_data,hypnogram, predicted_hypnogram)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T17:04:50.904879Z","iopub.execute_input":"2022-11-27T17:04:50.905209Z","iopub.status.idle":"2022-11-27T17:04:53.716066Z","shell.execute_reply.started":"2022-11-27T17:04:50.905173Z","shell.execute_reply":"2022-11-27T17:04:53.714993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"You've reached the end of this second tutorial.\n\nLet's go to the last part about deep learning methods.\nOpen the **Tutorial_Sleep_Staging_C** tutorial","metadata":{}}]}