data_exploration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Higgs Boson challenge\n",
    "\n",
    "# Exploratory data analysis - feature selection\n",
    "\n",
    "\n",
    "This dataset has been built from official ATLAS full-detector simulation, with \"Higgs to tautau\" events mixed with different backgrounds. The task is to classify events into \"tau tau decay of a Higgs boson\" versus \"background\".\n",
    "\n",
    "\n",
    "## Dataset characteristics\n",
    "\n",
    "[Kaggle Challenge Page](https://www.kaggle.com/c/higgs-boson)\n",
    "\n",
    "The dataset from Kaggle has 800000 events (195.5 MB in total):\n",
    "- Training set of 250000 events\n",
    "- Test set of 550000 events\n",
    "\n",
    "Training set has 30 feature columns, a weight column and a label column.\n",
    "Test set has 30 feature columns and a label column.\n",
    "\n",
    "\n",
    "### Feature characteristics\n",
    "\n",
    "- all variables are floating point, except PRI_jet_num which is integer\n",
    "- variables prefixed with PRI (for PRImitives) are “raw” quantities about the bunch collision as measured by the detector.\n",
    "- variables prefixed with DER (for DERived) are quantities computed from the primitive features, which were selected by  the physicists of ATLAS\n",
    "- it can happen that for some entries some variables are meaningless or cannot be computed; in this case, their value is −999.0, which is outside the normal range of all variables\n",
    "\n",
    "### Class distribution\n",
    "\n",
    "The class distribution of the training set is\n",
    "\n",
    "- b (background) : 164333 events (66%)\n",
    "- s (Higgs to tau tau):  85667 (34%)\n",
    "\n",
    "We can see there is a class imbalance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Usual imports\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "sns.set(rc={'figure.figsize':(9,6)})\n",
    "plt.rcParams['figure.figsize'] = (9,6)\n",
    "\n",
    "import IPython.display as ipd\n",
    "pd.options.display.max_columns = 35\n",
    "\n",
    "# random seed for reproducibility\n",
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"data/training.csv\",sep=',').astype({'Label': 'category'})\n",
    "print(\"Raw training data:\")\n",
    "ipd.display(df.tail(10))\n",
    "\n",
    "data = df.drop(['EventId', 'Label', 'Weight'], axis=1)\n",
    "kaggle_weights = df['Weight']\n",
    "# encode categorical class variable\n",
    "target = df['Label']\n",
    "# replace -990.0 values with NaNs\n",
    "data = data.replace(-999.0, np.NaN)\n",
    "weight = df['Weight']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of classes\n",
    "ipd.display(target.value_counts())\n",
    "ipd.display(data.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Missing values\n",
    "\n",
    "-990.0 values represent entries that have no meaning/that could not be calculated. The figure below shows the the percentage of those missing values per column.\n",
    "\n",
    "We can see than there are a lot of missing values: 72% of rows contain at least one missing value and  of some columns have up to 71% of them. So it is not viable to delete rows or columns, otherwise it will throw away too much data. Moreover the fact that a variable has no meaning or couldn't be calculated might be correlated with the label.\n",
    "\n",
    "The three possibilities we have to solve the problem of missing data are:\n",
    "- conserve -990.0 values: the value was chosen to be abnormal and distant from real data points; hence a model could be capable of identifying it and treat this value as missing.\n",
    "- replace with NaN values and use algorithms that manage missing values (for example XGBOOST)\n",
    "- impute those value swith the median (more robust than the man), this will smooth the data but may remove the explicative potential of those values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show percentage of nans per column\n",
    "perc_nans = data.isna().sum() * 100 / len(data)\n",
    "perc_rows_nans = (len(data) - len(data.dropna()))*100/(len(data))\n",
    "print(\"Percentage of rows containing at least one missing value:\", perc_rows_nans)\n",
    "plt.figure(figsize=(9,6))\n",
    "ax = perc_nans.plot.barh()\n",
    "ax.set_xlim(0, 100)\n",
    "ax.set_xticks(range(0, 101, 10))\n",
    "ax.tick_params(axis='y', labelsize=10)\n",
    "plt.title(\"Percentage of Missing Values Per Column\\n\", fontsize=15)\n",
    "#plt.savefig(\"figures/Percentage_of_Missing_Values_Per_Column.eps\", bbox_inches = \"tight\")\n",
    "plt.show()\n",
    "\n",
    "ipd.display(perc_nans)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Histograms\n",
    "\n",
    "To get a better understanding of the features distributions we plot the histograms for each features. \n",
    "\n",
    "Some variables seem to follow well known probability distributions:\n",
    "\n",
    "- Gaussians (ex. `PRI_tau_eta`, `PRI_lep_eta`)\n",
    "- Beta distribution with one parameter alpha (U distribution) (`DER_lep_eta_centrality`, `DER_met_phi_centrality`)\n",
    "- Uniform distributions (ex. `PRI_tau_phi`, `PRI_met_phi`)\n",
    "- Exponential distributions (ex. `DER_sum_pt`, `PRI_jet_leading_pt`)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show histograms of features\n",
    "data.hist(bins=30, figsize=(27,18))\n",
    "plt.suptitle(\"Features Histograms\", fontsize=30, y=0.94)\n",
    "#plt.savefig(\"figures/Histograms.svg\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Per-class histograms for primitive features and derived features\n",
    "\n",
    "Here, the features are separated into primitives (PRI) and derived (DER). For each feature we plot two histograms, one per class, to see any differences in distribution.\n",
    "\n",
    "Certain features, for ex. `DER_mass_MMC` and `DER_mass_transverse_met_lep` have significative differences in their distributions between classes, they will probably be discriminant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show histograms of primitives features\n",
    "\n",
    "pri = df.loc[:, df.columns.str.startswith(('PRI', 'Label'))]\n",
    "pri_columns = pri.drop('Label', axis=1).columns\n",
    "\n",
    "# show histograms of features\n",
    "#data.hist(bins=30, figsize=(27,18))\n",
    "#plt.show()\n",
    "\n",
    "def sephist(df, col):\n",
    "    s = df[df['Label'] == 's'][col]\n",
    "    b = df[df['Label'] == 'b'][col]\n",
    "    return s, b\n",
    "\n",
    "plt.figure(figsize=(27, 18),)\n",
    "\n",
    "for num, column in enumerate(pri_columns):\n",
    "    plt.subplot(4, 5, num+1)\n",
    "    sep = sephist(pri, column)\n",
    "    plt.hist(sep[0], bins=30, density=True, alpha=0.6, label='s')\n",
    "    plt.hist(sep[1], bins=30, density=True, alpha=0.6, label='b')\n",
    "    plt.legend(loc='best')\n",
    "    plt.title(column)\n",
    "\n",
    "# plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)\n",
    "plt.suptitle(\"Primitive Features Histograms\", fontsize=28, y=0.94)\n",
    "#plt.savefig(\"figures/Primitive_Features_Histograms.svg\")\n",
    "plt.show()\n",
    "\n",
    "# show histograms of primitives features\n",
    "\n",
    "der = df.loc[:, df.columns.str.startswith(('DER', 'Label'))]\n",
    "der_columns = der.drop('Label', axis=1).columns\n",
    "\n",
    "# show histograms of features\n",
    "#data.hist(bins=30, figsize=(27,18))\n",
    "#plt.show()\n",
    "\n",
    "def sephist(df, col):\n",
    "    s = df[df['Label'] == 's'][col]\n",
    "    b = df[df['Label'] == 'b'][col]\n",
    "    return s, b\n",
    "\n",
    "plt.figure(figsize=(21, 13))\n",
    "\n",
    "for num, column in enumerate(der_columns):\n",
    "    plt.subplot(4, 5, num+1)\n",
    "    sep = sephist(der, column)\n",
    "    plt.hist(sep[0], bins=30, density=True, alpha=0.6, label='s')\n",
    "    plt.hist(sep[1], bins=30, density=True, alpha=0.6, label='b')\n",
    "    plt.legend(loc='best')\n",
    "    plt.title(column)\n",
    "\n",
    "#plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)\n",
    "plt.suptitle(\"Derived Features Histograms\", fontsize=22, y=0.94)\n",
    "#plt.savefig(\"figures/Derived_Features_Histograms.svg\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlations matrices\n",
    "\n",
    "We compute the pairwise correlation of columns, excluding NaN values, with the Pearson's coefficient.  \n",
    "\n",
    "We can observe that several variables are highly correlated (for example `DER_sum_pt` with `PRI_jet_all_pt` (0.97)). As some variables are derived from others this is expected. There are also high correlations between primitives variables (ex. `PRI_met_sumet` with `PRI_jet_all_pt` (0.88))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation matrix\n",
    "corr = data.corr()\n",
    "corr.index = data.columns\n",
    "plt.figure(figsize=(30,30))\n",
    "sns.heatmap(corr, annot=True, cmap='RdYlGn', fmt='.2f')\n",
    "plt.title(\"Correlation Heatmap\\n\", fontsize=30)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pri = pri_data = data.loc[:, data.columns.str.startswith('PRI')]\n",
    "\n",
    "# Correlation matrix\n",
    "pri_corr = pri_data.corr()\n",
    "pri_corr.index = pri_data.columns\n",
    "plt.figure(figsize=(16,16))\n",
    "sns.heatmap(pri_corr, annot=True, cmap='RdYlGn', fmt='.2f')\n",
    "plt.title(\"Correlation Heatmap\\n\", fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis of the importance of features via a random forest classifier.\n",
    "\n",
    "To get a grasp of feature importances for classification, the best method we see is an \"a posteriori\" selection from the learning of a classifier capable of assigning importance to each of the features.\n",
    "We train a random forest with untuned hyperparams to get a scores of features importances.\n",
    "\n",
    "For each decision tree, the importance of a feature importance is computed as the (normalized) total reduction of the Grini criterion brought by that feature. It is also known as the Gini importance. It follows that the importances of features for the random forest are the mean of those importances for each tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.feature_selection import f_regression, mutual_info_regression\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "le = LabelEncoder()\n",
    "le.fit(target)\n",
    "target_enc = le.transform(target)\n",
    "\n",
    "#mi = mutual_info_regression(imp_data, target_enc)\n",
    "#mi /= np.max(mi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Forest Classifier to visualize feature importances\n",
    "\n",
    "By fitting classifier like a random forest on the data, we can observe the importance of features for the classification afterwards.\n",
    "\n",
    "For each decision tree, the importance of a feature importance is computed as the (normalized) total reduction of the criterion brought by that feature. It is also known as the Gini importance.\n",
    "It follows that the importances of features for the random forest are the mean of those importances for each tree. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.impute import SimpleImputer\n",
    "\n",
    "y = target.cat.codes\n",
    "imputer = SimpleImputer(strategy='median')\n",
    "X = imputer.fit_transform(data.values)\n",
    "\n",
    "clf = RandomForestClassifier(n_estimators=150, max_depth=10, n_jobs=-1)\n",
    "clf.fit(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "importances = clf.feature_importances_\n",
    "indices = np.argsort(importances)[::-1]\n",
    "std = np.std([tree.feature_importances_ for tree in clf.estimators_],\n",
    "             axis=0)[indices]\n",
    "# Print the feature ranking\n",
    "print(\"Feature ranking:\")\n",
    "importances = pd.Series(clf.feature_importances_, index=data.columns)[indices]\n",
    "\n",
    "# Plot the feature importances of the forest\n",
    "plt.figure(figsize=(9,6))\n",
    "plt.title(\"Feature Importances of Random Forest\", fontsize=15)\n",
    "\n",
    "importances.plot.bar(yerr=std, color='r')\n",
    "plt.gcf().subplots_adjust(bottom=0.45)\n",
    "\n",
    "#plt.savefig(\"figures/Feature_Importances.svg\")\n",
    "\n",
    "plt.show()\n",
    "\n",
    "#ax.tick_params(axis='y', labelsize=10)\n",
    "plt.show()"
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}