<jupyter_start><jupyter_text>IntroductionIn this section, we will proceed to explain in general terms the process followed for the project. The process is split into several stages to obtain as many insights as possible from the data:1.  **Problem definition:** An informal definition of the problem is detailed, along with the initial assumptions.2.  **Exploratory Data Analysis (EDA):** The data is analyzed using statistical methods and unsupervised modeling to gain insights that may help the posterior data preparation and modeling stages. This step is also useful to verify our initial assumptions.3.  **Data Preparation:** This step involves data transformation, imputation of missing values, etc.4.  **Modeling or Spot-Checking:** We evaluate different supervised models and select the best candidates among them.5.  **Final model:** We choose the final model and retrain it with all data available (train and test set).6.  **Results:** We provide a summary of the results obtained so far.  Problem DefinitionThe problem we want to tackle is a classification problem based on the Census Income dataset. This dataset shows different variables regarding individuals, such as _education_, _sex_, and other attributes. Our goal is to create a model to predict whether an individual earns more than 50k dollars per year or not. This dataset was created using Census data from American individuals in 1996. It has been an object of study and benchmark on multiple occasions.In this project, we will try to approach the problem using modern tools and domain knowledge.  Exploratory Data Analysis<jupyter_code>from project import *

pd.set_option('display.latex.repr', True)
configure_plotting()

train, test = load_dataset()
train.head()<jupyter_output><empty_output><jupyter_text>The table above represents a small sample of the dataset. As we can see, there are several variables which some of them are not self-explanatory. Because of that, we need to dig deeper into the dataset description to find out the meaning of these variables:-   **age:** Individual's age-   **workclass:** Type of job of the individual.-   Private, Self-emp-not-inc, Self-emp-inc, Federal-gov, Local-gov, State-gov, Without-pay, Never-worked.-   **fnlwgt**: This continuous variable represents the final weight, which is the number of units in the target population that the sample represents.-   **education**: Level of education of the individual.-   Bachelors, Some-college, 11th, HS-grad, Prof-school, Assoc-acdm, Assoc-voc, 9th, 7th-8th, 12th, Masters, 1st-4th, 10th, Doctorate, 5th-6th, Preschool.-   **education-num**: This numeric variable has the same meaning as the previous one, but with a numerical and ordered representation. That means that we can discard the variable _education_ in favor of this one.-   **marital-status**: Martial status of the individual:-   Married-civ-spouse, Divorced, Never-married, Separated, Widowed, Married-spouse-absent, Married-AF-spouse.-   **occupation**: Work sector of the individual:-   Tech-support, Craft-repair, Other-service, Sales, Exec-managerial, Prof-specialty, Handlers-cleaners, Machine-op-inspct, Adm-clerical, Farming-fishing, Transport-moving, Priv-house-serv, Protective-serv, Armed-Forces.-   **relationship**: Represents the relationship of the individual with its family.-   Wife, Own-child, Husband, Not-in-family, Other-relative, Unmarried.-   **race**: Individual's race:-   White, Asian-Pac-Islander, Amer-Indian-Eskimo, Other, Black.-   **sex**: Individual' sex:-   Male, Female-   **capital-gain**: Yearly capital gains.-   **capital-loss**: Yearly capital losses.-   **hours-per-week**: Work hours per week.-   **country**: Country of origin.-   **target**: The target variable. It is a categorical variable that indicates if an individual earns more than 50k dollars or not.As we can see, most of the variables that will give us some insights into the problem. However, it may seem that some important variables are not included but may be useful, such as the number of children. Presumably, we can think that people that have children earn higher salaries to cover the costs of raising a child. Another useful variable would be the GDP per capita of the state where the individual works, as we may find a disparity of salaries between states.<jupyter_code>_, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
ax1.set_title('Train set')
ax2.set_title('Test set')
sns.countplot(train['target'], ax=ax1)
sns.countplot(test['target'], ax=ax2)
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>As we face a classification problem, it is always good to split the data in a stratified way, so the ratio of samples for each class remains the same for both training and test set. The plot above shows the ratio of positive and negative samples for each dataset. We can also notice how the ratio of positive samples is much higher than the negative one.  Numerical Variables Analysis> The rest of the EDA will be carried on using the training set only to avoid introducing any bias in our modeling that may affect the generalization error.For this dataset, only 6 out of 14 variables are numeric. The following table describes those variables briefly:<jupyter_code>continuous_columns = ['age', 'fnlwgt', 'education-num','capital-gain',
                   'capital-loss', 'hours-per-week']
