{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Introduction & motivations\n\n* Background & context\n> Historically, underwriters within insurance companies have worked closely with pricing actuaries to assess and manage the risks associated with providing insurance coverage to their clients, as well as to determine a suitable price/rate (known as the pure risk premium) at which an agreed level of risk coverage can be provided - as per the insurance company's risk management policy. In the past, underwriters would have needed to manually review a lot of information before they could decide what price to charge for a policy - for instance, they would need to determine:\n>\n> * what type of peril/s the insurance company agrees to insure,\n> * what type of policy coverage can be offered, given both the insurance company's policies and the client's requirements, and\n> * what level of risk is posed against the insured, and whether there is an appropriate rate that can be charged.\n>\n>>\n> Nowadays, the underwriting process is becoming more and more automated, supported \"by a combination of machine and deep learning models built within companies' technology stacks\" ([McKinsey, Insurance 2030 — The impact of AI on the future of insurance](https://www.mckinsey.com/industries/financial-services/our-insights/insurance-2030-the-impact-of-ai-on-the-future-of-insurance)), which enables insurers to make rapid decisions regarding underwriting/pricing, as well as to extensively provide tailored coverage that is specific to each client's risk profile.\n\n* What is this project about?\n> The aim of the project is to provide a demonstration of how un/supervised ML techniques can be used to predict risk ratings for life insurance applicants, which could then in turn support underwriters in their decision-making process on how prospective new business should be valued. As a case study, we will use the dataset featured in the **Prudential Life Insurance Assessment** competition previously hosted on Kaggle, and showcase how ML classifiers can be used to quantitatively assess risk.\n>\n> More information on the dataset used can be found at the following page: [Prudential Life Insurance Assessment](https://www.kaggle.com/competitions/prudential-life-insurance-assessment/).\n\n* What will be discussed/shown in the project?\n> In this project, we will consider how to explore, pre-process and encode a dataset of life insurance applicants, how to select important risk features from the dataset, how to train/calibrate/test a range of ML classifiers to predict risk ratings based on these selected features, how to evaluate the performance of each ML classifier, and how to understand/interpret predictions generated by an ML classifier.","metadata":{}},{"cell_type":"markdown","source":"# Code initialisation","metadata":{}},{"cell_type":"code","source":"# Import key modules that will be used throughout the project.\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt # graphs/plotting\nimport seaborn as sns\n\n# Check to ensure that both CSV files are held in the correct (input) directory.\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-08T17:51:27.296038Z","iopub.execute_input":"2022-07-08T17:51:27.297045Z","iopub.status.idle":"2022-07-08T17:51:27.838742Z","shell.execute_reply.started":"2022-07-08T17:51:27.296945Z","shell.execute_reply":"2022-07-08T17:51:27.837408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Step 1: Import the datasets into dataframes","metadata":{}},{"cell_type":"code","source":"# Load the CSV data file into a Pandas dataframe.\nmain_data = pd.read_csv('../input/prudential-life-insurance-assessment/train.csv.zip')\n\n# Review the main characteristics of the Pandas dataframe.\nprint(main_data.dtypes)\nmain_data.describe()","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:51:27.844181Z","iopub.execute_input":"2022-07-08T17:51:27.844977Z","iopub.status.idle":"2022-07-08T17:51:29.145512Z","shell.execute_reply.started":"2022-07-08T17:51:27.844945Z","shell.execute_reply":"2022-07-08T17:51:29.144690Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> From the cell above, we can see that there are a variety of datatypes within our dataframe - any columns with `object` dtype contain non-numerical (character) data, which will need to be pre-processed in order for these to be machine-interpretable.\n>\n> This will be explained in further detail later on.","metadata":{}},{"cell_type":"code","source":"# Set 'id' as the index column.\nmain_data_index_set = main_data.set_index('Id')","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:51:29.146985Z","iopub.execute_input":"2022-07-08T17:51:29.147365Z","iopub.status.idle":"2022-07-08T17:51:29.166039Z","shell.execute_reply.started":"2022-07-08T17:51:29.147333Z","shell.execute_reply":"2022-07-08T17:51:29.164800Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"main_data_index_set.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:51:29.170193Z","iopub.execute_input":"2022-07-08T17:51:29.170636Z","iopub.status.idle":"2022-07-08T17:51:29.196800Z","shell.execute_reply.started":"2022-07-08T17:51:29.170600Z","shell.execute_reply":"2022-07-08T17:51:29.195754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> Here, the `Id` has been set as the index, and each row lists out the values of each column (risk feature) for the first 5 life insurance applicants within the dataset.\n> \n> Note above that there are a mixture of numeric-valued features that are both normalised and non-normalised - we will need to review and pre-process each of these columns (where required), before training/testing our models further on.","metadata":{}},{"cell_type":"code","source":"# Visualise the distribution of applicants' risk ratings in the full dataset.\nplt.hist(main_data_index_set['Response'], bins=sorted(main_data_index_set['Response'].unique()))\nplt.xlabel('Response')\nplt.ylabel('# of Applicants')\nplt.title('Response Distribution')","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:51:29.198559Z","iopub.execute_input":"2022-07-08T17:51:29.198939Z","iopub.status.idle":"2022-07-08T17:51:29.444247Z","shell.execute_reply.started":"2022-07-08T17:51:29.198901Z","shell.execute_reply":"2022-07-08T17:51:29.442915Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The histogram above shows the full dataset's distribution of applicants by their determined risk ratings (i.e. `Response`).\n> \n> It is important to note that the distribution is unbalanced and is skewed towards classes 6-8, although classes 1-2 also account for a notable proportion of the dataset. We will see later on as to how closely the ML models are able to 'mimic' this distribution, as an indication of whether they have been adequately fitted.","metadata":{}},{"cell_type":"markdown","source":"# Step 2: Perform EDA to visualise distributions/outliers/correlations/excess zeroes","metadata":{}},{"cell_type":"markdown","source":"### Exploratory Data Analysis (EDA)\n\n> In this section, we will perform EDA in order to visualise and understand the dataset that we are working on. This is essential in order to understand how our models will need to be trained/evaluated/understood, later on.\n>\n> The cell below contains separate lists that pertain to each of the sets of variables within the dataset. For instance, `ColSet1_ProdInfo` refers to the set of columns that are prefixed by \"**Product_Info_**\".","metadata":{}},{"cell_type":"code","source":"# Split out the full set of the main dataset's columns into separate lists for easier use.\n\nColSet1_ProdInfo = ['Product_Info_1','Product_Info_2','Product_Info_3','Product_Info_4','Product_Info_5','Product_Info_6','Product_Info_7']\nColSet2_ApplicantInfo = ['Ins_Age','Ht','Wt','BMI']\nColSet3_EmploymentInfo = ['Employment_Info_1','Employment_Info_2','Employment_Info_3','Employment_Info_4','Employment_Info_5','Employment_Info_6']\nColSet4_InsuredInfo = ['InsuredInfo_1','InsuredInfo_2','InsuredInfo_3','InsuredInfo_4','InsuredInfo_5','InsuredInfo_6','InsuredInfo_7']\nColSet5_InsuranceHistoryInfo = ['Insurance_History_1','Insurance_History_2','Insurance_History_3','Insurance_History_4','Insurance_History_5','Insurance_History_7','Insurance_History_8','Insurance_History_9']\nColSet6_FamilyHistoryInfo = ['Family_Hist_1','Family_Hist_2','Family_Hist_3','Family_Hist_4','Family_Hist_5']\n\nColSet7_MedicalHistoryInfo = ['Medical_History_1','Medical_History_2','Medical_History_3','Medical_History_4','Medical_History_5','Medical_History_6','Medical_History_7','Medical_History_8',\n                              'Medical_History_9','Medical_History_10','Medical_History_11','Medical_History_12','Medical_History_13','Medical_History_14','Medical_History_15',\n                              'Medical_History_16','Medical_History_17','Medical_History_18','Medical_History_19','Medical_History_20','Medical_History_21','Medical_History_22',\n                              'Medical_History_23','Medical_History_24','Medical_History_25','Medical_History_26','Medical_History_27','Medical_History_28','Medical_History_29',\n                              'Medical_History_30','Medical_History_31','Medical_History_32','Medical_History_33','Medical_History_34','Medical_History_35','Medical_History_36',\n                              'Medical_History_37','Medical_History_38','Medical_History_39','Medical_History_40','Medical_History_41']\n\nColSet8_MedicalKeywordInfo = ['Medical_Keyword_1','Medical_Keyword_2','Medical_Keyword_3','Medical_Keyword_4','Medical_Keyword_5','Medical_Keyword_6','Medical_Keyword_7','Medical_Keyword_8',\n                              'Medical_Keyword_9','Medical_Keyword_10','Medical_Keyword_11','Medical_Keyword_12','Medical_Keyword_13','Medical_Keyword_14','Medical_Keyword_15','Medical_Keyword_16',\n                              'Medical_Keyword_17','Medical_Keyword_18','Medical_Keyword_19','Medical_Keyword_20','Medical_Keyword_21','Medical_Keyword_22','Medical_Keyword_23','Medical_Keyword_24',\n                              'Medical_Keyword_25','Medical_Keyword_26','Medical_Keyword_27','Medical_Keyword_28','Medical_Keyword_29','Medical_Keyword_30','Medical_Keyword_31','Medical_Keyword_32',\n                              'Medical_Keyword_33','Medical_Keyword_34','Medical_Keyword_35','Medical_Keyword_36','Medical_Keyword_37','Medical_Keyword_38','Medical_Keyword_39','Medical_Keyword_40',\n                              'Medical_Keyword_41','Medical_Keyword_42','Medical_Keyword_43','Medical_Keyword_44','Medical_Keyword_45','Medical_Keyword_46','Medical_Keyword_47','Medical_Keyword_48']","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:51:29.445957Z","iopub.execute_input":"2022-07-08T17:51:29.446285Z","iopub.status.idle":"2022-07-08T17:51:29.459857Z","shell.execute_reply.started":"2022-07-08T17:51:29.446257Z","shell.execute_reply":"2022-07-08T17:51:29.458568Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Distribution plots\n> Here, we will generate Kernel Density Estimate (KDE) plots for each of the features - where `Response` is set as the hue of each curve - in order to compare the distributions across each risk rating, and to understand whether there are any trends/correlations within the data.\n> \n> To do this, we will use the `seaborn.kdeplot()` function as a high-level interface to plot these pairwise relationships in the `main_data_index_set` dataset.","metadata":{}},{"cell_type":"code","source":"# Set up a subplot grid.\nfig, axes = plt.subplots(nrows=2, ncols=3, figsize=(25,15))\n\n# Product_Info_2 has been excluded - as this has not yet been encoded into numeric values.\nColSet1_ProdInfo_kde = ['Product_Info_1','Product_Info_3','Product_Info_4','Product_Info_5','Product_Info_6','Product_Info_7']\n\n# Produce kernel density estimate plots for each set of columns.\nfor i, column in enumerate(main_data_index_set[ColSet1_ProdInfo_kde].columns):\n    sns.kdeplot(data=main_data_index_set,\n                x=column,\n                hue=\"Response\", fill=True, common_norm=True, alpha=0.05,\n                ax=axes[i//3,i%3])","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:51:29.461415Z","iopub.execute_input":"2022-07-08T17:51:29.462684Z","iopub.status.idle":"2022-07-08T17:51:33.079340Z","shell.execute_reply.started":"2022-07-08T17:51:29.462636Z","shell.execute_reply":"2022-07-08T17:51:33.078072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> These KDE plots show a number of distributions with varying modality, however the consistent trend to them is that they all closely overlap between each `Response` group/cohort of applicants, with no major difference in relative densities. As a result, any variation in these features is unlikely to individually help towards predicting an applicant's risk rating.","metadata":{}},{"cell_type":"code","source":"# Set up a subplot grid.\nfig, axes = plt.subplots(nrows=2, ncols=2, figsize=(25,15))\n\n# Produce kernel density estimate plots for each set of columns.\nfor i, column in enumerate(main_data_index_set[ColSet2_ApplicantInfo].columns):\n    sns.kdeplot(data=main_data_index_set,\n                x=column,\n                hue=\"Response\", fill=True, common_norm=True, alpha=0.05,\n                ax=axes[i//2,i%2])","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:51:33.080964Z","iopub.execute_input":"2022-07-08T17:51:33.081739Z","iopub.status.idle":"2022-07-08T17:51:35.506259Z","shell.execute_reply.started":"2022-07-08T17:51:33.081679Z","shell.execute_reply":"2022-07-08T17:51:35.505016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* **Ins_Age**: This KDE plot displays significant variation in each `Response` group distribution's composition and structure, where each peak exhibits both broadening and shouldering. Whilst the majority of the distribution density is spread between x=0 and x=0.9, each cohort's distribution does vary somewhat in shape and skew - for instance, class 8 features positive skew towards x=0.2 whereas most of the other classes are negatively skewed towards x=0.4.\n* **Ht**: This KDE plot displays significant variation in each `Response` group distribution's composition and structure, where each peak exhibits both broadening and shouldering. The majority of the distribution density is spread between x=0.6 and x=0.8, and each cohort's skew is similar to that shown in `Ins_Age` (i.e. positive for class 8, negative for most others).\n* **Wt**: This KDE plot displays significant variation in each `Response` group distribution's composition and structure, where each peak exhibits both broadening and shouldering. The majority of the distribution density is spread between x=0.2 and x=0.5. The distribution for class 8 is notably centred around low values of `Wt` (x=0.2), whereas the remainder of the population is mostly spread across the interval between x=0.2 and x=0.5.\n* **BMI**: This KDE plot displays significant variation in each `Response` group distribution's composition and structure, where each peak exhibits both broadening and shouldering. The majority of the distribution density is spread in a similar fashion to `Wt` - class 8's peak is centred at x=0.4 whereas most of the other classes' distributions are spread across the interval between x=0.4 and x=0.6. However, one of the \"medium\" risk-rating distributions features a notably stronger skew in its distribution compared to its peers, with significant negative skew towards x=0.6 rather than towards x=0.5.","metadata":{}},{"cell_type":"code","source":"# Set up a subplot grid.\nfig, axes = plt.subplots(nrows=2, ncols=3, figsize=(25,15))\n\n# Produce kernel density estimate plots for each set of columns.\nfor i, column in enumerate(main_data_index_set[ColSet3_EmploymentInfo].columns):\n    sns.kdeplot(data=main_data_index_set,\n                x=column,\n                hue=\"Response\", fill=True, common_norm=True, alpha=0.05,\n                ax=axes[i//3,i%3])","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:51:35.507864Z","iopub.execute_input":"2022-07-08T17:51:35.509101Z","iopub.status.idle":"2022-07-08T17:51:39.137229Z","shell.execute_reply.started":"2022-07-08T17:51:35.509050Z","shell.execute_reply":"2022-07-08T17:51:39.135980Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> These KDE plots show a number of distributions with varying modality, however the consistent trend to them is that they all closely overlap between each `Response` group/cohort of applicants, with no major difference in relative densities. As a result, any variation in these features is unlikely to individually help towards predicting an applicant's risk rating.","metadata":{}},{"cell_type":"code","source":"# Set up a subplot grid.\nfig, axes = plt.subplots(nrows=3, ncols=3, figsize=(25,20))\n\n# Produce kernel density estimate plots for each set of columns.\nfor i, column in enumerate(main_data_index_set[ColSet4_InsuredInfo].columns):\n    sns.kdeplot(data=main_data_index_set,\n                x=column,\n                hue=\"Response\", fill=True, common_norm=True, alpha=0.05,\n                ax=axes[i//3,i%3])\n\n# Delete any unused sets of axes in the subplot grid.    \nfig.delaxes(axes[2,1])\nfig.delaxes(axes[2,2])","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:51:39.138761Z","iopub.execute_input":"2022-07-08T17:51:39.139117Z","iopub.status.idle":"2022-07-08T17:51:43.154824Z","shell.execute_reply.started":"2022-07-08T17:51:39.139085Z","shell.execute_reply":"2022-07-08T17:51:43.153620Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> These KDE plots show a number of distributions with varying modality, however the consistent trend to them is that they all closely overlap between each `Response` group/cohort of applicants, with no major difference in relative densities. As a result, any variation in these features is unlikely to individually help towards predicting an applicant's risk rating.","metadata":{}},{"cell_type":"code","source":"# Set up a subplot grid.\nfig, axes = plt.subplots(nrows=3, ncols=3, figsize=(25,20))\n\n# Produce kernel density estimate plots for each set of columns.\nfor i, column in enumerate(main_data_index_set[ColSet5_InsuranceHistoryInfo].columns):\n    sns.kdeplot(data=main_data_index_set,\n                x=column,\n                hue=\"Response\", fill=True, common_norm=True, alpha=0.05,\n                ax=axes[i//3,i%3])\n    \n# Delete any unused sets of axes in the subplot grid.    \nfig.delaxes(axes[2,2])","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:51:43.156331Z","iopub.execute_input":"2022-07-08T17:51:43.156691Z","iopub.status.idle":"2022-07-08T17:51:47.652253Z","shell.execute_reply.started":"2022-07-08T17:51:43.156652Z","shell.execute_reply":"2022-07-08T17:51:47.650791Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> These KDE plots show a number of distributions with varying modality, however the consistent trend to them is that they all closely overlap between each `Response` group/cohort of applicants, with no major difference in relative densities. As a result, any variation in these features is unlikely to individually help towards predicting an applicant's risk rating.","metadata":{}},{"cell_type":"code","source":"# Set up a subplot grid.\nfig, axes = plt.subplots(nrows=2, ncols=3, figsize=(25,15))\n\n# Produce kernel density estimate plots for each set of columns.\nfor i, column in enumerate(main_data_index_set[ColSet6_FamilyHistoryInfo].columns):\n    sns.kdeplot(data=main_data_index_set,\n                x=column,\n                hue=\"Response\", fill=True, common_norm=True, alpha=0.05,\n                ax=axes[i//3,i%3])\n\n# Delete any unused sets of axes in the subplot grid.    \nfig.delaxes(axes[1,2])","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:51:47.654021Z","iopub.execute_input":"2022-07-08T17:51:47.654460Z","iopub.status.idle":"2022-07-08T17:51:49.707135Z","shell.execute_reply.started":"2022-07-08T17:51:47.654412Z","shell.execute_reply":"2022-07-08T17:51:49.705951Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> * **Family_Hist_1:** This KDE plot shows a bimodal distribution comprised of two peaks at x=2 and x=3 (along with a very low-density curve at x=1), with the most prominent peak at x=3. However, as this feature appears to show little variation in relative densities between each of the `Response` classes, any variation in this feature is unlikely to individually help towards predicting an applicant's risk rating.\n> \n> * **Family_Hist_2 - Family_Hist_5**: These KDE plots show a number of unimodal distributions which display some variation in terms of each `Response` group distribution's composition and structure, where each peak exhibits both broadening and shouldering. Whilst the majority of the distributions' densities are spread between x=0.2 and x=0.8, each cohort's distribution does vary somewhat in shape and skew/kurtosis - for instance, class 8 generally features more positive kurtosis than most of the other classes' distributions, which tend to show broader density plots.","metadata":{}},{"cell_type":"code","source":"# Set up a subplot grid.\nfig, axes = plt.subplots(nrows=11, ncols=4, figsize=(25,75))\n\n# Produce kernel density estimate plots for each set of columns.\nfor i, column in enumerate(main_data_index_set[ColSet7_MedicalHistoryInfo].columns):\n    sns.kdeplot(data=main_data_index_set,\n                x=column,\n                hue=\"Response\", fill=True, common_norm=True, alpha=0.05,\n                ax=axes[i//4,i%4])\n\n# Delete any unused sets of axes in the subplot grid.    \nfig.delaxes(axes[10,1])\nfig.delaxes(axes[10,2])\nfig.delaxes(axes[10,3])","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:51:49.712398Z","iopub.execute_input":"2022-07-08T17:51:49.712788Z","iopub.status.idle":"2022-07-08T17:52:12.149524Z","shell.execute_reply.started":"2022-07-08T17:51:49.712756Z","shell.execute_reply":"2022-07-08T17:52:12.148320Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A majority of these KDE plots feature distributions that closely overlap between each `Response` group/cohort of applicants, and hence any variation in the underlying features is unlikely to individually lend any predictive power for determining an applicant's risk rating. However, there are a handful of notable exceptions:\n\n* **Medical_History_2/15/24**: These KDE plots appears to show multimodal distributions that features some degree of predictive distinction in terms of variance, as each Response group distribution's peaks exhibit different levels of broadening and shouldering. However, it is important to note the scales of the y-axes here - the densities for each underlying distribution are very small and provide little in the way of allowing each feature to individually help in distinguishing between each Response group.\n* **Medical_History_10**: This KDE plot seems to display some interesting features. For the low-risk applicant cohorts, the plots appear to show distributions that are bimodal, however at higher risk levels the distribution tends towards having a single peak instead. Unfortunately, as will be shown in further detail later on, this column features a very high proportion of missing values and thus its distribution/s should not be misconstrued as highly predictive.\n* **Medical_History_23**: This KDE plot helps to illustrate that this feature *does* show some potential. As the `Response` value/risk rating increases, each peak in the bimodal distribution becomes sharper as the degree of kurtosis becomes more positive. In simpler terms, having a value (within this column) that is further away from the peaks' centres would tend to correlate with having a lower risk rating, whereas values that overlap closely with the peaks tend to represent applicants with higher risk ratings.","metadata":{}},{"cell_type":"code","source":"# Set up a subplot grid.\nfig, axes = plt.subplots(nrows=12, ncols=4, figsize=(25,50))\n\n# Produce kernel density estimate plots for each set of columns.\nfor i, column in enumerate(main_data_index_set[ColSet8_MedicalKeywordInfo].columns):\n    sns.kdeplot(data=main_data_index_set,\n                x=column,\n                hue=\"Response\", fill=True, common_norm=True, alpha=0.05,\n                ax=axes[i//4,i%4])","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:12.151246Z","iopub.execute_input":"2022-07-08T17:52:12.151623Z","iopub.status.idle":"2022-07-08T17:52:41.054332Z","shell.execute_reply.started":"2022-07-08T17:52:12.151589Z","shell.execute_reply":"2022-07-08T17:52:41.053085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> These KDE plots show a number of distributions with varying modality, however the consistent trend to them is that they all closely overlap between each `Response` group/cohort of applicants, with no major difference in relative densities. As a result, any variation in these features is unlikely to individually help towards predicting an applicant's risk rating.","metadata":{}},{"cell_type":"markdown","source":"### Correlation plots/heatmaps\n\n> Next, we will generate a correlation heatmap plot for all of the features - in order to better understand the correlations between each pair of features, and help to uncover whether there are any possible interactions within the data.\n> \n> To do this, we will use the `seaborn.heatmap()` function as a high-level interface to plot these pairwise relationships in the `main_data_index_set` dataset.","metadata":{}},{"cell_type":"code","source":"# Produce a correlation matrix of the dataset - then, create a mask to hide the upper-right half of the matrix.\ncorrs = main_data_index_set.corr()\nmask = np.zeros_like(corrs)\nmask[np.triu_indices_from(mask)] = True\n\n# Convert the correlation matrix into a heatmap using Seaborn.\nplt.figure(figsize=(24,16))\nsns.heatmap(corrs, cmap='RdBu_r', mask=mask)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:41.056249Z","iopub.execute_input":"2022-07-08T17:52:41.056989Z","iopub.status.idle":"2022-07-08T17:52:46.064820Z","shell.execute_reply.started":"2022-07-08T17:52:41.056918Z","shell.execute_reply":"2022-07-08T17:52:46.063592Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From the chart above, we can infer (based on each column set):\n\n* **Column Set 1 - Product Info:** These appear to show little interaction/correlation with a majority of the other feature sets, with the exception of `Employment_Info_1`, `Employment_Info_5` and `Insured_Info_6` - these columns may be directly correlated as, to give an example, an applicant's employment/financial status will have an impact on what type of policy/product the applicant is applying for.\n* **Column Set 2 - Applicant Info:** These columns show a varying range of interactions with the other feature sets; the strongest anti-/correlations (excluding those within the same column set) are between a handful of the `Family_Hist` columns as well as `Insured_Info_6`.\n* **Column Set 3 - Employment Info:** With the exception of two strong anti-correlations - between `Employment_Info_2` and `Employment_Info_3`, plus between `Employment_Info_5` and `Product_Info_3` - and a couple of moderate interactions between `Employment_Info_6` and `Family_Hist_2` / `Family_Hist_4`, this column set does not interact very strongly with the rest of the features.\n* **Column Set 4 - Insured Info:** The column `InsuredInfo_2` shows a fairly strong correlation with `InsuredInfo_7`, and also a strong anticorrelation with some of the Applicant Info columns; otherwise, this column set does not interact very much with the rest of the other features.\n* **Column Set 5 - Insurance History Info:** This feature set exhibits several strong inter-correlations with other `Insurance_History` columns, but does not interact very much with the rest of the features.\n* **Column Set 6 - Family History Info:** The columns `Family_Hist_2` and `Family_Hist_4` show a very strong positive correlation with `Ins_Age`, and also with `Medical_History_10` and `Medical_History_15` to a lesser degree.\n* **Column Set 7 - Medical History Info:** This column set shows a number of correlation hotspots against several `Medical_Keyword` columns, as well as against `Ins_Age` and some of the `Family_Hist` columns.\n* **Column Set 8 - Medical Keyword Info:** These columns show a number of correlation hotspots against several `Medical_History` columns, but do not otherwise show any notable interactions with the rest of the features.","metadata":{}},{"cell_type":"markdown","source":"### Missing values/excess zeroes\n\n> It is important to review the completeness of our dataset, in order to make sure that any correlations/trends can be substantiated with sufficient evidence - in other words, that there are enough data-points to validate any inferences/predictions that we may make from the dataset.\n> \n> For demonstration purposes, the following checks below are aimed at displaying where the entire dataset is incomplete.","metadata":{}},{"cell_type":"code","source":"# Determine which columns contain nulls/missing values.\ncols_with_missing = [col for col in main_data_index_set.columns\n                     if main_data_index_set[col].isnull().any()]\n\n# Summarise how many missing values are present in each column.\nmain_data_index_set[cols_with_missing].isna().sum()","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:46.066224Z","iopub.execute_input":"2022-07-08T17:52:46.066593Z","iopub.status.idle":"2022-07-08T17:52:46.107893Z","shell.execute_reply.started":"2022-07-08T17:52:46.066561Z","shell.execute_reply":"2022-07-08T17:52:46.106929Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Calculate the proportion of zeroes relative to non-zero values.\nfor col in cols_with_missing:\n    sum = main_data_index_set[col].isna().sum()\n    length = len(main_data_index_set[col].index)\n    ratio = sum/length\n    print('Proportion of zeroes in', col, 'is: ', round(ratio*100,2), '%.')","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:46.109318Z","iopub.execute_input":"2022-07-08T17:52:46.109617Z","iopub.status.idle":"2022-07-08T17:52:46.123790Z","shell.execute_reply.started":"2022-07-08T17:52:46.109589Z","shell.execute_reply":"2022-07-08T17:52:46.122195Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> In the code above, we have checked through the entire dataset to determine how many columns contain missing values, and how many missing values there are in each of these columns.\n> \n> However, as we wish to avoid test-data leakage later on during our model validation/testing process, we need to ensure that we are only using a \"training\" set of data in order to inform our view of the completeness of our data. Hence, we will need to split the full dataset into training/validation/testing sets *and then* review the training set for missing values afterwards.","metadata":{}},{"cell_type":"markdown","source":"# Step 3: Perform a train-validation-test split of the dataset\n\n> Here, we will create separate dataframes that will store the features and target variables.\n>\n> These are then supplied to the `sklearn` function `train_test_split()` in order to split the data into training/validation/test subsets.","metadata":{}},{"cell_type":"code","source":"# Assign the features to their own dataframe.\nX_full = main_data_index_set.drop(['Response'], axis=1)\n\n# Assign the target variable to its own dataframe.\ny_full = main_data_index_set.Response\n\n# Perform a train-test split to obtain the training, validation and test data as separate dataframes.\nfrom sklearn.model_selection import train_test_split\n\n# Split out test/holdout set from full dataset.\n# We will set the size of the X/y test datasets to be 20% of the original (full) X/y datasets, via the train_size/test_size parameters.\nX_rem, X_test, y_rem, y_test = train_test_split(X_full, y_full, train_size=0.8, test_size=0.2, random_state=0, stratify=y_full)\n\n# Split remaining portion into training/validation sets.\n# We will set the size of the X/y train datasets to be 60% of the original (full) X/y datasets, via the train_size/test_size parameters.\nX_train, X_valid, y_train, y_valid = train_test_split(X_rem, y_rem, train_size=0.75, test_size=0.25, random_state=0, stratify=y_rem)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:46.125325Z","iopub.execute_input":"2022-07-08T17:52:46.126286Z","iopub.status.idle":"2022-07-08T17:52:46.446494Z","shell.execute_reply.started":"2022-07-08T17:52:46.126246Z","shell.execute_reply":"2022-07-08T17:52:46.445219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Step 4: Review the training dataset's columns and handle them accordingly","metadata":{}},{"cell_type":"markdown","source":"### Locating/handling excess zeroes\n> As detailed above, we will now review the training subset in order to understand where there are missing values.\n>\n> Statistical tests can be performed to assess whether values in these columns are either missing completely at random (MCAR), missing at random (MAR) or missing not at random (MNAR) - this will indicate whether or not our remaining data is likely to still be representative of the general population. In the interest of brevity, we assume here that any missing data are MCAR and that any subsequent analysis/imputation is not subject to implicit bias.","metadata":{}},{"cell_type":"code","source":"# Determine which columns contain nulls/missing values.\nX_train_cols_with_missing = [col for col in X_train.columns\n                     if X_train[col].isnull().any()]\n\n# Summarise how many missing values are present in each column.\nX_train[X_train_cols_with_missing].isna().sum()","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:46.447779Z","iopub.execute_input":"2022-07-08T17:52:46.448160Z","iopub.status.idle":"2022-07-08T17:52:46.491252Z","shell.execute_reply.started":"2022-07-08T17:52:46.448125Z","shell.execute_reply":"2022-07-08T17:52:46.489977Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Calculate the proportion of zeroes relative to non-zero values.\nfor col in X_train_cols_with_missing:\n    sum = X_train[col].isna().sum()\n    length = len(X_train[col].index)\n    ratio = sum/length\n    print('Proportion of zeroes in', col, 'is: ', round(ratio*100,2), '%.')","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:46.492637Z","iopub.execute_input":"2022-07-08T17:52:46.492985Z","iopub.status.idle":"2022-07-08T17:52:46.506603Z","shell.execute_reply.started":"2022-07-08T17:52:46.492954Z","shell.execute_reply":"2022-07-08T17:52:46.504684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> A number of columns have been identified in the code above, with different proportions of missing values - some are still in workable condition and can be preprocessed via imputation methods in order to provide machine-interpretable inputs for our models. However, there are still a handful of columns which are highly incomplete - any attempt to perform imputation would likely introduce significant bias/skew into these features' distributions.\n>\n> Under these conditions, it may be safer to simply remove the columns altogether - intuitively, this makes sense as the number of *non-blank* values already represents a very small portion of these columns. Here, we have elected to delete columns where their proportions of missing values in the *training* subset are **greater than 40%**, although any other sensible threshold could be set instead.","metadata":{}},{"cell_type":"code","source":"# These columns have been selected as they contain a high proportion of blanks/missing values (deemed here as >40%) in the TRAINING dataset.\ncols_to_delete_due_to_missing_data = ['Insurance_History_5',\n                                      'Family_Hist_2', 'Family_Hist_3', 'Family_Hist_5',\n                                      'Medical_History_10', 'Medical_History_15', 'Medical_History_24', 'Medical_History_32']\n\n# Delete columns from ALL datasets where the proportion of zeroes in the TRAINING dataset exceeds a stipulated threshold.\nX_train = X_train.drop(cols_to_delete_due_to_missing_data, axis=1)\nX_valid = X_valid.drop(cols_to_delete_due_to_missing_data, axis=1)\nX_test = X_test.drop(cols_to_delete_due_to_missing_data, axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:46.508479Z","iopub.execute_input":"2022-07-08T17:52:46.509290Z","iopub.status.idle":"2022-07-08T17:52:46.545099Z","shell.execute_reply.started":"2022-07-08T17:52:46.509240Z","shell.execute_reply":"2022-07-08T17:52:46.543795Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Handling columns with missing values - via iterative imputation\n> Here, we perform iterative imputation on the remaining columns with missing values, which involves estimating feature values as a function of all other features, such that any gaps are instead replaced with predicted values.\n>\n> This means that all rows/columns within the dataset will now contain machine-interpretable values, ready for use in feature scaling/selection as well as model fitting.","metadata":{}},{"cell_type":"code","source":"# These columns still contain missing values, and require imputation before they can be supplied as inputs to each ML classifier.\ncols_to_impute = ['Employment_Info_1', 'Employment_Info_4', 'Employment_Info_6',\n                 'Family_Hist_4',\n                 'Medical_History_1']","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:46.546597Z","iopub.execute_input":"2022-07-08T17:52:46.547449Z","iopub.status.idle":"2022-07-08T17:52:46.553009Z","shell.execute_reply.started":"2022-07-08T17:52:46.547406Z","shell.execute_reply":"2022-07-08T17:52:46.551804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Import the IterativeImputer class from sklearn - NOTE: enable_iterative_imputer also needs to be imported as this is an experimental feature. \nfrom sklearn.experimental import enable_iterative_imputer\nfrom sklearn.impute import IterativeImputer\n\n# Take a copy of each dataset before transforming.\ncopy_X_train = X_train.copy()\ncopy_X_valid = X_valid.copy()\ncopy_X_test = X_test.copy()\n\n# Filter the splits down to the columns that require imputation.\nX_train_pre_impute = copy_X_train[cols_to_impute]\nX_valid_pre_impute = copy_X_valid[cols_to_impute]\nX_test_pre_impute = copy_X_test[cols_to_impute]\n\n# Save the other columns into separate dataframes, for re-joining later on.\nX_train_no_impute = copy_X_train.drop(cols_to_impute, axis=1)\nX_valid_no_impute = copy_X_valid.drop(cols_to_impute, axis=1)\nX_test_no_impute = copy_X_test.drop(cols_to_impute, axis=1)\n\n# Initialise the IterativeImputer transformer.\nX_imputer = IterativeImputer(random_state=0)\n\n# Transform the train/val/test datasets using iterative imputation.\nX_train_post_impute = pd.DataFrame(X_imputer.fit_transform(X_train_pre_impute), columns=X_train_pre_impute.columns)\nX_valid_post_impute = pd.DataFrame(X_imputer.transform(X_valid_pre_impute), columns=X_valid_pre_impute.columns)\nX_test_post_impute = pd.DataFrame(X_imputer.transform(X_test_pre_impute), columns=X_train_pre_impute.columns)\n\n# Reset the indexes of each dataset, as they are dropped during imputation.\nX_train_post_impute.index = X_train_pre_impute.index\nX_valid_post_impute.index = X_valid_pre_impute.index\nX_test_post_impute.index = X_test_pre_impute.index\n\n# Re-join the imputed columns with the remaining columns in each dataset.\nX_train_imputed = pd.concat([X_train_no_impute, X_train_post_impute], axis=1)\nX_valid_imputed = pd.concat([X_valid_no_impute, X_valid_post_impute], axis=1)\nX_test_imputed = pd.concat([X_test_no_impute, X_test_post_impute], axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:46.555061Z","iopub.execute_input":"2022-07-08T17:52:46.555790Z","iopub.status.idle":"2022-07-08T17:52:47.314090Z","shell.execute_reply.started":"2022-07-08T17:52:46.555744Z","shell.execute_reply":"2022-07-08T17:52:47.312658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Checks - before/after iterative imputation\n\n> We will now review the dataset before and after iterative imputation, in order to understand how each column's distribution has been affected.","metadata":{}},{"cell_type":"code","source":"# KDE Plots - Before imputation.\n\nfig, axes = plt.subplots(nrows=2, ncols=3, figsize=(25,10))\n\nfor i, column in enumerate(main_data_index_set[X_train_post_impute.columns].columns):\n    sns.kdeplot(data=main_data_index_set[X_train_post_impute.columns],\n                x=column,\n                fill=True, common_norm=True, alpha=0.05,\n                ax=axes[i//3,i%3])\n    \nfig.delaxes(axes[1,2])","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:47.316200Z","iopub.execute_input":"2022-07-08T17:52:47.317121Z","iopub.status.idle":"2022-07-08T17:52:49.241936Z","shell.execute_reply.started":"2022-07-08T17:52:47.317070Z","shell.execute_reply":"2022-07-08T17:52:49.240868Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# KDE Plots - After imputation.\n\nfig, axes = plt.subplots(nrows=2, ncols=3, figsize=(25,10))\n\nfor i, column in enumerate(X_train_post_impute.columns):\n    sns.kdeplot(data=X_train_post_impute,\n                x=column,\n                fill=True, common_norm=True, alpha=0.05,\n                ax=axes[i//3,i%3])\n    \nfig.delaxes(axes[1,2])","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:49.243498Z","iopub.execute_input":"2022-07-08T17:52:49.244394Z","iopub.status.idle":"2022-07-08T17:52:51.264382Z","shell.execute_reply.started":"2022-07-08T17:52:49.244354Z","shell.execute_reply":"2022-07-08T17:52:51.263212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> As shown above, three of the five columns' distributions appear to be mostly unchanged - however, we have introduced additional probability density/peak splitting into `Family_Hist_4` (at approx. x=0.4, x=0.6) and `Medical_History_1` (at approx. x=10).\n>\n> We will need to take care later on, when evaluating the importance of these two features and explaining how we have generated predictions for our chosen model.","metadata":{}},{"cell_type":"markdown","source":"# Step 5: Review and consider any other sources of data leakage\n\n> This dataset is comprised of over a hundred variables describing various keywords/attributes of life insurance applicants (i.e. information that they would declare when applying for life insurance). As all of this information will be available at the time of predicting the `Response` variable, the dataset does not feature any sources of target leakage.\n> \n> Furthermore, we have used `train_test_split()` in order to split our dataset into training/validation/test samples, so that our models can be assembled/evaluated in a fair manner that does not invoke any form of data snooping/train-test contamination.","metadata":{}},{"cell_type":"markdown","source":"# Step 6: Perform feature engineering (using training data)\n\n> In this section, we will implement both supervised (e.g. categorical encoding) and unsupervised (e.g. K-means clustering) learning techniques in order to create new features within our dataset that our models can be fitted to.","metadata":{}},{"cell_type":"code","source":"# Create a clone copy of each imputed dataset to avoid changing any original data.\ncopy_X_train_imputed = X_train_imputed.copy()\ncopy_X_valid_imputed = X_valid_imputed.copy()\ncopy_X_test_imputed = X_test_imputed.copy()","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:51.265820Z","iopub.execute_input":"2022-07-08T17:52:51.266160Z","iopub.status.idle":"2022-07-08T17:52:51.336468Z","shell.execute_reply.started":"2022-07-08T17:52:51.266130Z","shell.execute_reply":"2022-07-08T17:52:51.335107Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### One-Hot Encoding\n\n> We will one-hot encode the `Product_Info_2` column, such that any categorical inputs are converted to numerical, machine-interpretable values that can be supplied to each classification model.","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import OneHotEncoder\n\n# Initialise a one-hot encoder to columns that contain categorical data.\nOH_encoder = OneHotEncoder(handle_unknown='ignore', sparse=False)\nOH_col = ['Product_Info_2']\n\n## We set handle_unknown='ignore' to avoid errors when the validation data contains classes that aren't represented\n## in the training data, and setting sparse=False ensures that the encoded columns are returned as a numpy array\n## (instead of a sparse matrix).\n\n# Use the one-hot encoder to transform the categorical data columns. \nOH_col_train = pd.DataFrame(OH_encoder.fit_transform(copy_X_train_imputed[OH_col]))\nOH_col_valid = pd.DataFrame(OH_encoder.transform(copy_X_valid_imputed[OH_col]))\nOH_col_test = pd.DataFrame(OH_encoder.transform(copy_X_test_imputed[OH_col]))\n\n# One-hot encoding removes the index; re-assign the original index.\nOH_col_train.index = copy_X_train_imputed.index\nOH_col_valid.index = copy_X_valid_imputed.index\nOH_col_test.index = copy_X_test_imputed.index\n\n# Add column-labelling back in, using the get_feature_names() function. \nOH_col_train.columns = OH_encoder.get_feature_names(OH_col)\nOH_col_valid.columns = OH_encoder.get_feature_names(OH_col)\nOH_col_test.columns = OH_encoder.get_feature_names(OH_col)\n\n# Create dataframes that only include the numerical features/columns (these will be concatenated with the one-hot encoded dataframes).\ncopy_X_train_imputed_no_OH_col = copy_X_train_imputed.drop(OH_col, axis=1)\ncopy_X_valid_imputed_no_OH_col = copy_X_valid_imputed.drop(OH_col, axis=1)\ncopy_X_test_imputed_no_OH_col = copy_X_test_imputed.drop(OH_col, axis=1)\n\n# Concatenate the one-hot encoded columns with the existing numerical features/columns.\nX_train_enc = pd.concat([copy_X_train_imputed_no_OH_col, OH_col_train], axis=1)\nX_valid_enc = pd.concat([copy_X_valid_imputed_no_OH_col, OH_col_valid], axis=1)\nX_test_enc = pd.concat([copy_X_test_imputed_no_OH_col, OH_col_test], axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:51.338188Z","iopub.execute_input":"2022-07-08T17:52:51.338543Z","iopub.status.idle":"2022-07-08T17:52:51.443002Z","shell.execute_reply.started":"2022-07-08T17:52:51.338509Z","shell.execute_reply":"2022-07-08T17:52:51.441679Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Scaling/normalisation\n> Next, we perform min-max scaling on the encoded datasets, such that all features lie between 0 and 1 - this is so that, when training any of the classification models, all features will have variances with the same order of magnitude as each other. Thus, no single feature will dominate the objective function and prohibit the model from learning from other features correctly as expected.\n","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import MinMaxScaler\n\n# Initialise the MinMaxScaler model, then fit it to the (encoded) training dataset.\nMM_scaler = MinMaxScaler()\nMM_scaler.fit(X_train_enc)\n\n# Then, normalise/transform the training, validation and test datasets.\nX_train_scale = pd.DataFrame(MM_scaler.transform(X_train_enc),\n                             index=X_train_enc.index,\n                             columns=X_train_enc.columns)\n\nX_valid_scale = pd.DataFrame(MM_scaler.transform(X_valid_enc),\n                             index=X_valid_enc.index,\n                             columns=X_valid_enc.columns)\n\nX_test_scale = pd.DataFrame(MM_scaler.transform(X_test_enc),\n                             index=X_test_enc.index,\n                             columns=X_test_enc.columns)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:51.444513Z","iopub.execute_input":"2022-07-08T17:52:51.444896Z","iopub.status.idle":"2022-07-08T17:52:51.620966Z","shell.execute_reply.started":"2022-07-08T17:52:51.444860Z","shell.execute_reply":"2022-07-08T17:52:51.619037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### K-means clustering\n> We will now use an unsupervised learning technique known as K-means clustering in order to group together applicants based on their commonalities. In order to do this, we will train our K-means clustering algorithm using the training subset, before predicting cluster groups on all three subsets.\n>\n> These cluster labels will be incorporated as an additional feature in our datasets, and may in fact prove to be useful in helping to understand applicants' risk rating assignments later on. ","metadata":{}},{"cell_type":"code","source":"# Create copies of the scaled datasets, prior to performing K-Means clustering.\ncopy_X_train_scale = X_train_scale.copy()\ncopy_X_valid_scale = X_valid_scale.copy()\ncopy_X_test_scale = X_test_scale.copy()\n\nfrom sklearn.cluster import KMeans","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:51.625198Z","iopub.execute_input":"2022-07-08T17:52:51.625657Z","iopub.status.idle":"2022-07-08T17:52:51.770503Z","shell.execute_reply.started":"2022-07-08T17:52:51.625617Z","shell.execute_reply":"2022-07-08T17:52:51.768865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> We will now use the elbow-analysis method for determining the optimal number of clusters to use.","metadata":{}},{"cell_type":"code","source":"# Determine the optimal number of clusters.\n# Method: Cluster the dataset into k clusters, then calculate the inertia/sum of squared distances.\n# Repeat this by looping through k=1 to k=30.\n\nSum_of_squared_distances = []\nK = range(1,30)\nfor k in K:\n    km = KMeans(n_clusters=k)\n    km = km.fit(copy_X_train_scale)\n    Sum_of_squared_distances.append(km.inertia_)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:52:51.772468Z","iopub.execute_input":"2022-07-08T17:52:51.773094Z","iopub.status.idle":"2022-07-08T17:55:35.955086Z","shell.execute_reply.started":"2022-07-08T17:52:51.773052Z","shell.execute_reply":"2022-07-08T17:55:35.953788Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create a plot of K-values versus their respective inertias/sums of squared distances.\nplt.plot(K, Sum_of_squared_distances, 'bx-')\nplt.xlabel('k')\nplt.ylabel('Sum_of_squared_distances')\nplt.title('Elbow Method For Optimal k')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:55:35.956733Z","iopub.execute_input":"2022-07-08T17:55:35.957318Z","iopub.status.idle":"2022-07-08T17:55:36.161449Z","shell.execute_reply.started":"2022-07-08T17:55:35.957281Z","shell.execute_reply":"2022-07-08T17:55:36.160165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> As shown above, the \"elbow\" of the curve begins to form at *k*=15.\n> \n> Increasing *k* beyond this value does not yield a significant benefit in the rate of reduction in the training dataset's inertia - hence, we will set *k*=15 when initialising the `KMeans()` clustering algorithm below.","metadata":{}},{"cell_type":"code","source":"# Set n_clusters=15, as derived from the elbow-method analysis above:\nkmeans = KMeans(n_clusters=15, n_init=10, random_state=0)\n\n# Fit the K-Means clustering algorithm to the training dataset, then predict the train/valid/test datasets.\ncopy_X_train_scale[\"Cluster\"] = kmeans.fit_predict(copy_X_train_scale)\ncopy_X_valid_scale[\"Cluster\"] = kmeans.predict(copy_X_valid_scale)\ncopy_X_test_scale[\"Cluster\"] = kmeans.predict(copy_X_test_scale)\n\n# Convert the cluster labels into one-hot encoded variants.\nX_train_cluster_OH_enc = pd.get_dummies(copy_X_train_scale.Cluster).add_prefix('KMeansCluster_')\nX_valid_cluster_OH_enc = pd.get_dummies(copy_X_valid_scale.Cluster).add_prefix('KMeansCluster_')\nX_test_cluster_OH_enc = pd.get_dummies(copy_X_test_scale.Cluster).add_prefix('KMeansCluster_')\n\n# Re-join the K-Means clustering labels onto the original dataframes.\nX_train_KMeans = pd.concat([copy_X_train_scale, X_train_cluster_OH_enc], axis=1)\nX_valid_KMeans = pd.concat([copy_X_valid_scale, X_valid_cluster_OH_enc], axis=1)\nX_test_KMeans = pd.concat([copy_X_test_scale, X_test_cluster_OH_enc], axis=1)\n\n# Remove the initially derived \"Cluster\" columns from each dataset.\nX_train_KMeans = X_train_KMeans.drop(['Cluster'], axis=1)\nX_valid_KMeans = X_valid_KMeans.drop(['Cluster'], axis=1)\nX_test_KMeans = X_test_KMeans.drop(['Cluster'], axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:55:36.162808Z","iopub.execute_input":"2022-07-08T17:55:36.163115Z","iopub.status.idle":"2022-07-08T17:55:43.020006Z","shell.execute_reply.started":"2022-07-08T17:55:36.163085Z","shell.execute_reply":"2022-07-08T17:55:43.018934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Step 7: Perform feature selection (using validation data)\n\n> In this section, we will use a number of different methods to determine which features within our dataset will be most important for fitting each of the classification models - this is done in order to prevent overfitting.\n>\n> At this stage, it is important to now switch over to using the validation dataset for feature selection/model refinement, so that we do not continually rely on the training dataset and risk invoking data leakage into our model generation process.","metadata":{}},{"cell_type":"markdown","source":"### Method 1: Mutual Information (MI)\n\n> Here, we will calculate the MI scores of the validation dataset, using two custom functions that rely on the `mutual_info_classif()` function available within `sklearn`.\n>\n> This will help us to understand whether there are any useful features in our dataset that should be preserved, during feature selection.","metadata":{}},{"cell_type":"code","source":"# Import the mutual_info_classif() class from sklearn.\nfrom sklearn.feature_selection import mutual_info_classif\n\n# Define a custom function that calculates Mutual Information (MI) scores for a given dataset.\ndef make_mi_scores(X, y):\n    mi_scores = mutual_info_classif(X, y)\n    mi_scores = pd.Series(mi_scores, name=\"MI Scores\", index=X.columns)\n    mi_scores = mi_scores.sort_values(ascending=False)\n    return mi_scores\n\n# Define a custom function that plots MI scores in descending order (i.e. most important to least important).\ndef plot_mi_scores(scores):\n    scores = scores.sort_values(ascending=True)\n    width = np.arange(len(scores))\n    ticks = list(scores.index)\n    plt.barh(width, scores)\n    plt.yticks(width, ticks)\n    plt.title(\"Mutual Information Scores\")","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:55:43.021883Z","iopub.execute_input":"2022-07-08T17:55:43.022631Z","iopub.status.idle":"2022-07-08T17:55:43.046295Z","shell.execute_reply.started":"2022-07-08T17:55:43.022580Z","shell.execute_reply":"2022-07-08T17:55:43.045029Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate MI scores on the validation dataset.\nmi_scores_X_valid = make_mi_scores(X_valid_KMeans, y_valid)\n\n# Plot the MI scores obtained from the validation dataset.\nplt.figure(dpi=100, figsize=(20,50))\nplot_mi_scores(mi_scores_X_valid)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:55:43.048937Z","iopub.execute_input":"2022-07-08T17:55:43.050296Z","iopub.status.idle":"2022-07-08T17:55:54.576958Z","shell.execute_reply.started":"2022-07-08T17:55:43.050196Z","shell.execute_reply":"2022-07-08T17:55:54.575597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> As per the chart above, we can see that the top 5 ranked features (in descending order) are `BMI`, `Wt`, `Product_Info_4`, `Medical_Keyword_15`, and `Medical_History_23`. This means that these features have strong statistical dependences with the `Response` variable, i.e. that they contribute significantly to reducing uncertainty in the value of the `Response` variable (given a known value of the feature).\n>\n> Features that have low/zero MI scores indicate that they do not significantly contribute towards reducing this uncertainty, and are hence less useful for guiding our predictions.","metadata":{}},{"cell_type":"markdown","source":"### Method 2: Multicollinearity analysis\n\n> We will now review our dataset via Variance Inflation Factor analysis, which is critical for detecting the presence of multicollinearity (i.e. where several independent variables in a model are highly correlated - hence resulting in less reliable statistical inferences).\n>\n> This will help us to understand whether there are any redundant features in our dataset that *should* not be kept, during feature selection.","metadata":{}},{"cell_type":"code","source":"# Import variance_inflation_factor from statsmodels\nfrom statsmodels.stats.outliers_influence import variance_inflation_factor\n\n# Define a custom function that calculates variance inflation factor (VIF) scores - for determining multicollinearity.\ndef calc_vif(X):\n    vif = pd.DataFrame()\n    vif[\"Variables\"] = X.columns\n    vif[\"VIF\"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n    return(vif)\n\n# Calculate VIF scores on the validation dataset.\nvif_scores = calc_vif(X_valid_KMeans)\n# \"RuntimeWarning: divide by zero\" can be safely ignored as this is caused by perfectly correlated dummy variables - see link below.","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:55:54.578569Z","iopub.execute_input":"2022-07-08T17:55:54.579122Z","iopub.status.idle":"2022-07-08T17:56:46.336494Z","shell.execute_reply.started":"2022-07-08T17:55:54.579065Z","shell.execute_reply":"2022-07-08T17:56:46.335155Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot the VIF scores obtained from the validation dataset. \nvif_scores.plot()\n\n# Display all columns with VIF scores > 10.\nvif_scores.loc[vif_scores['VIF'] > 10]","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:46.342924Z","iopub.execute_input":"2022-07-08T17:56:46.347141Z","iopub.status.idle":"2022-07-08T17:56:46.566439Z","shell.execute_reply.started":"2022-07-08T17:56:46.347065Z","shell.execute_reply":"2022-07-08T17:56:46.564456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The features listed above have very high VIF scores, which indicate a high level of multicollinearity. However, in the edge cases where some of these values tend towards infinity, these can be discounted as they represent dummy variables that are perfectly anti-correlated (e.g. one applicant/row in the dataset can only belong to a single K-Means cluster).\n> \n> More information on when it is safe to ignore multicollinearity can be found here: [Multicollinearity - Statistical Horizons](https://statisticalhorizons.com/multicollinearity/).","metadata":{}},{"cell_type":"markdown","source":"### Method 3: Principal Component Analysis\n\n> Next, we will perform Principal Component Analysis in order to understand the most significant sources of variation within our dataset.\n>\n> This will help us to understand whether there are particularly useful features in our dataset that *should* be preserved, during feature selection.","metadata":{}},{"cell_type":"code","source":"# Import/initialise key modules that will be used for visualising the Principal Component Analysis.\nfrom IPython.display import display\n\nplt.style.use(\"seaborn-whitegrid\")\nplt.rc(\"figure\", autolayout=True)\nplt.rc(\n    \"axes\",\n    labelweight=\"bold\",\n    labelsize=\"large\",\n    titleweight=\"bold\",\n    titlesize=14,\n    titlepad=10,\n)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:46.567945Z","iopub.execute_input":"2022-07-08T17:56:46.569028Z","iopub.status.idle":"2022-07-08T17:56:46.576336Z","shell.execute_reply.started":"2022-07-08T17:56:46.568982Z","shell.execute_reply":"2022-07-08T17:56:46.575031Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Define a custom function that plots the explained/cumulative variances for each Principal Component (PC) of a given dataset.\n\ndef plot_variance(pca, width=8, dpi=100):\n    fig, axs = plt.subplots(1, 2)\n    n = pca.n_components_\n    grid = np.arange(1, n + 1)\n    \n    evr = pca.explained_variance_ratio_\n    axs[0].bar(grid, evr)\n    axs[0].set(xlabel=\"Component\",\n               title=\"% Explained Variance\",\n               ylim=(0.0, 0.2))\n    \n    cv = np.cumsum(evr)\n    axs[1].plot(np.r_[0, grid], np.r_[0, cv], \"o-\")\n    axs[1].set(xlabel=\"Component\",\n               title=\"% Cumulative Variance\",\n               ylim=(0.0, 1.0))\n    \n    fig.set(figwidth=8, dpi=100)\n    return axs","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:46.578229Z","iopub.execute_input":"2022-07-08T17:56:46.578558Z","iopub.status.idle":"2022-07-08T17:56:46.594186Z","shell.execute_reply.started":"2022-07-08T17:56:46.578527Z","shell.execute_reply":"2022-07-08T17:56:46.593006Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.decomposition import PCA\n\n# Initialise the Principal Component Analysis (PCA) algorithm.\npca = PCA()\n\n# Fit the PCA algorithm to the validation dataset, and generate its corresponding PCs.\nX_valid_pca = pca.fit_transform(X_valid_KMeans)\n\n# Create a list of labels for each PC, equal in length to the number of columns in the validation dataset.\nX_valid_component_names = [f\"PC{i+1}\" for i in range(X_valid_pca.shape[1])]\n\n# Create a dataframe that contains the PCs generated, along with their respective labels.  \nX_valid_pca = pd.DataFrame(X_valid_pca, columns=X_valid_component_names)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:46.595781Z","iopub.execute_input":"2022-07-08T17:56:46.596105Z","iopub.status.idle":"2022-07-08T17:56:46.760891Z","shell.execute_reply.started":"2022-07-08T17:56:46.596078Z","shell.execute_reply":"2022-07-08T17:56:46.759523Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Use the custom function to plot the explained/cumulative variances, for each PC, within the validation dataset.\nplot_variance(pca)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:46.777900Z","iopub.execute_input":"2022-07-08T17:56:46.778995Z","iopub.status.idle":"2022-07-08T17:56:47.398055Z","shell.execute_reply.started":"2022-07-08T17:56:46.778939Z","shell.execute_reply":"2022-07-08T17:56:47.396771Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate the cumulative sum of the explained variation ratios for the first 40 PCs in the validation dataset.\npca.explained_variance_ratio_[:40].cumsum()","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:47.399663Z","iopub.execute_input":"2022-07-08T17:56:47.400184Z","iopub.status.idle":"2022-07-08T17:56:47.409413Z","shell.execute_reply.started":"2022-07-08T17:56:47.400139Z","shell.execute_reply":"2022-07-08T17:56:47.408045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> As shown in the cells above, the first 40 PCs contain just over 80% of the cumulative variance in the validation dataset.\n>\n> This means that we can still capture a significant majority of the dataset's cumulative variance, were we to use a lower dimensionality feature-space instead, rather than simply using all features together.","metadata":{}},{"cell_type":"code","source":"# Create a dataframe which displays each principal component's loading/s on each original feature in the dataset.\n\nX_valid_loadings = pd.DataFrame(\n    pca.components_.T,  # We need to transpose the matrix of loadings.\n    columns=X_valid_component_names,  # Columns are set as the principal components.\n    index=X_valid_KMeans.columns,  # Rows are set as the original features.\n)\nX_valid_loadings","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:47.411396Z","iopub.execute_input":"2022-07-08T17:56:47.411867Z","iopub.status.idle":"2022-07-08T17:56:47.447394Z","shell.execute_reply.started":"2022-07-08T17:56:47.411823Z","shell.execute_reply":"2022-07-08T17:56:47.446274Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> We can also display how each principal component is comprised, in terms of the original features' loadings.\n>\n> Hence, when looking in the top 40 principal components, we can see which of the original features contribute most strongly towards the cumulative variance of the validation dataset.","metadata":{}},{"cell_type":"code","source":"X_valid_loadings_first40 = X_valid_loadings.iloc[:, :40]\nX_valid_loadings_first40","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:47.449028Z","iopub.execute_input":"2022-07-08T17:56:47.449677Z","iopub.status.idle":"2022-07-08T17:56:47.483624Z","shell.execute_reply.started":"2022-07-08T17:56:47.449631Z","shell.execute_reply":"2022-07-08T17:56:47.482556Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"useful_cols = []\n\n# Visualise which columns in the top 40 PCs contain notable variance (i.e. where the absolute value of the columns's PC loading is > 0.25).\nfor col in X_valid_loadings_first40.columns:\n    cols = X_valid_loadings_first40[col].loc[abs(X_valid_loadings_first40[col]) > 0.25]\n    cols_df = pd.DataFrame(cols)\n    useful_cols.append(cols_df)\n    \nuseful_cols","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:47.484868Z","iopub.execute_input":"2022-07-08T17:56:47.485164Z","iopub.status.idle":"2022-07-08T17:56:47.565045Z","shell.execute_reply.started":"2022-07-08T17:56:47.485138Z","shell.execute_reply":"2022-07-08T17:56:47.563897Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The cell above indicates which columns/features appear most commonly throughout the top 40 PCs (which contain over 80% of the cumulative variance in the dataset).\n>\n> The more frequently a column name appears above, the more indicative this is that the column is useful for capturing significant variance within the dataset. It is thus more important for the models to place greater importance on these features later on.","metadata":{}},{"cell_type":"markdown","source":"### Method 4: Lasso (L1) Regularisation\n\n> Finally, we will perform feature selection via L1 (lasso) regularisation using functions that are readily available in `sklearn` - for more information, see this link: [Feature Selection via scikit-learn](https://scikit-learn.org/stable/modules/feature_selection.html#l1-based-feature-selection).","metadata":{}},{"cell_type":"code","source":"# Import key modules in order to perform LASSO-based feature selection.\nfrom sklearn.svm import LinearSVC\nfrom sklearn.feature_selection import SelectFromModel\n\n# Establish the Lasso (L1) Regularisation model that will perform feature selection.\nlinearsvc = LinearSVC(penalty=\"l1\", dual=False, tol=1e-3, C=1e-2, random_state=0).fit(X_valid_KMeans, y_valid)\nmodel = SelectFromModel(linearsvc, prefit=True)\n\n# Reduce the dataset to the most important features, using the regularisation model above.\nX_valid_L1 = model.transform(X_valid_KMeans)\n\n# Convert the transformed dataset into a dataframe with the same size/shape as the original dataset.\n# For features that were previously removed, this dataset will now include zeroes instead of their original values.\nselected_features = pd.DataFrame(model.inverse_transform(X_valid_L1),\n                                 index=X_valid_KMeans.index,\n                                 columns=X_valid_KMeans.columns)\n\n# Drop columns from the dataframe where features are deemed unimportant in capturing the dataset's variance.\n# To achieve this, we selectively drop columns where their variance is equal to 0 (i.e. where a column only contains zeroes).\nX_valid_L1reg = selected_features.drop(selected_features.columns[selected_features.var() == 0], axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:47.566700Z","iopub.execute_input":"2022-07-08T17:56:47.567350Z","iopub.status.idle":"2022-07-08T17:56:50.558173Z","shell.execute_reply.started":"2022-07-08T17:56:47.567313Z","shell.execute_reply":"2022-07-08T17:56:50.556938Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Determine the number of columns that are kept after L1 regularisation.\n\nlen(selected_features.columns[selected_features.var() != 0])","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.559345Z","iopub.execute_input":"2022-07-08T17:56:50.559650Z","iopub.status.idle":"2022-07-08T17:56:50.592883Z","shell.execute_reply.started":"2022-07-08T17:56:50.559623Z","shell.execute_reply":"2022-07-08T17:56:50.591853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> We have now reduced our dataset down to 57 features, from a starting value of 126.\n>\n> Next, we will transform our training/test subsets so that they are now constrained to the same set of columns as derived above.","metadata":{}},{"cell_type":"code","source":"# Produce a separate list containing the columns preserved after L1 regularisation.\nselected_columns = selected_features.columns[selected_features.var() != 0]\n\n# Reduce the training/test datasets to the same set of columns.\nX_train_L1reg = X_train_KMeans[selected_columns]\nX_test_L1reg = X_test_KMeans[selected_columns]","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.594269Z","iopub.execute_input":"2022-07-08T17:56:50.594671Z","iopub.status.idle":"2022-07-08T17:56:50.633264Z","shell.execute_reply.started":"2022-07-08T17:56:50.594640Z","shell.execute_reply":"2022-07-08T17:56:50.632132Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Step 8: Define the classifiers/models used\n\n> In this project, we aim to predict the target (`Response`) as well as the likelihood of belonging to the predicted group using the following approaches:\n> \n> 1. Logistic/Softmax Regression\n> \n>> **Softmax Regression** is a generalised form of the standard logistic regression model. Given an instance $X$, input features are supplied to the model in order to compute a score for each class - then, the probability of $X$ belonging to each class is estimated by applying the softmax function (which calculates the exponential of every score, then normalises them). Then, the algorithm returns a prediction as the class with the highest estimated probability (simply the class with the highest score).\n> \n> 2. Gaussian Naive Bayes\n> \n>> **Gaussian Naive Bayes** works as an extension of Naive Bayes (which is termed as \"naive\" due to its main assumption that all features are independent), by instead assuming that each feature within the data can be represented by a normal/Gaussian distribution. As a result, we only need to calculate the mean/standard deviation of each Gaussian distribution in order to predict the likelihood of an instance belonging to a given class, reducing computational complexity as compared to using other Bayesian methods.\n> \n> 3. Support Vector Machines\n> \n>> **Support Vector Machines** are capable of performing linear or nonlinear classification, and are particularly well-suited for classification of complex small- and medium-sized datasets/segmentation of high dimensionality feature spaces. They essentially work by fitting the widest possible margin between each of the classes, as a function of each of the dataset's features. As can be imagined, this technique is sensitive to feature scaling, which is why min-max scaling was performed earlier above.\n> \n> 4. Decision Trees/Random Forests\n> \n>> **Decision Trees** are powerful and versatile algorithms that are capable of fitting to/classifying complex datasets. They are generated by splitting a training dataset into two subsets recursively, based on a single feature and a corresponding threshold value (which are chosen as the pair that produces subsets with the lowest possible Gini impurity).\n>>\n>> **Random Forest Classification** works by training multiple decision trees, based on the random sampling (with replacement) of a training dataset. Input features from an unseen dataset can then be supplied to each trained decision tree in order to generate a prediction, which is subsequently averaged across all predictions to produce a final classification output; averaging across all predictions has the benefit of reducing overfitting to any given random sample within the training set.\n>\n> 5. Gradient Boosting Classifiers\n> \n>> **Boosting** refers to any ensemble method that combines several weak \"learners\" into a strong learner. Usually, this is accomplished by training predictors sequentially, where each tries to correct its predecessor.\n>>\n>> With **AdaBoost classifiers**, a base classifier (such as a Decision Tree) is trained and used to make predictions on the training set. The algorithm then increases the relative weight of misclassified training instances, and then trains a second classifier using the updated weights, before making a new set of predictions. This loop is usually repeated up until there is no significant improvement in the last round's performance, or until the maximum number of iterations has been reached.\n>>\n>> In the case of **Gradient Boosting classifiers**, the algorithm works very similarly to that used for training AdaBoost classifiers, which is by sequentially adding predictors to an ensemble - except for that, instead of updating instance weights at every iteration, the next predictor is instead fitted to the residual errors made by the previous predictor. Over the course of several iterations, this can lead to a highly refined model which is capable of recognising complex nuances in the training dataset and providing accurate classifications, however care must also be taken to avoid overfitting.","metadata":{}},{"cell_type":"code","source":"# Import the classification models from sklearn/xgboost.\n\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.naive_bayes import GaussianNB\nfrom sklearn.svm import SVC\nfrom sklearn.tree import DecisionTreeClassifier\nfrom sklearn.ensemble import RandomForestClassifier, AdaBoostClassifier, GradientBoostingClassifier\nfrom xgboost import XGBClassifier","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.635784Z","iopub.execute_input":"2022-07-08T17:56:50.636530Z","iopub.status.idle":"2022-07-08T17:56:50.819193Z","shell.execute_reply.started":"2022-07-08T17:56:50.636497Z","shell.execute_reply":"2022-07-08T17:56:50.818016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Step 9: Perform hyperparameter optimisation via grid-search methods (using validation data)\n\n> Here, we will utilise exhaustive grid-search methods in order to optimise several of the models' hyperparameters. These are:\n>\n>> `LogisticRegression`\n>> * `tol` represents the tolerance for the stopping criteria, i.e. the minimum step size of the loss function's improvement.\n>> * `C` represents the inverse of the regularisation strength, thus determining how strongly the model fits to the data.\n>>\n>>`GaussianNB`\n>> * `var_smoothing` represents the degree to which we widen/\"smooth\" our Gaussian distributions, to account for additional samples that are further away from the distribution mean.\n>> \n>>`SVC(kernel=\"linear\")`\n>> * `C` represents the inverse of the regularisation strength, thus determining how strongly the model fits to the data.\n>> * `tol` represents the tolerance for the stopping criteria, i.e. the minimum step size of the loss function's improvement.\n>> \n>>`LinearSVC`\n>> * `tol` represents the tolerance for the stopping criteria, i.e. the minimum step size of the loss function's improvement.\n>> * `C` represents the inverse of the regularisation strength, thus determining how strongly the model fits to the data.\n>> \n>>`SVC(kernel=\"poly\")`\n>> * `C` represents the inverse of the regularisation strength, thus determining how strongly the model fits to the data.\n>> * `degree` represents the degree of the polynomial kernel function that we want the model to fit.\n>> * `tol` represents the tolerance for the stopping criteria, i.e. the minimum step size of the loss function's improvement.\n>> \n>>`SVC(kernel=\"rbf\")`\n>> * `C` represents the inverse of the regularisation strength, thus determining how strongly the model fits to the data.\n>> * `tol` represents the tolerance for the stopping criteria, i.e. the minimum step size of the loss function's improvement.\n>> \n>>`SVC(kernel=\"sigmoid\")`\n>> * `C` represents the inverse of the regularisation strength, thus determining how strongly the model fits to the data.\n>> * `tol` represents the tolerance for the stopping criteria, i.e. the minimum step size of the loss function's improvement.\n>> \n>>`DecisionTreeClassifier`\n>> * `max_depth` represents the maximum depth of the tree, thus determining how strongly the model fits to the data.\n>> * `max_features` represents the number of features to consider, when looking for the optimum split.\n>> \n>>`RandomForestClassifier`\n>> * `n_estimators` represents the number of decision trees that are implemented by the random forest classifier.\n>> * `max_depth` represents the maximum depth of the tree, thus determining how strongly the model fits to the data.\n>> * `max_features` represents the number of features to consider, when looking for the optimum split.\n>> \n>>`AdaBoostClassifier`\n>> * `n_estimators` represents the maximum number of boosting rounds/decision trees to be used, at which point the boosting process is terminated.\n>> * `learning_rate` represents the contribution weighting applied to each decision tree, at each boosting iteration.\n>> \n>>`GradientBoostingClassifier`\n>> * `learning_rate` represents the contribution weighting applied to each decision tree, at each boosting iteration.\n>> * `n_estimators` represents the maximum number of boosting rounds/decision trees to be used, at which point the boosting process is terminated.\n>> \n>>`XGBClassifier`\n>> * `n_estimators` represents the maximum number of boosting rounds/decision trees to be used, at which point the boosting process is terminated.\n>> * `learning_rate` represents the contribution weighting applied to each decision tree, at each boosting iteration.\n>>\n>\n> More information regarding exhaustive grid-search methods can be found at the following page: [Exhaustive grid-search methods for tuning hyperparameters - scikit-learn](https://scikit-learn.org/stable/modules/grid_search.html#exhaustive-grid-search).","metadata":{}},{"cell_type":"code","source":"# Import the GridSearchCV function from sklearn.\n\nfrom sklearn.model_selection import GridSearchCV","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.820960Z","iopub.execute_input":"2022-07-08T17:56:50.821416Z","iopub.status.idle":"2022-07-08T17:56:50.827168Z","shell.execute_reply.started":"2022-07-08T17:56:50.821371Z","shell.execute_reply":"2022-07-08T17:56:50.826079Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Initialise each model\n\n> We set up each of the classifiers as baseline models, by initialising each model as a new object.","metadata":{}},{"cell_type":"code","source":"## Initialise each classifier - these are capable of providing predictions as well as their associated probabilities/likelihoods/confidences.\n\n### Logistic/softmax regressors\nModel1_Base = LogisticRegression(random_state=0,\n                                 solver='liblinear')\n\n### Naive Bayes classifiers\nModel2_Base = GaussianNB()\n\n### Support Vector Machines (linear/non-linear)\nModel3_Base = SVC(kernel='linear',\n                  probability=True,\n                  max_iter=1000,\n                  random_state=0)\n\nModel4_Base = LinearSVC(dual=False,\n                        random_state=0,\n                        max_iter=1000)\n\nModel5_Base = SVC(kernel='poly',\n                  probability=True,\n                  max_iter=1000,\n                  random_state=0)\n\nModel6_Base = SVC(kernel='rbf',\n                  probability=True,\n                  max_iter=1000,\n                  random_state=0)\n\nModel7_Base = SVC(kernel='sigmoid', \n                  probability=True,\n                  max_iter=1000,\n                  random_state=0)\n\n### Decision Trees\nModel8_Base = DecisionTreeClassifier(random_state=0)\n\n### Random Forests\nModel9_Base = RandomForestClassifier(n_jobs=-1,\n                                     random_state=0)\n\n### Gradient Boosting Machines\nModel10_Base = AdaBoostClassifier(random_state=0)\n\nModel11_Base = GradientBoostingClassifier(random_state=0)\n\nModel12_Base = XGBClassifier(random_state=0,\n                             n_jobs=-1,\n                             eval_metric=\"merror\")","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.828673Z","iopub.execute_input":"2022-07-08T17:56:50.829062Z","iopub.status.idle":"2022-07-08T17:56:50.841018Z","shell.execute_reply.started":"2022-07-08T17:56:50.829023Z","shell.execute_reply":"2022-07-08T17:56:50.839791Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Declare the hyperparameter grids for each model\n\n> Next, we declare the models' parameter grids that we require the grid-searches to be performed over.","metadata":{}},{"cell_type":"code","source":"# Config for model 1 - LogisticRegression.\nparam_grid_model1 = {'tol': [0.1, 0.01, 0.001, 0.0001, 0.00001, 0.000001],\n                     'C': [100.0, 10.0, 1.0, 0.1, 0.01, 0.001]}\n\n\n# Config for model 2 - GaussianNB.\nparam_grid_model2 = {'var_smoothing': [1e-01, 1e-02, 1e-03, 1e-04, 1e-05, 1e-06, 1e-07, 1e-08, 1e-09, 1e-10]}\n\n\n# Config for model 3 - SVC (linear kernel).\nparam_grid_model3 = {'C': [100, 10, 1.0, 0.1, 0.01, 0.001],\n                     'tol': [1e-01, 1e-02, 1e-03, 1e-04, 1e-05, 1e-06]}\n\n\n# Config for model 4 - LinearSVC.\nparam_grid_model4 = {'tol': [1e-01, 1e-02, 1e-03, 1e-04, 1e-05, 1e-06],\n                     'C': [100, 10, 1.0, 0.1, 0.01, 0.001]}\n\n\n# Config for model 5 - SVC (polynomial kernel).\nparam_grid_model5 = {'C': [100, 10, 1.0, 0.1, 0.01, 0.001],\n                     'degree': [0, 1, 2, 3, 4, 5, 6],\n                     'tol': [1e-01, 1e-02, 1e-03, 1e-04, 1e-05, 1e-06]}\n\n\n# Config for model 6 - SVC (rbf kernel).\nparam_grid_model6 = {'C': [100, 10, 1.0, 0.1, 0.01, 0.001],\n                     'tol': [1e-01, 1e-02, 1e-03, 1e-04, 1e-05, 1e-06]}\n\n\n# Config for model 7 - SVC (sigmoid kernel).\nparam_grid_model7 = {'C': [100, 10, 1.0, 0.1, 0.01, 0.001],\n                     'tol': [1e-01, 1e-02, 1e-03, 1e-04, 1e-05, 1e-06]}\n\n\n# Config for model 8 - DecisionTreeClassifier.\nparam_grid_model8 = {'max_depth': [5, 10, 15, 20, 25, 30, 35, 40],\n                     'max_features': [0.2, 0.4, 0.6, 0.8, 1.0]}\n\n\n# Config for model 9 - RandomForestClassifier.\nparam_grid_model9 = {'n_estimators': [10, 100, 250, 500, 1000],\n                     'max_depth': [5, 10, 15, 20, 25, 30, 35, 40],\n                     'max_features': [0.2, 0.4, 0.6, 0.8, 1.0]}\n\n\n# Config for model 10 - AdaBoostClassifier.\nparam_grid_model10 = {'n_estimators': [10, 100, 250, 500, 1000],\n                      'learning_rate': [1e-0, 1e-01, 1e-02, 1e-03, 1e-04, 1e-05]}\n\n\n# Config for model 11 - GradientBoostingClassifier.\nparam_grid_model11 = {'learning_rate': [1e-0, 1e-01, 1e-02, 1e-03, 1e-04, 1e-05],\n                      'n_estimators': [10, 100, 250, 500, 1000]}\n\n\n# Config for model 12 - XGBClassifier.\nparam_grid_model12 = {'n_estimators': [10, 100, 250, 500, 1000],\n                      'learning_rate': [1e-0, 1e-01, 1e-02, 1e-03, 1e-04, 1e-05]}\n\n# Create a list containing each of the baseline models.\nBaseModels = [Model1_Base, Model2_Base, Model3_Base, \n              Model4_Base, Model5_Base, Model6_Base,\n              Model7_Base, Model8_Base, Model9_Base,\n              Model10_Base, Model11_Base, Model12_Base]\n\n# Create a list containing each of the model's parameter-grid dictionaries.\nParamGrids = [param_grid_model1, param_grid_model2, param_grid_model3, \n              param_grid_model4, param_grid_model5, param_grid_model6,\n              param_grid_model7, param_grid_model8, param_grid_model9,\n              param_grid_model10, param_grid_model11, param_grid_model12]","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.842612Z","iopub.execute_input":"2022-07-08T17:56:50.842998Z","iopub.status.idle":"2022-07-08T17:56:50.861841Z","shell.execute_reply.started":"2022-07-08T17:56:50.842965Z","shell.execute_reply":"2022-07-08T17:56:50.860730Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Define a custom gridsearch function to time each experiment\n\n> In order to run each grid-search consecutively in a single loop, we define a custom function that performs a grid-search optimisation for a given model, and returns an object containing a set of optimised hyperparameters that can be easily retrieved afterwards.","metadata":{}},{"cell_type":"code","source":"from datetime import datetime\nimport pytz\n\n# Define a custom function that performs a grid-search for a given model and parameter grid.\ndef GridSearcher(model, param_grid):\n    ## Startup\n    timezone = pytz.timezone('Europe/London')\n    start_time = datetime.now(timezone)\n    print(\"Running GridSearchCV for:\", str(model))\n    print(\"Starting at:\",start_time.strftime(\"%H:%M:%S\"))\n    \n    ## Perform the grid-search\n    GridSearcher = GridSearchCV(model, param_grid, scoring='balanced_accuracy', refit=False, error_score='raise', verbose=2)\n    GridSearcher.fit(X_valid_L1reg, y_valid)\n    # HOLDOUT n_jobs=-1\n    \n    ## Finish\n    finish_time = datetime.now(timezone)\n    print(\"Finished at:\", finish_time.strftime(\"%H:%M:%S\"))\n    duration = finish_time - start_time\n    dur = divmod(duration.seconds, 60)\n    print(\"Duration: \", dur[0], 'minutes', dur[1], 'seconds')\n    return GridSearcher","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.864309Z","iopub.execute_input":"2022-07-08T17:56:50.864737Z","iopub.status.idle":"2022-07-08T17:56:50.877777Z","shell.execute_reply.started":"2022-07-08T17:56:50.864685Z","shell.execute_reply":"2022-07-08T17:56:50.876786Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Run grid-searches for all models\n\n> For each model, we loop through the custom grid-search function and retrieve its optimised hyperparameters.\n>\n> This section is commented out as it takes several hours to run.","metadata":{}},{"cell_type":"code","source":"##     tz = pytz.timezone('Europe/London')\n##     start = datetime.now(tz)\n##     print(\"GridSearchCV runs starting at:\",start.strftime(\"%H:%M:%S\"))\n\n##     for i in range(0, 12):\n##         GridSearchResult = GridSearcher(BaseModels[i], ParamGrids[i])\n##         print(GridSearchResult.best_params_)\n##         GridSearchResults.append(GridSearchResult)\n    \n##     finish = datetime.now(tz)\n##     print(\"GridSearchCV runs finished at:\", finish.strftime(\"%H:%M:%S\"))\n##     duration = finish - start\n##     dur = divmod(duration.seconds, 60)\n##     print(\"Total Duration: \", dur[0], 'minutes', dur[1], 'seconds')\n\n##     for i in GridSearchResults:\n##         print(i.best_params_)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.879041Z","iopub.execute_input":"2022-07-08T17:56:50.879518Z","iopub.status.idle":"2022-07-08T17:56:50.890493Z","shell.execute_reply.started":"2022-07-08T17:56:50.879489Z","shell.execute_reply":"2022-07-08T17:56:50.889307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The following sets of optimised hyperparameters were obtained via the cell above.\n\nModel1_Opt_Params = {'C': 100.0, 'tol': 0.01}\nModel2_Opt_Params = {'var_smoothing': 0.0001}\nModel3_Opt_Params = {'C': 1.0, 'tol': 0.1}\nModel4_Opt_Params = {'C': 100, 'tol': 0.001}\nModel5_Opt_Params = {'C': 0.1, 'degree': 6, 'tol': 0.1}\nModel6_Opt_Params = {'C': 10, 'tol': 0.01}\nModel7_Opt_Params = {'C': 1.0, 'tol': 0.01}\nModel8_Opt_Params = {'max_depth': 10, 'max_features': 0.8}\nModel9_Opt_Params = {'max_depth': 40, 'max_features': 0.6, 'n_estimators': 1000}\nModel10_Opt_Params = {'learning_rate': 1.0, 'n_estimators': 100}\nModel11_Opt_Params = {'learning_rate': 0.1, 'n_estimators': 250}\nModel12_Opt_Params = {'learning_rate': 1.0, 'n_estimators': 10}\n\n# Create a list containing each of the models' optimised hyperparameters.\nOptimisedParams = [Model1_Opt_Params, Model2_Opt_Params, Model3_Opt_Params,\n                   Model4_Opt_Params, Model5_Opt_Params, Model6_Opt_Params,\n                   Model7_Opt_Params, Model8_Opt_Params, Model9_Opt_Params,\n                   Model10_Opt_Params, Model11_Opt_Params, Model12_Opt_Params]","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.891593Z","iopub.execute_input":"2022-07-08T17:56:50.892233Z","iopub.status.idle":"2022-07-08T17:56:50.907999Z","shell.execute_reply.started":"2022-07-08T17:56:50.892197Z","shell.execute_reply":"2022-07-08T17:56:50.906791Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Update all baseline models to use optimised hyperparameters\n\n> Finally, we update the baseline models declared earlier above, by passing the optimised hyperparameters to each model using the `.set_params(**kwargs)` method.","metadata":{}},{"cell_type":"code","source":"OptimisedModels = []\n\n# Update each of the baseline models to include their (respective) optimised hyperparameters.\nfor i in range(0, 12):\n    OptimisedModel = BaseModels[i].set_params(**OptimisedParams[i])\n    OptimisedModels.append(OptimisedModel)\n    \nOptimisedModels","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.910151Z","iopub.execute_input":"2022-07-08T17:56:50.910706Z","iopub.status.idle":"2022-07-08T17:56:50.931276Z","shell.execute_reply.started":"2022-07-08T17:56:50.910660Z","shell.execute_reply":"2022-07-08T17:56:50.929836Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Step 10: Train each model (using training data)\n\n> In this section, we will now train our optimised models on the encoded training dataset.","metadata":{}},{"cell_type":"code","source":"# Create a duplicate copy of the OptimisedModels list, for use in model training.\nTrainedModels = OptimisedModels","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.933093Z","iopub.execute_input":"2022-07-08T17:56:50.934036Z","iopub.status.idle":"2022-07-08T17:56:50.939802Z","shell.execute_reply.started":"2022-07-08T17:56:50.933989Z","shell.execute_reply":"2022-07-08T17:56:50.938387Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> XGBoost's classifiers require labels to be supplied in a zero-based fashion (i.e. classes 1-8 must instead be converted into labels 0-7 prior to training/prediction). Previously, this would have been handled automatically via the use of the `use_label_encoder` keyword argument - however, we must now explicitly label-encode our dataset prior to fitting/predicting as this has since become deprecated.\n>\n> Hence, we will establish a label encoder that converts the classes into compatible labels.","metadata":{}},{"cell_type":"code","source":"# Label-encode each dataset for compatibility with the XGBoost classifier (model 12) - \"use_label_encoder\" has since become deprecated.\nfrom sklearn.preprocessing import LabelEncoder\nle = LabelEncoder()\n\n# Convert each set of class labels (1-8) into encoded labels (0-7).\ny_train = le.fit_transform(y_train)\ny_valid = le.transform(y_valid)\ny_test = le.transform(y_test)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.941657Z","iopub.execute_input":"2022-07-08T17:56:50.942835Z","iopub.status.idle":"2022-07-08T17:56:50.954971Z","shell.execute_reply.started":"2022-07-08T17:56:50.942797Z","shell.execute_reply":"2022-07-08T17:56:50.953705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> We can now train all of our models using the encoded/transformed labels.","metadata":{}},{"cell_type":"code","source":"# Train each of the optimised models, using the training dataset.\nfor model in TrainedModels:\n    tz = pytz.timezone('Europe/London')\n    start = datetime.now(tz)\n    print(\"Model training started at:\",start.strftime(\"%H:%M:%S\"))\n    model.fit(X_train_L1reg, y_train)\n    finish = datetime.now(tz)\n    print(\"Model training finished at:\",finish.strftime(\"%H:%M:%S\"))\n    duration = finish - start\n    dur = divmod(duration.seconds, 60)\n    print(\"Total Duration for model training: \", dur[0], 'minutes', dur[1], 'seconds')\n    print(\"Model \"+str(TrainedModels.index(model)+1)+\" - \"+model.__class__.__name__+\" has been trained.\")","metadata":{"execution":{"iopub.status.busy":"2022-07-08T17:56:50.956399Z","iopub.execute_input":"2022-07-08T17:56:50.956780Z","iopub.status.idle":"2022-07-08T18:19:33.917282Z","shell.execute_reply.started":"2022-07-08T17:56:50.956741Z","shell.execute_reply":"2022-07-08T18:19:33.916316Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Step 11: Perform probability calibration on each model (using validation data)\n\n> Our models are now capable of providing predictions as well as some initial probability estimates - however, we wish to obtain true probabilities as currently, by default, ours do not yet factor in the expected distributions of the predicted classes. To this end, we will need to perform probability calibration for all of our models, so that they are capable of providing true likelihoods of belonging to each respective class.\n>\n> More information on this topic can be found at the following page: [Probability Calibration - scikit-learn](https://scikit-learn.org/stable/modules/calibration.html).","metadata":{}},{"cell_type":"code","source":"from sklearn.calibration import CalibratedClassifierCV","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:19:33.918739Z","iopub.execute_input":"2022-07-08T18:19:33.919541Z","iopub.status.idle":"2022-07-08T18:19:33.928032Z","shell.execute_reply.started":"2022-07-08T18:19:33.919502Z","shell.execute_reply":"2022-07-08T18:19:33.926759Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calibrate each classifier using the validation dataset.\n\nModel1_Calibrated = CalibratedClassifierCV(base_estimator=TrainedModels[0], method=\"isotonic\", cv=\"prefit\").fit(X_valid_L1reg, y_valid)\nModel2_Calibrated = CalibratedClassifierCV(base_estimator=TrainedModels[1], method=\"isotonic\", cv=\"prefit\").fit(X_valid_L1reg, y_valid)\nModel3_Calibrated = CalibratedClassifierCV(base_estimator=TrainedModels[2], method=\"isotonic\", cv=\"prefit\").fit(X_valid_L1reg, y_valid)\nModel4_Calibrated = CalibratedClassifierCV(base_estimator=TrainedModels[3], method=\"isotonic\", cv=\"prefit\").fit(X_valid_L1reg, y_valid)\nModel5_Calibrated = CalibratedClassifierCV(base_estimator=TrainedModels[4], method=\"isotonic\", cv=\"prefit\").fit(X_valid_L1reg, y_valid)\nModel6_Calibrated = CalibratedClassifierCV(base_estimator=TrainedModels[5], method=\"isotonic\", cv=\"prefit\").fit(X_valid_L1reg, y_valid)\nModel7_Calibrated = CalibratedClassifierCV(base_estimator=TrainedModels[6], method=\"isotonic\", cv=\"prefit\").fit(X_valid_L1reg, y_valid)\nModel8_Calibrated = CalibratedClassifierCV(base_estimator=TrainedModels[7], method=\"isotonic\", cv=\"prefit\").fit(X_valid_L1reg, y_valid)\nModel9_Calibrated = CalibratedClassifierCV(base_estimator=TrainedModels[8], method=\"isotonic\", cv=\"prefit\").fit(X_valid_L1reg, y_valid)\nModel10_Calibrated = CalibratedClassifierCV(base_estimator=TrainedModels[9], method=\"isotonic\", cv=\"prefit\").fit(X_valid_L1reg, y_valid)\nModel11_Calibrated = CalibratedClassifierCV(base_estimator=TrainedModels[10], method=\"isotonic\", cv=\"prefit\").fit(X_valid_L1reg, y_valid)\nModel12_Calibrated = CalibratedClassifierCV(base_estimator=TrainedModels[11], method=\"isotonic\", cv=\"prefit\").fit(X_valid_L1reg, y_valid)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:19:33.930126Z","iopub.execute_input":"2022-07-08T18:19:33.930919Z","iopub.status.idle":"2022-07-08T18:21:15.088512Z","shell.execute_reply.started":"2022-07-08T18:19:33.930859Z","shell.execute_reply":"2022-07-08T18:21:15.087509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Use each model to generate predictions (as probability estimates of belonging to each Response class) from the validation dataset.\n\nModel1_Valid_PredProb = Model1_Calibrated.predict_proba(X_valid_L1reg)\nModel2_Valid_PredProb = Model2_Calibrated.predict_proba(X_valid_L1reg)\nModel3_Valid_PredProb = Model3_Calibrated.predict_proba(X_valid_L1reg)\nModel4_Valid_PredProb = Model4_Calibrated.predict_proba(X_valid_L1reg)\nModel5_Valid_PredProb = Model5_Calibrated.predict_proba(X_valid_L1reg)\nModel6_Valid_PredProb = Model6_Calibrated.predict_proba(X_valid_L1reg)\nModel7_Valid_PredProb = Model7_Calibrated.predict_proba(X_valid_L1reg)\nModel8_Valid_PredProb = Model8_Calibrated.predict_proba(X_valid_L1reg)\nModel9_Valid_PredProb = Model9_Calibrated.predict_proba(X_valid_L1reg)\nModel10_Valid_PredProb = Model10_Calibrated.predict_proba(X_valid_L1reg)\nModel11_Valid_PredProb = Model11_Calibrated.predict_proba(X_valid_L1reg)\nModel12_Valid_PredProb = Model12_Calibrated.predict_proba(X_valid_L1reg)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:21:15.090193Z","iopub.execute_input":"2022-07-08T18:21:15.090694Z","iopub.status.idle":"2022-07-08T18:22:58.861969Z","shell.execute_reply.started":"2022-07-08T18:21:15.090647Z","shell.execute_reply":"2022-07-08T18:22:58.860594Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Convert the probability estimates into label-based predictions, then label-encode for compatibility with XGBoost.\n\n# Model 1 - converting probabilities into label predictions\nModel1_Valid_Preds = np.argmax(Model1_Valid_PredProb, axis=1)\nModel1_Valid_Classes = Model1_Calibrated.classes_\nModel1_Valid_Preds = [Model1_Valid_Classes[i] for i in Model1_Valid_Preds]\nModel1_Valid_Preds = le.inverse_transform(Model1_Valid_Preds)\n\n# Model 2 - converting probabilities into label predictions\nModel2_Valid_Preds = np.argmax(Model2_Valid_PredProb, axis=1)\nModel2_Valid_Classes = Model2_Calibrated.classes_\nModel2_Valid_Preds = [Model2_Valid_Classes[i] for i in Model2_Valid_Preds]\nModel2_Valid_Preds = le.inverse_transform(Model2_Valid_Preds)\n\n# Model 3 - converting probabilities into label predictions\nModel3_Valid_Preds = np.argmax(Model3_Valid_PredProb, axis=1)\nModel3_Valid_Classes = Model3_Calibrated.classes_\nModel3_Valid_Preds = [Model3_Valid_Classes[i] for i in Model3_Valid_Preds]\nModel3_Valid_Preds = le.inverse_transform(Model3_Valid_Preds)\n\n# Model 4 - converting probabilities into label predictions\nModel4_Valid_Preds = np.argmax(Model4_Valid_PredProb, axis=1)\nModel4_Valid_Classes = Model4_Calibrated.classes_\nModel4_Valid_Preds = [Model4_Valid_Classes[i] for i in Model4_Valid_Preds]\nModel4_Valid_Preds = le.inverse_transform(Model4_Valid_Preds)\n\n# Model 5 - converting probabilities into label predictions\nModel5_Valid_Preds = np.argmax(Model5_Valid_PredProb, axis=1)\nModel5_Valid_Classes = Model5_Calibrated.classes_\nModel5_Valid_Preds = [Model5_Valid_Classes[i] for i in Model5_Valid_Preds]\nModel5_Valid_Preds = le.inverse_transform(Model5_Valid_Preds)\n\n# Model 6 - converting probabilities into label predictions\nModel6_Valid_Preds = np.argmax(Model6_Valid_PredProb, axis=1)\nModel6_Valid_Classes = Model6_Calibrated.classes_\nModel6_Valid_Preds = [Model6_Valid_Classes[i] for i in Model6_Valid_Preds]\nModel6_Valid_Preds = le.inverse_transform(Model6_Valid_Preds)\n\n# Model 7 - converting probabilities into label predictions\nModel7_Valid_Preds = np.argmax(Model7_Valid_PredProb, axis=1)\nModel7_Valid_Classes = Model7_Calibrated.classes_\nModel7_Valid_Preds = [Model7_Valid_Classes[i] for i in Model7_Valid_Preds]\nModel7_Valid_Preds = le.inverse_transform(Model7_Valid_Preds)\n\n# Model 8 - converting probabilities into label predictions\nModel8_Valid_Preds = np.argmax(Model8_Valid_PredProb, axis=1)\nModel8_Valid_Classes = Model8_Calibrated.classes_\nModel8_Valid_Preds = [Model8_Valid_Classes[i] for i in Model8_Valid_Preds]\nModel8_Valid_Preds = le.inverse_transform(Model8_Valid_Preds)\n\n# Model 9 - converting probabilities into label predictions\nModel9_Valid_Preds = np.argmax(Model9_Valid_PredProb, axis=1)\nModel9_Valid_Classes = Model9_Calibrated.classes_\nModel9_Valid_Preds = [Model9_Valid_Classes[i] for i in Model9_Valid_Preds]\nModel9_Valid_Preds = le.inverse_transform(Model9_Valid_Preds)\n\n# Model 10 - converting probabilities into label predictions\nModel10_Valid_Preds = np.argmax(Model10_Valid_PredProb, axis=1)\nModel10_Valid_Classes = Model10_Calibrated.classes_\nModel10_Valid_Preds = [Model10_Valid_Classes[i] for i in Model10_Valid_Preds]\nModel10_Valid_Preds = le.inverse_transform(Model10_Valid_Preds)\n\n# Model 11 - converting probabilities into label predictions\nModel11_Valid_Preds = np.argmax(Model11_Valid_PredProb, axis=1)\nModel11_Valid_Classes = Model11_Calibrated.classes_\nModel11_Valid_Preds = [Model11_Valid_Classes[i] for i in Model11_Valid_Preds]\nModel11_Valid_Preds = le.inverse_transform(Model11_Valid_Preds)\n\n# Model 12 - converting probabilities into label predictions\nModel12_Valid_Preds = np.argmax(Model12_Valid_PredProb, axis=1)\nModel12_Valid_Classes = Model12_Calibrated.classes_\nModel12_Valid_Preds = [Model12_Valid_Classes[i] for i in Model12_Valid_Preds]\nModel12_Valid_Preds = le.inverse_transform(Model12_Valid_Preds)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:22:58.865376Z","iopub.execute_input":"2022-07-08T18:22:58.865834Z","iopub.status.idle":"2022-07-08T18:22:58.974825Z","shell.execute_reply.started":"2022-07-08T18:22:58.865794Z","shell.execute_reply":"2022-07-08T18:22:58.973229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Step 12: Evaluate the performance of each model (using validation data)\n\n> In a real-world situation, other factors - such as training time, model explainability and costs/difficulties associated with model deployment/maintenance - will also have an impact on which model should be chosen for live use.\n>\n> However, for the purposes of this project we are simply aiming to compare each of the machine learning algorithms used thus far, solely based in terms of their predictive performances, in order to select the \"best\" model. In the code below, a number of techniques are demonstrated in order to display how each model has performed in terms of classifying applicants into each Response group.\n>\n> We will use a combination of functions across both `sklearn` as well as `scikitplot` in order to visualise each model's performance.","metadata":{}},{"cell_type":"code","source":"# Import key modules/functions for evaluating model performance.\nfrom scikitplot.metrics import plot_roc, plot_confusion_matrix\nfrom sklearn.metrics import classification_report, balanced_accuracy_score","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:22:58.976713Z","iopub.execute_input":"2022-07-08T18:22:58.977542Z","iopub.status.idle":"2022-07-08T18:22:59.006863Z","shell.execute_reply.started":"2022-07-08T18:22:58.977479Z","shell.execute_reply":"2022-07-08T18:22:59.004398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Convert the validation dataset's encoded labels (0-7) back to the original set of classes (1-8), for clearer reviews of model performance.\ny_valid = le.inverse_transform(y_valid)\n\n# Create a list containing each model's predicted labels.\nValid_Preds = [Model1_Valid_Preds, Model2_Valid_Preds, Model3_Valid_Preds,\n               Model4_Valid_Preds, Model5_Valid_Preds, Model6_Valid_Preds,\n               Model7_Valid_Preds, Model8_Valid_Preds, Model9_Valid_Preds,\n               Model10_Valid_Preds, Model11_Valid_Preds, Model12_Valid_Preds]\n\n# Create a list containing each model's prediction probabilities.\nValid_PredProbs = [Model1_Valid_PredProb, Model2_Valid_PredProb, Model3_Valid_PredProb,\n                   Model4_Valid_PredProb, Model5_Valid_PredProb, Model6_Valid_PredProb,\n                   Model7_Valid_PredProb, Model8_Valid_PredProb, Model9_Valid_PredProb,\n                   Model10_Valid_PredProb, Model11_Valid_PredProb, Model12_Valid_PredProb]","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:22:59.009165Z","iopub.execute_input":"2022-07-08T18:22:59.009797Z","iopub.status.idle":"2022-07-08T18:22:59.023159Z","shell.execute_reply.started":"2022-07-08T18:22:59.009681Z","shell.execute_reply":"2022-07-08T18:22:59.021161Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### ROC Curves\n\n> The plots below display the ROC Curves (i.e. TP rate vs. FP rate) for each of the 12 models.","metadata":{}},{"cell_type":"code","source":"# For each model, plot an Receiver Operating Characteristic (ROC) curve to display model performance - TP rate vs FP rate.\n\nfor y_valid_predprobs in Valid_PredProbs:\n    plot_roc(y_valid, y_valid_predprobs)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:22:59.024416Z","iopub.execute_input":"2022-07-08T18:22:59.024990Z","iopub.status.idle":"2022-07-08T18:23:04.603251Z","shell.execute_reply.started":"2022-07-08T18:22:59.024935Z","shell.execute_reply":"2022-07-08T18:23:04.601893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> As shown in the charts above, each model has a wide range of AUC values against each class, indicating a varying degree of sensitivity/specificity across the dataset. The two highest average AUC values were achieved by the gradient boosting classifiers (model **11**: macro-average AUC=0.84, and model **12**: macro-average AUC=0.83).\n>\n> It is also worth noting that both of these models performed best when predicting applicants for classes 3/4/8, as indicated by the high AUC scores (~0.9) for their respective ROC curves, whereas they perform slightly worse when classifying applicants into the remaining groups.","metadata":{}},{"cell_type":"markdown","source":"### Classification Reports (Precisions/Recalls/F1-Scores)","metadata":{}},{"cell_type":"code","source":"# For each model, print a Classification Report to display model performance - Precision/Recall/F1-score/Accuracy.\n\nfor y_valid_preds in Valid_Preds:\n    print(classification_report(y_valid, y_valid_preds))","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:04.605338Z","iopub.execute_input":"2022-07-08T18:23:04.606121Z","iopub.status.idle":"2022-07-08T18:23:04.976259Z","shell.execute_reply.started":"2022-07-08T18:23:04.606059Z","shell.execute_reply":"2022-07-08T18:23:04.974541Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> As shown in the summary tables above, each model has a wide range of precision/recall values against each class, indicating a varying degree of accuracy across the dataset.\n>\n> In terms of precision, each of the models display somewhat mediocre performance for each class, with the highest value consistently belonging to class 8. Each model also consistently shows very high recall for class 8, and poorer values against the other groups. The same trend can also be observed for F1-score, indicating that the models are strongly fitted towards predicting applicants in class 8 as we would expect (given that they represent the largest class in our dataset).\n> \n> It is interesting to observe that, for class 3, model **11**'s precision/recall/F1-score values are equal to 0, which means that it did not correctly predict any applicants for this risk rating. On the other hand, model **12** showed a precision of 0.27, a recall of 0.05 and an F1-score of 0.08. However, as will be discussed further in the next section, this does not necessarily mean that model 11 is wholly inaccurate and should not be trusted altogether.\n>\n> Finally, it is important not to judge the models too severely based on their precision/recall/F1-scores alone. These are typically harsh metrics that can indicate where models are strongly under/overfitting as they evaluate accuracy in a \"one vs. rest\" fashion - if the set of possible `Response` values was much smaller instead, e.g. by grouping together classes 1-8 into Low/Med/High-risk, then each model's performance would appear to significantly improve.","metadata":{}},{"cell_type":"markdown","source":"### Confusion Matrices","metadata":{}},{"cell_type":"code","source":"# For each model, plot a Confusion Matrix to display model performance - True labels vs False labels.\n\nfor y_valid_preds in Valid_Preds:\n    plot_confusion_matrix(y_valid, y_valid_preds)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:04.977819Z","iopub.execute_input":"2022-07-08T18:23:04.978531Z","iopub.status.idle":"2022-07-08T18:23:12.016589Z","shell.execute_reply.started":"2022-07-08T18:23:04.978491Z","shell.execute_reply":"2022-07-08T18:23:12.015106Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The confusion matrices above have been assembled in order to visualise how each model's predictions map to the true set of labels available in the test dataset. These plots help to illustrate as to whether the model has become overfitted to any particular classes in the dataset, or whether they remain sufficiently generalised and capable of sensibly classifying applicants.\n>\n> From the plots above, most appear to be fairly well-fitted to the data, in that they are capable of replicating the distribution of `Response` values (shown earlier above in Step 1) reasonably well. We can also see that a majority of the missed cases are in fact relatively close to where they should be (e.g. predicted as class 8, was actually class 7). However, there are some notably underperforming models which appear to be highly overfitted towards predicting applicants as belonging to class 8 (which represents the largest proportion across all applicants in the dataset). The confusion matrix for model **7** (SVC with sigmoid kernel) demonstrates this trait extremely well, which shows that only a tiny proportion of predictions were made in any other class apart from 8. Models **3** and **5** (SVCs with linear and polynomial kernels, respectively) also demonstrate this trend to a lesser extent as well.\n> \n> Models **11** and **12**, on the other hand, appear to generalise reasonably well to the dataset, mimicking the distribution plot shown in Step 1. Furthermore, the only notable drawback of these models that can be observed is that they tend to predict some low-risk applicants (e.g. classes 1-2) as high-risk applicants (e.g. classes 6-8) with a higher than expected frequency - see the top-right corners of each plot. In terms of using either of these models in a real-world business scenario, this would only mean that more effort/time is potentially wasted on scrutinising low-risk applicants further before offering a policy, rather than treating high-risk applicants with a light touch and introducing unnecessary risk into the insurer's portfolio.","metadata":{}},{"cell_type":"markdown","source":"### Balanced Accuracy Scores\n\n> Here, we calculate the balanced accuracy score of each model, by comparing their predictions against the known set of results in the validation dataset. This can be derived as the arithmetic mean of the recalls of each class.\n>\n> It is important to note that, following similar logic to what has already been mentioned above, this metric is also harsh in that it categorises the predictions in a 'one-vs-rest' fashion - this does not adequately reflect the level of difficulty in correctly classifying each applicant with the exact rating they were given in the validation dataset.\n>\n> We can derive the balanced accuracy scores of each model as follows:","metadata":{}},{"cell_type":"code","source":"BalancedAccScores = []\n\n# Calculate the balanced accuracy score for each model's set of predictions.\nfor y_valid_preds in Valid_Preds:\n    BalancedAccScore = balanced_accuracy_score(y_valid, y_valid_preds)\n    BalancedAccScores.append(BalancedAccScore)\n    \n# Collect all Balanced Accuracy scores in a single dictionary.    \nBalancedAccScoreResults_y_valid = {'Model1_Calibrated': BalancedAccScores[0],\n                                   'Model2_Calibrated': BalancedAccScores[1],\n                                   'Model3_Calibrated': BalancedAccScores[2],\n                                   'Model4_Calibrated': BalancedAccScores[3],\n                                   'Model5_Calibrated': BalancedAccScores[4],\n                                   'Model6_Calibrated': BalancedAccScores[5],\n                                   'Model7_Calibrated': BalancedAccScores[6],\n                                   'Model8_Calibrated': BalancedAccScores[7],\n                                   'Model9_Calibrated': BalancedAccScores[8],\n                                   'Model10_Calibrated': BalancedAccScores[9],\n                                   'Model11_Calibrated': BalancedAccScores[10],\n                                   'Model12_Calibrated': BalancedAccScores[11]}\n\n# Select the model with the highest Balanced Accuracy.\nModel_LowestBalancedAccScore_Valid = max(BalancedAccScoreResults_y_valid, key=BalancedAccScoreResults_y_valid.get)\nprint(Model_LowestBalancedAccScore_Valid)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:12.018877Z","iopub.execute_input":"2022-07-08T18:23:12.019371Z","iopub.status.idle":"2022-07-08T18:23:12.203566Z","shell.execute_reply.started":"2022-07-08T18:23:12.019322Z","shell.execute_reply":"2022-07-08T18:23:12.202410Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Step 13: Assess & determine the best performing model, based on the metrics discussed above\n\n> In the section above, we considered a variety of performance-based factors that are commonly involved in choosing the best model to use for generating a final set of predictions.\n>\n> Based on the discussions above RE: ROC Curve/Confusion Matrix/Balanced Accuracy Score analyses, we select Model 12 (`XGBClassifier`) as the best performing model to be used for testing purposes.","metadata":{}},{"cell_type":"code","source":"# Designate model 12 (XGBClassifier) as the best performing model of the cohort.\nBestModel = Model12_Calibrated","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:12.205183Z","iopub.execute_input":"2022-07-08T18:23:12.206130Z","iopub.status.idle":"2022-07-08T18:23:12.211148Z","shell.execute_reply.started":"2022-07-08T18:23:12.206088Z","shell.execute_reply":"2022-07-08T18:23:12.209539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Step 14: Generate predictions/probability estimates using the best model (using test data)\n\n> Using the best performing model as determined above, we can now generate a set of test predictions which will be used to evaluate/review the model's performance on unseen data.","metadata":{}},{"cell_type":"code","source":"# Use the best performing model to generate predictions (as probability estimates) on the test dataset.\nBestModel_Test_PredProb = BestModel.predict_proba(X_test_L1reg)\n\n## Convert the probability estimates into label-encoded predictions, then convert back to the original set of classes.\nBestModel_Test_Preds = np.argmax(BestModel_Test_PredProb, axis=1)\nBestModel_Test_Classes = BestModel.classes_\nBestModel_Test_Preds = [BestModel_Test_Classes[i] for i in BestModel_Test_Preds]\nBestModel_Test_Preds = le.inverse_transform(BestModel_Test_Preds)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:12.213011Z","iopub.execute_input":"2022-07-08T18:23:12.213412Z","iopub.status.idle":"2022-07-08T18:23:12.279433Z","shell.execute_reply.started":"2022-07-08T18:23:12.213377Z","shell.execute_reply":"2022-07-08T18:23:12.277874Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Step 15: Evaluate the performance of the final model (using test data)\n\n> In the section below, we will review how the selected model has performed in terms of classifying applicants into each Response group, based on the test dataset.","metadata":{}},{"cell_type":"code","source":"# Convert the test dataset's encoded labels (0-7) back to the original set of classes (1-8), for clearer reviews of model performance.\ny_test = le.inverse_transform(y_test)\n\n# Create a duplicate copy of the model's predicted labels.\nTest_Preds = BestModel_Test_Preds\n\n# Create a duplicate copy of the model's prediction probabilities.\nTest_PredProbs = BestModel_Test_PredProb","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:12.281211Z","iopub.execute_input":"2022-07-08T18:23:12.281768Z","iopub.status.idle":"2022-07-08T18:23:12.290413Z","shell.execute_reply.started":"2022-07-08T18:23:12.281673Z","shell.execute_reply":"2022-07-08T18:23:12.289100Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### ROC Curve","metadata":{}},{"cell_type":"code","source":"# Plot the ROC curve to display the selected model's performance - TP rate vs FP rate.\nplot_roc(y_test, Test_PredProbs)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:12.292283Z","iopub.execute_input":"2022-07-08T18:23:12.293378Z","iopub.status.idle":"2022-07-08T18:23:12.647484Z","shell.execute_reply.started":"2022-07-08T18:23:12.293323Z","shell.execute_reply":"2022-07-08T18:23:12.646163Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> As shown in the chart above, the `XGBClassifier` model features a similarly broad range of AUC values against each class, indicating a varying degree of sensitivity/specificity across the dataset. Furthermore, the model has also performed relatively strongly when predicting applicants for classes 3/4/8, as indicated by the high AUC scores (~0.9) for their respective ROC curves - this is a good sign that our model is still able to generalise well, even when handling previously unseen data.","metadata":{}},{"cell_type":"markdown","source":"### Classification Report (i.e. Precision/Recall/F1-Score)","metadata":{}},{"cell_type":"code","source":"# Print a Classification Report to display the selected model's performance - Precision/Recall/F1-score/Accuracy.\n\nprint(classification_report(y_test, Test_Preds))","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:12.649215Z","iopub.execute_input":"2022-07-08T18:23:12.650555Z","iopub.status.idle":"2022-07-08T18:23:12.689538Z","shell.execute_reply.started":"2022-07-08T18:23:12.650508Z","shell.execute_reply":"2022-07-08T18:23:12.688512Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> As shown in the classification report above, the model shows a broad range of precision/recall/F1-score values against each class, indicating a varying degree of accuracy across the dataset and emphasising how the model has been generally well-fitted towards the most important `Response` groups. The model still appears to perform best when predicting for class 8 applicants, as evidenced by the relatively higher precision/recall/F1-score values for this group.","metadata":{}},{"cell_type":"markdown","source":"### Confusion Matrix","metadata":{}},{"cell_type":"code","source":"# Plot a Confusion Matrix to display the selected model's performance - True labels vs False labels.\n\nplot_confusion_matrix(y_test, Test_Preds)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:12.691185Z","iopub.execute_input":"2022-07-08T18:23:12.691604Z","iopub.status.idle":"2022-07-08T18:23:13.274929Z","shell.execute_reply.started":"2022-07-08T18:23:12.691567Z","shell.execute_reply":"2022-07-08T18:23:13.273940Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> From the confusion matrix above, it is clear that the model has performed on the test dataset very similarly to the validation scenario discussed earlier above. Furthermore, the model has also been relatively successful in replicating the `Response` distribution shown earlier in Step 1, \n>\n> Any off-diagonal hotspots (representing notable proportions of incorrect predictions) have already been revealed as part of the model validation process earlier above - so this should not come as a huge surprise. The important thing to note is that, with the exception of the observations in the top-right corner (where the model has predicted higher risk levels than required), the model appears to be well-fitted towards handling a multitude of applicants with varying levels of risk.","metadata":{}},{"cell_type":"markdown","source":"### Balanced Accuracy Score\n\n> Finally, we can calculate the balanced accuracy score of the selected model, using the test dataset:","metadata":{}},{"cell_type":"code","source":"# Calculate the balanced accuracy score for the selected model.\nBalancedAccScore_y_test = balanced_accuracy_score(y_test, Test_Preds)\n\nprint(BalancedAccScore_y_test)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:13.276445Z","iopub.execute_input":"2022-07-08T18:23:13.277392Z","iopub.status.idle":"2022-07-08T18:23:13.300786Z","shell.execute_reply.started":"2022-07-08T18:23:13.277325Z","shell.execute_reply":"2022-07-08T18:23:13.299646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Step 16: Review model performance in the context of \"feature importance\" and explainability\n\n> In the section above, we considered how well the model has performed against a number of key metrics such as ROC-AUC scores and precisions/recalls. However, we also need to think about how the model's decisions/behaviour can be explained in a clear and consistent manner - this is important in a real-world scenario from both a legal and ethical perspective, as we are assessing the prospective risk levels of numerous applicants for life insurance policies. Agreeing/rejecting applicants on insufficient/incorrect grounds could otherwise lead to reputational/financial damages.\n>\n> In this section, we will now consider:\n> * what features have the biggest impact on the `XGBClassifier` model's predictions, and\n> * how the `XGBClassifier` model has generated its predictions based on the values of the most important features.","metadata":{}},{"cell_type":"markdown","source":"### Feature Importances\n\n> Firstly, we will take a look at the `XGBClassifier` a bit more closely, in order to understand what features the model values most highly.","metadata":{}},{"cell_type":"code","source":"# Review the feature importances of the XGBoost classifier.\n\nFeatureImportances_BestModel = BestModel.base_estimator.feature_importances_\nFeatureImportances_BestModel","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:13.302490Z","iopub.execute_input":"2022-07-08T18:23:13.303170Z","iopub.status.idle":"2022-07-08T18:23:13.313326Z","shell.execute_reply.started":"2022-07-08T18:23:13.303131Z","shell.execute_reply":"2022-07-08T18:23:13.312322Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from xgboost import plot_importance\n\n# Generate a Feature Importance plot using the selected model.\n\nfig, ax = plt.subplots(figsize=(10, 10))\nplot_importance(BestModel.base_estimator,\n                importance_type=\"gain\",\n                xlabel=\"Gain\",\n                show_values=False,\n                ax=ax)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:13.315181Z","iopub.execute_input":"2022-07-08T18:23:13.315602Z","iopub.status.idle":"2022-07-08T18:23:14.884680Z","shell.execute_reply.started":"2022-07-08T18:23:13.315565Z","shell.execute_reply":"2022-07-08T18:23:14.883120Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The Gain chart above demonstrates the relative contribution of each feature to the model - which is calculated by taking each feature's contribution for each decision tree that is used in the model.\n>\n> This plot shows that the top five most important features (in terms of Gain/model contribution) are:\n>* `Medical_History_23`, \n>* `Medical_History_4`, \n>* `Medical_Keyword_3`, \n>* `Medical_Keyword_15`, and\n>* `BMI`.\n>\n> Interestingly, these five features were also highly ranked within the Mutual Information score chart (see Step 7 earlier above), which lends some credibility to the view that these features may be closely involved in governing what risk rating an applicant should be assigned. \n>\n> The bottom five (i.e. least important) features in terms of Gain are:\n>* `Medical_History_34`,\n>* `Product_Info_2_E1`,\n>* `KMeansCluster_4`,\n>* `Insurance_History_8`, and\n>* `Medical_History_41`.\n>\n> Four of these five features were also poorly ranked within the MI score chart shown in Step 7, however `KMeansCluster_4` was instead moderately ranked (residing within the top 15 features). This implies that, whilst `KMeansCluster_4` was initially deemed to show some potential in terms of predictive power, the `XGBClassifier` does not value this feature as highly when generating predictions for the test dataset.","metadata":{}},{"cell_type":"markdown","source":"### Permutation Importance\n\n> Next, we can review how the model handles random column-wise reordering of data - performing what is more commonly known as permutation importance analysis - to see which features increase the volatility of our model's predictions upon random shuffling (and hence, which features the model relies on most heavily for generating predictions).\n>\n> To do this, we will calculate the permutation importance weights using our test dataset, as shown below:","metadata":{}},{"cell_type":"code","source":"from eli5.sklearn import PermutationImportance\nfrom eli5 import show_weights\n\n# Calculate the Permutation Importances of the selected model.\n\nperm = PermutationImportance(BestModel.base_estimator, random_state=0).fit(X_test_L1reg, y_test)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:14.886422Z","iopub.execute_input":"2022-07-08T18:23:14.886862Z","iopub.status.idle":"2022-07-08T18:23:36.851886Z","shell.execute_reply.started":"2022-07-08T18:23:14.886824Z","shell.execute_reply":"2022-07-08T18:23:36.850754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Show the Permutation Importance weights (plus errors) of the top 20 features.\n\nshow_weights(perm, feature_names=X_test_L1reg.columns.tolist())","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:36.853811Z","iopub.execute_input":"2022-07-08T18:23:36.855385Z","iopub.status.idle":"2022-07-08T18:23:36.914782Z","shell.execute_reply.started":"2022-07-08T18:23:36.855324Z","shell.execute_reply":"2022-07-08T18:23:36.912854Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The values towards the top of the table are the most important, and those towards the bottom matter the least.\n>\n> Here, the top five features that appear to have the strongest impact on the model's performance/accuracy when shuffled randomly are:\n>* `Medical_Keyword_3`\n>* `Medical_History_39`\n>* `InsuredInfo_5`\n>* `Product_Info_2_D1`\n>* `Medical_History_17`\n>\n> Most of these features were also highlighted as having the highest feature importances in the section above. Interestingly, `Medical_History_23` (the highest scored feature in terms of Feature Importance) does not show up in the top 20 permutation importance weightings - this could imply that this feature does not have any meaningful causational relationship with `Response`, and has unintentionally been given higher importance due to possible overfitting.\n>\n> More information on the topic of comparing tree-based feature importances against permutation importances can be found at the following page: [Relation to impurity-based importance in trees - scikit-learn](https://scikit-learn.org/stable/modules/permutation_importance.html#relation-to-impurity-based-importance-in-trees).","metadata":{}},{"cell_type":"markdown","source":"### SHAP Values\n\n> The `SHAP` framework (acronym derived from **SH**apley **A**dditive ex**P**lanations) is excellent for breaking down predictions in order to show the impact of each feature. This is especially useful for when we want to explain how/why we have classified certain applicants who may have demonstrated high/low risk potential.\n>\n> The advantage that this approach has over calculating permutation importances alone, is that `SHAP` values can express whether a feature has a broad effect across all predictions, or whether its effect is more localised for a handful of predictions and negligible in general; permutation importance simply captures the \"average\" impact of each feature.\n>\n> Here, we will calculate the SHAP values for the test dataset to visualise how the `XGBClassifier` model has behaved, based on the values of the most important features.","metadata":{}},{"cell_type":"code","source":"import shap\n\n# Initialise the explainer object in order to calculate SHAP values.\nexplainer = shap.TreeExplainer(BestModel.base_estimator)\n\n# Calculate SHAP values for the whole test dataset.\nshap_values = explainer.shap_values(X_test_L1reg)\n\n# Create a summary plot of the SHAP values.\nshap.summary_plot(shap_values[1], X_test_L1reg)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:36.916430Z","iopub.execute_input":"2022-07-08T18:23:36.916853Z","iopub.status.idle":"2022-07-08T18:23:49.309843Z","shell.execute_reply.started":"2022-07-08T18:23:36.916817Z","shell.execute_reply":"2022-07-08T18:23:49.308841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The SHAP summary plot above provides us with a top-down view of feature importance, for the top 20 features as calculated via the `SHAP` framework. The colour of each dot represents whether that feature was high or low (for that row in the dataset), and its horizontal location shows whether the effect of that value caused a higher (towards 8) or lower (towards 1) prediction.\n>\n> For instance, the `BMI` feature clearly expresses that as the BMI of the applicant increases, then their predicted risk rating also increases strongly. As another example, `Product_Info_2_A6` - i.e. whether the applicant's selection for `Product_Info_2` is equal to **A6** - shows a clear negative correlation with the predicted risk rating.\n>\n> However, there are also a couple of other interesting cases to note - for example, `Ins_Age` shows a negative correlation when away from the baseline value (class 5) but is more mixed as the value approaches the baseline. `Family_Hist_4` appears to show a mixed effect on the value of `Response` regardless of the choice of class.","metadata":{}},{"cell_type":"markdown","source":"### SHAP Dependence Plots\n\n> We can also review individual features more closely, by creating SHAP dependence contribution plots.\n>\n> These are extremely helpful for displaying what the distribution of effects is, and linking this to a model's predictions.","metadata":{}},{"cell_type":"code","source":"# Create a SHAP Dependence plot of the highest ranked feature - \"BMI\".\n\nshap.dependence_plot(\"rank(0)\", shap_values[1], X_test_L1reg)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:49.311165Z","iopub.execute_input":"2022-07-08T18:23:49.312085Z","iopub.status.idle":"2022-07-08T18:23:51.494076Z","shell.execute_reply.started":"2022-07-08T18:23:49.312039Z","shell.execute_reply":"2022-07-08T18:23:51.492700Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> Here, we have generated a SHAP dependence contribution plot for `BMI` as a function of `Ins_Age`. The horizontal location of each dot represents the actual value from the test dataset, and the vertical location represents the impact (on the prediction) of having that value.\n>\n> The shape of the plot indicates that, once `BMI` exceeds 0.7, the distribution slopes sharply upwards - inferring that the likelihood of being assigned to a higher risk rating greatly increases after this point. However, as `BMI` approaches 1.0, the distribution then quickly slopes downwards. The drastic increase in the distribution's spread as the feature value approaches 1.0 implies that other features begin to interact more strongly with `BMI` at this upper limit.\n>\n> Upon reviewing the interaction between `BMI` and `Ins_Age` more closely at this upper limit, it appears that on average, younger applicants with high BMIs tend to be assigned higher risk ratings by the model than high-BMI applicants that are more elderly, however more data would be required in order to validate this hypothesis, before any relevant business decisions could be made. ","metadata":{}},{"cell_type":"code","source":"# Create a SHAP Dependence plot of the 2nd highest ranked feature - \"Product_Info_4\".\n\nshap.dependence_plot(\"rank(1)\", shap_values[1], X_test_L1reg)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:51.495712Z","iopub.execute_input":"2022-07-08T18:23:51.496418Z","iopub.status.idle":"2022-07-08T18:23:53.663111Z","shell.execute_reply.started":"2022-07-08T18:23:51.496359Z","shell.execute_reply":"2022-07-08T18:23:53.661743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The shape of this plot indicates a broad range of interaction between `Product_Info_4` and `Response` within the test dataset. However, the general trend is that at low values of `Product_Info_4`, the likelihood of being assigned to a higher risk rating is increased - however a couple of notable exceptions to this are when the distribution's spread is at its widest, where `Product_Info_4` is close to either 0.0, 0.5 or 1.0.\n>\n> The interaction between `Product_Info_4` and `Medical_History_13` seems to vary quite significantly - however a notable trend that can be observed is that, for a given value of `Product_Info_4`, applicants with lower values of `Medical_History_13` tend to be rated as less risky than their counterparts with higher values of `Medical_History_13`. This trend is especially notable when `Product_Info_4` approaches either 0.0, 0.5 or 1.0.","metadata":{}},{"cell_type":"code","source":"# Create a SHAP Dependence plot of the 3rd highest ranked feature - \"Product_Info_2_A6\".\n\nshap.dependence_plot(\"rank(2)\", shap_values[1], X_test_L1reg)","metadata":{"execution":{"iopub.status.busy":"2022-07-08T18:23:53.665292Z","iopub.execute_input":"2022-07-08T18:23:53.666112Z","iopub.status.idle":"2022-07-08T18:23:55.919477Z","shell.execute_reply.started":"2022-07-08T18:23:53.666058Z","shell.execute_reply":"2022-07-08T18:23:55.918234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> This plot displays that, when the applicant's selection for `Product_Info_2` is equal to **A6**, the predicted risk rating is generally much lower on average (as already mentioned in the summary plot), whereas if it is not equal to **A6** then the risk rating is extremely likely to stay at the baseline (class 5).\n>\n> Although it is difficult to spot visually, `Product_Info_2_A6` does also show some interaction with `Medical_Keyword_3`. In cases where `Product_Info_2_A6` feature is equal to 1.0, applicants with lower values of `Medical_Keyword_3` tend to shift the model's predictions away from class 1, and slightly closer towards class 2.","metadata":{}},{"cell_type":"markdown","source":"# Step 17: Review for areas of improvement\n\n> Whilst this project has aimed to showcase an implementation of the DS/ML workflow in an insurance-based context, it does not aim to provide a comprehensive \"A-Z\" approach on how to apply ML techniques for the risk classification of insurance applicants. There are, however, a couple of areas where this project could be improved.\n>\n> For example, in terms of **model performance**:\n> * The aim of this project has been to establish an ML model that is capable of accurately predicting risk ratings (split between 8 different categories) for insurance applicants. However, as this involves working with a fairly complex dataset with a wide range of features, the model has not been fully able to capture every single nuance within the dataset, and has instead generalised to predict across the entire distribution to a reasonable degree of accuracy. Therefore, the model may actually prove to be more useful for business purposes if it was instead assigned to predict based on broader sets of risk bandings (e.g. classes 1-3 could be grouped together into Low, classes 4-6 into Medium, and classes 7-9 into High). Then, further refinement of these initial groupings could be made by assembling another layer of classifiers that are specifically trained on datasets that are representative of each risk banding (e.g. another model specifically trained for performing further segmentation of low-risk applicants).\n>\n> In terms of **feature engineering**:\n> * As this dataset's feature names have been anonymised, it has not been possible to incorporate \"business logic\"-driven feature engineering directly into the model preparation process. However, potential interactions between features could still be studied further during EDA in order to produce new features that take advantage of statistical correlations within the dataset. With the assistance of an \"unanonymised\" dataset as well as the input of subject matter expertise, it may be possible to create more predictive features and validate them for use in production.\n>\n> Lastly, in terms of **feature selection**:\n> * More stringent and quantitative-based approaches could have been used, in order to combine the methods employed for selecting the list of features that would be supplied to each model during the training process. For clarity, whilst we used Mutual Information/Variance Inflation Factor/Principal Component analysis techniques to display how the importance of each feature varied throughout the dataset, we only relied on Lasso regularisation in order to produce the final cut of selected features. It may be worth considering the use of a \"voting system\" in order to select features based on how they are ranked across each of the techniques used, so that one technique alone does not provide a final view.","metadata":{}},{"cell_type":"markdown","source":"# Acknowledgements\n\n> In addition to the sources cited earlier above, I would also like to thank **Prudential** for publishing the life insurance applicant dataset that was considered in this project.\n> \n> Additional information regarding the dataset that was used throughout this project, as well as the affiliated 2016 competition \"**Prudential Life Insurance Assessment**\" hosted on Kaggle, can be found [here](https://www.kaggle.com/competitions/prudential-life-insurance-assessment/).","metadata":{}}]}