Introduction

In this project, I will apply several supervised learning algorithms of my choice to accurately model individuals' income using data collected from the 1994 U.S. Census. Then I will choose the best candidate algorithm from preliminary results and further optimize this algorithm to model the data in a best way possible. My goal with this implementation is to construct a model that accurately predicts whether an individual makes more than $50,000. This sort of task can arise in a non-profit setting, where organizations survive on donations. Understanding an individual's income can help a non-profit better understand how large of a donation to request, or whether or not they should reach out to begin with. While it can be difficult to determine an individual's general income bracket directly from public sources, we can (as we will see) infer this value from other publicly available features.

The dataset for this project originates from the UCI Machine Learning Repository. The dataset was donated by Ron Kohavi and Barry Becker, after being published in the article "Scaling Up the Accuracy of Naive-Bayes Classifiers: A Decision-Tree Hybrid". You can find the article by Ron Kohavi online. The data I will be investigating here consists of small changes to the original dataset, such as removing the 'fnlwgt' feature and records with missing or ill-formatted entries.


Part 1: Importing Libraries

In [1]:
# Import libraries necessary for this project
import numpy as np
import pandas as pd
from time import time
from IPython.display import display # Allows the use of display() for DataFrames

# Import visualisation libraries
import matplotlib.pyplot as plt
import seaborn as sns
import visuals as vs

# Allow pretty display for Jupyter notebooks
%matplotlib inline

# Disable warnings
import warnings
warnings.filterwarnings("ignore")

Part 2: Importing Dataset

In [2]:
data = pd.read_csv("census.csv")
test_data = pd.read_csv("test_census.csv")

Part 3: Data Exploration

In [12]:
data.head()
----------------------------------------
Out[12]:
Unnamed: 0 age workclass education_level education-num marital-status occupation relationship race sex capital-gain capital-loss hours-per-week native-country
0 0 21.0 Private 10th 6.0 Married-civ-spouse Craft-repair Husband White Male 0.0 0.0 40.0 United-States
1 1 49.0 Private Bachelors 13.0 Married-civ-spouse Adm-clerical Wife White Female 0.0 0.0 40.0 United-States
2 2 44.0 Self-emp-not-inc Assoc-acdm 12.0 Married-civ-spouse Other-service Wife White Female 0.0 0.0 99.0 United-States
3 3 34.0 Private Bachelors 13.0 Married-civ-spouse Sales Husband White Male 7298.0 0.0 46.0 United-States
4 4 24.0 Private HS-grad 9.0 Married-civ-spouse Machine-op-inspct Husband White Male 0.0 0.0 40.0 United-States
In [13]:
test_data.head()
Out[13]:
Unnamed: 0 age workclass education_level education-num marital-status occupation relationship race sex capital-gain capital-loss hours-per-week native-country
0 0 21.0 Private 10th 6.0 Married-civ-spouse Craft-repair Husband White Male 0.0 0.0 40.0 United-States
1 1 49.0 Private Bachelors 13.0 Married-civ-spouse Adm-clerical Wife White Female 0.0 0.0 40.0 United-States
2 2 44.0 Self-emp-not-inc Assoc-acdm 12.0 Married-civ-spouse Other-service Wife White Female 0.0 0.0 99.0 United-States
3 3 34.0 Private Bachelors 13.0 Married-civ-spouse Sales Husband White Male 7298.0 0.0 46.0 United-States
4 4 24.0 Private HS-grad 9.0 Married-civ-spouse Machine-op-inspct Husband White Male 0.0 0.0 40.0 United-States
In [15]:
data.info()
print("-"*40)
test_data.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 45222 entries, 0 to 45221
Data columns (total 14 columns):
age                45222 non-null int64
workclass          45222 non-null object
education_level    45222 non-null object
education-num      45222 non-null float64
marital-status     45222 non-null object
occupation         45222 non-null object
relationship       45222 non-null object
race               45222 non-null object
sex                45222 non-null object
capital-gain       45222 non-null float64
capital-loss       45222 non-null float64
hours-per-week     45222 non-null float64
native-country     45222 non-null object
income             45222 non-null object
dtypes: float64(4), int64(1), object(9)
memory usage: 4.8+ MB
----------------------------------------
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 45222 entries, 0 to 45221
Data columns (total 14 columns):
Unnamed: 0         45222 non-null int64
age                45205 non-null float64
workclass          45200 non-null object
education_level    45202 non-null object
education-num      45208 non-null float64
marital-status     45201 non-null object
occupation         45200 non-null object
relationship       45206 non-null object
race               45203 non-null object
sex                45203 non-null object
capital-gain       45207 non-null float64
capital-loss       45205 non-null float64
hours-per-week     45209 non-null float64
native-country     45206 non-null object
dtypes: float64(5), int64(1), object(8)
memory usage: 4.8+ MB
In [16]:
data.describe()
Out[16]:
age education-num capital-gain capital-loss hours-per-week
count 45222.000000 45222.000000 45222.000000 45222.000000 45222.000000
mean 38.547941 10.118460 1101.430344 88.595418 40.938017
std 13.217870 2.552881 7506.430084 404.956092 12.007508
min 17.000000 1.000000 0.000000 0.000000 1.000000
25% 28.000000 9.000000 0.000000 0.000000 40.000000
50% 37.000000 10.000000 0.000000 0.000000 40.000000
75% 47.000000 13.000000 0.000000 0.000000 45.000000
max 90.000000 16.000000 99999.000000 4356.000000 99.000000