get_summary(train[continuous_columns])<jupyter_output><empty_output><jupyter_text>In terms of nullity, no continuous variable contains null values, so we can skip imputation for those columns.The **Age** variable ranges between 17-90. It is highly likely that both extremes are outliers, or the distribution has long tails, as the first and third quartiles are 28 and 48, respectively.The **education-num** variable is, in fact, an ordinal variable codified as numbers, so we will analyze the **education** variable as the labels are self-explanatory. As far as **capital-gain** and **capital-loss** are concerned, the majority of values are zero. It does not mean that they contain null values, but that most of the people do not gain or lose capital. Presumably, we may say that only self-employed and business owners can suffer from capital changes.Finally, the **hours-per-week** variable ranges from 1 to 99. Both extremes could be outliers, as both values are far from Q1 and Q2. In fact, 99 hours per week is more than 14 hours per day, from Monday to Sunday. An amount that is far from being legal for employees and too high for self-employed people as well.The best way to interpret numerical variables is by plotting histograms:<jupyter_code>fig, axes = plt.subplots(3, 2, figsize=(13, 10))
for ax, col in zip(axes.flatten(), continuous_columns):
    sns.distplot(train[col], kde=False, bins=20, ax=ax)
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>As we can see, no variable follows a priori, a Gaussian distribution. Besides, as we commented before, most individuals have no losses or capital gains. Moreover, it is clear how more than 60 hours of work per week can consider an outlier. On the other hand, we see how age, although it is not a pure Gaussian distribution, is quite close.We will now see how the distributions of each variable differ depending on the labels.<jupyter_code>fig, axes = plt.subplots(3, 2, figsize=(12, 12))
for ax, col in zip(axes.flatten(), continuous_columns):
    sns.boxplot(x='target', y=col, data=train, ax=ax)
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>As one could sense from the plots, the **capital-loss** and **capital-gain** variables do not provide much information. Nor do they differ depending on income. Those who earn more than 50k still have a capital loss and capital gain of 0 in most cases. For this reason, we will eliminate both variables. Also, **fnlwgt** practically does not differ in relation to the labels, so we will eliminate it as well, as it does not provide information to the problem.Although these variables do not provide much information, some other variables help us understand the problem. For example, we can see how people who earn more are usually older. That can be mainly caused by older people having more experience, and therefore they tend to have higher salaries. Usually, juniors start with low salaries, and as the years go by and they become seniors, the salaries grow proportionally. We also see how there are many points of older people with salaries above and below 50k. Since these people are probably not working, the data has been collected based on the last salary received, the retirement pension, or their pension plan. In any case, we do not consider them outliers.We can also see how high levels of education and weekly hourly rates often result in higher salaries. Let us look at the relationship between education level and age (for visualization purposes, we have taken a random sample of 10,000 individuals, rather than the entire training set).<jupyter_code>sns.relplot(y='education-num', x='age', hue='target', data=train.sample(10000))
plt.show()<jupyter_output><empty_output><jupyter_text>As we see in the scatter plot, middle-aged adults with high educationare most likely to have salaries above 50k, while young people or adults who are not highly educated tend to have salaries below that threshold.As a final part of our continuous variables analysis, we will see the different correlations between the variables to discover whether there is colinearity or not. Because, as we know, if two variables have an almost perfect correlation, we can eliminate one of them and not lose much information.<jupyter_code>plot_correlation_matrix(train[continuous_columns], figsize=(12, 6), annot_size=15)
plt.show()<jupyter_output><empty_output><jupyter_text>Figure X shows no strong correlation between the data, so we do not have to remove any features due to multicollinearity. We have also calculated Spearman's correlation coefficients. Since no variable is Gaussian, it makes much more sense to perform a non-parametric test based on ranking than a parametric one, which assumes normally distributed variables.After carrying out the analysis of continuous variables, we will keep the variables age, education-num, and hours-per-week exclusively due to the information provided by each one of them.  Categorical Variables Analysis Most of the variables for our data set are categorical. Before converting them to numerical variables,we will analyze them deeply to increase our understanding of the problem.As we did with the continuous variables, we will see if there are null values. The following table shows the percentage of null values for each variable:<jupyter_code>categorical_columns = ['education', 'marital-status', 'education','occupation',
                       'relationship', 'race', 'sex', 'country']
