{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":50160,"databundleVersionId":7921029,"sourceType":"competition"}],"dockerImageVersionId":30665,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Introduction\n\nThis notebooks is comprised by\n\n* an exploratory data analysis, that describes the importance of some features on the probability of credit default by the applicant;\n\n* the training of a Gradient Boosted Decision Tree from [LightGBM](https://lightgbm.readthedocs.io/en/stable/), and tuning of its parameters through [scikit-learn](https://scikit-learn.org/stable/)'s [GridSearchCV](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html). Note that in order to avoid getting a model that overlooks the minority class (settled credit) while training, it was decided to balance the classes by applying scaling weights to the loss functions;\n\n* the computation of [Shapley values](https://en.wikipedia.org/wiki/Shapley_value) of some features using [SHAP](https://shap.readthedocs.io/en/latest/#) - to ascertain their relevance to the model.","metadata":{}},{"cell_type":"markdown","source":"# Import Python modules","metadata":{}},{"cell_type":"code","source":"# ---> General\n# warnings - to manage irrelevant warnings when runnig the code\nimport warnings\n# Ignore (that is, do not print) warnings\nwarnings.filterwarnings(\"ignore\")\n# matplotlib's pyplot - for general plotting\nimport matplotlib.pyplot as plt\n# Markdown (to output Markdown from Python code) \nfrom IPython.display import Markdown\n# time - to compute elapsed running time\nimport time\n\n# ---> Manage data\n# glob - to search for files whose names follow some pattern\nimport glob\n# NumPy - for basic mathematical operations\nimport numpy as np\n# polars - to efficiently manage (better than pandas) large quantities of data and to\n# deal with tables\nimport polars as pl\n# pandas - to convert polars dataframes to pandas after feature engineering (almost all\n# models still do not support polars)\nimport pandas as pd\n\n# ---> Exploratory Data Analysis (EDA)\n# Searborn - to plot statistical data\n# (the alias \"sns\" stands for \"Samuel Norman Seaborn\")\nimport seaborn as sns\n# collections - for alternative containers. Useful for counting\nimport collections\n\n# ---> Modelling\n# LightGBM - for applying gradient-boosting algorithms\nimport lightgbm as lgb\nimport xgboost as xgb\n# scikit-learn's GridSearchCV for hyperparameter tuning (using Grid Search and\n# cross-validation) and StratifiedGroupKFold for creating stratified grouped subfolds\n# for cross-validation\nfrom sklearn.model_selection import GridSearchCV\n# Import some scikit-learn's metrics \nfrom sklearn.metrics import (\n    accuracy_score,\n    roc_curve,\n    roc_auc_score,\n    classification_report,\n    confusion_matrix,\n    ConfusionMatrixDisplay,\n)\n\n# ---> Study on the Shapley values\nimport shap\n\n\nfrom sklearn.tree import DecisionTreeClassifier\nfrom sklearn.ensemble import RandomForestClassifier","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:44:14.48155Z","iopub.execute_input":"2024-05-05T16:44:14.481921Z","iopub.status.idle":"2024-05-05T16:44:14.65786Z","shell.execute_reply.started":"2024-05-05T16:44:14.481896Z","shell.execute_reply":"2024-05-05T16:44:14.656969Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Competition data in the kaggle kernel\n\nIf the notebook was created from the [competition page](https://www.kaggle.com/competitions/home-credit-credit-risk-model-stability/code), the respective data should be automatically loaded into the kaggle kernel. If other way was taken, the user would need to add the data through the button \"Add Input\" on the right-hand side pane, and choose the dataset for the project \"Home Credit - Credit Risk Model Stability\". The content may then be found in the directory `/kaggle/input/home-credit-credit-risk-model-stability/`.","metadata":{}},{"cell_type":"code","source":"# Path to data directory in the kaggle kernel\nPATH_DATA_ROOT = \"/kaggle/input/home-credit-credit-risk-model-stability/\"\n\n# List of names for the data batches\nbatches = [\"train\", \"test\"]\n\n# Dictionary of paths to training and test Parquet data directory in the kaggle kernel\n# [NOTE: it is usually much faster to handle Parquet files than CSV ones.]\nPATH_DATA = {\n    batch: f\"{PATH_DATA_ROOT}parquet_files/{batch}/\"\n    for batch in batches\n}","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:19.928217Z","iopub.execute_input":"2024-05-05T14:51:19.928826Z","iopub.status.idle":"2024-05-05T14:51:19.933896Z","shell.execute_reply.started":"2024-05-05T14:51:19.928798Z","shell.execute_reply":"2024-05-05T14:51:19.93288Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create polars dataframes from selected training and test data\n\nDataframes are found to be the most suitable objects for storing and managing table data in Python. The tables issued in this competition are quite huge. The popular [pandas](https://pandas.pydata.org/) package is not that efficient in dealing with such amount of data. Fortunately, there is another package, called [polars](https://pandas.pydata.org/) which is said to be [\"blazingly fast\"](https://docs.pola.rs/). Therefore, the polars package will be herein used in place of pandas' to create dataframes from the training and test data.\n\nThe data files are available in both CSV and Parquet formats. Since it is much faster to handle Parquet files than CSV ones, it was decided to take the former in place of the latter.\n\nThe information on the data files is posted on the [competition Data page](https://www.kaggle.com/competitions/home-credit-credit-risk-model-stability/data). It is also found in the file `feature_definitions.csv` of the input directory in this very kaggle kernel.\n\n## The data files\n\nThe Parquet files considered for the training and test dataset are described below. Furthermore, tables detailing the data are presented, having columns \"Used?\" and \"As?\" stating if the respective field was indeed used in the notebook or not, and how they were used (if the answer to the former is \"Yes\"), respectively.\n\nFor the sake of simplicity not all files were considered.\n\n* base files, `csv_files/train/train_base.csv` for the case of training and `csv_files/test/test_base.csv` for the case of testing. These describe the basic details of each credit case (a credit case corresponds to the credit contract and its state at the instant `WEEK_NUM`). The column `case_id` is the primary key of the table, uniquely identifying the credit case in this and all the other tables;\n        \n| Column          | Type | Description | Used? | As? |\n| --------------- | ---- | ----------- | ----- | --- |\n| `case_id`       | int  | A unique identifier for each credit case. | Yes | Case identifier |\n| `date_decision` | str  | Date of decision for the approval of the credit (format \"YYYY-MM-DD\"). | Yes | Auxiliar feature component, $x^{(j)}$ |\n| `MONTH`         | int  | Month of `date_decision` (format \"format YYYYMM\"). | No | --- |\n| `WEEK_NUM`      | int  | Number of weeks passed since `date_decision`. | Yes | Performance metric parameter |\n| `target`        | int  | Target (label) value ($1$ if the credit contract was defaulted, or $0$ if not). | Yes | Label, $y$ |\n\n> **_NOTE:_** in the case of the test file, there is no `target` column - test labels are not issued since these are to be predicted by the constructed model.\n\n<br>\n\n* static internal files, `csv_files/train/train_static_0_0.csv` and `csv_files/train/train_static_0_1.csv` for the case of training, and `csv_files/test/test_static_0_0.csv`, `csv_files/test/test_static_0_1.csv` and `csv_files/test/test_static_0_2.csv` for the case of testing. The files are static (without data being added as time runs) and they are internal (obtained from the credit contract itself). Note that the names of the files follow the format `*_<table_depth>_<table_part>`. `<table_depth>` is the depth of the table relationship tree (e.g. $2$ would mean that there would be a range of records (rows) for each `case_id` and another range of records for each of the records mentioned before). `<table_part>` is the part (set of records) of the whole table that is stored into a specific file (e.g. $2$ would mean \"part $2$\" of the whole table, implying that there would also be parts $0$ and $1$ in other files, necessarily). The static internal files have feature components of all types of transforms except \"T\" (see [competition Data page](https://www.kaggle.com/competitions/home-credit-credit-risk-model-stability/data) to know more about these transforms). All of them were used. In particular, three notable coloumns of the \"L\" type are herein detailed;\n\n| Column          | Type  | Description | Used? | As? |\n| --------------- | ----- | ----------- | ----- | --- |\n| `case_id`       | int   | Same as base files' `case_id`. | Yes | Case identifier |\n| `*P`     | float   | Feature component of P-type (\"transform days past due\"). | Yes | Feature component, $x^{(j)}$ |\n| `*M`     | str   | Feature component of M-type (\"masking categories\"). | Yes | Feature component, $x^{(j)}$ |\n| `*A`     | float | Feature component of A-type (\"transform amount\"). | Yes | Feature component, $x^{(j)}$ |\n| `*D`     | str   | Feature component of D-type (\"transform dates\"). | Yes | Feature component, $x^{(j)}$ |\n| `cntpmts24_3658933L`     | int | Number of monthly payments done in the last $24$ months and in the current one. | Yes | Auxiliar feature component, $x^{(j)}$ |\n| `mobilephncnt_593L` | int | Number of persons of the same contract using the same mobile phone number. | Yes | Feature component, $x^{(j)}$ |\n| `pmtnum_254L` | int | Total number of payments made by the applicant. | Yes | Feature component, $x^{(j)}$ |\n| `other *L`     | ---   | Feature component of L-type (\"unspecified transform\"). | Yes | Feature component, $x^{(j)}$ |\n\n<br>\n\n* static external (from a Credit Bureau) files, `csv_files/train/train_static_cb_0.csv` for the case of training, and `csv_files/test/test_static_cb_0.csv` for the case of testing. These files have feature components of all types of transforms except \"P\". All of them were considered;\n        \n| Column          | Type  | Description | Used? | As? |\n| --------------- | ----- | ----------- | ----- | --- |\n| `case_id`       | int   | Same as base files' `case_id`. | Yes | Case identifier |\n| `*M`     | str   | Feature component of M-type (\"masking categories\"). | Yes | Feature component, $x^{(j)}$ |\n| `*A`     | float | Feature component of A-type (\"transform amount\"). | Yes | Feature component, $x^{(j)}$ |\n| `*D`     | str   | Feature component of D-type (\"transform dates\"). | Yes | Feature component, $x^{(j)}$ |\n| `*L`     | ---   | Feature component of L-type (\"unspecified transform\"). | Yes | Feature component, $x^{(j)}$ |\n| `*T`     | ---   | Feature component of T-type (\"unspecified transform\"). | Yes | Feature component, $x^{(j)}$ |\n\n<br>\n\n* contract persons' data files, `csv_files/train/train_person_1.csv` for the case of training, and `csv_files/test/test_person_1.csv` for the case of testing. These files have feature components of all types of transforms except \"P\". All of them were considered. In particular, six feature components are herein detailed. Furthermore, for compactness reasons (note that the depth of this table is greater than $0$), the entries of these features were grouped by `case_id` and \"aggregated\" through a convenient operation to obtain a singular value from each group, as described in the table below;\n\n| Column                  | Type  | Description | Used?| As? |\n| ----------------------- | ----- | ----------- | ---- | --- |\n| `case_id`               | int   | Same as base files' `case_id`. | Yes | Case identifier |\n| `num_group1`            | int   | Identifier of the person in the respective credit case ($0$ identifies the applicant). | Yes | Person identifier |\n| `mainoccupationinc_max_A` <br> (derived) | float | Maximum income (from main occupation) of the set of people associated with each `case_id` group. This column derives from the `max` operation on the grouped values of the original column `mainoccupationinc_384A`. | Yes | Feature component, $x^{(j)}$ |\n| `anyselfemployed_T` <br> (derived)       | bool  | Boolean that holds `True` if any of the set of people of each `case_id` group is self-employed. This column derives from the `any` operation on the outcome of `== \"SELFEMPLOYED\"` applied to the grouped values of the original column `incometype_1044T`. | Yes | Feature component, $x^{(j)}$ |\n| `housetype_applicant_905L` <br> (derived)   | str   | House type of the applicant of each `case_id` group. This column derives from the values of the original column `housetype_905L` for which the values of the column `num_group1` are $0$.| Yes | Feature component, $x^{(j)}$ |\n| `sex_applicant_738L` <br> (derived)   | str   | Gender of the applicant of each `case_id` group. This column derives from the values of the original column `sex_738L` for which the values of the column `num_group1` are $0$.| Yes | Feature component, $x^{(j)}$ |\n| `birth_applicant_259D` <br> (derived)   | str   | Birth date of the applicant of each `case_id` group. This column derives from the values of the original column `birth_259D` for which the values of the column `num_group1` are $0$.| Yes | Auxiliar feature component, $x^{(j)}$ |\n| `empl_employedfrom_applicant_271D` <br> (derived)   | str   | Starting employment date of the applicant of each `case_id` group. This column derives from the values of the original column `empl_employedfrom_271D` for which the values of the column `num_group1` are $0$.| Yes | Auxiliar feature component, $x^{(j)}$ |\n| `*_applicant_*M` <br> (derived)    | str   | Feature component of M-type (\"masking categories\") associated with the applicant (that is, for which `num_group1` is $0$). | Yes | Feature component, $x^{(j)}$ |\n| `other *_applicant_*A` <br> (derived)   | float | Feature component of A-type (\"transform amount\") associated with the applicant (that is, for which `num_group1` is $0$). | Yes | Feature component, $x^{(j)}$ |\n| `other *_applicant_*D` <br> (derived)    | str   | Feature component of D-type (\"transform dates\") associated with the applicant (that is, for which `num_group1` is $0$). | Yes | Feature component, $x^{(j)}$ |\n| `other *_applicant_*L` <br> (derived)    | ---   | Feature component of L-type (\"unspecified transform\") associated with the applicant (that is, for which `num_group1` is $0$). | Yes | Feature component, $x^{(j)}$ |\n| `other *_applicant_*T` <br> (derived)    | ---   | Feature component of T-type (\"unspecified transform\") associated with the applicant (that is, for which `num_group1` is $0$). | Yes | Feature component, $x^{(j)}$ |    \n<br>\n    \n* Credit Bureau B's data files, `csv_files/train/train_credit_bureau_b_2.csv` for the case of training, and `csv_files/test/test_credit_bureau_b_2.csv` for the case of testing. These files have feature components of all types of transforms. As for the case of the contract persons' data files, the entries of some of the columns were grouped by `case_id` and \"aggregated\".\n\n| Column                    | Type  | Description | Used? | As? |\n| ------------------------- | ----- | ----------- | ----- | --- |\n| `case_id`                 | int   | Same as base files' `case_id`. | Yes | Case identifier |\n| `num_group1`              | int   | Identifier of the contract in the respective credit case. | No | --- |\n| `num_group2`              | int   | Identifier of the payment in the respective credit case. | No | --- |\n| `pmts_date_1107D`         | date  | Payment date. | No | --- |\n| `pmts_dpdvalue_anyover31_P` <br> (derived) | bool  | Boolean that holds `True` if there is any number of days of overdue payment greater than $31$ for each `case_id` group. This column derives from the `any` operation on the outcome of `> 31` applied to the grouped values of the original column `pmts_dpdvalue_108P`. | Yes | Feature component, $x^{(j)}$ |\n| `pmts_pmtsoverdue_max_A` <br> (derived) | float | Maximum value of the number of overdue payments for each `case_id` group. This columns derives from the `max` operation on the grouped values of the original column `pmts_pmtsoverdue_635A`. | Yes | Feature component, $x^{(j)}$ |","metadata":{}},{"cell_type":"code","source":"# ---> Function for setting polars dataframes dtypes\ndef set_pl_dtypes(df):\n    for col in df.columns:\n        # If column is associated with P or A-type transforms, set its dtype to to\n        # Float64\n        if col[-1] in (\"P\", \"A\"):\n            df = df.with_columns(pl.col(col).cast(pl.Float64))\n        # If column is associated with M-type transform, set its dtype to to Categorical\n        # [NOTE: the ordering style is set to \"lexical\" so that categories are sorted\n        # according to the string value aka \"alphabet\", instead of order of appearance\n        # in the dataframe (\"physical\").]\n        if col[-1] in (\"M\"):\n            df = df.with_columns(pl.col(col).cast(pl.Categorical(\"lexical\")))\n        # If column is associated with D-type transform, set its dtype to to Date\n        if col[-1] in (\"D\"):\n            df = df.with_columns(pl.col(col).cast(pl.Date))\n        # If column is associated with L or T-type transform, and it is not of numeric\n        # dtype, set it to Categorical\n        if col[-1] in (\"L\", \"T\") and not df[col].is_numeric:\n            df = df.with_columns(pl.col(col).cast(pl.Categorical(\"lexical\")))\n    return df","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:19.93496Z","iopub.execute_input":"2024-05-05T14:51:19.935209Z","iopub.status.idle":"2024-05-05T14:51:19.949745Z","shell.execute_reply.started":"2024-05-05T14:51:19.935186Z","shell.execute_reply":"2024-05-05T14:51:19.94891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Create polars dataframes with training and test data\n\n# Dictionary of polars dataframes\ndt_data = {\n    batch: {\n        # Define polars dataframe containing data from base Parquet file\n        # [NOTE: column \"date_decision\" is converted to polars' Date type.]\n        \"base\": (pl.read_parquet(f\"{PATH_DATA[batch]}{batch}_base.parquet\")\\\n                 .with_columns(pl.col(\"date_decision\").cast(pl.Date))),\n        # Define polars dataframe containing data from static internal Parquet files\n        \"static\": pl.concat([pl.read_parquet(PATH_DATA_STATIC)\\\n                             .pipe(set_pl_dtypes) for PATH_DATA_STATIC in\n                             glob.glob(f\"{PATH_DATA[batch]}{batch}_static_0*.parquet\")],\n                            # Vertically concatenate, while conveniently redefining\n                            # column's data types to support the data from common\n                            # columns\n                            how=\"vertical_relaxed\"),\n        # Define polars dataframe containing data from static external (from a Credit\n        # Bureau) Parquet files\n        \"static_cb\": (pl.read_parquet(f\"{PATH_DATA[batch]}{batch}_static_cb_0.parquet\")\\\n                      .pipe(set_pl_dtypes)),\n        # Define polars dataframe containing data from the persons' Parquet file of\n        # depth 1\n        \"person_1\": (pl.read_parquet(f\"{PATH_DATA[batch]}{batch}_person_1.parquet\")\\\n                     .pipe(set_pl_dtypes)),\n        # Define polars dataframe containing data from Credit Bureau B's Parquet file of\n        # depth 2\n        \"credit_bureau_b_2\": (pl.read_parquet(f\"{PATH_DATA[batch]}{batch}_credit_bureau_b_2.parquet\")\\\n                              .pipe(set_pl_dtypes))\n\n    }   \n    for batch in batches\n}","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:19.952279Z","iopub.execute_input":"2024-05-05T14:51:19.952695Z","iopub.status.idle":"2024-05-05T14:51:31.981024Z","shell.execute_reply.started":"2024-05-05T14:51:19.952663Z","shell.execute_reply":"2024-05-05T14:51:31.98004Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Filter training and test dataframes\n\nfor batch in batches:\n\n    # Group person_1 dataframe by case_id and aggregate the maximum income (from main\n    # occupation) of the set of people associated with each group, as well as a boolean\n    # that holds true if any of the set of people of each group is self-employed.\n    df_person_1_1 = dt_data[batch][\"person_1\"].group_by(\"case_id\").agg(\n        pl.col(\"mainoccupationinc_384A\").max().alias(\"mainoccupationinc_max_A\").cast(pl.Float64),\n        (pl.col(\"incometype_1044T\") == \"SELFEMPLOYED\").any().alias(\"anyselfemployed_T\").cast(pl.Boolean)\n    )\n    # Select from person_1 dataframe all the other columns. Keep solely rows associated\n    # with applicants (num_group1 = 0), drop the column \"num_group1\" and rename all\n    # columns except \"case_id\" to now refer to the applicant\n    df_person_1_2 = dt_data[batch][\"person_1\"].drop(\n        [\"mainoccupationinc_384A\", \"incometype_1044T\"])\\\n        .filter(pl.col(\"num_group1\") == 0).drop(\"num_group1\")\n    df_person_1_2 = df_person_1_2.rename(\n        {col: col.rsplit(\"_\",1)[0] + \"_applicant_\" + col.rsplit(\"_\",1)[1] for\n         col in df_person_1_2.drop(\"case_id\").columns}\n    )\n    # Redefine person_1 dataframe as the auxiliar df_person_1_1 and df_person_1_2 (left)\n    # joined together\n    dt_data[batch][\"person_1\"] = (df_person_1_1\\\n        .join(other=df_person_1_2,\n              how=\"left\",\n              on=\"case_id\")\n    )\n    \n    # Group credit_bureau_b_2 dataframe by case_id and aggregate maximum value of the\n    # number of overdue payments (pmts_pmtsoverdue_635A) and a boolean that holds true\n    # if there is any number of days of overdue payment greater than 31.\n    dt_data[batch][\"credit_bureau_b_2\"] = dt_data[batch][\"credit_bureau_b_2\"]\\\n        .group_by(\"case_id\").agg(\n            pl.col(\"pmts_pmtsoverdue_635A\").max().alias(\"pmts_pmtsoverdue_max_A\").cast(pl.Float64),\n            (pl.col(\"pmts_dpdvalue_108P\") > 31).any().alias(\"pmts_dpdvalue_anyover31_P\").cast(pl.Boolean)\n    )","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:31.983186Z","iopub.execute_input":"2024-05-05T14:51:31.983508Z","iopub.status.idle":"2024-05-05T14:51:34.48008Z","shell.execute_reply.started":"2024-05-05T14:51:31.983481Z","shell.execute_reply":"2024-05-05T14:51:34.479248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Collapse training and test dataframes through left join of their subdataframes\ndt_data = {\n    batch: dt_data[batch][\"base\"]\\\n    .join(dt_data[batch][\"static\"],\n          how=\"left\",\n          on=\"case_id\")\\\n    .join(dt_data[batch][\"static_cb\"],\n          how=\"left\",\n          on=\"case_id\")\\\n    .join(dt_data[batch][\"person_1\"],\n          how=\"left\",\n          on=\"case_id\")\\\n    .join(dt_data[batch][\"credit_bureau_b_2\"],\n          how=\"left\",\n          on=\"case_id\")\n    for batch in batches\n}","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:34.48128Z","iopub.execute_input":"2024-05-05T14:51:34.481607Z","iopub.status.idle":"2024-05-05T14:51:36.458121Z","shell.execute_reply.started":"2024-05-05T14:51:34.48158Z","shell.execute_reply":"2024-05-05T14:51:36.45732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Display training and test dataframes","metadata":{}},{"cell_type":"code","source":"# Display first five entries of the training dataframe\ndt_data[\"train\"].head()","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:36.459205Z","iopub.execute_input":"2024-05-05T14:51:36.459504Z","iopub.status.idle":"2024-05-05T14:51:36.47954Z","shell.execute_reply.started":"2024-05-05T14:51:36.45948Z","shell.execute_reply":"2024-05-05T14:51:36.478491Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Display first five entries of the test dataframe\ndt_data[\"test\"].head()","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:36.48073Z","iopub.execute_input":"2024-05-05T14:51:36.481386Z","iopub.status.idle":"2024-05-05T14:51:36.501543Z","shell.execute_reply.started":"2024-05-05T14:51:36.48136Z","shell.execute_reply":"2024-05-05T14:51:36.50051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Exploratory Data Analysis (EDA)","metadata":{}},{"cell_type":"markdown","source":"## Verify if the training data is balanced","metadata":{}},{"cell_type":"code","source":"# ---> Check if the training data is balanced\n# Dictionary of counts of the different classes in the training data\ndt_N_y_train = dict(collections.Counter(dt_data[\"train\"][\"target\"]))\n\n# Ratio between numbers of settled (y=0) and defaulted (y=1) credit cases\ndt_N_y_train[\"ratio\"] = dt_N_y_train[0] / dt_N_y_train[1]\n\n# Initialise figure and axes\nplt.figure(figsize=(6.4, 4.8))\nax = plt.axes()\nplt.title(\"Number of counts per class\", pad=20)\n\n# Define count plot for the occurrence of the different digits in the training data    \nsns.countplot(\n    ax=ax,\n    x=dt_data[\"train\"][\"target\"].to_numpy(),\n    color=\"blue\",\n    alpha=0.5,\n    edgecolor=\"black\",\n    linewidth=1.0,\n    width=0.075,\n    hatch=\"////\",\n    zorder=2\n)\n\n# Define axes labels                                \nax.set_xlabel(r\"Class, $y$\", fontdict={\"fontsize\": 10})\nax.set_ylabel(r\"Counts\", fontdict={\"fontsize\": 10})\n\n# Enable axes' minor ticks\nax.minorticks_on()\n\n# Define grid\nax.grid(\n    visible=True,\n    which=\"major\",\n    color=\"lightgray\",\n    linestyle=\"solid\",\n    linewidth=0.5\n)\nax.grid(\n    visible=True,\n    which=\"minor\",\n    color=\"lightgray\",\n    linestyle=\"dotted\",\n    linewidth=0.5\n)\n\n# Show plot\nplt.show() \n\n# Display ratio between numbers of settled (y=0) and defaulted (y=1) credit cases\nprint()\ndisplay(pd.DataFrame(data={\"$n^-/n^+$\": dt_N_y_train[\"ratio\"]},\n                     index=[0])\\\n        .style\\\n        .format({\"$n^-/n^+$\": \"{:.2f}\"})\\\n        .set_caption(\"Ratio between numbers of settled ($y=0$) and defaulted ($y=1$) credit contract cases\")\\\n        .set_table_styles([\n                # Column width\n                {\"selector\": \"th.col_heading,td\",\n                 \"props\": [(\"width\", \"300px\")]\n                 },\n                # Caption style\n                {\"selector\": \"caption\",\n                 \"props\": [(\"font-size\", \"16px\"),\n                           (\"font-weight\", \"bold\"),\n                           (\"font-style\", \"italic\")]\n                 }\n            ]))","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:36.502595Z","iopub.execute_input":"2024-05-05T14:51:36.502857Z","iopub.status.idle":"2024-05-05T14:51:37.338669Z","shell.execute_reply.started":"2024-05-05T14:51:36.502834Z","shell.execute_reply":"2024-05-05T14:51:37.33774Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As expected the training data is imbalanced - there are much less cases of \"defaulted credit contracts\" than \"settled credit contracts\". The number of the latter is $30.81$ times higher than the number of the former.\n\nSince the training data is imbalanced, if training is done in the \"standard\" way, the respective obtained model could exhibit biased performance tending to favor the majority class and overlook the minority one. An approach for making the training of the model to identically regard the majority and minority classes is to consider weights for the loss functions of the minority class. The weights would need to be such that the number of training points associated with the minority class scaled by these weights coincides witht he number of training points associated with the majority class.\n\nLet's mathematicaly describe this balancing approach. Let $y_i$ be the label of the $i$-th training point. The values $y=0$ and $y=1$ represent the minority and majority classes (\"settled credit\" and \"defaulted credit\"), respectively. Also, let $L_i$ be the loss function associated with the $i$-th point, and $w_i$ the respective weight. The cost function $C$ would then correspond to\n\n$$C = \\sum_{i=1}^{n} w_i \\cdot L_i\\text{ ,}$$\n\nwhere $n$ is the number of training points. And the weight $w_i$ would be defined as\n\n$$\nw_i=\n\\begin{cases}\n1\\text{, } & \\text{ if }y_i=0\\\\\n\\frac{n^-}{n^+}\\text{, } & \\text{ if }y_i=1\n\\end{cases}\n$$\n\nwhere $n^-$ and $n^+ = n - n^-$ are the numbers of training points of the classes \"settled credit\" and \"defaulted credit\", respectively.\n\nTraining is said to be \"balanced\" if the part of the cost function associated with one class is equal to the part associated with the other in the case of the loss function being uniform (with $L_i=L,\\,\\forall i=1,\\,\\dots,\\,n$). The usage of the above-mentioned weights indeed produces such result:\n\n$$\nC = \\sum_{i=1}^{n} w_i \\cdot L_i = \\left(\\sum_{i=0 \\, \\land \\, y_i=0} w_i \\cdot L_i\\right) +  \\left(\\sum_{i=0 \\, \\land \\, y_i=1} w_i \\cdot L_i\\right) = \\underbrace{\\left(\\sum_{i=0 \\, \\land \\, y_i=0} 1 \\right)}_{=n^-} \\cdot L + \\underbrace{\\left(\\sum_{i=0 \\, \\land \\, y_i=1} 1 \\right)}_{=n^+} \\cdot \\frac{n^-}{n^+} \\cdot L =\n$$\n\n$$\n= \\underbrace{\\left(n^- \\cdot L\\right)}_{\\text{Majority class cost}} + \\underbrace{\\left(n^- \\cdot L\\right)}_{\\text{Minority class cost}}\\text{ .}\\quad\\quad\\quad \\scriptsize{■}\n$$","metadata":{}},{"cell_type":"code","source":"# ---> Add a sample weight column to the training dataframe\ndt_data[\"train\"] = dt_data[\"train\"].with_columns(\n    pl.when(pl.col(\"target\") == 1).then(dt_N_y_train[\"ratio\"]).otherwise(1)\\\n    .alias(\"sample_weight\").cast(pl.Float64)\n)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:37.342636Z","iopub.execute_input":"2024-05-05T14:51:37.342925Z","iopub.status.idle":"2024-05-05T14:51:37.360233Z","shell.execute_reply.started":"2024-05-05T14:51:37.3429Z","shell.execute_reply":"2024-05-05T14:51:37.359503Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Check columns emptiness\n\nLet one get an histogram of counts of columns distributed along bins of emptiness value for the training dataset. Emptiness is herein defined as the ratio between the number of [missing values](https://docs.pola.rs/user-guide/expressions/missing-data/) (polars' `null` and `NaN` values) and the number of entries in the column.\n\nColumns with emptiness higher than $99.5\\%$ are considered to be \"almost empty\" and, therefore, dropped. The same ones in the test dataset are identically dropped.","metadata":{}},{"cell_type":"code","source":"# ---> Get columns' emptiness value in the training set\n\n# Dictionary describing columns emptiness\ndt_emp = {\n    # Columns' emptiness value\n    \"emptiness\": {col: dt_data[\"train\"][col].null_count() / len(dt_data[\"train\"][col]) for\n                  col in dt_data[\"train\"].columns},\n    # Columns' almost emptiness indicator: if true, more than 99.5 % of the entries are\n    # empty\n    \"almost_empty\": {col: dt_data[\"train\"][col].null_count() / len(dt_data[\"train\"][col]) > 0.995 for\n                     col in dt_data[\"train\"].columns}\n}","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:37.361322Z","iopub.execute_input":"2024-05-05T14:51:37.361682Z","iopub.status.idle":"2024-05-05T14:51:37.372038Z","shell.execute_reply.started":"2024-05-05T14:51:37.361651Z","shell.execute_reply":"2024-05-05T14:51:37.37117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Plot emptiness distribution of the training set\n\n# Initialise figure and axes\nplt.figure(figsize=(6.4, 4.8))\nax = plt.axes()\nplt.title(\"Number of counts per emptiness range in the training dataset\", pad=20)\n\n# Define histogram plot for counts (of columns) versus emptiness range\nsns.histplot(\n    ax=ax,\n    data=dt_emp[\"emptiness\"].values(),\n    stat=\"count\",\n    bins=20,\n    binrange=(0, 1),\n    palette=[\"blue\"],\n    alpha=0.5,\n    edgecolor=\"black\",\n    linewidth=1.0,\n    hatch=\"////\",\n    zorder=2,\n    legend=False\n)\n\n# Plot verical line emptiness = 0.995\nplt.axvline(\n    x=0.995,\n    color=\"black\",\n    linewidth=2,\n    alpha=1,\n    linestyle=\"dashed\",\n    label=\"$\\mathrm{emptiness}=0.995$\"\n)\n\n# Define axes labels                                \nax.set_xlabel(r\"Emptiness\", fontdict={\"fontsize\": 10})\nax.set_ylabel(r\"Counts\", fontdict={\"fontsize\": 10})\n\n# Define axes limits\nax.set_xlim(\n    left=0,\n    right=1\n)\n    \n# Enable axes' minor ticks\nax.minorticks_on()\n\n# Define grid\nax.grid(\n    visible=True,\n    which=\"major\",\n    color=\"lightgray\",\n    linestyle=\"solid\",\n    linewidth=0.5\n)\nax.grid(\n    visible=True,\n    which=\"minor\",\n    color=\"lightgray\",\n    linestyle=\"dotted\",\n    linewidth=0.5\n)\n\n# Define legend\nax.legend(fontsize=8)\n\n# Show plot\nplt.show() \n\n# ---> Display number of \"almost empty\" columns\nprint()\ndisplay(pd.DataFrame(data={\"Counts ($\\mathrm{emptiness}>0.995$)\":\n                           sum(dt_emp[\"almost_empty\"].values())},\n                     index=[0])\\\n        .style\\\n        .format({\"Counts ($\\mathrm{emptiness}>0.995$)\": \"{:d}\"})\\\n        .set_caption(\"Number of almost empty columns in the training dataset\")\\\n        .set_table_styles([\n                # Column width\n                {\"selector\": \"th.col_heading,td\",\n                 \"props\": [(\"width\", \"300px\")]\n                 },\n                # Caption style\n                {\"selector\": \"caption\",\n                 \"props\": [(\"font-size\", \"16px\"),\n                           (\"font-weight\", \"bold\"),\n                           (\"font-style\", \"italic\")]\n                 }\n            ]))","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:37.37334Z","iopub.execute_input":"2024-05-05T14:51:37.37361Z","iopub.status.idle":"2024-05-05T14:51:37.887985Z","shell.execute_reply.started":"2024-05-05T14:51:37.373587Z","shell.execute_reply":"2024-05-05T14:51:37.887174Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Drop almost empty columns\ncols_drop = [key for key in dt_emp[\"almost_empty\"].keys() if dt_emp[\"almost_empty\"][key] == True]\ndt_data[\"train\"] = dt_data[\"train\"].drop(cols_drop)\ndt_data[\"test\"] = dt_data[\"test\"].drop(cols_drop)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:37.889234Z","iopub.execute_input":"2024-05-05T14:51:37.889563Z","iopub.status.idle":"2024-05-05T14:51:37.895386Z","shell.execute_reply.started":"2024-05-05T14:51:37.889539Z","shell.execute_reply":"2024-05-05T14:51:37.894538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Check categorical columns' number of categories\n\nLet one plot the distribution of the number of categories in the categorical columns for the training dataset. Categorical columns with solely $1$ category are redundant. And categorical columns with too many categories (e.g more than $1000$) might be hard to handle. These will be, therefore, dropped, as well as the respective ones in the test dataset.","metadata":{}},{"cell_type":"code","source":"# ---> Get categorical columns' number of categories in the training set\n\n# List of categorical columns\ncols_cat = dt_data[\"train\"].select(pl.col(pl.Categorical)).columns\n\n# Dictionary describing categorical columns' number of categories \ndt_N_cat = {\n    # Columns' numbers of categories\n    \"N_cat\": {col: len(dt_data[\"train\"][col].cat.get_categories()) for\n              col in cols_cat},\n}\n\ndt_N_cat.update({\n    # Columns' indicator of single category: if true, the number of categories is 1\n    \"N_cat==1\": {col: dt_N_cat[\"N_cat\"][col] == 1 for col in cols_cat},\n    # Columns' indicator of admissible number of categories: if true, the number of\n    # categories is greater than 1 and smaller or equal to 1000\n    \"N_cat>1 and N_cat<=1\": {col: dt_N_cat[\"N_cat\"][col] > 1 and\n                             dt_N_cat[\"N_cat\"][col] <= 1000 for col in cols_cat},\n    # Columns' indicator of too many categories: if true, the number of categories is\n    # greater than 100\n    \"N_cat>1000\": {col: dt_N_cat[\"N_cat\"][col] > 1000 for col in cols_cat}\n})","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:37.896658Z","iopub.execute_input":"2024-05-05T14:51:37.896979Z","iopub.status.idle":"2024-05-05T14:51:37.908586Z","shell.execute_reply.started":"2024-05-05T14:51:37.896949Z","shell.execute_reply":"2024-05-05T14:51:37.907802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Plot number of categories distribution of the training set\n\n# Bar points' unique abscissas\nx = list(range(len([\"N_cat==1\", \"N_cat>1 and N_cat<=1\", \"N_cat>1000\"]) + 1))\n\n# Indices of x\ni_x = list(range(len(x)))\n# Bar points' abscissas\nx_bar_points = sum(\n    [[x[i]] * \n     (2 if i in (0, i_x[-1]) else 4)\n     for i in i_x],\n    []\n)\n\n# Bar points' ordinates\ny_bar_points = sum(\n    [[0, sum(dt_N_cat[key].values()), sum(dt_N_cat[key].values()), 0]\n     for key in [\"N_cat==1\", \"N_cat>1 and N_cat<=1\", \"N_cat>1000\"]],\n    []\n)\n\n# x-axis' tick positions\nx_ticks = [(x[i] + x[i + 1]) / 2 for i in i_x[:-1]]\n\n# Initialise figure and axes\nplt.figure(figsize=(6.4, 4.8))\nax = plt.axes()\nplt.title(\"Number of counts per number of categories range in the training dataset\", pad=20)\n\n# Plot bar lines\nax.plot(\n    x_bar_points,\n    y_bar_points,\n    color=\"black\",\n    linewidth=1.0,\n    zorder=2\n)\n\n# Fill bars\nax.fill_between(\n    x=x_bar_points,\n    y1=np.zeros(len(x_bar_points)),\n    y2=y_bar_points,\n    facecolor=\"blue\",\n    edgecolor=\"black\",\n    hatch=\"////\",\n    linewidth=1.0,\n    alpha=0.5,\n    zorder=2\n)\n        \n# Define axes labels                                \nax.set_xlabel(r\"Number of categories, $N_{\\mathrm{cat}}$\", fontdict={\"fontsize\": 10})\nax.set_ylabel(r\"Counts\", fontdict={\"fontsize\": 10})\n\n# Define axes limits\nax.set_xlim(\n    left=0,\n    right=3\n)\nax.set_ylim(\n    bottom=0\n)\n\n# Define x axis' ticks\nax.set_xticks(x_ticks, labels=[\"$1$\", \"$]1,\\,1000]$\", f\"$>1000$\"])\n\n# Define grid\nax.grid(\n    visible=True,\n    which=\"major\",\n    color=\"lightgray\",\n    linestyle=\"solid\",\n    linewidth=0.5\n)\nax.grid(\n    visible=True,\n    which=\"minor\",\n    color=\"lightgray\",\n    linestyle=\"dotted\",\n    linewidth=0.5\n)\n\n# Show plot\nplt.show() \n\n# ---> Display number of \"almost empty\" columns\nprint()\ndisplay(pd.DataFrame(data={\"Counts ($N_{\\mathrm{cat}}=1$)\": sum(dt_N_cat[\"N_cat==1\"].values()),\n                           \"Counts ($N_{\\mathrm{cat}}>1000$)\": sum(dt_N_cat[\"N_cat>1000\"].values())},\n                     index=[0])\\\n        .style\\\n        .format({\"N\": \"{:d}\"})\\\n        .set_caption(\"Numbers of columns of single and of too many categories in the training dataset\")\\\n        .set_table_styles([\n                # Column width\n                {\"selector\": \"th.col_heading,td\",\n                 \"props\": [(\"width\", \"300px\")]\n                 },\n                # Caption style\n                {\"selector\": \"caption\",\n                 \"props\": [(\"font-size\", \"16px\"),\n                           (\"font-weight\", \"bold\"),\n                           (\"font-style\", \"italic\")]\n                 }\n            ]))","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:37.909988Z","iopub.execute_input":"2024-05-05T14:51:37.910306Z","iopub.status.idle":"2024-05-05T14:51:38.198076Z","shell.execute_reply.started":"2024-05-05T14:51:37.910277Z","shell.execute_reply":"2024-05-05T14:51:38.19714Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Drop single and too many categories columns\ncols_drop = [col for col in cols_cat if\n             True in (dt_N_cat[\"N_cat==1\"][col], dt_N_cat[\"N_cat>1000\"])]\n\ndt_data[\"train\"] = dt_data[\"train\"].drop(cols_drop)\ndt_data[\"test\"] = dt_data[\"test\"].drop(cols_drop)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:38.199555Z","iopub.execute_input":"2024-05-05T14:51:38.200243Z","iopub.status.idle":"2024-05-05T14:51:38.207527Z","shell.execute_reply.started":"2024-05-05T14:51:38.200207Z","shell.execute_reply":"2024-05-05T14:51:38.206614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Get probability distributions associated with notable features\n\nAccording to @[paddykb](https://www.kaggle.com/paddykb)'s [study](https://www.kaggle.com/competitions/home-credit-credit-risk-model-stability/discussion/473950#2644886), some of the most impactful (on the contract state) features are\n\n* contract's monthly payment value (`annuity_780A` of `static_0.parquet` files);\n\n* applicant's age at decision date (inferred from`birth_259D` (birth date) of `person_1.parquet` files, and `date_decision` of `base.parquet` files);\n\n* applicant's employment time at decision date (inferred from `empl_employedfrom_271D` of `person_1.parquet` files, and `date_decision` of `base.parquet` files);\n\n* number of months with missing payment in the last 24 months and current one (inferred from `cntpmts24_3658933L` of `static_0.parquet`);\n\n* applicant's gender (`sex_738L` of `person_1.parquet` files);\n\nas well as\n\n* the number of persons of the same contract using the same mobile phone number (`mobilephncnt_593L` of `static_0.parquet` files);\n\n* weekday of decision for the approval of the credit (inferred from`date_decision` (decision date) of `base.parquet` files);\n\n* total number of payments made by the applicant (`pmtnum_254L` of `static_0.parquet` files).\n\nLet one plot the probability of credit default $P:=P(y=1)$ for different intervals of values of these features $x^{(j)}$.","metadata":{}},{"cell_type":"code","source":"# ---> Function for distribution plotting\ndef dist_plot(df, col, col_fancy, binning=True, N_bins=10, xticklabels_rotation=0):\n\n    # ---------------------------------- Prepare --------------------------------- #\n    if binning == True:\n        # Array of bins' limits\n        x_bin_limit = np.linspace(df[col].min(), df[col].max(), N_bins + 1)\n\n        # List of indices of the array of bins' limits\n        i_x_bin_limit = range(len(x_bin_limit))\n\n        # Polars dataframe with distributions of counts and probabilities\n        df_dist = pl.DataFrame(\n            schema=[\n                \"Bin's $x_{\\mathrm{min}}$\",\n                \"Bin's $x_{\\mathrm{av}}$\",\n                \"Bin's $x_{\\mathrm{max}}$\",\n                \"Counts ($y=0$)\",\n                \"Counts ($y=1$)\",\n                \"$P(y=1)\\,[\\%]$\"\n            ]\n        )\n\n        # Polars dataframe with distributions of probabilities for plotting in a bar style\n        df_dist_plot = pl.DataFrame(\n            schema=[\n                \"x\",\n                \"y\"\n            ]\n        )\n\n        # Construct distribution dataframes\n        # [NOTE: the dataframes are constructed in a way that bins with no counts are not\n        # accounted.]\n        for i in i_x_bin_limit[:-1]:\n            df_dist_i = (\n                df.select(col, \"target\")\\\n                    # Group by bin range\n                    .groupby(\n                        ((pl.col(col) >= x_bin_limit[i]) &\n                         ((pl.col(col) < x_bin_limit[i + 1]) if \n                          i != i_x_bin_limit[-2] else \n                          (pl.col(col) <= x_bin_limit[i + 1])))\n                    )\\\n                    # Aggregate to compute counts\n                    .agg(\n                        (pl.col(\"target\") == 0).sum().alias(\"Counts ($y=0$)\"),\n                        (pl.col(\"target\") == 1).sum().alias(\"Counts ($y=1$)\")\n                    )\\\n                    # Add probabilities P(y=1) column\n                    .with_columns(\n                        (pl.col(\"Counts ($y=1$)\") /\n                         (pl.col(\"Counts ($y=0$)\") + pl.col(\"Counts ($y=1$)\")) *\n                         100).alias(\"$P(y=1)\\,[\\%]$\")\n                    )\\\n                    # Keep solely row for which the bin range condition when grouping holds true\n                    .filter(pl.col(col) == True)\\\n                    # Drop col column\n                    .drop(col)\n                    # Add columns for the bin range limits\n                    .with_columns(\n                        pl.lit(x_bin_limit[i]).alias(\"Bin's $x_{\\mathrm{min}}$\"),\n                        pl.lit((x_bin_limit[i] + x_bin_limit[i + 1]) / 2).alias(\"Bin's $x_{\\mathrm{av}}$\"),\n                        pl.lit(x_bin_limit[i + 1]).alias(\"Bin's $x_{\\mathrm{max}}$\"),\n                    )\n                    # Rearrange columns\n                    .select(\n                        \"Bin's $x_{\\mathrm{min}}$\",\n                        \"Bin's $x_{\\mathrm{av}}$\",\n                        \"Bin's $x_{\\mathrm{max}}$\",\n                        \"Counts ($y=0$)\",\n                        \"Counts ($y=1$)\",\n                        \"$P(y=1)\\,[\\%]$\"\n                    )\n            )\n\n            df_dist = pl.concat([\n                df_dist,\n                df_dist_i\n            ],\n                how=\"vertical_relaxed\"\n            )\n\n            df_dist_plot = pl.concat([\n                df_dist_plot,\n                df_dist_i.select(pl.col(\"Bin's $x_{\\mathrm{min}}$\").alias(\"x\")).with_columns(pl.lit(0).alias(\"y\")),\n                df_dist_i.select(pl.col(\"Bin's $x_{\\mathrm{min}}$\").alias(\"x\"), pl.col(\"$P(y=1)\\,[\\%]$\").alias(\"y\")),\n                df_dist_i.select(pl.col(\"Bin's $x_{\\mathrm{max}}$\").alias(\"x\"), pl.col(\"$P(y=1)\\,[\\%]$\").alias(\"y\")),\n                df_dist_i.select(pl.col(\"Bin's $x_{\\mathrm{max}}$\").alias(\"x\")).with_columns(pl.lit(0).alias(\"y\"))\n            ],\n                how=\"vertical_relaxed\"\n            )\n    else:\n        # Polars dataframe with distributions of counts and probabilities\n        df_dist = (df.select(col, \"target\")\\\n            .groupby(col)\\\n            # Aggregate to compute counts\n            .agg((pl.col(\"target\") == 0).sum().alias(\"Counts ($y=0$)\"),\n                 (pl.col(\"target\") == 1).sum().alias(\"Counts ($y=1$)\"))\\\n            # Add probabilities P(y=1) column\n            .with_columns(\n            (pl.col(\"Counts ($y=1$)\") /\n             (pl.col(\"Counts ($y=0$)\") + pl.col(\"Counts ($y=1$)\")) *\n             100).alias(\"$P(y=1)\\,[\\%]$\"))\\\n            # Rename feature column\n            .rename({col: col_fancy.capitalize()})\n        )\n\n    # ----------------------------------- Plot ----------------------------------- #\n\n    # Display horizontal line\n    print()\n    display(Markdown('---'))\n    print()\n\n    # Initialise figure and axes\n    plt.figure(figsize=(6.4, 4.8))\n    ax = plt.axes()\n    plt.title(f\"Relation between probability of credit default and {col_fancy}\", pad=20)\n\n    if binning == True:\n        # Plot bar lines\n        ax.plot(\n            df_dist_plot[\"x\"],\n            df_dist_plot[\"y\"],\n            color=\"blue\",\n            linewidth=1.0,\n            zorder=2\n        )\n\n        # Fill bars\n        ax.fill_between(\n            x=df_dist_plot[\"x\"],\n            y1=np.zeros(len(df_dist_plot[\"x\"])),\n            y2=df_dist_plot[\"y\"],\n            facecolor=\"blue\",\n            edgecolor=\"black\",\n            hatch=\"////\",\n            linewidth=1.0,\n            alpha=0.5,\n            zorder=2\n        )\n\n        # Annotate probability values\n        for row in df_dist.rows(named=True):\n            ax.annotate(\n                text=\"${:.2f}\\,\\%$\".format(row[\"$P(y=1)\\,[\\%]$\"]),\n                xy=(row[\"Bin's $x_{\\mathrm{av}}$\"], \n                    row[\"$P(y=1)\\,[\\%]$\"]),\n                ha=\"left\",\n                va=\"bottom\",\n                size=9,\n                xytext=(-7.5, 2),\n                textcoords=\"offset points\",\n                rotation=70\n            )\n    else:\n        # Define bar plot   \n        sns.barplot(\n            ax=ax,\n            data=df_dist.to_pandas(),\n            x=col_fancy.capitalize(),\n            y=\"$P(y=1)\\,[\\%]$\",\n            color=\"blue\",\n            edgecolor=\"black\",\n            hatch=\"////\",\n            linewidth=1.0,\n            alpha=0.5,\n            width=0.200,\n            zorder=2\n        )\n\n        # Annotate probability values\n        for bar in ax.patches:\n            ax.annotate(\n                text=\"${:.2f}\\,\\%$\".format(bar.get_height()),\n                xy=(bar.get_x() + bar.get_width() / 2, \n                    bar.get_height()),\n                ha=\"left\",\n                va=\"bottom\",\n                size=9,\n                xytext=(-5, 2),\n                textcoords=\"offset points\",\n                rotation=70\n            )\n\n    # Define axes labels                                \n    ax.set_xlabel(col_fancy.capitalize(), fontdict={\"fontsize\": 10})\n    ax.set_ylabel(r\"$P(y=1)\\,[\\%]$\", fontdict={\"fontsize\": 10})\n\n    # Axes ticks\n    ax.tick_params(axis=\"x\", rotation=xticklabels_rotation)\n\n    # Axes limits\n    if binning == True:\n        ax.set_xlim(\n            left=min(df_dist_plot[\"x\"]),\n            right=max(df_dist_plot[\"x\"])\n        )\n    ax.set_ylim(bottom=0, top=1.125 * ax.get_ylim()[1])\n\n    # Enable axes' minor ticks\n    ax.minorticks_on()\n\n    # Define grid\n    ax.grid(\n        visible=True,\n        which=\"major\",\n        color=\"lightgray\",\n        linestyle=\"solid\",\n        linewidth=0.5\n    )\n    ax.grid(\n        visible=True,\n        which=\"minor\",\n        color=\"lightgray\",\n        linestyle=\"dotted\",\n        linewidth=0.5\n    )\n\n    # Show plot\n    plt.show()\n\n    # ------------------------- Display summary dataframe ------------------------ #\n    print()\n    display(df_dist.to_pandas()\\\n            .style\\\n            .format({\n                \"P_1\": \"{:.2f}\",\n                \"Bin's $x_{\\mathrm{min}}$\": \"{:.2f}\",\n                \"Bin's $x_{\\mathrm{av}}$\": \"{:.2f}\",\n                \"Bin's $x_{\\mathrm{max}}$\": \"{:.2f}\"\n            } if binning==True else\n            {\"P_1\": \"{:.2f}\"})\\\n            .set_caption(f\"{col_fancy.capitalize()} - counts and probabilities\")\\\n            .set_table_styles([\n                    # Column width\n                    {\"selector\": \"th.col_heading,td\",\n                     \"props\": [(\"width\", \"150px\")]\n                     },\n                    # Caption style\n                    {\"selector\": \"caption\",\n                     \"props\": [(\"font-size\", \"16px\"),\n                               (\"font-weight\", \"bold\"),\n                               (\"font-style\", \"italic\")]\n                     }\n                ]))","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:38.209251Z","iopub.execute_input":"2024-05-05T14:51:38.209838Z","iopub.status.idle":"2024-05-05T14:51:38.242436Z","shell.execute_reply.started":"2024-05-05T14:51:38.209805Z","shell.execute_reply":"2024-05-05T14:51:38.241651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Plot for contract's monthly payment\n\n# Rename column \"annuity_780A\" of the training and test dataframes to \"monthly_payment_A\"\ndt_data = {\n    batch: dt_data[batch].rename({\"annuity_780A\": \"monthly_payment_A\"})\n    for batch in batches\n}\n\n# Plot\ndist_plot(\n    df=dt_data[\"train\"],\n    col=\"monthly_payment_A\",\n    col_fancy=\"contract's monthly payment\",\n    binning=True,\n    N_bins=10,\n    xticklabels_rotation=0\n)\n\n# ---> Plot for applicant's age at decision time\n\n# Compute applicant's age at decision date and add it as column \"age_applicant_A\" of the\n# training and test dataframes\ndt_data = {\n    batch: dt_data[batch].with_columns(\n        (pl.col(\"date_decision\").sub(pl.col(\"birth_applicant_259D\"))\\\n         .alias(\"age_applicant_A\").dt.days() / 365.25).cast(pl.Float64)\n    ).drop(\"birth_applicant_259D\")\n    for batch in batches\n}\n\n# Plot\ndist_plot(\n    df=dt_data[\"train\"],\n    col=\"age_applicant_A\",\n    col_fancy=\"applicant's age at decision time\",\n    binning=True,\n    N_bins=20,\n    xticklabels_rotation=0\n)\n\n# ---> Applicant's employment time at decision date\n\n# Compute applicant's employment time (in years) and add it as column\n# \"employment_time_applicant_A\" of the training and test dataframes\ndt_data = {\n    batch: dt_data[batch].with_columns(\n        (pl.col(\"date_decision\").sub(pl.col(\"empl_employedfrom_applicant_271D\"))\\\n         .alias(\"employment_time_applicant_A\").dt.days() / 365.25).cast(pl.Float64)\n    ).drop(\"empl_employedfrom_applicant_271D\")\n    for batch in batches\n}\n\n# Plot\ndist_plot(\n    df=dt_data[\"train\"],\n    col=\"employment_time_applicant_A\",\n    col_fancy=\"applicant's employment time (in years) at decision date\",\n    binning=True,\n    N_bins=12,\n    xticklabels_rotation=0\n)\n\n# ---> Number of missing payments in the last 24 months (and current one)\n\n# Compute number of missing payments in the last 24 months (and current one) and add it\n# as column \"N_missing_payments_24_L\" of the training and test dataframes\ndt_data = {\n    batch: dt_data[batch].with_columns(\n        (25 - pl.col(\"cntpmts24_3658933L\"))\\\n         .alias(\"N_missing_payments_24_L\").cast(pl.Int64)\n    ).drop(\"cntpmts24_3658933L\")\n    for batch in batches\n}\n\n# Plot\ndist_plot(\n    df=dt_data[\"train\"],\n    col=\"N_missing_payments_24_L\",\n    col_fancy=\"number of missing payments in the last $24$ months (and current one)\",\n    binning=True,\n    N_bins=10,\n    xticklabels_rotation=0\n)\n\n# ---> Applicant's gender\n\n# Plot\ndist_plot(\n    df=dt_data[\"train\"],\n    col=\"sex_applicant_738L\",\n    col_fancy=\"applicant's gender\",\n    binning=False,\n    xticklabels_rotation=0\n)\n\n# ---> Number of persons of the same contract using the same mobile phone number \n\n# Redefine name and data type of column \"mobilephncnt_593L\"\ndt_data = {\n    batch: dt_data[batch].with_columns(\n        pl.col(\"mobilephncnt_593L\")\\\n        .alias(\"n_persons_same_phone_A\").cast(pl.Int64)\n    ).drop(\"mobilephncnt_593L\")\n    for batch in batches\n}\n\n# Plot\ndist_plot(\n    df=dt_data[\"train\"],\n    col=\"n_persons_same_phone_A\",\n    col_fancy=\"number of persons using the same phone number\",\n    binning=True,\n    N_bins=10,\n    xticklabels_rotation=0\n)\n\n# ---> Weekday of contract's decision date\n\n# Compute weekday of contract's decision date and add it as column\n# \"weekday_date_decision_A\"\n# For the sakeness of completeness, also compute month number and year\n# [NOTE: method dt.weekday() applied to a date returns the weekday in its number format.\n# The numbers are sorted according to the weekday in a week (1: \"Monday\", ..., 7:\n# \"Sunday\").]\ndt_data = {\n    batch: dt_data[batch].with_columns(\n        pl.col(\"date_decision\").dt.weekday().cast(pl.Int64)\\\n        .alias(\"weekday_date_decision_A\"),\n        pl.col(\"date_decision\").dt.month().cast(pl.Int64)\\\n        .alias(\"month_date_decision_A\"),\n        pl.col(\"date_decision\").dt.month().cast(pl.Int64)\\\n        .alias(\"year_date_decision_A\"),\n    )\n    for batch in batches\n}\n\n# Plot\ndist_plot(\n    df=dt_data[\"train\"],\n    col=\"weekday_date_decision_A\",\n    col_fancy=\"weekday of contract's decision date\",\n    binning=False,\n    xticklabels_rotation=45\n)\n\n# ---> Total number of payments made by the applicant\n\n# Rename column \"pmtnum_254L\" of the training and test dataframes to \"n_payments_A\" and\n# redefine type as int\ndt_data = {\n    batch: dt_data[batch].with_columns(\n        pl.col(\"pmtnum_254L\")\\\n        .alias(\"n_payments_A\").cast(pl.Int64)\n    ).drop(\"pmtnum_254L\")\n    for batch in batches\n}\n\n# Plot\ndist_plot(\n    df=dt_data[\"train\"],\n    col=\"n_payments_A\",\n    col_fancy=\"number of done payments\",\n    binning=True,\n    N_bins=12,\n    xticklabels_rotation=0\n)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:38.243662Z","iopub.execute_input":"2024-05-05T14:51:38.244069Z","iopub.status.idle":"2024-05-05T14:51:47.430255Z","shell.execute_reply.started":"2024-05-05T14:51:38.244037Z","shell.execute_reply":"2024-05-05T14:51:47.4293Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"One may conclude that\n\n* contracts with the highest monthly payments are less prone to be defaulted by the applicant. However, note that the number of counts at such range is quite low;\n\n* younger applicants have a higher probability of defaulting the credit contract than older ones. This could be justified by the fact that younger people usually do not have a so strong enconomic stability and also that they might take riskier decisions when asking for loans. The trend is such that the probability of credit default decreases with age until around $65$, from which it then starts to increase, however, never really attaining the very high values that are associated with the youth. This also evidences that the complications of aging after retirement might reduce economic stability;\n\n* the probability of credit default decreases with applicant's employment time, which may be explained by an associated increase in economic stability;\n\n* as one would expect, the probability of credit default increases with the number of missing payments;\n\n* males are approximately $1.50$ times more prone to default on their credit contract than females. This might evidence that males take riskier decisions when asking for loans.\n\n* The probability of credit default increases quite significantly with the number of persons (of the same contract) with the same phone number except at a very high values. However, this decreasing behaviour at high values is supported by very little amount of data;\n\n* The probability of credit default does depend on the weekday of contract's decisiond date: it increases from Sunday to Wednesday and decreases from Wednesday to Sunday. This might evidence that days of less work avail reflection, and, therefore, more careful decisions;\n\n* The probability of credit default does vary a lot with the number of done payments. But it is in average smaller for smaller numbers of payments and higher for higher ones. The smallest probability occurs at around $15$ payments.","metadata":{}},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Make the data ready for training\n\nIn order to train a model from [LightGBM](https://lightgbm.readthedocs.io/en/stable/), the respective dataset needs to be trasnformed into a [pandas](https://pandas.pydata.org/)' dataframes since the methods are currently not supported for [polars](https://pola.rs/)'.\n\nFurthermore, to better organise the data, each dataframe will be subdivided into three other: one with feature components $x^{(j)}$, one with labels $y$ and the third (a base) with auxiliary data for training and computation of performance metrics.\n\nDate columns are transformed into [ordinals](https://pandas.pydata.org/docs/reference/api/pandas.Timestamp.toordinal.html), allowing the respective entries to be handled like numbers. The ordinals are the number of days associated with the dates. The $1^{\\mathrm{st}}$ of January of year $1$ would be associated with the ordinal $1$. Such transformation would make huge sets of categories (the particular dates) to be represented by a single continuous variable.\n\nFor compatibility reasons, categories (unique values of object columns) of the feature test dataframe that do not pertain to the feature training dataframe are disregarded.\n\nAlso, since [SHAP](https://shap.readthedocs.io/en/latest/) cannot handle categorical columns, and Shapley values are to be computed later, these categorical columns are converted to integer ones using categories' codes. The respective code-category maps are added to the datasets. Furthermore, [SHAP](https://shap.readthedocs.io/en/latest/) also cannot handle boolean columns. These are converted to the integer type.","metadata":{}},{"cell_type":"code","source":"# ---> Shuffle training data (important when doing GridSearchCV to avoid geeting\n# \"biased\" splits)\n# [NOTE: a seed number is given for the sake of reproducibility: the shuffling result\n# would be the same in future calls.]\n# [NOTE: the rows are first ordered by \"case_id\" before shuffling to make sure that the\n# same dataframe is obtained from run to run, as the previous reading and the operations\n# on the data files might not result on the same order.]\ndt_data[\"train\"] = dt_data[\"train\"].sort(\"case_id\")\ndt_data[\"train\"] = dt_data[\"train\"].sample(fraction=1, shuffle=True, seed=42)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:47.431342Z","iopub.execute_input":"2024-05-05T14:51:47.431652Z","iopub.status.idle":"2024-05-05T14:51:51.232377Z","shell.execute_reply.started":"2024-05-05T14:51:47.431626Z","shell.execute_reply":"2024-05-05T14:51:51.231431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Feature columns\n\n# List of labels of the columns of the dataframes which are associated with features\n# [NOTE: features are such that the labels of the respective columns have lower case\n# characters in their whole extent except in the last position.]\ncols_x = []\nfor col in dt_data[\"train\"].columns:\n    if col[-1].isupper() and col[:-1].islower():\n        cols_x.append(col)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:51.233686Z","iopub.execute_input":"2024-05-05T14:51:51.234143Z","iopub.status.idle":"2024-05-05T14:51:51.24051Z","shell.execute_reply.started":"2024-05-05T14:51:51.234102Z","shell.execute_reply":"2024-05-05T14:51:51.239528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Define auxiliary functions\n\n# Function for converting object columns of pandas dataframes to category columns\n# [NOTE: if there is a column of the tuple of dataframes that is type \"object\", the\n# respective columns in the that and the other dataframes are converted to the type\n# category.]\n# [NOTE: a column of type \"object\" is such that it has entries of solely type \"str\" or\n# entries of multiple different types (e.g. \"NoneType\" and \"float\").]\n# [NOTE: category columns are columns whose entries pertain to a finite list of text\n# values.]\n# [NOTE: the category \"Unknown\" is added to the list of categories to support the case\n# in which validation and test dataframes have exclusive categories not pertaining to\n# the training dataframe - these exclusive categories which are not supported by the\n# trained model should be replaced by the \"Unknown\" category.]\ndef convert_cols_obj_to_cols_cat(*dfs):\n    # List of columns of the tuple of dataframes that are of type \"object\" in at least\n    # one of the dataframes\n    cols_object = list(set().union(*(df.select_dtypes(include=[\"object\"]).columns for df in dfs)))\n    # For each column of dtype \"object\"\n    for col in cols_object:\n        # For each dataframe of the tuple\n        for df in dfs:\n            # Convert current column to dtype \"category\"\n            df[col] = df[col].astype(\"category\")\n            # New categorical dtype whose categories correspond to the ones of the\n            # current column and the category \"Unknown\", being ordered\n            new_dtype = pd.CategoricalDtype(categories=df[col].cat.categories.to_list() +\n                                            [\"Unknown\"],\n                                            ordered=True)\n            # Assign new dtype to current column\n            df[col] = df[col].astype(new_dtype)\n    return dfs\n\n# Function for converting entries of date columns of a pandas dataframe to ordinals,\n# that is, the number of days, being day 1 of January of year 1 corresponding to the\n# number 1\n# [NOTE: this would make dates to be handled like numbers.]\ndef convert_cols_date_to_cols_ord(df):\n    # List of columns of date type\n    cols_date = df.select_dtypes(include=[\"datetime64\"]).columns\n    # For each column of date type\n    for col in cols_date:\n        # Convert each date in the current column to the respective ordianal\n        # [NOTE: if the entry is missing, let one associate a negative value (-1000) as\n        # ordinal.]\n        df[col] = df[col].apply(lambda x: x.toordinal() if not\n                                pd.isnull else -1000)\n    return df\n\n\n# Function for converting boolean columns of pandas dataframes to integer columns\ndef convert_cols_bool_to_cols_int(df):\n    # List of columns of the pandas dataframe that are of boolean type\n    cols_bool = df.select_dtypes(include=[\"bool\"]).columns\n    # For each column of type bool\n    for col in cols_bool:\n        # Convert column to int\n        df[col] = df[col].astype(\"int64\")\n    return df\n\n# Function for making categories of some pandas dataframe that do not pertain to a\n# reference one be replaced by the category \"Unknown\"\ndef make_cat_excl_unknown(df, df_ref):\n    # For each categorical column of the reference pandas dataframe\n    for col in df_ref.select_dtypes(include=[\"category\"]).columns:\n        # List of categories in the reference pandas dataframe\n        cat_ref = df_ref[col].cat.categories.to_list()\n        # List of categories in the pandas dataframe of interest\n        cat = df[col].cat.categories.to_list()\n        # List of common categories\n        cat_common = list(set(cat).intersection(cat_ref))\n        # List of exclusive categories\n        cat_exc = list(set(cat).difference(cat_common))\n        # New categorical dtype whose categories correspond to the common ones\n        new_dtype = pd.CategoricalDtype(categories=cat_common,\n                                        ordered=True)\n        # Replace current column's entries associated with exclusive categories as\n        # \"Unknown\"\n        df[col] = df[col].replace(to_replace=cat_exc, value=\"Unknown\")\n        # Assign the new dtype to the current column\n        df[col] = df[col].astype(new_dtype)\n    return df\n\n# Convert categorical columns to code ones while saving the respective mapping\ndef covert_cols_cat_to_cols_code(df):\n    # Dictionary dictionaries of category codes (keys) and respective category names\n    # (values) (a dictionary per feature categorical column)\n    # [NOTE: these may be used later to convert the codes to the category names through\n    # pandas Series' map method.]\n    dt_map = {}\n    \n    # For each categorical column of the pandas dataframe\n    for col in df.select_dtypes(include=[\"category\"]).columns:\n        \n        # Update dictionary of maps\n        dt_map.update({col:\n            dict(enumerate(df[col].cat.categories))\n        })\n        \n        # Convert categorical column to a code one\n        df[col] = df[col].cat.codes\n\n    return (df, dt_map)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:51.241618Z","iopub.execute_input":"2024-05-05T14:51:51.241898Z","iopub.status.idle":"2024-05-05T14:51:51.258524Z","shell.execute_reply.started":"2024-05-05T14:51:51.241874Z","shell.execute_reply":"2024-05-05T14:51:51.257624Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Define dictionaries of pandas dataframes associated with the training and test\n# data\ndt_data[\"train\"] = {\n    # Base dataframe with auxiliar data for the computation of the performance metrics\n    \"base\": (dt_data[\"train\"][[\"case_id\", \"date_decision\", \"WEEK_NUM\", \"target\", \"sample_weight\"]]\\\n             .rename({\"target\": \"y\"}).to_pandas()),\n    # Features' dataframe\n    \"x\": dt_data[\"train\"][cols_x].to_pandas(),\n    # Labels' dataframe\n    \"y\": dt_data[\"train\"][\"target\"].to_pandas(),\n}\ndt_data[\"test\"] = {\n    \"base\": dt_data[\"test\"][[\"case_id\", \"date_decision\", \"WEEK_NUM\"]].to_pandas(),\n    \"x\": dt_data[\"test\"][cols_x].to_pandas()\n}\n\n# Convert object columns of feature pandas dataframes to category columns and also add\n# the category \"Unknown\"\n# [NOTE: if there is a column of the tuple of dataframes that is type \"object\", the\n# respective columns in the that and the other dataframes are converted to the type\n# category.]\n(dt_data[\"train\"][\"x\"], dt_data[\"test\"][\"x\"]) = convert_cols_obj_to_cols_cat(\n    dt_data[\"train\"][\"x\"], dt_data[\"test\"][\"x\"]\n)\n\n# Convert date columns of feature pandas dataframes to ordinal ones\ndt_data[\"train\"][\"x\"] = convert_cols_date_to_cols_ord(dt_data[\"train\"][\"x\"])\ndt_data[\"test\"][\"x\"] = convert_cols_date_to_cols_ord(dt_data[\"test\"][\"x\"])\n\n# For compatibility reasons, make categories of the feature test dataframes which do not\n# pertain to the the training dataframe be replaced by the category \"Unknown\"\ndt_data[\"test\"][\"x\"] = make_cat_excl_unknown(\n    df=dt_data[\"test\"][\"x\"], df_ref=dt_data[\"train\"][\"x\"]\n)\n\n# Convert feature categorical columns to integer ones using categories' codes\n# [NOTE: category codes are the indices of the categories in column's categories array.]\n# [NOTE: such conversion is required to compute Shapley values when using SHAP. SHAP\n# cannot handle categorical columns.]\n(dt_data[\"train\"][\"x\"], dt_data[\"train\"][\"cat_map\"]) = covert_cols_cat_to_cols_code(dt_data[\"train\"][\"x\"])\n(dt_data[\"test\"][\"x\"], dt_data[\"test\"][\"cat_map\"]) = covert_cols_cat_to_cols_code(dt_data[\"test\"][\"x\"])\n\n# Convert feature boolean columns to integer ones\n# [NOTE: such conversion is required to compute Shapley values when using SHAP. SHAP\n# cannot handle boolean columns.]\ndt_data[\"train\"][\"x\"] = convert_cols_bool_to_cols_int(dt_data[\"train\"][\"x\"])\ndt_data[\"test\"][\"x\"] = convert_cols_bool_to_cols_int(dt_data[\"test\"][\"x\"])","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:51:51.259721Z","iopub.execute_input":"2024-05-05T14:51:51.26005Z","iopub.status.idle":"2024-05-05T14:52:50.379384Z","shell.execute_reply.started":"2024-05-05T14:51:51.259988Z","shell.execute_reply":"2024-05-05T14:52:50.378478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Theoretical background on LightGBM's Gradient Boosting method","metadata":{}},{"cell_type":"markdown","source":"\n\nIn this work it was decided to train a Gradient Boosting Decision Tree using the [LightGBM](https://lightgbm.readthedocs.io/en/stable/) package. This package issues a training method, [`lightgbm.train()`](https://lightgbm.readthedocs.io/en/latest/pythonapi/lightgbm.train.html), and [parameters](https://lightgbm.readthedocs.io/en/latest/Parameters.html).\n\n## Gradient boosting\n\nSection $10.10.2$ of the book [\"The Elements of Statistical Learning\"](https://hastie.su.domains/Papers/ESLII.pdf#page=378) by Trevor Hastie et al. is worth reading.\n\nThe [gradient boosting](https://en.wikipedia.org/wiki/Gradient_boosting) algorithm is in some way analogous to the one of  [gradient descent](https://en.wikipedia.org/wiki/Gradient_descent).\n\nLet $\\{x_i\\}$ and $y_i$ be the feature vector and label of the $i$-th training data point, with $i=1,\\,\\dots,\\,n$ where $n$ is the number of points. Let $\\hat{y}_i=f(\\{x_i\\})$ be the model estimate for the label associated with the feature vector $\\{x_i\\}$. And let $L(y_i,\\,\\hat{y}_i) = L(y_i,\\,f(\\{x_i\\}))$ be the loss function for the prediction of that label. The cost function would be the sum of all loss functions:\n\n$$C = \\sum_{i=1}^n L\\left(y_i,\\,f(\\{x_i\\})\\right)\\text{ .}$$\n\n### The gradient descent approach\n\nClassically, function $f$ would be described by a set of parameters $\\{\\theta\\}$, as $f_{\\{\\theta\\}}$, which would need to be fitted in order to minimise the cost function. The classical problem would then be to find the optimum set of parameters,\n\n$$\\{\\theta^*\\}=\\underset{\\{\\theta\\}}{\\mathrm{argmin}}\\sum_{i=1}^n L\\left(y_i,\\,f_{\\{\\theta\\}}(\\{x_i\\})\\right)\\text{ .}$$\n\nIn the gradient descent approach, an initial guess $\\{\\theta_0\\}$ would be given to the set of parameters, and these would be successively updated with negatively scaled gradients of the cost $C$ with respect to $\\{\\theta\\}$ at the previous value of the set of parameters, until convergence. The $m$-th iteration of the set of parameters would then be given by\n\n$$\\{\\theta_m\\} = \\{\\theta_{m-1}\\} - \\beta_m \\left\\{\\nabla_{\\{\\theta\\}}C\\right\\}(\\{\\theta_{m-1}\\})\\text{ ,}$$\n\nwhere $\\beta_m>0$ is the learning rate for that iteration.\n\n### The gradient boosting approach\n\nThe gradient boosting approach follows a particular motivation: let one consider for a moment that function $f$ is infinitely flexible, so that its evaluation at each data feature vector $\\{x_i\\}$ may be regarded as a variable. Let $\\{f\\}=[f(\\{x_1\\}),\\,\\dots,\\,f(\\{x_n\\})]^T$ be such vector of variables. The component $f^{(i)}=f(\\{x_i\\})$ would correspond to the variable associated with the $i$-th data point. The optimisation problem would then be to find the optimum vector of variables,\n\n$$\\{f^*\\} = \\underset{\\{f\\}}{\\mathrm{argmin}}\\sum_{i=1}^n L\\left(y_i,\\,f_i\\right)\\text{ .}$$\n\nBy analogy to the gradient descent approach, the variable $f^{(i)}$ would be regarded as a parameter and the set $\\{f\\}$ would be updated according to the recursion rule\n\n$$\n\\begin{aligned}\n\\{f_m\\} & = \\{f_{m-1}\\} - \\beta_m \\left\\{\\nabla_{\\{f\\}}C\\right\\}(\\{f_{m-1}\\})\\\\\n& = \\{f_{m-1}\\} - \\beta_m \\sum_{i=1}^n \\left[\\frac{\\partial\\,L\\left(y_i,\\,f\\right)}{\\partial\\, f}\\right]_{f=f_{m-1}^{(i)}}\\{e_i\\}\n\\text{ .}\n\\end{aligned}\n$$\n\nwhere $\\{e_i\\}$ is the versor in the $i$-th direction.\n\nFurthermore, one could be more \"ambitious\" and choose $\\beta_m$ as the learning rate that would minimise the cost function resulting from the above rule:\n\n$$\n\\begin{aligned}\n\\beta_m & = \\underset{\\beta}{\\mathrm{argmin}}\\sum_{i=1}^n L\\left(y_i,\\,f_{m,\\,\\beta}^{(i)}\\right)\\\\\n& = \\underset{\\beta}{\\mathrm{argmin}}\\sum_{i=1}^n L\\left(y_i,\\,f_{m-1}^{(i)} - \\beta \\sum_{i=1}^n \\left[\\frac{\\partial\\,L\\left(y_i,\\,f\\right)}{\\partial\\, f}\\right]_{f=f_{m-1}^{(i)}}\\right)\\text{ .}\n\\end{aligned}\n$$\n\nAnd by replacing $f^{(i)}$ by its definition, $f^{(i)}=f(\\{x_i\\})$, in the above equations, one gets\n\n$$\nf_m(\\{x_i\\}) = f_{m-1}(\\{x_i\\}) - \\beta_m \\left[\\frac{\\partial\\,L\\left(y_i,\\,f\\right)}{\\partial\\, f}\\right]_{f=f_{m-1}(\\{x_i\\})}\n$$\n\nand\n\n$$\n\\beta_m = \\underset{\\beta}{\\mathrm{argmin}}\\sum_{i=1}^n L\\left(y_i,\\,f_{m-1}(\\{x_i\\}) - \\beta \\left[\\frac{\\partial\\,L\\left(y_i,\\,f\\right)}{\\partial\\, f}\\right]_{f=f_{m-1}(\\{x_i\\})}\\right)\\text{ .}\n$$\n\nHowever, a bearable (in the sense of not being too complex) function cannot be infinitely flexible. It needs to have some structure, and this structure would support the whole feature vector $\\{x\\}$ domain. It would be impossible to define some bearable function that could independently adjust itselft to each of the data points. \n\nThe pragmatic approach would be to define a simple function $f_0(\\{x\\})$ as an initial estimate for $f(\\{x\\})$, and then use an update rule of the form \n\n$$ f_m(\\{x\\}) = f_{m-1}(\\{x\\}) + \\beta_m\\,g_m(\\{x\\})\\text{ ,}$$\n\nwhere $g_m(\\{x\\})$ would be another simple function - termed \"weak learner\" - that better reproduces the symmetric of the derivatives of the loss function derivatives at the training points, $-\\left[\\frac{\\partial\\,L\\left(y_i,\\,f\\right)}{\\partial\\, f}\\right]_{f=f_{m-1}(\\{x_i\\})}$. Let\n\n$$r_{m,i}=-\\left[\\frac{\\partial\\,L\\left(y_i,\\,f\\right)}{\\partial\\, f}\\right]_{f=f_{m-1}(\\{x_i\\})}\\text{ ,}$$\n\nbe the so-called \"pseudo-residual\" associated with the $m$-th iteration and $i$-th data point. The optimum weak learner $g_m(\\{x\\})$ would be the one that better fits into the points $\\left\\{\\left(\\{x_i\\},\\,r_{m,i}\\right)\\right\\}_{i=1}^n$. In the case of Gradient Boosting Decision Tree, the weak learners would correspond to decision trees of not so large depth.\n\nFor the new updating rule, the optimum learning rate would correspond to\n\n$$\n\\beta_m = \\underset{\\beta}{\\mathrm{argmin}}\\sum_{i=1}^n L\\left(y_i,\\,f_{m-1}(\\{x_i\\}) + \\beta \\, g_m(\\{x_i\\})\\right)\\text{ .}\n$$\n\n\nThis is the gradient boosting approach that encompasses the idea of gradient descent and the incremental usage of weak learners to obtain a \"strong learner\" - a well behaved model that accurately predicts the training data points as well as others out of the domain of the former.\n","metadata":{}},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Tune and fit LightGBM's Gradient Boosting Decision Tree with GridSearchCV","metadata":{}},{"cell_type":"markdown","source":"\nLightGBM's Gradient Boosting Decision Tree has several hyperparameters. Two of them are the [maximum depth](https://lightgbm.readthedocs.io/en/latest/Parameters.html#max_depth) of the weak decision trees, $\\mathrm{d}_{\\mathrm{max}}$, and the [number of such trees](https://lightgbm.readthedocs.io/en/latest/Parameters.html#num_iterations), $N_{\\mathrm{trees}}$. The default value for the maximum depth is `max_depth=-1`, which according to the documentation means \"no limit\" (as it is non-positive). And the default value of the number of trees is `n_estimators=100`.\n\nThe hyperparameters `max_depth` and`n_estimators` may be tuned by using grid-search together with cross-validation (CV). skit-learn's [GridSearchCV](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html#sklearn.model_selection.GridSearchCV) may be herein applied for such purpose.\n\nThe grid-search method consists in defining a matrix of hyperparameters $[\\Gamma]$ with its rows and columns associated with the hyperparameters' values and variables respectively. The distribution of values would be \"cartesian\", that is, each value of a hyperparamenter variable would combine with every single one of the other hyperparameter variables. Consider, for instance, that there are $3$ hyperparameter variables, $\\gamma^{(1)}$, $\\gamma^{(2)}$ and $\\gamma^{(3)}$, and that these may take values $\\{\\gamma_1^{(1)}\\}$, $\\{\\gamma_1^{(2)},\\,\\gamma_2^{(2)},\\,\\gamma_3^{(2)}\\}$ and $\\{\\gamma_1^{(3)},\\,\\gamma_2^{(3)}\\}$, respectively. The hyperparameter matrix would have $1 \\times 3 \\times 2 = 6$ rows and $3$ columns:\n\n\\begin{equation*}\n[\\Gamma] =\n\\begin{bmatrix}\n\\gamma_1^{(1)} & \\gamma_1^{(2)} & \\gamma_1^{(3)} \\\\\n\\gamma_1^{(1)} & \\gamma_1^{(2)} & \\gamma_2^{(3)} \\\\\n\\gamma_1^{(1)} & \\gamma_2^{(2)} & \\gamma_1^{(3)} \\\\\n\\gamma_1^{(1)} & \\gamma_2^{(2)} & \\gamma_2^{(3)} \\\\\n\\gamma_1^{(1)} & \\gamma_3^{(2)} & \\gamma_1^{(3)} \\\\\n\\gamma_1^{(1)} & \\gamma_3^{(2)} & \\gamma_2^{(3)} \\\\\n\\end{bmatrix}\\text{ .}\n\\end{equation*}\n\nEach combination of hyperparameter values, that is, each row of $[\\Gamma]$, would need to be tried to obtain a model from the training data and some metric score evaluated. This may be done through $K$-Fold Cross-Validation. In this approach, the training set is divided into $k$ folds of identical size. Then, for each combination of hyperparameter values, the model is trained using $k-1$ folds, validated using the unpicked fold and a performance score on the validation fold is computed. And still for each combination the process is repeated until all $k$ folds are used as validation ones. The combination of hyperparameter values for which the mean value of the validation score is the highest would be the one to be picked. Note that the mean value of the validation score for some combination would correspond to the artithmetic mean of the performance scores on each of the $k$ validation folds using the models trained on the $k-1$ other and with this very combination of hyperparamater values.\n\nA typical performance score in classification is accuracy. However, the one considered in this competition for the evaluation of the submitted restults directly depends on the area under the ROC Curve (that is, AUC), and, therefore, it was decided to herein take it as GridSearchCV's validation score.\n\nRegarding GridSearchCV's parameters, one should mention:\n\n* If GridSearchCV's `cv` parameter is set to an integer, let this be $k$, the training set is split into $k$ [stractified folds](https://scikit-learn.org/stable/modules/cross_validation.html#stratified-k-fold). Stractified folds preserve the proportion between counts of the differents classes in them, that is, they have the same proportion as the original training set. By default, `cv = 5`. Note that if an integer is considered, shuffling will not be formely applied to to the data before splitting it.\n\n* GridSearchCV's `scoring` corresponds to the strategy used for evaluating the performance of the cross-validated model in each training and validation fold. For the [case of classifiers](https://datascience.stackexchange.com/a/94263/158543) one has by default `scoring=\"accuracy\"`.\n\n* GridSearchCV's `return_train_score` parameter is a boolean that defines if the training score (besides the validation one) of each trial is to be reported in a returned GridSearchCV's `cv_results_` attribute. By setting `return_train_score=True`, the training score is reported. By default, `return_train_score = False`.\n\n* GridSearchCV's `n_jobs` parameter corresponds to the number of jobs to run in parallel. By setting `n_jobs = -1`, [all CPUs of the kaggle kernel is used in parallel for doing the computations](https://scikit-learn.org/stable/glossary.html#term-n_jobs). By default, `n_jobs = 1` except in a so-called \"[joblib.parallel_backend](https://joblib.readthedocs.io/en/latest/generated/joblib.parallel_backend.html#joblib.parallel_backend)\" context.\n\n* GridSearchCV's `refit` parameter is a boolean that defines if at the end of the procedure the best found estimator should be refitted to the whole training set. It is then allocated in GridSearchCV's `best_estimator_` attribute. By default, `refit = True`.\n\n* GridSearchCV's `verbose` parameter is an integer that controls the verbosity: the higher, the more messages. The default value is `verbose = 0`.","metadata":{}},{"cell_type":"code","source":"# ---> LightGBM's Gradient Boosting Decision Tree (the estimator template)\n\n# Dictionary of parameters for training the LightGBM estimator\nparams_lgb = {\n    # Gradient Boosting type\n    # [NOTE: possible values:\n    #  * \"gbdt\", standing for \"Gradient Boosting Decision Tree\",\n    #  * \"dart\", standing for \"Dropouts meet multiple Additive Regression Trees\",\n    #  * \"rf\", standing for \"Random Forest\",\n    # .],\n    \"boosting_type\": \"gbdt\",\n    # Learning objective\n    \"objective\": \"binary\",\n    # List of metrics to be evaluated while training\n    # [NOTE: \"auc\" stands for \"Area Under the Curve\" - this curve is the \"Receiver\n    # Operating Characterisitc\" (ROC).]\n    \"metric\": \"auc\",\n    # Maximum number of leave nodes in a weak tree learner\n    \"num_leaves\": 31,\n    # Learning rate\n    \"learning_rate\": 0.05,\n    # Fraction of feature components to use in each Gradient Boosting iteration\n    # [NOTE: LightGBM will randomly select a subset of feature components on each\n    # iteration (tree). The dimensionality of this subset is a fraction feature_fraction\n    # of the dimensionality of the whole feature vector.]\n    # [NOTE: feature_fraction must pertain to the interval ]0, 1].]\n    # [NOTE: the default value is 1.]\n    # [NOTE: feature_fraction may be used to speed up training and to avoid getting an\n    # overfitting model.]\n    \"feature_fraction\": 1,\n    # Fraction of data to use in each Gradient Boosting iteration\n    # [NOTE: LightGBM will randomly select a subset of training points on each\n    # iteration (tree). The cardinality of this subset is a fraction bagging_fraction\n    # of the cardinality of the whole training set.]\n    # [NOTE: to consider bagging, the bagging frequency parameter (bagging_freq) must be\n    # set to a non-null value (by default it is null).]\n    # [NOTE: bagging_fraction must pertain to the interval ]0, 1].]\n    # [NOTE: the default value is 1.]\n    # [NOTE: bagging_fraction may be used to speed up training and to avoid getting an\n    # overfitting model.]\n    \"bagging_fraction\": 1,\n    # Bagging frequency\n    # [NOTE: A bagging frequency of k means that bagging is done at each k iterations\n    # (trees) and that resultant subset of training points is used in the current\n    # iteration and in the next k-1.]\n    # [NOTE: the default value is 0.]\n    \"bagging_freq\": 0,\n    # Maximum number of weak learners (the same as num_iterations)\n    \"n_estimators\": 1200,\n    # Max depth of the weak tree learners\n    # [NOTE: default value is -1. Any value less or equal to 0 means no limit.]\n    \"max_depth\": 3,\n    # Number of parallel threads\n    # [NOTE: -1 means using all threads.]\n    \"n_jobs\": -1,\n    # Random state seed number, for the sake of reproducibility\n    # [NOTE: actually even with a given seed number, reproducibility is not guaranteed.\n    # Using multiple threads (as defined by n_jobs) or using GPU breaks reproducibility\n    # since each processor treats its own sub-batch data, and the way the data is\n    # distributed to the processers is not deterministic. It has been reported that\n    # bagging with multi-threading may be one of the principal causes for\n    # non-reproducibility. Fortunately, in the latest versions (starting at v3.0.0rc1)\n    # of LightGBM, bagging was redefined so that its behaviour does not depend on the\n    # number of threads. See FAQ: https://lightgbm.readthedocs.io/en/latest/FAQ.html]\n    \"random_state\": 42,\n    # Type of device to perform the computations (\"cpu\", \"gpu\" or \"cuda\")\n    # [NOTE: to use a GPU in kaggle, this needs to be chosen through the button\n    # \"Accelerator\" in the panel \"Session options\" at the right-hand side. A device such\n    # as \"GPU P100\" could be chosen.]\n    \"device_type\": \"gpu\",\n    # Use double precision (64-bit float point) instead of single (32-bit) in the case\n    # of a GPU device being picked\n    # [NOTE: this may help to get reproducible results (though without guarantees).\n    # However it may slow down training.]\n    \"gpu_use_dp\": True,\n    # Verbosity level of LightGBM\n    # [NOTE: \"-1\" means solely \"fatal errors\".]\n    \"verbose\": -1,\n}\n\n# LightGBM's estimator\nlgb_estimator = lgb.LGBMClassifier(**params_lgb)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:52:50.380665Z","iopub.execute_input":"2024-05-05T14:52:50.380965Z","iopub.status.idle":"2024-05-05T14:52:50.390997Z","shell.execute_reply.started":"2024-05-05T14:52:50.380942Z","shell.execute_reply":"2024-05-05T14:52:50.390139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Perform GridSearchCV on LightGBM's estimator\n\n# Dictionary of parameters' lists of values to be tried in GridSearchCV\n# [NOTE: the keys of this dictionary should coincide with the names of the estimator\n# arguments that represent such hyperparamenters.]\nparam_grid_lgb = {\n    \"n_estimators\": [600, 800, 1000],\n    \"max_depth\": [8, 9, 10]\n}\n\n# Create a GridSearchCV object\ngs_lgb = GridSearchCV(\n    estimator=lgb_estimator, \n    param_grid=param_grid_lgb,\n    cv=5,\n    scoring=\"roc_auc\",\n    return_train_score=True,\n    n_jobs=1,\n    refit=True,\n    verbose=0\n)\n\n# Initial instant of time [s]\nt_i = time.perf_counter()\n\n# Perform GridSearchCV on the training data\n# [NOTE: to perform a balanced training a sample weight array is issued. Note that\n# GridSearchCV uses it only in the cost function and not in the scorer. Anyway, the ROC\n# curve is actually invariable to class weights. Sample weighting could then in this\n# case be disregarded in scoring.]\ngs_lgb.fit(\n    X=dt_data[\"train\"][\"x\"],\n    y=dt_data[\"train\"][\"y\"],\n    sample_weight=dt_data[\"train\"][\"base\"][\"sample_weight\"]\n)\n\n# Final instant of time [s]\nt_f = time.perf_counter()\n\n# Print elapsed time\nprint()\nprint(f\"Running time: {(t_f - t_i)/3600:.2f} h\")","metadata":{"execution":{"iopub.status.busy":"2024-05-05T14:52:50.392117Z","iopub.execute_input":"2024-05-05T14:52:50.392374Z","iopub.status.idle":"2024-05-05T16:26:00.981114Z","shell.execute_reply.started":"2024-05-05T14:52:50.392352Z","shell.execute_reply":"2024-05-05T16:26:00.980345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Display summary on results of the performed GridSearchCV\n\n# Get cross-validation results\ncv_results_lgb = pd.DataFrame(gs_lgb.cv_results_)\n\n# Replace \"score\" substring in the column names by \"auc\"\ncv_results_lgb.columns = cv_results_lgb.columns.str.replace(\"score\", \"auc\")\n\n# Save to CSV file\ncv_results_lgb.to_csv(\"cv_results_lgb.csv\", index=True)\n\n# Display mean training and validation AUCs as well as respective standard devation for\n# each grid point\nprint()\ndisplay(cv_results_lgb[cv_results_lgb.filter(like=\"param_\", axis=\"columns\").columns.values.tolist() + \n                   [\"mean_train_auc\", \"std_train_auc\", \"mean_test_auc\", \"std_test_auc\"]]\\\n        .style.set_caption(\"Mean AUCs and respective standard devations\")\\\n        .set_table_styles(\n            [\n                # Column width\n                {\"selector\": \"th.col_heading,td\",\n                 \"props\": [(\"width\", \"100px\")]\n                 },\n                # Caption style\n                {\"selector\": \"caption\",\n                 \"props\": [(\"font-size\", \"16px\"),\n                           (\"font-weight\", \"bold\"),\n                           (\"font-style\", \"italic\")]\n                 }\n            ]))\n\n# Tuple of cross-validation results associated with best training and best validation\n# mean AUCs\n(cv_results_best_train_auc_lgb, cv_results_best_valid_auc_lgb) = (\n    pd.DataFrame(\n        data=(cv_results_lgb[cv_results_lgb.filter(like=\"param_\", axis=\"columns\").columns.values.tolist() + \n                         [\"mean_train_auc\", \"std_train_auc\", \"mean_test_auc\", \"std_test_auc\"]]\\\n              .loc[cv_results_lgb[f\"mean_{batch}_auc\"].idxmax()].to_dict()),\n        index=[cv_results_lgb[f\"mean_{batch}_auc\"].idxmax()]\n    ) for batch in [\"train\", \"test\"]\n)\n\n# Display\nfor (df, batch_fancy) in ((cv_results_best_train_auc_lgb, \"training\"),\n                          (cv_results_best_valid_auc_lgb, \"validation\")):\n    print()\n    display(df\\\n        .style.set_caption(f\"Mean AUCs and respective standard devations at best {batch_fancy} mean AUC\")\\\n        .set_table_styles(\n            [\n                # Column width\n                {\"selector\": \"th.col_heading,td\",\n                 \"props\": [(\"width\", \"100px\")]\n                 },\n                # Caption style\n                {\"selector\": \"caption\",\n                 \"props\": [(\"font-size\", \"16px\"),\n                           (\"font-weight\", \"bold\"),\n                           (\"font-style\", \"italic\")]\n                 }\n            ])\n    )","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:39:17.795723Z","iopub.execute_input":"2024-05-05T16:39:17.796092Z","iopub.status.idle":"2024-05-05T16:39:17.830964Z","shell.execute_reply.started":"2024-05-05T16:39:17.796062Z","shell.execute_reply":"2024-05-05T16:39:17.830086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Visualise dependence of the mean training and validation AUCs on the respective\n# hyperparameters\n\n# Array of unique n_estimators values\nn_estimators_lgb = cv_results_lgb[\"param_n_estimators\"].unique()\n# Array of unique max_depth values\nmax_depth_lgb = cv_results_lgb[\"param_max_depth\"].unique()\n# Array of arrays of mean_test_auc values (first dimension (rows) associated with\n# n_estimators and second (columns) with max_depth)\nmean_test_auc_lgb = (cv_results_lgb.groupby(\"param_n_estimators\")\\\n                 .agg(list)[\"mean_test_auc\"].to_list())\n\n# Initialise figure and axes\nplt.figure(figsize=(6.4, 4.8))\nax = plt.axes()\nplt.title(\n    r\"Dependence of $\\mathrm{AUC}_{\\mathrm{valid},\\,\\mathrm{av}}$ on $N_{\\mathrm{trees}}$ and $\\mathrm{d}_{\\mathrm{max}}$\",\n    pad=20\n)\n\n# Arrays of coordinates\nX = n_estimators_lgb\nY = max_depth_lgb\nZ = np.transpose(mean_test_auc_lgb)\n\n# Unfilled contour plot\ncs = ax.contour(X, Y, Z, levels=10, colors=\"white\", linestyles=\"solid\", linewidths=0.5)\n\n# Filled contour plot\ncsf = ax.contourf(X, Y, Z, levels=10, cmap=\"plasma\")\n\n# Define contour line labels\nax.clabel(CS=cs, fmt=\"%1.3f\", inline=True, fontsize=8)\n\n# Define axes labels                                \nax.set_xlabel(r\"Number of trees, $N_{\\mathrm{trees}}$\", fontdict={\"fontsize\": 10})\nax.set_ylabel(r\"Maximum tree depth, $\\mathrm{d}_{\\mathrm{max}}$\", fontdict={\"fontsize\": 10})\n\n# Colour bar\ncbar = plt.colorbar(mappable=csf, label=\"$\\mathrm{AUC}_{\\mathrm{valid},\\,\\mathrm{av}}$\", ax=ax)\ncbar.add_lines(cs)\n\n# Show plot\nplt.show() ","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:55:08.350394Z","iopub.execute_input":"2024-05-05T16:55:08.351425Z","iopub.status.idle":"2024-05-05T16:55:08.90817Z","shell.execute_reply.started":"2024-05-05T16:55:08.351389Z","shell.execute_reply":"2024-05-05T16:55:08.907264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For the considered range of hyperparameters values, the average validation $\\mathrm{AUC}$ actually does not change too much. The best value was obtained for a maximum tree depth of $\\mathrm{d}_{\\mathrm{max}}=8$ of and a number of trees of $N_{\\mathrm{trees}}=600$. But the best training $\\mathrm{AUC}$ was actually obtained for the greatest hyperparameters values, $\\mathrm{d}_{\\mathrm{max}}=10$ and $N_{\\mathrm{trees}}=1000$. This shows that overfitting occurs in that extreme regime.","metadata":{}},{"cell_type":"code","source":"# ---> XGBoost的梯度提升决策树（估计器模板）\n\n# 用于训练XGBoost估计器的参数字典\nparams_xgb = {\n    \"boosting_type\": \"gbdt\",\n    \"objective\": \"binary\",\n    \"metric\": \"auc\",\n    \"num_leaves\": 31,\n    \"learning_rate\": 0.05,\n    \"feature_fraction\": 1,\n    \"bagging_fraction\": 1,\n    \"bagging_freq\": 0,\n    \"n_estimators\": 1200,\n    \"max_depth\": 3,\n    \"n_jobs\": -1,\n    \"random_state\": 42,\n    \"device_type\": \"gpu\",\n    \"gpu_use_dp\": True,\n    \"verbose\": -1,\n}\n\n# XGBoost估计器\nxgb_estimator = xgb.XGBClassifier(**params_xgb)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:45:01.043376Z","iopub.execute_input":"2024-05-05T16:45:01.044114Z","iopub.status.idle":"2024-05-05T16:45:01.050219Z","shell.execute_reply.started":"2024-05-05T16:45:01.044083Z","shell.execute_reply":"2024-05-05T16:45:01.049205Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Perform GridSearchCV on LightGBM's estimator\n\n# Dictionary of parameters' lists of values to be tried in GridSearchCV\n# [NOTE: the keys of this dictionary should coincide with the names of the estimator\n# arguments that represent such hyperparamenters.]\nparam_grid_xgb = {\n    \"n_estimators\": [600, 800, 1000],\n    \"max_depth\": [8, 9, 10]\n}\n\n# Create a GridSearchCV object\ngs_xgb = GridSearchCV(\n    estimator=xgb_estimator,\n    param_grid=param_grid_xgb,\n    cv=5,  # 交叉验证的折数\n    scoring=\"roc_auc\",  # 评分标准\n    n_jobs=1,  # 并行工作的进程数\n    refit=True,  # 使用整个训练集对找到的最佳模型进行拟合\n    verbose=0\n)\n\n# Initial instant of time [s]\nt_i = time.perf_counter()\n\n# Perform GridSearchCV on the training data\n# [NOTE: to perform a balanced training a sample weight array is issued. Note that\n# GridSearchCV uses it only in the cost function and not in the scorer. Anyway, the ROC\n# curve is actually invariable to class weights. Sample weighting could then in this\n# case be disregarded in scoring.]\ngs_xgb.fit(\n    X=dt_data[\"train\"][\"x\"],\n    y=dt_data[\"train\"][\"y\"],\n    sample_weight=dt_data[\"train\"][\"base\"][\"sample_weight\"]\n)\n\n# Final instant of time [s]\nt_f = time.perf_counter()\n\n# Print elapsed time\nprint()\nprint(f\"Running time: {(t_f - t_i)/3600:.2f} h\")","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:56:06.814206Z","iopub.execute_input":"2024-05-05T16:56:06.814928Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Model performance assessment on the whole training data","metadata":{}},{"cell_type":"markdown","source":"## Confusion matrix\n\nLet one now compute the confusion matrix of the best estimator for the whole original training set. Note that since GridSearchCV's `refit` parameter was set to `True`, the best estimator was refitted to the original training set after Grid-search Cross-Validation being performed, corresponding to GridSearchCV's `best_estimator_` attribute.","metadata":{}},{"cell_type":"code","source":"# ---> Get best estimator\n# Best estimator (which has been refitted to the whole training data)\nbest_estimator_lgb = gs_lgb.best_estimator_","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.105261Z","iopub.status.idle":"2024-05-05T16:26:01.105626Z","shell.execute_reply.started":"2024-05-05T16:26:01.105435Z","shell.execute_reply":"2024-05-05T16:26:01.105449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Compute predicted probabilities and predicted labels for the whole original\n# training set\n          \n# Add predicted probabilities for label corresponding to 1 (credit default case)\n# to base pandas dataframe of the dataset dictionary\n# [NOTE: the method predict_proba returns a numpy.ndarray of how many columns there are\n# classes. The columns are ordered according to the attribute classes_ of the\n# classifier.]\ndt_data[\"train\"][\"base\"][\"P_pred_lgb\"] = best_estimator_lgb.predict_proba(\n    X=dt_data[\"train\"][\"x\"]\n)[:, 1]\n\n# Add predicted labels for the case of a threshold of 0.5 to base pandas dataframe of\n# the dataset dictionary\ndt_data[\"train\"][\"base\"][\"y_pred_lgb\"] = ((dt_data[\"train\"][\"base\"][\"P_pred\"] >= 0.5)\\\n                                      .astype(dtype=\"int32\"))","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.107007Z","iopub.status.idle":"2024-05-05T16:26:01.107339Z","shell.execute_reply.started":"2024-05-05T16:26:01.107177Z","shell.execute_reply":"2024-05-05T16:26:01.10719Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Compute and plot confusion matrix for the whole original training set\n\n# Confusion matrix associated with the predictions\ncm = confusion_matrix(\n    y_true=dt_data[\"train\"][\"y\"],\n    y_pred=dt_data[\"train\"][\"base\"][\"y_pred_lgb\"]\n)\n\n# Plot confusion matrix\ndisp = ConfusionMatrixDisplay(cm).plot(cmap=\"cividis\")\ndisp.ax_.set_title(\"LGBM Confusion matrix\", pad=20);","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.108532Z","iopub.status.idle":"2024-05-05T16:26:01.108867Z","shell.execute_reply.started":"2024-05-05T16:26:01.108701Z","shell.execute_reply":"2024-05-05T16:26:01.108714Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The $(i,\\,j)$-the entry of the confusion matrix corresponds to the occurrence of the frequency of predictions of the $j$-th class when the $i$-th class is the actual (correct) one. The confusion matrix of a good estimator would have its diagonal values much greater than the off-diagonal ones. And that's indeed the case of the present model.\n\nThe obtained number false positives ($\\mathrm{FP}$) is remarkably high - much higher (by a factor of $9.66$) than the number of true positives ($\\mathrm{TP}$). However, note that the number of actual negatives is much higher (by a factor of $30.81$) than the number of actual positives - the number of samples incorrectly labelled as positives would then be arguably higher when considering balancing sample weights in the loss functions.","metadata":{}},{"cell_type":"markdown","source":"## Precisions, recalls, $F_1$-score and accuracy\n\nOne should now compute performance metrics of the best estimator on the whole original training set, namely precisions, recalls, $F_1$ scores and accuracy.\n\nLet $f_{ij}$ be the $(i,\\,j)$ entry of the confusion matrix.\n\nThe precision of the estimator on the prediction of the $i$-th class corresponds to the ratio between the number of correct predictions of the $i$-th class and the total number of predictions of that class:\n\n\\begin{equation}\n\\mathrm{precision}_i = \\frac{f_{ii}}{\\sum_{j}f_{ji}} = \\frac{\\mathrm{TP}_i}{\\mathrm{TP}_i + \\mathrm{FP}_i}\\text{ .}\n\\label{eq:precisioni}\n\\end{equation}\n\nThe $i$-th precision is then the ratio between the $i$-th diagonal entry and the sum of the $i$-th column's entries of the confusion matrix. It may be understood as the ratio between true positives ($\\mathrm{TP}_i$) and all positives ($\\mathrm{TP}_i + \\mathrm{FP}_i$) of that class. Herein, a \"positive\" of the $i$-th class corresponds to its prediction regardless of being correct (\"true positive\") or not (\"false positive\"). In binary classification, usually solely precision of the positive class ($+1$) is considered:\n\n\n\\begin{equation}\n\\mathrm{precision}:= \\mathrm{precision}_{+}= \\frac{\\mathrm{TP}_+}{\\mathrm{TP}_+ + \\mathrm{FP}_+} = \\frac{\\mathrm{TP}}{\\mathrm{TP} + \\mathrm{FP}}\\text{ .}\n\\label{eq:precision}\n\\end{equation}\n\nA value of $1$ would mean no false positives, though, not necessarily no false negatives. In the current problem, no false positives would correspond to no falsely predicted credit defaults.\n\nThe recall of the estimator on the prediction of the $i$-th class corresponds to the ratio between the number of correct predictions of the $i$-th class and the total number of actual occurrences of that class:\n\n\\begin{equation}\n\\mathrm{recall}_i = \\frac{f_{ii}}{\\sum_{j}f_{ij}} = \\frac{\\mathrm{TP}_i}{\\mathrm{TP}_i + \\mathrm{FN}_i}\\text{ .}\n\\label{eq:recalli}\n\\end{equation}\n\nThe $i$-th recall is then the ratio between the $i$-th diagonal entry and the sum of the $i$-th row's entries of the confusion matrix. It may be understood as the ratio between true positives ($\\mathrm{TP}_i$) and all actual occurrences ($\\mathrm{TP}_i + \\mathrm{FN}_i$). The term \"recall\" comes from ability of the estimator that the quantity represents: recalling (remembering) actual occurrences of some class. A value of $1$ would mean no false negatives, though, not necessarily no false positives. In binary classification, as in the case of precision, usually solely recall of the positive class ($+1$) is considered:\n\n\\begin{equation}\n\\mathrm{recall} := \\mathrm{recall}_+ = \\frac{\\mathrm{TP}_+}{\\mathrm{TP}_+ + \\mathrm{FN}_+} = \\frac{\\mathrm{TP}}{\\mathrm{TP} + \\mathrm{FN}}\\text{ .}\n\\label{eq:recall}\n\\end{equation}\n\nIn the current problem no false negatives would mean no falsely predicted credit non-defaults.\n\nThe $F_1$ score of the estimator on the prediction of the $i$-th class is simply the harmonic mean (the reciprocal of the arithmetic mean of the reciprocals) of precision and recall of the $i$-th class:\n\n\\begin{equation}\nF_{1,i} = \\frac{2}{\\mathrm{precision}_i^{-1} + \\mathrm{recall}_i^{-1}}\\text{ .}\n\\label{eq:F1i}\n\\end{equation}\n\nAnd again, in the case of binary classification, one usually just computes the $F_1$ score of the positive class:\n\n\\begin{equation}\nF_{1} := F_{1,+} = \\frac{2}{\\mathrm{precision}_+^{-1} + \\mathrm{recall}_+^{-1}} = \\frac{2}{\\mathrm{precision}^{-1} + \\mathrm{recall}^{-1}}\\text{ .}\n\\label{eq:F1}\n\\end{equation}\n\n\nFinally, the accuracy of the estimator corresponds to the ratio between the number of correct predictions and the total number of predictions:\n\n\\begin{equation}\n\\mathrm{accuracy} = \\frac{\\sum_{i}f_{ii}}{\\sum_{i}\\sum_{j}f_{ij}} = \\frac{\\sum_i \\left(\\mathrm{TP}_i + \\mathrm{TN}_i\\right)}{\\sum_i \\left(\\mathrm{TP}_i + \\mathrm{TN}_i + \\mathrm{FP}_i + \\mathrm{FN}_i\\right)}\\text{ .}\n\\label{eq:accuracy}\n\\end{equation}\n\nIt is then the ratio between the sum of all diagonal entries and the sum of all entries of the confusion matrix. And it may be understood as the ratio between the sum of trues positives of all classes ($\\sum_i \\mathrm{TP}_i$) and the sum of true positives and the average of falses of all classes $\\left(\\sum_i \\mathrm{TP}_i + (\\mathrm{FP}_i + \\mathrm{FN}_i)/2\\right)$. In the case of binary classification, this may be expressed as\n\n\\begin{equation}\n\\mathrm{accuracy} = \\frac{\\mathrm{TP} + \\mathrm{TN}}{\\mathrm{TP} + \\mathrm{TN} + \\mathrm{FP} + \\mathrm{FN}}\\text{ .}\n\\label{eq:accuracy_bin}\n\\end{equation}","metadata":{}},{"cell_type":"code","source":"# ---> Compute precisions, recalls, F1-score and accuracy on the whole original training\n# set\n\n# pandas' DataFrame with the performance metrics on the whole original training set\n# [NOTE: the sample_weight argument is such that the confusion matrix is transformed\n# into another using its values as weights. When defining and entry of the matrix, the\n# respective occurrence frequency becomes the sum of the weights of the samples whose\n# labels are the same as the one of the row (the true label) and whose predicted labels\n# are the same as the one of the column.]\ncr_lgb = pd.DataFrame(\n    classification_report(y_true=dt_data[\"train\"][\"y\"],\n                          y_pred=dt_data[\"train\"][\"base\"][\"y_pred_lgb\"],\n                          sample_weight=dt_data[\"train\"][\"base\"][\"sample_weight\"],\n                          output_dict=True)).transpose()\n\n# Making pandas' DataFrame to properly present the accuracy score\n# [NOTE: see https://stackoverflow.com/questions/39662398/scikit-learn-output-metrics-classification-report-into-csv-tab-delimited-format#comment125983216_53780589]\ncr_lgb.loc[\"accuracy\"] = [\"---\", \"---\", cr_lgb.loc[\"accuracy\", \"precision\"], cr_lgb.loc[\"macro avg\", \"support\"]]\n\n# Format \"support\" column so that the values correspond to integers\n# [NOTE: the support value associated with the i-th class corresponds to the number of\n# actual occurrences of that class in the data. This is equivalent to the sum of i-th\n# row's entries of the confusion matrix.]\ncr_lgb[\"support\"] = cr[\"support\"].astype(int) \n\n# Add positive class precision, recall, F1-score and the accuracy to the dataset\n# dictionary\ndt_data[\"train\"][\"LGBM precision\"] = cr_lgb[\"precision\"].loc[\"1\"]\ndt_data[\"train\"][\"LGBM recall\"] = cr_lgb[\"recall\"].loc[\"1\"]\ndt_data[\"train\"][\"LGBM f1-score\"] = cr_lgb[\"f1-score\"].loc[\"1\"]\ndt_data[\"train\"][\"LGBM accuracy\"] = cr_lgb[\"f1-score\"].loc[\"accuracy\"]\n\n# ---> Display performance metric values\nprint()\ndisplay(cr.style.set_caption(\"Classification report\")\\\n        .set_table_styles(\n            [\n                # Column width\n                {\"selector\": \"th.col_heading,td\",\n                 \"props\": [(\"width\", \"100px\")]\n                 },\n                # Caption style\n                {\"selector\": \"caption\",\n                 \"props\": [(\"font-size\", \"16px\"),\n                           (\"font-weight\", \"bold\"),\n                           (\"font-style\", \"italic\")]\n                 }\n            ]))","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.110115Z","iopub.status.idle":"2024-05-05T16:26:01.110449Z","shell.execute_reply.started":"2024-05-05T16:26:01.11027Z","shell.execute_reply":"2024-05-05T16:26:01.110304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the table above, the row `macro avg` has the arithmetic means of the precisions, recalls and $F_1$ scores associated with each class. They are simply the sums of the class-wise entries of the columns above divided by the number of classes. For instance, in the case of the macro-averaged precision, one gets\n\n\\begin{equation}\n\\mathrm{precision}_{\\text{macro-avg}} = \\frac{\\sum_i \\mathrm{precision}_i}{\\sum_i 1}\\text{ .}\n\\label{eq:macro avg}\n\\end{equation}\n\nThe row `weighted avg` has the support-weighted averages. The support associated with the $i$-th class corresponds to the sample weight-scaled number of actual occurrences of such class in the data. This is equivalent to the sum of $i$-th row's entries of a sample weight-scaled confusion matrix. The support-weighted arithmetic means correspond to the class-wise entries of the columns of the classification report multiplied by the respective support value in the same row and divided by the sample weight-scaled number of samples. For instance, the support-weighted average of precision would be given by\n\n\\begin{equation}\n\\mathrm{precision}_{\\text{weighted avg}} = \\frac{\\sum_i \\mathrm{precision}_i\\,\\cdot\\,\\sum_j f_{ij}}{\\sum_i\\sum_j f_{ij}}\\text{ .}\n\\label{eq:weighted avg}\n\\end{equation}\n\nNote that here $f_{ij}$ are sample weight-scaled occurrence frequencies.\n\nIn the case of imbalanced datasets, without balancing sample weights, the numbers of actual occurrences of the minority classes would be much lower than the ones of the majority classes which means that the weights of the minority classes for the support-weighted averaging would be negligible and their scores irrelevant. Therefore, in the case of imbalanced datasets, the macro-averaging might be more useful than the support-weighted averaging.\n\nIn the table above one finds that both average precision and recall are identical, regardless of the averaging method, having a value $\\sim 78 \\%$. Since the $F_1$ score is a harmonic mean of precision and recall, and the averages of these are identical, the average $F_1$ score is also identical to them. Coincidentally, the accuracy is found to be arround this same value. This is not surprising since one may show that accuracy is equivalent to the support-weighted average recall:\n\n\\begin{equation}\n\\mathrm{accuracy} = \\frac{\\sum_i f_{ii}}{\\sum_i\\sum_j f_{ij}} = \\frac{\\sum_i \\mathrm{recall}_i \\sum_jf_{ij}}{\\sum_i\\sum_j f_{ij}} =: \\mathrm{recall}_{\\text{weighted avg}}\\text{ .}\n\\label{eq:acc_recall}\n\\end{equation}\n\nIn a previous experiment it was found that by disregarding sample weights in training and in the evaluation metrics, the recall and $F_1$-score of the positive class achieve very low values ($<1\\%$), revealing that the number of false negatives ($\\mathrm{FN}$) gets to be much higher than the number of true positives ($\\mathrm{TP}$). Such result shows how the minority class (the positive one) is overlooked when sample weights are not considered.","metadata":{}},{"cell_type":"markdown","source":"## The ROC curve, the area under it (AUC) and the Gini coefficient in binary classification\n\nThe implemented Gradient Boosting Decision Tree returns the predicted probability of a given feature $\\{x\\}$ being associated with the label $y=1$ (the case of credit default), that is, $\\hat{P}:=\\hat{P}\\left(y(\\{x\\}) = 1\\right)$. One may then argue which threshold value, $p$, for the predicted probability $\\hat{P}$ should be considered for labelling the feature $\\{x\\}$ with $y=1$. It would be reasonable to consider the threshold to be simply $50\\,\\%$, that is, associating the result $\\hat{P}\\ge 0.5$ with the label $y=1$ label and the result $\\hat{P}< 0.5$ with $y=0$. On the other hand, it would be also wise to choose the threshold value as the one that gives the best result for some performance metric. The receiver operating characteristic (ROC) curve may be useful for identifying this value.\n\n### The ROC curve\n\nThe ROC curve corresponds to the relation between the true positive rate ($\\mathrm{TPR}$) and the false positive rate ($\\mathrm{FPR}$) at different threshold values. The true positive rate is defined as the ratio between true positives ($\\mathrm{TP}$) and actual positives (the sum of true positives ($\\mathrm{TP}$) and false negatives ($\\mathrm{FN}$)):\n\n$$\n\\mathrm{TPR}=\\frac{\\mathrm{TP}}{\\mathrm{TP} + \\mathrm{FN}}\\text{ .}\n$$\n\nAnd the false positive rate ($\\mathrm{FPR}$) is defined as the ratio between false positives ($\\mathrm{FP}$) and actual negatives (the sum of false positives ($\\mathrm{TP}$) and true negatives ($\\mathrm{TN}$)):\n\n$$\n\\mathrm{FPR}=\\frac{\\mathrm{FP}}{\\mathrm{FP} + \\mathrm{TN}}\\text{ .}\n$$\n\nFrom the definitions above, one understands that the true positive and false positive rates are constrained to the range $[0,\\,1]$.\n\nAn ideal model would be one with no false negatives ($\\mathrm{FN}$) neither false positives ($\\mathrm{FP}$), which would mean $\\mathrm{TPR}=1$ and $\\mathrm{FPR}=0$.\n\nThe ROC curve is plotted on a graph of $x$ and $y$-axes associated with $\\mathrm{FPR}$ and $\\mathrm{TPR}$, respectively. The ideal point would then correspond to the upper left corner.\n\nThe first threshold value corresponds to $p=\\infty$, and since by definition, $0\\le\\hat{P}\\le 1$, the condition $\\hat{P}\\ge \\infty$ is impossible. No cases would be positively labelled (with $y=1$), meaning that $\\mathrm{FP}=0$ and $\\mathrm{TP}=0$, and subsequently, $\\mathrm{FPR}=0$ and $\\mathrm{TPR}=0$. The respective ROC point would necessarily be at the origin (lower left corner), that is, $(\\mathrm{TPR},\\,\\mathrm{FPR})_{p=\\infty}=(0,\\,0)$.\n\nThe following threshold value would correspond to the highest obtained predicted probability $\\hat{P}_{\\mathrm{max}}$, meaning that solely the cases with such probability value (which in principle would not be that many) would be positively labelled. \n\nA reasonably accurate model would be such that high predicted probabilities would be solely associated with true positives. Therefore, for $p\\lesssim 1$, no false positives would be expected, meaning $\\mathrm{FP}=0$ and, therefore, $\\mathrm{FPR}=0$. However, such threshold is so high that there would be a lot of false negatives, making the denominator in $\\mathrm{TPR}$ to be quite large and $\\mathrm{TPR}\\gtrsim 0$. One could expect the ROC curve to be close-to-vertical at threshold values $p\\lesssim 1$.\n\nTo obtain the ROC curve one should define the threshold values after the first one as the unique values of the predicted probabilities $\\hat{P}$ in decreasing order. Consider, for instance, that a threshold value between two consecutive unique values was considered, that is, with $p_{i+1}<p<p_{i}$ where $p_{i}$ is the $i$-th greatest unique predicted probability and $p_{i+1}$ is the following smaller one. If there is no unique predicted probability in the range $[p,\\,p_{i}[$, the cases satisfying $\\hat{P}\\ge p$ would be the same as the ones satisfying $\\hat{P}\\ge p_{i}$, and, therefore, there would be no change in the respective predictions, that is, the number of true positives ($\\mathrm{TP}$), true negatives ($\\mathrm{TN}$), false positives ($\\mathrm{FP}$) and false negatives ($\\mathrm{FN}$) would remain the same. And the respective ROC point would coincide with the one for the previous threshold, $p=p_{i}$.\n\nFor the smallest threshold value, $p=\\hat{P}_{\\mathrm{min}}$, the condition $\\hat{P}\\ge p$ always holds, and, therefore, all cases would be positively labelled. Without any negativelly labelled cases, $\\mathrm{FN}=0$ and $\\mathrm{TN}=0$, meaning that $(\\mathrm{TPR},\\,\\mathrm{FPR})_{p=\\hat{P}_{\\mathrm{min}}}=(1,\\,1)$, and the ROC point would be at the upper right corner.\n\n### The ROC curve is non-decreasing\n\nOne may show that the ROC curve cannot decrease, that is, $\\mathrm{TPR}_{p=p_i}\\ge\\mathrm{TPR}_{p=p_{i + 1}}\\,\\forall\\,i$. When going from the threshold value $p_i$ to $p_{i+1}$, cases with predicted probabilities that were lower than the previous threshold, $p_i$, and that are now higher than the current threshold, $p_{i+1}$, were previously negativelly labelled but are now positively labelled. If these cases are truly positive, then, such change would transform false negatives ($\\mathrm{FN}$) into true positives ($\\mathrm{TP}$), decreasing and increasing their numbers, respectively. If the cases are actually truly negative, then, such change would transform true negatives ($\\mathrm{TN}$) into false positives ($\\mathrm{FP}$). Well, the first circumstance would make $\\mathrm{TPR}$ to increase while the second circumstance would not affect $\\mathrm{TPR}$, meaning that $\\mathrm{TPR}$ cannot decrease.\n\n### The ROC curve of a random classifier is the diagonal line $\\mathrm{TPR}=\\mathrm{FPR}$\n\nOne may show that for the case of a random classifier, $\\mathrm{TPR}=\\mathrm{FPR}$. Let $n^+$ the number of actual positives, and $n^-$ the number of actual negatives. Also, let the threshold $p$ be regarded as the probability of the random classifier labelling a case (any, that is, regardless of its feature vector, $\\{\nx\\}$) as positive.\n\nOne would expect that a fraction $p$ of actual positives would be correctly labelled, and a fraction $p$ of actual negatives would be incorrectly labelled as positive. Therefore,\n\n$$\n\\begin{aligned}\n\\mathrm{TP} & = p\\cdot n^+\\text{ ,}\\\\ \\\\\n\\mathrm{FN} & = n^+ - \\mathrm{TP} = (1-p)\\cdot n^+\\text{ ,}\\\\ \\\\\n\\mathrm{FP} & = p\\cdot n^-\\text{ ,}\\\\ \\\\\n\\mathrm{TN} & = n^- - \\mathrm{FP} = (1-p)\\cdot n^-\\text{ .}\n\\end{aligned}\n$$\n\nAnd,\n\n$$\n\\begin{aligned}\n\\mathrm{TPR} & = \\frac{\\mathrm{TP}}{\\mathrm{TP} + \\mathrm{FN}} = \\frac{p\\cdot n^+}{p\\cdot n^+ + (1-p)\\cdot n^+}=p\\text{ ,}\\\\ \\\\\n\\mathrm{FPR} & = \\frac{\\mathrm{FP}}{\\mathrm{FP} + \\mathrm{TN}} = \\frac{p\\cdot n^-}{p\\cdot n^- + (1-p)\\cdot n^-}=p=\\mathrm{TPR}\\text{ .}\\quad\\quad\\quad \\scriptsize{■}\\\\\n\\end{aligned}\n$$\n\n### The ROC curve is invariable to balancing class weights\n\nFor the case of unbalanced datasets, someone agnostic to the structure of ROC curves may argue that it would be more reasonable to compute a \"weighted\" curve in which counts associated with the actual labels are multiplied by their respective balancing weights. Let the balancing weights for the actual negatives and actual positives be $w^-$ and $w^+$, respectively. These weights could be defined as $w^-=1$ and $w^+ = n^-/n^+$. True positives ($\\mathrm{TP}$) and false negatives ($\\mathrm{FN}$) are associated with actual positives while true negatives ($\\mathrm{TN}$) and false positives ($\\mathrm{FP}$) with actual negatives. When \"weighting\" the ROC curve, the terms $\\mathrm{TP}$ and $\\mathrm{FN}$ and the terms $\\mathrm{TN}$ and $\\mathrm{FP}$ involved in the definition of the ROC points $(\\mathrm{FPR},\\,\\mathrm{TPR}$) are multiplied by $w^+$ and $w^-$, respectively. However, one may note that $\\mathrm{FPR}$ and $\\mathrm{TPR}$ are invariable to such transformations, and so are also the ROC points. Indeed, \n\n$$\n\\begin{aligned}\n\\mathrm{TPR} & = \\frac{w^+\\cdot\\mathrm{TP}}{w^+\\cdot\\mathrm{TP} + w^+\\cdot\\mathrm{FN}} = \\frac{\\mathrm{TP}}{\\mathrm{TP} + \\mathrm{FN}}\\text{ ,}\\\\ \\\\\n\\mathrm{FPR} & = \\frac{w^-\\cdot\\mathrm{FP}}{w^-\\cdot\\mathrm{FP} + w^-\\cdot\\mathrm{TN}} =  \\frac{\\mathrm{FP}}{\\mathrm{FP} + \\mathrm{TN}}\\text{ .}\\quad\\quad\\quad \\scriptsize{■}\\\\\n\\end{aligned}\n$$\n\nTherefore, it is irrelevant to consider balancing weights for the ROC curve when the dataset in unbalanced.\n\n### The Area Under the ROC Curve (AUC)\n\nThe Area Under the ROC Curve (AUC) is literally what the term states. Since the ROC curve is constrained to a $1$ by $1$ square, AUC is within the range $[0,\\,1]$. The higher the value of AUC is, the closer the model would be to a perfect classifier. A random classifier would have a $\\mathrm{AUC}$ value of $0.5$. Therefore, values below $0.5$ are not desirable.\n\nThe Gini coefficient is defined as the area under the ROC curve of the model and above the ROC curve of a random classifier, scaled by a factor of $2$:\n\n$$\nG = 2 \\left(\\mathrm{AUC} - 0.5\\right) = 2\\cdot\\mathrm{AUC} - 1\\text{ .}\n$$\n\nSince $\\mathrm{AUC}$ is within $[0,\\,1]$, the Gini coefficient is within $[-1,\\,1]$. A random classifier would have a Gini coefficient of $0$ and a perfect classifier would have a Gini coefficient of $1$. The sign of the Gini coefficient would determine if the model is indeed better (being positive) or not (being negative) than random.","metadata":{}},{"cell_type":"code","source":"# ---> Compute ROC curve, AUC and Gini coefficient for the traning dataset\n\n# Add ROC's FPR and TPR to the dataset dictionary\n(dt_data[\"train\"][\"LGBM FPR\"], dt_data[\"train\"][\"LGBM TPR\"], _) = roc_curve(\n    y_true=dt_data[\"train\"][\"y\"],\n    y_score=dt_data[\"train\"][\"base\"][\"P_pred_lgb\"]\n)\n\n# Add AUC to the dataset dictionary\ndt_data[\"train\"][\"LGBM auc\"] = roc_auc_score(\n    y_true=dt_data[\"train\"][\"y\"],\n    y_score=dt_data[\"train\"][\"base\"][\"P_pred_lgb\"]\n)\n\n# Add Gini coefficient to the dataset dictionary\ndt_data[\"train\"][\"LGBM g\"] = 2 * dt_data[\"train\"][\"LGBM auc\"] - 1","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.111818Z","iopub.status.idle":"2024-05-05T16:26:01.112123Z","shell.execute_reply.started":"2024-05-05T16:26:01.111975Z","shell.execute_reply":"2024-05-05T16:26:01.111987Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Plot ROC curve for the training data\n\n# Initialise figure and axes\nplt.figure(figsize=(4.8, 4.8))\nax = plt.axes()\nplt.title(rf\"ROC curves, AUC and Gini coefficient of LGBM\", pad=20)\n\n# Plot ROC\nax.plot(dt_data[\"train\"][\"LGBM FPR\"],\n        dt_data[\"train\"][\"LGBM TPR\"],\n        linestyle=\"solid\",\n        color=\"black\",\n        alpha=1,\n        linewidth=1.5,\n        label=\"Estimator's ROC\")\n\n# Plot AUC\nax.fill_between(\n    x=dt_data[\"train\"][\"LGBM FPR\"],\n    y1=np.zeros(len(dt_data[\"train\"][\"LGBM FPR\"])),\n    y2=dt_data[\"train\"][\"LGBM TPR\"],\n    edgecolor=\"midnightblue\",\n    facecolor=\"midnightblue\",\n    linewidth=0,\n    # hatch = \"///\",\n    alpha=0.65,\n    label=(r\"Estimator's $\\mathrm{AUC}$ ($\\mathrm{AUC} = \" +\n           rf\"{dt_data['train']['LGBM auc']:.3f}$)\"))\n\n# Plot 1/2 Gini\nax.fill_between(\n    x=dt_data[\"train\"][\"LGBM FPR\"],\n    y1=dt_data[\"train\"][\"LGBM FPR\"],\n    y2=dt_data[\"train\"][\"LGBM TPR\"],\n    edgecolor=\"yellow\",\n    facecolor=\"None\",\n    linewidth=0,\n    hatch = \"xxx\",\n    alpha=1,\n    label=(r\"Estimator's $1/2$ Gini ($G = \" +\n           rf\"{dt_data['train']['LGBM g']:.3f}$)\"))\n\n# Plot random classifier's ROC\nax.plot(\n    [0, 1],\n    [0, 1],\n    linestyle=\"dashed\",\n    color=\"black\",\n    alpha=1,\n    linewidth=1.25,\n    label=\"Random estimator's ROC\")\n\n# Define axes labels                                \nax.set_xlabel(r\"$\\mathrm{FPR}$\", fontdict={\"fontsize\": 10})\nax.set_ylabel(r\"$\\mathrm{TPR}$\", fontdict={\"fontsize\": 10})\n\n# Define axes limits\nax.set_xlim([0, 1])\nax.set_ylim([0, 1])\n\n# Set aspect ratio\n# [NOTE: \"equal\" means that the scales for x and y-axes are equals.]\nax.set_aspect(\"equal\")\n\n# Legend\nax.legend(loc=\"upper left\", fontsize=8)\n\n# Show plot\nplt.show() ","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.113509Z","iopub.status.idle":"2024-05-05T16:26:01.113939Z","shell.execute_reply.started":"2024-05-05T16:26:01.113717Z","shell.execute_reply":"2024-05-05T16:26:01.113735Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## The stability metric\n\nOne may order and group the data points by week number, and then compute a Gini coefficient for each group. Let $G_i$ be the Gini coefficient for the $i$-th group, and $n$ the number of groups.\n\nQualitatively, the stability score considered in this competition is such that it is larger for models that give \n\n* a larger average Gini coefficient,\n    $$\\overline{G}=\\frac{1}{n}\\sum_{i=1}^n G_i\\text{ ;}$$\n\n* whose Gini coefficient points $(i,\\, G_i)$ are closer to its fitting straight line, that is, the root mean square deviation of the linear regression Gini coefficients from the actual ones,\n    $$\\mathrm{RMSD}=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\left(G_{\\mathrm{fit},\\,i} - G_i\\right)^2}\\text{ ,}$$\n    is smaller;\n    \n* and whose slope of such line, $a$, is non-negative or, if negative, it is smaller (in absolute value).\n\nQuantitatively, the stability score is given by\n\n$$ \\mathrm{score} = \\underbrace{\\overline{G}}_{\\text{Base score}} - \\underbrace{0.5\\cdot \\mathrm{RMSD}}_{\\text{Oscillation penalty}}+ \\underbrace{88.0\\cdot \\mathrm{min}(0,\\,a)}_{-\\left(\\text{Negative slope penalty}\\right)}\\text{ .}$$\n\nBasically, it corresponds to the average Gini coefficient subtracted by two penalties: one proportional (by a factor of $0.5$) to the root mean square deviation of the (straight line) approximation Gini coefficients from the actual ones, $\\mathrm{RMSD}$, and another proportional (by a factor of $88.0$) to the negativity of the approximating line slope, $\\mathrm{max}(0,\\,-a) = -\\mathrm{min}(0,\\,a)$.\n\n> **_NOTE:_** The [\"Overview\" page of the competition](https://www.kaggle.com/competitions/home-credit-credit-risk-model-stability/overview) mentions \"standard deviation of the residuals\",\n>\n> $$\\mathrm{s}\\left(\\{\\Delta G\\}\\right)=\\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\left(\\underbrace{\\Delta G_{i}}_{\\text{i-th residual}} - \\underbrace{\\overline{\\Delta G}}_{\\text{Average residual}}\\right)^2} = \\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\left(G_{\\mathrm{fit},\\,i} - G_i - \\overline{\\Delta G}\\right)^2}\\text{ ,}$$\n>\n> instead of the root mean square deviation (RMSD). However, these are actually equivalent since in a linear regression the average residual $\\overline{\\Delta G}$ is null. Indeed, the straight line bias $b$ which minimises the residual sum of squares $\\mathrm{RSS} := \\sum_{i=1}^n\\left(G_{\\mathrm{fit},i} - G_i\\right)^2 = \\sum_{i=1}^n\\left(a\\cdot i + b - G_i\\right)^2$, is defined as\n>\n> $$ b = b^*:\\quad \\frac{\\partial\\,\\mathrm{RSS}}{\\partial\\,b}(b^*) = 0 \\Leftrightarrow 2\\left[a\\cdot\\left(\\sum_{i=1}^n i\\right) + b\\underbrace{\\left(\\sum_{i=1}^n 1\\right)}_{=n} - \\left(\\sum_{i=1}^n G_i\\right)\\right]=0 \\Leftrightarrow$$\n>\n> $$\\Leftrightarrow b = \\frac{1}{n}\\left(\\sum_{i=1}^n G_i\\right) - a\\cdot\\frac{1}{n}\\left(\\sum_{i=1}^n i\\right) \\Leftrightarrow$$\n>\n> $$\\Leftrightarrow b = \\overline{G} - a\\cdot\\overline{i}\\text{ ,}$$\n>\n> which implies that\n>\n> $$\\overline{\\Delta G} = \\overline{\\left(G_{\\mathrm{fit},\\,i} - G_i\\right)} = \\overline{G}_{\\mathrm{fit}} - \\overline{G}= a\\cdot\\overline{i} + b - \\overline{G} = a\\cdot\\overline{i} + \\overline{G} - a\\cdot\\overline{i} - \\overline{G} = 0\\text{ .}$$","metadata":{}},{"cell_type":"code","source":"# ---> Compute stability scores for the training data\n \n# Function for computing stability score's elements\ndef get_stability_score_lgb(\n    # Base pandas dataframe of the dataset dictionary\n    df_base,\n    # Weight for average (in week number) of Gini coefficient\n    w_G_av=1,\n    # Weight for the slope of the Gini coefficient (if negative)\n    w_a=88.0,\n    # Weight for the root mean square deviation of the linear regression Gini\n    # coefficients from the actual ones\n    w_RMSD=-0.5\n):\n    # List of Gini coefficients - one for each week number\n    # [NOTE: the base pandas dataframe is sorted and grouped by WEEK_NUM. The respective\n    # lists of labels (y) and predicted probabilities (P_pred) for each week number are\n    # taken and a respective Gini coefficient is computed.]\n    G_lgb = df_base[[\"WEEK_NUM\", \"y\", \"P_pred_lgb\"]]\\\n        .sort_values(by=\"WEEK_NUM\")\\\n        .groupby(by=\"WEEK_NUM\")[[\"y\", \"P_pred_lgb\"]]\\\n        .apply(lambda x:\n               2 * roc_auc_score(x[\"y\"], x[\"P_pred_lgb\"]) - 1).tolist()\n    \n    # Average (in week number) Gini coefficient\n    G_av_lgb = np.mean(G_lgb)\n\n    # Array of indices for the Gini coefficients\n    i = np.arange(len(G_lgb))\n    \n    # Weight (a) and bias (_) of the linear regression\n    [a, b] = np.polyfit(x=i, y=G_lgb, deg=1)\n    \n    # Array of fit Gini coefficients\n    G_fit = a * i + b\n    \n    # Root mean square deviation of the fit Gini values from the actual ones \n    RMSD_lgb = np.sqrt(np.mean((G_fit - G_lgb)**2))\n\n    # Stability score\n    stability_score_lgb = w_G_av * G_av + w_a * min(0, a) + w_RMSD * RMSD_lgb\n    \n    # Dictionary of stability score elements\n    dt = {\n        \"LGBM_g_week\": G_lgb,\n        \"a\": a,\n        \"b\": b,\n        \"LGBM_RMSD\": RMSD_lgb,\n        \"LGBM_stability_score\": stability_score_lgb\n    }\n    \n    return dt\n\n\n# Add stability score elemements to the dataset dictionary\ndt_data[\"train\"].update(get_stability_score(dt_data[\"train\"][\"LGBM base\"]))","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.117026Z","iopub.status.idle":"2024-05-05T16:26:01.117505Z","shell.execute_reply.started":"2024-05-05T16:26:01.117245Z","shell.execute_reply":"2024-05-05T16:26:01.117263Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Plot stability score's elements for the training data\n\n# Initialise figure and axes\nplt.figure(figsize=(6.4, 4.8))\nax = plt.axes()\nplt.title(rf\"Stability score's elements of LGBM\", pad=20)\n\n# Plot Gini coefficients\nax.plot(range(len(dt_data[\"train\"][\"LGBM_g_week\"])),\n        dt_data[\"train\"][\"LGBM_g_week\"],\n        linestyle=\"None\",\n        marker=\"o\",\n        markeredgecolor=\"black\",\n        markerfacecolor=\"None\",\n        markersize=3,\n        alpha=1,\n        label=\"Points ($G_i$)\")\n\n# Plot linear regression line\ni_lgb = np.array(range(len(dt_data[\"train\"][\"LGBM_g_week\"])))[[0, -1]]\nax.plot(i_lgb,\n        dt_data[\"train\"][\"a\"] * i + dt_data[\"train\"][\"b\"],\n        linestyle=\"solid\",\n        color=\"blue\",\n        alpha=1,\n        linewidth=1,\n        label=(\"Fitting line ($G_{\\mathrm{fit}}(i) = a \\cdot i +b$, with \" +\n               f\"$a=${dt_data['train']['a']:.2e})\"))\n\n# Define axes labels                                \nax.set_xlabel(r\"Week number, $i$\", fontdict={\"fontsize\": 10})\nax.set_ylabel(r\"Gini coefficient, $G$\", fontdict={\"fontsize\": 10})\n\n# Enable axes' minor ticks\nax.minorticks_on()\n\n# Define grid\nax.grid(visible=True, which=\"major\", color=\"lightgray\", linestyle=\"solid\",\n        linewidth=0.5)\nax.grid(visible=True, which=\"minor\", color=\"lightgray\", linestyle=\"dotted\",\n        linewidth=0.5)\n\n# Legend\nax.legend(fontsize=8)\n\n# Show plot\nplt.show() ","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.1186Z","iopub.status.idle":"2024-05-05T16:26:01.118918Z","shell.execute_reply.started":"2024-05-05T16:26:01.118755Z","shell.execute_reply":"2024-05-05T16:26:01.118768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Summary table on performance metrics","metadata":{}},{"cell_type":"code","source":"# Summary dictionary for display\ndt_summary_LGBM = {\n    \"dataset\": \"train\",\n    \"precision\": dt_data[\"train\"][\"LGBM precision\"],\n    \"recall\": dt_data[\"train\"][\"LGBM recall\"],\n    \"f1-score\": dt_data[\"train\"][\"LGBM f1-score\"],\n    \"accuracy\": dt_data[\"train\"][\"LGBM accuracy\"],\n    \"auc\": dt_data[\"train\"][\"LGBM auc\"],\n    \"g\": dt_data[\"train\"][\"LGBM g\"],\n    \"stability_score\": dt_data[\"train\"][\"LGBM_stability_score\"]\n}\n\n# Display summary\nprint()\ndisplay(pd.DataFrame(data=dt_summary_LGBM, index=[0])\\\n        .style.set_caption(\"LGBM Performance metrics for the training data\")\\\n        .set_table_styles(\n            [\n                # Column width\n                {\"selector\": \"th.col_heading,td\",\n                 \"props\": [(\"width\", \"100px\")]\n                 },\n                # Caption style\n                {\"selector\": \"caption\",\n                 \"props\": [(\"font-size\", \"16px\"),\n                           (\"font-weight\", \"bold\"),\n                           (\"font-style\", \"italic\")]\n                 }\n            ]))","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.120384Z","iopub.status.idle":"2024-05-05T16:26:01.120735Z","shell.execute_reply.started":"2024-05-05T16:26:01.120574Z","shell.execute_reply":"2024-05-05T16:26:01.120588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Submit results","metadata":{}},{"cell_type":"markdown","source":"As mentioned in the [\"Overview\" page of the competition](https://www.kaggle.com/competitions/home-credit-credit-risk-model-stability/overview), one needs to compute probabilies of credit default for the issued test dataset and submit a csv (Comma-Separated Values) file with the first column having the credit case and the second column having the probability values. Note that the respective column headers (\"case_id\" and \"score\") need to be defined. Also, note that the actual labels of the test dataset are not issued and, therefore, there is no way to know the performance of the model before submitting the file. \n\nThe file is submitted by going to the \"Submit to competition\" section on the right-hand side pane and clicking on the \"Submit\" button. The stability score should be shown on that section after submitting.","metadata":{}},{"cell_type":"code","source":"# ---> Add predicted probabilities for label corresponding to 1 (credit default case)\n# to base pandas dataframe of the test dataset\n# [NOTE: the method predict_proba returns a numpy.ndarray of how many columns there are\n# classes. The columns are ordered according to the attribute classes_ of the\n# classifier.]\ndt_data[\"test\"][\"base\"][\"P_pred\"] = best_estimator.predict_proba(\n    X=dt_data[\"test\"][\"x\"]\n)[:, 1]","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.122055Z","iopub.status.idle":"2024-05-05T16:26:01.122357Z","shell.execute_reply.started":"2024-05-05T16:26:01.122208Z","shell.execute_reply":"2024-05-05T16:26:01.12222Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Create submission pandas dataframe\n# Submission dataframe\nsubmission = pd.DataFrame({\n    \"case_id\": dt_data[\"test\"][\"base\"][\"case_id\"].to_numpy(),\n    \"score\": dt_data[\"test\"][\"base\"][\"P_pred\"]\n})\n\n# Display dataframe\nprint()\ndisplay(submission.style.set_caption(\"Submission dataframe\")\\\n        .set_table_styles(\n            [\n                # Column width\n                {\"selector\": \"th.col_heading,td\",\n                 \"props\": [(\"width\", \"100px\")]\n                 },\n                # Caption style\n                {\"selector\": \"caption\",\n                 \"props\": [(\"font-size\", \"16px\"),\n                           (\"font-weight\", \"bold\"),\n                           (\"font-style\", \"italic\")]\n                 }\n            ]))","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.12371Z","iopub.status.idle":"2024-05-05T16:26:01.124045Z","shell.execute_reply.started":"2024-05-05T16:26:01.123886Z","shell.execute_reply":"2024-05-05T16:26:01.123899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Create the csv submission file\n# [NOTE: after running the command below, the file \"submission.csv\" will appear on the\n# right-hand side pane, at the \"Output\" section.]\nsubmission.to_csv(\"submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.125066Z","iopub.status.idle":"2024-05-05T16:26:01.125386Z","shell.execute_reply.started":"2024-05-05T16:26:01.125228Z","shell.execute_reply":"2024-05-05T16:26:01.125241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 决策树\n\n# 假设 X_train, y_train 是你的训练数据和标签\n# 创建决策树分类器的实例\ndt_classifier = DecisionTreeClassifier(random_state=42)\n\n# 定义要尝试的参数的字典\nparam_grid = {\n    \"criterion\": [\"gini\", \"entropy\"],  # 决策树的划分标准\n    \"max_depth\": [None, 10, 20, 30],    # 树的最大深度\n    \"min_samples_split\": [2, 5, 10],    # 分裂内部节点所需的最小样本数\n    \"min_samples_leaf\": [1, 2, 4]       # 叶节点所需的最小样本数\n}\n\n# 创建 GridSearchCV 对象\ngrid_search = GridSearchCV(\n    estimator=dt_classifier,\n    param_grid=param_grid,\n    cv=5,  # 交叉验证的折数\n    scoring=\"accuracy\",  # 评分标准\n    n_jobs=1,  # 并行工作的进程数\n    refit=True,  # 使用整个训练集对找到的最佳模型进行拟合\n    verbose=0  # 输出的详细程度\n)\n\n# 记录初始时间\nt_i = time.perf_counter()\n\n# 在训练数据上执行 GridSearchCV\ngrid_search.fit(X=dt_data[\"train\"][\"x\"], y=dt_data[\"train\"][\"y\"])\n\n# 记录结束时间\nt_f = time.perf_counter()\n\n# 输出运行时间\nprint()\nprint(f\"Running time: {(t_f - t_i)/3600:.2f} h\")\n\n# 输出最佳参数和最佳模型的得分\nprint(\"Best parameters found: \", grid_search.best_params_)\nprint(\"Best cross-validated score: \", grid_search.best_score_)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.127073Z","iopub.status.idle":"2024-05-05T16:26:01.1274Z","shell.execute_reply.started":"2024-05-05T16:26:01.127241Z","shell.execute_reply":"2024-05-05T16:26:01.127255Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 假设 X_train, y_train 是你的训练数据和标签\n\n# 创建随机森林分类器的实例\nrf_classifier = RandomForestClassifier(random_state=42)\n\n# 定义要尝试的参数的字典\nparam_grid = {\n    \"n_estimators\": [100, 200, 300],  # 决策树的数量\n    \"max_depth\": [None, 10, 20, 30],   # 树的最大深度\n    \"min_samples_split\": [2, 5, 10],   # 分裂内部节点所需的最小样本数\n    \"min_samples_leaf\": [1, 2, 4],     # 叶节点所需的最小样本数\n    \"bootstrap\": [True, False]          # 是否使用bootstrap样本\n}\n\n# 创建 GridSearchCV 对象\ngrid_search = GridSearchCV(\n    estimator=rf_classifier,\n    param_grid=param_grid,\n    cv=5,  # 交叉验证的折数\n    scoring=\"accuracy\",  # 评分标准\n    n_jobs=-1,  # 使用所有可用的CPU核心进行并行处理\n    refit=True,  # 使用整个训练集对找到的最佳模型进行拟合\n    verbose=0  # 输出的详细程度\n)\n\n# 记录初始时间\nt_i = time.perf_counter()\n\n# 在训练数据上执行 GridSearchCV\ngrid_search.fit(X=dt_data[\"train\"][\"x\"], y=dt_data[\"train\"][\"y\"])\n\n# 记录结束时间\nt_f = time.perf_counter()\n\n# 输出运行时间\nprint()\nprint(f\"Running time: {(t_f - t_i)/3600:.2f} h\")\n\n# 输出最佳参数和最佳模型的得分\nprint(\"Best parameters found: \", grid_search.best_params_)\nprint(\"Best cross-validated score: \", grid_search.best_score_)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.129077Z","iopub.status.idle":"2024-05-05T16:26:01.129386Z","shell.execute_reply.started":"2024-05-05T16:26:01.129231Z","shell.execute_reply":"2024-05-05T16:26:01.129243Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> XGBoost的梯度提升决策树（估计器模板）\n​\n# 用于训练XGBoost估计器的参数字典\nparams = {\n    \"boosting_type\": \"gbdt\",\n    \"objective\": \"binary\",\n    \"metric\": \"auc\",\n    \"num_leaves\": 31,\n    \"learning_rate\": 0.05,\n    \"feature_fraction\": 1,\n    \"bagging_fraction\": 1,\n    \"bagging_freq\": 0,\n    \"n_estimators\": 1200,\n    \"max_depth\": 3,\n    \"n_jobs\": -1,\n    \"random_state\": 42,\n    \"device_type\": \"gpu\",\n    \"gpu_use_dp\": True,\n    \"verbose\": -1,\n}\n​\n# XGBoost估计器\nxgb_estimator = xgb.XGBClassifier(**params)\nadd Codeadd Markdown\n# ---> Perform GridSearchCV on LightGBM's estimator\n​\n# Dictionary of parameters' lists of values to be tried in GridSearchCV\n# [NOTE: the keys of this dictionary should coincide with the names of the estimator\n# arguments that represent such hyperparamenters.]\nparam_grid = {\n    \"n_estimators\": [600, 800, 1000],\n    \"max_depth\": [8, 9, 10]\n}\n​\n# Create a GridSearchCV object\ngs = GridSearchCV(\n    estimator=xgb_classifier,\n    param_grid=param_grid,\n    cv=5,  # 交叉验证的折数\n    scoring=\"roc_auc\",  # 评分标准\n    n_jobs=1,  # 并行工作的进程数\n    refit=True,  # 使用整个训练集对找到的最佳模型进行拟合\n    verbose=0\n)\n​\n# Initial instant of time [s]\nt_i = time.perf_counter()\n​\n# Perform GridSearchCV on the training data\n# [NOTE: to perform a balanced training a sample weight array is issued. Note that\n# GridSearchCV uses it only in the cost function and not in the scorer. Anyway, the ROC\n# curve is actually invariable to class weights. Sample weighting could then in this\n# case be disregarded in scoring.]\ngs.fit(\n    X=dt_data[\"train\"][\"x\"],\n    y=dt_data[\"train\"][\"y\"],\n    sample_weight=dt_data[\"train\"][\"base\"][\"sample_weight\"]\n)\n​\n# Final instant of time [s]\nt_f = time.perf_counter()\n​\n# Print elapsed time\nprint()\nprint(f\"Running time: {(t_f - t_i)/3600:.2f} h\")","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.131129Z","iopub.status.idle":"2024-05-05T16:26:01.131588Z","shell.execute_reply.started":"2024-05-05T16:26:01.131343Z","shell.execute_reply":"2024-05-05T16:26:01.131361Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Appendix: Shapley values","metadata":{}},{"cell_type":"markdown","source":"The Shapley value of the feature component $x^{(j)}$ with respect to the $i$-th data point, herein denoded as $\\phi_{i}^{j}$, corresponds to a measure of the contribution of such component to a shift on the predicted probability at the $i$-th point when the component is included in the model.\n\nLet $\\mathcal{N}={1,\\,\\dots,\\,d}$ be the whole set of indices of the feature components, where $d$ is their number. And let $\\mathcal{S}\\subseteq\\mathcal{N}$ be a subset of that set.\n\nFurthermore, let $p(\\{x\\})$ be the density probability function of the whole feature vector (such function may be inferred from the distribution of feature values in the dataset). And let $P(\\{x\\})$ be the predicted probability of the feature vector $\\{x\\}$ being associated with the positive label $y=1$. The mean predicted probability is defined as\n\n$$\n\\bar{P}=\\mathbb{E}_{\\{x\\}}\\left[P\\right]=\\int P\\left(\\{x\\}\\right)\\cdot p\\left(\\{x\\}\\right)\\,d\\{x^{}\\}\\text{ .}\n$$\n\nThe partial dependence of the predicted probability on the subset of component features $\\{x^{(\\mathcal{S})}\\}\\subseteq \\{x\\}$, or in other words, the marginal predicted probability of variable $\\{x^{(\\mathcal{S})}\\}$ corresponds to a mean predicted probability in which the mean is done solely on the complementary subset, $\\{x^{(\\mathcal{C})}\\}=\\{x\\}\\setminus \\{x^{(\\mathcal{S})}\\}$:\n\n$$\nP\\left(\\{x^{(\\mathcal{S})}\\}\\right)=\\mathbb{E}_{\\{x^{(\\mathcal{C})}\\}}\\left[P\\right]=\\int P\\left(\\{x\\}\\right)\\cdot p\\left(\\{x\\}\\right)\\,d\\{x^{(\\mathcal{S})}\\}\\text{ .}\n$$\n\nThe sensible marginal probability of variable $\\{x^{(\\mathcal{S})}\\}$ would be the difference between the marginal probability and the mean probability:\n\n$$\nv\\left(\\{x^{(\\mathcal{S})}\\}\\right) = P\\left(\\{x^{(\\mathcal{S})}\\}\\right) - \\bar{P}\\text{ .}\n$$\n\nThe Shapley value associated with the component $x^{(j)}$ for the $i$-th data point is then defined as\n\n$$\n\\begin{align}\n\\phi_{i}^{(j)} & = \\overbrace{\\frac{1}{d}\\underbrace{\\sum_{m=0}^{d-1}}_{\\text{$n$ terms: the sizes} \\\\ \\text{of the subsets $\\mathcal{S}$}}}^{\\text{Averaging in all subset sizes}}\\overbrace{\\frac{1}{\\binom{d-1}{n-m-1}}\\underbrace{\\sum_{\\mathcal{S}\\subseteq \\mathcal{N}\\setminus\\{j\\}\\\\|\\mathcal{S}|=m}}_{\\text{$\\binom{d-1}{d-m-1}$ terms: all}\\\\\\text{subsets $\\mathcal{S}$ of size $m$}}}^{\\text{Averaging in all $m$-size subsets $\\mathcal{S}$}} \\underbrace{v\\left(\\{x^{(\\mathcal{S}\\,\\cup\\,\\{j\\})}\\}\\right) - v\\left(\\{x^{(\\mathcal{S})}\\}\\right)}_{\\text{Shift in sensible marginal probability} \\\\ \\text{when going from the index subset $\\mathcal{S}$} \\\\ \\text{of size $m$ to $\\mathcal{S}$ plus the $j$-th index}}\\\\\\\\\n& = \\frac{1}{d} \\sum_{\\mathcal{S}\\subseteq \\mathcal{N}\\setminus\\{j\\}} \\binom{d-1}{d-|\\mathcal{S}|-1}^{-1} v\\left(\\{x^{(\\mathcal{S}\\,\\cup\\,\\{j\\})}\\}\\right) - v\\left(\\{x^{(\\mathcal{S})}\\}\\right)\n\\text{ .}\n\\end{align}\n$$\n\nThe term $\\binom{d-1}{d-|\\mathcal{S}|-1} = \\frac{\\left(d-1\\right)!}{\\left(d-|\\mathcal{S}|-1\\right)!|\\mathcal{S}|!}$ is a binomial coefficient which gives the number of $|\\mathcal{S}|$-size subsets $\\mathcal{S}$ that may be created from the $d$ elements except $j$ of the set $\\mathcal{N}$. Such number is equivalent to distributing $|\\mathcal{S}|$ undistinguishable positions (the ones in the subset $\\mathcal{S}$, undisguishable because the order is irrelevant) to $d-1$ distinguishable indices (the ones except $j$ in the set $\\mathcal{N}$) without repetition. Indeed, $\\left(d-1\\right)\\cdot\\,\\dots\\,\\cdot\\left(d-|\\mathcal{S}|\\right)=\\frac{\\left(n-1\\right)!}{\\left(n-|\\mathcal{S}|-1\\right)!}$ is the number of ways of placing $|\\mathcal{S}|$  indices from the $d-1$ distinguishable ones into $|\\mathcal{S}|$ hypothetically distinguishable positions. And $|\\mathcal{S}|\\cdot\\dots\\cdot1=|\\mathcal{S}|!$ is the number of ways of ordering the $|\\mathcal{S}|$ distinguishable distinguishable positions. Since the positions are actually undistinguishable, the so-wanted actual number is the former divided by the latter, $\\frac{\\left(d-1\\right)!}{\\left(d-|\\mathcal{S}|-1\\right)!|\\mathcal{S}|!} = \\binom{d-1}{d-|\\mathcal{S}|-1}$.\n\nA positive Shapley value would mean that the feature component has a positive contribution for the predicted probability of the respective point. In the current problem that would imply that the existence of such component and its respective value increases the probability of credit default. A negative Shapley value would mean the converse.\n\nShapley values may be computed in this notebook by using the methods of [SHAP](https://shap.readthedocs.io/en/latest/#)'s [Explainer](https://shap.readthedocs.io/en/latest/generated/shap.Explainer.html) class.\n\nThe most important feature components would have the highest absolute Shapley values. One could define the [importance of the $j$-th feature component](https://christophm.github.io/interpretable-ml-book/shap.html#shap-feature-importance) as the arithmetic mean of the respective absolute Shapley values:\n\n$$\nI^{(j)} = \\frac{1}{n}\\sum_{i=1}^{n}|\\phi_i^{(j)}|\\text{ ,}\n$$\n\nwhere $n$ is the number of data points.","metadata":{}},{"cell_type":"code","source":"# ---> Define SHAP's Explainer object\n\n# SHAP's masker - an object used by SHAP's Explainer for constructing the probability\n# density function of the data and the output function of the the model, based on the\n# given data\nmasker = shap.maskers.Independent(\n    # The base data\n    data=dt_data[\"train\"][\"x\"],\n    # The maximum number of samples to taken from the base data\n    max_samples=100)\n\n# SHAP's Explainer object for the probability prediction method of the built model\n# [NOTE: predict_proba returns both probabilities of deafault and non-default. That\n# would make SHAP to also compute Shapley values to both cases. However, solely the ones\n# for the latter are to be considered.]\nexplainer = shap.Explainer(model=best_estimator.predict_proba, masker=masker)","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.133047Z","iopub.status.idle":"2024-05-05T16:26:01.133508Z","shell.execute_reply.started":"2024-05-05T16:26:01.133254Z","shell.execute_reply":"2024-05-05T16:26:01.133272Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Compute Shapley values\n\n# Initial instant of time [s]\nt_i = time.perf_counter()\n\n# Shapley values for the training data\n# [NOTE: not all data was taken as it would require too much computational time.]\nshapley = explainer(dt_data[\"train\"][\"x\"].sample(frac=0.005, random_state=42))\n\n# Final instant of time [s]\nt_f = time.perf_counter()\n\n# Print elapsed time\nprint()\nprint(f\"Running time: {(t_f - t_i)/3600:.2f} h\")","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.134615Z","iopub.status.idle":"2024-05-05T16:26:01.134928Z","shell.execute_reply.started":"2024-05-05T16:26:01.134772Z","shell.execute_reply":"2024-05-05T16:26:01.134784Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Plot features' importance\n# [NOTE: in shap.plots.bar the features are by default sorted in decreasing value of the\n# absolute value of the given shap_values (in this case, this corresponds to the feature\n# importance defined above).]\nshap.plots.bar(\n    shap_values = shapley[:,:,1].abs.mean(0),\n    # How many top features to include in the plot \n    max_display=11,\n    # False, to not immediately plot, and later use matplotlib to configure the plot\n    show=False\n)\n\n# Set plot title\nplt.title(\n    label=\"Feature importance, $I^{(j)}$\\n (the mean absolute Shapley value, $|\\phi^{(j)}|_{\\mathrm{av}}$)\",\n    loc=\"center\",\n    pad=20,\n    fontdict={\n        \"fontsize\": 15,\n        \"verticalalignment\": \"baseline\",\n        \"horizontalalignment\": \"center\"})\n\n# Show plot\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.136272Z","iopub.status.idle":"2024-05-05T16:26:01.136606Z","shell.execute_reply.started":"2024-05-05T16:26:01.136423Z","shell.execute_reply":"2024-05-05T16:26:01.136436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Plot Shapley values\n# [NOTE: in shap.plots.beeswarm, the features are by default sorted in decreasing value\n# of the mean of the given shap_values (in this case, this corresponds to the feature\n# importance defined above).]\nshap.plots.beeswarm(\n    shap_values = shapley[:,:,1],\n    # How many top features to include in the plot \n    max_display=11,\n    # Colour map\n    color=shap.plots.colors.red_blue,\n    # False, to not immediately plot, and later use matplotlib to configure the plot\n    show=False\n)\n\n# Set plot title\nplt.title(\n    label=\"Features and respective Shapley values,\\n sorted by their importance ($I^{(j)}$)\",\n    loc=\"center\",\n    pad=20,\n    fontdict={\n        \"fontsize\": 15,\n        \"verticalalignment\": \"baseline\",\n        \"horizontalalignment\": \"center\"})\n\n# Show plot\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-05-05T16:26:01.137796Z","iopub.status.idle":"2024-05-05T16:26:01.138129Z","shell.execute_reply.started":"2024-05-05T16:26:01.13797Z","shell.execute_reply":"2024-05-05T16:26:01.137984Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"One may then infer that according to the Shapley values, the $10$ most important features are\n\n* `sex_applicant_738L`, (studied before) applicant's gender - males associated with higher probabilities of defaulting;\n\n* `n_persons_same_phone_A`, (studied before) number of persons of the same contract using the same mobile phone number - larger values associated with higher probabilities of defaulting;\n\n* `age_applicant_A`, (studied before) applicant's age - smaller values associated with with higher probabilities of defaulting;\n\n* `n_payments_A`, (studied before) total number of payments made by the applicant - larger values associated with higher probabilities of defaulting;\n\n* `monthly_payment_A`, (studied before) contract's monthly payment value - larger values associated with higher probabilities of defaulting;\n\n* `requesttype_4525192L`, tax authority request type;\n\n* `avgdpdtolclosure24_3658938P`, average number of days past due with tolerance withing the last $24$ months from the maximum closure date, assuming that the contract is finished (if the contract is ongoing, the calculation is based on the current date) - larger values associated with higher probabilities of defaulting;\n\n* `eir_270L`, interest rate - larger values associated with higher probabilities of defaulting;\n\n* `pmtssum_45A`, sum of tax deductions for the client - smaller values associated with higher probabilities of defaulting;\n\n* `totaldebt_9A`, total amount of debt - larger values associated with higher probabilities of defaulting.\n\n<!-- * `employment_time_applicant_A`, (studied before) applicant's employment time at decision date - smaller values associated with higher probabilities of defaulting; -->\n\n<!-- * `homephncnt_628L`, number of distinct home phones on client's application - smaller values associated with higher probabilities of defaulting. -->\n","metadata":{}}]}