A cursory investigation of the dataset will determine how many individuals fit into either group, and will tell us about the percentage of these individuals making more than $50,000.

Let's take a look at the following:

  1. The total number of records, "n_records".
  2. The number of individuals making more than \$50,000 annually, "n_greater_50k".
  3. The number of individuals making at most \$50,000 annually, "n_at_most_50k".
  4. The percentage of individuals making more than \$50,000 annually, "greater_percent".
In [17]:
# Total number of records
n_records = len(data)

# Number of records where individual's income is more than $50,000
n_greater_50k = len(data[data["income"] == ">50K"])

# Number of records where individual's income is at most $50,000
n_at_most_50k = len(data[data["income"] == "<=50K"])

# Percentage of individuals whose income is more than $50,000
greater_percent = 100 * n_greater_50k / n_records

# Print the results
print("Total number of records: {}".format(n_records))
print("Individuals making more than $50,000: {}".format(n_greater_50k))
print("Individuals making at most $50,000: {}".format(n_at_most_50k))
print("Percentage of individuals making more than $50,000: {:.2f}%".format(greater_percent))
Total number of records: 45222
Individuals making more than $50,000: 11208
Individuals making at most $50,000: 34014
Percentage of individuals making more than $50,000: 24.78%

Thanks to Seaborn, it is also possible to visualize the relationship between different features of a given individual and their incomes.

Let's see the breakdown of the counts of people earning more or less than \$50,000 annually, based on their gender and education levels.

In [18]:
sns.set(style="whitegrid", color_codes=True)
sns.factorplot("sex", col="education_level", data=data, hue="income", kind="count", col_wrap=4);
In [19]:
# Check how many unique values each feature has
for column in data.columns:
    print(column, len(data[column].unique()))
age 74
workclass 7
education_level 16
education-num 16
marital-status 7
occupation 14
relationship 6
race 5
sex 2
capital-gain 121
capital-loss 97
hours-per-week 96
native-country 41
income 2

Age, capital_gain, capital_loss, hours_per_week and education_num are continuous features and the rest are categorical features.

In [20]:
continuous = ["age", "capital-gain", "capital-loss", "hours-per-week", "education-num"]
categorical = ["workclass", "education_level", "marital-status", "occupation", "relationship", "race", "sex", "native-country"]
In [21]:
# For each categorical features print the unique values:
for column in categorical:
    print(column, data[column].unique())
workclass [' State-gov' ' Self-emp-not-inc' ' Private' ' Federal-gov' ' Local-gov'
 ' Self-emp-inc' ' Without-pay']
education_level [' Bachelors' ' HS-grad' ' 11th' ' Masters' ' 9th' ' Some-college'
 ' Assoc-acdm' ' 7th-8th' ' Doctorate' ' Assoc-voc' ' Prof-school'
 ' 5th-6th' ' 10th' ' Preschool' ' 12th' ' 1st-4th']
marital-status [' Never-married' ' Married-civ-spouse' ' Divorced'
 ' Married-spouse-absent' ' Separated' ' Married-AF-spouse' ' Widowed']
occupation [' Adm-clerical' ' Exec-managerial' ' Handlers-cleaners' ' Prof-specialty'
 ' Other-service' ' Sales' ' Transport-moving' ' Farming-fishing'
 ' Machine-op-inspct' ' Tech-support' ' Craft-repair' ' Protective-serv'
 ' Armed-Forces' ' Priv-house-serv']
relationship [' Not-in-family' ' Husband' ' Wife' ' Own-child' ' Unmarried'
 ' Other-relative']
race [' White' ' Black' ' Asian-Pac-Islander' ' Amer-Indian-Eskimo' ' Other']
sex [' Male' ' Female']
native-country [' United-States' ' Cuba' ' Jamaica' ' India' ' Mexico' ' Puerto-Rico'
 ' Honduras' ' England' ' Canada' ' Germany' ' Iran' ' Philippines'
 ' Poland' ' Columbia' ' Cambodia' ' Thailand' ' Ecuador' ' Laos'
 ' Taiwan' ' Haiti' ' Portugal' ' Dominican-Republic' ' El-Salvador'
 ' France' ' Guatemala' ' Italy' ' China' ' South' ' Japan' ' Yugoslavia'
 ' Peru' ' Outlying-US(Guam-USVI-etc)' ' Scotland' ' Trinadad&Tobago'
 ' Greece' ' Nicaragua' ' Vietnam' ' Hong' ' Ireland' ' Hungary'
 ' Holand-Netherlands']