get_missing_values(train[categorical_columns])<jupyter_output><empty_output><jupyter_text>As we can see, in this case, there are null values, specifically in the **occupation** and **country** variables. That will make us have to impute them or eliminate these entries directly later.To see how each variable can influence the predictions, we will analyze the number of examples of each one in relation to the labels.<jupyter_code>plot_count_by_income('sex', train)
plt.show()<jupyter_output><empty_output><jupyter_text>The graph above shows that the proportion of women and men in each class differs considerably. Overall, there are more men than women in the data set, but the difference varies between classes. In this case, we see that the highest salaries are usually those of men, while women's proportion is more noticeable for low salaries. For this reason, we may consider this feature relevant to predictive modeling.<jupyter_code>plot_count_by_income('marital-status', train)
plt.show()<jupyter_output><empty_output><jupyter_text>For **marital status**, the distribution also varies between classes. For the *<50k* class, most individuals are single. We should also take into account that correlation does not mean causality. Because we know that youngsters tend to have lower salaries because of inexperience, and those same people are not usually married. However, it may also be possible that people with low income decide not to get married to save money. In either case, we will use this variable because it will most certainly help predict classes.<jupyter_code>g = plot_count_by_income('education', train)
plt.show()<jupyter_output><empty_output><jupyter_text>The figure above shows that the proportion of people with a college education is higher for salaries above the threshold. In contrast, the proportion of those who have only studied up to high school is more significant for low salaries. If we remember the previous analysis, the education level and age give us a good indicator of the target variable. So we can consider this feature relevant, although we will use education-num to avoid preprocessing.<jupyter_code>g = plot_count_by_income('occupation', train)
plt.show()
g = plot_count_by_income('relationship', train)
plt.show()<jupyter_output><empty_output><jupyter_text>As we see in these two graphs above, both **occupation** and **relationship** variables will help in the prediction. In both cases, the distribution changes between classes. In the first case, we see how few are the cleaners, livestock/fishers, and clergy people, who have a high salary. Most are below the $50K threshold.The second variable, we could say that is related to **marital-status**. Since the predominant classes in both variables are single and married. In the case of **relationship,** it corresponds to no family and husband/wife, respectively. Thus, we will transform this variable by another binary variable called **own-child**, which is the only value not represented in marital-status. We could even perform a chi-square test at 95% confidence on both variables' contingency table to see if a relationship exists.<jupyter_code>from scipy.stats import chi2_contingency
import statsmodels.api as sm

contingency_table = sm.stats.Table.from_data(
    train[["relationship", "marital-status"]])
