{"cells":[{"metadata":{"id":"EIMxxH2V9xMz"},"cell_type":"markdown","source":"# **Categorical Feature Encoding Challenge II**\n\nby Robin Munier\n\n---\n\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp;The challenge is to achieve a binary classification with every feature a categorical one. Besides, no information is available about what the target and the features represent (except what we can deduce from the values they take). Thus, the variables names are anonymized. The performance measure to maximize on this competition is the area under the ROC curve.\n\n&nbsp;&nbsp;&nbsp;&nbsp;Comparing to the first Categorical feature encoding challenge, this one adds the additional complexity of feature interactions, as well as missing data.\n\n&nbsp;&nbsp;&nbsp;&nbsp;The purpose of this competition is to develop its encoding and preprocessing skills. It is the opportunity to try different encoding schemes for different algorithms and to compare how they perform. Last but not least, of course, it is also an occasion to share our knowledge, experience and thoughts with members of the Kaggle community. So here is my proposal on this competition."},{"metadata":{"id":"hPgKQidb8pG6"},"cell_type":"markdown","source":"## **Table of content**\n\n[Introduction](#introduction_id)<br>\n\n[I. Preliminaries](#I)<br>\n\n[II. Missing data processing](#II)<br>\n&nbsp;&nbsp;&nbsp;&nbsp; [II.1 The missing data of our data set](#II.1)<br>\n&nbsp;&nbsp;&nbsp;&nbsp; [II.2 Different missing data mechanisms](#II.2)<br>\n&nbsp;&nbsp;&nbsp;&nbsp; [II.3 Missing data exploration strategy](#II.3)<br>\n&nbsp;&nbsp;&nbsp;&nbsp; [II.4 Tests of independence](#II.4)<br>\n&nbsp;&nbsp;&nbsp;&nbsp; [II.5 Relationships between more than 2 features](#II.5)<br>\n&nbsp;&nbsp;&nbsp;&nbsp; [II.6 Strategy adopted](#II.6)<br>\n\n[III. Encoding, modeling and interpretation](#III)<br>\n&nbsp;&nbsp;&nbsp;&nbsp; [III.1 First strategy](#III.1)<br>\n&nbsp;&nbsp;&nbsp;&nbsp; [III.2 Second strategy](#III.2)<br>\n&nbsp;&nbsp;&nbsp;&nbsp; [III.3 Third strategy](#III.3)<br>\n&nbsp;&nbsp;&nbsp;&nbsp; [III.4 Ensembling](#III.4)<br>\n\n[IV. Conclusion and potential extentions](#IV)<br>\n\n[V. Bibliographic references](#V)<br>"},{"metadata":{"id":"FqoeSjllAWlJ"},"cell_type":"markdown","source":"## **Introduction** <a id='introduction_id'></a>\n\n&nbsp;&nbsp;&nbsp;&nbsp;This notebook is an overview of the useful part my work on this challenge. <br>\n&nbsp;&nbsp;&nbsp;&nbsp;It does not include a particular part on exploratory data visualization. For clarity, the visualizations are instead included along the reasoning thread when they bring something useful.<br>\n&nbsp;&nbsp;&nbsp;&nbsp;Some notebooks present vast EDA to finally move on to the modeling part without getting anything out of it. I don't want to make visualizations just to decorate. So, for concision's sake, this notebook does not show all the exploration and visualizations realized, just the ones which guided some useful choices. <br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp;This work has the ambition of being beginner friendly while being enough elaborated to interest more advance readers. As rigor and clarity are not synonymous, I have tried to strike a good balance in the textual parts. In particular, I have chosen to spend less time to explain the modeling than the missing values management, the encoding and the interpretation of the results. Indeed, after having read many notebooks and discussion threads on Kaggle, I feel the contents dealing with the preprocessing are rarer than the ones dealing with modeling. But, as stated Larry Wasserman in the preface of his book All of Statistics (Springer) : \"Using fancy tools like neural nets, boosting, and support vector machines without understanding basic statistics is like doing brain surgery before knowing how to use a band-aid\".<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp;Some bibliographic references are available at the end. They contains ressources dealing with the concepts I evoke. In this way, if some of my explanations are unclear to you, you can check them and maybe they will suit you more. The outputs (the dataframes, the models, the submission files, etc.) of the notebook have been saved and are available in the part \"Output Files\" below.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp;Feel free to share all your feedbacks and remarks, publicly or privately. I am very open to criticism and I take it as an opportunity of improvement. What's more, I am not a native speaker of english, so I hope that I have not mistreated this language.<br>\n&nbsp;&nbsp;&nbsp;&nbsp;If this notebook pleased you in any way, or learned you anything, you can upvote it (don't be timid, it is free). As I plan to write other notebooks in the future, you can follow me to be informed when that will happen."},{"metadata":{"id":"yqk06XgzNTN5"},"cell_type":"markdown","source":"## **I. Preliminaries** <a id='I'></a>\n<br>\nInstallation of some libraries."},{"metadata":{"id":"X-g3sGa8pzTw","outputId":"69d82551-739c-4ea9-99d7-97e5e6858eae","trusted":true},"cell_type":"code","source":"!pip install --upgrade xgbfir","execution_count":null,"outputs":[]},{"metadata":{"id":"q-D2LpemRwXV"},"cell_type":"markdown","source":"Import of some libraries."},{"metadata":{"id":"gbqE6nJBXemp","outputId":"691598ce-62c1-4655-fff8-ce84da0a3a6a","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport sklearn as sk\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport seaborn as sbn\nimport sys\nfrom datetime import datetime\nimport category_encoders as ce\nfrom scipy.stats import chi2_contingency\nfrom statsmodels.stats.multitest import fdrcorrection, multipletests\nfrom sklearn.model_selection import train_test_split, StratifiedKFold\nfrom sklearn.metrics import roc_auc_score, roc_curve, recall_score, accuracy_score, precision_recall_curve\nfrom sklearn.preprocessing import OneHotEncoder, StandardScaler, LabelEncoder\nimport xgboost as xgb\nimport xgbfir\nimport string\nimport joblib\nfrom category_encoders import TargetEncoder, OrdinalEncoder\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.naive_bayes import CategoricalNB","execution_count":null,"outputs":[]},{"metadata":{"id":"Q2ZCCMIaR483"},"cell_type":"markdown","source":"Let's display the versions for reproducibility."},{"metadata":{"id":"qbcxdat7R4Aq","outputId":"c54f8031-2f18-40db-c687-65d7727784fe","trusted":true},"cell_type":"code","source":"print(\"python\", sys.version)\nfor module in np, pd, mpl, sbn, ce, sk, xgb, joblib:\n    print(module.__name__, module.__version__)","execution_count":null,"outputs":[]},{"metadata":{"id":"f61FQMeLSVi-"},"cell_type":"markdown","source":"Initialisation of the random number generators (scikit-learn uses the numpy seed)."},{"metadata":{"id":"WLJnrKoPfCdp","trusted":true},"cell_type":"code","source":"np.random.seed(0)","execution_count":null,"outputs":[]},{"metadata":{"id":"l_ZBY1fTSSYM"},"cell_type":"markdown","source":"Data import."},{"metadata":{"id":"jmqa4Jud0Oi-","trusted":true},"cell_type":"code","source":"df = pd.read_csv('../input/cat-in-the-dat-ii/train.csv').drop('id',axis=1)\ntest = pd.read_csv(r'../input/cat-in-the-dat-ii/test.csv', index_col=0)\nsample_submission = pd.read_csv(r'../input/cat-in-the-dat-ii/sample_submission.csv', index_col=0)","execution_count":null,"outputs":[]},{"metadata":{"id":"dW49LtQfKtLF"},"cell_type":"markdown","source":"The submission file contains a dataframe with an id key for each test observations and only one columns: the target. It is filled with a constant value of 0.5. On each row, we have to replace it with our estimated probability for the sample to be of class 1."},{"metadata":{"id":"EBM3xTZsKNLs","outputId":"797d8ace-40f1-432a-c0a7-ff4d44ff1a58","trusted":true},"cell_type":"code","source":"sample_submission.head()","execution_count":null,"outputs":[]},{"metadata":{"id":"acptuBpGM4Y7"},"cell_type":"markdown","source":"The train and test files contain respectively the training and test sets. The second one contains the observations for which the probabilities must be estimated and submitted on Kaggle to participate to the challenge.<br>\n\nLet's have a first look at the data."},{"metadata":{"id":"UL7EvZv_zSuF","outputId":"227ec79f-369e-4588-ba82-241a01ae9039","trusted":true},"cell_type":"code","source":"df.head().T","execution_count":null,"outputs":[]},{"metadata":{"id":"6abmMExOMvVN","outputId":"bb8910a0-c7cb-4e68-ac08-8188e8c1f6f8","trusted":true},"cell_type":"code","source":"test.head()","execution_count":null,"outputs":[]},{"metadata":{"id":"J70BNGVN0jXn","outputId":"bb008fa0-36bc-46da-83c1-b5f88e8b9fcc","trusted":true},"cell_type":"code","source":"print('Shape of the training set: ', df.shape)\nprint('Shape of the test set: ', test.shape)","execution_count":null,"outputs":[]},{"metadata":{"id":"c3VKj89PtfYL"},"cell_type":"markdown","source":"## II. Missing data processing.<a id='II'></a> <br>\n\n### II.1 The missing data of our data set.<a id='II.1'></a> <br>\n\nAll the features are from the expected dtype. Each one has approximately 18000 missing values."},{"metadata":{"id":"KzX8OvmJQVWO","outputId":"ab7b0f18-7110-4818-c7df-24e24d701e22","trusted":true},"cell_type":"code","source":"df.info()","execution_count":null,"outputs":[]},{"metadata":{"id":"HpvW64ho1M1z"},"cell_type":"markdown","source":"About half of the rows hold at least one missing value."},{"metadata":{"id":"e0H3lRNtzkIn","outputId":"6466329a-8a13-4b1e-daed-d6dd79f07d74","trusted":true},"cell_type":"code","source":"(df.isna().sum(axis=1)==0).sum()/df.shape[0]","execution_count":null,"outputs":[]},{"metadata":{"id":"IlmNqKhtzrbJ","outputId":"c675f54b-3355-43d2-bccd-c92f551e7b1b","trusted":true},"cell_type":"code","source":"(test.isna().sum(axis=1)==0).sum()/test.shape[0]","execution_count":null,"outputs":[]},{"metadata":{"id":"Akhmy-3fMdY5"},"cell_type":"markdown","source":"When some rows have too many missing values, we can consider deleting them. Here, every row contains less than 8 missing values. So we would probably discard more information than noise by deleting it. Moreover, to drop rows can bias the dataset if the values are not randomly missing."},{"metadata":{"id":"Q2myTvpuMGfw","outputId":"b6cd42f6-2329-412c-ae9e-3e20d130543c","trusted":true},"cell_type":"code","source":"df.isna().sum(axis=1).value_counts()","execution_count":null,"outputs":[]},{"metadata":{"id":"Y_oCX5CodNAn"},"cell_type":"markdown","source":"### II.2 Different missing data mechanisms.<a id='II.2'></a> <br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; Understanding the reasons why data are missing helps to handle them. Therefore, let's remind the generally considered missing data mechanisms.<br>\n\n*   Missingness completely at random (MCAR): the values of a variable are missing completely at random if the propensity of missingness is the same for each value. When data are MCAR, the analysis performed on the non-missing data are unbiased. Unfortunately, data are rarely MCAR.<br>\n\n*   Missing at random (MAR) occurs when the probability a variable is missing depends only on available information, that is to say on the other variables of the dataset (the expression \"at random\" refers to the random impact of the unavailable information on the missingness). <br>\n\n*   Missing not at random (MNAR): this mechanism includes the cases where the missingness depends on unobserved predictors or on the missing value itself (whose a special case is censoring : it happens for example if people with a salary superior to a certain threshold refuse to communicate on its value).<br>\n\n\n&nbsp;&nbsp;&nbsp;&nbsp; Of course, the missing data mechanisms can be different for different variables of our dataset. We can never be sure of the mechanism by which the data are missing. Furthermore, data are never MCAR or even MAR. If a data is missing, there is necessarily a cause (or even several) as nothing happens by chance. This cause is in fact a variable linked to the missingness. But we are looking for a reasonable (in the sense that it is coherent and it conforms well to our observations) assumption more than the strict reality.\n<br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; We should usually try to include new features in our dataset if possible. Generally, more we add variables in our dataset, more reasonable is the MAR assumption. Even if these features are not necessarily efficient for the final task (whatever it is, prediction or clustering), they can be useful to handle the NaN (it is always possible to drop them once the NaN have been managed).\n<br>\n<br>\n\n### II.3 Missing data exploration strategy.<a id='II.3'></a><br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; To identify the mechanisms involved, we can use knowledge about the data (their meaning, how they were collected, etc.) and investigate them to find patterns. Here, as mentionned previously, we have no knowledge about data. So we are bound to investigate patterns. <br>\n<br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; The exploration of the relationship between two categorical variables is usually done with contingency tables. For example, below is the contingency table between nom_1 and nom_2. As all the categories are not equally represented, each value in the table is divided by the number of observations in the category corresponding to the row (so we can compare the values between the rows)."},{"metadata":{"id":"y9fOa2NCtfYE","outputId":"f170515e-24f8-4f57-a03d-04c8552b5c0b","trusted":true},"cell_type":"code","source":"cont_table = pd.crosstab(df.nom_1, df.nom_2, normalize='index')\ncont_table","execution_count":null,"outputs":[]},{"metadata":{"id":"5bx57yRYBFze"},"cell_type":"markdown","source":"Sometimes, a heatmap of the contingency table is clearer (maybe not in that case)."},{"metadata":{"id":"MEnDi9IjZtp6","outputId":"29f6a251-a4f0-40d5-d5d4-c0ff572b31e6","trusted":true},"cell_type":"code","source":"plt.matshow(cont_table, cmap=plt.cm.gray)","execution_count":null,"outputs":[]},{"metadata":{"id":"kKZL-DJ6tfYJ"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; Here, we observe that nom_2 takes more often the value Dog when nom_1 takes the value Polygon than when nom_1 takes the value Square. But how to know if this observation is enough to conclude to a link between the two variables ? <br>\n&nbsp;&nbsp;&nbsp;&nbsp; We can't do it properly just by visualizing the data. Instead, we need to quantify the pattern observed. The classical way to check (or at least to estimate properly) if a pattern is due to sampling fluctuation or to an actual link is to perform a statistical test. Thus, we will not display a myriad of contingency tables. It would be useless. Instead, we will perform statistical tests.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; The famous Little's test of MCAR is not suitable for categorical variables. So, for each feature, we will firstly define an indicator variable taking the value 1 if the feature is missing and 0 otherwise and then we will perform one chi-squared test of independence between this indicator variable and each other available feature (and also one with the target as it can reveal useful patterns). <br>\n&nbsp;&nbsp;&nbsp;&nbsp; A chi-squared test makes it possible to check the absence of statistical link between two variables $X$ and $Y$. They are said to be independent when there is no statistical link between them. In other words, when knowing X does not offer any information about Y. <br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; If these statistical tests of independence suggest the missingness of one variable is not independent from the values taken by at least one of the others, then we can argue this variable is not MCAR.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; However, the missingness of a variable can depend on the interaction between several other variables. These multiple interactions are not considered by the pairwise tests of independence evoked. Thus, to take them into account, we will next use gradient boosted trees (as they allow interactions) to predict the missingness of each variable thanks to the others.<br>\nIf the variables present in our dataset predict well enough the missingness, we could argue it is because the missingness depends strongly on available information (so the data can be considered MAR). <br>\nIf the prediction is poor, we can support the opposite. Nevertheless, the cause could be different (maybe the model was not well suited). When the data are not MAR, it is hard to check the MCAR and MNAR assumptions.<br>"},{"metadata":{"id":"ohyHIfAyBFzl"},"cell_type":"markdown","source":"### II.4 Tests of independence.<a id='II.4'></a><br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; In the general framework, we have two variables $X$ and $Y$ which are supposed to take a finite number of possible values, $I$ for $X$, $J$ for $Y$. We have a sample of $N$ data. Consider testing $H0: X\\: and\\: Y\\: are\\: independent$ versus $H1: X\\: and\\: Y\\: are\\: not\\: independent$.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Let's denote $p(i)$, $p(j)$ and $p(i,j)$ the respective theoretical probabilities that $X$ equals $i$, $Y$ equals $j$ and $(X,Y)$ equals $(i,j)$ (that is to say $X$ equals $i$ and $Y$ equals $j$ simultaneously). Now the hypothesis can be rewritten as $H0: p(i,j)=p(i) \\times p(j)$ versus $H1: p(i,j) \\neq p(i) \\times p(j)$.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Let's denote by $O(i)$ the observed number of data for which $X$ takes the value $i$, by $O(j)$ the observed number of data for which $Y$ takes the value $j$, and by $O(i,j)$ the observed number of data for which $(X,Y)$ takes the value $(i,j)$. Besides, if we denote by $E(i,j)$, $E(i)$ and $E(j)$ the expectations of $O(i,j)$, $O(i)$ and $O(j)$ respectively, then we have: $E(i) = p(i) \\times N$, $E(j)=p(j) \\times N$, and $E(i,j)=p(i,j) \\times N$.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Under $H0$, $p(i,j)=p(i)*p(j)$, so the maximum likelihood estimator of $p(i,j)$ is $\\hat{p}(i,j)=\\hat{p}(i)*\\hat{p}(j)= O(i)/N * O(j)/N$. So the maximum likelihood estimator of $E(i,j)$ is $\\hat{E}(i,j)=\\hat{p}(i,j) \\times N=O(i)*O(j)/N$.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; We calculate the distance $T$ between the observed values $O(i,j)$ (or empirical values) and the theoretical values $\\hat{E}(i,j)$ (that is to say our estimation of the expected values if there was independence) using the formula:<br>\n<br>\n$$T=\\sum_{i,j} \\frac{(O(i,j)-\\hat{E}(i,j))^2}{\\hat{E}(i,j)}$$<br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; Under $H0$, the chi-squared statistic $T$ follows asymptotically a chi-squared distribution with $(I – 1)*(J – 1)$ degrees of freedom. This fact can be used to approximate $T$ with a chi-squared distribution and thus to define a rejection region. An approximate level $\\alpha$ test is obtained by rejecting $H0$ if and only if $T$ is superior to the $\\alpha$-quantile of the chi-squared distribution with $(I – 1)(J – 1)$ degree of freedom.<br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; As $T$ follows a chi-squared law asymptotically, larger is the dataset, more acceptable is the approximation of $T$ with a chi-squared distribution.\nFor the test to be valid, we need the values of the contingency table between $X$ and $Y$ to be sufficiently large, meaning the dataset contains enough value of each category pair (larger is the dataset, better are the estimations, more reliable is the result). Diverse arbitrary criteria are usually used to impose constraints on the number of observations. For the vastly used Cochran's criterion to be satisfied, we need: $\\forall\\:(i,j)$, $\\hat{E}(i,j)>0$ and for at least 80% of the category pairs $(i,j)$, $\\hat{E}(i,j)>4$.<br>\n<br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; For each variable, we will apply this test between the missingness indicator of that variable and each other variable in our dataset. The null hypothesis will be rejected when the p-value is under $\\alpha$ = 5% (the threshold generally used).<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; In some cases, we can exclude by hand some pairs of variable that must not be linked (but for which a test could possibly return a low p-value by fluctuation of sampling). As in this competition, the meaning of the variables is hidden, we will roughly perform the tests with all of them.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; When the hypothesis to test are suggested by a given data set, they must be tested with another data set (to avoid data snooping). Here, as we will tests all possible pair dependencies, the choice of the tests is not based on a data exploration. So we don't have this problem.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Finally, we should keep in mind that the results of the tests are not an absolute truth. In particular, as we will perform each test at level 5%, then for each one we have 5% chance to reject the independence hypothesis when it is actually true (what is called the type I error). Moreover, we can also fail to reject false null hypothesis (but these type II errors are harder to estimate).<br>\n&nbsp;&nbsp;&nbsp;&nbsp; As we will conduct 529 tests (23 for each of the 23 features, as the target has no missing value), the probability that at least one is wrongly rejected is much larger than 5% (this is known as the multiple testing problem). As it depends on the number of true null hypothesis, we cannot know this probability. <br>\n&nbsp;&nbsp;&nbsp;&nbsp; However, we can control the false discovery rate with the Benjamini-Hochberg method. The false discovery proportion is the proportion of incorrect rejections among all the rejections (so it is the number of incorrect rejections divided by the number of rejections, which is different from the porportion of incorrect rejections among the true null hypothesis). The false discovery rate is the expectation of the false discovery proportion. Then, the Benjamini-Hochberg method consists in following steps:\n- Firstly, we order the p-values $P{(i)}$ we obtained by performing all the tests, as follows: $P_{(1)} < P_{(2)} < ... < P_{(m)}$ (where $m$=529 is the number of tests conducted).<br>\n- Then, for $i$ going from 1 to m, we define $l_{(i)} = \\frac{i}{m} \\times \\alpha$ (where $\\alpha$ is 5% is our case) and $R = \\max{\\{i: P_{(i)} < l_{(i)}\\}}$\n- The BH (Benjamini-Hochberg) threshold $T$ is defined by: $T=P_{(R)}$. In other words, T is the biggest p_value $P_{(i)}$ such that $P_{(i)} < l_{(i)}$.\n- Finally, each null hypothesis $H0_{(i)}$ (whose the p-value is $P_{(i)}$) is rejected if and only if $P_{(i)} \\leq T$.\n<br>\n\nWith this method, the false discovery rate is inferior to $\\alpha$ (whatever the number of actual true null hypothesis and the distribution of the p-values when the null hypothesis is false).<br>\nHence, next we will conduct the tests with the Benjamini-Hochberg method.<br>\n<br>\nWe first define a function that takes a contingency table and return a boolean indicating if the Cochran's criterion is satisfied for this contingency table."},{"metadata":{"id":"lT_1YqDJP4dp","trusted":true},"cell_type":"code","source":"def cochran_criterion(crosstab):\n  criterion = True\n  E = []\n  N = df.shape[0]\n  for O_i in crosstab.sum(axis=1):\n    for O_j in crosstab.sum():\n      if O_i*O_j/N==0:\n        critetion = False\n      E.append(O_i*O_j/N > 4)\n  criterion = criterion & (np.mean(E)>0.8)\n  return criterion","execution_count":null,"outputs":[]},{"metadata":{"id":"xoGoLObNm0U1"},"cell_type":"markdown","source":"Here, we create a function taking two variables as argument and returning a list holding the four following outputs :\n<br>\n\n*   A boolean indicating if the contingency table between the two input variables satisfies the Cochran's criterion.\n<br>\n*   The chi-squared statistic T calculated between the two input variables.\n<br>\n*   The p-value of the chi-squared test.\n<br>\n*   The significance of the chi-squared test alone (without considering the Benjamini-Hochberg method), namely a boolean equals to True if and only if the Cochran's criterion is satisfied and the p-value is under 5%.\n"},{"metadata":{"id":"LQRZVKlvQXlF","trusted":true},"cell_type":"code","source":"def chi2_test(X,Y):\n  crosstab = pd.crosstab(X, Y)\n  criterion = cochran_criterion(crosstab)\n  chi2, p = chi2_contingency(crosstab)[:2]\n  return [criterion, chi2, p]","execution_count":null,"outputs":[]},{"metadata":{"id":"Rd-MI38Noj3F"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; Here, we will perform the tests on the full dataset, that is to say the concatenation of the training and test sets (and not only on the training set). In the general case where we want to produce the model that will perform best when deployed on new data sets, it can be a cause of data snooping. Indeed, let's assume we construct a classification model, train it on the training data, test it on a test set we used (with the training set) for the statistical tests, and then test it again on new test sets we have never seen before. As some information from the first test set would be used for the statistical tests, the performance obtained on this test set would be an overestimation of the performance we would get on new test sets.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; But we are not in this case on this challenge. In fact, our performance on this classification competition is measured (by the ROC AUC score) only on the test set provided by Kaggle. So we don't care whether our performance on this test set is an overestimation of the perfomance of our model. Therefore, in this particular case, it is not problematic to use the test set for the independence tests.<br>\nBesides, as mentionned before, larger the dataset on which the tests are performed, more reliable are the results. Consequently, we realized the tests with all the data from both the training dataset and the test set. <br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; We create a dataframe containing the outputs of the chi2_test function for each pair of variables. As the target variable is only present on the training set, only the training set is used for the tests involving the target.<br>"},{"metadata":{"id":"qI6S3mhp7YAV","outputId":"645e5ba5-cd81-41ea-e4aa-529fde75e92b","trusted":true},"cell_type":"code","source":"start = datetime.now()\n\ndf_total = pd.concat([df.drop('target', axis=1), test])\nchi2_missingness = pd.DataFrame(columns=['Cochran_criterion', 'Chi2', 'p_value'])\n\nfor col1 in df_total.columns:\n  for col2 in df_total.columns.drop(col1):\n    missingness_indicator = df_total[col1].isna()\n    other_variable = df_total[col2]\n    chi2_missingness.loc[col1 + '_' + col2] = chi2_test(missingness_indicator, other_variable)\n\nfor col in df.columns.drop('target'):\n  missingness_indicator = df[col].isna()\n  chi2_missingness.loc[col + '_target'] = chi2_test(missingness_indicator, df.target)\n\nruntime = datetime.now() - start\nprint('Runtime ', runtime)","execution_count":null,"outputs":[]},{"metadata":{"id":"8oSxuaDsBFz1"},"cell_type":"markdown","source":"Now we use the Benjamini-Hochberg method and add the results to the dataframe."},{"metadata":{"id":"fRm8GL0uBFz2","trusted":true},"cell_type":"code","source":"chi2_missingness.sort_values('p_value', inplace=True)\nreject, pvalue_corrected = fdrcorrection(pvals=chi2_missingness.p_value, alpha=0.05, method='indep', is_sorted=True)\nchi2_missingness['p_value_corrected'] = pvalue_corrected\nchi2_missingness['BH_reject'] = chi2_missingness['Cochran_criterion'] & reject\nchi2_missingness.to_csv('chi2_tests_missingness')","execution_count":null,"outputs":[]},{"metadata":{"id":"f_sM8JvMBFz6"},"cell_type":"markdown","source":"Let's take a look at the dataframe."},{"metadata":{"id":"gmVSVx_FB9Z7","outputId":"e08573f0-ca28-435c-e1d8-7dbc6e72d934","trusted":true},"cell_type":"code","source":"chi2_missingness.head()","execution_count":null,"outputs":[]},{"metadata":{"id":"3c3ubbgMBF0E"},"cell_type":"markdown","source":"There is no significant test."},{"metadata":{"id":"bEmGlBZM5jib","outputId":"8c6e8b4b-5099-4d87-d5b3-ddefee056e87","trusted":true},"cell_type":"code","source":"significant_tests = chi2_missingness[chi2_missingness['BH_reject']==True]\nprint('Number of statistically significant tests: ', significant_tests.shape[0])","execution_count":null,"outputs":[]},{"metadata":{"id":"o9K9wa3mRDX5"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; It could be tempting to conclude no variable is MAR. Nonetheless, we just looked at interactions between couples of variables. Maybe the missingness of a variable is linked to the observed values of several other variables.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; So we could introduce interactions terms in the dataframe and then perform again chi-squared tests. That would need to carry out numerous tests.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; A more pragmatic approach is to try to predict the missingness of a variable from the others. Then we should support the data is MAR if and only if we achieve good results. But if we don't, we can't say if it is because the data is MCAR or MNAR. Indeed, whether the missingness of a variable is random or mainly depending on unobserved data or on itself, it is unpredictable.\n<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; But before doing it, let's clarify some points. Everytime we tested the independence between a variable and the missingness indicator of another variable, the rows containing missing data for the first variable were dropped. We could have added a category for the missingness in every variable before the tests. In this way, we would have considered the relationship between the missingness of different variables.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; However, if the missingness indicator of a variable is linked to the observed values of a second variable but not to the missingness of that second variable, then the significance of the test would be reduced by adding a new category. Indeed, by adding a missing category, we add degrees of freedom to the chi-squared law asymptotically followed by the $T$ statistic. But if the added category is not linked to the missingness of the other variable, then the $T$ statistic is barely higher. Thus, the test would have less propensity to be significant.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; One solution could be to test the independence between the missingness indicators of two different variables. In this way, the exploration of the joint missingness dependence would not affect the study of the dependence between the missingness of a variable and the observed values of another.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; If you have spare time, you can try it. Maybe it can reveal something that I missed. I did not want to spend too much time on this.<br>\n<br>\nAfter this clarification, let's move on to the prediction of the missingness.<br>\n<br>\n\n### II.5 Relationships between more than 2 features.<a id='II.5'></a><br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; Below, we create a function to forecast the missingness of each feature from the others.<br>\nFor this process, to use the target would be a case of data leakage as the target is unavailable for the test set (so we could not impute the NaN of the test set this way). Nonetheless, we use the entire dataset (training set plus test set) for the same reason we did it for the chi-squared tests.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; We will use a tree-based gradient boosting model to take into account the interactions between the predictors. The features are first target encoded (with a 5-fold strategy for regularization). As we are looking for dependencies more than the perfect model, here we use a quick and simple encoding (just to have numerical values to feed to the model) with a model whose hyperparameters are not optimized but just set to classical values.<br>"},{"metadata":{"id":"bFbflxtsX73n","trusted":true},"cell_type":"code","source":"def predict_missingness(col):\n\n  X = df_total.drop(col, axis=1)\n  y = df_total[col].isna()\n  cols = X.columns.tolist()\n\n  X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, stratify=y, random_state=0)\n  X_train.reset_index(drop=True, inplace=True)\n  y_train.reset_index(drop=True, inplace=True)\n  X_train_te = X_train.copy()\n\n  # Target encoding\n  for index_fit, index_transform in StratifiedKFold(n_splits=5, random_state=0, shuffle=True).split(X_train, y_train):\n    target_enc = TargetEncoder(cols = cols, smoothing=0.5, min_samples_leaf=5)\n    target_enc.fit(X_train.iloc[index_fit,:], y_train.iloc[index_fit])\n    X_train_te.loc[index_transform,:] = target_enc.transform(X_train.iloc[index_transform, :])\n    \n  X_train_te = X_train_te.astype(float)\n\n  target_enc = TargetEncoder(cols = cols, smoothing=0.5, min_samples_leaf=5)\n  target_enc.fit(X_train, y_train)\n  X_val_te = target_enc.transform(X_val)\n\n  weight = sum(y_train==0)/sum(y_train==1)\n\n  # Gradient boosted trees\n  gbt = xgb.XGBClassifier(n_estimators=3000, random_state=0, objective='binary:logistic', scale_pos_weight=weight,\n                          learning_rate=0.2, max_depth=4, subsample=0.8, colsample_bytree = 0.8, min_child_weight=1000)\n  \n  gbt.fit(X_train_te, y_train, eval_set= [(X_train_te, y_train), (X_val_te, y_val)], eval_metric=['auc'], early_stopping_rounds=50, \n          verbose=False)\n\n  # Performance\n  y_pred = gbt.predict_proba(X_train_te)\n  train_score = roc_auc_score(y_train, y_pred[:,1])\n  y_pred = gbt.predict_proba(X_val_te)\n  val_score = roc_auc_score(y_val, y_pred[:,1])\n\n  return(train_score, val_score)","execution_count":null,"outputs":[]},{"metadata":{"id":"EemkIZxg5EwW","outputId":"f2c622ec-f54f-42f8-f544-60578aec14a6","trusted":true},"cell_type":"code","source":"start = datetime.now()\n\npred_missingness = pd.DataFrame(columns=['Feature', 'train_score', 'val_score'])\n\nfor col in df_total.columns:\n  train_score, val_score = predict_missingness(col)\n  pred_missingness.loc[col] = [col, train_score, val_score]\n\npred_missingness.sort_values('val_score', ascending=False, inplace=True)\npred_missingness.to_csv('pred_missingness.csv', index=False)\n\nruntime = datetime.now() - start\nprint('Runtime ', runtime)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_missingness.head()","execution_count":null,"outputs":[]},{"metadata":{"id":"FjTLEbsjOzRi"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; The results are very poor. We could have the same performances by imputing values randomly. That confirms the results of the tests.<br>\nEven by optimizing the classifier for this task, it is unlikely that the validation performances (the area under the ROC curve) increase from 0.5 to 0.75 (even the training performance is weak).\nWe should not obtain such poor results under the MAR assumption.<br> Consequently, we effectively answered the question we had asked ourselves: the MAR assumption is dismissed for all the features. So we do not conclude about their missingness mechanism as we cannot actually distinguish the MCAR from the MNAR. <br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; We can try to predict the values taken by each feature with a model trained on the others. Let's assume we achieve this well for a feature. Then if the MCAR assumption holds for this feature, the prediction of the actual values of the NaN should be as efficient as the predictions of the non-missing values (by definition of the MCAR mechanism). So we would then just impute the NaN with their predicted values and see if this provides better results for the prediction of the target than other imputation methods. If the MCAR assumption holds, then it should be the case. If the missingness depends too strongly on unavailable information, then it should not be the case.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Let us try it with a naive bayesian classifier, for the same reasons as previously: it is quick and simple.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; This classifier just uses the Bayes theorem (making the naive assumption the features are conditionally independent given the target) to estimate the probability of belonging to each class given the values taken by the categorical features.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; The scikit-learn's version of this model needs the input data to be ordinal encoded. So we first apply an ordinal encoding to the input. Nonetheless each ordinal value is actually interpreted as a modality of a categorical feature, not as a number.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; The alpha hyperparameter of the naive bayes classifier is a smoothing parameter which makes the calculated probabilities closer to equiprobability (the case where the probability of pertaining to each category of a given feature is the same, which can be thought of our a priori). It enables to avoid the classifier to estimate the probabilities by the exact frequencies in the training set, as it could lead to overfitting and a bad generalization of the observed patterns. Here, the default value of 1 is used.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; This model does not actually consider the interactions because of the assumption of conditional independence of the features. Instead, it measures an aggregation (precisely the product) of the individual impacts of each feature.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; For the prediction of each feature, the model's performance is measured as follows. When the feature is binary, it is measured by the area under the ROC curve (a.k.a AUROC). When the feature contains more than two classes, a ROC curve is computed for each class (by considering all the other classes as a single class, \"the rest\"). Then the performance is measured by the average of all the AUROC computed.<br>\nThe prediction of the high cardinality features is set aside as they are not well suited to this model. We will not try to predict them."},{"metadata":{"id":"YivkzfZEFoAb","trusted":true},"cell_type":"code","source":"def predict_feature(col):\n\n  df_new = df_total[df_total[col].isna()==False].astype(str) # We remove the NaN of the variable to predict.\n\n  X = df_new.drop(col, axis=1)\n  y = df_new[col]\n\n  le = LabelEncoder() # The label encoding of the variable that we want to predict is needed to use this classifier.\n  y = le.fit_transform(y)\n\n  train_scores = []\n  val_scores = []\n\n  for index_train, index_val in StratifiedKFold(n_splits=5, random_state=42, shuffle=True).split(X, y):\n          \n      X_train = X.iloc[index_train]\n      y_train = y[index_train]\n      X_valid = X.iloc[index_val]\n      y_valid = y[index_val]\n\n      # Ordinal encoding.\n      oe = OrdinalEncoder()\n      X_train_oe = oe.fit_transform(X_train.astype(str))\n      X_val_oe = oe.transform(X_valid.astype(str))\n\n      X_train_oe += 2 # we add 2 to avoid negative values resulting from the ordinal encoding \n      X_val_oe += 2   # (the classifier only accept positive values)\n\n      # Training.\n      nb = CategoricalNB()\n      nb.fit(X_train_oe, y_train)\n\n      # We stock the results.\n      if np.max(y)==1:  \n        y_pred = nb.predict_proba(X_train_oe)\n        train_score = roc_auc_score(y_train, y_pred[:,1])\n        train_scores.append(train_score)\n\n        y_pred = nb.predict_proba(X_val_oe)\n        val_score = roc_auc_score(y_valid, y_pred[:,1])\n        val_scores.append(val_score)\n      else:\n        y_pred = nb.predict_proba(X_train_oe)\n        train_score = roc_auc_score(y_train, y_pred, multi_class='ovr') \n        train_scores.append(train_score)                                \n\n        y_pred = nb.predict_proba(X_val_oe)\n        val_score = roc_auc_score(y_valid, y_pred, multi_class='ovr')\n        val_scores.append(val_score)\n\n  return(train_scores, val_scores)","execution_count":null,"outputs":[]},{"metadata":{"id":"HGeBRU6AQB38","outputId":"5e1105ee-4804-4d89-9e4f-c9d65ad3cb0f","trusted":true},"cell_type":"code","source":"start = datetime.now()\n\npred_features = pd.DataFrame(columns=['Feature', 'train_score', 'val_score', 'std_val_score'])\n\nhigh_cardinality_features = ['nom_' + str(i) for i in range(5,10)] + ['ord_5']\nfor col in df_total.columns.drop(high_cardinality_features):\n  train_score, val_score = predict_feature(col)\n  pred_features.loc[col] = [col, np.mean(train_score), np.mean(val_score), np.std(val_score)]\n\npred_features.sort_values('val_score', ascending=False, inplace=True)\npred_features.to_csv('pred_features.csv', index=False)\n\nruntime = datetime.now() - start\nprint('Runtime ', runtime)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_features.head()","execution_count":null,"outputs":[]},{"metadata":{"id":"LCIZI942Sevu"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; The results are as poor as for the prediction of the missingness (close to a random prediction). It would be unlikely to obtain such poor results if an efficient prediction were possible (even if the models were not optimized for these predictions). So it seems to be a bad idea to use models to try to predict the actual values of the NaN.<br>\n<br>\n\n### II.6 Strategy adopted.<a id='II.6'></a><br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; Some missing data imputaion method should not be used when the MCAR asumption is not satisfied. For instance, we should not use listwise deletion (a.k.a complete-case analysis) since it would bias the dataset (in addition to the fact this approach drops information and shorten the dataset, making the model less efficient).<br>\nSimilarly, when the data is not MCAR, we are not justified to think the mode of the observed values is also the mode of the NaN (even if it were the case, the mode imputation would still bias the distribution of the values toward the mode). <br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; If the missingess depends on unavailable information, to impute the NaN would result in loosing this information that we have no other way to get. So to create a new category for the NaN is prone to be efficient by bringing new information.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; We can see the mode imputation as a regularization technique compared to the creation of a new category. In fact, the mode being the most common category of a feature, to learn an observation has taken the mode for a given variable does not gives us much information. On the other hand, to know an observation has taken the less represented value is likely to have much more impact on the prediction (as an illustration, in NLP, we sometimes take advantage of this fact by giving less importance to the most frequent words of a text).<br>\nThus, any model will associate a middle target value to an observation taking the mode (other features being equal). But when a category is created for the missingness, it can possibly be associated with a more discriminant value. When the missingness category is poorly represented (what is generally the case), that discriminant value can be a cause of overfitting as the model has few cases to learn the impact of the missingness. Moreover, the missing values can be very different from each other, what makes it difficult to seize the effect of the missingness.\nConsequently, to impute the missing values with the mode can help to regularize a model.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Besides, one other thing to consider is the number of caterogies of each feature (a.k.a the cardinality). For instance, the bin_4 variable takes the value N for 312 444 rows, the value Y for 269 609 rows and has 18047 missing values."},{"metadata":{"id":"G-fSTVYNQeix","outputId":"fd9784dd-dec5-43d5-f8b1-9321853eb7d8","trusted":true},"cell_type":"code","source":"df.bin_4.fillna('NaN').value_counts()","execution_count":null,"outputs":[]},{"metadata":{"id":"ZAS0LLdQNVX2"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; Therefore, to consider the missingness as a new category of the bin_4 feature would create a category with just 18 047 observations. But if we consider the missingness of nom_5 as a new category, that new category would be its most common category with 17 778 observations against only 977 for its second most common category.\n"},{"metadata":{"id":"31acytyTLLNP","outputId":"0c0ed1fb-208c-4738-f04a-40e6bdf65d8f","trusted":true},"cell_type":"code","source":"df.nom_5.fillna('NaN').value_counts()","execution_count":null,"outputs":[]},{"metadata":{"id":"JMjhUtwnTkgx"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; Thus, to impute a non-representative value to the NaN of nom_5 would bias the most common category, in addition to waste all potential information contained in the missingness of that variable. On the other hand, bin_4 has roughly 15 times less missing values than its number of values in any other category. So the imputation of the missing values of bin_4 would probably have less impact on bin_4 than imputation of the ones of nom_5 would have on nom_5.<br>\n<br>\nThe missingness will be considered as a new category for the following predictors: <br>\n*   nom_1, nom_2, nom_3 and nom_4, day and month as the number of elements per category of these variables is rather adapted to their number of missing values.<br>\n*   nom_5, nom_6, nom_7, nom_8 and nom_9 as we can't ignore their most common category.<br>\n<br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; The missingness will be imputed with the mode for the predictors: bin_0, bin_1, bin_2, bin_3, bin_4 and nom_0 as their number of missing values is small in comparison to the number of observed cases in each category.<br>\nTo avoid overfitting, the models we trained needed whether to be very simple (resulting in more bias but less variance) whether to be particularly regularized. Consequently, to perform a mode imputation of some features seemed to be a reasonable choice. That was another argument to the mode imputation we realized.<br>\n<br>\n\nThe imputations of the missing values of ord_3, ord_4 and ord_5 are kept for later."},{"metadata":{"id":"ZPE44ruLZPOk"},"cell_type":"markdown","source":"## III. Encoding, modeling and interpretation.<a id='III'></a> <br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; The winners of Kaggle competitions usually use large ensembles of models. Indeed, when we pool together the predictions of enough different and approximately equally efficient models, the resulting prediction is prone to be better than the prediction of the best single model of the ensemble. The reason is that different models fail differently. Therefore, it is unlikely that the majority of the models of the ensemble fail in the same way on a prediction. Thus, when one model fail, the others will balanced this error and the ensemble will probably still have a good performance (unless the models perform poorly, of course), that is to say, in our case, predict the right class.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Consequently, in order to get competitive results on this challenge, we will use ensembling. But for our ensembling to be useful, we need to produce different models. This can be done by ensembling different kind of models (for example tree-based models with neural networks), or it can be done by training models with different preprocessing of the data. We will do both to optimize our ensembling.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; In this way, we will produce several representations of our data, and then we will train diverse kind of models on them. Here, our goal is not to find the best possible single encoding (that is to say the one that performs best, the one that is the most theoretically justifiable or the most interpretable). We are not looking for the best possible single model either. Instead, we aim to use different encoding to produce different good models and get the most of them.<br>\nImproving the performance of each model is only an instrumental objective. The final objective is to obtain the best possible ranking, which may involve producing models that are individually less efficient but form a better ensemble.<br>\n<br>\n\n### III.1 First strategy.<a id='III.1'></a> <br>\n\n\n&nbsp;&nbsp;&nbsp;&nbsp; As we imputed all their missing values, the features bin_0, bin_1, bin_2, bin_3 and bin_4 are binary variables (hence their names) taking values 0 and 1 (for bin_0, bin_1 and bin_2) or F and T (for bin_3) or N and Y (for bin_4). Therefore, bin_0, bin_1 and bin_2 are kept in this form. The values of bin_3 and bin_4 are just mapped with 0 (for F and N respectively) and 1 (for T and Y respectively) to be in an appropriate format.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; The predictors nom_0, nom_1, nom_2, nom_3 and nom_4 have a priori no order and have a cardinality low enough to be one-hot encoded. So we will one-hot encode them.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; The predictors nom_5, nom_6, nom_7, nom_8, nom_9 have a priori no order either. Albeit their values look like hexadecimal numbers, to encode them this way was not efficient. But, as they have a high cardinality, we will not one-hot encode them (neither binary encode them or base-N encode them). We won't use ordinal encoding as it would create an artificial order in these nominal features. Instead, we will apply a target encoding in order to create an order that has sense: each modality will be encoded with the probability that the target equals 1 given that modality (in fact this probability is regularized by a weighted averaging with the a priori probability the target equals 1, the weights being tuned by a smoothing parameter). <br>\n&nbsp;&nbsp;&nbsp;&nbsp; Nonetheless, target encoding is more prone to overfitting than ordinal encoding as target encoding use information about the variable we want to predict. To avoid overfitting, one approach is to add random noise (and sometimes to exclude the current row’s target when calculating the mean target for a level to reduce the effect of outliers). Instead, here, we use a stratified 5-fold splitting strategy in the same way as for 5-fold cross-validation: we train a target encoder on 80% of the training set and use it to encode the remaining 20%. Thus, with 5 splits we encode all the training set. The test set is then encoded with a target encoder trained on all the training set. <br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Although ord_0 take three ordered values (1, 2 and 3), as we don't know the meaning of that variable, there is no particular reason to think the distance between the three values should be the same. I tried to use some transformations (taking the square, the square root, etc.) of ord_0 but the default order worked better, so we will keep it. <br>\n&nbsp;&nbsp;&nbsp;&nbsp; To impute its missing values with the mode (which is 1) would bias not only its distribution (because it would over-represent the value 1) but also its relation with the target. Indeed, since the test of independence between the missingness of ord_0 and the target was not significant, there is no reason to think the target more often equals to 1 when ord_0 is missing. But, as shown in the graph below, the modality 1 of ord_0 is associated with a higher propensity of the target to equal 0 (as the target is positively correlated to ord_0 and 1 is the minimum of ord_0).<br>"},{"metadata":{"id":"TM3ux0JI6xDD","outputId":"285d7c0b-3cab-4c9b-99ec-969c5ec43e0d","trusted":true},"cell_type":"code","source":"def plot_hist_bar(col, dtype=\"numeric\", dataset=df, xlabels=True, histogram=True):\n\n    \"\"\" Function to plot a barplot and optionally an histogram of one chosen feature from the data.\n\n        Parameters\n        ----------\n        col : String\n            The feature whose the function will create graphics from.\n        dtype : String (default=\"numeric\")\n            If \"str\", then the feature will be treated as a string feature and so will be ordinal encoded before plotting\n            the charts. Otherwise, the feature is considered to be already numeric.\n        dataset : DataFrame (default=df, namely the training set)\n            The dataframe from containing the variables used.\n        xlabels : Boolean (default=True)\n            Whether to keep the xlabels.\n        histogram : Boolean (default=True)\n            Whether to plot an histogram and a barplot or just a barplot.\n        \n        Display\n        -------\n        \n        - Optionally an histogram given the distribution of the feature col.\n        - A barplot providing the target frequency for each value of the feature col.\n        - A dictionary mapping the encoded values on the x axis of the precedent charts to their original values. For all \n        of them, the missing values of the feature col are filled with the value -1 and represented with a red bar.\n        \n        \"\"\"\n\n    if (dtype=='str'):\n      map_categories = {cat:i for i,cat in enumerate(df.dropna()[col].unique())}\n    else:\n      map_categories = {}\n  \n    map_categories[np.nan] = -1\n    dataset = dataset.replace({col: map_categories})\n    dataset.loc[col] = dataset[col].astype(int)\n\n    fig = plt.figure()\n    palette = dict((value,'grey') if (value != -1) else (value,'red') for value in dataset[col].unique())\n\n    if histogram==True:\n\n      plt.subplot(1,3,1)\n      data = dataset[col].value_counts().reset_index()\n      data.columns = ['Values','Count']\n      ax1 = sbn.barplot(x = 'Values', y = 'Count', data = data, palette=palette)\n      plt.title('Histogram')\n      if xlabels==False:\n        frame1 = plt.gca()\n        frame1.axes.xaxis.set_ticklabels([])\n\n      plt.subplot(1,3,3)\n      ax2 = sbn.barplot(x = col, y = 'target', data = dataset, palette=palette)\n      ax2.set(xlabel='Values', ylabel='Target frequency')\n      if xlabels==False:\n        frame2 = plt.gca()\n        frame2.axes.xaxis.set_ticklabels([])\n      plt.title('Barplot')\n\n    else:\n      ax = sbn.barplot(x = col, y = 'target', data = dataset, palette=palette)\n      ax.set(xlabel='Values', ylabel='Target frequency')\n      if xlabels==False:\n        frame = plt.gca()\n        frame.axes.xaxis.set_ticklabels([])\n\n    fig.suptitle(\"Feature: \" + col, fontsize=16)\n    plt.show()\n  \n    if xlabels==True:\n      print('\\033[30m'' Relation between the numbers on the x axis and the categories they represent:', '\\n', map_categories)","execution_count":null,"outputs":[]},{"metadata":{"id":"hH3Ucjj79ZIC","outputId":"7a3785de-54e8-47d4-9d17-b6fa995799b4","trusted":true},"cell_type":"code","source":"plot_hist_bar('ord_0')","execution_count":null,"outputs":[]},{"metadata":{"id":"3Xp8ufaV9GS6"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; Thus, we will imput the missing values with the value 2. That choice corresponds to a median imputation (whose interpretation is close to the one of the mode imputation for categorical variables). This preserves the relationship between ord_0 and the target.<br>\nWe could also imput the missing values with the mean of ord_0 which approximately equals 1.9. This value would also keep the relationship between ord_0 and the target. Here, the mean imputation also enables the tree-based models to consider its missingness. Indeed, with this imputation, a tree can split the data set in order to separate either the samples which take the value 1 from the others, the samples wich take the value 1 or NaN from the others or the samples wich take the value 3 from the others. With a median imputation, the second splitting could not be done by a tree. To add this possibility could led to a better fitting or to an overfitting. My trials have proven better results with the median imputation.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; As their name state it, ord_1 and ord_2 are also ordinal features. A one-hot encoding (as the cardinality is low) would get over the question of the distances between the values at the cost of a higher dimensionality and a potential loss of the effect of the interactions between the features. One more time, my trials gave me better results by using an ordinal encoding with equal distances. <br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; We tried several mappings between alphabet letters and numbers in order to encode ord_3 and ord_4. One has proven better results with the diverse classifiers we trained. It consists in mapping the capital letters A to Z with the numbers 0 to 25 and the lowercase letters from a to z with the numbers 26 to 51. Thus, as shown in the graph below, the target is an increasing function of ord_3 and ord_4.<br>"},{"metadata":{"id":"xiCamOq36eUb","outputId":"f0844e44-3891-444a-b0f6-229288e7d1b7","trusted":true},"cell_type":"code","source":"# Creation of a dictionnary to encode the features.\nreplace_letters =  {j:i for i,j in enumerate(list(string.ascii_uppercase + string.ascii_lowercase))}\n\n# We plot a barplot representing the target frequency of each modality of ord_3 after encoding.\n# For more clarity, the x-axis is made invisible. The red bar corresponds to the NaN.\ncol ='ord_3'\nplot_hist_bar(col, dataset = df.replace({col:replace_letters}), xlabels=False, histogram=False)","execution_count":null,"outputs":[]},{"metadata":{"id":"149Ew4xR9uck","outputId":"bcc854a4-6005-4054-9331-b430584bafc3","trusted":true},"cell_type":"code","source":"# Same chart for ord_4.\ncol = 'ord_4'\nplot_hist_bar(col, dataset = df.replace({col:replace_letters}), xlabels=False, histogram=False)","execution_count":null,"outputs":[]},{"metadata":{"id":"fiXVXbOP-J7P"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; Note that, as the values taken by ord_3 and ord_4 are respectively lowercase letters and uppercase letters, we could have encoded each letter with its rank in the alphabet, without distinguishing both cases. We still treated both separately because we're going to use the same dictionnary replace_letters to encode ord_5.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; The values taken by ord_3 range from 26 to 51, but if they were going from 0 to 25, the results would be inchanged because we will only use this strategy with tree-based methods or after standardization (which means to substract the mean value and then to divide by the standard deviation). Indeed, just the order between the values influence the construction of a tree, so only the order matters for the tree-based methods. Furthermore, to standardize a variable before feeding it to a model makes it invariant by translation.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; This encoding (of ord_3 and ord_4) is similar to a target encoding as the target is an increasing function of the resulting features. However, a target encoding would either add noise (by a 5-fold spliting strategy or by directly adding random noise, as explained before) or grow the risk of overfitting by inducing distances between the values depending exclusively on the target. Here, we found a simple and straight encoding strategy that seems to make sense in view of its interpretation (each letter is associated with its rank in the alphabet and the uppercases correspond to larger numbers than the lowercases). As we did not deduce this order from the target values but for our interpretation of the feature itself, there is no reason to think we will overfit by using it, even if its link with the target is monotone. So adding noise would probably worsen the relationship we found. <br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; The above graphs highlight that imputing the missing values of ord_3 and ord_4 with the median of the encoded values does not break the monotone relationship with the target. So we will handle them this way.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Next, ord_5 is an ordinal predictors whose values are combinations of two letters (in uppercase or lowercase). After having tortured that variable, I have come to conclusion that a relevant way to use it is to just keep the first letter of each value and encode it in the same way that ord_3 and ord_4 (I have nothing personal against the second letter, it is just feature engineering). One more time, the target is an increasing function of the resulting feature."},{"metadata":{"id":"tUoI2-USPJXQ","outputId":"ce8bb0d7-422e-430a-e3c3-2c3ee57c4458","trusted":true},"cell_type":"code","source":"df_new=df.copy()\n# We fill the NaN with 0 and keep only the first character of ord_5.\ndf_new.ord_5 = df_new.ord_5.fillna('0').astype(str).str[0] \ndf_new.ord_5 = df_new.ord_5.replace({'0':np.nan}) # We transform back the 0 to NaN. \n#This trick allowed to conserve the NaN even when keeping only the first character.\nplot_hist_bar(col, dataset = df_new.replace({col:replace_letters}), xlabels=False, histogram=False)","execution_count":null,"outputs":[]},{"metadata":{"id":"wT8tTDnee33M"},"cell_type":"markdown","source":"The missing values are also replaced by the median, for the exact same reason (underscored by the above chart) as for ord_3 and ord_4.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Although the predictors day and month are cyclical, they are provided by Kaggle encoded as ordinal features. <br>\nMonday is the day after Sunday. So, in a sense, Sunday and Monday are closer than Sunday and Wednesday. Nevertheless, as the day variable is numbered from 1 to 7 (to represent the days of the week), Monday equals 1, Wednesday 3 and Sunday 7. In this way, Sunday and Wednesday are closer than Sunday and Monday. The month variable, taking values going from 1 to 12 to number the months of the year, has the same peculiarity.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; To overcome this issue, day and month could be one-hot encoded. In order to restrain the dimensionality, we rather choose to project each of these features in two dimensions using a classical mapping. The idea is simple. The values of the day feature range from 1 to 7. So they are first multiplied by $\\frac{2\\pi}{7}$ in order to belongs to the interval $[0;2\\pi]$. Then, they are wrapped around the unit circle (that is to say the circle centered in zero, the origin of a two-dimensional space, and with radius 1) whose perimeter equals $2\\pi$. Thus, the values don't overlap and are all equidistant on the circle. An illustration is available below."},{"metadata":{"id":"aS3HtfhQ-E5j","outputId":"76495073-6fd8-4245-d4bd-daa7cb1c93b4","trusted":true},"cell_type":"code","source":"def wrap_around_unit_circle(n, with_nan=False, title=None):\n\n    x = np.cos(2*np.pi*np.arange(1,n+1)/n)\n    y = np.sin(2*np.pi*np.arange(1,n+1)/n)\n    \n    unit_circle = plt.Circle((0, 0), 1, color='grey', fill=False)\n    \n    fig, ax = plt.subplots()\n    ax.scatter(x, y, s=10, color='black', marker='x')\n    ax.set_xlim((-1.1,1.1))\n    ax.add_artist(unit_circle)\n    \n    for i in range(1,n+1):\n        if (i==n)&(with_nan==True):\n            ax.annotate(\"NaN\", (x[i-1]+0.02, y[i-1]), color='red')\n        else:\n            ax.annotate(i, (x[i-1]+0.02, y[i-1]), color='black')\n\n    if title!=None:\n        plt.title(title, color='white')\n\nwrap_around_unit_circle(7, title=\"Wrapping of the days around the unit circle\")","execution_count":null,"outputs":[]},{"metadata":{"id":"fdZvOTrXFFBu"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; Thus, closer are days in the calendar, closer they are on the circle. For instance, Sunday and Monday are closer than Sunday and Wednesday. To feed the model with the information of the spatial localisation of each day on that circle, we just need to provide it the coordinates of each day on this two-dimensional space. Fortunately, these coordinates are nothing but the cosinus and the sinus of the days values after having multiplied them by $\\frac{2\\pi}{7}$. <br> \n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; To handle the missing days, we will just refine a little bit this procedure. We will firstly impute the missing days with the value 8. Consequently, next, instead of multiply the values of the day feature by $\\frac{2\\pi}{7}$, we multiply them by $\\frac{2\\pi}{8}$ (as there is now 8 day values). Then, we wrap them around the unit circle. Thus, one more time, the values do not overlap and are all equidistant on the circle. Moreover, as the missing days take the biggest value on the abscisse (wich is 1), a tree can easily split the dataset along the missingness of the day feature. As we will use a tree-based method with this preprocessing, there is no need to create an indicator variable for the missingness of the day feature.<br>\nAn illustration of this procedure is available below."},{"metadata":{"id":"rv-4kYJOFFBv","outputId":"1ec0425a-ebc8-4d6b-b7e9-81e619a1e2d5","trusted":true},"cell_type":"code","source":"wrap_around_unit_circle(8, with_nan=True, title=\"Wrapping of the days (including the missing ones) around the unit circle\")","execution_count":null,"outputs":[]},{"metadata":{"id":"zqx7REqLFFB2"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; The same rationale applies to the month feature. We therefore encode it by firstly imputing the missing months with the value 13, then multiplying all its values by $\\frac{2\\pi}{13}$ and then replacing it with its cosinus and its sinus.<br>\nWe illustrate the result below."},{"metadata":{"id":"_FYmh0CQFFB4","outputId":"c431f472-0017-4ad8-c551-b59f9a79653e","trusted":true},"cell_type":"code","source":"wrap_around_unit_circle(13, with_nan=True, title=\"Wrapping of the months (including the missing ones) around the unit circle\")","execution_count":null,"outputs":[]},{"metadata":{"id":"lNZsIuVXFFB9"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; Brief remark. When using one-hot encoding, we create a dummy (binary) variable for each modality of the encoded feature. Hence, we introduce multicollinearity as we can deduce the belonging to one category from the knowledge of the belonging to the others. <br>\nTo be rigorous, it is not the case if we did not create a dummy variable for the missing or unknown values (unknown values are values present in a test set and not in the training set). Thus, an observation with a missing or unknown value would take zero for every dummy variable.<br>\nHere, as we created dummy variables to represent the missing values, our dummy variables are in fact multicolinear. That issue (also known as the dummy variable trap) is less problematic with tree-based models (what we will utilize) than with others. Nevertheless, to be rigorous, after each one-hot encoding, we will delete one dummy variable. In particular, every time we will delete the one associated with the most common modality.<br>\nThe choice of the deleted dummy variable impact the interpretation. With the choice we made, the samples taking 0 for all the dummy variables associated with a feature are the ones taking its mode. We can then interpret them as the \"normal\" samples as they are the most numerous.<br>\n<br>\nLet's define a function to preprocess the data."},{"metadata":{"id":"-yNjL4FJVc8O","trusted":true},"cell_type":"code","source":"def preprocessing_1(df_train, df_test):\n\n  \"\"\" Function to preprocess the data.\n\n        Parameters\n        ----------\n        df_train : DataFrame, shape = (600000, 24)\n            Training set to preprocess, containing the features plus the target.\n        df_test : DataFrame, shape = (400000, 23)\n            Test set to preprocess.\n        \n        Returns\n        -------\n        X_train_encoded :  DataFrame, shape = (600000, 44)\n            Training set preprocessed, without the target variable.\n        y_train :  Series, shape = (600000,)\n            Training's target variable.\n        X_test_encoded :  DataFrame, shape = (400000, 44)\n            Test set preprocessed.\n\n        \"\"\"\n\n  # Extraction of the target from the training set.    \n  y_train = df_train['target'].copy()\n  df = pd.concat([df_train.drop('target', axis=1), df_test])\n\n  # Filling of the missing values.\n  df = df.fillna(value={'bin_0': df.bin_0.mode()[0], 'bin_1':df.bin_1.mode()[0], 'bin_2':df.bin_2.mode()[0], 'bin_3':df.bin_3.mode()[0], \n                        'bin_4':df.bin_4.mode()[0], 'nom_0':df.nom_0.mode()[0], 'nom_1':'nan', 'nom_2':'nan', 'nom_3':'nan', 'nom_4':'nan',\n                        'nom_5':'nan', 'nom_6':'nan', 'nom_7':'nan', 'nom_8':'nan', 'nom_9':'nan', 'ord_0':df.ord_0.median(), 'ord_5':'0',\n                        'day':8, 'month':13})\n\n  # Dropping of the second letter of ord_5.\n  df['ord_5'] = df['ord_5'].astype(str).str[0] # If we had not filled the NaN with zero, they would have resulted in the \"n\" value when taking the first letter.\n  \n  # Creation of a dictionnary for replacing the letters of ord_3, ord_4 and ord_5.\n  replace_letters = {j:i for i,j in enumerate(list(string.ascii_uppercase + string.ascii_lowercase))}\n  replace_letters['0'] = np.nan # This trick enables to later fill the NaN with the median of the encoded non-missing data.\n\n  # Encoding of bin_3, bin_4, ord_1, ord_2, ord_3, ord_4 and ord_5.\n  df = df.replace({'bin_3':{'F':0, 'T':1}, 'bin_4':{'N':0, 'Y':1},'ord_1': {'Novice':1, 'Contributor':2, 'Expert':3, 'Master':4, 'Grandmaster':5},\n                   'ord_2':{'Freezing': 1, 'Cold':2, 'Warm':3, 'Hot':4, 'Boiling Hot': 5, 'Lava Hot':6},'ord_3': replace_letters, 'ord_4':replace_letters,\n                   'ord_5':replace_letters})\n  \n  # The missing values of ord_1, ord_2, ord_3, ord_4 and ord_5 are imputed with the median of the values after encoding. \n  df = df.fillna(value={'ord_1':df['ord_1'].median(), 'ord_2':df['ord_2'].median(), 'ord_3':df['ord_3'].median(), 'ord_4':df['ord_4'].median(), 'ord_5':df['ord_5'].median()})\n\n  # One hot encoding.\n  one_hot_cols = ['nom_0','nom_1','nom_2','nom_3','nom_4']\n  X_cat = df[one_hot_cols].astype(str).copy()\n  one_hot_enc = OneHotEncoder()\n  X_cat_1hot = one_hot_enc.fit_transform(X_cat)\n  X_cat_1hot = pd.DataFrame(X_cat_1hot.toarray(), index=X_cat.index, columns = one_hot_enc.get_feature_names())\n  X_cat_1hot.drop(['x'+str(i)+'_' + X_cat['nom_' + str(i)].mode()[0] for i in range(5)], axis=1, inplace=True)\n  X = pd.concat([df.drop(one_hot_cols, axis=1), X_cat_1hot], axis=1)\n\n  # Two-dimensional projection of day and month.\n  X['sin_day'] = np.sin(2*np.pi*X.day/8)\n  X['cos_day'] = np.cos(2*np.pi*X.day/8)\n  X['sin_month'] = np.sin(2*np.pi*X.month/13)\n  X['cos_month'] = np.cos(2*np.pi*X.month/13)\n  X.drop(['day', 'month'], axis=1, inplace=True)\n\n  # To apply a target encoding, we need to separate the training set from the test set.\n  X_train = X.iloc[:600000]\n  X_test = X.iloc[600000:1000000]\n\n  # Target encoding.\n  target_enc_cols = ['nom_5', 'nom_6', 'nom_7', 'nom_8', 'nom_9']\n  X_train_encoded = X_train.copy()\n\n  for index_fit, index_transform in StratifiedKFold(n_splits=5, random_state=42, shuffle=True).split(X_train, y_train):\n    target_enc = TargetEncoder(cols = target_enc_cols, smoothing=0.5)\n    target_enc.fit(X_train.iloc[index_fit,:], y_train.iloc[index_fit])\n    X_train_encoded.loc[index_transform,:] = target_enc.transform(X_train.iloc[index_transform, :])\n  \n  target_enc = TargetEncoder(cols = target_enc_cols, smoothing=0.5)\n  target_enc.fit(X_train, y_train)\n  X_test_encoded = target_enc.transform(X_test)\n\n  # Data types management.\n  X_train_encoded.loc[:,target_enc_cols] = X_train_encoded[target_enc_cols].astype(float)\n  X_train_encoded.loc[:, X_train_encoded.columns.drop(target_enc_cols).to_list()] = X_train_encoded.loc[:, X_train_encoded.columns.drop(target_enc_cols).to_list()].astype(float)\n  X_test_encoded.loc[:,'ord_3'] = X_test_encoded.loc[:,'ord_3'].astype(float)\n  X_test_encoded.loc[:,'ord_4'] = X_test_encoded.loc[:,'ord_4'].astype(float)\n  \n  return(X_train_encoded, y_train, X_test_encoded)","execution_count":null,"outputs":[]},{"metadata":{"id":"Hy06QoNBt84G"},"cell_type":"markdown","source":"Let's use it to produce a first preprocessing of the data."},{"metadata":{"id":"7oHTSXlm9kiR","outputId":"1bbefdba-3de9-4f5d-e460-b0f92dc65851","trusted":true},"cell_type":"code","source":"start = datetime.now()\nX, y, X_test = preprocessing_1(df, test)\nend = datetime.now()\nprint('Duration: ', end-start)","execution_count":null,"outputs":[]},{"metadata":{"id":"TJFe2qFK9kiY","outputId":"c0740973-ac87-4ea0-e193-1942582edfbc","scrolled":true,"trusted":true},"cell_type":"code","source":"print('X shape :', X.shape)\nprint('y shape :', y.shape)\nprint('X_test shape :', X_test.shape)","execution_count":null,"outputs":[]},{"metadata":{"id":"BRE9sMcVt84W"},"cell_type":"markdown","source":"Quick visualization of the resulting training set."},{"metadata":{"id":"2V20brsy9kic","outputId":"b2731113-e0d2-4e6f-b8b7-3e390e8d63d4","trusted":true},"cell_type":"code","source":"X.head().T","execution_count":null,"outputs":[]},{"metadata":{"id":"ralGLsQ3RiP1","outputId":"04da53cb-ffc7-47ae-acf8-da5ad2cfc5e3","trusted":true},"cell_type":"code","source":"# Saving of the dataframes.\nX.to_csv('X_1st_strategy.csv', index=False)\ny.to_csv('y.csv', index=False, header=True)\nX_test.to_csv('X_test_1st_strategy.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"id":"r62ZSQK1t84x"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; We will train a gradient boosted random forest on this dataframe. For that, we will use the XGBClassifier from the library XGBoost. Although it is more common to present firstly the simplest model (as a logistic regression), we begin with this one as it is trained on the preprocessing strategy that seems the more natural to me (given the previous analysis). <br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Gradient boosted random forests are similar to gradient boosted trees. The difference is that at each iteration, a random forest is trained instead of a single tree. Here, I first trained and tuned gradient boosted trees. Once I got a good performance, I tried to add trees at each iteration to see if it was worthy. That's why I used just 5 trees by iteration. My hyperparameters were already tuned for gradient boosted trees. To add more trees may provide a better performance if the hyperparameters are modified consequently.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; This classifier can take many hyperparameters. The values we used has been choosen by cross-validation with a gridsearch strategy. The hyperparameter values could probably still be optimized. Indeed, I use Kaggle to practice data science on free data and to exchange with the great Kaggle community, not to run hyperparameter search during hours. I completely agree with the words of François Chollet in his book Deep Learning With Python (Manning): \"it shouldn't be your job as a human to fiddle with hyperparemeters all day -- that is better left to a machine\".<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; The hyperparameter values selected are quite classical, except for colsample_bytree. <br>\nThe colsample_bytree hyperparameter determines the proportion of the features that will randomly be choosen for the construction of each tree. It is set to 0.2. Usually, the used values are larger. This little value allows to regularize the model.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Let's define two functions to plot ROC curves and precision and recall curves. Both have been taken (and slightly modified) from the book Hands-On Machine Learning with Scikit-Learn and Tensorflow by Aurélien Géron (O'Reilly)."},{"metadata":{"id":"Dopu_B1Ftfcx","trusted":true},"cell_type":"code","source":"# Function to plot the ROC curve of all predictions on validation sets.\ndef plot_roc_curve(fpr, tpr):\n    plt.plot(fpr, tpr, linewidth=2)\n    plt.plot([0,1],[0,1],'k--')\n    plt.axis([0,1,0,1])\n    plt.xlabel('Taux de faux positifs')\n    plt.ylabel('Taux de vrais positifs')\n    plt.suptitle('ROC curve of all predictions on validation sets', color=\"white\")\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"id":"d9m39ug4tfc0","trusted":true},"cell_type":"code","source":"# Precision and recall as a function of the decision threshold.\ndef plot_precision_recall_vs_threshold(y, y_pred):\n    precisions, recalls, thresholds = precision_recall_curve(y, y_pred)\n    plt.plot(thresholds, precisions[:-1],'b--', label='Precision')\n    plt.plot(thresholds, recalls[:-1], 'g-', label='Recall')\n    plt.xlabel('Thresholds')\n    plt.legend(loc='center left')\n    plt.ylim([0,1])\n    plt.suptitle('Precision and recall versus thresholds of all predictions on validation sets', color=\"white\")\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"id":"axTQtuTc0RVl"},"cell_type":"markdown","source":"Let's define another function to display the results of the cross-validation."},{"metadata":{"id":"ycv3MaRg0RVm","trusted":true},"cell_type":"code","source":"def display_scores(train_scores, val_scores, y_validations, y_predictions):\n\n    # Printing of the scores.\n    print('Training scores: ', train_scores)\n    print('Mean training score: ', np.mean(train_scores))\n    print('Standard deviation of the training scores: ', np.std(train_scores), '\\n')\n\n    print('Validation scores: ', val_scores)\n    print('Mean validation score: ', np.mean(val_scores))\n    print('Standard deviation of the validation scores: ', np.std(val_scores), '\\n\\n')\n\n    # Precision-Recall versus decision thresholds.\n    y_valid = np.concatenate(tuple(y_validations))\n    y_pred = np.concatenate(tuple(y_predictions))\n\n    plot_precision_recall_vs_threshold(y_valid, y_pred)\n\n    print('\\n')\n\n    # ROC curve.\n    fpr, tpr, threshold = roc_curve(y_valid, y_pred)\n    plot_roc_curve(fpr,tpr)","execution_count":null,"outputs":[]},{"metadata":{"id":"T8bFOulLcoDK","outputId":"ac431762-9ee2-4d99-af2b-ff6953ede5d5","trusted":true},"cell_type":"code","source":"start = datetime.now()\n\ntrain_scores = []\nval_scores = []\ny_validations = []\ny_predictions = []\n\nfor index_train, index_val in StratifiedKFold(n_splits=5, random_state=42, shuffle=True).split(X, y):\n        \n    X_train = X.iloc[index_train]\n    y_train = y.iloc[index_train]\n    X_valid = X.iloc[index_val]\n    y_valid = y.iloc[index_val]\n\n    weights = sum(y_train.values==0)/sum(y_train.values==1)\n\n    \n    xgb_gbrf = xgb.XGBClassifier(n_estimators=3000, random_state=0, objective='binary:logistic', scale_pos_weight=weights, \n                                    learning_rate=0.15, max_depth=2, subsample=0.7, min_child_weight=500,  \n                                    colsample_bytree = 0.2, reg_lambda = 3.5, reg_alpha=1.5, num_parallel_tree= 5)\n    \n    xgb_gbrf.fit(X_train, y_train, eval_set= [(X_train, y_train), (X_valid, y_valid)], eval_metric=['auc'], \n                    early_stopping_rounds=100, verbose=False)\n    \n    # We stock the results.\n    y_pred = xgb_gbrf.predict_proba(X_train)\n    train_score = roc_auc_score(y_train, y_pred[:,1])\n    train_scores.append(train_score)\n\n    y_pred = xgb_gbrf.predict_proba(X_valid)\n    val_score = roc_auc_score(y_valid, y_pred[:,1])\n    val_scores.append(val_score)\n    \n    y_validations.append(y_valid)\n    y_predictions.append(y_pred[:,1])\n    \n    \ndisplay_scores(train_scores, val_scores, y_validations, y_predictions)\n\nruntime = datetime.now() - start\nprint('\\n\\nRuntime ', runtime)","execution_count":null,"outputs":[]},{"metadata":{"id":"EqLex66Wt85T"},"cell_type":"markdown","source":"This classifier got a mean area under the ROC curve equals to 0.7880 on the 5 fold cross-validation.<br>\n<br>\nXGBoost offers three ways of measuring the importance of the features. The three kinds of importance are the weight (or frequency), the cover and the gain.<br>\n<br>\nThe weight importance type returns, for each features, the number of times it is used in all the trees of the model, as a percentage of the total sum of this value across all the features."},{"metadata":{"id":"cYc8UIw06fqN","outputId":"8960a498-9bed-40c5-c125-0bc563cac457","trusted":true},"cell_type":"code","source":"xgb.plot_importance(xgb_gbrf, importance_type='weight', max_num_features=15)","execution_count":null,"outputs":[]},{"metadata":{"id":"ttUoXwXlE_Ct"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; Consequently, the predictors with the highest cardinality tend to have the largest importance as measured by the weight. Indeed, higher is the cardinality of a variable, larger the number of times we need to split it in order to extract its whole information.<br> \n&nbsp;&nbsp;&nbsp;&nbsp; To illustrate this fact, let's assume we want to separate the samples according to the modality they take for a particular feature. If the feature is binary, the model just need to split the samples once along this feature. But if the feature takes for example 100 values, the model needs to split the dataset 99 times along this feature. So the predictors with a high cardinality are expected to be used more often than the ones with a low cardinality. Hence the observed results.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; The cover importance type returns, for each feature, an indicator of the number of observations it splitted, that is to say an indicator of the number of observations whose the decision of the leaf node is made using this feature (for all trees)."},{"metadata":{"id":"--0rTdRb7Yci","outputId":"e032f466-e438-4b4c-ac2b-a021c57d38b9","trusted":true},"cell_type":"code","source":"xgb.plot_importance(xgb_gbrf, importance_type='cover', max_num_features=15)","execution_count":null,"outputs":[]},{"metadata":{"id":"VEokmXDLLbtA"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; Here, the predictors with a lower cardinality tend to have higher importance. Indeed, the fact the binary features have higher values for this importance type is coherent as a single split along them is sufficient to extract all the information they include. \n<br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; The gain importance type returns, for each features, the reduction of the loss function brought by the splits along this feature. So, more than the other importance types, it actually indicates the contribution of each feature to the prediction."},{"metadata":{"id":"PVDspxLj27lJ","outputId":"750a1bc4-d785-4ebe-9b0c-2f46df4cb577","trusted":true},"cell_type":"code","source":"xgb.plot_importance(xgb_gbrf, importance_type='gain', max_num_features=15)","execution_count":null,"outputs":[]},{"metadata":{"id":"KzRhfoM19jt1"},"cell_type":"markdown","source":"Below is the complete ranking of the feature importances measured by the gain (in the below table, it is normalized so that the sum of all the gains equals 1)."},{"metadata":{"id":"DFj5a2tutwOz","trusted":true},"cell_type":"code","source":"gbrf_feature_importances = pd.DataFrame({'Predictors': X.columns, 'Normalized gains': xgb_gbrf.feature_importances_}) \ngbrf_feature_importances.sort_values('Normalized gains',ascending=False, inplace=True)\ngbrf_feature_importances.reset_index(inplace=True, drop=True)\ngbrf_feature_importances.to_csv('fi_gbrf_1st_strategy')","execution_count":null,"outputs":[]},{"metadata":{"id":"CMm2gmdNtwO7","outputId":"21fcd729-3803-4a29-d55f-b6972fc4e63b","trusted":true},"cell_type":"code","source":"gbrf_feature_importances","execution_count":null,"outputs":[]},{"metadata":{"id":"cUug0wFr9jt7"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; Although nom_7 has been the most used feature and x1_Trapezoid the one wich impacted the classification of the greatest number of samples, it is ord_3 which allowed the largest reduction of the loss function. Generally, the ordinal featured were the most important. The target encoded and the one-hot encoded features were not so discriminant.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; To be rigorous, the importances given by the gains must not be taken literally. In fact, the classifier has been trained with the log loss (a.k.a logistic loss or cross-entropy loss). So the gain importance type has been computed as the reduction of the log loss allowed by each feature on the training set. But the metric we actually want to optimize is the area under the roc curve on the test set, and the ranking of the features according to their impact on the optimization of this metric may be different from their ranking on the precedent table. However, this effect should be slight. <br>\n<br>    \n&nbsp;&nbsp;&nbsp;&nbsp; Now let's train this classifier on the whole training set. For this, we will use 818 iterations as it is the average best number of iterations on the cross-validation (the median were 759). The number of iterations is provided when the verbosity of the fit method of the classifier is set to True. Here, I set it to False to shorten the notebook."},{"metadata":{"id":"OSozbRQ50ZQR","outputId":"7656d256-de7d-4c17-ade4-cecb29108575","scrolled":true,"trusted":true},"cell_type":"code","source":"np.mean([749,738,759,943,901])","execution_count":null,"outputs":[]},{"metadata":{"id":"23XXvOT6IXXQ","trusted":true},"cell_type":"code","source":"#### Modelisation with XGBoost Gradient Boosting trees Classifier.\n\nweight = sum(y.values==0)/sum(y.values==1)\n\n\nxgb_gbrf = xgb.XGBClassifier(n_estimators=818, random_state=0, objective='binary:logistic', scale_pos_weight=weight, \n                             learning_rate=0.15, max_depth=2, subsample=0.7, min_child_weight=500,  colsample_bytree = 0.2,\n                             reg_lambda = 3.5, reg_alpha=1.5, num_parallel_tree= 5)","execution_count":null,"outputs":[]},{"metadata":{"id":"eUm5JfihfyKZ","outputId":"c88ce711-7ffa-4142-a444-25e981443700","trusted":true},"cell_type":"code","source":"xgb_gbrf.fit(X, y, eval_set= [(X, y)], eval_metric=['auc'], verbose=False)","execution_count":null,"outputs":[]},{"metadata":{"id":"MindLp0LUqkx"},"cell_type":"markdown","source":"Saving of the model."},{"metadata":{"id":"J3QGBuhxGztw","outputId":"0832b91a-45a6-4ab6-d19b-9eecdfd68fe9","trusted":true},"cell_type":"code","source":"joblib.dump(xgb_gbrf, 'model_gbrf_1st_strategy')","execution_count":null,"outputs":[]},{"metadata":{"id":"-Muimkh8UvHL"},"cell_type":"markdown","source":"We can now use the model to perform a prediction on the test set."},{"metadata":{"id":"uUDRatxlizau","outputId":"a14a85ea-af38-4ad0-bde7-db7f961e5447","trusted":true},"cell_type":"code","source":"y_test = xgb_gbrf.predict_proba(X_test)[:,1]\noutput = pd.DataFrame({'id':sample_submission.index, 'target':y_test})\noutput","execution_count":null,"outputs":[]},{"metadata":{"id":"OvJkaoeLlHLI"},"cell_type":"markdown","source":"Some descriptive statistics about the predictions."},{"metadata":{"id":"ydOVUfjgizbJ","outputId":"fed8b426-3dc7-457d-afdd-7eced79321bf","trusted":true},"cell_type":"code","source":"output.describe()","execution_count":null,"outputs":[]},{"metadata":{"id":"TAvtrIVdVAVy"},"cell_type":"markdown","source":"Saving of the predictions."},{"metadata":{"id":"rnSMRH1FizbN","trusted":true},"cell_type":"code","source":"output.to_csv('output_gbrf_1st_strategy.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"id":"wRbKw1kttfe0"},"cell_type":"markdown","source":"This classifier got a public leaderboard score of 0.78560. That is a good performance as it corresponds roughly to the first third of the leaderboard (when I write these lines).<br>\nAs say the competition description about the public leaderboard: \"This leaderboard is calculated with approximately 25% of the test data. The final results will be based on the other 75%, so the final standings may be different.\" Thereafter, we will use the terms \"private leaderboard\" to refer to the 75% of the test data used to calculate the final results.<br>\nOn the private leaderboard, this classifer (alone) got an AUROC equal to 0.78713."},{"metadata":{"id":"tzotynpusOIF"},"cell_type":"markdown","source":"### III.2 Second strategy.<a id='II.2'></a> <br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; The second classifier of our ensemble is a logistic regression. The trick we used to encode the features day and month is not suited for this model. Indeed, a linear model cannot take into account the interactions between sin_day and cos_day and between sin_month and cos_month. It just associate a coefficient to each feature, modelizing its direct impact on the target. Thus, before training a logistic regression, we remove the features sin_day, cos_day, sin_month, and cos_month and replace them with a one-hot encoding of the original features day and month. A one-hot encoding do not break the cyclical aspect of these features. What is more, we keep the dummy variable indicating the missingness for both day and month (to remove it would be equivalent to a mode imputation).<br>\n&nbsp;&nbsp;&nbsp;&nbsp; The introduction of an indicator variable of the missingness of the binary variables bin_0, bin_1, bin_2, bin_3 and bin_4 has been considered as an hyperparameter. The creation of such indicator variables for the missingness of bin_2 and bin_4 afforded slightly better results and allowed to add another difference with the first model. <br>"},{"metadata":{"id":"MZ5cKMiNFFHI","trusted":true},"cell_type":"code","source":"# Creation of the training set used with the logistic regression.\ndf_total = pd.concat([df.drop('target', axis=1), test])\n\none_hot_cols = ['day','month']\nX_cat = df_total[one_hot_cols].astype(str).copy()\n\none_hot_enc = OneHotEncoder()\nX_cat_1hot = one_hot_enc.fit_transform(X_cat)\nX_cat_1hot = pd.DataFrame(X_cat_1hot.toarray(), index=X_cat.index, \n                          columns = pd.Series(one_hot_enc.get_feature_names()).str.replace('x0', 'day').str.replace('x1', 'month'))\nX_cat_1hot.drop(['day_' + X_cat['day'].mode()[0], 'month_' + X_cat['month'].mode()[0]], axis=1, inplace=True)\n\nX_lr_total = pd.concat([X, X_test])\nX_lr_total = pd.concat([X_lr_total.drop(['sin_day', 'cos_day', 'sin_month', 'cos_month'], axis=1), X_cat_1hot], axis=1)\n\nX_lr_total['bin_2_nan'] = df_total.bin_2.isna().astype(int)\nX_lr_total['bin_4_nan'] = df_total.bin_4.isna().astype(int)\n\nX_lr = X_lr_total.iloc[:600000]\nX_test_lr = X_lr_total.iloc[600000:1000000]","execution_count":null,"outputs":[]},{"metadata":{"id":"LVlRxzvgZtwJ"},"cell_type":"markdown","source":"Quick visualization of the new training set."},{"metadata":{"id":"mNipo_PvBUVq","outputId":"826e9df5-bafe-4a3b-87b7-afa60cbfe164","trusted":true},"cell_type":"code","source":"X_lr.shape","execution_count":null,"outputs":[]},{"metadata":{"id":"4z6N0ftSR0tk","outputId":"a0dbf751-eb96-447b-e38b-3fa26400de70","trusted":true},"cell_type":"code","source":"X_lr.head().T","execution_count":null,"outputs":[]},{"metadata":{"id":"1ngJSD4fxAQx"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; A logistic regression is a very simple model, which is more prone to underfitting than to overfitting in our case. To add interactions terms allow to consider more complex patterns. Nevertheless, there are a lot of possible interactions to consider. So let's look at the interactions exploited by a tree-based model to have an idea of the profitable ones.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; The gradient boosted random forest we trained in the precedent part just used first order interactions (as the depth of its trees was limited to 2). To consider more complex interactions, let's train gradient boosted trees with a depth up to 4. <br>\n&nbsp;&nbsp;&nbsp;&nbsp; Here, to save time, we do not use a k-fold cross-validation. Instead, we just split the training set once at random into training and validation sets (what is known as the holdout method). <br>"},{"metadata":{"outputId":"4ca068fa-d513-450d-fd5c-e500cbbd8930","id":"EVJNJgWXMrZK","trusted":true},"cell_type":"code","source":"start = datetime.now()\n\nX_train, X_valid, y_train, y_valid = train_test_split(X_lr, y, test_size=0.2, stratify=y, random_state=42)\nweight = sum(y_train.values==0)/sum(y_train.values==1)\n    \ngbrf = xgb.XGBClassifier(n_estimators=3000, random_state=0, objective='binary:logistic', scale_pos_weight=weight, \n                                learning_rate=0.10, max_depth=4, subsample=0.7, min_child_weight=500,  \n                                colsample_bytree = 0.2, reg_lambda = 3.5, reg_alpha=1.5, num_parallel_tree= 5)\n    \ngbrf.fit(X_train, y_train, eval_set= [(X_train, y_train), (X_valid, y_valid)], eval_metric=['auc'], \n                    early_stopping_rounds=100, verbose=False)\n    \ny_pred = gbrf.predict_proba(X_train)\nprint(\"Train AUROC = \", roc_auc_score(y_train, y_pred[:,1]))\n\ny_pred = gbrf.predict_proba(X_valid)\nprint(\"Validation AUROC = \", roc_auc_score(y_valid, y_pred[:,1]))\n\nruntime = datetime.now() - start\nprint('\\n\\nRuntime ', runtime)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"joblib.dump(gbrf, 'model_gbrf_2nd_strategy')","execution_count":null,"outputs":[]},{"metadata":{"id":"BT4o_GjOvH8l"},"cell_type":"markdown","source":"Now Xgbfir is used to produce a file containing informations about the importance of each interaction."},{"metadata":{"id":"xYyZCwOOqv-a","trusted":true},"cell_type":"code","source":"xgbfir.saveXgbFI(gbrf, feature_names=X_lr.columns, MaxInteractionDepth=3, OutputXlsxFile='fi_gbrf_2nd_strategy.xlsx')","execution_count":null,"outputs":[]},{"metadata":{"id":"C66TDVALwH_P","trusted":true},"cell_type":"code","source":"xls = pd.ExcelFile('fi_gbrf_2nd_strategy.xlsx')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here, we are interested by the columns \"Gain\" of the following tables (the gain has the same sense as when we interpreted the feature importances of the first classifier). The meaning of the other columns can be found in the Xgbfir's github page, available [here](https://github.com/limexp/xgbfir)."},{"metadata":{"id":"cVg9n4oYvqrL","trusted":true},"cell_type":"code","source":"fi_depth_1 = pd.read_excel(xls, 'Interaction Depth 1')\nfi_depth_2 = pd.read_excel(xls, 'Interaction Depth 2')\nfi_depth_3 = pd.read_excel(xls, 'Interaction Depth 3')","execution_count":null,"outputs":[]},{"metadata":{"id":"KEh2gWwUvzWd","outputId":"3c4496a8-cce6-40a6-d56f-e41a63cb2ef3","trusted":true},"cell_type":"code","source":"fi_depth_1[:10]","execution_count":null,"outputs":[]},{"metadata":{"id":"ePyyt64UInLY","outputId":"54639b28-0cc4-4bb3-8c06-e35a3d5e8820","trusted":true},"cell_type":"code","source":"fi_depth_2[:10]","execution_count":null,"outputs":[]},{"metadata":{"id":"AeF_FlOrIvtY","outputId":"e47f078f-fec4-48a2-9a50-025817bbc649","trusted":true},"cell_type":"code","source":"fi_depth_3[:10]","execution_count":null,"outputs":[]},{"metadata":{"id":"L1T_zF0HxdCN"},"cell_type":"markdown","source":"This information allows to save time by considering the inclusion of only a few interaction terms. <br>\n<br>\nMost of the interactions bring no improvement to the logistic regression. In the case where the interactions between the features do not impact the target, the trees just split the data set along the most gainful features regardless of the previous splits. The precedent visualizations showed that the interactions which bring the greatest gains are those which relate to the most influential predictors. So these results are consistent with the assumption of a low impact of the interactions. <br>\nTo include the product of ord_0, ord_3 and ord_5 still improves the performance.<br>\n<br>\nFurthermore, some predictors have been succesfully removed based on the coefficients of the linear combination used by the logistic regresssion (after having iteratively trained the classifier with different features). This concerns bin_3 and the dummy variables corresponding to the modality Square of nom_1, the modality Cat of nom_2, the modality Oboe of nom_4, and the modality 5 of day. <br>"},{"metadata":{"id":"xcmqoSaHxZ7g","outputId":"63072c6a-966a-4b2e-e3b3-3db78e1ca5fa","trusted":true},"cell_type":"code","source":"X_lr_total['int_1'] = X_lr_total['ord_0'] * X_lr_total['ord_3'] * X_lr_total['ord_5']\n\nX_lr_total.drop(['bin_3', 'x1_Square', 'x2_Cat', 'x4_Oboe', 'day_5.0'], axis=1, inplace=True)\n\nX_lr = X_lr_total.iloc[:600000]\nX_test_lr = X_lr_total.iloc[600000:1000000]","execution_count":null,"outputs":[]},{"metadata":{"id":"wQRm6SpQZtwS"},"cell_type":"markdown","source":"The features are standardized before training in order to make coefficients of the classifier more representative of the actual impact of the features. We use again a 5-fold cross-validation strategy."},{"metadata":{"id":"iZkNvPBgZtwT","outputId":"43d79ef0-d72a-428a-bf9e-8a5c9061d301","trusted":true},"cell_type":"code","source":"start = datetime.now()\n\ntrain_scores = []\nval_scores = []\ny_validations = []\ny_predictions = []\n\nfor index_train, index_val in StratifiedKFold(n_splits=5, random_state=42, shuffle=True).split(X_lr, y):\n    X_train = X_lr.iloc[index_train]\n    y_train = y.iloc[index_train]\n    X_valid = X_lr.iloc[index_val]\n    y_valid = y.iloc[index_val]\n    \n    # Logistic regression.\n    lr = LogisticRegression(random_state=0, max_iter=10000, fit_intercept=True, class_weight='balanced', \n                            penalty='l2', C=0.4, verbose=0)\n\n    scaler = StandardScaler()\n    X_train_scaled = scaler.fit_transform(X_train)\n    X_valid_scaled = scaler.transform(X_valid)\n\n    lr.fit(X_train_scaled, y_train)\n\n    # We stock the results.\n    y_pred_train = lr.predict_proba(X_train_scaled)[:,1]\n    train_score = roc_auc_score(y_train, y_pred_train)\n    train_scores.append(train_score)\n\n    y_pred_valid = lr.predict_proba(X_valid_scaled)[:,1]\n    val_score = roc_auc_score(y_valid, y_pred_valid)\n    val_scores.append(val_score)\n    \n    y_validations.append(y_valid)\n    y_predictions.append(y_pred_valid)\n    \ndisplay_scores(train_scores, val_scores, y_validations, y_predictions)\n\nruntime = datetime.now() - start\nprint('\\n\\nRuntime ', runtime)","execution_count":null,"outputs":[]},{"metadata":{"id":"yS-QT12SRHSc"},"cell_type":"markdown","source":"We got an average AUROC of roughly 0.7877 over the cross-validation.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; As the logistic regression just compute the estimated probabilities by applying the logistic function to a linear combination of the standardized features, then the coefficients of this linear combination are generally representative of the importances of the features in the classification (the presence of multicolinearities can yet cause misleading coefficients). Larger is the absolute value of a coefficient, greater is its influence on this classification. The sign of the coefficients indicate how an increase of the features impact the estimated probabilities (positively or negatively).<br>"},{"metadata":{"id":"hRc-NBEr5VZD","trusted":true},"cell_type":"code","source":"lr_coefficients = pd.DataFrame({'Predictors': X_lr.columns.to_list() + ['Intercept'], \n                                'Coefficients': list(lr.coef_.reshape(-1)) + list(lr.intercept_),\n                                'Direction': np.sign(list(lr.coef_.reshape(-1)) + list(lr.intercept_)),\n                                'Importances': np.abs(list(lr.coef_.reshape(-1)) + list(lr.intercept_))})\n\nlr_coefficients.sort_values('Importances',ascending=False, inplace=True)\nlr_coefficients.reset_index(inplace=True, drop=True)\nlr_coefficients.to_csv('coefficients_lr_2nd_strategy')","execution_count":null,"outputs":[]},{"metadata":{"id":"LWyneT0M5VZH","outputId":"0589afd4-8656-425f-9723-e8f0f481d109","trusted":true},"cell_type":"code","source":"lr_coefficients[:50]","execution_count":null,"outputs":[]},{"metadata":{"id":"CkjK9j1N_nMA","outputId":"0c560e3e-8d25-4e30-d8aa-a0c904a35b92","trusted":true},"cell_type":"code","source":"lr_coefficients[50:]","execution_count":null,"outputs":[]},{"metadata":{"id":"2tqHqxgAEZmB"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; The features' importances are similar to the ones rendered by the gradient boosted random forest we trained before, despite the greater influence of the target encoded variables (especially nom_8 and nom_9) on the logistic regression. As the trees were allowed to use first order interactions, that suggests their impact was small (otherwise, the features' importances should be more different because the logistic regression considered just one interaction).<br>\n&nbsp;&nbsp;&nbsp;&nbsp; As a remark, the precedent statement could not be deduced from the independence of the features. Indeed, I previously reported that it was unlikely that I achieve to predict efficiently a feature from the others (given the results of a basic naive bayesian classifier). From that, we can deduce the features were not much linked. However, the target variable could have been affected by the interactions between unlinked features.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; As the Kaggle's description of the competition evoked the presence of interactions, it is surprising that they are so weak. Thus, I am not completely satisfied that I did not find stronger interactions. <br>\n\nLet's now train the logistic regression on the entire training set and submit its predictions on the test set."},{"metadata":{"id":"DEO-1bH7Pxab","outputId":"c8507c76-cb08-4449-f257-70a75b6ef6c8","trusted":true},"cell_type":"code","source":"# Full training.\nlr = LogisticRegression(random_state=0, max_iter=10000, fit_intercept=True, class_weight='balanced', penalty='l2', \n                        C=0.4, verbose=0)\n\nscaler = StandardScaler()\nX_lr_scaled = scaler.fit_transform(X_lr)\nX_test_lr_scaled = scaler.transform(X_test_lr)\n\nlr.fit(X_lr_scaled, y)","execution_count":null,"outputs":[]},{"metadata":{"id":"_B5Xj4lO5VYu","outputId":"02509d01-7c75-4999-aae8-2d38e3731b36","trusted":true},"cell_type":"code","source":"joblib.dump(lr,'model_lr_2nd_strategy')","execution_count":null,"outputs":[]},{"metadata":{"id":"Zt7h3DHf5VZK","outputId":"36c1e349-c938-48d0-fce8-655e608d780e","trusted":true},"cell_type":"code","source":"#### Submissions.\ny_test = lr.predict_proba(X_test_lr_scaled)[:,1]\noutput = pd.DataFrame({'id':sample_submission.index, 'target':y_test})\noutput","execution_count":null,"outputs":[]},{"metadata":{"id":"i0cNqHJS5VZh","outputId":"852bc551-91ce-40de-d487-2f23d6757b45","trusted":true},"cell_type":"code","source":"output.describe()","execution_count":null,"outputs":[]},{"metadata":{"id":"09GtmWPi5VZk","trusted":true},"cell_type":"code","source":"output.to_csv('output_lr_2nd_strategy.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This logistic regression reached scores equal to 0.78552 on the public leaderboard and to 0.78702 on the private leaderboard. Once again, the difference with the gradient boosted random forest is small."},{"metadata":{"id":"weLD-zcXtfge"},"cell_type":"markdown","source":"### III.3 Third strategy.<a id='III.3'></a> <br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; In this section, we will train a naive bayesian classifier. As this model does not consider interactions, the choice of the features to include must be done based on the relationships between each one and the target. Whereof, a simple way to get an idea of the useless features is to use chi-squared tests of independence between each feature and the target.<br>\n&nbsp;&nbsp;&nbsp;&nbsp; In contrast, I did not conducted tests with the encoded features used for the logistic regression as the coefficients of the linear combination already gave an idea of the features' impact on this classifier, allowing to select the useful features iteratively."},{"metadata":{"outputId":"97158a5b-599e-4035-d8f5-5ae50d03f275","id":"sRfPc1yPS3hH","trusted":true},"cell_type":"code","source":"start = datetime.now()\n\nchi2_target = pd.DataFrame(columns=['Cochran_criterion', 'Chi2', 'p_value'])\nfor col in df.columns.drop('target'):\n  feature = df[col]\n  chi2_target.loc[col + '_target'] = chi2_test(feature, df.target)\n\nchi2_target['rejected'] = chi2_target['Cochran_criterion'] & (chi2_target['p_value']<0.05)\nchi2_target.sort_values('p_value', ascending=True, inplace=True)\nchi2_target.to_csv('chi2_target_tests.csv', index=False)\nprint(chi2_target, '\\n')\n\nruntime = datetime.now() - start\nprint('Runtime ', runtime)","execution_count":null,"outputs":[]},{"metadata":{"id":"S3y7r8hfZtxH"},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; Only the test between bin_3 and the target was not rejected. Unsurprisingly, to remove bin_3 improved the model.<br>\nI considered as an hyper-parameter the inclusion of other features (particularly the ones whose the link with the target as measures by the chi-squared statistic is weak). To drop nom_6 happened to be profitable. That was not the case of the other variables I tried to drop. <br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Here, for each feature, the missing values are considered as a particular category. As approximately 18000 values are missing for each feature, it is enough to estimate reasonably their probabilities. Moreover, as stated before, to handle the missing values differently increases the propensity of the different classifiers to learn different aspects of the data and so to form an efficient ensemble.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; With a naive bayesian classifier, the little frequent categories are a source of overfitting. Indeed, the estimated probabilities of belonging to these categories are calculated with few examples, so they are potentially inaccurate. Consequently, we delete categories with less than 100 observations and consider them as missing values.<br>"},{"metadata":{"trusted":true},"cell_type":"code","source":"for col in ['nom_'+str(i) for i in range(5,10)]:\n    rare_categories = df[col].value_counts().iloc[np.where(df[col].value_counts()<100)].index\n    print(\"The predictor \" + col + \" has \", rare_categories.size, \" rare categories among its \", df[col].unique().size,\n         \" categories.\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"&nbsp;&nbsp;&nbsp;&nbsp; The alpha hyperparameter (tackled previously) has been set to 3.4 (thanks to a gridsearch) to reduce overfitting. <br>\nFor each feature, all the categories present in the test set and absent from the training set are considered as a single new category. As this category is new, the classifier would assign it a probability of 0 without the smoothing parameter. So here is another interest of the alpha hyperparameter. However, in our case, that concerns only one value of nom_6. <br>"},{"metadata":{"id":"IKeSUHOHOVcc","outputId":"41a73fd2-cd19-448d-863b-29ba3abe9028","trusted":true},"cell_type":"code","source":"for col in df.columns.drop('target'):\n  unknowns = set(test[col].astype(str).unique()) - set(df[col].astype(str).unique())\n  print('The predictor ', col, ' contains ', len(unknowns), ' value present only in the test set.')","execution_count":null,"outputs":[]},{"metadata":{"id":"toKSYLU33WxB"},"cell_type":"markdown","source":"As nom_6 has been dropped, it is not an issue.<br>\n<br>\nFor the classifier to consider the interactions between two features, we can drop the two features and replace them with a new variable taking a different values for each pair of modalities a sample can take for these two features. It has been tried for some pairs of features, but it was vain (the classifier is already overfitting).<br>\n<br>\nBelow is the cross-validation."},{"metadata":{"id":"Vyuwye_TLVTU","outputId":"9177e8de-c256-440f-c5aa-30fa23f53ec8","trusted":true},"cell_type":"code","source":"start = datetime.now()\n\ntrain_scores = []\nval_scores = []\ny_validations = []\ny_predictions = []\n\nX_nb = df.drop(['bin_3', 'nom_6', 'target'], axis=1).astype(str)\n\nfor index_train, index_val in StratifiedKFold(n_splits=5, random_state=42, shuffle=True).split(X_nb, y):\n        \n    X_train = X_nb.iloc[index_train].copy()\n    y_train = y.iloc[index_train].copy()\n    X_valid = X_nb.iloc[index_val].copy()\n    y_valid = y.iloc[index_val].copy()\n\n    for col in ['nom_5', 'nom_7', 'nom_8', 'nom_9']:\n      rare_categories = X_train[col].value_counts().iloc[np.where(X_train[col].value_counts()<100)].index\n      X_train.loc[:,col] = X_train.loc[:,col].replace(rare_categories, 'nan')\n      X_valid.loc[:,col] = X_valid.loc[:,col].replace(rare_categories, 'nan')\n\n    # Ordinal encoding.\n    oe = OrdinalEncoder()\n    X_train_oe = oe.fit_transform(X_train)\n    X_val_oe = oe.transform(X_valid)\n\n    X_train_oe += 2 # we add 2 to avoid negative values resulting from the ordinal encoding\n    X_val_oe += 2\n\n    # Training.\n    nb = CategoricalNB(alpha = 3.4)\n    nb.fit(X_train_oe, y_train)\n\n    # We stock the results.\n    y_pred = nb.predict_proba(X_train_oe)\n    train_score = roc_auc_score(y_train, y_pred[:,1])\n    train_scores.append(train_score)\n\n    y_pred = nb.predict_proba(X_val_oe)\n    val_score = roc_auc_score(y_valid, y_pred[:,1])\n    val_scores.append(val_score)\n\n    y_validations.append(y_valid)\n    y_predictions.append(y_pred[:,1])\n\ndisplay_scores(train_scores, val_scores, y_validations, y_predictions)\n\nruntime = datetime.now() - start\nprint('\\n\\nRuntime ', runtime)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We got a validation score equal to 0.7864.<br>\n<br>\nLet's preprocess the entire training set and the test set."},{"metadata":{"id":"fLfjCng87pqs","trusted":true},"cell_type":"code","source":"# Full preprocessing\noe = OrdinalEncoder()\n\nX_nb = df.drop(['bin_3', 'nom_6', 'target'], axis=1).astype(str)\nX_test_nb = test.drop(['bin_3', 'nom_6'], axis=1).astype(str)\n\nfor col in ['nom_5', 'nom_7', 'nom_8', 'nom_9']:\n  rare_categories = X_nb[col].value_counts().iloc[np.where(X_nb[col].value_counts()<100)].index\n  X_nb[col] = X_nb[col].replace(rare_categories, 'nan')\n  X_test_nb[col] = X_test_nb[col].replace(rare_categories, 'nan')\n\nX_oe = oe.fit_transform(X_nb)\nX_test_oe = oe.transform(X_test_nb)\n  \nX_oe += 2 # we add 2 to avoid negative values resulting from the ordinal encoding\nX_test_oe += 2","execution_count":null,"outputs":[]},{"metadata":{"id":"GbEGhY7O7IIt"},"cell_type":"markdown","source":"To finish, here are the full training and the predictions."},{"metadata":{"id":"BW4dvaOKSFRw","outputId":"49d0793d-0dfe-4fb0-c4ba-4dd72f232c26","trusted":true},"cell_type":"code","source":"start = datetime.now()\n\n# Full training\nnb = CategoricalNB(alpha = 3.4) \nnb.fit(X_oe, df.target)\n\ny_pred = nb.predict_proba(X_oe)\ntrain_score = roc_auc_score(df.target, y_pred[:,1])\nprint('AUROC on the entire training set: ', train_score)\n\nruntime = datetime.now() - start\nprint('\\n\\nRuntime ', runtime)","execution_count":null,"outputs":[]},{"metadata":{"id":"01k-Mz-0TLEG","outputId":"25650c6b-eb6d-4a1e-9a3b-17abb35bc568","trusted":true},"cell_type":"code","source":"joblib.dump(nb, 'model_nb_3rd_strategy')","execution_count":null,"outputs":[]},{"metadata":{"id":"3_y9pxGYTS5R","outputId":"460f0910-1a2a-46d1-cb9a-7e83f2d5d8de","trusted":true},"cell_type":"code","source":"y_test = nb.predict_proba(X_test_oe)[:,1]\noutput = pd.DataFrame({'id':sample_submission.index, 'target':y_test})\noutput","execution_count":null,"outputs":[]},{"metadata":{"id":"iGbFMjKITiWa","outputId":"7ed5ba1b-bf05-439f-f80e-50a8c9a9b8e9","trusted":true},"cell_type":"code","source":"output.describe()","execution_count":null,"outputs":[]},{"metadata":{"id":"vA8iK5eRTIv6","trusted":true},"cell_type":"code","source":"output.to_csv('output_nb_3rd_strategy.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"id":"bImfhZWfApQZ"},"cell_type":"markdown","source":"This classifier got a public leaderboard score of 0.78477. It is smaller than with the precedent models. Nevertheless, as this classifier is different from what we have done so far (except the target encoding wich is a bayesian encoding method), it is expected to improve our ensembling.<br>\nIts private leaderboard score is equal to 0.78621.\n\n### III.4 Ensembling.<a id='III.4'></a> <br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; Now we compute a weighted average of the estimated probabilities of our different classifiers. The weights of each classifier is just its public leaderboard score (as the private leaderboard score was not yet available when I wrote this code). Thus, the more efficient a model, the more weight it is given."},{"metadata":{"id":"gcVpMzd0tfhk","trusted":true},"cell_type":"code","source":"output_gbrf = pd.read_csv('output_gbrf_1st_strategy.csv')\noutput_lr = pd.read_csv('output_lr_2nd_strategy.csv')\noutput_nb = pd.read_csv('output_nb_3rd_strategy.csv')","execution_count":null,"outputs":[]},{"metadata":{"id":"l4IXWe-W_lfh","trusted":true},"cell_type":"code","source":"mean = output_nb\nmean.target = (output_gbrf.target * 0.78560 + output_lr.target * 0.78552 + output_nb.target * 0.78477)/(0.78560 + 0.78552 + 0.78477)","execution_count":null,"outputs":[]},{"metadata":{"id":"clu0poDC_w70","trusted":true},"cell_type":"code","source":"mean.to_csv('ouput_ensemble.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This ensembling got a public leaderboard score of 0.78605 and a private leaderboard score of 0.78757"},{"metadata":{"id":"Yz4lzqMOtfhp"},"cell_type":"markdown","source":"## IV. Conclusion and potential extensions.<a id='IV'></a> <br>\n\n&nbsp;&nbsp;&nbsp;&nbsp; This work allowed me to reach the top 21% (with the rank 240 among 1161 participants) of the final standing (computed with the AUROC on the private leaderboard). Nevertheless, I see numerous way of improvement (my spare time and computing power were limited, so I did what I could). <br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; As the Kaggle's description of the competition mentionned the presence of interactions, I am not completely satisfied that I did not find stronger interactions (as evoked before). Indeed, to increase the depth of the tree-based methods that I used worsen them (because of overfitting), the inclusion of most of the possible interactions worsen the logistic regression, and the naive bayesian classifier does not consider interactions at all. Moreover, the features' importances measured by the coefficients of the logistic regression and by the gain of the gradient boosted random forest are similar, suggesting a small impact of these first order interactions. With more spare time, I would try gain to unearth them. Maybe tree-based methods with deeper trees could succeed to discover them. But, to override overfitting, extra-trees could be more suited that random forest and gradient boosted trees (I had begun to use extra-trees without success).<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; Besides, the scikit-learn's implementation of the naive bayesian classifier I used do not enable to use a different alpha hyperparameter for each feature. However, high cardinality features would certainly benefit from a higher value of alpha than the low cardinality features. Indeed, when a feature have numerous categories, each one cover few samples. But for the estimated probability of belonging to a categories to be efficient, we need to have enough samples from that category. Contrariwise, the categories of a low cardinality feature cover more samples. So the resulting estimated probabilites of belonging to each category are better and less regularization should provide better result.<br>\nConsequently, to create my own naive bayesian classifier (or find one well suited to my needs) could have provided a better performance.<br>\n<br>\n&nbsp;&nbsp;&nbsp;&nbsp; I could have tried more models with more hyperparameters combination when I tried to predict either the missingness or the values taken by the features. Maybe I missed something.<br>\nThe preprocessings (encoding and missing value handling) of the predictors were almost the same for the first two strategies. To include to the ensembling models trained on differently preprocessed data is another opportunity of improvement.<br>\n<br>\n## V. Bibliographic references. <a id='V'></a><br>\n<br>\n\nThis bibliography aims to point the reader to appropriate sources covering the concepts used in this notebook. They are provided with links to simplify the potential searches of the reader.<br>\n<br>\n\n[All of statistics, by Larry Wasserman (Springer).](https://www.amazon.com/All-Statistics-Statistical-Inference-Springer/dp/0387402721/ref=sr_1_1?ie=UTF8&qid=1249141007&sr=8-1) This book is an introductory book on statistics which cover a broad range of topics. In particular, it covers the chi-squared test of independence I used (with examples) and many others.<br>\n<br>\n[Hands-on Machine Learning with Scikit-Learn and Tensorflow, by Aurélien Géron (O'Reilly).](Link : https://www.oreilly.com/library/view/hands-on-machine-learning/9781492032632/) <br>\n<br>\n[Deep Learning with Python, by François Chollet (Manning).](https://www.manning.com/books/deep-learning-with-python)<br>\n<br>\n[Data Analysis Using Regression and Multilevel/Hierarchical Models by Andrew Gelman and Jennifer Hill (Cambridge University Press).](http://www.stat.columbia.edu/~gelman/arm/) <br>\nThe chapter 5 of this book covers the most used methods to handle the missing values.<br>\n<br>\nSimply put: here are the wikipedia articles on [Missing data](https://en.wikipedia.org/wiki/Missing_data), [Chi-squared tests](https://en.wikipedia.org/wiki/Pearson%27s_chi-squared_test), [Data leakage](https://en.wikipedia.org/wiki/Leakage_%28machine_learning%29), and [Cross-validation](https://en.wikipedia.org/wiki/Cross-validation_%28statistics%29). <br>\n<br>\n[Here](https://towardsdatascience.com/smarter-ways-to-encode-categorical-data-for-machine-learning-part-1-of-3-6dca2f71b159) is a clear article dealing with the mainly used feature encoding methods. <br>\n<br>\n[This link](https://ianlondon.github.io/blog/encoding-cyclical-features-24hour-time/) points to an article dealing with the projection of cyclical features in two dimensions. <br>\n<br>\n[Here](https://www.kaggle.com/caesarlupum/2020-20-lines-target-encoding) is a simple and straight notebook using a 5-fold splitting strategy for target encoding.  <br>\n<br>\n[This](https://github.com/limexp/xgbfir) is a link to the Xgbfi's github page.<br>\n<br>\nAs I was brief when dealing whith the gradient boosted random forest, here are two links from the XGBoost's documentation: [An introduction to boosted trees](https://xgboost.readthedocs.io/en/latest/tutorials/model.html) and [A guide to the available parameters](https://xgboost.readthedocs.io/en/latest/parameter.html) <br>\n<br>\nThe scikit-learn library contains great user guides available [here](https://scikit-learn.org/stable/user_guide.html). In particular, [this one](https://scikit-learn.org/stable/modules/naive_bayes.html) deals with the naive bayesian methods available in scikit-learn.<br>\n<br>\n<br>\n<br>\nThanks for reading. I hope you have taken as pleasure by reading this notebook as I have taken by writing it.<br>\nI will conclude by quoting the words of a brilliant data scientist in the introduction of this notebook (I love to quote myself): \"If this notebook pleased you in any way, or learned you anything, you can upvote it (don't be timid, it is free). As I plan to write other notebooks in the future, you can follow me to be informed when that will happen\"."}],"metadata":{"accelerator":"GPU","colab":{"name":"Categorical_Feature_Encoding_Challenge_II_notebook_6.ipynb","provenance":[],"collapsed_sections":[],"machine_shape":"hm"},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.3"}},"nbformat":4,"nbformat_minor":4}