3.1. Summary of the Data

  • age: continuous.
  • workclass: Private, Self-emp-not-inc, Self-emp-inc, Federal-gov, Local-gov, State-gov, Without-pay, Never-worked.
  • education: 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: continuous.
  • marital-status: Married-civ-spouse, Divorced, Never-married, Separated, Widowed, Married-spouse-absent, Married-AF-spouse.
  • occupation: 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: Wife, Own-child, Husband, Not-in-family, Other-relative, Unmarried.
  • race: Black, White, Asian-Pac-Islander, Amer-Indian-Eskimo, Other.
  • sex: Female, Male.
  • capital-gain: continuous.
  • capital-loss: continuous.
  • hours-per-week: continuous.
  • native-country: United-States, Cambodia, England, Puerto-Rico, Canada, Germany, Outlying-US(Guam-USVI-etc), India, Japan, Greece, South, China, Cuba, Iran, Honduras, Philippines, Italy, Poland, Jamaica, Vietnam, Mexico, Portugal, Ireland, France, Dominican-Republic, Laos, Ecuador, Taiwan, Haiti, Columbia, Hungary, Guatemala, Nicaragua, Scotland, Thailand, Yugoslavia, El-Salvador, Trinadad&Tobago, Peru, Hong, Holand-Netherlands.

Part 4: Preparing the Data

Before data can be used as input for machine learning algorithms, it often must be cleaned, formatted, and restructured — this process is typically known as preprocessing. Fortunately, for this dataset, there are no invalid or missing entries we must deal with, however, there are some qualities about certain features that must be adjusted. This preprocessing can help tremendously with the outcome and predictive power of nearly all learning algorithms.

4.1. Transforming Skewed Continuous Features

A dataset may sometimes contain at least one feature whose values tend to lie near a single number, but will also have a non-trivial number of vastly larger or smaller values than that single number. Algorithms can be sensitive to such distributions of values and can underperform if the range is not properly normalized. With the census dataset two features fit this description: "capital-gain" and "capital-loss".

Let’s plot a histogram of these two features and see how they are distributed.

In [22]:
# Split the data into features and target label
income_raw = data["income"]
features_raw = data.drop("income", axis=1)

# Visualize skewed continuous features of original data
vs.distribution(data)

For highly-skewed feature distributions such as "capital-gain" and "capital-loss", it is a common practice to apply a logarithmic transformation on the data so that the very large and very small values do not negatively affect the performance of a learning algorithm. Using a logarithmic transformation significantly reduces the range of values caused by outliers. However, we must be very careful when we are applying this transformation because the logarithm of 0 is undefined, so we must translate the values by a small amount above 0 to be able to apply the logarithm successfully.

In [23]:
# Log-transform the skewed features
skewed = ['capital-gain', 'capital-loss']
features_log_transformed = pd.DataFrame(data = features_raw)
features_log_transformed[skewed] = features_raw[skewed].apply(lambda x: np.log(x + 1))

# Visualize the new log distributions
vs.distribution(features_log_transformed, transformed=True)

4.2. Normalizing Numerical Features

In addition to performing transformations on features that are highly skewed, it is often good practice to perform some type of scaling on numerical features as well. Applying a scaling to the data does not change the shape of each feature's distribution (such as "capital-gain" or "capital-loss" above); however, normalization ensures that each feature is treated equally when applying supervised learners. Note that once scaling is applied, observing the data in its raw form will no longer have the same original meaning, as you will see below.

I will use sklearn.preprocessing.MinMaxScaler for this step.

In [13]:
# Import sklearn.preprocessing.MinMaxScaler
from sklearn.preprocessing import MinMaxScaler

# Initialize a scaler, then apply it to the features
scaler = MinMaxScaler()
numerical = ["age", "education-num", "capital-gain", "capital-loss", "hours-per-week"]

features_log_minmax_transform = pd.DataFrame(data = features_log_transformed)
features_log_minmax_transform[numerical] = scaler.fit_transform(features_log_transformed[numerical])

# Show an example of a record with scaling applied
display(features_log_minmax_transform.head(n = 5))
age workclass education_level education-num marital-status occupation relationship race sex capital-gain capital-loss hours-per-week native-country
0 0.301370 State-gov Bachelors 0.800000 Never-married Adm-clerical Not-in-family White Male 0.667492 0.0 0.397959 United-States
1 0.452055 Self-emp-not-inc Bachelors 0.800000 Married-civ-spouse Exec-managerial Husband White Male 0.000000 0.0 0.122449 United-States
2 0.287671 Private HS-grad 0.533333 Divorced Handlers-cleaners Not-in-family White Male 0.000000 0.0 0.397959 United-States
3 0.493151 Private 11th 0.400000 Married-civ-spouse Handlers-cleaners Husband Black Male 0.000000 0.0 0.397959 United-States
4 0.150685 Private Bachelors 0.800000 Married-civ-spouse Prof-specialty Wife Black Female 0.000000 0.0 0.397959 Cuba

4.3. One-hot Encoding

In our dataset, there are several features for each record that are non-numeric. Typically, learning algorithms expect numeric input, which requires that non-numeric features (called categorical variables) to be converted. One popular way to convert categorical variables is by using the one-hot encoding scheme. One-hot encoding creates a "dummy" variable for each possible category of each non-numeric feature. For example, assume someFeature has three possible entries: A, B, or C. We then encode this feature into someFeature_A, someFeature_B and someFeature_C.

Additionally, as with the non-numeric features, we need to convert the non-numeric target label, "income" to numerical values for the learning algorithm to work. Since there are only two possible categories for this label ("<=50K" and ">50K"), we can avoid using one-hot encoding and simply encode these two categories as 0 and 1, respectively.