chi_test_result = contingency_table.test_nominal_association()
print('CHI^2 STATISTICAL TEST RESULTS:')
print(chi_test_result)<jupyter_output>CHI^2 STATISTICAL TEST RESULTS:
df          30
pvalue      0.0
statistic   38740.875018132836<jupyter_text>As we see, the evidence suggests that we reject the null hypothesis that the data are independent of each other<jupyter_code>g = plot_count_by_income('race', train)
plt.show()
_, (ax1, ax2) = plt.subplots(1, 2)
ax1.set_title('Income <= 50K')
ax2.set_title('Income > 50K')
plot_pie(train[train['target'] == '<=50K'], 'country', ax1)
plot_pie(train[train['target'] == '>50K'], 'country', ax2)
plt.show()<jupyter_output><empty_output><jupyter_text>On the other hand, we will discard the **country** and **race** variables.Mainly because they do not provide much information. In the case ofrace, the proportions are similar between classes. Moreover, thecountry of origin variable is dominated by the United Statesfor both classes, making the variable irrelevant.<jupyter_code>g = plot_count_by_income('workclass', train)
plt.show()<jupyter_output><empty_output><jupyter_text>As we see, most individuals work for the private sector. Few (or none) have never worked or are not paid for their work. Moreover, the proportions are similar between both classes, so this feature does not provide relevant information that makes us discriminate between classes.With this, we conclude the exploratory analysis and remember that we are going to select only the following variables: **age, education-num, hours-per-week, sex, marital status, and occupation**. Additionally, we will transform the **relationship** variable into a binary (own-child) one.It is worth pointing out that we carried out an in-depth analysis of each variable because the data's dimensionality has allowed us to do so. That has made it possible for us to understand each variable's meaning and relationship with the target and reduce the number of relevant characteristics. There are situations where the sample size or the number of characteristics make the data untreatable for this type of analysis. In these cases, we would have to use automatic feature selection, reduction of dimensionality, and other methods.  Preprocessing In this section, we will apply some transformations to the data. The first transformation will consist of filtering the dataset to keep only a subset of the features, as we stated before. We will also substitute the relationship feature with a so-called own-child feature, which will be a binary feature that represents whether an individual has children or not. We will also encode the target variable to a binary variable (0 and 1).Subsequently, and depending on the type of data of each feature, different transformations will be applied:+ If the data is numerical, it will be scaled so that none of the numerical variables has more weight than others. This scaling will also consider the outliers so that it will scale according to the interquartile range (first and third quartile). As we saw in the EDA, the hours-per-week variable had values that reached 99, a clear outlier, and we must be careful.    + If the data is of a categorical type, a transformation pipeline of the following style will be applied:    1.  The missing data is imputed with the most frequent value. There are only two variables with null values: **occu pation** and **country**. As the latter has been discarded, this method only applies to the **occupation** variable.    2.  One-Hot encoding is applied to nominal features. In our case, all the categorical variables are nominal except for education. But as we know, this variable is already numerically coded in education-num. ![](assets/preprocessing.png)<jupyter_code>X_train, y_train = split_input_labels(train)
y_train = encode_labels(y_train)<jupyter_output><empty_output><jupyter_text>Spot-Checking In this step, we will evaluate different algorithms. To do so, we will follow the following steps:1.  **Choosing a set of algorithms**: We will include some linear models seen in class in addition to other more complex models to make a more interesting comparison. Also, it is interesting to test several configurations for each algorithm by modifying the hyperparameters. It is not about searching through the hyperparameter space but allowing different models to behave well even if they are not optimal. That allows us to quickly recognize what kind of algorithms behave best for this problem to explore them further later. As we have a large sample, we have to be very critical when selecting the models, especially when evaluating more complex proposals such as RF or MLP. Knowing this, we have to reduce the set of distinct configurations using our knowledge about the problem, which we have acquired in the previous steps.2.  **Evaluate all selected models using cross-validation**. To evaluate the models, we will use a **k-fold** CV schema. That ensures that even using the training set, we can accurately estimate the error out-of-sample. There are different ways to perform cross-validation. In this case, we will use 5-fold CV mainly because we have a large sample and different models to evaluate. If we were to use 10-fold, for example, we would have to fit ten different models for each configuration, and for each fit, 90% of the data would be used. Our sample size allows us to discard 20% for evaluation while limiting the number of models to fit. Also, the partitions should be stratified.3.  **Interpret the results**. To do this, we will generate different box-plots for each configuration, and from there, select the best models based on their average score and their estimated variance.A crucial part of this analysis is selecting the metric with which we will evaluate the different algorithms. As we saw at the beginning of the EDA, the dataset is not balanced. For this reason, using a naive metric like accuracy can lead to an erroneous interpretation of the results, and therefore an incorrect selection of the final model. Since a priori, there are no discrepancies between false positives and false negatives; we will use ROC AUC to evaluate the models.The complete list of algorithms to try is the following:- **Linear**: Ridge, Lasso, and SVM models with a linear kernel (the latter will be explained in the SVM section below). Both Ridge and Lasso are implemented using SGD, as it is faster than the analytical solution due to the dataset's size. In the case of Ridge and Lasso, we have implemented three different configurations for each one, in which we modify the value of the hyperparameter $\alpha$. This parameter controls the level of regularization: $\alpha=0.1$ (little regularization), $\alpha=0.5$ (moderate regularization) and $\alpha=1$ (maximum regularization). We have added polynomial transformations to the grade 2 and 3 data for these last two $\alpha$ values, respectively. The idea is to increase the regularization as we increase the number of features, to have models balanced in power and complexity.- **Neural Networks**: 4 MLP architectures have been selected with the following topologies: _(20), (50, 20), (25, 25, 25)_. Each component of the tuple represents a hidden layer and the number of units for that layer. Therefore, we have three architectures of different sizes with 1 to 3 hidden layers and 20 to 50 units per layer. These architectures have been chosen based on the dimensionality of the input data and the training sample size. We could have selected architectures with many more units per layer (or more hidden layers). However, we have less than 30 features, making it unnecessary to use architectures with thousands of neurons. Every architecture uses *RELU* as activation function, which solves some inherent problems from the *sigmoid* and *tanh* functions, e.g., vanishing gradient. On the other hand, the training algorithm we will use is Adam, which has a fast convergence and few hyperparameters to modify. It usually works very well on medium to large datasets like this, providing good solutions in a fraction of the computing time. The hyperparameters of the training algorithm, specifically _beta1, beta2_, and the initial learning rate, will be the default ones. Later, we will adjust more precisely each of the parameters.- **SVM**: SVM family will be used as a linear and non-linear model given the flexibility provided by kernels' use. In the first case, we will fit a linear SVM (without kernel) with two different configurations in which we set the hyperparameter $ C $ to 0.5 and 1 to control the hyperplane margin. The higher the value of $C$, the more importance will be given to misclassified examples instead of achieving a wide margin between classes. The non-linear models have two variants depending on the type of kernel: polynomial (with grade 2) and Gaussian. For both cases, $C=1$. Both configurations are chosen based on their linear counterparts. In this case, we relax the margins ($C=1$), and that is why we have selected the grade 2 polynomial kernel, not the grade 3, as in the case of linear models.- **Random Forest**: For RF, we have added two configurations varying the number of estimators, 100 and 500 trees. In both cases, the Gini criterion is used to evaluate the impurity, while for the depth, no specific value has been specified. The nodes will expand as long as their leaf nodes are not pure or have more than two samples.- **AdaBoost**: This boosting model is composed of decision tree estimators with a maximum depth of 1. The number of iterations used is 100, but early stopping is applied.- **Gradient Boosting:** We have also considered it interesting to include a Gradient Boosting model. Gradient Boosting is a generalization of Boosting that uses gradient-based methods to adjust the model. Among this model's benefits, we can find its robustness against outliers, flexibility in working with heterogeneous data, and excellent predictive power . Apart from the difficulty in its parallelization, the main disadvantage of this algorithm is that it is quite complex to adjust its hyperparameters. For this reason, we have only modified the number of trees or rounds.- **Naive Bayes:** We are going to make a comparison with a probability-based model. We will use the Gaussian variant because, despite used One-Hot encoding, we also have continuous variables, which makes Bernouilli's version unfeasible. As can be seen, the assumption that $P(x_i \mid y) = \frac{1}{\sqrt{2\pi\sigma^2_y}} \exp\left(-\frac{(x_i - \mu_y)^2}{2\sigma^2_y}\right)$ is quite strong (and false in many cases). Despite this, this model tends to give good results in some cases and is pretty fast to fit.The following sections show the results of the evaluation of the algorithms divided into groups. These groups have been created based on the family or nature of each algorithm.  Linear Models<jupyter_code>linear_models = get_linear_models()
linear_results = evaluate_models(X_train,
                                 y_train,
                                 linear_models,
                                 metric='roc_auc',
                                 cv=5,
                                 verbose=2)
