{"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":"none","dataSources":[{"sourceId":50160,"databundleVersionId":7921029,"sourceType":"competition"}],"dockerImageVersionId":30646,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"<h1 style=\"font-size:26pt;\"> \n    <center>\n        The Starter Notebook - thoroughly explained\n    </center>\n</h1>\n\n# Introduction\n\nThis work thoroughly explains [@danielherman](https://www.kaggle.com/jetakow)'s amazing [Home Credit 2024 Starter Notebook](https://www.kaggle.com/code/jetakow/home-credit-2024-starter-notebook). Even though such original notebook is a good starter, I think that it may be hard for the less experienced users to grasp it, as the code is not that much commented, there is no introduction to the theoretical background, neither some useful visualisations. I intended in this new version to add to [@danielherman](https://www.kaggle.com/jetakow)'s awesome work some of these lacking \"features\".","metadata":{}},{"cell_type":"markdown","source":"***","metadata":{}},{"cell_type":"markdown","source":"# Import Python modules","metadata":{}},{"cell_type":"code","source":"# ---> General\n# warnings - to manage irrelevant warnings when running 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\n# ---> Manage data\n# NumPy - for basic mathematical operations\nimport numpy as np\n# polars - to efficiently manage large quantities of data and to deal with tables\nimport polars as pl\n# Pandas - to deal with tables \nimport pandas as pd\n\n# ---> Modelling\n# LightGBM - for applying gradient-boosting algorithms\nimport lightgbm as lgb\n# Import some scikit-learn's metrics \nfrom sklearn.metrics import (\n    accuracy_score,\n    roc_curve,\n    roc_auc_score\n)","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:05.538115Z","iopub.execute_input":"2024-03-14T11:05:05.539267Z","iopub.status.idle":"2024-03-14T11:05:11.573948Z","shell.execute_reply.started":"2024-03-14T11:05:05.539174Z","shell.execute_reply":"2024-03-14T11:05:11.572695Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"***","metadata":{}},{"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 added to 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 = \"/kaggle/input/home-credit-credit-risk-model-stability/\"","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:11.576046Z","iopub.execute_input":"2024-03-14T11:05:11.576526Z","iopub.status.idle":"2024-03-14T11:05:11.582520Z","shell.execute_reply.started":"2024-03-14T11:05:11.576486Z","shell.execute_reply":"2024-03-14T11:05:11.581082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"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 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 data in this very kaggle kernel.\n\nFor the sake of simplicity not all issued data files were used in this notebook, neither all of their fields (columns). \n\n## The data files\n\nThe CSV 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\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| `data_decision` | str  | Date of decision for the approval of the credit (format \"YYYY-MM-DD\"). | No | --- |\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 `<string>_<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). But solely \"M\" and \"A\" types were used in this notebook;\n\n| Column          | Type  | Description | Used? | As? |\n| --------------- | ----- | ----------- | ----- | --- |\n| `case_id`       | int   | Same as base files' `case_id`. | Yes | Case identifier |\n| `<string>M`     | str   | Feature component of M-type (\"masking categories\"). | Yes | Feature component, $x^{(j)}$ |\n| `<string>A`     | float | Feature component of A-type (\"transform amount\"). | Yes | Feature component, $x^{(j)}$ |\n| `<other>`       |       | Feature components of other types. | No | --- |\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\". But as in the case of the static internal files, solely \"M\" and \"A\" types were used in this notebook;\n        \n| Column          | Type  | Description | Used? | As? |\n| --------------- | ----- | ----------- | ----- | --- |\n| `case_id`       | int   | Same as base files' `case_id`. | Yes | Case identifier |\n| `<string>M`     | str   | Feature component of M-type (\"masking categories\"). | Yes | Feature component $x^{(j)}$ |\n| `<string>A`     | float | Feature component of A-type (\"transform amount\"). | Yes | Feature component $x^{(j)}$ |\n| `<other>`       |       | Feature components of other types. | No | --- |\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\". But solely three feature components were used. Furthermore, for compactness reasons, 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_L` <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| `<other>`               |       | Other feature components. | No | --- |\n    \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. 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        if col[-1] in (\"M\"):\n            df = df.with_columns(pl.col(col).cast(pl.Categorical))\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    return df","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:11.584728Z","iopub.execute_input":"2024-03-14T11:05:11.585145Z","iopub.status.idle":"2024-03-14T11:05:11.599463Z","shell.execute_reply.started":"2024-03-14T11:05:11.585106Z","shell.execute_reply":"2024-03-14T11:05:11.597993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Create polars dataframes for training data\n\n# Define polars dataframe containing data from base training CSV file\ndf_train_base = pl.read_csv(path_data + \"csv_files/train/train_base.csv\")\n\n# Define polars dataframe containing data from static internal training CSV files 0 and\n# 1\ndf_train_static = pl.concat([pl.read_csv(path_data + \"csv_files/train/train_static_0_0.csv\").pipe(set_pl_dtypes),\n                             pl.read_csv(path_data + \"csv_files/train/train_static_0_1.csv\").pipe(set_pl_dtypes)],\n                            # vertically concatenate, while conveniently redefining\n                            # column's data types to support the data from common\n                            # columns\n                            how=\"vertical_relaxed\")\n\n# Define polars dataframe containing data from static external (from a Credit Bureau)\n# training CSV files\ndf_train_static_cb = pl.read_csv(path_data + \"csv_files/train/train_static_cb_0.csv\").pipe(set_pl_dtypes)\n\n# Define polars dataframe containing data from the persons' training CSV file of depth 1\ndf_train_person_1 = pl.read_csv(path_data + \"csv_files/train/train_person_1.csv\").pipe(set_pl_dtypes)\n\n# Define polars dataframe containing data from Credit Bureau B's training CSV file of\n# depth 2\ndf_train_credit_bureau_b_2 = pl.read_csv(path_data + \"csv_files/train/train_credit_bureau_b_2.csv\").pipe(set_pl_dtypes)","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:11.603290Z","iopub.execute_input":"2024-03-14T11:05:11.603809Z","iopub.status.idle":"2024-03-14T11:05:38.874497Z","shell.execute_reply.started":"2024-03-14T11:05:11.603763Z","shell.execute_reply":"2024-03-14T11:05:38.873087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Create polars dataframes for test data\n\n# Define polars dataframe containing data from base test CSV file\ndf_test_base = pl.read_csv(path_data + \"csv_files/test/test_base.csv\")\n\n# Define polars dataframe containing data from static internal test CSV files 0 and\n# 1\ndf_test_static = pl.concat(\n    [pl.read_csv(path_data + \"csv_files/test/test_static_0_0.csv\").pipe(set_pl_dtypes),\n     pl.read_csv(path_data + \"csv_files/test/test_static_0_1.csv\").pipe(set_pl_dtypes),\n     pl.read_csv(path_data + \"csv_files/test/test_static_0_2.csv\").pipe(set_pl_dtypes)],\n    # vertically concatenate, while conveniently redefining column's data types to\n    # support the data from common columns\n    how=\"vertical_relaxed\")\n# Define polars dataframe containing data from static external (from a Credit Bureau)\n# test CSV files\ndf_test_static_cb = pl.read_csv(path_data + \"csv_files/test/test_static_cb_0.csv\").pipe(set_pl_dtypes)\n\n# Define polars dataframe containing data from the persons' test CSV file of depth 1\ndf_test_person_1 = pl.read_csv(path_data + \"csv_files/test/test_person_1.csv\").pipe(set_pl_dtypes)\n\n# Define polars dataframe containing data from Credit Bureau B's test CSV file of depth\n# 2\ndf_test_credit_bureau_b_2 = pl.read_csv(path_data + \"csv_files/test/test_credit_bureau_b_2.csv\").pipe(set_pl_dtypes)","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:38.876147Z","iopub.execute_input":"2024-03-14T11:05:38.876603Z","iopub.status.idle":"2024-03-14T11:05:38.973626Z","shell.execute_reply.started":"2024-03-14T11:05:38.876566Z","shell.execute_reply":"2024-03-14T11:05:38.972254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Define selecting columns for the static dataframes\n\n# List of labels of the selected columns of training and test static internal dataframes\n# [NOTE: for the sake of simplicity, solely A and M-types are herein selected. These\n# stand for \"transform amount\" and \"masking categories\" of the groups of transforms,\n# with capital letters \"A\" and \"M\" as the last characters of the labels, respectively.]\ncols_static_selec = []\nfor col in df_train_static.columns:\n    if col[-1] in (\"A\", \"M\"):\n        cols_static_selec.append(col)\n\n# List of labels of the selected columns of training and test static Credit Bureau\n# dataframes\n# [NOTE: for the sake of simplicity, solely A and M-types are herein selected.]\ncols_static_cb_selec = []\nfor col in df_train_static_cb.columns:\n    if col[-1] in (\"A\", \"M\"):\n        cols_static_cb_selec.append(col)","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:38.975563Z","iopub.execute_input":"2024-03-14T11:05:38.975997Z","iopub.status.idle":"2024-03-14T11:05:38.984455Z","shell.execute_reply.started":"2024-03-14T11:05:38.975961Z","shell.execute_reply":"2024-03-14T11:05:38.982966Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Filter training dataframes\n\n# Group df_train_person_1 by case_id (that is, collapse rows with same case_id value,\n# which with an additional \"aggregation\" (to be performed later) makes the entries of\n# the chosen columns to accomodate lists with the respective collapsed values (or\n# operations on them, as wished))\ndf_train_person_1_group = df_train_person_1.group_by(\"case_id\")\n\n\n# Aggregate to the group, the maximum income (from main occupation) of the set of people\n# associated with each group, as well as a boolean that holds true if any of the set of\n# people of each group is self-employed\ndf_train_person_1_1 = df_train_person_1_group.agg(\n    pl.col(\"mainoccupationinc_384A\").max().alias(\"mainoccupationinc_max_A\"),\n    (pl.col(\"incometype_1044T\") == \"SELFEMPLOYED\").any().alias(\"anyselfemployed_T\")\n)\n\n# Select from client df_train_person_1 columns \"case_id\", \"num_group1\" and\n# \"housetype_905L\". Keep solely rows associated with applicants (num_group1 = 0), drop\n# the column \"num_group1\" and rename \"housetype_905L\" to \"housetype_applicant_L\", since\n# it refers now to the house type (e.g. \"owned\") of the applicants.\ndf_train_person_1_2 = df_train_person_1.select(\n    [\"case_id\", \"num_group1\", \"housetype_905L\"]\n).filter(pl.col(\"num_group1\") == 0).drop(\"num_group1\").rename(\n    {\"housetype_905L\": \"housetype_applicant_L\"}\n)\n\n# Group df_train_credit_bureau_b_2 by case_id and aggregate maximum value of the number\n# of overdue payments (pmts_pmtsoverdue_635A) and a boolean that holds true if there is\n# any number of days of overdue payment greater than 31.\ndf_train_credit_bureau_b_2 = df_train_credit_bureau_b_2.group_by(\"case_id\").agg(\n    pl.col(\"pmts_pmtsoverdue_635A\").max().alias(\"pmts_pmtsoverdue_max_A\"),\n    (pl.col(\"pmts_dpdvalue_108P\") > 31).any().alias(\"pmts_dpdvalue_anyover31_P\")\n)\n\n# Training dataframe defined by a left join of all training dataframes.\ndf_train = df_train_base\\\n.join(df_train_static.select([\"case_id\"] + cols_static_selec),\n      how=\"left\",\n      on=\"case_id\")\\\n.join(df_train_static_cb.select([\"case_id\"] + cols_static_cb_selec),\n      how=\"left\",\n      on=\"case_id\")\\\n.join(df_train_person_1_1,\n      how=\"left\",\n      on=\"case_id\")\\\n.join(df_train_person_1_2,\n      how=\"left\",\n      on=\"case_id\")\\\n.join(df_train_credit_bureau_b_2,\n      how=\"left\",\n      on=\"case_id\")","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:38.986653Z","iopub.execute_input":"2024-03-14T11:05:38.987885Z","iopub.status.idle":"2024-03-14T11:05:43.898489Z","shell.execute_reply.started":"2024-03-14T11:05:38.987832Z","shell.execute_reply":"2024-03-14T11:05:43.897382Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Filter test dataframes\n\n# Group df_test_person_1 by case_id (that is, collapse rows with same case_id value,\n# which with an additional \"aggregation\" (to be performed later) makes the entries of\n# the chosen columns to accomodate lists with the respective collapsed values (or\n# operations on them, as wished))\ndf_test_person_1_group = df_test_person_1.group_by(\"case_id\")\n\n\n# Aggregate to the group, the maximum income (from main occupation) of the set of people\n# associated with each group, as well as a boolean that holds true if any of the set of\n# people of each group is self-employed\ndf_test_person_1_1 = df_test_person_1_group.agg(\n    pl.col(\"mainoccupationinc_384A\").max().alias(\"mainoccupationinc_max_A\"),\n    (pl.col(\"incometype_1044T\") == \"SELFEMPLOYED\").any().alias(\"anyselfemployed_T\")\n)\n\n# Select from client df_test_person_1 columns \"case_id\", \"num_group1\" and\n# \"housetype_905L\". Keep solely rows associated with applicants (num_group1 = 0), drop\n# the column \"num_group1\" and rename \"housetype_905L\" to \"housetype_applicant_L\", since\n# it refers now to the house type (e.g. \"owned\") of the applicants.\ndf_test_person_1_2 = df_test_person_1.select(\n    [\"case_id\", \"num_group1\", \"housetype_905L\"]\n).filter(pl.col(\"num_group1\") == 0).drop(\"num_group1\").rename(\n    {\"housetype_905L\": \"housetype_applicant_L\"}\n)\n\n# Group df_test_credit_bureau_b_2 by case_id and aggregate maximum value of the number\n# of overdue payments (pmts_pmtsoverdue_635A) and a boolean that holds true if there is\n# any number of days of overdue payment greater than 31.\ndf_test_credit_bureau_b_2 = df_test_credit_bureau_b_2.group_by(\"case_id\").agg(\n    pl.col(\"pmts_pmtsoverdue_635A\").max().alias(\"pmts_pmtsoverdue_max_A\"),\n    (pl.col(\"pmts_dpdvalue_108P\") > 31).any().alias(\"pmts_dpdvalue_anyover31_P\")\n)\n\n# Submission dataframe defined by a left join of all test dataframes.\n# [NOTE: the label \"submission\" is used in place of \"test\" as another test dataset is to\n# be created from the training one.]\ndf_submission = df_test_base\\\n.join(df_test_static.select([\"case_id\"] + cols_static_selec),\n      how=\"left\",\n      on=\"case_id\")\\\n.join(df_test_static_cb.select([\"case_id\"] + cols_static_cb_selec),\n      how=\"left\",\n      on=\"case_id\")\\\n.join(df_test_person_1_1,\n      how=\"left\",\n      on=\"case_id\")\\\n.join(df_test_person_1_2,\n      how=\"left\",\n      on=\"case_id\")\\\n.join(df_test_credit_bureau_b_2,\n      how=\"left\",\n      on=\"case_id\")","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:43.899931Z","iopub.execute_input":"2024-03-14T11:05:43.900344Z","iopub.status.idle":"2024-03-14T11:05:43.921177Z","shell.execute_reply.started":"2024-03-14T11:05:43.900309Z","shell.execute_reply":"2024-03-14T11:05:43.919935Z"},"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\ndf_train.head()","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:43.923481Z","iopub.execute_input":"2024-03-14T11:05:43.923947Z","iopub.status.idle":"2024-03-14T11:05:43.954872Z","shell.execute_reply.started":"2024-03-14T11:05:43.923909Z","shell.execute_reply":"2024-03-14T11:05:43.953719Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Display first five entries of the test dataframe\ndf_submission.head()","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:43.959267Z","iopub.execute_input":"2024-03-14T11:05:43.959979Z","iopub.status.idle":"2024-03-14T11:05:43.976746Z","shell.execute_reply.started":"2024-03-14T11:05:43.959933Z","shell.execute_reply":"2024-03-14T11:05:43.975484Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"***","metadata":{}},{"cell_type":"markdown","source":"# Define training, validation and test datasets from the original training one\n\nOne may consider a validation dataset for performing [early stopping](https://en.wikipedia.org/wiki/Early_stopping) when training the model - the score on the validation dataset is evaluated and if not improving within some tolerance the training is stopped. If early stopping was not considered, overfitting could happen, and even though the model could predict quite well the training data, it might not for the case of the validation and test ones.\n\nSince the labels of the validation dataset would need to be known, the latter would need to be obtained from the original training dataset.\n\nAlso, [@danielherman](https://www.kaggle.com/jetakow) proposes considering an additional test dataset, obtained from the original training dataset that allows one to test the model even before submitting. Note that this test dataset would differ from the validation one since it would not participate in any way in the training of the model. Note also, that one should henceforward refer to the original test dataset as \"submission dataset\" to avoid ambiguity with this derived dataset.\n\nThe considered datasets, their origin, and split fractions are described in the table below.\n\n| Dataset     | Taken from                | Split fraction |\n| ----------- | ------------------------- | -------------- |\n| Training    | Original training dataset | $0.6$          |\n| Validation  | Original training dataset | $0.2$          |\n| Test        | Original training dataset | $0.2$          |\n| Submission  | Original test dataset     | $1$            |","metadata":{}},{"cell_type":"code","source":"# ---> Info on the split of the original training dataset\n\n# Dictionary of parameters for splitting the original training dataset into a new\n# training, validation and test datasets\ntrain_valid_test_split = {\n    # Fraction of the original training dataset to be used on validation\n    \"valid_frac\": 0.2,\n    # Fraction of the original training dataset to be used on testing\n    \"test_frac\": 0.2\n}","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:43.978121Z","iopub.execute_input":"2024-03-14T11:05:43.978861Z","iopub.status.idle":"2024-03-14T11:05:43.994694Z","shell.execute_reply.started":"2024-03-14T11:05:43.978820Z","shell.execute_reply":"2024-03-14T11:05:43.993064Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Define polars series of case_id for the training, validation, test and submission\n# dataframes\n\n# polars series of unique values of case_id in the original training dataframe, suffle\n# with seed 42\n# [NOTE: unique values are taken in order to disregard possible duplicates.]\n# [NOTE: shuffle is done to randomly distribute the training, validation and test\n# datasets at an upcoming moment.]\n# [NOTE: by issuing a seed number to the shuffle random process, later calls of such\n# process would produce the same result, ensuring reproducibility of this work.]\ncase_id_train_old = df_train[\"case_id\"].unique().shuffle(seed=42)\n\n# Number of rows in the original training dataset\nN_train_old = len(case_id_train_old)\n\n# Number of rows in the new original training dataset\nN_train = int((1 - train_valid_test_split[\"valid_frac\"] -\n               train_valid_test_split[\"test_frac\"]) * N_train_old)\n\n# Number of rows in the new validation dataset\nN_valid = int(train_valid_test_split[\"valid_frac\"] * N_train_old)\n\n# polars series of case_id in the new training dataframe\ncase_id_train = case_id_train_old.head(N_train)\n\n# polars series of case_id in the new validation dataframe\ncase_id_valid = case_id_train_old.tail(-N_train).head(N_valid)\n\n# polars series of case_id in the new test dataframe\ncase_id_test = case_id_train_old.tail(-N_train).tail(-N_valid)\n\n# polars series of case_id in the submission dataframe\ncase_id_submission = df_submission[\"case_id\"].unique()","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:43.996587Z","iopub.execute_input":"2024-03-14T11:05:43.997027Z","iopub.status.idle":"2024-03-14T11:05:44.151719Z","shell.execute_reply.started":"2024-03-14T11:05:43.996991Z","shell.execute_reply":"2024-03-14T11:05:44.149917Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 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 df_train.columns:\n    if col[-1].isupper() and col[:-1].islower():\n        cols_x.append(col)","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:44.153466Z","iopub.execute_input":"2024-03-14T11:05:44.153887Z","iopub.status.idle":"2024-03-14T11:05:44.161527Z","shell.execute_reply.started":"2024-03-14T11:05:44.153854Z","shell.execute_reply":"2024-03-14T11:05:44.159759Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Define auxiliary functions\n\n# Function that returns a dictionary of pandas dataframes from given polars dataframe\n# and list of case indices\n# [NOTE: lightgbm (the package used for training the model) only supports pandas\n# dataframes, hence the need to convert the polars dataframes to pandas'.]\ndef from_polars_to_pandas(df, case_id, cols_x):\n    # polars dataframe corresponding to df for the given case_id\n    df_case_id = df.filter(pl.col(\"case_id\").is_in(case_id))\n    \n    # Dictionary of pandas dataframes\n    dt = {\n        # Base dataframe with case_id and columns that are required for the computation\n        # of the Gini coefficient\n        \"base\": df_case_id[[\"case_id\", \"WEEK_NUM\"]].to_pandas(),\n        # Features' dataframe\n        \"x\": df_case_id[cols_x].to_pandas()\n    }\n    \n    # If the polars dataframe has labels (column \"target\") add them to the dictioanry\n    if \"target\" in df_case_id.columns:\n        dt[\"base\"].insert(loc=dt[\"base\"].shape[1], column=\"y\", value=df_case_id[\"target\"])\n        dt[\"y\"] = df_case_id[\"target\"].to_pandas()\n\n    return dt\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 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","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:44.164053Z","iopub.execute_input":"2024-03-14T11:05:44.164495Z","iopub.status.idle":"2024-03-14T11:05:44.188262Z","shell.execute_reply.started":"2024-03-14T11:05:44.164461Z","shell.execute_reply":"2024-03-14T11:05:44.186492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Define dictionaries of pandas dataframes associated with the training,\n# validation, test and submission data\n\n# Dictionary of pandas dataframes associated with the training data\ndt_train = from_polars_to_pandas(df_train, case_id_train, cols_x)\n# Dictionary of pandas dataframes associated with the validation data\ndt_valid = from_polars_to_pandas(df_train, case_id_valid, cols_x)\n# Dictionary of pandas dataframes associated with the test data\ndt_test = from_polars_to_pandas(df_train, case_id_test, cols_x)\n# Dictionary of pandas dataframes associated with the submission data\ndt_submission = from_polars_to_pandas(df_submission, case_id_submission, cols_x)\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_train[\"x\"], dt_valid[\"x\"], dt_test[\"x\"], dt_submission[\"x\"]) = convert_cols_obj_to_cols_cat(\n    dt_train[\"x\"], dt_valid[\"x\"], dt_test[\"x\"], dt_submission[\"x\"]\n)\n# For compatibility reasons, make categories of the feature validation, test and\n# submission dataframes which do not pertain to the the training dataframe be replaced\n# by the category \"Unknown\"\ndt_valid[\"x\"] = make_cat_excl_unknown(df=dt_valid[\"x\"], df_ref=dt_train[\"x\"])\ndt_test[\"x\"] = make_cat_excl_unknown(df=dt_test[\"x\"], df_ref=dt_train[\"x\"])\ndt_submission[\"x\"] = make_cat_excl_unknown(df=dt_submission[\"x\"], df_ref=dt_train[\"x\"])\n\n# Display shapes of feature pandas dataframes\ndisplay(pd.DataFrame(data=\n                     {\"Feature dataset\": [\"train\", \"valid\", \"test\", \"submission\"],\n                      \"N_rows\": [dt[\"x\"].shape[0]for\n                                          dt in (dt_train, dt_valid, dt_test, dt_submission)],\n                      \"N_cols\": [dt[\"x\"].shape[1]for\n                                          dt in (dt_train, dt_valid, dt_test, dt_submission)]\n                     }))","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:44.190195Z","iopub.execute_input":"2024-03-14T11:05:44.190686Z","iopub.status.idle":"2024-03-14T11:05:46.246584Z","shell.execute_reply.started":"2024-03-14T11:05:44.190637Z","shell.execute_reply":"2024-03-14T11:05:46.245281Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"***","metadata":{}},{"cell_type":"markdown","source":"# Training LightGBM\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\n","metadata":{}},{"cell_type":"code","source":"# LightGBM's training dataset\nds_train = lgb.Dataset(\n    data=dt_train[\"x\"],\n    label=dt_train[\"y\"]\n)\n\n# LightGBM's validation dataset\n# [NOTE: by setting reference=ds_train, the defined dataset is regarded as a validation\n# one for ds_train.]\nds_valid = lgb.Dataset(\n    data=dt_valid[\"x\"],\n    label=dt_valid[\"y\"],\n    reference=ds_train\n)\n\n# Dictionary of parameters for training the LightGBM model\nparams = {\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    # [NOTE: \"binary\" stands for binary classification.]\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    # [NOTE: \"cross_entropy\" is the cross-entropy cost function.]\n    \"metric\": [\"auc\", \"cross_entropy\"],\n    # Max depth of the weak tree learners\n    \"max_depth\": 3,\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\": 0.9,\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\": 0.8,\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\": 5,\n    # Maximum number of weak learners (the same as num_iterations)\n    \"n_estimators\": 1000,\n    # Verbosity level of LightGBM\n    # [NOTE: \"-1\" means solely \"fatal errors\".]\n    \"verbose\": -1\n}\n\n\n# Dictionary of evaluation results of the model (to be used later for plotting)\nmodel_eval_result = {}\n\n# Trained Booster model\nmodel = lgb.train(\n    params=params,\n    train_set=ds_train,\n    # List of LightGBM datasets whose points are used to evaluate the model while\n    # training\n    # [NOTE: since ds_train was already used in the argument of train_set, lgb.train\n    # will recognise it as the training dataset and early stopping would only be applied\n    # on the validation one, ds_valid. Anyway, if the user wants to get lgb.train to\n    # compute the metrics for the training dataset, it needs to be included in the\n    # valid_sets argument.]\n    valid_sets=[ds_valid, ds_train],\n    # List of names associated with the list valid_sets\n    valid_names=[\"valid\", \"train\"],\n    # List of callback functions that are applied at each iteration\n    callbacks=[\n        # Report evaluation results at each 50 iterations\n        lgb.log_evaluation(period=50),\n        # Enable early stopping (the model will train until the validation score doesn’t\n        # improve by at least min_delta in the last stopping_rounds-1 iterations and in\n        # the current one)\n        lgb.early_stopping(\n            # Boolean that if set to True makes the trainer to solely use the first\n            # metric in the list of metrics for performing early stopping\n            first_metric_only=True,\n            stopping_rounds=10,\n            verbose=True,\n            min_delta=0),\n        # Save evaluation results into a dictionary\n        lgb.record_evaluation(eval_result=model_eval_result)\n    ]\n)","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:05:46.250378Z","iopub.execute_input":"2024-03-14T11:05:46.251433Z","iopub.status.idle":"2024-03-14T11:07:31.555383Z","shell.execute_reply.started":"2024-03-14T11:05:46.251354Z","shell.execute_reply":"2024-03-14T11:07:31.553465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"***","metadata":{}},{"cell_type":"markdown","source":"# Model Performance Assessment\n\n## Evolution curves of the metrics while training","metadata":{}},{"cell_type":"code","source":"# ---> Get numbers of trees\n\n# List of numbers of trees considered while training\nN_trees = list(range(1, len(model_eval_result[\"train\"][\"cross_entropy\"]) + 1))\n\n# Best number of trees (that is, it is the one that gives the highest validation AUC)\nN_trees_best = model.best_iteration\n\n# Index of the list of numbers of trees for which the respective entry is the best\ni_N_trees_best = model.best_iteration - 1","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:07:31.557111Z","iopub.execute_input":"2024-03-14T11:07:31.557579Z","iopub.status.idle":"2024-03-14T11:07:31.566036Z","shell.execute_reply.started":"2024-03-14T11:07:31.557537Z","shell.execute_reply":"2024-03-14T11:07:31.564322Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Plot cross-entropy cost and AUC curves obtained while training\n\nfor (key_metric, key_metric_display) in [(\"cross_entropy\", \"Cross-entropy cost\"),\n                                         (\"auc\", \"$\\mathrm{AUC}$\")]:\n    # Initialise figure and axes\n    plt.figure(figsize=(6.4, 4.8))\n    ax = plt.axes()\n    plt.title(rf\"{key_metric_display} curves\")\n\n    # Plot training metric                  \n    ax.plot(N_trees,\n            model_eval_result[\"train\"][key_metric],\n            linestyle=\"solid\",\n            color=\"blue\",\n            alpha=1,\n            linewidth=1,\n            marker=\"None\",\n            markeredgecolor=\"blue\",\n            markerfacecolor=\"None\",\n            markersize=3,\n            label=f\"Training curve\")\n\n    # Plot validation metric                   \n    ax.plot(N_trees,\n            model_eval_result[\"valid\"][key_metric],\n            linestyle=\"solid\",\n            color=\"red\",\n            alpha=1,\n            linewidth=1,\n            marker=\"None\",\n            markeredgecolor=\"red\",\n            markerfacecolor=\"None\",\n            markersize=3,\n            label=\"Validation curve\")\n\n    # Plot training metric at best validation AUC\n    ax.plot(N_trees_best,\n            model_eval_result[\"train\"][key_metric][i_N_trees_best],\n            linestyle=\"None\",\n            color=\"blue\",\n            alpha=1,\n            linewidth=1,\n            marker=\"o\",\n            markeredgecolor=\"blue\",\n            markerfacecolor=\"None\",\n            markersize=5,\n            label=\"Training point at best validation AUC\")\n\n    # Plot validation metric at best validation AUC              \n    ax.plot(N_trees_best,\n            model_eval_result[\"valid\"][key_metric][i_N_trees_best],\n            linestyle=\"None\",\n            color=\"red\",\n            alpha=1,\n            linewidth=1,\n            marker=\"o\",\n            markeredgecolor=\"red\",\n            markerfacecolor=\"None\",\n            markersize=5,\n            label=\"Validation point at best validation AUC\")\n\n    # Plot vertical line x = N_trees_best\n    ax.axvline(x=N_trees_best,\n               color=\"black\", linestyle=\"dashed\", linewidth=0.75, zorder=0)\n\n    # Define axes labels                                \n    ax.set_xlabel(r\"Number of trees\", fontdict={\"fontsize\": 10})\n    ax.set_ylabel(rf\"{key_metric_display}\", fontdict={\"fontsize\": 10})\n\n    # Enable axes' minor ticks\n    ax.minorticks_on()\n\n    # Define grid\n    ax.grid(visible=True, which=\"major\", color=\"lightgray\", linestyle=\"solid\",\n            linewidth=0.5)\n    ax.grid(visible=True, which=\"minor\", color=\"lightgray\", linestyle=\"dotted\",\n            linewidth=0.5)\n\n    # Legend\n    ax.legend(fontsize=8)\n\n    # Show plot\n    plt.show() ","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:07:31.567751Z","iopub.execute_input":"2024-03-14T11:07:31.568162Z","iopub.status.idle":"2024-03-14T11:07:32.923340Z","shell.execute_reply.started":"2024-03-14T11:07:31.568129Z","shell.execute_reply":"2024-03-14T11:07:32.921862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Display number of trees at best validation AUC\nprint()\ndisplay(pd.DataFrame(data={\"N_trees_best\": N_trees_best},\n                     index=[0]).style.set_caption(\"Number of trees at best validation AUC\")\\\n        .set_table_styles(\n            [\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-03-14T11:07:32.925334Z","iopub.execute_input":"2024-03-14T11:07:32.926215Z","iopub.status.idle":"2024-03-14T11:07:33.014917Z","shell.execute_reply.started":"2024-03-14T11:07:32.926171Z","shell.execute_reply":"2024-03-14T11:07:33.013620Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## The ROC curve, the area under it (AUC) and the Gini coefficient\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 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\n### The Gini coefficient\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 accuracies, ROC curves, AUCs and Gini coefficients (at best validation\n# AUC) for the training, validation and test sets\n\nfor dt in (dt_train, dt_valid, dt_test):\n    # Add predicted probabilities for label corresponding to 1 (credit default case)\n    # to base pandas dataframe of the dataset dictionary\n    # [NOTE: model.predict() returns a numpy.ndarray.]\n    dt[\"base\"][\"P_pred\"] = model.predict(\n        data=dt[\"x\"],\n        # Number of iterations (trees) to be used\n        # [NOTE: model.best_iteration returns the best found number of iterations.]\n        num_iteration=model.best_iteration\n    )\n    # Add predicted labels for the case of a threshold of 0.5\n    dt[\"base\"][\"y_pred\"] = (dt[\"base\"][\"P_pred\"] >= 0.5).astype(dtype=\"int32\")\n    \n    # Add accuracy to the dataset dictionary\n    dt[\"accuracy\"] = accuracy_score(\n        y_true=dt[\"y\"],\n        y_pred=dt[\"base\"][\"y_pred\"]\n    )\n    \n    # Add ROC's FPR and TPR to the dataset dictionary\n    (dt[\"FPR\"], dt[\"TPR\"], _) = roc_curve(\n        y_true=dt[\"y\"],\n        y_score=dt[\"base\"][\"P_pred\"]\n    )\n\n    # Add AUC to the dataset dictionary\n    dt[\"auc\"] = roc_auc_score(\n        y_true=dt[\"y\"],\n        y_score=dt[\"base\"][\"P_pred\"]\n    )\n    # Add Gini coefficient to the dataset dictionary\n    dt[\"g\"] = 2 * dt[\"auc\"] - 1","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:07:33.016661Z","iopub.execute_input":"2024-03-14T11:07:33.017223Z","iopub.status.idle":"2024-03-14T11:07:52.609368Z","shell.execute_reply.started":"2024-03-14T11:07:33.017171Z","shell.execute_reply":"2024-03-14T11:07:52.607933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Plot ROC curves at best validation AUC\n\nfor (dt, label_dt_display) in [(dt_train, \"Training\"),\n                               (dt_valid, \"Validation\"),\n                               (dt_test, \"Test\")]:\n    \n    # Initialise figure and axes\n    plt.figure(figsize=(4.8, 4.8))\n    ax = plt.axes()\n    plt.title(rf\"{label_dt_display}\")\n\n    # Plot model's ROC\n    ax.plot(dt[\"FPR\"],\n            dt[\"TPR\"],\n            linestyle=\"solid\",\n            color=\"black\",\n            alpha=1,\n            linewidth=1.5,\n            label=\"Model's ROC\")\n    \n    # Plot model's AUC\n    ax.fill_between(\n        x=dt[\"FPR\"],\n        y1=np.zeros(len(dt[\"FPR\"])),\n        y2=dt[\"TPR\"],\n        edgecolor=\"blue\",\n        facecolor=\"blue\",\n        linewidth=0,\n        # hatch = \"///\",\n        alpha=0.5,\n        label=(r\"Model's $\\mathrm{AUC}$ ($\\mathrm{AUC} = \" +\n               rf\"{dt['auc']:.3f}$)\"))\n\n    # Plot model's 1/2 Gini\n    ax.fill_between(\n        x=dt[\"FPR\"],\n        y1=dt[\"FPR\"],\n        y2=dt[\"TPR\"],\n        edgecolor=\"yellow\",\n        facecolor=\"None\",\n        linewidth=0,\n        hatch = \"xxx\",\n        alpha=1,\n        label=(r\"Model's $1/2$ Gini ($G = \" +\n               rf\"{dt['g']:.3f}$)\"))\n\n    # Plot random classifier's ROC\n    ax.plot(\n        [0, 1],\n        [0, 1],\n        linestyle=\"dashed\",\n        color=\"black\",\n        alpha=1,\n        linewidth=1.25,\n        label=\"Random classifier's ROC\")\n\n    # Define axes labels                                \n    ax.set_xlabel(r\"$\\mathrm{FPR}$\", fontdict={\"fontsize\": 10})\n    ax.set_ylabel(r\"$\\mathrm{TPR}$\", fontdict={\"fontsize\": 10})\n\n    # Define axes limits\n    ax.set_xlim([0, 1])\n    ax.set_ylim([0, 1])\n\n    # Set aspect ratio\n    # [NOTE: \"equal\" means that the scales for x and y-axes are equals.]\n    ax.set_aspect(\"equal\")\n    \n    # Legend\n    ax.legend(loc=\"upper left\", fontsize=8)\n\n    # Show plot\n    plt.show() ","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:07:52.611345Z","iopub.execute_input":"2024-03-14T11:07:52.611919Z","iopub.status.idle":"2024-03-14T11:07:53.945233Z","shell.execute_reply.started":"2024-03-14T11:07:52.611866Z","shell.execute_reply":"2024-03-14T11:07:53.943907Z"},"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\n \n# Function for computing stability score's elements\ndef get_stability_score(\n    # Base pandas dataframe of the dataset dictionary\n    dt_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 = dt_base[[\"WEEK_NUM\", \"y\", \"P_pred\"]]\\\n        .sort_values(by=\"WEEK_NUM\")\\\n        .groupby(by=\"WEEK_NUM\")[[\"y\", \"P_pred\"]]\\\n        .apply(lambda x:\n               2 * roc_auc_score(x[\"y\"], x[\"P_pred\"]) - 1).tolist()\n    \n    # Average (in week number) Gini coefficient\n    G_av = np.mean(G)\n\n    # Array of indices for the Gini coefficients\n    i = np.arange(len(G))\n    \n    # Weight (a) and bias (_) of the linear regression\n    [a, b] = np.polyfit(x=i, y=G, 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 = np.sqrt(np.mean((G_fit - G)**2))\n\n    # Stability score\n    stability_score = w_G_av * G_av + w_a * min(0, a) + w_RMSD * RMSD\n    \n    # Dictionary of stability score elements\n    dt = {\n        \"g_week\": G,\n        \"a\": a,\n        \"b\": b,\n        \"RMSD\": RMSD,\n        \"stability_score\": stability_score\n    }\n    \n    return dt\n\nfor dt in (dt_train, dt_valid, dt_test):\n    # Add stability score elemements to the dataset dictionary\n    dt.update(get_stability_score(dt[\"base\"]))\n","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:07:53.946692Z","iopub.execute_input":"2024-03-14T11:07:53.947658Z","iopub.status.idle":"2024-03-14T11:07:55.274848Z","shell.execute_reply.started":"2024-03-14T11:07:53.947615Z","shell.execute_reply":"2024-03-14T11:07:55.273206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Plot stability scores' elements\n\nfor (dt, label_dt_display) in [(dt_train, \"Training\"),\n                               (dt_valid, \"Validation\"),\n                               (dt_test, \"Test\")]:\n    \n    # Initialise figure and axes\n    plt.figure(figsize=(6.4, 4.8))\n    ax = plt.axes()\n    plt.title(rf\"{label_dt_display}\")\n\n    # Plot Gini coefficients\n    ax.plot(range(len(dt[\"g_week\"])),\n            dt[\"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\n    i = np.array(range(len(dt[\"g_week\"])))[[0, -1]]\n    ax.plot(i,\n            dt[\"a\"] * i + dt[\"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['a']:.2e})\"))\n\n    # Define axes labels                                \n    ax.set_xlabel(r\"Week number, $i$\", fontdict={\"fontsize\": 10})\n    ax.set_ylabel(r\"Gini coefficient, $G$\", fontdict={\"fontsize\": 10})\n\n    # Enable axes' minor ticks\n    ax.minorticks_on()\n\n    # Define grid\n    ax.grid(visible=True, which=\"major\", color=\"lightgray\", linestyle=\"solid\",\n            linewidth=0.5)\n    ax.grid(visible=True, which=\"minor\", color=\"lightgray\", linestyle=\"dotted\",\n            linewidth=0.5)\n\n    # Legend\n    ax.legend(fontsize=8)\n\n    # Show plot\n    plt.show() ","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:07:55.277288Z","iopub.execute_input":"2024-03-14T11:07:55.278446Z","iopub.status.idle":"2024-03-14T11:07:57.155027Z","shell.execute_reply.started":"2024-03-14T11:07:55.278369Z","shell.execute_reply":"2024-03-14T11:07:57.154015Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Metric summary table","metadata":{}},{"cell_type":"code","source":"# Summary dictionary for display\ndt_summary = {\"dataset\": [\"train\", \"valid\", \"test\"]}\ndt_summary.update({key_metric: [dt[key_metric] for dt in\n                                (dt_train, dt_valid, dt_test)] for\n                   key_metric in [\"accuracy\", \"auc\", \"g\", \"stability_score\"]})\n\n# Display summary\nprint()\ndisplay(pd.DataFrame(data=dt_summary)\\\n        .style.set_caption(\"Metrics at best validation 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            ]))","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:07:57.156532Z","iopub.execute_input":"2024-03-14T11:07:57.157116Z","iopub.status.idle":"2024-03-14T11:07:57.174149Z","shell.execute_reply.started":"2024-03-14T11:07:57.157082Z","shell.execute_reply":"2024-03-14T11:07:57.172794Z"},"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 submission dataset dictionary\ndt_submission[\"base\"][\"P_pred\"] = model.predict(\n    data=dt_submission[\"x\"],\n    # Number of iterations (trees) to be used\n    # [NOTE: model.best_iteration returns the best found number of iterations.]\n    num_iteration=model.best_iteration\n)","metadata":{"execution":{"iopub.status.busy":"2024-03-14T11:07:57.175632Z","iopub.execute_input":"2024-03-14T11:07:57.176894Z","iopub.status.idle":"2024-03-14T11:07:57.212890Z","shell.execute_reply.started":"2024-03-14T11:07:57.176846Z","shell.execute_reply":"2024-03-14T11:07:57.211510Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---> Create submission pandas dataframe\n# Submission dataframe\nsubmission = pd.DataFrame({\n    \"case_id\": dt_submission[\"base\"][\"case_id\"].to_numpy(),\n    \"score\": dt_submission[\"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-03-14T11:07:57.214301Z","iopub.execute_input":"2024-03-14T11:07:57.214675Z","iopub.status.idle":"2024-03-14T11:07:57.233237Z","shell.execute_reply.started":"2024-03-14T11:07:57.214644Z","shell.execute_reply":"2024-03-14T11:07:57.231910Z"},"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-03-14T11:07:57.234563Z","iopub.execute_input":"2024-03-14T11:07:57.235046Z","iopub.status.idle":"2024-03-14T11:07:57.250571Z","shell.execute_reply.started":"2024-03-14T11:07:57.235004Z","shell.execute_reply":"2024-03-14T11:07:57.249157Z"},"trusted":true},"execution_count":null,"outputs":[]}]}