04_validation_nsrr.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NSRR testing set validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import pingouin as pg\n",
    "import sklearn.metrics as skm\n",
    "from tqdm.notebook import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "sns.set(style=\"ticks\", font_scale=1.1)\n",
    "\n",
    "from helper_functions import NUM2STR, STR2NUM\n",
    "from helper_functions import mean_std, median_iqr, perc_transition\n",
    "\n",
    "# Define paths\n",
    "wdir = \"output/cv/\"\n",
    "wdir_demo = \"output/demo/\"\n",
    "outdir = \"output/plots/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Choose model\n",
    "model = \"eeg+eog+emg+demo\"\n",
    "\n",
    "feat_files = glob.glob(wdir + \"%s/cv_loo_nsrr_*.parquet\" % model)\n",
    "\n",
    "df = []\n",
    "for f in feat_files:\n",
    "    df.append(pd.read_parquet(f))\n",
    "    \n",
    "df = pd.concat(df)\n",
    "df['subj'] = df['subj'].astype(str)\n",
    "df['dataset'] = df['dataset'].str.upper()\n",
    "print(df['subj'].nunique(), 'subjects')\n",
    "print(df.shape)\n",
    "df.head().round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demographics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, let's merge with the main demographics\n",
    "df_demo = pd.read_csv(wdir_demo + \"demo_nsrr_all.csv\")\n",
    "df_demo.set_index(['subj', 'dataset'], inplace=True)\n",
    "df_demo.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hours of data\n",
    "df.shape[0] / 120"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Distribution of hours of data per subject\n",
    "(df.groupby('subj')['model'].count() / 120).hist();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Per each night\n",
    "df_scores = []\n",
    "\n",
    "labels = ['N1', 'N2', 'N3', 'R', 'W']\n",
    "\n",
    "for sub in tqdm(df['subj'].unique(), leave=False):\n",
    "    df_sub = df[df['subj'] == sub]\n",
    "    yt = df_sub['y_true']\n",
    "    yp = df_sub['y_pred']\n",
    "    n = yt.shape[0]\n",
    "\n",
    "    sub_scores = {\n",
    "        'dataset': df_sub['dataset'].iloc[0],\n",
    "        # Accuracy\n",
    "        'accuracy': 100 * skm.accuracy_score(yt, yp),\n",
    "        'kappa': 100 * skm.cohen_kappa_score(yt, yp, labels=labels),\n",
    "        'mcc': 100 * skm.matthews_corrcoef(yt, yp),\n",
    "        'f1_macro': 100 * skm.f1_score(yt, yp, labels=labels, average='macro', zero_division=1),\n",
    "        # % Transitions\n",
    "        'dur_hours': (yt.size / 2) / 60,\n",
    "        'perc_trans_true': perc_transition(yt),\n",
    "        'perc_trans_pred': perc_transition(yp),\n",
    "        # Confidence\n",
    "        'avg_confidence': 100 * df_sub['confidence'].mean()\n",
    "    }\n",
    "\n",
    "    # F1 for each stage\n",
    "    f1 = 100 * skm.f1_score(yt, yp, average=None, labels=labels, zero_division=1)\n",
    "    for f, l in zip(f1, labels):\n",
    "        sub_scores['f1_' + l] = f\n",
    "        \n",
    "    # In the case of an error, is the second-highest probability typically the correct one?\n",
    "    error = (yp != yt).to_numpy()\n",
    "    sub_scores['accuracy_second'] = 100 * skm.accuracy_score(yt[error], df_sub['y_pred_second'][error])\n",
    "\n",
    "    # Proportion of each stage (NaN = 0)\n",
    "    prop_true = (yt.value_counts() / n).add_prefix('perc_').add_suffix('_true')\n",
    "    prop_pred = (yp.value_counts() / n).add_prefix('perc_').add_suffix('_pred')\n",
    "    sub_scores.update(prop_true.to_dict())\n",
    "    sub_scores.update(prop_pred.to_dict())\n",
    "\n",
    "    # Append to main dataframe\n",
    "    df_scores.append(pd.DataFrame(sub_scores, index=[sub]))\n",
    "\n",
    "df_scores = pd.concat(df_scores)\n",
    "df_scores.sort_index(axis=1, inplace=True)\n",
    "df_scores.index.name = 'subj'\n",
    "df_scores.set_index('dataset', append=True, inplace=True)\n",
    "\n",
    "# Join with demographics\n",
    "df_scores = df_scores.join(df_demo, how=\"left\")\n",
    "\n",
    "df_scores.round(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fill the NaN in perc_XX by zero\n",
    "# df_scores.isna().sum(0)\n",
    "df_scores.fillna(0, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save to Excel for visual inspection\n",
    "df_scores.round(3).to_excel(wdir + model + \"/df_scores.xlsx\", freeze_panes=(1, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the median\n",
    "df_scores._get_numeric_data().agg(median_iqr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the median per dataset\n",
    "df_scores._get_numeric_data().groupby('dataset').agg(median_iqr).T.round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define colors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = list(sns.color_palette(\"Blues\", n_colors=10, as_cmap=False, desat=1))\n",
    "color_pred = cmap[-1]\n",
    "color_ref = \"tab:orange\"\n",
    "cmap_stages = ['#99d7f1', '#009DDC', 'xkcd:twilight blue', 'xkcd:rich purple', 'xkcd:sunflower']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boxplots\n",
    "\n",
    "#### Accuracy, MCC and F1-macro"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Boxplot of accuracy\n",
    "fig, ax = plt.subplots(1, 1, figsize=(4.5, 4.5), dpi=100)\n",
    "order = df_scores.groupby(\"dataset\")['accuracy'].median().sort_values(ascending=False).index\n",
    "\n",
    "ax1 = sns.boxplot(data=df_scores.reset_index(), y='accuracy', x='dataset', color=color_pred, saturation=1,\n",
    "                  order=order, fliersize=0, width=0.6, notch=True, linewidth=1.5)\n",
    "                  # boxprops=dict(edgecolor=\"k\"))\n",
    "\n",
    "# plt.title(\"Performance on testing set (n=%i)\" % df_scores.shape[0])\n",
    "plt.xlabel(\"Dataset (~100 test nights each)\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.ylim(50, 100)\n",
    "\n",
    "ax.annotate(\"Accuracy (median): %.2f\\nKappa (median): %.2f\" % \n",
    "            (df_scores['accuracy'].median(), df_scores['kappa'].median()),\n",
    "            xy=(5, 55), ha=\"right\", fontstyle=\"italic\", fontweight=\"semibold\")\n",
    "\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_boxplot_accuracy_datasets.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Distribution of F1-score per stage\n",
    "\n",
    "For N3, one issue is that we have a lot of participants with only a few epochs of N3 sleep across the night. If the algorithm misses some or all of these epochs, the resulting F1-score can be very low, or even zero.\n",
    "\n",
    "One solution is to set the F1-scores values to NaN whenever less than 5 min of the specified stage is present in the ground-truth scoring."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_f1 = df_scores[['f1_N1', 'f1_N2', 'f1_N3', 'f1_R', 'f1_W']].copy()\n",
    "df_f1.columns = df_f1.columns.str.split('_').str.get(1)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(4.5, 4.5), dpi=100)\n",
    "sns.boxplot(data=df_f1, palette=cmap_stages, fliersize=0, ax=ax, saturation=1, notch=True)\n",
    "plt.xlabel(\"Stage\")\n",
    "plt.ylabel(\"F1-score\")\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distribution of confidence levels across sleep stages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(4.5, 4.5), dpi=100)\n",
    "sns.boxplot(data=df, x=\"y_true\", y=\"confidence\", palette=cmap_stages, fliersize=0, ax=ax, saturation=1, notch=False)\n",
    "plt.xlabel(\"Stage\")\n",
    "plt.ylabel(\"Confidence\")\n",
    "plt.axhline(0.8, ls=\":\", lw=1, color=\"tab:grey\", zorder=-20)\n",
    "plt.ylim(0, 1.05)\n",
    "sns.despine()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_stage_confidence.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Average confidence per stage\n",
    "df.groupby('y_true')['confidence'].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Percent of epochs flagged as high-confidence per stage\n",
    "df['highconf'] = df['confidence'] >= 0.8\n",
    "(df.groupby('y_true')['highconf'].agg(['sum', 'count', 'mean']).round(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Impact of the EEG reference (mastoids vs Fpz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_scores['ref_mastoids'] = np.isin(df_scores.index.get_level_values(-1), ['SHHS', 'MESA'])\n",
    "display(df_scores.groupby(\"ref_mastoids\")['accuracy'].agg(mean_std))\n",
    "\n",
    "# Regression controlling for age and AHI\n",
    "pg.linear_regression(X=df_scores[[\"ref_mastoids\", \"age\"]].astype(int), y=df_scores['accuracy']).round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Percent transitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Effect size of the percentage of transitions\n",
    "display(df_scores[['perc_trans_pred', 'perc_trans_true']].agg(mean_std))\n",
    "pg.ttest(df_scores['perc_trans_pred'], df_scores['perc_trans_true'], paired=True).round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stage discrepancies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_prop_pred = df_scores.filter(like=\"_pred\").iloc[:, :-1].melt(var_name=\"Stage\", value_name=\"Proportion\", ignore_index=False)\n",
    "df_prop_true = df_scores.filter(like=\"_true\").iloc[:, :-1].melt(var_name=\"Stage\", value_name=\"Proportion\", ignore_index=False)\n",
    "\n",
    "df_prop_pred['Stage'] = df_prop_pred['Stage'].str.split('_').str.get(1)\n",
    "df_prop_true['Stage'] = df_prop_true['Stage'].str.split('_').str.get(1)\n",
    "\n",
    "df_prop_pred['Scoring'] = 'Predicted'\n",
    "df_prop_true['Scoring'] = 'Reference'\n",
    "\n",
    "df_prop = pd.concat((df_prop_pred.reset_index(), df_prop_true.reset_index()))\n",
    "df_prop = df_prop.sort_values(by=['subj', 'Stage', 'Scoring']).reset_index(drop=True)\n",
    "\n",
    "df_prop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the effect size\n",
    "ptest = df_prop.pairwise_ttests(dv=\"Proportion\", within=['Stage', \"Scoring\"], subject=\"subj\", effsize=\"cohen\").iloc[11:, :].round(3)\n",
    "ef = ptest.loc[:, ['Stage', 'cohen']].set_index(\"Stage\").abs()\n",
    "display(ef)\n",
    "\n",
    "# Boxplot\n",
    "fig, ax = plt.subplots(1, 1, figsize=(4.5, 4.5), dpi=100)\n",
    "\n",
    "sns.boxplot(y=df_prop['Proportion'] * 100, x=df_prop['Stage'], hue=df_prop['Scoring'],\n",
    "            hue_order=['Reference', 'Predicted'], \n",
    "            palette=[color_ref, color_pred], \n",
    "            saturation=1, width=0.6, fliersize=0, linewidth=1.5, notch=True)\n",
    "\n",
    "plt.ylim(0, 80)\n",
    "plt.yticks([0, 20, 40, 60, 80])\n",
    "plt.legend(frameon=False, loc=\"upper right\")\n",
    "plt.ylabel(\"Proportion of total recording time (%)\");\n",
    "\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_stage_proportion.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bland-Altman plots\n",
    "\n",
    "See https://sri-human-sleep.github.io/sleep-trackers-performance/AnalyticalPipeline_v1.0.0.html#24_bland-altman_plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 4, figsize=(16, 4), dpi=80, sharex=False, sharey=True)\n",
    "\n",
    "ax = ax.flatten()\n",
    "\n",
    "for i, l in enumerate(labels[:-1]):\n",
    "    ax_i = pg.plot_blandaltman(\n",
    "        x=df_scores['perc_%s_pred' % l] * 100, \n",
    "        y=df_scores['perc_%s_true' % l] * 100,\n",
    "        xaxis=\"y\",\n",
    "        annotate=False,\n",
    "        scatter_kws={'color': color_pred, 'alpha': 0.3, 's': 20},\n",
    "        confidence=None, ax=ax[i])\n",
    "    if l == 'R':\n",
    "        l = \"REM\"\n",
    "    ax_i.set_title(\"%s sleep\" % l)\n",
    "    ax_i.set_xlabel(\"Reference %TIB\")\n",
    "    ax_i.set_ylabel(\"\")\n",
    "    \n",
    "plt.subplots_adjust(wspace=0.2)\n",
    "\n",
    "ax[0].set_ylabel(\"Predicted - reference\\n%TIB\")\n",
    "sns.despine()\n",
    "\n",
    "plt.savefig(outdir + \"cv_baplots.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export for further analyses in R\n",
    "df[['subj', 'y_true', 'y_pred']].to_csv(\"BAplots/predictions.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Confusion matrices\n",
    "\n",
    "The normalized confusion matrices show the sensitivity (= recall)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = 100 * skm.confusion_matrix(df['y_true'], df['y_pred'], labels=labels, normalize=\"true\")\n",
    "cm = pd.DataFrame(cm, index=labels, columns=labels)\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(1, 1, dpi=100, figsize=(4.5, 4.5))\n",
    "sns.heatmap(cm, annot=True, vmin=0, vmax=100, cmap=\"Blues\", square=True, cbar=False, fmt=\".1f\")\n",
    "plt.ylabel(\"Reference\")\n",
    "plt.xlabel(\"Predicted\")\n",
    "# plt.title(\"Recall\", y=1.02, fontweight=\"bold\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_confusion_recall.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Precision \n",
    "# skm.precision_recall_fscore_support(df['y_true'], df['y_pred'], labels=labels)\n",
    "cm = 100 * skm.confusion_matrix(df['y_true'], df['y_pred'], labels=labels, normalize=\"pred\")\n",
    "cm = pd.DataFrame(cm, index=labels, columns=labels).round(1)\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(1, 1, dpi=100, figsize=(4.5, 4.5))\n",
    "sns.heatmap(cm, annot=True, vmin=0, vmax=100, cmap=\"Blues\", square=True, cbar=False, fmt=\".1f\")\n",
    "plt.ylabel(\"Reference\")\n",
    "plt.xlabel(\"Predicted\")\n",
    "# plt.title(\"Precision\", y=1.02, fontweight=\"bold\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_confusion_prec.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2, 3, dpi=90, figsize=(13, 9), sharex=False, sharey=False)\n",
    "\n",
    "ax = ax.flatten()\n",
    "\n",
    "for i, dset in enumerate(tqdm(df['dataset'].unique(), leave=False)):\n",
    "    df_dset = df[df['dataset'] == dset]\n",
    "    cm = 100 * skm.confusion_matrix(df_dset['y_true'], df_dset['y_pred'], labels=labels, normalize=\"true\")\n",
    "    cm = pd.DataFrame(cm, index=labels, columns=labels).round(1)\n",
    "\n",
    "    # Plot\n",
    "    sns.heatmap(cm, annot=True, vmin=0, vmax=100, cmap=\"Blues\", square=False, cbar=False, ax=ax[i])\n",
    "    ax[i].set_title(\"%s\" % dset, fontweight=\"bold\", y=1.02)\n",
    "    ax[i].set_ylabel(\"Reference\")\n",
    "    ax[i].set_xlabel(\"Predicted\")\n",
    "\n",
    "plt.subplots_adjust(hspace=0.4, wspace=0.4)\n",
    "plt.savefig(outdir + \"cv_confusion_datasets.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*************\n",
    "\n",
    "## Moderator analyses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First we define a color palette that we'll use in each individual plot\n",
    "# cmap = list(sns.color_palette(\"magma_r\", n_colors=8, desat=1))\n",
    "# sns.color_palette(\"magma_r\", n_colors=8, desat=1)\n",
    "\n",
    "# From: https://carto.com/carto-colors/\n",
    "# cmap = ['#E58606', '#5D69B1', '#CC61B0', '#24796C', \n",
    "#         '#DAA51B', '#2F8AC4', '#764E9F', '#ED645A', '#CC3A8E', '#A5AA99']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Accuracy as a function of overall confidence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_conf = df_scores.pairwise_corr(['avg_confidence', 'accuracy']).round(3)\n",
    "display(corr_conf)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(4.5, 4.5), dpi=100)\n",
    "sns.regplot(data=df_scores, x=\"avg_confidence\", y=\"accuracy\", truncate=True, order=1,\n",
    "            scatter_kws={\"s\": 20, \"alpha\": .2, \"lw\": 1},\n",
    "            line_kws={\"color\": \"k\", \"lw\": 3}, \n",
    "            color=color_pred, ax=ax)\n",
    "plt.xlim(60, 100)\n",
    "plt.ylim(50, 100)\n",
    "plt.xlabel(\"Average confidence\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "\n",
    "plt.annotate(\"r=%.2f\" % corr_conf.loc[0, 'r'], (0.6, 0.1), xycoords=\"axes fraction\", fontstyle=\"italic\")\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_accuracy_confidence.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Accuracy as a function of age / sex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_scores['Sex'] = df_scores['male'].map({0: \"Female\", 1: \"Male\"})\n",
    "\n",
    "display(df_scores['Sex'].value_counts())\n",
    "\n",
    "# T-test\n",
    "display(df_scores.pairwise_ttests(dv=\"accuracy\", between=\"Sex\").round(3))\n",
    "\n",
    "# Plot\n",
    "sns.set(style=\"ticks\", font_scale=1.1)\n",
    "fig, ax = plt.subplots(1, 1, figsize=(2.5, 3.5), dpi=100)\n",
    "\n",
    "sns.boxplot(data=df_scores, x=\"Sex\", y=\"accuracy\", width=0.7, fliersize=0,\n",
    "            saturation=1, color=color_pred, ax=ax)\n",
    "\n",
    "# sns.barplot(data=df_scores, x=\"Sex\", y=\"accuracy\", ci=\"sd\", capsize=.1, errwidth=1, saturation=1, color=cmap[1], ax=ax)\n",
    "\n",
    "# Add \"n=\" on top\n",
    "for i, ac in enumerate([\"Female\", \"Male\"]):\n",
    "    n = df_scores['Sex'].value_counts()[ac]\n",
    "    plt.annotate(\"n=%i\" % n, (i, 98), ha=\"center\", fontsize=11, fontstyle=\"italic\")\n",
    "\n",
    "plt.ylim(50, 100)\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.xlabel(\"Sex\")\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_accuracy_male.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test for raw bivariate correlations between age and accuracy / F1-macro\n",
    "corr_age = df_scores.pairwise_corr([['age'], ['accuracy', 'MCC', 'f1_macro']]).round(3)\n",
    "display(corr_age)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(3.5, 3.5), dpi=100)\n",
    "sns.regplot(data=df_scores, x=\"age\", y=\"accuracy\", truncate=True, order=1,\n",
    "            scatter_kws={\"s\": 20, \"alpha\": .2, \"lw\": 1},\n",
    "            line_kws={\"color\": \"k\", \"lw\": 3}, \n",
    "            color=color_pred, ax=ax)\n",
    "plt.ylim(50, 100)\n",
    "plt.xlim(0, 100)\n",
    "plt.xlabel(\"Age\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "\n",
    "plt.annotate(\"r=%.2f\" % corr_age.loc[0, 'r'], (0.6, 0.1), xycoords=\"axes fraction\", fontstyle=\"italic\")\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_accuracy_agecorr.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create age categories\n",
    "df_scores['age_cut'] = pd.cut(\n",
    "    df_scores['age'], [0, 10, 20, 60, 75, 100],\n",
    "    labels=['<10', '10-20', '20-60', '60-75', '≥75'],\n",
    "    right=False, include_lowest=True)\n",
    "display(df_scores['age_cut'].value_counts(sort=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ANOVA\n",
    "display(df_scores.anova(dv=\"accuracy\", between=\"age_cut\").round(3))\n",
    "display(df_scores.pairwise_tukey(dv=\"accuracy\", between=\"age_cut\").round(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "fig, ax = plt.subplots(1, 1, figsize=(4, 3.5), dpi=100)\n",
    "\n",
    "sns.boxplot(data=df_scores, x=\"age_cut\", y=\"accuracy\", width=0.7, fliersize=0,\n",
    "            saturation=1, color=color_pred, ax=ax)\n",
    "\n",
    "# Add \"n=\" on top\n",
    "for i, ac in enumerate(df_scores['age_cut'].cat.categories):\n",
    "    n = df_scores['age_cut'].value_counts()[ac]\n",
    "    # yloc = df_scores.groupby(['age_cut'])['accuracy'].quantile(0.99)[ac] + 0.03\n",
    "    plt.annotate(\"n=%i\" % n, (i, 98), ha=\"center\", fontsize=11, fontstyle=\"italic\")\n",
    "\n",
    "plt.ylim(50, 100)\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.xlabel(\"Age\")\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_accuracy_age.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Accuracy as a function of BMI / AHI / Ethnicity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation with BMI\n",
    "corr_bmi = df_scores.pairwise_corr([['bmi'], ['accuracy', 'MCC', 'f1_macro']]).round(3)\n",
    "display(corr_bmi)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(3.5, 3.5), dpi=100)\n",
    "sns.regplot(data=df_scores, x=\"bmi\", y=\"accuracy\", truncate=True,\n",
    "            scatter_kws={\"s\": 20, \"alpha\": .2, \"lw\": 1},\n",
    "            line_kws={\"color\": \"k\", \"lw\": 3}, \n",
    "            color=color_pred, ax=ax)\n",
    "plt.ylim(50, 100)\n",
    "plt.xlim(10, 60)\n",
    "plt.xlabel(\"BMI\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "\n",
    "plt.annotate(\"r=%.2f\" % corr_bmi.loc[0, 'r'], (0.6, 0.1), xycoords=\"axes fraction\", fontstyle=\"italic\")\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_accuracy_bmi.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation with AHI - same result when adjusting for age!\n",
    "corr_ahi = df_scores.pairwise_corr([['ahi'], ['accuracy', 'MCC', 'f1_macro']]).round(3)  # covar=['age']\n",
    "display(corr_ahi)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(3.5, 3.5), dpi=100)\n",
    "sns.regplot(data=df_scores, x=\"ahi\", y=\"accuracy\", truncate=True,\n",
    "            scatter_kws={\"s\": 20, \"alpha\": .2, \"lw\": 1},\n",
    "            line_kws={\"color\": \"k\", \"lw\": 3}, \n",
    "            color=color_pred, ax=ax)\n",
    "plt.ylim(50, 100)\n",
    "plt.xlim(0, 110)\n",
    "plt.xlabel(\"AHI\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "\n",
    "plt.annotate(\"r=%.2f\" % corr_ahi.loc[0, 'r'], (0.6, 0.1), xycoords=\"axes fraction\", fontstyle=\"italic\")\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_accuracy_ahi.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Does AHI correlates with the number of stage transitions?\n",
    "df_scores.pairwise_corr(['perc_trans_true', 'ahi']).round(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Is the overall confidence correlated with AHI?\n",
    "df_scores.pairwise_corr(['avg_confidence', 'ahi']).round(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ethnicity - no sig difference\n",
    "display(df_scores.anova(dv=\"accuracy\", between=\"ethnicity\").round(3))\n",
    "display(df_scores.pairwise_tukey(dv=\"accuracy\", between=\"ethnicity\").round(3))\n",
    "\n",
    "df_scores['Race/Ethnicity'] = df_scores['ethnicity'].replace({\n",
    "    'caucasian': \"White\",\n",
    "    'african': \"Black\",\n",
    "    'hispanic': \"Hispanic\",\n",
    "    'other': \"Other\",\n",
    "})\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(1, 1, figsize=(3.5, 3.5), dpi=100)\n",
    "sns.boxplot(data=df_scores, x=\"Race/Ethnicity\", y=\"accuracy\", width=0.7, fliersize=0,\n",
    "            saturation=1, color=color_pred, ax=ax)\n",
    "plt.ylim(50, 100)\n",
    "\n",
    "# Add \"n=\" on top\n",
    "for i, ac in enumerate(df_scores['Race/Ethnicity'].unique()):\n",
    "    n = df_scores['Race/Ethnicity'].value_counts(dropna=False)[ac]\n",
    "    # yloc = df_scores.groupby(['age_cut'])['accuracy'].quantile(0.99)[ac] + 0.03\n",
    "    plt.annotate(\"n=%i\" % n, (i, 98), ha=\"center\", fontsize=11, fontstyle=\"italic\")\n",
    "\n",
    "plt.ylabel(\"Accuracy\")\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_accuracy_ethnicity.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stage transitions and confidence\n",
    "\n",
    "Here, we use PSG to define the transitions between stages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_trans = []\n",
    "\n",
    "for sub in tqdm(df['subj'].unique(), leave=False):\n",
    "    df_sub = df[df['subj'] == sub]\n",
    "    yt = df_sub['y_true']\n",
    "    yp = df_sub['y_pred']\n",
    "\n",
    "    # Identify stable periods, i.e. the 3 epochs before / after are similar (3 minutes window)\n",
    "    first_ep, last_ep = yt.iloc[0], yt.iloc[-1]\n",
    "    stable = np.logical_and.reduce((\n",
    "        yt.shift(1, fill_value=first_ep) == yt,  # = same as previous one\n",
    "        yt.shift(-1, fill_value=last_ep) == yt, # = same as next one\n",
    "        yt.shift(2, fill_value=first_ep) == yt,\n",
    "        yt.shift(-2, fill_value=last_ep) == yt,\n",
    "        yt.shift(3, fill_value=first_ep) == yt,\n",
    "        yt.shift(-3, fill_value=last_ep) == yt,\n",
    "    ))\n",
    "    \n",
    "    # Confidence of the algorithm\n",
    "    highconf = (df_sub['confidence'] >= 0.8).to_numpy()\n",
    "\n",
    "    # Append to main dict\n",
    "    sub_scores = {\n",
    "        # Stage transition\n",
    "        'p_stable': len(stable[stable]) / n,\n",
    "        'p_trans': len(stable[~stable]) / n,\n",
    "        'acc_stable': skm.accuracy_score(yt[stable], yp[stable]),\n",
    "        'acc_trans': skm.accuracy_score(yt[~stable], yp[~stable]),\n",
    "        # Confidence\n",
    "        'p_highconf': len(highconf[highconf]) / n,\n",
    "        'p_lowconf': len(highconf[~highconf]) / n,\n",
    "        'acc_highconf': skm.accuracy_score(yt[highconf], yp[highconf]),\n",
    "        'acc_lowconf': skm.accuracy_score(yt[~highconf], yp[~highconf]),\n",
    "    }\n",
    "\n",
    "    # Append to main dataframe\n",
    "    df_trans.append(100 * pd.DataFrame(sub_scores, index=[sub]))\n",
    "\n",
    "df_trans = pd.concat(df_trans)\n",
    "df_trans.sort_index(axis=1, inplace=True)\n",
    "df_trans.index.name = 'subj'\n",
    "df_trans.round(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Accuracy x Confidence\n",
    "display(df_trans[['acc_highconf', 'acc_lowconf']].apply(mean_std))\n",
    "pg.ttest(df_trans['acc_highconf'], df_trans['acc_lowconf'], paired=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Accuracy x Stage Transitions\n",
    "display(df_trans[['acc_stable', 'acc_trans']].apply(mean_std))\n",
    "pg.ttest(df_trans['acc_stable'], df_trans['acc_trans'], paired=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "sns.set(style=\"ticks\", font_scale=1.1)\n",
    "fig, ax = plt.subplots(1, 1, figsize=(2.5, 3.5), dpi=100)\n",
    "sns.boxplot(data=df_trans[['acc_stable', 'acc_trans']], width=0.7, fliersize=0,\n",
    "            saturation=1, color=color_pred, ax=ax)\n",
    "\n",
    "# sns.barplot(data=df_trans[['acc_stable', 'acc_trans']], ci=\"sd\", capsize=.1, errwidth=1, saturation=1, color=cmap[6], ax=ax)\n",
    "\n",
    "plt.ylim(50, 101)\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.xlabel(\"Epoch\")\n",
    "plt.xticks([0, 1], ['Stable', 'Transition'])\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_accuracy_transitions.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "******\n",
    "\n",
    "## Contributors of accuracy\n",
    "\n",
    "Do low and high accuracy sessions differ in some aspects?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictors = [\n",
    "    'age', 'male', 'bmi', 'ahi', 'ethnicity', 'depression', 'diabete',\n",
    "    'hypertension', 'insomnia', 'dur_hours', 'perc_trans_true', \n",
    "    'perc_N1_true', 'perc_N2_true', 'perc_N3_true', 'perc_R_true', 'perc_W_true', \n",
    "]\n",
    "\n",
    "X = pd.get_dummies(df_scores[predictors])\n",
    "\n",
    "# pg.linear_regression(X=df_scores[predictors], y=df_scores['accuracy'], relimp=True).round(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_scores.pairwise_corr(columns=[[\"accuracy\"], predictors]).round(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "\n",
    "rf = RandomForestRegressor(n_estimators=100, random_state=42)\n",
    "rf.fit(X, y=df_scores['accuracy'])\n",
    "\n",
    "# Percent feature importance\n",
    "relcont = 100 * pd.DataFrame({\"Importance\": rf.feature_importances_}, index=X.columns)\n",
    "relcont"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# relcont.sort_values(\"Importance\", ascending=False).to_clipboard(excel=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*******\n",
    "\n",
    "## Features importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dic_features = {\n",
    "    \"eog_abspow\": \"EOG: Absolute power\",\n",
    "    \"eeg_petrosian\": \"EEG: Petrosian fractal dimension\",\n",
    "    \"eeg_abspow\": \"EEG: Absolute power\",\n",
    "    \"eeg_beta\": \"EEG: Beta power\",\n",
    "    \"eeg_petrosian_c7min_norm\": \"EEG: Petrosian fractal dimension (7.5 min)\",\n",
    "    \"eeg_fdelta\": \"EEG: Fast delta power\",\n",
    "    \"eog_fdelta_c7min_norm\": \"EOG: Fast delta power (7.5 min)\",\n",
    "    \"eeg_db\": \"EEG: Delta/beta power ratio\",\n",
    "    \"eeg_perm\": \"EEG: Permutation entropy\",\n",
    "    \"eog_petrosian\": \"EOG: Petrosian fractal dimension\",\n",
    "    \"eeg_abspow_c7min_norm\": \"EEG: Absolute power (7.5 min)\",\n",
    "    \"eog_std\": \"EOG: Standard deviation\",\n",
    "    \"eeg_theta_c7min_norm\": \"EEG: Theta power (7.5 min)\",\n",
    "    \"eog_skew_c7min_norm\": \"EOG: Skewness (7.5 min)\",\n",
    "    \"eog_higuchi_c7min_norm\": \"EOG: Higuchi fractal dimension (7.5 min)\",\n",
    "    \"eog_hcomp_c7min_norm\": \"EOG: Hjorth complexity (7.5 min)\",\n",
    "    \"eeg_beta_p2min_norm\": \"EEG: Beta power (2 min)\",\n",
    "    \"emg_std_p2min_norm\": \"EMG: Standard deviation (2 min)\",\n",
    "    \"eeg_higuchi\": \"EEG: Higuchi fractal dimension\",\n",
    "    \"time_norm\": \"Time elapsed from beginning of recording (norm.)\"\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fimp = 'output/classifiers/clf_%s_lgb_gbdt_custom_shap.csv' % model\n",
    "df_fimp = pd.read_csv(fimp).head(20)\n",
    "df_fimp['Features'] = df_fimp['Features'].replace(dic_features)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(10, 7))\n",
    "sns.barplot(data=df_fimp, y=\"Features\", x=\"Importance\", palette=\"magma\", saturation=1)\n",
    "\n",
    "plt.ylabel(\"\")\n",
    "plt.xlabel(\"Importance (SHAP values)\")\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.savefig(outdir + \"cv_fimp.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*********\n",
    "\n",
    "## Single-night probabilities & spectrogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import yasa\n",
    "# from mne.io import read_raw_edf\n",
    "\n",
    "# eeg_file = \"data/cfs-visit5-802348.edf\"\n",
    "# hypno_file = \"data/cfs-visit5-802348-profusion.xml\"\n",
    "\n",
    "# include = ['C4', 'LOC', 'EMG1']\n",
    "\n",
    "# sf = 100\n",
    "# path_to_model = 'output/classifiers/clf_eeg+eog+emg+demo_lgb_gbdt_custom.joblib'\n",
    "\n",
    "# raw = read_raw_edf(eeg_file, preload=False, verbose=0)\n",
    "# raw = read_raw_edf(eeg_file, preload=True, \n",
    "#                     exclude=np.setdiff1d(raw.info['ch_names'], include), \n",
    "#                     verbose=0)\n",
    "\n",
    "# # Resample and low-pass filter\n",
    "# raw.resample(sf, npad=\"auto\")\n",
    "    \n",
    "# # LOAD HYPNOGRAM\n",
    "# hypno, sf_hyp = yasa.load_profusion_hypno(hypno_file)\n",
    "# hypno_up = yasa.hypno_upsample_to_data(hypno, sf_hyp, raw)\n",
    "  \n",
    "# # PREDICT SLEEP STAGES\n",
    "# metadata = dict(age=23, male=False)\n",
    "# sls = yasa.SleepStaging(raw, eeg_name=include[0], eog_name=include[1], \n",
    "#                         emg_name=include[2], metadata=metadata)\n",
    "\n",
    "# pred = sls.predict(path_to_model)\n",
    "# pred = pd.Series(pred).replace({'W': 0, 'N1': 1, 'N2': 2, 'N3': 3, 'R': 4}).to_numpy()\n",
    "# pred_up = yasa.hypno_upsample_to_data(pred, sf_hyp, raw)\n",
    "\n",
    "# # Show the accuracy score\n",
    "# skm.accuracy_score(hypno, pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Spectrogram with predicted stage\n",
    "# eeg_data = np.squeeze(raw.pick_channels([include[0]]).get_data() * 1e6)\n",
    "# fig = yasa.plot_spectrogram(eeg_data, sf, hypno=pred_up)\n",
    "# fig.savefig(outdir + \"cv_spec_pred.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Spectrogram with true stage\n",
    "# fig = yasa.plot_spectrogram(eeg_data, sf, hypno=hypno_up)\n",
    "# fig.savefig(outdir + \"cv_spec_true.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Predict proba\n",
    "\n",
    "# def plot_proba(proba):\n",
    "#     proba.index = proba.index / 120\n",
    "#     ax = proba.plot(kind='area', color=cmap_stages, figsize=(12, 4), alpha=.8,\n",
    "#                     stacked=True, lw=0)\n",
    "#     # Add smoothed confidence\n",
    "#     confidence = proba.max(1).rolling(3, center=True).mean()\n",
    "#     ax.plot(confidence, lw=2, color='k', ls='-', alpha=0.9, label='Confidence')\n",
    "#     ax.set_xlim(0, proba.index[-1])\n",
    "#     ax.set_ylim(0, 1)\n",
    "#     ax.set_ylabel(\"Probability\")\n",
    "#     ax.set_xlabel(\"Time [hrs]\")\n",
    "#     return ax\n",
    "\n",
    "# proba = sls.predict_proba(path_to_model)\n",
    "# ax = plot_proba(proba)\n",
    "# plt.savefig(outdir + \"cv_predict_proba_conf.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*********\n",
    "\n",
    "## Features across the night"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Features across the night\n",
    "# sns.set(style=\"darkgrid\", font_scale=1.2)\n",
    "\n",
    "# # Get features and apply smoothing\n",
    "# features = sls.get_features()\n",
    "# features = features.rolling(5, win_type=\"triang\", center=True, min_periods=1).mean()\n",
    "\n",
    "# features['eeg_db'] = np.log(features['eeg_db'])\n",
    "\n",
    "# fig, axes = plt.subplots(4, 1, figsize=(12, 9), sharex=True, sharey=False)\n",
    "\n",
    "# time = np.arange(features.shape[0]) / 120\n",
    "\n",
    "# plt.subplots_adjust(hspace=0.35)\n",
    "\n",
    "# params = dict(color=\"#009DDC\", lw=2)\n",
    "\n",
    "# feat_to_plot = ['eeg_std', 'eeg_db', 'eeg_nzc', 'eeg_higuchi']\n",
    "\n",
    "# axes[0].plot(time, features[feat_to_plot[0]], label=\"Standard deviation\", **params)\n",
    "# axes[1].plot(time, features[feat_to_plot[1]], label=\"Delta/beta ratio (log)\", **params)\n",
    "# axes[2].plot(time, features[feat_to_plot[2]], label=\"Number of zero-crossings\", **params)\n",
    "# axes[3].plot(time, features[feat_to_plot[3]], label=\"Higuchi fractal dimension\", **params)\n",
    "\n",
    "# palette = ['xkcd:sunflower', 'xkcd:twilight blue', 'xkcd:rich purple']\n",
    "\n",
    "# # Mask WAKE, DEEP and REM sleep\n",
    "# for j, st in enumerate([0, 3, 4]):\n",
    "#     mask = np.vstack([hypno == st] * features.shape[1]).T\n",
    "#     features_mask = features.where(mask)\n",
    "#     for i, ax in enumerate(axes):\n",
    "#         axes[i].plot(time, features_mask[feat_to_plot[i]], color=palette[j], lw=2)  \n",
    "\n",
    "# for i, ax in enumerate(axes):\n",
    "#     ax.set_xticks([])\n",
    "#     ax.set_xlim(0, time[-1])\n",
    "#     if i <= 1:\n",
    "#         loc = \"upper right\"\n",
    "#     else:\n",
    "#         loc = \"lower right\"\n",
    "#     ax.legend(loc=loc, frameon=False, fontsize=\"medium\", handlelength=0, handletextpad=0)\n",
    "#     sns.despine(bottom=False, ax=ax)\n",
    "    \n",
    "# axes[0].set_ylim(0, 75)\n",
    "# axes[1].set_ylim(0, 7.5)\n",
    "# axes[1].set_yticks([0, 2.5, 5, 7.5])\n",
    "# axes[2].set_ylim(0, 600)\n",
    "# axes[2].set_yticks([0, 200, 400, 600])\n",
    "# axes[3].set_ylim(1, 1.9)\n",
    "# axes[3].set_yticks([1, 1.3, 1.6, 1.9])\n",
    "    \n",
    "# axes[-1].set_xticks([0, 1, 2, 3, 4, 5, 6])\n",
    "# axes[-1].set_xlabel(\"Time [hrs]\")\n",
    "\n",
    "# sns.despine(bottom=False, ax=axes[-1])\n",
    "\n",
    "# plt.savefig(outdir + \"features_night.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}