summarize_results(linear_results)<jupyter_output>Evaluating:  Logistic-Lasso-0.1<jupyter_text>The figure X, shows how SVM and Ridge algorithms have very similar metrics that never go below 0.85 for any configuration. On the other hand, Lasso is quite behind with much worse accuracy in any of its versions. In choosing between SVM and Ridge, we will choose SVM due to its greater generalization and optimality capacity.  Non-lineal Models<jupyter_code>non_linear_models = get_non_linear_models()
non_linear_results = evaluate_models(X_train,
                                     y_train,
                                     non_linear_models,
                                     metric='roc_auc',
                                     cv=5,
                                     verbose=2)
summarize_results(non_linear_results)<jupyter_output>Evaluating:  SVM-PolynomialKernel<jupyter_text>Among non-linear models, the best configuration is KNN with 32 neighbors. For this reason, we will use KNN and find the optimal number of neighbors.  Ensembles<jupyter_code>ensemble_models = get_ensemble_models()
ensemble_results = evaluate_models(X_train,
                                   y_train,
                                   ensemble_models,
                                   metric='roc_auc',
                                   cv=5,
                                   verbose=2)
summarize_results(ensemble_results)<jupyter_output>Evaluating:  AdaBoost-100<jupyter_text>As depicted in the figure above, the results are quite favorable for boosting algorithms, where both configurations obtain very similar results. Random Forest is quite behind, although we have two configurations with different numbers of trees. Therefore, we will discard it as a possible final model. Gradient Boosting seems to give better results with the same number of trees than AdaBoost, so we have decided to choose it as a candidate for the final model.  Neural Networks<jupyter_code>nn_models = get_neural_network_models()
nn_results = evaluate_models(X_train, y_train, nn_models, metric='roc_auc', cv=5, verbose=2)
summarize_results(nn_results)<jupyter_output>Rank=1, Name=mlp-20, Score=0.891 (+/- 0.004)
Rank=2, Name=mlp-25-25-25, Score=0.885 (+/- 0.003)
Rank=3, Name=mlp-50-20, Score=0.882 (+/- 0.004)<jupyter_text>The best MLP architecture a priori is the most simple, only one hidden layer of 20 units. Besides being the architecture with the best average AUC, it is also the one with the smallest VC dimension, so the generalization will be better than in the rest of the architectures. In the next step, we will find the optimum number of units for the hidden layer, and adjust some values of the training algorithm. To finish this section of Spot-Checking, we will show the results of the models selected so far:<jupyter_code>all_results = {}
all_results.update(nn_results)
all_results.update(ensemble_results)
all_results.update(linear_results)
all_results.update(non_linear_results)
top_results = {k: all_results[k] for k in ('LinearSVM-0.5', 'KNN-32', 'GradientBoosting-100', 'mlp-20')}
summarize_results(top_results)<jupyter_output>Rank=1, Name=mlp-20, Score=0.891 (+/- 0.004)
Rank=2, Name=GradientBoosting-100, Score=0.891 (+/- 0.003)
Rank=3, Name=LinearSVM-0.5, Score=0.880 (+/- 0.004)
Rank=4, Name=KNN-32, Score=0.880 (+/- 0.004)<jupyter_text>This plot allows us to take perspective and see which algorithms start with an advantage for the next step, although it is worth noting that the advantage of some over others is of hundredths. We will now adjust each model's hyperparameters and see if the trend is maintained or if some model outperforms the others.  Final Model In the previous step, we selected one model of each type. Now we are going to find the optimal hyperparameters for each model and compare them. Finally, we will pick the final model based on the selected metrics and the generalization capacity. That is, we will look for a balance between bias and variance. We will also take into account other factors such as interpretability or prediction time.To adjust the hyperparameters, we will do it in two ways. We will use Grid Search for the algorithms that are not very expensive to train. While for more expensive algorithms such as neural networks, we will use Random Search. Random search allows us to find near-optimal solutions with much less effort. We started with SVM linear. To find an optimal model based on this algorithm, we will adjust the parameter $C$. As we know, it is in charge of controlling the margin size. As it is the most important parameter of this algorithm, we will Random Search for the value of $C$ between 0 and 1.<jupyter_code># TODO: Modelo lineal
linear_param_space = {
    'predictor__C': np.linspace(0, 1, 101)
}