In [14]:
# One-hot encode the 'features_log_minmax_transform' data using pandas.get_dummies()
features_final = pd.get_dummies(features_log_minmax_transform)

# Encode the 'income_raw' data to numerical values
income = income_raw.apply(lambda x: 1 if x == '>50K' else 0)

# Print the number of features after one-hot encoding
encoded = list(features_final.columns)
print("{} total features after one-hot encoding.".format(len(encoded)))
103 total features after one-hot encoding.

4.4. Shuffle and Split the Dataset

Now that we have all of our categorical variables have been converted into numerical features, and all of our numerical features have been normalized. As always, we will now split the data (both the features and their labels) into training and testing sets. I will use 80% of the data for training and 20% for testing.

In [15]:
# Import train_test_split
from sklearn.model_selection import train_test_split

# Split the 'features' and 'income' data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(features_final, income, test_size = 0.2, random_state = 0)

# Show the results of the split
print("Training set has {} samples.".format(X_train.shape[0]))
print("Testing set has {} samples.".format(X_test.shape[0]))
Training set has 36177 samples.
Testing set has 9045 samples.

Part 5: Evaluating Model Performance

In this section, I will investigate four different algorithms, and determine which one is best for my use case. Three of these algorithms will be supervised learners, and the fourth algorithm will be a naive predictor.

CharityML, equipped with their research, knows individuals that make more than $50,000 are most likely to donate to their charity. Because of this, CharityML is particularly interested in predicting who makes more than $50,000 accurately. It would seem that using accuracy as a metric for evaluating a particular model's performace would be appropriate. Additionally, identifying someone that does not make more than $50,000 as someone who does would be detrimental to CharityML, since they are looking to find individuals willing to donate. Therefore, a model's ability to precisely predict those who make more than $50,000 is more important than the model's ability to recall those individuals. We can use F-beta score as a metric that considers both precision and recall:

$$ F_{\beta} = (1 + \beta^2) \cdot \frac{precision \cdot recall}{\left( \beta^2 \cdot precision \right) + recall} $$

In particular, when $\beta = 0.5$, more emphasis is placed on precision. This is called the F$_{0.5}$ score (or F-score for simplicity).

Looking at the distribution of classes (those who make at most $50,000, and those who make more), it's clear that most individuals do not make more than $50,000. This can greatly affect accuracy, since we could simply say "this person does not make more than \$50,000" and generally be right, without ever looking at the data! Making such a statement would be called naive, since we have not considered any information to substantiate the claim. It is always important to consider the naive prediction for our data, to help establish a benchmark for whether a model is performing well. That being said, using that prediction would be pointless: If we predicted all people made less than $50,000, CharityML would identify no one as donors.

5.1. A Quick Recap of Accuracy, Precision and Recall

Accuracy measures how often the classifier makes the correct prediction. It is the ratio of the number of correct predictions to the total number of predictions (the number of test data points).

Precision tells us what proportion of messages we classified as spam, actually were spam. It is a ratio of true positives (words classified as spam, and which are actually spam) to all positives (all words classified as spam, irrespective of whether that was the correct classification), in other words it is the ratio of: [True Positives/(True Positives + False Positives)]

Recall (or sensitivity) tells us what proportion of messages that actually were spam were classified by us as spam. It is a ratio of true positives (words classified as spam, and which are actually spam) to all the words that were actually spam, in other words it is the ratio of: [True Positives/(True Positives + False Negatives)]

For classification problems that are skewed in their classification distributions like in our case, for example if we had a 100 text messages and only 2 were spam and the rest 98 weren't, accuracy by itself is not a very good metric. We could classify 90 messages as not spam(including the 2 that were spam but we classify them as not spam, hence they would be false negatives) and 10 as spam(all 10 false positives) and still get a reasonably good accuracy score. For such cases, precision and recall come in very handy. These two metrics can be combined to get the F1 score, which is weighted average (harmonic mean) of the precision and recall scores. This score can range from 0 to 1, with 1 being the best possible F1 score (we take the harmonic mean as we are dealing with ratios).

5.2. Naive Predictor Performance

The purpose of generating a naive predictor is simply to show what a base model without any intelligence would look like. In the real world, ideally our base model would be either the results of a previous model or could be based on a research paper upon which we are looking to improve. When there is no benchmark model set, getting a result better than random choice is a place we could start from.

What if we chose a model that always predicted an individual made more than $50,000, what would that model’s accuracy and F-score be on this dataset?

In [16]:
# Calculate accuracy
accuracy = n_greater_50k / n_records

# Calculate precision
precision = n_greater_50k / (n_greater_50k + n_at_most_50k)

# Calculate recall
recall = n_greater_50k / (n_greater_50k + 0)

# Calculate F-score using the formula above for beta = 0.5
fscore =  (1  + (0.5*0.5)) * ( precision * recall / (( 0.5*0.5 * (precision))+ recall))

# Print the results 
print("Naive Predictor: [Accuracy score: {:.4f}, F-score: {:.4f}]".format(accuracy, fscore))
Naive Predictor: [Accuracy score: 0.2478, F-score: 0.2917]

