{"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":30673,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"##
Population stability index for credit scoring models
\n \nThe main aim of this notebook is to provide a basic understanding of the population stability index (**PSI**) and its potential applications in monitoring the model's stability.\n \n**Key observations**:\n\n* PSI score can be used to evaluate whatever training features are stable over time examples are shown in [7.1 Numerical - low number of unique values](#7.1), [7.2 Numerical - more than 10 unique values](#7.2) and - [7.3 Categorical - low number of unique values](#7.3) sub-sections.\n* PSI score can be also used as an additional validation metric for the scoring model (see [9.2 PSI validation](#9.2) and [9.3 Risk differentiation](#9.3) sub-sections).\n* **Experiment #1** - trained log. regression model using only 5 low weekly PSI numerical features without missing values and had which had less than 10 unique values (**nb v5**). Models's weekly GINI score hovered around 0.115 on the validation sample. A strong negative correlation was observed between weekly GINI score and number of observations, i.e.weeks with a small number of observations tend to have large GINI score variance. Linear curve fit `y = k * a + b` on weekly GINI score was very small `k = 0.001`. The public **LB** score for this model was 0.101.\n* **Experiment #2** - trained LightGBM model with additional numerical features with high weekly PSI (**nb v6**). The average GINI score on the validation sample increased from 0.12 to 0.185 (~55% improved as compared to training with 5 features from **Experiment #1**). A similar improvement of the public **LB** score of 0.16 was observed when compared to the previous experiment. Interestingly, adding high PSI features did not result in increased linear dependence on weekly GINI score and increased standard deviation of GINI score (compare [9. Validation](#9) charts between different notebook versions). \n* **Experiment #3** - numerical features without missing values were cut into 10 equal buckets using the `qcut` function from the pandas package. All 5 features were rather stable as the maximum weekly PSI score was around 0.1. The average weekly GINI score error distribution has a larger mean of ~50% as compared to previous experiments and a relatively smaller variance. While validation score increased, the public **LB** score was sligthly smaller comapred to previous experiment.\n* **Experiment #4** - trained model on stable features selected in experiments 1 and 3. The standard deviation of average weekly GINI scores decreased compared to previous experiments. Combining features from two experiments resulted in increased **LB** score of 0.214.\n* **Experiment #5** - trained model on all features investigated with PSI. Using features from all experiments resulted in increased **LB** score of 0.243.\n* **Experiment #6** - trained model on 4 categorical features with less than 10 unique values. One of those features had stable weekly PSI. The average cross-validation GINI score was similiar to **Experiment #2**.\n* **Experiment #7** - trained model on categorical features with more than 10 unique values. Predictions could not be split into 5 equal buckets. In addition, risk differentiation is also poor for this model (see [9.3 Risk differentiation](#9.3)). As expected **LB** score was very low (0.094).\n* **Experiment #8** - trained model on all features which were investigated using PSI in section [7. Features (no NaN values)](#7). The average ROC AUC score on the validation sample increased to 0.7 during ensembling training.\n \n \n##
Table of contents
\n\n- [1. Prepare project](#1)\n- [2. Theory](#2)\n- [3. Helper functions](#3)\n- [4. Load data](#4)\n- [5. Feature engineering](#5)\n- [6. Describe features](#6) \n - [6.1 Numerical features](#6.1)\n - [6.2 Categorical features](#6.2) \n- [7. Features (no NaN values)](#7)\n - [7.1 Numerical - low number of unique values](#7.1)\n - [7.2 Numerical - more than 10 unique values](#7.2)\n - [7.3 Categorical - low number of unique values](#7.3)\n - [7.4 Categorical - more than 10 unique values](#7.4)\n- [8. Modeling](#8)\n- [9. Validation](#9)\n - [9.1 Classical](#9.1)\n - [9.2 PSI validation](#9.2)\n - [9.3 Risk differentiation](#9.3)\n- [10. Submit predictions](#10)\n- [11. Changelog](#11)\n \n\n# 1. Prepare project","metadata":{}},{"cell_type":"code","source":"# suppress enjoying warnings from seaborn\nimport warnings\nwarnings.simplefilter(action='ignore', category=FutureWarning)\n\n# for fiding file names\nfrom pathlib import Path\nfrom glob import glob\nimport gc\n\n# data processing libraries\nimport polars as pl\nimport numpy as np\nimport pandas as pd\n\n# for visualizing data\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\n# sklearn functions\nfrom sklearn.model_selection import train_test_split, StratifiedGroupKFold\nfrom sklearn.preprocessing import LabelEncoder\n# for validating models\nfrom sklearn.metrics import roc_auc_score, roc_curve\n\n# LightGBM modeling\nimport lightgbm as lgb\n\n# define default colors for plots in notebook\nfrom matplotlib import cycler\nfrom matplotlib.colors import LinearSegmentedColormap\nCOLORS = [\"#068D9D\", \"#53599A\", \"#607BB0\", \"#6D9DC5\", \"#77BECF\", \"#80DED9\", \"#AEECEF\"]\nplt.rc('axes', facecolor='#E6E6E6', edgecolor='none', axisbelow=True, grid=True, prop_cycle=cycler('color', COLORS))\n\n# project CONSTANTS\nROOT = Path(\"/kaggle/input/home-credit-credit-risk-model-stability\")\nTRAIN_DIR = ROOT / \"parquet_files\" / \"train\"\nTEST_DIR = ROOT / \"parquet_files\" / \"test\"\nSEED = 42","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:00.271533Z","iopub.execute_input":"2024-04-08T15:08:00.272512Z","iopub.status.idle":"2024-04-08T15:08:00.286409Z","shell.execute_reply.started":"2024-04-08T15:08:00.272456Z","shell.execute_reply":"2024-04-08T15:08:00.285151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# 2. Theory\n\n**PSI** is a vital tool in credit risk model performance monitoring. It ensures the stability of model predictions over time by comparing input data used to build the model and a validation/in-production dataset. A low PSI indicates that validation/in-production data is very similar to the development data. If the **PSI** is high, it suggests that the model’s input data varies significantly between development and tested datasets, indicating instability. This could be a sign that the model needs recalibration or redevelopment. By calculating and monitoring **PSI**:\n\n* we can estimate how stable data is in development data set through time, e.g. weeks/months;\n* we can expect to decrease in the predictive and discriminatory power of the model if input data changes significantly compared to development data.\n\n**PSI** is a statistical measure with a basis in information theory that quantifies the difference between one probability distribution from a reference probability distribution. A discrete form of PSI:\n\n$$PSI = \\sum_{}^{}\\left((Actual,\\% - Expected,\\%) \\times \\ln{\\left(\\frac{Actual,\\%}{Expected,\\%} \\right)} \\right),$$\n\nwhere $Actual,\\%$ and $Expected,\\%$ are the shares of observations in bucket $i$.\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"def PSI_custom(array_x, array_y):\n \"\"\"\n TODO: update docstring\n \n Calculates PSI between two vectors\n Output:\n return pandas DataFrame and PSI score, int\n \"\"\"\n \n # calculate shares\n share_x = pd.DataFrame({\"bucket\": array_x}).value_counts().reset_index()\n share_y = pd.DataFrame({\"bucket\": array_y}).value_counts().reset_index()\n share_x[\"proportion_x\"] = share_x[\"count\"] / share_x[\"count\"].sum() \n share_y[\"proportion_y\"] = share_y[\"count\"] / share_y[\"count\"].sum() \n \n # add small constant of 1e-6 for buckets with no observations in one of the arrays\n df_shares = share_x.merge(share_y, on=\"bucket\", how = \"outer\").fillna(1e-6)\n df_shares['diff'] = df_shares['proportion_x'] - df_shares['proportion_y']\n df_shares['psi'] = (df_shares['proportion_x'] - df_shares['proportion_y']) * np.log(df_shares['proportion_x'] / df_shares['proportion_y'])\n \n # return results\n return df_shares, df_shares['psi'].sum()\n\n\ndef plot_PSI(df, psi_score):\n \"\"\"\n TODO: update docstring\n \n Plot distrubutions in raw counts and in proportions\n \"\"\"\n fig, ax = plt.subplots(1, 2, figsize=(16, 4))\n\n df_left = df[['bucket', 'count_x', 'count_y']].melt(id_vars='bucket', var_name=\"array\", value_name='count')\n df_right = df[['bucket', 'proportion_x', 'proportion_y']].melt(id_vars='bucket', var_name=\"array\", value_name='proportion %')\n df_right['proportion %'] = df_right['proportion %'] * 100\n \n df_left['array'] = df_left['array'].replace({'count_x': \"array x\", 'count_y': \"array y\"})\n df_right['array'] = df_right['array'].replace({'proportion_x': \"array x\", 'proportion_y': \"array y\"})\n\n sns.barplot(df_left, x=\"bucket\", y=\"count\", hue=\"array\", ax=ax[0])\n sns.barplot(df_right, x=\"bucket\", y=\"proportion %\", hue=\"array\", ax=ax[1])\n \n ax[1].set_title(f\"PSI score: {psi_score:.3f}\")\n\n for i in range(2):\n # add white background to legend\n legend = ax[i].legend(frameon=1)\n frame = legend.get_frame()\n frame.set_facecolor('w')\n\n plt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:00.289282Z","iopub.execute_input":"2024-04-08T15:08:00.290384Z","iopub.status.idle":"2024-04-08T15:08:00.311569Z","shell.execute_reply.started":"2024-04-08T15:08:00.290336Z","shell.execute_reply":"2024-04-08T15:08:00.310000Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Below is an example of two distributions with the same number of observations.","metadata":{}},{"cell_type":"code","source":"# example 1\narray_x = [1, 2, 1, 1, 1, 3, 2]\narray_y = [2, 1, 1, 1, 2, 2, 3]\n\ndf_shares, psi_score = PSI_custom(array_x, array_y)\nplot_PSI(df_shares, psi_score)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:00.315141Z","iopub.execute_input":"2024-04-08T15:08:00.315607Z","iopub.status.idle":"2024-04-08T15:08:01.123934Z","shell.execute_reply.started":"2024-04-08T15:08:00.315574Z","shell.execute_reply":"2024-04-08T15:08:01.122760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The PSI score doesn't depend on the number of observations. THe same exaple as above, just the number of observations in array x was increased 5 times.","metadata":{}},{"cell_type":"code","source":"# example 2\narray_x = [1, 2, 1, 1, 1, 3, 2] * 5\narray_y = [2, 1, 1, 1, 2, 2, 3]\n\ndf_shares, psi_score = PSI_custom(array_x, array_y)\nplot_PSI(df_shares, psi_score)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:01.129258Z","iopub.execute_input":"2024-04-08T15:08:01.129916Z","iopub.status.idle":"2024-04-08T15:08:01.780178Z","shell.execute_reply.started":"2024-04-08T15:08:01.129883Z","shell.execute_reply":"2024-04-08T15:08:01.778657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For cases when observation is not found in one of the distributions, zeros are substituted with `1e-6`. Such cases will always result to large PSI scores.","metadata":{}},{"cell_type":"code","source":"# example 3\narray_x = [1, 2, 2, 3, 3, 3] * 3\narray_y = [1, 1, 2, 3, 4, 4]\n\ndf_shares, psi_score = PSI_custom(array_x, array_y)\nplot_PSI(df_shares, psi_score)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:01.782335Z","iopub.execute_input":"2024-04-08T15:08:01.782823Z","iopub.status.idle":"2024-04-08T15:08:02.440213Z","shell.execute_reply.started":"2024-04-08T15:08:01.782780Z","shell.execute_reply":"2024-04-08T15:08:02.439012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# 3. Helper functions\n\nMost of the functions borrowed from the notebook https://www.kaggle.com/code/daviddirethucus/home-credit-risk-lightgbm .\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"class Pipeline:\n \"\"\"\n Helper class taken from notebook:\n https://www.kaggle.com/code/daviddirethucus/home-credit-risk-lightgbm\n \"\"\"\n def set_table_dtypes(df):\n for col in df.columns:\n if col in [\"case_id\", \"WEEK_NUM\", \"num_group1\", \"num_group2\"]:\n df = df.with_columns(pl.col(col).cast(pl.Int64))\n elif col in [\"date_decision\"]:\n df = df.with_columns(pl.col(col).cast(pl.Date))\n elif col[-1] in (\"P\", \"A\"):\n df = df.with_columns(pl.col(col).cast(pl.Float64))\n elif col[-1] in (\"M\",):\n df = df.with_columns(pl.col(col).cast(pl.String))\n elif col[-1] in (\"D\",):\n df = df.with_columns(pl.col(col).cast(pl.Date))\n return df\n\n def handle_dates(df):\n for col in df.columns:\n if col[-1] in (\"D\",):\n df = df.with_columns(pl.col(col) - pl.col(\"date_decision\")) #!!?\n df = df.with_columns(pl.col(col).dt.total_days()) # t - t-1\n df = df.drop(\"date_decision\", \"MONTH\")\n return df\n\n def filter_cols(df):\n for col in df.columns:\n if col not in [\"target\", \"case_id\", \"WEEK_NUM\"]:\n isnull = df[col].is_null().mean()\n if isnull > 0.7:\n df = df.drop(col)\n \n for col in df.columns:\n if (col not in [\"target\", \"case_id\", \"WEEK_NUM\"]) & (df[col].dtype == pl.String):\n freq = df[col].n_unique()\n if (freq == 1) | (freq > 200):\n df = df.drop(col)\n \n return df","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:02.442073Z","iopub.execute_input":"2024-04-08T15:08:02.443337Z","iopub.status.idle":"2024-04-08T15:08:02.460999Z","shell.execute_reply.started":"2024-04-08T15:08:02.443269Z","shell.execute_reply":"2024-04-08T15:08:02.459446Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Aggregator:\n \"\"\"\n Helper class taken from notebook:\n https://www.kaggle.com/code/daviddirethucus/home-credit-risk-lightgbm\n \"\"\" \n def num_expr(df):\n cols = [col for col in df.columns if col[-1] in (\"P\", \"A\")]\n expr_max = [pl.max(col).alias(f\"max_{col}\") for col in cols]\n return expr_max\n \n def date_expr(df):\n cols = [col for col in df.columns if col[-1] in (\"D\")]\n expr_max = [pl.max(col).alias(f\"max_{col}\") for col in cols]\n return expr_max\n \n def str_expr(df):\n cols = [col for col in df.columns if col[-1] in (\"M\",)]\n expr_max = [pl.max(col).alias(f\"max_{col}\") for col in cols]\n return expr_max\n \n def other_expr(df):\n cols = [col for col in df.columns if col[-1] in (\"T\", \"L\")]\n expr_max = [pl.max(col).alias(f\"max_{col}\") for col in cols]\n return expr_max\n \n def count_expr(df):\n cols = [col for col in df.columns if \"num_group\" in col]\n expr_max = [pl.max(col).alias(f\"max_{col}\") for col in cols] # max & replace col name\n return expr_max\n \n def get_exprs(df):\n exprs = Aggregator.num_expr(df) + \\\n Aggregator.date_expr(df) + \\\n Aggregator.str_expr(df) + \\\n Aggregator.other_expr(df) + \\\n Aggregator.count_expr(df)\n\n return exprs","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:02.463416Z","iopub.execute_input":"2024-04-08T15:08:02.464057Z","iopub.status.idle":"2024-04-08T15:08:02.481127Z","shell.execute_reply.started":"2024-04-08T15:08:02.464018Z","shell.execute_reply":"2024-04-08T15:08:02.479445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_file(path, depth=None):\n \"\"\"\n Helper function taken from notebook:\n https://www.kaggle.com/code/daviddirethucus/home-credit-risk-lightgbm\n \"\"\"\n df = pl.read_parquet(path)\n df = df.pipe(Pipeline.set_table_dtypes)\n if depth in [1,2]:\n df = df.group_by(\"case_id\").agg(Aggregator.get_exprs(df)) \n return df","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:02.482997Z","iopub.execute_input":"2024-04-08T15:08:02.483438Z","iopub.status.idle":"2024-04-08T15:08:02.501605Z","shell.execute_reply.started":"2024-04-08T15:08:02.483399Z","shell.execute_reply":"2024-04-08T15:08:02.499728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_files(regex_path, depth=None):\n \"\"\"\n Helper function taken from notebook:\n https://www.kaggle.com/code/daviddirethucus/home-credit-risk-lightgbm\n \"\"\"\n chunks = []\n \n for path in glob(str(regex_path)):\n df = pl.read_parquet(path)\n df = df.pipe(Pipeline.set_table_dtypes)\n if depth in [1, 2]:\n df = df.group_by(\"case_id\").agg(Aggregator.get_exprs(df))\n chunks.append(df)\n \n df = pl.concat(chunks, how=\"vertical_relaxed\")\n df = df.unique(subset=[\"case_id\"])\n return df","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:02.503368Z","iopub.execute_input":"2024-04-08T15:08:02.503798Z","iopub.status.idle":"2024-04-08T15:08:02.517353Z","shell.execute_reply.started":"2024-04-08T15:08:02.503753Z","shell.execute_reply":"2024-04-08T15:08:02.515167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def feature_eng(df_base, depth_0, depth_1, depth_2):\n \"\"\"\n Helper function taken from notebook:\n https://www.kaggle.com/code/daviddirethucus/home-credit-risk-lightgbm\n \"\"\"\n df_base = (\n df_base\n .with_columns(\n month_decision = pl.col(\"date_decision\").dt.month(),\n weekday_decision = pl.col(\"date_decision\").dt.weekday(),\n )\n )\n for i, df in enumerate(depth_0 + depth_1 + depth_2):\n df_base = df_base.join(df, how=\"left\", on=\"case_id\", suffix=f\"_{i}\")\n df_base = df_base.pipe(Pipeline.handle_dates)\n return df_base","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:02.523894Z","iopub.execute_input":"2024-04-08T15:08:02.524911Z","iopub.status.idle":"2024-04-08T15:08:02.535725Z","shell.execute_reply.started":"2024-04-08T15:08:02.524855Z","shell.execute_reply":"2024-04-08T15:08:02.533767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def to_pandas(df_data, cat_cols=None):\n \"\"\"\n Helper function taken from notebook:\n https://www.kaggle.com/code/daviddirethucus/home-credit-risk-lightgbm\n \"\"\"\n df_data = df_data.to_pandas()\n if cat_cols is None:\n cat_cols = list(df_data.select_dtypes(\"object\").columns)\n df_data[cat_cols] = df_data[cat_cols].astype(\"category\")\n return df_data, cat_cols","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:02.537947Z","iopub.execute_input":"2024-04-08T15:08:02.538482Z","iopub.status.idle":"2024-04-08T15:08:02.548034Z","shell.execute_reply.started":"2024-04-08T15:08:02.538447Z","shell.execute_reply":"2024-04-08T15:08:02.546432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def reduce_mem_usage(df, df_name, verbose=True):\n \"\"\"\n This function changes pandas numerical dtypes (reduces bit size if possible)\n to reduce memory usage\n \n :param df: pandas DataFrame\n :param df_name: str, name of DataFrame\n :param verbose: bool, if True prints out message of how much memory usage was reduced\n \n :return: pandas DataFrame \n \"\"\"\n numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n start_mem = df.memory_usage().sum() / 1024**2\n for col in df.columns:\n col_type = df[col].dtypes\n if col_type in numerics:\n c_min = df[col].min()\n c_max = df[col].max()\n if str(col_type)[:3] == 'int':\n if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n df[col] = df[col].astype(np.int16)\n elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n df[col] = df[col].astype(np.int32)\n elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n df[col] = df[col].astype(np.int64)\n elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n df[col] = df[col].astype(np.int64)\n else:\n if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\n df[col] = df[col].astype(np.float32)\n elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n df[col] = df[col].astype(np.float64)\n else:\n df[col] = df[col].astype(np.float64)\n # calculate memory after reduction\n end_mem = df.memory_usage().sum() / 1024**2\n if verbose:\n # reduced memory usage in percent\n diff_pst = 100 * (start_mem - end_mem) / start_mem\n msg = f'{df_name} mem. usage decreased to {end_mem:5.2f} Mb ({diff_pst:.1f}% reduction)'\n print(msg)\n return df","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:02.550096Z","iopub.execute_input":"2024-04-08T15:08:02.550520Z","iopub.status.idle":"2024-04-08T15:08:02.568296Z","shell.execute_reply.started":"2024-04-08T15:08:02.550486Z","shell.execute_reply":"2024-04-08T15:08:02.567049Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# 4. Load data\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"%%time\ndata_train = {\n \"df_base\": read_file(TRAIN_DIR / \"train_base.parquet\"),\n \"depth_0\": [\n read_file(TRAIN_DIR / \"train_static_cb_0.parquet\"),\n read_files(TRAIN_DIR / \"train_static_0_*.parquet\"),\n ],\n \"depth_1\": [\n read_files(TRAIN_DIR / \"train_applprev_1_*.parquet\", 1),\n read_file(TRAIN_DIR / \"train_tax_registry_a_1.parquet\", 1),\n read_file(TRAIN_DIR / \"train_tax_registry_b_1.parquet\", 1),\n read_file(TRAIN_DIR / \"train_tax_registry_c_1.parquet\", 1),\n read_files(TRAIN_DIR / \"train_credit_bureau_a_1_*.parquet\", 1),\n read_file(TRAIN_DIR / \"train_credit_bureau_b_1.parquet\", 1),\n read_file(TRAIN_DIR / \"train_other_1.parquet\", 1),\n read_file(TRAIN_DIR / \"train_person_1.parquet\", 1),\n read_file(TRAIN_DIR / \"train_deposit_1.parquet\", 1),\n read_file(TRAIN_DIR / \"train_debitcard_1.parquet\", 1),\n ],\n \"depth_2\": [\n read_file(TRAIN_DIR / \"train_credit_bureau_b_2.parquet\", 2),\n read_files(TRAIN_DIR / \"train_credit_bureau_a_2_*.parquet\", 2),\n ]\n}","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:08:02.570043Z","iopub.execute_input":"2024-04-08T15:08:02.570521Z","iopub.status.idle":"2024-04-08T15:11:27.821753Z","shell.execute_reply.started":"2024-04-08T15:08:02.570488Z","shell.execute_reply":"2024-04-08T15:11:27.818958Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndata_test = {\n \"df_base\": read_file(TEST_DIR / \"test_base.parquet\"),\n \"depth_0\": [\n read_file(TEST_DIR / \"test_static_cb_0.parquet\"),\n read_files(TEST_DIR / \"test_static_0_*.parquet\"),\n ],\n \"depth_1\": [\n read_files(TEST_DIR / \"test_applprev_1_*.parquet\", 1),\n read_file(TEST_DIR / \"test_tax_registry_a_1.parquet\", 1),\n read_file(TEST_DIR / \"test_tax_registry_b_1.parquet\", 1),\n read_file(TEST_DIR / \"test_tax_registry_c_1.parquet\", 1),\n read_files(TEST_DIR / \"test_credit_bureau_a_1_*.parquet\", 1),\n read_file(TEST_DIR / \"test_credit_bureau_b_1.parquet\", 1),\n read_file(TEST_DIR / \"test_other_1.parquet\", 1),\n read_file(TEST_DIR / \"test_person_1.parquet\", 1),\n read_file(TEST_DIR / \"test_deposit_1.parquet\", 1),\n read_file(TEST_DIR / \"test_debitcard_1.parquet\", 1),\n ],\n \"depth_2\": [\n read_file(TEST_DIR / \"test_credit_bureau_b_2.parquet\", 2),\n read_files(TEST_DIR / \"test_credit_bureau_a_2_*.parquet\", 2),\n ]\n}","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:11:27.825373Z","iopub.execute_input":"2024-04-08T15:11:27.825888Z","iopub.status.idle":"2024-04-08T15:11:28.554483Z","shell.execute_reply.started":"2024-04-08T15:11:27.825850Z","shell.execute_reply":"2024-04-08T15:11:28.552752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get column names of original/raw features\nFEATS_ORIG = []\n\n# get column names\nfor _key in data_train.keys():\n if isinstance(data_train[_key], list):\n for _df in data_train[_key]:\n FEATS_ORIG += _df.columns\n \n# leave only unique values\nFEATS_ORIG = list(set(FEATS_ORIG))\n# drop case_id column\nFEATS_ORIG.remove(\"case_id\")","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:11:28.556616Z","iopub.execute_input":"2024-04-08T15:11:28.557052Z","iopub.status.idle":"2024-04-08T15:11:28.565140Z","shell.execute_reply.started":"2024-04-08T15:11:28.557020Z","shell.execute_reply":"2024-04-08T15:11:28.563814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# 5. Feature engineering\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"%%time\ndf_train = feature_eng(**data_train)\nprint(\"train data shape:\\t\", df_train.shape)\n# clean memory\ndel data_train\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:11:28.568220Z","iopub.execute_input":"2024-04-08T15:11:28.569236Z","iopub.status.idle":"2024-04-08T15:11:43.117510Z","shell.execute_reply.started":"2024-04-08T15:11:28.569179Z","shell.execute_reply":"2024-04-08T15:11:43.116371Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf_test = feature_eng(**data_test)\nprint(\"test data shape:\\t\", df_test.shape)\n# clean memory\ndel data_test\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:11:43.119025Z","iopub.execute_input":"2024-04-08T15:11:43.119450Z","iopub.status.idle":"2024-04-08T15:11:43.286720Z","shell.execute_reply.started":"2024-04-08T15:11:43.119418Z","shell.execute_reply":"2024-04-08T15:11:43.285283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get column names of new created features\nFEATS_NEW = set(df_train.columns)\n# remove speific columns\nFEATS_NEW = FEATS_NEW.difference({'case_id', 'WEEK_NUM', 'target'})\nFEATS_NEW = FEATS_NEW.difference(set(FEATS_ORIG))\nFEATS_NEW = list(FEATS_NEW)\n\nprint(f\"Before feature engineering: {len(FEATS_ORIG)} featues.\")\nprint(f\"{len(FEATS_NEW)} new features created.\")","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:11:43.288525Z","iopub.execute_input":"2024-04-08T15:11:43.289342Z","iopub.status.idle":"2024-04-08T15:11:43.297414Z","shell.execute_reply.started":"2024-04-08T15:11:43.289283Z","shell.execute_reply":"2024-04-08T15:11:43.296365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Drop the insignificant features\ndf_train = df_train.pipe(Pipeline.filter_cols)\ndf_test = df_test.select([col for col in df_train.columns if col != \"target\"])\n\nprint(\"train data shape: \", df_train.shape)\nprint(\"test data shape: \", df_test.shape)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:11:43.298623Z","iopub.execute_input":"2024-04-08T15:11:43.298993Z","iopub.status.idle":"2024-04-08T15:11:46.701496Z","shell.execute_reply.started":"2024-04-08T15:11:43.298962Z","shell.execute_reply":"2024-04-08T15:11:46.700261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get column names of the remaining features\nFEATS_REMAIN = set(df_train.columns)\n# remove speific columns\nFEATS_REMAIN = FEATS_REMAIN.difference({'case_id', 'WEEK_NUM', 'target'})\n# remaining in original features\n_1 = [_ for _ in FEATS_REMAIN if _ in FEATS_ORIG]\n# remaining of the new features\n_2 = [_ for _ in FEATS_REMAIN if _ in FEATS_NEW]\n\nprint(\"After removing insignificant features:\")\nprint(f\"{len(_1)} original features left\")\nprint(f\"{len(_2)} new features left\")","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:11:46.702879Z","iopub.execute_input":"2024-04-08T15:11:46.703243Z","iopub.status.idle":"2024-04-08T15:11:46.712843Z","shell.execute_reply.started":"2024-04-08T15:11:46.703212Z","shell.execute_reply":"2024-04-08T15:11:46.711701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# convert back to pandas\ndf_train, cat_cols = to_pandas(df_train)\ndf_test, cat_cols = to_pandas(df_test, cat_cols)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:11:46.714469Z","iopub.execute_input":"2024-04-08T15:11:46.715633Z","iopub.status.idle":"2024-04-08T15:12:11.074118Z","shell.execute_reply.started":"2024-04-08T15:11:46.715572Z","shell.execute_reply":"2024-04-08T15:12:11.072703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# reduce memory usage if available\ndf_train = reduce_mem_usage(df_train, \"df_train\")\ndf_test = reduce_mem_usage(df_test, \"df_test\")","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:11.076780Z","iopub.execute_input":"2024-04-08T15:12:11.077470Z","iopub.status.idle":"2024-04-08T15:12:18.357000Z","shell.execute_reply.started":"2024-04-08T15:12:11.077420Z","shell.execute_reply":"2024-04-08T15:12:18.355808Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# 6. Describe features\n\n\n## 6.1. Numerical features\n\nFind features with no missing values in both `train` and `test` data sets.\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"%%time\n# select numeric columns\nnum_cols = df_train.iloc[:, 3:].select_dtypes(include='number').columns\n# get column names without NaN values\n_ = df_train[num_cols].isna().sum()\n# get numerical column names without NaN\nnum_cols_no_nan = _[_==0].index.tolist()","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:18.358618Z","iopub.execute_input":"2024-04-08T15:12:18.358974Z","iopub.status.idle":"2024-04-08T15:12:28.181012Z","shell.execute_reply.started":"2024-04-08T15:12:18.358943Z","shell.execute_reply":"2024-04-08T15:12:28.179378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get train data set stats\n_df_1 = df_train[num_cols_no_nan].describe().T\n_df_1.drop(columns=['count'], inplace=True)\n_df_1.columns = [f\"train {_}\" for _ in _df_1.columns]\n# get test data set stats\n_df_2 = df_test[num_cols_no_nan].describe().T\n_df_2.drop(columns=['count'], inplace=True)\n_df_2.columns = [f\"test {_}\" for _ in _df_2.columns]\npd.merge(_df_1, _df_2, left_index=True, right_index=True, how=\"outer\").sort_values(by=\"train mean\")","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:28.182860Z","iopub.execute_input":"2024-04-08T15:12:28.183321Z","iopub.status.idle":"2024-04-08T15:12:30.783634Z","shell.execute_reply.started":"2024-04-08T15:12:28.183262Z","shell.execute_reply":"2024-04-08T15:12:30.782153Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 6.2. Categorical features\n\nFind features with no missing values in both `train` and `test` data sets.\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"%%time\n# select numeric columns\ncat_cols = df_train.iloc[:, 3:].select_dtypes(include='category').columns\n# get column names without NaN values\n_ = df_train[cat_cols].isna().sum()\n# get numerical column names without NaN\ncat_cols_no_nan = _[_==0].index.tolist()","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:30.785320Z","iopub.execute_input":"2024-04-08T15:12:30.785813Z","iopub.status.idle":"2024-04-08T15:12:32.804887Z","shell.execute_reply.started":"2024-04-08T15:12:30.785773Z","shell.execute_reply":"2024-04-08T15:12:32.803409Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get train data set stats\n_df_1 = df_train[cat_cols_no_nan].describe().T\n_df_1.drop(columns=['count'], inplace=True)\n_df_1.columns = [f\"train {_}\" for _ in _df_1.columns]\n# get test data set stats\n_df_2 = df_test[cat_cols_no_nan].describe().T\n_df_2.drop(columns=['count'], inplace=True)\n_df_2.columns = [f\"test {_}\" for _ in _df_2.columns]\npd.merge(_df_1, _df_2, left_index=True, right_index=True, how=\"outer\").sort_values(by=\"train unique\")","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:32.808339Z","iopub.execute_input":"2024-04-08T15:12:32.808867Z","iopub.status.idle":"2024-04-08T15:12:33.035937Z","shell.execute_reply.started":"2024-04-08T15:12:32.808824Z","shell.execute_reply":"2024-04-08T15:12:33.034659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# 7. Features (no NaN values)\n\n\n## 7.1. Numerical - low number of unique values\n\nStart with numerical features which have less than 10 and more than 1 unique observations. In addition, these features should not have a signle bucket with more that 95% of values in train data set.\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"# columns to select for analysis\nnum_cols_no_nan_low_unique = []\nfor _ in num_cols_no_nan:\n if 1 < df_train[_].nunique() < 10: \n num_cols_no_nan_low_unique.append(_)\n \nprint(f\"{len(num_cols_no_nan_low_unique)} numerical features were found with less than 10 unique values.\")\n\nnum_cols_no_nan_low_unique_2 = []\n# check whatever 1 value doesn't represent 95% of the data\nfor _ in num_cols_no_nan_low_unique:\n _df = df_train[_].value_counts(normalize=True)\n if _df.max() < 0.95:\n num_cols_no_nan_low_unique_2.append(_)\n \nprint(f\"{len(num_cols_no_nan_low_unique_2)}/{len(num_cols_no_nan_low_unique)} numerical features don't have > 95% of their data in one bucket.\")","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:33.038128Z","iopub.execute_input":"2024-04-08T15:12:33.038563Z","iopub.status.idle":"2024-04-08T15:12:33.861635Z","shell.execute_reply.started":"2024-04-08T15:12:33.038528Z","shell.execute_reply":"2024-04-08T15:12:33.860005Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_distributions(df, col_names, nrows = 5, rotation = 0):\n fig, ax = plt.subplots(int(np.ceil(len(col_names) / nrows)), nrows, figsize=(16, 8))\n ax = ax.flatten()\n \n # iterate over columns\n for i, col in enumerate(col_names):\n # count values\n _df = df[col].value_counts().reset_index()\n sns.barplot(data = _df, x = col, y = \"count\", ax = ax[i])\n if rotation:\n ax[i].set_xticklabels(ax[i].get_xticklabels(), rotation=rotation)\n \n plt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:33.871355Z","iopub.execute_input":"2024-04-08T15:12:33.871783Z","iopub.status.idle":"2024-04-08T15:12:33.882262Z","shell.execute_reply.started":"2024-04-08T15:12:33.871754Z","shell.execute_reply":"2024-04-08T15:12:33.881159Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_distributions(df_train, num_cols_no_nan_low_unique_2)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:33.884278Z","iopub.execute_input":"2024-04-08T15:12:33.884681Z","iopub.status.idle":"2024-04-08T15:12:36.726510Z","shell.execute_reply.started":"2024-04-08T15:12:33.884651Z","shell.execute_reply":"2024-04-08T15:12:36.724842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# mannualy calculating PSI is x3 faster than using function\ndata = list()\n\n# calculate weekly PSI scores\nfor col in num_cols_no_nan_low_unique_2:\n # all development data\n array_x = df_train.loc[:, col].values\n share_x = pd.DataFrame({\"bucket\": array_x}).value_counts().reset_index()\n share_x[\"proportion_x\"] = share_x[\"count\"] / share_x[\"count\"].sum()\n \n # iterate over weeks\n for week in df_train[\"WEEK_NUM\"].unique():\n array_y = df_train.loc[df_train['WEEK_NUM'] == week, col].values\n share_y = pd.DataFrame({\"bucket\": array_y}).value_counts().reset_index()\n share_y[\"proportion_y\"] = share_y[\"count\"] / share_y[\"count\"].sum()\n \n # calculate PSI\n df_shares = share_x.merge(share_y, on=\"bucket\", how = \"outer\").fillna(1e-6)\n df_shares['psi'] = (df_shares['proportion_x'] - df_shares['proportion_y']) * np.log(df_shares['proportion_x'] / df_shares['proportion_y'])\n \n data.append([col, week, df_shares['psi'].sum()])\n \ndf_plot = pd.DataFrame(data, columns=['Feature', \"WEEK_NUM\", \"PSI\"])","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:36.728569Z","iopub.execute_input":"2024-04-08T15:12:36.728943Z","iopub.status.idle":"2024-04-08T15:12:44.742257Z","shell.execute_reply.started":"2024-04-08T15:12:36.728914Z","shell.execute_reply":"2024-04-08T15:12:44.740149Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_lines(df, ncols = 2):\n fig, ax = plt.subplots(int(np.ceil(df['Feature'].nunique() / ncols)), ncols, figsize=(16, df['Feature'].nunique() / ncols * 2))\n ax = ax.flatten()\n \n # iterate over columns\n for i, col in enumerate(df['Feature'].unique()):\n # count values\n _df = df[df['Feature'] == col]\n ax[i].plot(_df['WEEK_NUM'], _df['PSI'], label=col, linewidth=2.5)\n ax[i].set_xlabel(\"Week num.\")\n ax[i].set_ylabel(\"PSI\")\n \n # add white background to legend\n legend = ax[i].legend(frameon=1)\n frame = legend.get_frame()\n frame.set_facecolor('w')\n \n \n plt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:44.744443Z","iopub.execute_input":"2024-04-08T15:12:44.744901Z","iopub.status.idle":"2024-04-08T15:12:44.754641Z","shell.execute_reply.started":"2024-04-08T15:12:44.744867Z","shell.execute_reply":"2024-04-08T15:12:44.753397Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_lines(df_plot)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:44.756057Z","iopub.execute_input":"2024-04-08T15:12:44.756461Z","iopub.status.idle":"2024-04-08T15:12:47.008504Z","shell.execute_reply.started":"2024-04-08T15:12:44.756431Z","shell.execute_reply":"2024-04-08T15:12:47.007506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# find TOP 10 most stable, i.e. lowest weekly PSI variance, features \n_df = df_plot.groupby(\"Feature\").agg({\"PSI\": [\"mean\", \"median\", \"max\", \"min\", \"std\"]}).round(3).sort_values(by=(\"PSI\", \"std\"))\n_df","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:47.009992Z","iopub.execute_input":"2024-04-08T15:12:47.010641Z","iopub.status.idle":"2024-04-08T15:12:47.051371Z","shell.execute_reply.started":"2024-04-08T15:12:47.010608Z","shell.execute_reply":"2024-04-08T15:12:47.050055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 7.2. Numerical - more than 10 unique values\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"# columns to select for analysis\nnum_cols_no_nan_high_unique = []\nfor _ in num_cols_no_nan:\n if 1 < df_train[_].nunique() >= 10: \n num_cols_no_nan_high_unique.append(_)\n \n# exclude month decision as there is no point in q cutting it\nnum_cols_no_nan_high_unique.remove('month_decision')\n \nprint(f\"{len(num_cols_no_nan_high_unique)} numerical features were found with more than 10 unique values.\")\n\nnum_cols_no_nan_high_unique_2 = []\n# check whatever 1 value doesn't represent 95% of the data\nfor _ in num_cols_no_nan_high_unique:\n _df = df_train[_].value_counts(normalize=True)\n if _df.max() < 0.95:\n num_cols_no_nan_high_unique_2.append(_)\n \nprint(f\"{len(num_cols_no_nan_high_unique_2)}/{len(num_cols_no_nan_high_unique)} numerical features don't have > 95% of their data in represented by one numerical value.\")","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:47.052858Z","iopub.execute_input":"2024-04-08T15:12:47.053228Z","iopub.status.idle":"2024-04-08T15:12:48.277857Z","shell.execute_reply.started":"2024-04-08T15:12:47.053197Z","shell.execute_reply":"2024-04-08T15:12:48.276727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# number of splits\nN = 10\n\ndf_train_qcuts = pd.DataFrame()\ndf_test_qcuts = pd.DataFrame()\n\n# create quantiles based on whole train data set\nfor _ in num_cols_no_nan_high_unique_2:\n ser, bins = pd.qcut(df_train[_].round(1), N, duplicates='drop', retbins=True)\n # exclude features with low number of unique quantile cuts\n if len(set(bins)) >= 10:\n df_train_qcuts[f\"{_}_qcut\"] = ser\n df_test_qcuts[f\"{_}_qcut\"] = pd.cut(df_test[_], bins=bins, include_lowest=True)\n \n# col names for plotting\nCOLS_10_bins = df_train_qcuts.columns","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:48.279334Z","iopub.execute_input":"2024-04-08T15:12:48.279697Z","iopub.status.idle":"2024-04-08T15:12:50.045716Z","shell.execute_reply.started":"2024-04-08T15:12:48.279666Z","shell.execute_reply":"2024-04-08T15:12:50.044535Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# join tables\ndf_train = pd.concat([df_train, df_train_qcuts], axis=1)\ndf_test = pd.concat([df_test, df_test_qcuts], axis=1)\n\ndel df_train_qcuts, df_test_qcuts\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:50.047288Z","iopub.execute_input":"2024-04-08T15:12:50.048666Z","iopub.status.idle":"2024-04-08T15:12:53.427406Z","shell.execute_reply.started":"2024-04-08T15:12:50.048613Z","shell.execute_reply":"2024-04-08T15:12:53.425658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_distributions(df_train, COLS_10_bins.values.tolist(), 3, 45)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:53.429119Z","iopub.execute_input":"2024-04-08T15:12:53.429915Z","iopub.status.idle":"2024-04-08T15:12:55.338507Z","shell.execute_reply.started":"2024-04-08T15:12:53.429879Z","shell.execute_reply":"2024-04-08T15:12:55.337401Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# mannualy calculating PSI is x3 faster than using function\ndata = list()\n\n# calculate weekly PSI scores\nfor col in COLS_10_bins:\n # all development data\n array_x = df_train.loc[:, col].values\n share_x = pd.DataFrame({\"bucket\": array_x}).value_counts().reset_index()\n share_x[\"proportion_x\"] = share_x[\"count\"] / share_x[\"count\"].sum()\n \n # iterate over weeks\n for week in df_train[\"WEEK_NUM\"].unique():\n array_y = df_train.loc[df_train['WEEK_NUM'] == week, col].values\n share_y = pd.DataFrame({\"bucket\": array_y}).value_counts().reset_index()\n share_y[\"proportion_y\"] = share_y[\"count\"] / share_y[\"count\"].sum()\n \n # calculate PSI\n df_shares = share_x.merge(share_y, on=\"bucket\", how = \"outer\")\n df_shares[['proportion_x', 'proportion_y']] = df_shares[['proportion_x', 'proportion_y']].fillna(1e-6)\n df_shares['psi'] = (df_shares['proportion_x'] - df_shares['proportion_y']) * np.log(df_shares['proportion_x'] / df_shares['proportion_y'])\n \n data.append([col, week, df_shares['psi'].sum()])\n \ndf_plot = pd.DataFrame(data, columns=['Feature', \"WEEK_NUM\", \"PSI\"])","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:55.340036Z","iopub.execute_input":"2024-04-08T15:12:55.341329Z","iopub.status.idle":"2024-04-08T15:12:59.650714Z","shell.execute_reply.started":"2024-04-08T15:12:55.341262Z","shell.execute_reply":"2024-04-08T15:12:59.649264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_lines(df_plot)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:12:59.652363Z","iopub.execute_input":"2024-04-08T15:12:59.652703Z","iopub.status.idle":"2024-04-08T15:13:02.483140Z","shell.execute_reply.started":"2024-04-08T15:12:59.652675Z","shell.execute_reply":"2024-04-08T15:13:02.479714Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get PSI variance stats \n_df = df_plot.groupby(\"Feature\").agg({\"PSI\": [\"mean\", \"median\", \"max\", \"min\", \"std\"]}).round(3).sort_values(by=(\"PSI\", \"std\"))\n_df","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:13:02.486839Z","iopub.execute_input":"2024-04-08T15:13:02.487678Z","iopub.status.idle":"2024-04-08T15:13:02.565274Z","shell.execute_reply.started":"2024-04-08T15:13:02.487631Z","shell.execute_reply":"2024-04-08T15:13:02.562504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 7.3. Categorical - low number of unique values\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"# columns to select for analysis\ncat_cols_no_nan_low_unique = []\nfor _ in cat_cols_no_nan:\n if 1 < df_train[_].nunique() < 10: \n cat_cols_no_nan_low_unique.append(_)\n \nprint(f\"{len(cat_cols_no_nan_low_unique)} categorical features were found with less than 10 unique values.\")\n\ncat_cols_no_nan_low_unique_2 = []\n# check whatever 1 value doesn't represent 95% of the data\nfor _ in cat_cols_no_nan_low_unique:\n _df = df_train[_].value_counts(normalize=True)\n if _df.max() < 0.95:\n cat_cols_no_nan_low_unique_2.append(_)\n \nprint(f\"{len(cat_cols_no_nan_low_unique_2)}/{len(cat_cols_no_nan_low_unique)} categorical features don't have > 95% of their data in one bucket.\")","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:13:02.569111Z","iopub.execute_input":"2024-04-08T15:13:02.569892Z","iopub.status.idle":"2024-04-08T15:13:02.888861Z","shell.execute_reply.started":"2024-04-08T15:13:02.569834Z","shell.execute_reply":"2024-04-08T15:13:02.887559Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_distributions(df_train, cat_cols_no_nan_low_unique_2, 2)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:13:02.890097Z","iopub.execute_input":"2024-04-08T15:13:02.890690Z","iopub.status.idle":"2024-04-08T15:13:05.149163Z","shell.execute_reply.started":"2024-04-08T15:13:02.890655Z","shell.execute_reply":"2024-04-08T15:13:05.146808Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# mannualy calculating PSI is x3 faster than using function\ndata = list()\n\n# calculate weekly PSI scores\nfor col in cat_cols_no_nan_low_unique_2:\n # all development data\n array_x = df_train.loc[:, col].values\n if col == 'max_incometype_1044T':\n array_x = array_x.astype(str)\n array_x[array_x == \"HANDICAPPED_2\"] = \"HANDICAPPED\"\n array_x[array_x == \"HANDICAPPED_3\"] = \"HANDICAPPED\"\n \n share_x = pd.DataFrame({\"bucket\": array_x}).value_counts().reset_index()\n share_x[\"proportion_x\"] = share_x[\"count\"] / share_x[\"count\"].sum()\n \n # iterate over weeks\n for week in df_train[\"WEEK_NUM\"].unique():\n array_y = df_train.loc[df_train['WEEK_NUM'] == week, col].values\n if col == 'max_incometype_1044T':\n array_y = array_y.astype(str)\n array_y[array_y == \"HANDICAPPED_2\"] = \"HANDICAPPED\"\n array_y[array_y == \"HANDICAPPED_3\"] = \"HANDICAPPED\"\n share_y = pd.DataFrame({\"bucket\": array_y}).value_counts().reset_index()\n share_y[\"proportion_y\"] = share_y[\"count\"] / share_y[\"count\"].sum()\n \n # calculate PSI\n df_shares = share_x.merge(share_y, on=\"bucket\", how = \"outer\")\n df_shares[[\"count_x\", \"count_y\"]] = df_shares[[\"count_x\", \"count_y\"]].fillna(0)\n df_shares[[\"proportion_x\", \"proportion_y\"]] = df_shares[[\"proportion_x\", \"proportion_y\"]].fillna(1e-6)\n df_shares['psi'] = (df_shares['proportion_x'] - df_shares['proportion_y']) * np.log(df_shares['proportion_x'] / df_shares['proportion_y'])\n \n data.append([col, week, df_shares['psi'].sum()])\n \ndf_plot = pd.DataFrame(data, columns=['Feature', \"WEEK_NUM\", \"PSI\"])","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:13:05.154817Z","iopub.execute_input":"2024-04-08T15:13:05.158088Z","iopub.status.idle":"2024-04-08T15:13:17.491171Z","shell.execute_reply.started":"2024-04-08T15:13:05.157966Z","shell.execute_reply":"2024-04-08T15:13:17.488789Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_lines(df_plot)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:13:17.493808Z","iopub.execute_input":"2024-04-08T15:13:17.494901Z","iopub.status.idle":"2024-04-08T15:13:20.487985Z","shell.execute_reply.started":"2024-04-08T15:13:17.494817Z","shell.execute_reply":"2024-04-08T15:13:20.483870Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get PSI variance stats \n_df = df_plot.groupby(\"Feature\").agg({\"PSI\": [\"mean\", \"median\", \"max\", \"min\", \"std\"]}).round(3).sort_values(by=(\"PSI\", \"std\"))\n_df","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:13:20.493850Z","iopub.execute_input":"2024-04-08T15:13:20.494949Z","iopub.status.idle":"2024-04-08T15:13:20.566641Z","shell.execute_reply.started":"2024-04-08T15:13:20.494852Z","shell.execute_reply":"2024-04-08T15:13:20.561825Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 7.4. Categorical - more than 10 unique values\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"# columns to select for analysis\ncat_cols_no_nan_high_unique = []\nfor _ in cat_cols_no_nan:\n if 1 < df_train[_].nunique() >= 10: \n cat_cols_no_nan_high_unique.append(_)\n \nprint(f\"{len(cat_cols_no_nan_high_unique)} categorical features were found with more than 10 unique values.\")\n\ncat_cols_no_nan_high_unique_2 = []\n# check whatever 1 value doesn't represent 95% of the data\nfor _ in cat_cols_no_nan_high_unique:\n _df = df_train[_].value_counts(normalize=True)\n if _df.max() < 0.95:\n cat_cols_no_nan_high_unique_2.append(_)\n \nprint(f\"{len(cat_cols_no_nan_high_unique_2)}/{len(cat_cols_no_nan_high_unique)} categorical features don't have > 95% of their data in one bucket.\")","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:13:20.573149Z","iopub.execute_input":"2024-04-08T15:13:20.574413Z","iopub.status.idle":"2024-04-08T15:13:20.948517Z","shell.execute_reply.started":"2024-04-08T15:13:20.574290Z","shell.execute_reply":"2024-04-08T15:13:20.944676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_distributions(df_train, cat_cols_no_nan_high_unique_2, 2)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:13:20.955700Z","iopub.execute_input":"2024-04-08T15:13:20.957864Z","iopub.status.idle":"2024-04-08T15:13:32.624462Z","shell.execute_reply.started":"2024-04-08T15:13:20.957695Z","shell.execute_reply":"2024-04-08T15:13:32.619925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# mannualy calculating PSI is x3 faster than using function\ndata = list()\n\n# calculate weekly PSI scores\nfor col in cat_cols_no_nan_high_unique_2:\n # all development data\n array_x = df_train.loc[:, col].values\n \n share_x = pd.DataFrame({\"bucket\": array_x}).value_counts().reset_index()\n share_x[\"proportion_x\"] = share_x[\"count\"] / share_x[\"count\"].sum()\n \n # iterate over weeks\n for week in df_train[\"WEEK_NUM\"].unique():\n array_y = df_train.loc[df_train['WEEK_NUM'] == week, col].values\n share_y = pd.DataFrame({\"bucket\": array_y}).value_counts().reset_index()\n share_y[\"proportion_y\"] = share_y[\"count\"] / share_y[\"count\"].sum()\n \n # calculate PSI\n df_shares = share_x.merge(share_y, on=\"bucket\", how = \"outer\")\n df_shares[[\"count_x\", \"count_y\"]] = df_shares[[\"count_x\", \"count_y\"]].fillna(0)\n df_shares[[\"proportion_x\", \"proportion_y\"]] = df_shares[[\"proportion_x\", \"proportion_y\"]].fillna(1e-6)\n df_shares['psi'] = (df_shares['proportion_x'] - df_shares['proportion_y']) * np.log(df_shares['proportion_x'] / df_shares['proportion_y'])\n \n data.append([col, week, df_shares['psi'].sum()])\n \ndf_plot = pd.DataFrame(data, columns=['Feature', \"WEEK_NUM\", \"PSI\"])","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:13:32.630003Z","iopub.execute_input":"2024-04-08T15:13:32.633160Z","iopub.status.idle":"2024-04-08T15:13:39.463407Z","shell.execute_reply.started":"2024-04-08T15:13:32.632995Z","shell.execute_reply":"2024-04-08T15:13:39.461673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# to many unique categories, i.e. feature is unstable over time\nplot_lines(df_plot)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:13:39.464881Z","iopub.execute_input":"2024-04-08T15:13:39.465283Z","iopub.status.idle":"2024-04-08T15:13:41.064615Z","shell.execute_reply.started":"2024-04-08T15:13:39.465245Z","shell.execute_reply":"2024-04-08T15:13:41.063395Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# 8. Modeling\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"%%time\n# qcut columns\nQCUT_COLS = ['annuity_780A_qcut', 'credamount_770A_qcut', 'disbursedcredamount_1113A_qcut', 'max_mainoccupationinc_384A_qcut', 'max_birth_259D_qcut']\n\n# use label encoded for qcut features\n# for more details refer here\n# https://stackoverflow.com/questions/54507269/valueerror-circular-reference-detected-in-lightgbm\nfor col_name in QCUT_COLS:\n le = LabelEncoder()\n le.fit(df_train[col_name].values)\n # re-transform features\n df_train.loc[df_train[col_name].notnull(), f\"{col_name}_new\"] = le.transform(df_train.loc[df_train[col_name].notnull(), col_name])\n df_test.loc[df_test[col_name].notnull(), f\"{col_name}_new\"] = le.transform(df_test.loc[df_test[col_name].notnull(), col_name])","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:13:41.066037Z","iopub.execute_input":"2024-04-08T15:13:41.066420Z","iopub.status.idle":"2024-04-08T15:14:35.550257Z","shell.execute_reply.started":"2024-04-08T15:13:41.066388Z","shell.execute_reply":"2024-04-08T15:14:35.548564Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"TRAIN_COLS = []\n# low weekly PSI numerical features with no NaN and less than 10 unique values\nTRAIN_COLS += ['clientscnt_533L', 'numactivecredschannel_414L', 'sellerplacecnt_915L', 'numactiverelcontr_750L', 'numactivecreds_622L']\n# remaining numerical features with no NaN and less than 10 unique values (higher PSI scores)\nTRAIN_COLS += ['max_persontype_1072L', 'max_persontype_792L', 'max_personindex_1023L', 'weekday_decision', 'max_num_group1_9']\n# numerical features without NaN values and large number of unique values\nTRAIN_COLS += [ f\"{_}_new\" for _ in QCUT_COLS]\n# categorical features with low number of unique values\nTRAIN_COLS += ['max_language1_981M', 'max_incometype_1044T', 'max_role_1084L', 'max_sex_738L']\n# categorical features with high number of unique values\nTRAIN_COLS += ['lastapprcommoditycat_1041M', 'lastcancelreason_561M', 'lastrejectcommoditycat_161M', 'lastrejectreason_759M', 'lastrejectreasonclient_4145040M']","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:14:35.552265Z","iopub.execute_input":"2024-04-08T15:14:35.552656Z","iopub.status.idle":"2024-04-08T15:14:35.561584Z","shell.execute_reply.started":"2024-04-08T15:14:35.552625Z","shell.execute_reply":"2024-04-08T15:14:35.560249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ncase_ids_train, case_ids_test = train_test_split(df_train['case_id'], train_size=0.8, random_state=SEED)\n\nX_train = df_train.loc[df_train['case_id'].isin(case_ids_train), TRAIN_COLS + [\"WEEK_NUM\"]]\nX_test = df_train.loc[df_train['case_id'].isin(case_ids_test), TRAIN_COLS + [\"WEEK_NUM\"] + ['target']]\n\ny_train = df_train.loc[df_train['case_id'].isin(case_ids_train), 'target']\ny_test = df_train.loc[df_train['case_id'].isin(case_ids_test), 'target']","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:14:35.563432Z","iopub.execute_input":"2024-04-08T15:14:35.563880Z","iopub.status.idle":"2024-04-08T15:14:36.314378Z","shell.execute_reply.started":"2024-04-08T15:14:35.563838Z","shell.execute_reply":"2024-04-08T15:14:36.313119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nNO_MODELS = 5\n\n# use week numbers for group spliting\nweeks = X_train['WEEK_NUM']\n\ncv = StratifiedGroupKFold(n_splits = NO_MODELS, shuffle=False)\n\nensemble_models = {}\n\n# COMPUTE CV SCORE WITH 5 GROUP K FOLD\nfor i, (train_index, test_index) in enumerate(cv.split(X_train, y_train, groups=weeks)):\n print('#'*25)\n print('### Fold',i+1)\n print('#'*25)\n \n # parameters were taken from this notebook:\n # https://www.kaggle.com/code/greysky/home-credit-baseline\n lgb_params = {\n \"boosting_type\": \"gbdt\",\n 'objective' : 'binary',\n 'metric' : 'auc',\n \"learning_rate\": 0.05,\n \"n_estimators\": 1500,\n \"colsample_bytree\": 0.8, \n \"colsample_bynode\": 0.8,\n \"verbose\": -1,\n \"reg_alpha\": 0.1,\n \"reg_lambda\": 10,\n \"extra_trees\":True,\n 'num_leaves':64,\n \"random_state\": SEED + i\n }\n \n # TRAIN DATA\n train_x = X_train.iloc[train_index][TRAIN_COLS]\n train_y = y_train.iloc[train_index]\n\n # VALID DATA\n valid_x = X_train.iloc[test_index][TRAIN_COLS]\n valid_y = y_train.iloc[test_index]\n \n # CONVERT TO DataSet\n train_dataset = lgb.Dataset(train_x, label=train_y)\n valid_dataset = lgb.Dataset(valid_x, label=valid_y, reference=train_dataset)\n \n # TRAIN MODEL\n model = lgb.train(\n lgb_params,\n train_dataset,\n valid_sets = valid_dataset,\n callbacks = [lgb.log_evaluation(100), lgb.early_stopping(50)]\n )\n \n # VALIDATE MODELS\n preds = model.predict(valid_x, num_iteration=model.best_iteration)\n \n # SCORE ON VALID SAMPLE\n auc = roc_auc_score(valid_y, preds)\n gini = 2 * auc - 1\n print(f\"\\nValidation AUC: {auc:.3f}\")\n print(f\"Validation GINI: {gini:.3f}\")\n\n # SAVE MODEL\n ensemble_models[f'model_{i}'] = model\n \n print()","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:14:36.316157Z","iopub.execute_input":"2024-04-08T15:14:36.316529Z","iopub.status.idle":"2024-04-08T15:22:09.161733Z","shell.execute_reply.started":"2024-04-08T15:14:36.316500Z","shell.execute_reply":"2024-04-08T15:22:09.159665Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf_preds = pd.DataFrame()\nfor i in range(NO_MODELS):\n model = ensemble_models[f'model_{i}'] \n df_preds[f'preds_{i}'] = model.predict(X_test[TRAIN_COLS], num_iteration=model.best_iteration)\n \n# add average predictions\nX_test['pred'] = df_preds.mean(axis = 1).values","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:22:09.164067Z","iopub.execute_input":"2024-04-08T15:22:09.164532Z","iopub.status.idle":"2024-04-08T15:22:46.745515Z","shell.execute_reply.started":"2024-04-08T15:22:09.164489Z","shell.execute_reply":"2024-04-08T15:22:46.742350Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# 9. Validation\n\n\n## 9.1 Classical\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"def validate_model(data, output = \"Chart\"):\n \"\"\"\n Input:\n preds, array like object\n target, array like object\n chart, boolean, if true return matplolib else table of summary statistics\n Output:\n pandas DataFrame or matplolib object\n \n Needs equal length prediction and taget column vectors\n \"\"\"\n # create DataFrame for plotting and summary stats\n _df = data.copy()\n _roc = data.loc[:, [\"WEEK_NUM\", \"target\", \"pred\"]]\\\n .sort_values(\"WEEK_NUM\").groupby(\"WEEK_NUM\")[[\"target\", \"pred\"]]\\\n .apply(lambda x: roc_auc_score(x[\"target\"], x[\"pred\"]))\n _df['ROC_AUC'] = _df['WEEK_NUM'].map(_roc)\n _df['GINI'] = 2 * _df['ROC_AUC'] - 1\n \n # stability score\n x = _df[['WEEK_NUM', 'GINI']].drop_duplicates()['WEEK_NUM']\n y = _df[['WEEK_NUM', 'GINI']].drop_duplicates()['GINI']\n a, b = np.polyfit(x, y, 1)\n y_hat = a*x + b\n residuals = y - y_hat\n res_std = np.std(residuals)\n avg_gini = _df[['WEEK_NUM', 'GINI']].drop_duplicates()['GINI'].mean()\n w_fallingrate = 88.0\n w_resstd = -0.5 \n # calculate competition score\n score = avg_gini + w_fallingrate * min(0, a) + w_resstd * res_std\n \n if output == \"Chart\":\n fig, ax = plt.subplots(2, 2, figsize=(16, 10))\n \n # temporal changes of gini scores\n _df_plot = _df.groupby('WEEK_NUM')['GINI'].mean().reset_index()\n ax[0, 0].plot(_df_plot['WEEK_NUM'], _df_plot['GINI'], \"o-\", label = 'Weekly GINI scores')\n # add linear fit\n ax[0, 0].plot(_df_plot['WEEK_NUM'], _df_plot['WEEK_NUM'] * a + b, \"r--\", label = f'linear fit, k={a:.3f}')\n ax[0, 0].set_ylabel(\"Gini score\")\n ax[0, 0].set_xlabel(\"Week number\")\n \n # GINI score boxplot\n sns.boxplot(data=_df_plot, x=\"GINI\", ax = ax[0, 1], color = COLORS[2])\n ax[0, 1].set_xlabel(\"Gini score\")\n ax[0, 1].set_title(f\"Q1: {np.quantile(_df_plot['GINI'].values, 0.25):.3f}, Q3: {np.quantile(_df_plot['GINI'].values, 0.75):.3f}\")\n \n # GINI score histogram\n ax[1, 1].hist(_df_plot['GINI'], bins=20, label=\"Weekly GINI scores\", color = COLORS[2])\n ax[1, 1].set_ylabel(\"Count\")\n ax[1, 1].set_xlabel(\"Gini score\")\n ax[1, 1].set_title(f\"Mean: {_df_plot['GINI'].mean():.3f}, STD: {_df_plot['GINI'].std():.3f}\")\n \n # GINI vs number of observations\n _df_plot = _df.groupby('WEEK_NUM').agg({\"GINI\": ['mean', \"count\"]})\n sns.regplot(x=_df_plot.iloc[:, 1], y=_df_plot.iloc[:, 0], ax=ax[1, 0], label=f'Corr. {_df_plot.corr().values[0, 1]:.3f}', color = COLORS[1])\n ax[1, 0].set_ylabel(\"Weekly Gini score\")\n ax[1, 0].set_xlabel(\"Number of applications\")\n \n # add white backgrounds to legends\n for i in range(2):\n for ii in range(2):\n legend = ax[i, ii].legend(frameon=1)\n frame = legend.get_frame()\n frame.set_facecolor('w')\n \n # return DataFrame with summary stats\n elif output == \"Stats\":\n _df_stats = pd.DataFrame({\"dataset\": [\"test\"]})\n _df_stats['Avg. weekly ROC AUC'] = _df[['WEEK_NUM', 'ROC_AUC']].drop_duplicates()['ROC_AUC'].mean().round(3)\n _df_stats['Std. of weekly ROC AUC'] = _df[['WEEK_NUM', 'ROC_AUC']].drop_duplicates()['ROC_AUC'].std().round(3)\n _df_stats['Avg. weekly GINI'] = _df[['WEEK_NUM', 'GINI']].drop_duplicates()['GINI'].mean().round(3)\n _df_stats['Std. of weekly GINI'] = _df[['WEEK_NUM', 'GINI']].drop_duplicates()['GINI'].std().round(3)\n _df_stats['Stability score'] = round(score, 3)\n \n return _df_stats\n else:\n return _df","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:22:46.747775Z","iopub.execute_input":"2024-04-08T15:22:46.748337Z","iopub.status.idle":"2024-04-08T15:22:46.787808Z","shell.execute_reply.started":"2024-04-08T15:22:46.748271Z","shell.execute_reply":"2024-04-08T15:22:46.785475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# make predictions on test data set\nvalidate_model(X_test)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:22:46.790379Z","iopub.execute_input":"2024-04-08T15:22:46.790974Z","iopub.status.idle":"2024-04-08T15:22:48.709099Z","shell.execute_reply.started":"2024-04-08T15:22:46.790916Z","shell.execute_reply":"2024-04-08T15:22:48.707487Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"validate_model(X_test, \"Stats\")","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:22:48.710991Z","iopub.execute_input":"2024-04-08T15:22:48.711394Z","iopub.status.idle":"2024-04-08T15:22:49.162443Z","shell.execute_reply.started":"2024-04-08T15:22:48.711361Z","shell.execute_reply":"2024-04-08T15:22:49.160896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 9.2 PSI validation\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"N = 5\nser = pd.qcut(X_test['pred'], N, duplicates='drop')\nle = LabelEncoder()\nX_test['pred_bins'] = le.fit_transform(ser.values)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:22:49.164174Z","iopub.execute_input":"2024-04-08T15:22:49.165506Z","iopub.status.idle":"2024-04-08T15:22:50.269729Z","shell.execute_reply.started":"2024-04-08T15:22:49.165468Z","shell.execute_reply":"2024-04-08T15:22:50.268217Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# mannualy calculating PSI is x3 faster than using function\ndata = list()\n\n# all development data\narray_x = X_test.loc[:, 'pred_bins'].values\nshare_x = pd.DataFrame({\"bucket\": array_x}).value_counts().reset_index()\nshare_x[\"proportion_x\"] = share_x[\"count\"] / share_x[\"count\"].sum()\n\n# iterate over weeks\nfor week in X_test[\"WEEK_NUM\"].unique():\n array_y = X_test.loc[X_test['WEEK_NUM'] == week, 'pred_bins'].values\n share_y = pd.DataFrame({\"bucket\": array_y}).value_counts().reset_index()\n share_y[\"proportion_y\"] = share_y[\"count\"] / share_y[\"count\"].sum()\n\n # calculate PSI\n df_shares = share_x.merge(share_y, on=\"bucket\", how = \"outer\")\n df_shares[['proportion_x', 'proportion_y']] = df_shares[['proportion_x', 'proportion_y']].fillna(1e-6)\n df_shares['psi'] = (df_shares['proportion_x'] - df_shares['proportion_y']) * np.log(df_shares['proportion_x'] / df_shares['proportion_y'])\n\n data.append([col, week, df_shares['psi'].sum()])\n \ndf_plot = pd.DataFrame(data, columns=['Feature', \"WEEK_NUM\", \"PSI\"])","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:22:50.274542Z","iopub.execute_input":"2024-04-08T15:22:50.274941Z","iopub.status.idle":"2024-04-08T15:22:51.052353Z","shell.execute_reply.started":"2024-04-08T15:22:50.274912Z","shell.execute_reply":"2024-04-08T15:22:51.049420Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(1, 2, figsize=(16, 5))\n\n# quick data for plotting\n_df = X_test.groupby(['WEEK_NUM'])['pred_bins'].value_counts(True).reset_index()\n_df['proportion'] = _df['proportion'] * 100\n\ny_bot = np.zeros(_df['WEEK_NUM'].nunique())\n\nfor _bin in range(_df[['pred_bins']].nunique()[0]):\n _df_b = _df[_df['pred_bins'] == _bin]\n ax[0].bar(_df_b['WEEK_NUM'], _df_b['proportion'], label=_bin, bottom=y_bot)\n y_bot += _df_b['proportion'].values\n\nax[1].plot(df_plot['WEEK_NUM'], df_plot['PSI'], label=\"pred buckets\", linewidth=2.5)\nax[0].set_ylabel(\"Share, %\")\nax[1].set_ylabel(\"PSI\")\n\n# add white background to legend\nfor _ in range(2):\n ax[_].set_xlabel(\"Week num.\")\n legend = ax[_].legend(frameon=1)\n frame = legend.get_frame()\n frame.set_facecolor('w')\n\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:22:51.054226Z","iopub.execute_input":"2024-04-08T15:22:51.054906Z","iopub.status.idle":"2024-04-08T15:22:53.928125Z","shell.execute_reply.started":"2024-04-08T15:22:51.054869Z","shell.execute_reply":"2024-04-08T15:22:53.926061Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 9.3 Risk differentiation\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(1, 2, figsize=(16, 5))\n\n# quick data for plotting\n_df = X_test.groupby(['WEEK_NUM', 'pred_bins'])['target'].mean(True).reset_index()\n_df['target'] = _df['target'] * 100\n\nsns.scatterplot(data=_df, x=\"WEEK_NUM\", y=\"target\", hue=\"pred_bins\", ax = ax[0], palette=\"crest\")\nsns.boxplot(data=_df, y=\"target\", x=\"pred_bins\", ax = ax[1], palette=\"crest\")\n\nax[0].set_ylabel(\"ODF, %\")\nax[0].set_xlabel(\"Week num.\")\n\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:22:53.929948Z","iopub.execute_input":"2024-04-08T15:22:53.930376Z","iopub.status.idle":"2024-04-08T15:22:54.929802Z","shell.execute_reply.started":"2024-04-08T15:22:53.930334Z","shell.execute_reply":"2024-04-08T15:22:54.928391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"_df = validate_model(X_test, \"Raw\")\n_df = _df.groupby(['WEEK_NUM'])['GINI'].mean().reset_index()\n_df = pd.merge(df_plot, _df, left_on='WEEK_NUM', right_on='WEEK_NUM', how=\"inner\")\n\nfig, ax = plt.subplots(figsize=(16, 5))\n\nsns.scatterplot(data=_df, x=\"PSI\", y=\"GINI\", ax = ax)\n\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:22:54.931527Z","iopub.execute_input":"2024-04-08T15:22:54.931891Z","iopub.status.idle":"2024-04-08T15:22:55.790733Z","shell.execute_reply.started":"2024-04-08T15:22:54.931861Z","shell.execute_reply":"2024-04-08T15:22:55.789426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# 10. Submit predictions\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"%%time\ndf_preds = pd.DataFrame()\nfor i in range(NO_MODELS):\n model = ensemble_models[f'model_{i}'] \n df_preds[f'preds_{i}'] = model.predict(df_test[TRAIN_COLS], num_iteration=model.best_iteration)","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:22:55.792636Z","iopub.execute_input":"2024-04-08T15:22:55.793036Z","iopub.status.idle":"2024-04-08T15:22:55.873882Z","shell.execute_reply.started":"2024-04-08T15:22:55.793004Z","shell.execute_reply":"2024-04-08T15:22:55.872246Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sub = pd.read_csv(ROOT / \"sample_submission.csv\")\ndf_sub = df_sub.set_index(\"case_id\")\ndf_sub[\"score\"] = df_preds.mean(axis = 1).values","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:22:55.876068Z","iopub.execute_input":"2024-04-08T15:22:55.876989Z","iopub.status.idle":"2024-04-08T15:22:55.908971Z","shell.execute_reply.started":"2024-04-08T15:22:55.876943Z","shell.execute_reply":"2024-04-08T15:22:55.907409Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sub.to_csv(\"submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-04-08T15:22:55.910859Z","iopub.execute_input":"2024-04-08T15:22:55.911370Z","iopub.status.idle":"2024-04-08T15:22:55.923721Z","shell.execute_reply.started":"2024-04-08T15:22:55.911322Z","shell.execute_reply":"2024-04-08T15:22:55.921995Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# 11. Changelog\n\n* **v1** created a basic structure for the notebook (added table of contents). Loaded data and created new features using code from @daviddirethucus's notebook https://www.kaggle.com/code/daviddirethucus/home-credit-risk-lightgbm.\n* **v2** added custom PSI calculation function, reduced memory for train and test pandas DataFrames. Described numerical features without missing values and ploted an example of distribution for `weekday_decision` feature for the first 40 weeks. PSI was calcualted between first 2 weeks distributions and distribution of first week agains whole train data set.\n* **v3** added visualization for theoretical distributions and calculated PSI scores. Selected 10 numerical features with less than 10 unique values. Calculated weekly PSI scores against whole development data and provided a quick visualization.\n* **v4** trained baseline logistic regression models on numerical features that are stable across time, i.e. low median and standard variation PSI weekly scores. The average ROC AUC score on 5-fold cross-validation samples was around `~0.55`.\n* **v5** save submission DataFrame to .csv, which I forgot to do in last notebook version.\n* **v6** change baseline model from log. regression to LGB as it easier handles missing values(there are features which don't have missing values in train data set, but have missing values in train). When model was trained on all numerical features with low bumber of unique values, the average ROC AUC score on 5-fold cross-validation samples increased to around `~0.58`.\n* **v7** selected numerical features without missing values and put them into 10 equal buckets by percentile. The average ROC AUC score on 5-fold cross-validation samples increased to around `0.62` compared to previous notebook versions.\n* **v8** added features for training from notebook versions 5 and 7. The average ROC AUC score on 5-fold cross-validation samples increased to around `~0.64` compared to previous notebook versions.\n* **v9** trained model on all features and added additional PSI validation curve for the test set.\n* **v10** added additional validation charts for predictions on test data set.\n* **v11** minor text fixes.\n* **v12** added new sub-section where PSI score was calculated on categorical features with low number of unique values.\n* **v13** tried calculating PSI score for catergorical features with high number of unique values, but could not. Re-trained model only of those features.\n* **v14** trained model on all features which were investigated using PSI in section [7. Features (no NaN values)](#7).\n\n
\n
\nBack to top\n
","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}