best_linear = tune_hyperparameters(
    X_train,y_train,
    LinearSVC(random_state=SEED),
    linear_param_space,
    exhaustive=False,
    n_iter=20,
    verbose=1,
    random_state=1)

render_optimized_model_results(best_linear)<jupyter_output>Fitting 5 folds for each of 20 candidates, totalling 100 fits<jupyter_text>The best predictor of linear SVM is found with a relatively small $C$, which indicates wide margins. That will make the model have an excellent generalization level. Furthermore, we are not using any kernel that projects into a higher space, nor are we using any polynomial transformations.Now we move on to optimize the KNN model. For this, we are going to find an optimal number of neighbors. As fitting a KNN is slow and the sample we have is quite large, the amount of configurations we will try is going to be limited. Specifically, we will try seven configurations, with the following number of neighbors: 32, 40, 48, 56, 64, 72, 80. These numbers have been selected based on the previous step, where we saw that increasing the number of neighbors increased the ROC-AUC metric.<jupyter_code>non_linear_param_space = {
    'predictor__n_neighbors': np.arange(32, 88, 8)
}

best_non_linear = tune_hyperparameters(
    X_train,y_train,
    KNeighborsClassifier(),
    non_linear_param_space,
    exhaustive=True,
    verbose=1)

render_optimized_model_results(best_non_linear)<jupyter_output>Fitting 5 folds for each of 7 candidates, totalling 35 fits<jupyter_text>The best model we have found out of 7 candidates is a configuration with 72 neighbors - more than twice as many as in the model from the previous step. The main reason is that, as with the Bagging models, for example, by increasing the number of neighbors, the variance in the model is reduced.Now we move on to the ensembles. Remember that we have selected GradientBoosting. This type of ensemble, like the rest, usually uses decision trees as weak classifiers. So, we will focus on optimizing the trees' properties that compose the ensemble besides the number of estimators. Specifically, we will optimize the trees' depth, one of the most important hyperparameters in the bias-variance tradeoff. We will also optimize the percentage of samples that each classifier uses for training and the learning rate, which controls each tree's contribution to the total ensemble. Finally, the search space would be as follows:- Number of trees: $\{50, 150, \dots, 500\}$- Learning rate: $[0.1, 0.9)$- Max depth: $\{2, 3, \dots, 10\}$- Subsample ratio (ratio of the training set used to fit each tree): $[0.6, 1)$Due to time constraints, and because this model with 100 estimators gives good results, hyperparameter space has been limited to 500 estimators. Therefore, for the learning rate $\alpha$, we do not have any knowledge about which values are optimal, so we tried out a large space. The maximum depth is limited to 10 in order to avoid having estimators with much variance. Finally, the sample has been selected with a minimum of 60% to have enough samples to train each estimator. The rest of the hyperparameters that can be adjusted have are the default values.<jupyter_code>from sklearn.ensemble import GradientBoostingClassifier