5.3. Supervised Learning Models

The following are some of the supervised learning models that are currently available in scikit-learn that we can choose from:

  • Gaussian Naive Bayes (GaussianNB)
  • Decision Trees
  • Ensemble Methods (Bagging, AdaBoost, Random Forest, Gradient Boosting)
  • K-Nearest Neighbors (KNeighbors)
  • Stochastic Gradient Descent Classifier (SGDC)
  • Support Vector Machines (SVM)
  • Logistic Regression

Now I will pick three supervised learning models that are listed above which are appropriate for our problem, and test them on the census data.

5.3.1. Decision Trees

title

  • Real world applications: Decision Trees, or in general, Classification and Regression Trees have a lot of real world application areas. Some of them are:

    • Variable Selection: The number of variables that are routinely monitored in clinical settings has increased dramatically with the introduction of electronic data storage. Many of these variables are of marginal relevance and, thus, should probably not be included in data mining exercises. Like stepwise variable selection in regression analysis, decision tree methods can be used to select the most relevant input variables that should be used to form decision tree models, which can subsequently be used to formulate clinical hypotheses and inform subsequent research.
    • Assessing the relative importance of variables: Once a set of relevant variables is identified, researchers may want to know which variables play major roles. Generally, variable importance is computed based on the reduction of model accuracy (or in the purities of nodes in the tree) when the variable is removed. In most circumstances the more records a variable have an effect on, the greater the importance of the variable.
    • Handling of missing values: A common - but incorrect - method of handling missing data is to exclude cases with missing values; this is both inefficient and runs the risk of introducing bias in the analysis. Decision tree analysis can deal with missing data in two ways: it can either classify missing values as a separate category that can be analyzed with the other categories or use a built decision tree model which set the variable with lots of missing value as a target variable to make prediction and replace these missing ones with the predicted value. [Source]
  • Strengths:

    • They are able to handle both categorical and numerical data.
    • They do not require an immense process of much data preprocessing. One example to this is that they can handle un-normalized data.
    • They are rather intuitive to us and how we make our decisions in a daily life setting. Because of that, they are easy to understand and interpret.
    • The number of hyper-parameters to be tuned is almost null.
  • Weaknesses:

    • High probability of overfitting because decision trees who are rather complex do not generalize well enough to the data.
    • Information gain in a decision tree with categorical variables gives a biased response for attributes with greater number of categories.
    • They are rather unstable because tiny variations in the data can result in a completely different decision tree. Because of that, they are usually packed with an ensemble, like Random Forests, to build robustness. [Source]
  • Candidacy:

    • Since our dataset consists of both categorical and numerical data, decision trees can come in really handy in our use case. Also, they are easy to interpret, so we are going to have more chance when it comes to understanding what is going on under the hood.

5.3.2. Support Vector Machines (SVMs)

title

  • Real world applications:

    • SVMs have a lot of application areas. Some great examples of their use cases are: image classification and image segmentation, i.e., face detection in an image. [Source]
  • Strengths:

    • They are incredibly effective in high dimensional spaces, or simply, when we have a lot of features in our data.
    • Their kernel functions often come in really handy. They can be used to adapt to different use cases, or can be completely customized if needed. Therefore, they are versatile.
  • Weaknesses:

    • They do not perform well with large datasets.
    • They have several key parameters that need to be set correctly to achieve the best classification results for any given problem. Parameters that may result in an excellent classification accuracy for problem A, may result in a poor classification accuracy for problem B. Therefore we may have to experiment with a number of different parameter settings in order to achieve a satisfactory result. The main parameters that we should experiment with are the SVM kernel type, the SVM type and the kernel-specific parameters (such as gamma, degree, nu, etc.). [Source]
  • Candidacy:

    • I included SVMs for this project because I want to incorporate their effectiveness given high dimensionality. Considering that we used dummy variables to turn our categorical features into numerical values, we have more than 100 features in our dataset now. But that should not be a problem for SVMs. Well, one can only hope. :)

5.3.3. AdaBoost

title

  • Real world applications:

  • Strengths:

    • Like AdaBoost, all ensemble methods are more robust than single estimators.
    • They have improved generalizability.
    • They are prone to overfitting.
  • Weaknesses:

    • It can be sensitive to noisy data and outliers.
  • Candidacy:

    • Ensemble methods are considered to be high quality classifiers, and AdaBoost is one of most popular boosting algorithms thanks to its practical success with applications in a wide variety of fields, such as biology, computer vision, and speech processing. Also, considering that we have a class imbalance in our dataset, we could really use the boosting methods' robustness.

Part 6: Creating a Training and Predicting Pipeline

To properly evaluate the performance of each model we have chosen, it is important that we create a training and predicting pipeline that allows us to quickly and effectively train models using various sizes of training data and perform predictions on the testing data.

In [17]:
# Import two metrics from sklearn: fbeta_score and accuracy_score
from sklearn.metrics import fbeta_score, accuracy_score

