bioresponse_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "Aim of the project is to predict wheter a molecule elicit a biological response or not. The project first appeared as Kaggle [Bioresponse Competition](https://www.kaggle.com/c/bioresponse).\n",
    "It is a Binary Classification problem; the target variable is `Activity`: an experimentally measured biological response to molecules in the dataset. `Activity` is `1` if there is a measured biological response, otherwise it is `0`.\n",
    "Available features are 1776 Molecular Descriptors (columns `D1` - `D1776`). \n",
    "\n",
    "The \"physical meaning\" of the available Molecular Descriptors, as well as the specific Biological Response which was measured, are not made known to competition partecipants. No \"domain knowledge\" is applicable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Library Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from time import time\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.ensemble import GradientBoostingClassifier\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "from sklearn.gaussian_process import GaussianProcessClassifier\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.metrics import roc_auc_score\n",
    "from sklearn.metrics import log_loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploratory Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"input/bioresponse/train.csv\")\n",
    "#df = df.sample(100, random_state=1)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No missing values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.isnull().any().any()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_name = \"Activity\"\n",
    "y = df.loc[:, target_name]\n",
    "X = df.drop(target_name, axis=\"columns\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dataset is balanced:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many input features, so it is pointless to visualize each one.\n",
    "We are informed by [Data Description](https://www.kaggle.com/c/bioresponse/data) that input features are already normalized. Can we visually confirm that information, by looking at Violin Plots of Maximum, Minimum and Mean of Features? What about Variance?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stat_functions = [\"min\", \"mean\", \"max\", \"var\"]\n",
    "\n",
    "stats = X.apply(stat_functions).T\n",
    "stats[\"Descriptor\"] = stats.index\n",
    "stats.reset_index(inplace=True, drop=True)\n",
    "stats.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for fun in stat_functions:\n",
    "    sns.violinplot(data=stats, x=fun)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Minimum is peaked around 0, but some values are much higher.\n",
    "Maximum is peaked around 1, but some values are much lower.\n",
    "Mean distribution is simular to Minimum distribution, but the peak is less pronounced.\n",
    "\n",
    "\n",
    "It is confirmed that Normalization, i.e. \"Min-Max Scaling\", was applied to the Descriptors. But predictive models may still benefit from other preliminary scaling techniques.\n",
    "\n",
    "Variance is around 0 for most of features, but some values are higher, up to 0.3 .\n",
    "\n",
    "Are Minimum, Maximum, Mean, and Variance, of Descriptors equally distributed between active and inactive compounds?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_active = df.query(f\"{target_name} == 1\").drop(target_name, axis=\"columns\")\n",
    "X_inactive = df.query(f\"{target_name} == 0\").drop(target_name, axis=\"columns\")\n",
    "\n",
    "stats_active = X_active.apply(stat_functions).T\n",
    "stats_active[target_name] = \"active\"\n",
    "\n",
    "stats_inactive = X_inactive.apply(stat_functions).T\n",
    "stats_inactive[target_name] = \"inactive\"\n",
    "\n",
    "stats = pd.concat([stats_active, stats_inactive])\n",
    "stats[\"Descriptor\"] = stats.index\n",
    "stats.reset_index(inplace=True, drop=True)\n",
    "stats.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for fun in stat_functions:\n",
    "    sns.violinplot(data=stats, x=fun, y=target_name)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can confirm that Descriptor statistics are equally distributed between active and inactive compounds.\n",
    "\n",
    "\n",
    "Are the Descriptors correlated?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr = X.corr()\n",
    "\n",
    "mask = np.zeros_like(corr)\n",
    "mask[np.triu_indices_from(mask)] = True\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True, figsize=(10, 10))\n",
    "sns.heatmap(corr,\n",
    "            vmin=-1,\n",
    "            vmax=1,\n",
    "            cmap=\"PRGn\",\n",
    "            mask=mask,\n",
    "            square=True,\n",
    "            ax=ax,\n",
    "            );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some features are positively correlated. A lower number of features are negatively correlated. Most features have small correlation. Still, the dataset may benefit from some kind of dimensionality reduction, such as PCA.\n",
    "\n",
    "Try PCA without scaling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = PCA()\n",
    "pcs = pca.fit_transform(X)\n",
    "\n",
    "pc_names = [f\"PC{i}\" for i in range(1, pcs.shape[1] + 1)]\n",
    "\n",
    "pcs = pd.DataFrame(pcs, columns=pc_names)\n",
    "pcs[target_name] = y.astype(\"category\")\n",
    "\n",
    "expl_ratios = pca.explained_variance_ratio_\n",
    "expl_ratios = pd.Series(expl_ratios, index=pc_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pc_x = 1\n",
    "pc_y = 2\n",
    "expl_ratio_x = expl_ratios.at[f\"PC{pc_x}\"]\n",
    "expl_ratio_y = expl_ratios.at[f\"PC{pc_y}\"]\n",
    "\n",
    "sns.lmplot(data=pcs,\n",
    "           x=f\"PC{pc_x}\",\n",
    "           y=f\"PC{pc_y}\",\n",
    "           hue=target_name,\n",
    "           fit_reg=False,\n",
    "           scatter_kws={\"alpha\": 0.1},\n",
    "           )\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.set_xlabel(f\"PC{pc_x}: {expl_ratio_x * 100:0.1f}%\")\n",
    "ax.set_ylabel(f\"PC{pc_y}: {expl_ratio_y * 100:0.1f}%\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Three clusters are visible using PC1 and PC2. They do not separe active from inactive compounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cumsum_expl_ratios = expl_ratios.cumsum()\n",
    "cumsum_expl_ratios.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tresholds_cumsum_expl_ratios = [0.7, 0.8, 0.9]\n",
    "\n",
    "\n",
    "n_pcs_for_tresholds_unscaled = [(cumsum_expl_ratios <= tresh).argmin() + 1\n",
    "                                for tresh in tresholds_cumsum_expl_ratios]\n",
    "\n",
    "for tresh, n_pcs in zip(tresholds_cumsum_expl_ratios, n_pcs_for_tresholds_unscaled):\n",
    "    print(f\"{n_pcs} components are needed to have {tresh * 100}% variance\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, without scaling, 242 components are needed to have 90% variance!\n",
    "\n",
    "Try PCA with scaling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = StandardScaler()\n",
    "X_scaled = scaler.fit_transform(X)\n",
    "\n",
    "\n",
    "pca = PCA()\n",
    "pcs = pca.fit_transform(X_scaled)\n",
    "\n",
    "pc_names = [f\"PC{i}\" for i in range(1, pcs.shape[1] + 1)]\n",
    "\n",
    "pcs = pd.DataFrame(pcs, columns=pc_names)\n",
    "pcs[target_name] = y.astype(\"category\")\n",
    "\n",
    "expl_ratios = pca.explained_variance_ratio_\n",
    "expl_ratios = pd.Series(expl_ratios, index=pc_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pc_x = 1\n",
    "pc_y = 2\n",
    "expl_ratio_x = expl_ratios.at[f\"PC{pc_x}\"]\n",
    "expl_ratio_y = expl_ratios.at[f\"PC{pc_y}\"]\n",
    "\n",
    "sns.lmplot(data=pcs,\n",
    "           x=f\"PC{pc_x}\",\n",
    "           y=f\"PC{pc_y}\",\n",
    "           hue=target_name,\n",
    "           fit_reg=False,\n",
    "           scatter_kws={\"alpha\": 0.1},\n",
    "           )\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.set_xlabel(f\"PC{pc_x}: {expl_ratio_x * 100:0.1f}%\")\n",
    "ax.set_ylabel(f\"PC{pc_y}: {expl_ratio_y * 100:0.1f}%\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two clusters are visible when plotting PC1 vs PC2, but they do not separe active from inactive compounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cumsum_expl_ratios = expl_ratios.cumsum()\n",
    "cumsum_expl_ratios.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_pcs_for_tresholds_scaled = [(cumsum_expl_ratios <= tresh).argmin() + 1\n",
    "                              for tresh in tresholds_cumsum_expl_ratios]\n",
    "\n",
    "for tresh, n_pcs in zip(tresholds_cumsum_expl_ratios, n_pcs_for_tresholds_scaled):\n",
    "    print(f\"{n_pcs} components are needed to have {tresh * 100}% variance\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, with scaling, 350 components will produce a dataset with 90% of total variance! 350 components are less than one third of total number of original features."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Algorithm Selection\n",
    "First, create train and test datasets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X,\n",
    "                                                    y,\n",
    "                                                    test_size=0.2,\n",
    "                                                    stratify=y,\n",
    "                                                    random_state=1,\n",
    "                                                    )\n",
    "\n",
    "X_train.shape, X_test.shape, y_train.shape, y_test.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thanks to option `stratify=y` of `train_test_split` function, both\n",
    "train and test contain a balanced amount of active and inactive compounds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train.mean(), y_test.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's experiment with some models before running GridSearchCV:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_time = time()\n",
    "model = make_pipeline(PCA(n_components=10), GradientBoostingClassifier())\n",
    "#model = GradientBoostingClassifier()\n",
    "model.fit(X_train, y_train)\n",
    "end_time = time()\n",
    "tot_time = end_time - start_time\n",
    "tot_time / 3600"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define base hyperparameters that will be fixed in all model pipelines of each type of model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kwargs_pca_unscaled = dict(n_components=n_pcs_for_tresholds_unscaled[1])\n",
    "kwargs_pca_scaled = dict(n_components=n_pcs_for_tresholds_scaled[1])\n",
    "kwargs_l1 = dict(penalty=\"l1\", solver=\"liblinear\", random_state=1)\n",
    "kwargs_l2 = dict(penalty=\"l2\", solver=\"liblinear\", random_state=1)\n",
    "kwargs_rf = dict(random_state=1)\n",
    "kwargs_gb = dict(random_state=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define hyperparameter grids:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#grid_pca_unscaled = {\"pca__n_components\": n_pcs_for_tresholds_unscaled}\n",
    "#grid_pca_scaled = {\"pca__n_components\": n_pcs_for_tresholds_scaled}\n",
    "grid_l1 = {\"logisticregression__C\": [0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50, 100, 500, 1000]}\n",
    "grid_l2 = grid_l1\n",
    "grid_rf = {\"randomforestclassifier__n_estimators\": [100, 200],\n",
    "           \"randomforestclassifier__max_features\": [\"auto\", \"sqrt\", 0.33],\n",
    "           \"randomforestclassifier__min_samples_leaf\": [1, 2, 3, 5, 10],\n",
    "           }\n",
    "grid_gb = {\"gradientboostingclassifier__n_estimators\": [100, 200],\n",
    "           \"gradientboostingclassifier__learning_rate\": [0.03, 0.05, 0.1, 0.2],\n",
    "           \"gradientboostingclassifier__max_depth\": [1, 2, 3, 5, None],\n",
    "           }\n",
    "\n",
    "hyperparameters = {\"l1\": grid_l1,\n",
    "                   \"l2\": grid_l2,\n",
    "                   \"rf\": grid_rf,\n",
    "                   \"gb\": grid_gb,\n",
    "                   }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define model pipelines: for each model type, different kinds of preprocessing will be tried. The model will be applied on a number of PCS corresponding to the previously defined tresholds, and PCA will be applied on the dataset as-is and on the scaled dataset\n",
    "\n",
    "First, the models without scaling will be fitted:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pipelines_unscaled = {\"l1\": make_pipeline(PCA(**kwargs_pca_unscaled),\n",
    "                                          LogisticRegression(**kwargs_l1)),\n",
    "                      \"l2\": make_pipeline(PCA(**kwargs_pca_unscaled),\n",
    "                                          LogisticRegression(**kwargs_l2)),\n",
    "                      \"rf\": make_pipeline(PCA(**kwargs_pca_unscaled),\n",
    "                                          RandomForestClassifier(**kwargs_rf)),\n",
    "                      \"gb\": make_pipeline(PCA(**kwargs_pca_unscaled),\n",
    "                                          GradientBoostingClassifier(**kwargs_gb)),\n",
    "                      }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitted_models_unscaled = {}\n",
    "model_short_names = [\"l1\", \"l2\", \"gb\", \"rf\"]\n",
    "#model_short_names = [\"l1\", \"l2\"]\n",
    "for name in model_short_names:\n",
    "    start_time = time()\n",
    "    \n",
    "    model = GridSearchCV(estimator=pipelines_unscaled.get(name),\n",
    "                         param_grid=hyperparameters.get(name),\n",
    "                         n_jobs=5,\n",
    "                         )\n",
    "    \n",
    "    model.fit(X_train, y_train),\n",
    "    \n",
    "    end_time = time()\n",
    "    tot_time = end_time - start_time\n",
    "    print(f\"{name} has been fitted in {tot_time / 3600:0.6f} hours!\")\n",
    "    \n",
    "    fitted_models_unscaled.update({name: model})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_results_unscaled = []\n",
    "for name, model in fitted_models_unscaled.items():\n",
    "    \n",
    "    best_score = model.best_score_\n",
    "    \n",
    "    pred_proba = model.predict_proba(X_test)[:, 1]\n",
    "    \n",
    "    auroc = roc_auc_score(y_test, pred_proba)\n",
    "    loss = log_loss(y_test, pred_proba)\n",
    "    \n",
    "    tmp = {\"model\": name,\n",
    "           \"accuracy_train\": best_score,\n",
    "           \"auroc_test\": auroc,\n",
    "           \"log_loss_test\": loss,\n",
    "           }\n",
    "    df_results_unscaled.append(tmp)\n",
    "\n",
    "df_results_unscaled = pd.DataFrame(df_results_unscaled)\n",
    "for score_name in [\"accuracy_train\", \"auroc_test\", \"log_loss_test\"]:\n",
    "    df_results_unscaled[f\"rank_{score_name}\"] = df_results_unscaled[score_name].rank(ascending=False).astype(int)\n",
    "\n",
    "df_results_unscaled"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}