ensemble_param_space = {
    'predictor__n_estimators' : np.arange(100, 500, 50),
    'predictor__learning_rate' : np.arange(0.1, 1, 0.1),
    'predictor__max_depth' : np.arange(2, 10, 1),
    'predictor__subsample' : np.arange(0.6, 1, 0.1)
}

best_ensemble = tune_hyperparameters(
    X_train, y_train,
    GradientBoostingClassifier(random_state=SEED),
    ensemble_param_space,
    exhaustive=False,
    n_iter=20,
    verbose=1,
    random_state=1)

render_optimized_model_results(best_ensemble)<jupyter_output>Fitting 5 folds for each of 20 candidates, totalling 100 fits<jupyter_text>As we can see, these parameters found for the ensemble make it have an excellent bias-variance relationship. The number of iterations makes the bias decrease, while the trees of depth six can increase the variance, but by using 70% of the data in each sample, we manage to reduce it.In the case of neural networks, let us remember that we have an architecture composed of a single hidden layer. For this reason, we are going to focus mainly on finding an optimal number of units for that layer. Also, we will optimize the parameters $\beta_1$ and $\beta_2$ of the Adam, as well as the regularization parameter $\alpha$. The search space is the following:- Hidden units $\in \{5,6,\dots, 40\}$- $\beta_1 \in [0.5, 0.99]$- $\beta_2 \in [0.9, 0.99]$- $\alpha \in [0, 1]$We have tried to select a wide interval of values, ranging from a few neurons to almost twice as many neurons as the input dimension. According to the literature and documentation of Sklearn, the values of $\beta_1$ and $\beta_2$ are recommended to be close to 1. We have no prior knowledge of the $\alpha$ parameter, so we choose a wide range of values.<jupyter_code>hidden_layer_sizes = list(map(lambda x: (x,), range(5, 40)))

nn_param_grid = {
    'predictor__hidden_layer_sizes': hidden_layer_sizes,
    'predictor__beta_1': np.linspace(0.8, 0.99, 10),
    'predictor__beta_2': np.linspace(0.8, 0.999, 10),
    'predictor__alpha': np.linspace(0, 1, 11)
}

best_nn = tune_hyperparameters(
    X_train,y_train,
    MLPClassifier(random_state=1),
    nn_param_grid,
    exhaustive=False,
    n_iter=20,
    verbose=1,
    random_state=1)