def train_predict(learner, sample_size, X_train, y_train, X_test, y_test):
    '''
    inputs:
       - learner: the learning algorithm to be trained and predicted on
       - sample_size: the size of samples (number) to be drawn from training set
       - X_train: features training set
       - y_train: income training set
       - X_test: features testing set
       - y_test: income testing set
    '''

    results = {}

    # Fit the learner to the training data using slicing with 'sample_size'
    start = time() # Get start time
    learner = learner.fit(X_train[:sample_size], y_train[:sample_size])
    end = time() # Get end time

    # Calculate the training time
    results['train_time'] = end - start

    # Get the predictions on the test set, then get predictions on the first 300 training samples
    start = time() # Get start time
    predictions_test = learner.predict(X_test)
    predictions_train = learner.predict(X_train[:300])
    end = time() # Get end time

    # Calculate the total prediction time
    results['pred_time'] = end - start

    # Compute accuracy on the first 300 training samples
    results['acc_train'] = accuracy_score(y_train[:300], predictions_train)

    # Compute accuracy on test set
    results['acc_test'] = accuracy_score(y_test, predictions_test)

    # Compute F-score on the the first 300 training samples
    results['f_train'] = fbeta_score(y_train[:300], predictions_train, 0.5)

    # Compute F-score on the test set
    results['f_test'] = fbeta_score(y_test,predictions_test, 0.5)

    # Success
    print ("{} trained on {} samples.".format(learner.__class__.__name__, sample_size))

    # Return the results
    return results

Part 7: Initial Model Evaluation

In [18]:
# Import the three supervised learning models we will be using
from sklearn.tree import DecisionTreeClassifier
from sklearn.svm import SVC
from sklearn.ensemble import AdaBoostClassifier

# Start timing
start = time()

# Initialize the models, the random states are set to 42 so we know how to reproduce the model later
clf_A = DecisionTreeClassifier(random_state = 42)
clf_B = SVC(random_state = 42)
clf_C = AdaBoostClassifier(random_state = 42)

# Calculate the number of samples for 1%, 10%, and 100% of the training data
samples_1 = int(round(len(X_train) / 100))
samples_10 = int(round(len(X_train) / 10))
samples_100 = len(X_train)

# Collect results on the learners
results = {}
for clf in [clf_A, clf_B, clf_C]:
    clf_name = clf.__class__.__name__
    results[clf_name] = {}
    for i, samples in enumerate([samples_1, samples_10, samples_100]):
        results[clf_name][i] = \
        train_predict(clf, samples, X_train, y_train, X_test, y_test)

# Run metrics visualization for the three supervised learning models chosen
vs.evaluate(results, accuracy, fscore)

# Display runtime
end = time()
runtime = end - start
print("Runtime:", runtime)
DecisionTreeClassifier trained on 362 samples.
DecisionTreeClassifier trained on 3618 samples.
DecisionTreeClassifier trained on 36177 samples.
SVC trained on 362 samples.
SVC trained on 3618 samples.
SVC trained on 36177 samples.
AdaBoostClassifier trained on 362 samples.
AdaBoostClassifier trained on 3618 samples.
AdaBoostClassifier trained on 36177 samples.
Runtime: 108.57590460777283

Let's print out the values in the visualizations above to examine the results in more detail.

In [19]:
# Printing out the values
for item in results.items():
    print(item[0])
    display(pd.DataFrame(item[1]).rename(columns={0:'1%', 1:'10%', 2:'100%'}))
DecisionTreeClassifier
1% 10% 100%
acc_test 0.753897 0.802322 0.818574
acc_train 1.000000 0.996667 0.970000
f_test 0.511395 0.595541 0.627948
f_train 1.000000 0.997191 0.963855
pred_time 0.016013 0.009005 0.009007
train_time 0.006003 0.026018 0.372834
SVC
1% 10% 100%
acc_test 0.756219 0.832615 0.837148
acc_train 0.760000 0.833333 0.853333
f_test 0.000000 0.671256 0.674477
f_train 0.000000 0.671296 0.719697
pred_time 0.195138 1.606148 14.265084
train_time 0.010507 0.758536 88.478011
AdaBoostClassifier
1% 10% 100%
acc_test 0.810282 0.849862 0.857601
acc_train 0.896667 0.840000 0.850000
f_test 0.610253 0.701882 0.724551
f_train 0.811688 0.680147 0.711538
pred_time 0.089563 0.088075 0.082070
train_time 0.049535 0.179626 1.640147

And also let's visualize the confusion matrices of the results above.

In [20]:
# Visualize the confusion matrix for each classifier
from sklearn.metrics import confusion_matrix

for i, model in enumerate([clf_A,clf_B,clf_C]):
    cm = confusion_matrix(y_test, model.predict(X_test))
    cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis] # normalize the data

    # Plot the heatmap
    plt.figure(i)
    sns.heatmap(cm, annot=True, annot_kws={"size":30}, cmap='Blues', square=True, fmt='.3f')
    plt.ylabel('True label')
    plt.xlabel('Predicted label')
    plt.title('Confusion matrix for:\n{}'.format(model.__class__.__name__));

Part 8: Choosing 'the' Model

Looking at the confusion matrices' heatmaps above, it is clear the AdaBoost is the model to go! First, it is by far the winner when it comes to performance on the testing data, in terms of both the accuracy and f-score, which is great. Second, comparing the SVM's whopping ~120 seconds runtime, it only takes a fraction of it. It is fast!

title

