notebooks/benchmark_feature_selection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmark feature selection\n",
    "\n",
    "This notebook contains the benchmarking of feature selection methods and models. The goal is to compare the performance of different feature selection methods and models on the same dataset. \n",
    "\n",
    "Firstly, we will create 5 feature sets based on general-purpose methods and our proposed multi-stage feature selection (look for its implementation in `multi_stage_feature_selection.ipynb` notebook). The following methods will be used:\n",
    "\n",
    "- **Spearman** which will be done with [`spearmanr`](`https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.spearmanr.html`) from `scipy.stats`.\n",
    "- **Lasso** which will be done with [`Lasso`](https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.RFECV.html) with cross-validation on top to select `alpha` parameter.\n",
    "- **SVM + Sequential Forward Selector** which will be done with [`SequentialForwardSelector`](https://rasbt.github.io/mlxtend/api_subpackages/mlxtend.feature_selection/#sequentialfeatureselector) in the `forward=True` manner. The `estimator` parameter underneath will be [`SVM.LinearSVR`](https://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVR.html#sklearn.svm.LinearSVR).\n",
    "- **SVM + Sequential Backward Selector** which will be done with [`SequentialBackwardSelector`](https://rasbt.github.io/mlxtend/api_subpackages/mlxtend.feature_selection/#sequentialfeatureselector) in the `forward=False` manner. The `estimator` parameter underneath will be [`SVM.LinearSVR`](https://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVR.html#sklearn.svm.LinearSVR).\n",
    "- **AIC / statistical p-value + Sequential Forward Selector** test which will be done with Sequential Forward Selection (SFS) and Akaike Information Criterion (AIC) using [`OLS`](https://www.statsmodels.org/stable/examples/notebooks/generated/ols.html).\n",
    "- **AIC / statistical p-value + Sequential Backward Selector** test which will be done with Sequential Backward Selection (SBS) and Akaike Information Criterion (AIC) using [`OLS`](https://www.statsmodels.org/stable/examples/notebooks/generated/ols.html).\n",
    "\n",
    "**Inputs**: \n",
    "\n",
    "Please specify `DATA_HOME` in the [\"Loading the data models\"](#Loading-the-data-models). The following files are needed:\n",
    "- Processed data `AD_processed_final.csv`, `RA_processed_final.csv`\n",
    "- Selected features from multi-stage feature selection `ad_multi_stage_features_{TOP_K}.txt`, `ra_multi_stage_features_{TOP_K}.txt`\n",
    "\n",
    "**Outputs**:\n",
    "\n",
    "- `ad_selected_features/` and `ra_selected_features/` directories with the selected features from each method.\n",
    "\n",
    "# Table of contents\n",
    "\n",
    " - [Benchmark feature selection](#Benchmark-feature-selection)<br>\n",
    " - [Table of contents](#Table-of-contents)<br>\n",
    " - [Preparation](#Preparation)<br>\n",
    " - [Loading the data models](#Loading-the-data-models)<br>\n",
    " - [Getting `train` and `test` data sets as separate dataframes](#Getting-%60train%60-and-%60test%60-data-sets-as-separate-dataframes)<br>\n",
    " - [Feature selection pipeline](#Feature-selection-pipeline)<br>\n",
    "   - [Get features based on multi-stage feature selection](#Get-features-based-on-multi-stage-feature-selection)<br>\n",
    "   - [Get features based on `spearmanr` test](#Get-features-based-on-%60spearmanr%60-test)<br>\n",
    "   - [Get features based on `Lasso`](#Get-features-based-on-%60Lasso%60)<br>\n",
    "   - [Get features based on `SFS+SVM.LinearSVR`](#Get-features-based-on-%60SFS%2BSVM.LinearSVR%60)<br>\n",
    "   - [Get features based on `SBS+SVM.LinearSVR`](#Get-features-based-on-%60SBS%2BSVM.LinearSVR%60)<br>\n",
    "   - [Get features based on `SFS+AIC`](#Get-features-based-on-%60SFS%2BAIC%60)<br>\n",
    "   - [Get features based on `SBS+AIC`](#Get-features-based-on-%60SBS%2BAIC%60)<br>\n",
    " - [Show the results](#Show-the-results)<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# System\n",
    "import os\n",
    "import sys\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "# Data science\n",
    "import pandas as pd\n",
    "\n",
    "# Adding the modules to the PYTHONPATH\n",
    "# add the path to your repository below\n",
    "REPO_PATH = \"\"\n",
    "sys.path.append(REPO_PATH)\n",
    "\n",
    "from src.utils import *\n",
    "\n",
    "# Reproducibility\n",
    "SEED = 42"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loading the data models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_HOME = \"\"\n",
    "\n",
    "# Loading the data\n",
    "ad_proc_data = pd.read_csv(os.path.join(DATA_HOME, \"AD_processed_final.csv\"))\n",
    "ra_proc_data = pd.read_csv(os.path.join(DATA_HOME, \"RA_processed_final.csv\"))\n",
    "\n",
    "# Setting up the target, patient ID column names\n",
    "ad_target_col, ad_patient_id = \"endpt_lb_easi1_total_score\", \"patient_id\"\n",
    "ra_target_col, ra_patient_id = \"endpt_lb_das28__crp_\", \"patient_id\"\n",
    "\n",
    "# Getting the feature sets for both use cases (dropping features with high missingness)\n",
    "ad_features = ad_proc_data.filter(regex=\"^ft\").dropna(axis=\"columns\").columns.to_list()\n",
    "ra_features = ra_proc_data.filter(regex=\"^ft_\").dropna(axis=\"columns\").columns.to_list()\n",
    "\n",
    "# Drop records where target is NaN\n",
    "ad_proc_data = ad_proc_data.dropna(subset=[ad_target_col]).reset_index(drop=True)\n",
    "ra_proc_data = ra_proc_data.dropna(\n",
    "    subset=[ra_target_col, \"ft_lb_c_reactive_protein__mg_l_\"]\n",
    ").reset_index(drop=True)\n",
    "\n",
    "# Creating block lists for both datasets based on clinical and technical teams\n",
    "ad_block_list = [\n",
    "    \"ft_sl_actarm\",\n",
    "    \"ft_sl_actarmcd\",\n",
    "    \"ft_sl_ageu\",\n",
    "    \"ft_sl_arm\",\n",
    "    \"ft_sl_armcd\",\n",
    "    \"ft_sl_domain\",\n",
    "    \"ft_sl_dthdtc\",\n",
    "    \"ft_sl_ethnic\",\n",
    "    \"ft_sl_invid\",\n",
    "    \"ft_sl_invnam\",\n",
    "    \"ft_sl_rfstdtc\",\n",
    "    \"ft_sl_rfendtc\",\n",
    "    \"ft_sl_rficdtc\",\n",
    "    \"ft_sl_rfpendtc\",\n",
    "    \"ft_sl_rfxendtc\",\n",
    "    \"ft_sl_rfxstdtc\",\n",
    "    \"ft_sl_rowid\",\n",
    "    \"ft_sl_dthfl\",\n",
    "    \"ft_sl_subjid\",\n",
    "    \"ft_sl_studyid\",\n",
    "    \"ft_cc_easi1_total_score\",\n",
    "] + ad_proc_data.filter(regex=\"^ft_cc\").columns.to_list()\n",
    "ra_block_list = (\n",
    "    [\n",
    "        \"ft_sl_actarm\",\n",
    "        \"ft_sl_actarmcd\",\n",
    "        \"ft_sl_ageu\",\n",
    "        \"ft_sl_arm\",\n",
    "        \"ft_sl_armcd\",\n",
    "        \"ft_sl_domain\",\n",
    "        \"ft_sl_dthdtc\",\n",
    "        \"ft_sl_dthfl\",\n",
    "        \"ft_sl_ethnic\",\n",
    "        \"ft_sl_invid\",\n",
    "        \"ft_sl_invnam\",\n",
    "        \"ft_sl_rfendtc\",\n",
    "        \"ft_sl_rficdtc\",\n",
    "        \"ft_sl_rfpendtc\",\n",
    "        \"ft_sl_rfstdtc\",\n",
    "        \"ft_sl_rfxendtc\",\n",
    "        \"ft_sl_rfxstdtc\",\n",
    "        \"ft_sl_studyid\",\n",
    "        \"ft_sl_subjid\",\n",
    "        \"ft_sl_siteid\",\n",
    "    ]\n",
    "    + ra_proc_data.filter(regex=\"ft_eff\").columns.to_list()\n",
    "    + [\"endpt_lb_das28__esr_\"]\n",
    ")\n",
    "\n",
    "# Getting features used for modeling\n",
    "ad_features_to_use = [f for f in ad_features if f not in ad_block_list]\n",
    "ra_features_to_use = [f for f in ra_features if f not in ra_block_list]\n",
    "\n",
    "print(\n",
    "    f\"AD dataset contains {ad_proc_data.shape[0]} samples and {ad_proc_data.shape[1]} columns..\"\n",
    ")\n",
    "print()\n",
    "print(\n",
    "    f\"RA dataset contains {ra_proc_data.shape[0]} samples and {ra_proc_data.shape[1]} columns..\"\n",
    ")\n",
    "print()\n",
    "print(\"AD dataset:\")\n",
    "display(ad_proc_data)\n",
    "print()\n",
    "print(\"RA dataset:\")\n",
    "display(ra_proc_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting `train` and `test` data sets as separate dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getting `X`s\n",
    "ad_X_train, ad_X_test = (\n",
    "    ad_proc_data.loc[ad_proc_data[\"split\"] == \"TRAIN\", ad_features_to_use],\n",
    "    ad_proc_data.loc[ad_proc_data[\"split\"] == \"TEST\", ad_features_to_use],\n",
    ")\n",
    "ra_X_train, ra_X_test = (\n",
    "    ra_proc_data.loc[ra_proc_data[\"split\"] == \"TRAIN\", ra_features_to_use],\n",
    "    ra_proc_data.loc[ra_proc_data[\"split\"] == \"TEST\", ra_features_to_use],\n",
    ")\n",
    "\n",
    "# Getting `groups`s\n",
    "ad_groups_train, ad_groups_test = (\n",
    "    ad_proc_data.loc[ad_proc_data[\"split\"] == \"TRAIN\", ad_patient_id],\n",
    "    ad_proc_data.loc[ad_proc_data[\"split\"] == \"TEST\", ad_patient_id],\n",
    ")\n",
    "ra_groups_train, ra_groups_test = (\n",
    "    ra_proc_data.loc[ra_proc_data[\"split\"] == \"TRAIN\", ra_patient_id],\n",
    "    ra_proc_data.loc[ra_proc_data[\"split\"] == \"TEST\", ra_patient_id],\n",
    ")\n",
    "\n",
    "# Creating `y`s\n",
    "ad_y_train, ad_y_test = (\n",
    "    ad_proc_data.loc[ad_proc_data[\"split\"] == \"TRAIN\", ad_target_col],\n",
    "    ad_proc_data.loc[ad_proc_data[\"split\"] == \"TEST\", ad_target_col],\n",
    ")\n",
    "ra_y_train, ra_y_test = (\n",
    "    ra_proc_data.loc[ra_proc_data[\"split\"] == \"TRAIN\", ra_target_col],\n",
    "    ra_proc_data.loc[ra_proc_data[\"split\"] == \"TEST\", ra_target_col],\n",
    ")\n",
    "\n",
    "print(\n",
    "    f\"AD dataset contains:\\n\\t\"\n",
    "    f\"TRAIN: {ad_X_train.shape[1]} features and {ad_X_train.shape[0]} samples!\\n\\t\"\n",
    "    f\"TEST: {ad_X_test.shape[1]} features and {ad_X_test.shape[0]} samples!\\n\\t\"\n",
    ")\n",
    "print(\n",
    "    f\"RA dataset contains:\\n\\t\"\n",
    "    f\"TRAIN: {ra_X_train.shape[1]} features and {ra_X_train.shape[0]} samples!\\n\\t\"\n",
    "    f\"TEST: {ra_X_test.shape[1]} features and {ra_X_test.shape[0]} samples!\\n\\t\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Feature selection pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get features based on multi-stage feature selection\n",
    "\n",
    "Let us first list the features selected by multi-stage feature selection (please look at `multi_stage_feature_selection.ipynb` for the implementation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# AD - read features from files\n",
    "fn = os.path.join(DATA_HOME, \"ad_selected_features\", \"ad_multi_stage_features_0.05.txt\")\n",
    "ad_multi_features_0_05 = read_list(fn)\n",
    "\n",
    "# RA - read features from files\n",
    "fn = os.path.join(DATA_HOME, \"ra_selected_features\", \"ra_multi_stage_features_0.05.txt\")\n",
    "ra_multi_features_0_05 = read_list(fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get features based on `spearmanr` test\n",
    "\n",
    "In this section we will select features based on Spearman test using [`spearmanr`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.spearmanr.html). We will take `pvalue=0.05` significance level to select features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import spearmanr\n",
    "\n",
    "significance_level = 0.05\n",
    "statistic_cutoff = 0.2\n",
    "\n",
    "# AD dataset\n",
    "print(\"AD dataset:\")\n",
    "ad_spearman_features = ad_X_train.columns[\n",
    "    [\n",
    "        spearmanr(ad_X_train[f], ad_y_train).pvalue <= significance_level\n",
    "        and spearmanr(ad_X_train[f], ad_y_train).statistic >= statistic_cutoff\n",
    "        for f in ad_X_train.columns\n",
    "    ]\n",
    "].to_list()\n",
    "print(\n",
    "    f\"The following features were selected ({len(ad_spearman_features)}): {ad_spearman_features}\"\n",
    ")\n",
    "print()\n",
    "\n",
    "# Save the features to a file in DATA_HOME\n",
    "fn = os.path.join(DATA_HOME, \"ad_selected_features\", \"ad_spearman_features.txt\")\n",
    "write_list(ad_spearman_features, fn)\n",
    "\n",
    "\n",
    "# RA dataset\n",
    "print(\"RA dataset:\")\n",
    "ra_spearman_features = ra_X_train.columns[\n",
    "    [\n",
    "        spearmanr(ra_X_train[f], ra_y_train).pvalue < significance_level\n",
    "        and spearmanr(ra_X_train[f], ra_y_train).statistic >= statistic_cutoff\n",
    "        for f in ra_X_train.columns\n",
    "    ]\n",
    "].to_list()\n",
    "print(\n",
    "    f\"The following features were selected ({len(ra_spearman_features)}): {ra_spearman_features}\"\n",
    ")\n",
    "# Save the features to a file in DATA_HOME\n",
    "fn = os.path.join(DATA_HOME, \"ra_selected_features\", \"ra_spearman_features.txt\")\n",
    "write_list(ra_spearman_features, fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get features based on `Lasso`\n",
    "\n",
    "In this section we will select our features based on [`LassoCV`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LassoCV.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LassoCV\n",
    "\n",
    "cv = 5\n",
    "\n",
    "# AD dataset\n",
    "print(\"AD dataset:\")\n",
    "ad_lasso = LassoCV(cv=cv, random_state=SEED, n_jobs=-1).fit(ad_X_train, ad_y_train)\n",
    "ad_lasso_features = ad_X_train.columns[ad_lasso.coef_ != 0].to_list()\n",
    "print(\n",
    "    f\"The following features were selected ({len(ad_lasso_features)}): {ad_lasso_features}\"\n",
    ")\n",
    "print()\n",
    "\n",
    "# Save the features to a file in DATA_HOME\n",
    "fn = os.path.join(DATA_HOME, \"ad_selected_features\", \"ad_lasso_features.txt\")\n",
    "write_list(ad_lasso_features, fn)\n",
    "\n",
    "# RA dataset\n",
    "print(\"RA dataset:\")\n",
    "ra_lasso = LassoCV(cv=cv, random_state=SEED).fit(ra_X_train, ra_y_train)\n",
    "ra_lasso_features = ra_X_train.columns[ra_lasso.coef_ != 0].to_list()\n",
    "print(\n",
    "    f\"The following features were selected ({len(ra_lasso_features)}): {ra_lasso_features}\"\n",
    ")\n",
    "# Save the features to a file in DATA_HOME\n",
    "fn = os.path.join(DATA_HOME, \"ra_selected_features\", \"ra_lasso_features.txt\")\n",
    "write_list(ra_lasso_features, fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get features based on `SFS+SVM.LinearSVR`\n",
    "\n",
    "In this section we will select features based on [`SequentialForwardSelector`](https://rasbt.github.io/mlxtend/api_subpackages/mlxtend.feature_selection/#sequentialfeatureselector) in the `forward=True` manner. The `estimator` parameter underneath will be [`SVM.LinearSVR`](https://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVR.html#sklearn.svm.LinearSVR). We are using [`mlxtend`](https://github.com/rasbt/mlxtend/tree/master) implementation instead of [`scikit-learn`](https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.SequentialFeatureSelector.html#sklearn.feature_selection.SequentialFeatureSelector) implementation because we can provide `k_features='best'` to return feature subset with the best cross-validation performance and have a really data-driven way to extract features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from mlxtend.feature_selection import SequentialFeatureSelector\n",
    "from sklearn.svm import LinearSVR\n",
    "\n",
    "cv = 5\n",
    "\n",
    "# AD dataset\n",
    "print(\"AD dataset:\")\n",
    "estimator = LinearSVR(dual=False, loss=\"squared_epsilon_insensitive\", random_state=SEED)\n",
    "ad_sfs = SequentialFeatureSelector(\n",
    "    estimator=estimator, cv=cv, k_features=(5, 100), n_jobs=-1, forward=True, verbose=1\n",
    ").fit(ad_X_train, ad_y_train, groups=ad_groups_train)\n",
    "ad_sfs_features = ad_X_train.iloc[:, list(ad_sfs.k_feature_idx_)].columns.to_list()\n",
    "print(\n",
    "    f\"The following features were selected ({len(ad_sfs_features)}): {ad_sfs_features}\"\n",
    ")\n",
    "print()\n",
    "\n",
    "# Save the features to a file in DATA_HOME\n",
    "fn = os.path.join(DATA_HOME, \"ad_selected_features\", \"ad_sfs_features.txt\")\n",
    "write_list(ad_sfs_features, fn)\n",
    "\n",
    "# RA dataset\n",
    "print(\"RA dataset:\")\n",
    "estimator = LinearSVR(dual=False, loss=\"squared_epsilon_insensitive\", random_state=SEED)\n",
    "ra_sfs = SequentialFeatureSelector(\n",
    "    estimator=estimator, cv=cv, k_features=(5, 100), n_jobs=-1, forward=True, verbose=1\n",
    ").fit(ra_X_train, ra_y_train, groups=ra_groups_train)\n",
    "ra_sfs_features = ra_X_train.iloc[:, list(ra_sfs.k_feature_idx_)].columns.to_list()\n",
    "print(\n",
    "    f\"The following features were selected ({len(ra_sfs_features)}): {ra_sfs_features}\"\n",
    ")\n",
    "# Save the features to a file in DATA_HOME\n",
    "fn = os.path.join(DATA_HOME, \"ra_selected_features\", \"ra_sfs_features.txt\")\n",
    "write_list(ra_sfs_features, fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get features based on `SBS+SVM.LinearSVR`\n",
    "\n",
    "In this section we will select features based on [`SequentialBackwardSelector`](https://rasbt.github.io/mlxtend/api_subpackages/mlxtend.feature_selection/#sequentialfeatureselector) in the `forward=False` manner. The `estimator` parameter underneath will be [`SVM.LinearSVR`](https://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVR.html#sklearn.svm.LinearSVR). We are using [`mlxtend`](https://github.com/rasbt/mlxtend/tree/master) implementation instead of [`scikit-learn`](https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.SequentialFeatureSelector.html#sklearn.feature_selection.SequentialFeatureSelector) implementation because we can provide `k_features='best'` to return feature subset with the best cross-validation performance and have a really data-driven way to extract features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mlxtend.feature_selection import SequentialFeatureSelector\n",
    "from sklearn.svm import LinearSVR\n",
    "\n",
    "cv = 5\n",
    "\n",
    "# AD dataset\n",
    "print(\"AD dataset:\")\n",
    "estimator = LinearSVR(dual=False, loss=\"squared_epsilon_insensitive\", random_state=SEED)\n",
    "ad_sbs = SequentialFeatureSelector(\n",
    "    estimator=estimator, cv=cv, k_features=(5, 100), n_jobs=-1, forward=False, verbose=5\n",
    ").fit(ad_X_train, ad_y_train, groups=ad_groups_train)\n",
    "ad_sbs_features = ad_X_train.iloc[:, list(ad_sbs.k_feature_idx_)].columns.to_list()\n",
    "print(\n",
    "    f\"The following features were selected ({len(ad_sbs_features)}): {ad_sbs_features}\"\n",
    ")\n",
    "print()\n",
    "\n",
    "# Save the features to a file in DATA_HOME\n",
    "fn = os.path.join(DATA_HOME, \"ad_selected_features\", \"ad_sbs_features.txt\")\n",
    "write_list(ad_sbs_features, fn)\n",
    "\n",
    "# RA dataset\n",
    "print(\"RA dataset:\")\n",
    "estimator = LinearSVR(dual=False, loss=\"squared_epsilon_insensitive\", random_state=SEED)\n",
    "ra_sbs = SequentialFeatureSelector(\n",
    "    estimator=estimator, cv=cv, k_features=(5, 100), n_jobs=-1, forward=False, verbose=5\n",
    ").fit(ra_X_train, ra_y_train, groups=ra_groups_train)\n",
    "ra_sbs_features = ra_X_train.iloc[:, list(ra_sbs.k_feature_idx_)].columns.to_list()\n",
    "print(\n",
    "    f\"The following features were selected ({len(ra_sbs_features)}): {ra_sbs_features}\"\n",
    ")\n",
    "# Save the features to a file in DATA_HOME\n",
    "fn = os.path.join(DATA_HOME, \"ra_selected_features\", \"ra_sbs_features.txt\")\n",
    "write_list(ra_sbs_features, fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get features based on `SFS+AIC`\n",
    "\n",
    "In this section we will select features based on Sequential Forward Selection (SFS) and Akaike Information Criterion (AIC) using [`OLS`](https://www.statsmodels.org/stable/examples/notebooks/generated/ols.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from joblib import Parallel, delayed\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "cv = 5\n",
    "max_features = 100\n",
    "\n",
    "\n",
    "def _calculate_aic(model, data, target):\n",
    "    \"\"\"AIC from here: https://machinelearningmastery.com/probabilistic-model-selection-measures/\"\"\"\n",
    "    n, num_params = data.shape[0], data.shape[1]\n",
    "    mse = mean_squared_error(target, model.predict(data))\n",
    "    aic = n * np.log(mse) + 2 * num_params\n",
    "\n",
    "    return aic\n",
    "\n",
    "\n",
    "def aic_iter(feature, curr_X_train, X_train, y_train, groups, cv):\n",
    "    temp_X_train = curr_X_train.join(X_train[feature])\n",
    "    curr_aic = cross_val_score(\n",
    "        LinearRegression(),\n",
    "        temp_X_train,\n",
    "        y_train,\n",
    "        groups=groups,\n",
    "        cv=cv,\n",
    "        n_jobs=-1,\n",
    "        scoring=_calculate_aic,\n",
    "    ).mean()\n",
    "    return curr_aic, feature\n",
    "\n",
    "\n",
    "def SFS_AIC(X_train, y_train, groups):\n",
    "    \"\"\"A helper function to select features based on AIC using SBS technique. Similar to\n",
    "    https://rasbt.github.io/mlxtend/api_subpackages/mlxtend.feature_selection/#sequentialfeatureselector\n",
    "    with `forward=True`.\"\"\"\n",
    "\n",
    "    aics, added_features = [], []\n",
    "\n",
    "    # Start with no features\n",
    "    curr_X_train = pd.DataFrame(index=X_train.index)\n",
    "\n",
    "    for _ in range(max_features):\n",
    "        best_aic, best_feature = float(\"inf\"), None\n",
    "\n",
    "        # Try adding each feature\n",
    "        results = Parallel(n_jobs=-1)(\n",
    "            delayed(aic_iter)(feature, curr_X_train, X_train, y_train, groups, cv)\n",
    "            for feature in X_train.columns\n",
    "            if feature not in curr_X_train.columns\n",
    "        )\n",
    "\n",
    "        # Extract best feature and best AIC from results according to first element of tuple\n",
    "        best_idx = np.argmin([result[0] for result in results])\n",
    "        best_aic, best_feature = results[best_idx]\n",
    "\n",
    "        # Add best feature\n",
    "        curr_X_train = curr_X_train.join(X_train[best_feature])\n",
    "\n",
    "        # Adding for output\n",
    "        aics.append(best_aic)\n",
    "        added_features.append(best_feature)\n",
    "\n",
    "    # Finding idx of best AIC\n",
    "    best_idx = np.argmin(aics)\n",
    "\n",
    "    # Get selected features\n",
    "    selected_features = added_features[: best_idx + 1]\n",
    "\n",
    "    return selected_features, aics, added_features\n",
    "\n",
    "\n",
    "# AD dataset\n",
    "print(\"AD dataset:\")\n",
    "ad_sfs_aic_features, ad_sfs_aics, ad_sfs_aic_added_features = SFS_AIC(\n",
    "    ad_X_train, ad_y_train, ad_groups_train\n",
    ")\n",
    "print(\n",
    "    f\"The following features were selected ({len(ad_sfs_aic_features)}): {ad_sfs_aic_features}\"\n",
    ")\n",
    "print()\n",
    "\n",
    "# Save the features to a file in DATA_HOME\n",
    "fn = os.path.join(DATA_HOME, \"ad_selected_features\", \"ad_sfs_aic_features.txt\")\n",
    "write_list(ad_sfs_aic_features, fn)\n",
    "\n",
    "# # RA dataset\n",
    "print(\"RA dataset:\")\n",
    "ra_sfs_aic_features, ra_sfs_aics, ra_sfs_aic_added_features = SFS_AIC(\n",
    "    ra_X_train, ra_y_train, ra_groups_train\n",
    ")\n",
    "print(\n",
    "    f\"The following features were selected ({len(ra_sfs_aic_features)}): {ra_sfs_aic_features}\"\n",
    ")\n",
    "print()\n",
    "\n",
    "# Save the features to a file in DATA_HOME\n",
    "fn = os.path.join(DATA_HOME, \"ra_selected_features\", \"ra_sfs_aic_features.txt\")\n",
    "write_list(ra_sfs_aic_features, fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get features based on `SBS+AIC`\n",
    "\n",
    "In this section we will select features based on Sequential Backward Selection (SBS) and Akaike Information Criterion (AIC) using [`OLS`](https://www.statsmodels.org/stable/examples/notebooks/generated/ols.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "cv = 5\n",
    "\n",
    "\n",
    "def _calculate_aic(model, data, target):\n",
    "    \"\"\"AIC from here: https://machinelearningmastery.com/probabilistic-model-selection-measures/\"\"\"\n",
    "    n, num_params = data.shape[0], data.shape[1]\n",
    "    mse = mean_squared_error(target, model.predict(data))\n",
    "    aic = n * np.log(mse) + 2 * num_params\n",
    "\n",
    "    return aic\n",
    "\n",
    "\n",
    "def aic_iter(feature, curr_X_train, y_train, groups, cv):\n",
    "    temp_X_train = curr_X_train.drop(columns=feature)\n",
    "    curr_aic = cross_val_score(\n",
    "        LinearRegression(),\n",
    "        temp_X_train,\n",
    "        y_train,\n",
    "        groups=groups,\n",
    "        cv=cv,\n",
    "        n_jobs=-1,\n",
    "        scoring=_calculate_aic,\n",
    "    ).mean()\n",
    "    return curr_aic, feature\n",
    "\n",
    "\n",
    "def SBS_AIC(X_train, y_train, groups):\n",
    "    \"\"\"A helper function to select features based on AIC using SBS technique. Similar to\n",
    "    https://rasbt.github.io/mlxtend/api_subpackages/mlxtend.feature_selection/#sequentialfeatureselector\n",
    "    with `forward=False`.\"\"\"\n",
    "\n",
    "    aics, removed_features = [], []\n",
    "\n",
    "    # Start with all features\n",
    "    curr_X_train = X_train.copy()\n",
    "\n",
    "    for _ in range(X_train.shape[1] - 1):\n",
    "        best_aic, best_feature = float(\"inf\"), None\n",
    "\n",
    "        # Try removing each feature\n",
    "        results = Parallel(n_jobs=-1)(\n",
    "            delayed(aic_iter)(feature, curr_X_train, y_train, groups, cv)\n",
    "            for feature in curr_X_train.columns\n",
    "        )\n",
    "\n",
    "        # Extract best feature and best AIC from results according to first element of tuple\n",
    "        best_idx = np.argmin([result[0] for result in results])\n",
    "        best_aic, best_feature = results[best_idx]\n",
    "\n",
    "        # Remove best feature\n",
    "        curr_X_train = curr_X_train.drop(columns=best_feature)\n",
    "\n",
    "        # Adding for output\n",
    "        aics.append(best_aic)\n",
    "        removed_features.append(best_feature)\n",
    "\n",
    "    # Finding idx of best AIC\n",
    "    best_idx = np.argmin(aics)\n",
    "\n",
    "    # Get selected features\n",
    "    selected_features = list(\n",
    "        set(X_train.columns) - set(removed_features[: best_idx + 1])\n",
    "    )\n",
    "\n",
    "    return selected_features, aics, removed_features\n",
    "\n",
    "\n",
    "# AD dataset\n",
    "print(\"AD dataset:\")\n",
    "ad_sbs_aic_features, ad_sbs_aics, ad_sbs_aic_added_features = SBS_AIC(\n",
    "    ad_X_train, ad_y_train, ad_groups_train\n",
    ")\n",
    "print(\n",
    "    f\"The following features were selected ({len(ad_sbs_aic_features)}): {ad_sbs_aic_features}\"\n",
    ")\n",
    "print()\n",
    "\n",
    "# Save the features to a file in DATA_HOME\n",
    "fn = os.path.join(DATA_HOME, \"ad_selected_features\", \"ad_sbs_aic_features.txt\")\n",
    "write_list(ad_sbs_aic_features, fn)\n",
    "\n",
    "# # RA dataset\n",
    "print(\"RA dataset:\")\n",
    "ra_sbs_aic_features, ra_sbs_aics, ra_sbs_aic_added_features = SBS_AIC(\n",
    "    ra_X_train, ra_y_train, ra_groups_train\n",
    ")\n",
    "print(\n",
    "    f\"The following features were selected ({len(ra_sbs_aic_features)}): {ra_sbs_aic_features}\"\n",
    ")\n",
    "print()\n",
    "\n",
    "# Save the features to a file in DATA_HOME\n",
    "fn = os.path.join(DATA_HOME, \"ra_selected_features\", \"ra_sbs_aic_features.txt\")\n",
    "write_list(ra_sbs_aic_features, fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Show the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the features and print them\n",
    "ad_selected_features = {\n",
    "    \"multi_stage_0_05\": ad_multi_features_0_05,\n",
    "    \"spearman\": ad_spearman_features,\n",
    "    \"lasso\": ad_lasso_features,\n",
    "    \"sfs\": ad_sfs_features,\n",
    "    \"sbs\": ad_sbs_features,\n",
    "    \"sfs_aic\": ad_sfs_aic_features,\n",
    "    \"sbs_aic\": ad_sbs_aic_features,\n",
    "}\n",
    "ra_selected_features = {\n",
    "    \"multi_stage_0_05\": ra_multi_features_0_05,\n",
    "    \"spearman\": ra_spearman_features,\n",
    "    \"lasso\": ra_lasso_features,\n",
    "    \"sfs\": ra_sfs_features,\n",
    "    \"sbs\": ra_sbs_features,\n",
    "    \"sfs_aic\": ra_sfs_aic_features,\n",
    "    \"sbs_aic\": ra_sbs_aic_features,\n",
    "}\n",
    "\n",
    "print(\"AD:\")\n",
    "print(ad_selected_features)\n",
    "print()\n",
    "print(\"RA:\")\n",
    "print(ra_selected_features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for tag, features in ad_selected_features.items():\n",
    "    print(f\"{tag} ({len(features)}):\")\n",
    "    for pref in [\"ft_sl\", \"ft_lb\", \"ft_mh\", \"ft_cm\"]:\n",
    "        curr_fts = [t.replace(f\"{pref}_\", \"\") for t in features if t.startswith(pref)]\n",
    "        print(f\"{pref} ({len(curr_fts)}): {', '.join(curr_fts)}\")\n",
    "        print()\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for tag, features in ra_selected_features.items():\n",
    "    print(f\"{tag} ({len(features)}):\")\n",
    "    for pref in [\"ft_sl\", \"ft_lb\", \"ft_mh\", \"ft_cm\"]:\n",
    "        curr_fts = [t.replace(f\"{pref}_\", \"\") for t in features if t.startswith(pref)]\n",
    "        print(f\"{pref} ({len(curr_fts)}): {', '.join(curr_fts)}\")\n",
    "        print()\n",
    "    print()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.18"
  },
  "vscode": {
   "interpreter": {
    "hash": "5546cda908818f46075379ccf437c01b71ba37d3306bb47673a7b21af4c78453"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}