render_optimized_model_results(best_nn)<jupyter_output>Fitting 5 folds for each of 20 candidates, totalling 100 fits<jupyter_text>The optimal architecture for MLP contains a hidden layer of 17 units, and as we mentioned earlier, $\beta_1$ and $\beta_2$ high. The model also requires a weak regularization, mainly because we are using a simple architecture.  ComparisonOnce we have found some optimal models for each category, we will compare them following different criteria to select the final model. We will start with the primary metric, ROC AUC. In general, all models have a similar AUC, between 0.88 and 0.9. Although the most performant model is MLP, the differences with Gradient Boosting are minimal.As we know, the generalization power of neural networks is not very high for medium-sized datasets. On the other hand, the best Gradient Boosting model is composed of trees of depth 3, which makes it a good option in terms of generalization, since it is composed of estimators with low variance and the bias has been reduced through the iterations. Finally, we have the SVM and KNN models, which have a similar metric close to 0.88. In this sense, SVM's generalization power is higher than KNN's, and the cost of making predictions is much lower since it is a parametric model.ROC AUC gives us a good estimate of the model's behavior, but we will also analyze the confusion matrices generated on the training data. Although it is hard to extrapolate from the training sample itself, we can see what kind of errors are more typical for each classifier.<jupyter_code>show_models_confusion_matrix(X_train, y_train, [best_linear, best_non_linear, 
                              best_nn, best_ensemble])<jupyter_output><empty_output><jupyter_text>In the image above, we can see the confusion matrix for each of the models. As we can see, for the positive class (> 50K), the SVM sorter behaves hardly better than a random selection. While for the ensemble and neural network, there is a higher proportion of examples correctly classified. Remember that we do not have a preference for any error type (Type I or Type II) because we do not know what use the model would have. In any case, both MLP and Gradient Boosting offer a sweet balance.Although we know that these confusion matrices have been generated from the training data, we also know that these two models we mentioned have the best metric calculated via CV. For this reason, either model generalizes quite well.Finally, after observing different points of view - metrics, confusion matrix,  and power of generalization of each model, we have decided to use the ensemble as the final model. The reason is that it offers a high metric and generalization power. Furthermore, it offers out-of-the-box the relative importance of each feature.<jupyter_code>final_model = best_ensemble.best_estimator_<jupyter_output><empty_output><jupyter_text>Results Once we have the final model, we can calculate the metric with the test sample. As we have not used the test sample for any procedure that might contaminate the results, we have made sure that we will calculate an unbiased estimator.The following table shows the results of the final model on the test set:<jupyter_code>show_test_results(test.iloc[1:, :-1], test.iloc[1:, -1], final_model)<jupyter_output>ROC AUC Score: 0.889038

              precision    recall  f1-score   support

           0       0.91      0.88      0.89     12916
           1       0.59      0.68      0.63      3365

   micro avg       0.84      0.84      0.84     16281
   macro avg       0.75      0.78      0.76     16281
weighted avg       0.85      0.84      0.84     16281<jupyter_text>The ROC AUC metric is very similar to the one found on previous steps when using cross-validation. That indicates that our model generalizes quite well. On the other hand, we see a breakdown of different metrics for each class. It is worth noting how the model can predict very accurately negative examples, while for the other class, it is slightly better than a random selection. As we said before, in this case, we do not consider it a problem. If we wanted to balance both accuracies and pay more attention to class 1, we could have reweighted the positive samples. However, the result is acceptable for our purposes.  Interpretation For the final part of our project, it would be interesting to assess the importance of each feature for the predictions. As we are using a tree-based ensemble model, it is quite simple to calculate it.<jupyter_code>plot_feature_importances(final_model)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0604897979521164, 'fraction_numerical': 0.014674007822815716, 'mean_word_length': 4.548030407740152, 'pattern_counts': {'":': 0, '<': 109, '<?xml version=': 0, '>': 110, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '29', 'dir': 'jupyter-structured-clean-dedup', 'id': '159543', 'license': 'MIT', 'n_tokens_mistral': 10629, 'n_tokens_neox': 9853, 'n_words': 5490, 'path': 'adult_dataset/notebook.ipynb', 'repo_name': 'antoniomdk/data_science_personal_projects', 'size': '40139'}