Also, AdaBoost uses a decision tree of depth 1 as its base classifier by default, which makes it really handy when it comes to handling both categorical and numerical data. Since the weak learners are faster to train, the size of our dataset shouldn't be a problem.

8.1. How does AdaBoost work?

AdaBoost, like any other boosting algorithm, does its magic by combining bunch of weak learners and makes a strong learner out of them. Let's create a real world scenario to demonstrate it better. Let's say that we want to create a Spam filter for our email service.

How would we classify an email as Spam or not? We would simply consider a couple of cases:

  1. Email only has one image file, which is highly likely to be a promotional image: SPAM
  2. Email has only link(s): SPAM
  3. Email body contains words like "Money, Win, Prize" etc.: SPAM
  4. Email body contains a lot of exclamation points: SPAM
  5. Email from my domain "www.mertozturk.org": NOT Spam
  6. Email from a known source: NOT Spam

Above, I have defined some rules to classify an email into spam or not spam. But, do we think these rules individually are strong enough to successfully classify an email? No.

Individually, these rules are not powerful enough to classify an email into spam or not spam. Therefore, these rules are called as weak learner.

To convert bunch of weak learners into a strong learner, we will combine the prediction of each weak learner using methods like:

  • Using average/ weighted average
  • Considering prediction has higher vote

For example: Above, we have defined 6 weak learners. Out of these 6, 4 are voted as SPAM and 2 are voted as Not SPAM. In this case, by default, we will consider an email as SPAM because we have higher vote for SPAM.

The AdaBoost algorithm will perform in exactly the same way:

  1. First it will combine several weak learners (in our case, decision trees) to create an ensemble of learners that can predict whether an individual earns more that 50K a year or not.
  2. Each of our weak learners, i.e., decision trees, have been created using the features we have about individuals in our dataset such as age, occupation, education, etc. Thanks to these weak learners, we now have a set of rules that we will use to predict an individual's income. (Since this project is on classification and not on regression, we will not predict the individual's income per se. We will try to classify whether the individual makes more than 50K a year or less.)
  3. During the training process, the AdaBoost algorithm looks at the cases where it predicted poorly, and prioritizes those cases to make a correct prediction on them in the next round of training.
  4. With each round, the algorithm finds the best weak learner to incorporate into the ensemble by repeating the process for the specified number of rounds, or until we can't improve the predictions any further.
  5. Then, all of the weak learners are combined to make a final ensembled model. After that, each of them vote to predict if a person earns more than 50K a year or not. Finally, we take the majority of the votes to make a final decision. Such a democratized model, huh!
  6. Using this process of ensembled learners on our census dataset, we can predict the same information for a potential new donor and predict if he/she earns more than 50K a year or not, and from there, we can make a decision on the probability of them donating to the CharityML or not.

Part 9: Optimizing our Model

title

In this step, I will be using grid search (GridSearchCV) with numerous different parameter/value combinations to tune our model for even better results.

I will tune the n_estimators and the learning_rate parameters for the AdaBoost. Given that our base classifier is a decision tree, I will also experiment with its parameters as well. Let's get started!

In [21]:
# Import 'GridSearchCV', 'make_scorer', and any other necessary libraries
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import make_scorer

# Start the timer
start = time()

# Initialize the classifier
clf = AdaBoostClassifier(base_estimator = DecisionTreeClassifier())

# Create the parameters list you wish to tune
parameters = {'n_estimators':[20, 50, 120],
              'learning_rate':[0.1, 0.5, 1.],
              'base_estimator__min_samples_split' : np.arange(2, 8, 2),
              'base_estimator__max_depth' : np.arange(1, 4, 1)
             }

# Make an fbeta_score scoring object
scorer = make_scorer(fbeta_score,beta=0.5)

# Perform grid search on the classifier using 'scorer' as the scoring method
grid_obj = GridSearchCV(clf, parameters, scorer)

# Fit the grid search object to the training data and find the optimal parameters
grid_fit = grid_obj.fit(X_train, y_train)

# Get the estimator
best_clf = grid_fit.best_estimator_

# Make predictions using the unoptimized and model
predictions = (clf.fit(X_train, y_train)).predict(X_test)
best_predictions = best_clf.predict(X_test)

# Report the before-and-afterscores
print("Unoptimized model\n------")
print("Accuracy score on testing data: {:.4f}".format(accuracy_score(y_test, predictions)))
print("F-score on testing data: {:.4f}".format(fbeta_score(y_test, predictions, beta = 0.5)))
print("\nOptimized Model\n------")
print("Final accuracy score on the testing data: {:.4f}".format(accuracy_score(y_test, best_predictions)))
print("Final F-score on the testing data: {:.4f}".format(fbeta_score(y_test, best_predictions, beta = 0.5)))
print(best_clf)

# Print runtime
end = time()
runtime = end - start
print("Runtime:", runtime)
Unoptimized model
------
Accuracy score on testing data: 0.8355
F-score on testing data: 0.6647

Optimized Model
------
Final accuracy score on the testing data: 0.8703
Final F-score on the testing data: 0.7529
AdaBoostClassifier(algorithm='SAMME.R',
          base_estimator=DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=3,
            max_features=None, max_leaf_nodes=None,
            min_impurity_decrease=0.0, min_impurity_split=None,
            min_samples_leaf=1, min_samples_split=6,
            min_weight_fraction_leaf=0.0, presort=False, random_state=None,
            splitter='best'),
          learning_rate=0.5, n_estimators=50, random_state=None)
Runtime: 669.4053871631622

Part 10: Final Model Evaluation

Results:

Metric Benchmark Predictor Unoptimized Model Optimized Model
Accuracy Score 0.2478 0.8355 0.8703
F-score 0.2917 0.6647 0.7529

As you can see above, our model has an accuracy of 0.8691 and F-score of 0.7494.

These scores are better than the unoptimized model, and also way better than the benchmark predictor.

Part 11: Feature Importance

An important task when performing supervised learning on a dataset like the census data we study here is determining which features provide the most predictive power. By focusing on the relationship between only a few crucial features and the target label we simplify our understanding of the phenomenon, which is most always a useful thing to do. In the case of this project, that means we wish to identify a small number of features that most strongly predict whether an individual makes at most or more than \$50,000.

11.1. Feature Relevance Observation

We know that there are thirteen available features for each individual on record in the census data. Of these thirteen records, we can guess which five features might be most important for prediction. Let’s dive right in.

If I were to guess what are the most important features for prediction, it would be:

  1. occupation: Every job has different paychecks. Some jobs pay significantly higher than others.
  2. education: People who have a higher level of education are usually better equipped to handle more technical/specialized jobs which happen to pay well.
  3. workclass: The working class that people belong to can also be correlated with how much money they make.
  4. age: Older people happen to accumulate a greater wealth.
  5. hours-per-week: If a person works more than usual, chances are he/she is making more money than average.

I ranked these features according to the impact I believe they have on an individual's income. I placed occupation as the most important feature as I believe it will be single-handedly the most correlated feature for our case. About the ordering of the other 4 features, I am not that certain.

11.2. Extracting Feature Importance

In [22]:
# Import a supervised learning model that has 'feature_importances_'
from sklearn.ensemble import AdaBoostClassifier

# Train the supervised model on the training set 
model = AdaBoostClassifier().fit(X_train,y_train)

# Extract the feature importances
importances = model.feature_importances_

# Plot
vs.feature_plot(importances, X_train, y_train)

Out of the 5 features that I predicted in the Section 11.1. above, 3 of them are included in the list of most important features: Age, hours-per-week, education-num.

I am quite surprised that occupation isn't in this list. And also I didn't consider capital-gain and capital-loss to be in this list, as I honestly didn't know what they meant. But now reading about them, it make sense that they are both in the list. People who have earned profits from sale of assets are definitely likelier to earn higher, while those who incurred losses are likely to have had lower total income.

11.3. Feature Selection

An interesting thing to think about here is that how does a model perform if we only use a subset of all the available features in the data? With less features required to train, the expectation is that training and prediction time is much lower, at the cost of performance metrics. From the visualization above, we see that the top five most important features contribute more than half of the importance of all features present in the data. This hints that we can attempt to reduce the feature space and simplify the information required for the model to learn. The code cell below will use the same optimized model we found earlier, and train it on the same training set with only the top five important features.

In [23]:
# Import functionality for cloning a model
from sklearn.base import clone

# Reduce the feature space
X_train_reduced = X_train[X_train.columns.values[(np.argsort(importances)[::-1])[:5]]]
X_test_reduced = X_test[X_test.columns.values[(np.argsort(importances)[::-1])[:5]]]

# Train on the "best" model found from grid search earlier
clf = (clone(best_clf)).fit(X_train_reduced, y_train)

# Make new predictions
reduced_predictions = clf.predict(X_test_reduced)

# Report scores from the final model using both versions of data
print("Final Model trained on full data\n------")
print("Accuracy on testing data: {:.4f}".format(accuracy_score(y_test, best_predictions)))
print("F-score on testing data: {:.4f}".format(fbeta_score(y_test, best_predictions, beta = 0.5)))
print("\nFinal Model trained on reduced data\n------")
print("Accuracy on testing data: {:.4f}".format(accuracy_score(y_test, reduced_predictions)))
print("F-score on testing data: {:.4f}".format(fbeta_score(y_test, reduced_predictions, beta = 0.5)))
Final Model trained on full data
------
Accuracy on testing data: 0.8703
F-score on testing data: 0.7529

Final Model trained on reduced data
------
Accuracy on testing data: 0.8437
F-score on testing data: 0.7065

11.4. Effects of Feature Selection

We can see that in the reduced dataset, the final model's accuracy and f-score values are still relatively good comparing to the full dataset.

The accuracy is ~2.6% lower, while the f-score is ~4.7% lower.

My conclusion is, depending on the person or the project, we might not always able to afford to get 2.6& or 4.7% less of an accuracy (or f-score) for the sake of reducing the training time. Also, AdaBoost is pretty fast, so I don't think I would choose to train my model with reduced data, instead, I would train it with the full data and wouldn't give up on accuracy or f-score. For a different model however, I could easily choose the other way around.

End

Thank you for your interest in this notebook. I completed this project on June 30, 2019 as a part of Udacity's Intro to Machine Learning Nanodegree.