Machine Learning and Modeling

What we will cover

In the previous chapters, we looked at the larger field of artificial intelligence which relates to automating intellectual tasks performed by humans. At one time, it was thought that all human decision making could be coded as a set of rules, which, if followed, would mimic intelligence. The idea was that while these rules could be extremely complex in terms of their length and count, and in the way these were nested with each other, but in the end a set of properly structured if-else rules held the key to creating an artificial mind.

Of course, we know now that is not accurate. Rule based systems cannot generalize from patterns like the human mind does, and tend to be brittle to the point that they can be practically unusable.

Machine learning algorithms attempt to identify patterns in the data with which they create a solution to solve problems that haven't been seen before. That is the topic for the discussion in this chapter.

Deep learning is a special case (or a subset) of machine learning where layers of data abstractions (called neural networks) are used. However, you may hear of a distinction being sometimes made between machine learning, also sometimes called 'shallow learning', from deep learning that we will cover in the next chapter. Machine learning is sometimes called 'shallow learning' because it is based on a single layer of data transformations. It is called so to distinguish it from 'deep learning' that relies upon multiple layers of data transformations, with each layer extracting a different elements of useful information from the input.

This is not to suggest that machine learning is less useful or less powerful than deep learning - on the contrary simpler algorithms regularly beat deep learning algorithms for certain kinds of tasks. The type of learning to use is driven by the use case, performance obtained, and the desired explainability.

Next, we will cover the key machine learning algorithms that are used for classification, regression and clustering. We will cover deep learning in the next chapter.

Agenda:

1. Decision Trees
2. Random Forest
3. XGBoost
4. Linear Discriminant Analysis
5. Support Vector Machines
6. Naïve Bayes
7. K-Nearest Neighbors
8. K-Means Clustering
9. Hierarchical Clustering

All our work will follow the ML workflow discussed earlier, and repeated below:

1. Prepare your data – cleanse, convert to numbers, etc
2. Split the data into training and test sets
   • Training sets are what algorithms learn from
   • Test sets are the ‘hold-out’ data on which model effectiveness is measured
   • No set rules, often a 80:20 split between train and test data suffices. If there is a lot of training data, you may keep a smaller number as the test set.
3. Fit a model.
4. Check model accuracy based on the test set.
5. Use for predictions.

What you need to think about
As we cover the algorithms, think about the below ideas for each.

1. What is the conceptual basis for the algorithm?
   • This will help you think about the problems the algorithm can be applied to,
   • You should also think about the parameters you can control in the model,
   • You should think about model explainability, how essential is it to your use case, and who your audience is.
2. Do you need to scale/standardize the data?
   • Or can you use the raw data as is?
3. Whether it can perform regression, classification or clustering
   • Regression models help forecast numeric quantities, while classification algorithms help determine class membership.
   • Some algorithms can only perform either regression or classification, while others can do both.
4. If it is a classification algorithm, does it provide just the class membership, or probability estimates
   • If reliable probability estimates are available from the model, you can perform more advanced model evaluations, and tweak the probability cut-off to obtain your desired True Positive/False Positive rates.

Some library imports first...

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
import seaborn as sns
from sklearn import tree
from sklearn.metrics import confusion_matrix, accuracy_score, classification_report
from sklearn.metrics import mean_absolute_error, mean_squared_error, ConfusionMatrixDisplay
from sklearn import metrics
# from sklearn.metrics import mean_absolute_percentage_error

from sklearn.preprocessing import StandardScaler, LabelEncoder
from sklearn.model_selection import train_test_split
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.ensemble import RandomForestClassifier
from sklearn import svm

import sklearn.preprocessing as preproc

Decision Trees

Decision Trees (DTs) are a non-parametric supervised learning method used for classification and regression. The goal is to create a model that predicts the value of a target variable by learning simple decision rules inferred from the data features. In simple words, decision trees are a collection of if/else conditions that are applied to data till a prediction is reached. Trees are constructed by splitting data by a variable, and then doing the same again and again till the desired level of accuracy is reached (or we run out of data). Trees can be visualized, and are therefore easier to interpret – and in that sense they are a ‘white box model’.

Trivial Example
Consider a made-up dataset where we know the employment and housing status of our customers, and whether they have paid back or defaulted on their loans. When a new customer requests a loan, can we use this data to decide the whether there is likely to be a default?

We would like the ‘leaf nodes’ to be pure, ie contain instances that tend to belong to the same class.

Several practical issues arise in the above example: - Attributes rarely neatly split a group. In the made up example, everything lined up neatly but rarely will in reality.
- How does one select what order to select attributes in? We could have started with housing instead of looking at whether a person was an employee or not.
- Many attributes will not be binary, may have multiple unique values.
- Some attributes may be numeric. For example, we may know their credit scores. In such a case, how do we split the nodes?
- Finally, how do we decide we are done? Should we keep going till we run out of variables, or till all leaf nodes are pure?

Measuring Purity - Entropy and Gini Impurity

Entropy
The most common splitting criterion is called information gain, and is based on a measure called entropy. Entropy is a measure of disorder that can be applied to a collection.
Disorder corresponds to how mixed (impure) the group is with respect to the properties of interest.

$\mbox{Entropy} = -p_1 log_2(p_1) -p_2 log_2(p_2) - ...$

A node is pure when entropy = 0. So we are looking for ways to minimize entropy.

Gini Impurity
Another measure of impurity is the Gini Impurity.
$\mbox{Gini Index} = 1 - p_1^2 - p_2^2$

Like entropy, the Gini Impurity has a minimum of 0. In a two class problem, the maximum value for the Gini Impurity will be 0.5. Both Entropy and the Gini Impurity behave similarly, the Gini Impurity is supposedly less computationally intensive.

With entropy as the measure of disorder, we calculate Information Gain offered by each attribute when used as the basis of segmentation.

Information gain is the reduction in entropy by splitting our data on the basis of a single attribute.

For our toy example, the entropy for the top parent node was 0.95. This was reduced to 0.41 at the next child node, calculated as $p(c_1) * entropy(c_1) + p(c_2) * entropy(c_2) + ...$. We start our segmentation with the attribute that provides the most information gain.

Fortunately, automated algorithms do this for us, so we do not have to calculate any of this. But the concept of information gain and how regression tree algorithms decide to split the data is important to be aware of.

Toy Example Continued
Let us continue the example introduced earlier.

# Entropy before the first split
entropy1 = -((6/16) * np.log2(6/16))-((10/16) * np.log2(10/16))
entropy1

0.954434002924965

# Entroy after the split
entropy2 = \
(8/16) * (-(8/8) * np.log2(8/8)) \
+ \
((8/16) * (- (2/8) * np.log2(2/8) - (6/8) * np.log2(6/8) ))
entropy2

0.4056390622295664

#Information Gain
entropy1 - entropy2

0.5487949406953987

Another simple example
We look at another example where we try to build a decision tree to predict whether a debt was written-off for a customer given other attributes.

df = pd.read_excel('write-off.xlsx')

df

Balance, Age and Employed are independent variables, and Write-off is the predicted variable. Of these, the Write-off and Employed columns are strings and have to be converted to numerical variables so they can be used in algorithms.

df['Write-off'] = df['Write-off'].astype('category') #convert to category
df['write-off-label'] = df['Write-off'].cat.codes #use category codes as labels
df = pd.get_dummies(df, columns=["Employed"]) #one hot encoding using pandas
df

type(df['Write-off'])

pandas.core.series.Series

df = df.iloc[:,[0,3,4,1,2,5,6]]
df

df.iloc[:, 2:]

# This below command is required only to get back to the home folder if you aren't there already

# import os
# os.chdir('/home/jovyan')

X = df.iloc[:,3:]
y = df.iloc[:,2]
clf = tree.DecisionTreeClassifier()
clf = clf.fit(X, y)

import graphviz 
dot_data = tree.export_graphviz(clf, out_file=None) 
graph = graphviz.Source(dot_data) 
# graph.render("df") 

dot_data = tree.export_graphviz(clf, out_file=None, 
                         feature_names=X.columns,  
                         class_names=['yes', 'no'],  # Plain English names for classes_
                         filled=True, rounded=True,  
                         special_characters=True)  
graph = graphviz.Source(dot_data)  
graph 

clf.classes_

array([0, 1], dtype=int8)

y

0    1
1    0
2    0
3    1
4    0
Name: write-off-label, dtype: int8

Iris Flower Dataset
We consider the Iris dataset, a multivariate data set introduced by the British statistician and biologist Ronald Fisher in a 1936 paper. The data set consists of 50 samples from each of three species of Iris (Iris setosa, Iris virginica and Iris versicolor).

Four features were measured from each sample: the length and the width of the sepals and petals, in centimeters.

Source: Wikipedia

Image Source/Attribution: https://commons.wikimedia.org/w/index.php?curid=248095

Difference between a petal and a sepal:

Scikit Learn’s decision tree classifier algorithm, combined with another package called graphviz, can provide decision trees together with good graphing capabilities.

Unfortunately, sklearn requires all data to be numeric and as numpy arrays. This creates practical problems for the data analyst – categorical variables have to be labeled or one-hot encoded, and their plain English meanings have to be tracked separately.

# Load the data
iris = sm.datasets.get_rdataset('iris').data
iris.sample(6)

Our task is: Based on these features, can we create a decision tree to distinguish between the three species of the Iris flower?

# Let us look at some basic descriptive stats for each of the flower species.  

iris.pivot_table(columns = ['Species'], aggfunc = [np.mean, min, max]).transpose()

# Next, we build the decision tree  

import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn import tree
iris = load_iris()
X, y = iris.data, iris.target

# Train-test split
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.20)

# Create the classifier and visualize the decision tree
clf = tree.DecisionTreeClassifier()
clf = clf.fit(X_train, y_train)

import graphviz 
dot_data = tree.export_graphviz(clf, out_file=None) 
graph = graphviz.Source(dot_data) 
graph.render("iris") 

dot_data = tree.export_graphviz(clf, out_file=None, 
                         feature_names=iris.feature_names,  
                         class_names=iris.target_names,  
                         filled=True, rounded=True,  
                         special_characters=True)  
graph = graphviz.Source(dot_data)  
graph 

# Another way to build the tree is something as simple as 
# typing `tree.plot_tree(clf)` but the above code gives 
# us much better results.

print(tree.plot_tree(clf))

[Text(0.4444444444444444, 0.9166666666666666, 'x[2] <= 2.5\ngini = 0.666\nsamples = 120\nvalue = [43, 38, 39]'), Text(0.3333333333333333, 0.75, 'gini = 0.0\nsamples = 43\nvalue = [43, 0, 0]'), Text(0.5555555555555556, 0.75, 'x[3] <= 1.75\ngini = 0.5\nsamples = 77\nvalue = [0, 38, 39]'), Text(0.3333333333333333, 0.5833333333333334, 'x[2] <= 5.35\ngini = 0.139\nsamples = 40\nvalue = [0, 37, 3]'), Text(0.2222222222222222, 0.4166666666666667, 'x[2] <= 4.95\ngini = 0.051\nsamples = 38\nvalue = [0, 37, 1]'), Text(0.1111111111111111, 0.25, 'gini = 0.0\nsamples = 35\nvalue = [0, 35, 0]'), Text(0.3333333333333333, 0.25, 'x[1] <= 2.45\ngini = 0.444\nsamples = 3\nvalue = [0, 2, 1]'), Text(0.2222222222222222, 0.08333333333333333, 'gini = 0.0\nsamples = 1\nvalue = [0, 0, 1]'), Text(0.4444444444444444, 0.08333333333333333, 'gini = 0.0\nsamples = 2\nvalue = [0, 2, 0]'), Text(0.4444444444444444, 0.4166666666666667, 'gini = 0.0\nsamples = 2\nvalue = [0, 0, 2]'), Text(0.7777777777777778, 0.5833333333333334, 'x[2] <= 4.85\ngini = 0.053\nsamples = 37\nvalue = [0, 1, 36]'), Text(0.6666666666666666, 0.4166666666666667, 'x[0] <= 5.95\ngini = 0.5\nsamples = 2\nvalue = [0, 1, 1]'), Text(0.5555555555555556, 0.25, 'gini = 0.0\nsamples = 1\nvalue = [0, 1, 0]'), Text(0.7777777777777778, 0.25, 'gini = 0.0\nsamples = 1\nvalue = [0, 0, 1]'), Text(0.8888888888888888, 0.4166666666666667, 'gini = 0.0\nsamples = 35\nvalue = [0, 0, 35]')]

# List categories in the classifier
iris.target_names

array(['setosa', 'versicolor', 'virginica'], dtype='<U10')

# Perform predictions
clf.predict(X_test)

array([1, 1, 1, 2, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 2, 1, 2, 2, 1, 2, 2,
       1, 0, 2, 0, 2, 0, 2, 0])

Confusion Matrix and Classification Report

We did not split the data into train/test sets. For now, we will evaluate the model based on the entire data set (ie, on the training set).

For this trivial example, the decision tree has done a perfect job of predicting flower species.

confusion_matrix(y_true = y_test, y_pred = clf.predict(X_test))

array([[ 7,  0,  0],
       [ 0, 12,  0],
       [ 0,  2,  9]], dtype=int64)

print(classification_report(y_true = y_test, y_pred = clf.predict(X_test)))

              precision    recall  f1-score   support

           0       1.00      1.00      1.00         7
           1       0.86      1.00      0.92        12
           2       1.00      0.82      0.90        11

    accuracy                           0.93        30
   macro avg       0.95      0.94      0.94        30
weighted avg       0.94      0.93      0.93        30

ConfusionMatrixDisplay.from_estimator(clf, X_test, y_test, display_labels=iris.target_names);

Confusion Matrix for All Data

# Just for the heck of it, let us predict the entire dataset using our model, 
# and check the results
print(classification_report(y_true = y, y_pred = clf.predict(X)))

              precision    recall  f1-score   support

           0       1.00      1.00      1.00        50
           1       0.96      1.00      0.98        50
           2       1.00      0.96      0.98        50

    accuracy                           0.99       150
   macro avg       0.99      0.99      0.99       150
weighted avg       0.99      0.99      0.99       150

ConfusionMatrixDisplay.from_estimator(clf, X, y, display_labels=iris.target_names);

Class probabilities with decision trees

Decision Trees do not do a great job of predicting the probability of belonging to a particular class, for example, when compared to Logistic Regression.

Probabilities for class membership are just the proportion of observations in a particular class in the appropriate leaf node. For a tree with unlimited nodes, we will always mostly have p=100% for most predictions.

Scikit Learn provides a method to predict probabilities, clf.predict_proba(). If we apply this to our decision tree (first five observations only), we get as below:

# As can be seen below, the model does not give class probabilities
clf.predict(X)

array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])

# Get class probabilities
clf.predict_proba(X[:5])

array([[1., 0., 0.],
       [1., 0., 0.],
       [1., 0., 0.],
       [1., 0., 0.],
       [1., 0., 0.]])

Predictions with decision trees

With this model, how do I predict if I have the measurements for a new flower?

Once a Decision Tree Classifier is built, new predictions can be obtained using the predict(X) method.

Imagine we have a new flower with dimensions 5, 3, 1 and 2 and need to predict its species. Since we have the featureset, we feed this information to the model and obtain the prediction.

Refer below for the steps in Python

# let us remind ourselves of what features need to predict a flower's species

iris.feature_names

['sepal length (cm)',
 'sepal width (cm)',
 'petal length (cm)',
 'petal width (cm)']

# let us also look at existing feature set

X[:4]

array([[5.1, 3.5, 1.4, 0.2],
       [4.9, 3. , 1.4, 0.2],
       [4.7, 3.2, 1.3, 0.2],
       [4.6, 3.1, 1.5, 0.2]])

# Next, the measurements for the new flower

new_flower = [[5,3,1,2]]

# Now the prediction

clf.predict(new_flower)

array([0])

# The above means it is the category at index 1 in the target
# Let us look at what the target names are.
# We see that the 'versicolor' is at index 1, so that is the prediction for the new flower

iris.target_names

array(['setosa', 'versicolor', 'virginica'], dtype='<U10')

# or, all of the above in one line

print(iris.target_names[clf.predict(new_flower)])

['setosa']

Decision Tree Regression

Decision trees can also be applied to estimating continuous values for the target variable.

They work in the same way as decision trees for classification, except that information gain is measured differently, eg by a reduction in standard deviation at the node level.

So splits for a node would be performed based on a variable/value that creates the maximum reduction in the standard deviation of the y values in the node.

The prediction is then the average of the observations in the leaf node.

As an example, let us consider the Boston House Price dataset that is built into sklearn. There are 506 rows × 14 variables

# Load the data

from sklearn.datasets import fetch_california_housing
housing = fetch_california_housing()

X = housing['data']
y = housing['target']
features = housing['feature_names']
DESCR = housing['DESCR']

cali_df = pd.DataFrame(X, columns = features)
cali_df.insert(0,'medv', y)

cali_df

# Let us look at the data dictionary

print(DESCR)

.. _california_housing_dataset:

California Housing dataset
--------------------------

**Data Set Characteristics:**

    :Number of Instances: 20640

    :Number of Attributes: 8 numeric, predictive attributes and the target

    :Attribute Information:
        - MedInc        median income in block group
        - HouseAge      median house age in block group
        - AveRooms      average number of rooms per household
        - AveBedrms     average number of bedrooms per household
        - Population    block group population
        - AveOccup      average number of household members
        - Latitude      block group latitude
        - Longitude     block group longitude

    :Missing Attribute Values: None

This dataset was obtained from the StatLib repository.
https://www.dcc.fc.up.pt/~ltorgo/Regression/cal_housing.html

The target variable is the median house value for California districts,
expressed in hundreds of thousands of dollars ($100,000).

This dataset was derived from the 1990 U.S. census, using one row per census
block group. A block group is the smallest geographical unit for which the U.S.
Census Bureau publishes sample data (a block group typically has a population
of 600 to 3,000 people).

A household is a group of people residing within a home. Since the average
number of rooms and bedrooms in this dataset are provided per household, these
columns may take surprisingly large values for block groups with few households
and many empty houses, such as vacation resorts.

It can be downloaded/loaded using the
:func:`sklearn.datasets.fetch_california_housing` function.

.. topic:: References

    - Pace, R. Kelley and Ronald Barry, Sparse Spatial Autoregressions,
      Statistics and Probability Letters, 33 (1997) 291-297

We can fit a decision tree regressor to the data.
1. First, we load the data.
1. Next, we split the data into train and test sets, keeping 20% for the test set.
1. Then we fit a model to the training data, and store the model object in the variable model.
1. Next we use the model to predict the test cases.
1. Finally, we evaluate the results.

# Train-test split
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.20)

# model = tree.DecisionTreeRegressor()
model = tree.DecisionTreeRegressor(max_depth=9)
model = model.fit(X_train, y_train)
model.predict(X_test)

array([3.0852    , 0.96175   , 0.96001341, ..., 2.26529224, 2.15065625,
       2.02124038])

print(model.tree_.max_depth)

9

y_pred  =  model.predict(X_test)

print('MSE = ', mean_squared_error(y_test,y_pred))
print('RMSE = ', np.sqrt(mean_squared_error(y_test,y_pred)))
print('MAE = ', mean_absolute_error(y_test,y_pred))

MSE =  0.38773905612449
RMSE =  0.622686964794101
MAE =  0.4204310076086696

# Just checking to see if we have everything working right

print('Count of predictions:', len(y_pred))
print('Count of ground truth labels:', len(y_test))

Count of predictions: 4128
Count of ground truth labels: 4128

# We plot the actual home prices vs the predictions in a scatterplot

plt.figure(figsize = (8,8))
plt.scatter(y_test, y_pred, alpha=0.5)
plt.title('Actual vs Predicted Home Value in $000s \n Closer to red line (identity) \
          means more accurate prediction')
plt.plot( [0,5],[0,5], color='red' )
plt.xlabel("Actual")
plt.ylabel("Predicted")

Text(0, 0.5, 'Predicted')

# Context for the RMSE.  What is the mean, min and max?

cali_df.medv.describe()

count    20640.000000
mean         2.068558
std          1.153956
min          0.149990
25%          1.196000
50%          1.797000
75%          2.647250
max          5.000010
Name: medv, dtype: float64

# R-squared calculation
pd.DataFrame({'actual':y_test, 'predicted':y_pred}).corr()**2

How well did my model generalize?

Let us see how my model did on the training data

# R-squared
pd.DataFrame({'actual':y_train, 'predicted':model.predict(X_train)}).corr()**2

# Calculate MSE, RMSE and MAE

y_pred  =  model.predict(X_train)

print('MSE = ', mean_squared_error(y_train,y_pred))
print('RMSE = ', np.sqrt(mean_squared_error(y_train,y_pred)))
print('MAE = ', mean_absolute_error(y_train,y_pred))

MSE =  0.27139484003226105
RMSE =  0.5209556987232802
MAE =  0.3574534397649355

# Scatterplot for actual vs predicted on TRAINING data

plt.figure(figsize = (8,8))
plt.scatter(y_train, model.predict(X_train), alpha=0.5)
plt.title('Actual vs Predicted Home Value in $000s \n Closer to red line (identity) means more accurate prediction\n TRAINING DATA')
plt.plot( [0,5],[0,5], color='red' )
plt.xlabel("Actual")
plt.ylabel("Predicted")

Text(0, 0.5, 'Predicted')

Addressing Overfitting in Decision Trees

The simplest way to address overfitting in decision trees is to limit the depth of the trees using the max_depth parameter when fitting the model. The depth of a decision tree is the length of the longest path from a root to a leaf.

Find out the current value of the max tree depth in the example (print(model.tree_.max_depth)), and change the max_depth parameter to see if you can reduce the RMSE for the test set.

You can also change the minimum count of samples required to be present in a leaf node (min_samples_leaf), and the minimum number of observations required before a node is allowed to split (min_samples_split).

Random Forest

A random forest fits a number of decision tree classifiers on various sub-samples of the dataset and uses averaging to improve the predictive accuracy and control over-fitting.

Random Forest almost always gives results superior to decision trees, and is therefore preferred over decision trees. However, because the results provided by random forest are the result of averaging multiple trees, explainability can become an issue.

Therefore decision trees may still be preferred over random forest in the interest of explainability.

At this point, it is important to introduce two new concepts: bootstrapping, and bagging.

Bootstrapping
In bootstrapping, you treat the sample as if it were the population, and draw repeated samples of equal size from it. The samples are drawn with replacement. Now think that for each of these new samples you calculate a population characteristic, say the median. Because you potentially have a very large number of samples (theoretically infinite), you can get a distribution of the median of the population from our original single sample.

If we hadn’t done bootstrapping (ie resample from the sample with replacement), we would have only one point estimate for the median.

Bootstrapping improves the estimation process and reduces variance.

Bagging (Bootstrap + Aggregation)
Bagging is a type of ensemble learning. Ensemble learning is where we combine multiple models to produce a better prediction or classification.

In bagging, we produce multiple different training sets (called bootstrap samples), by sampling with replacement from the original dataset. Then, for each bootstrap sample, we build a model.
The results in an ensemble of models, where each model votes with the equal weight. Typically, the goal of this procedure is to reduce the variance of the model of interest (e.g. decision trees).

The Random Forest algorithm is when the above technique is applied to decision trees.

Random Forests
Random forests are an example of ensemble learning, where multiple models are combined to produce a better prediction or classification.

Random forests are collections of trees. Predictions are equivalent to the average prediction of component trees.

Multiple decision trees are created from the source data using a technique called bagging. Multiple different training sets (called bootstrap samples) are created by sampling with replacement from the original dataset.

Then, for each bootstrap sample, we build a model. The results in an ensemble of models, where each model votes with the equal weight. Typically, the goal of this procedure is to reduce the variance of the model of interest.

When applied to decision trees, this becomes random forest.

Random Forest for Classification

# load the data

college = pd.read_csv('collegePlace.csv')
college.shape

(2966, 8)

college

# divide the dataset into train and test sets, separating the features and target variable

X = college[['Age', 'Internships', 'CGPA', 'Hostel', 'HistoryOfBacklogs']].values
y = college['PlacedOrNot'].values

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

# classify using random forest classifier

RandomForest = RandomForestClassifier()
model_rf = RandomForest.fit(X_train, y_train)
pred = model_rf.predict(X_test)

print(classification_report(y_true = y_test, y_pred = pred))
ConfusionMatrixDisplay.from_estimator(model_rf, X_test, y_test);

              precision    recall  f1-score   support

           0       0.79      0.96      0.87       279
           1       0.96      0.77      0.86       315

    accuracy                           0.86       594
   macro avg       0.87      0.87      0.86       594
weighted avg       0.88      0.86      0.86       594

# get probabilities for each observation in the test set

model_rf.predict_proba(X_test)

array([[0.92120988, 0.07879012],
       [0.79331614, 0.20668386],
       [0.        , 1.        ],
       ...,
       [0.01571429, 0.98428571],
       [0.        , 1.        ],
       [0.73565005, 0.26434995]])

y_test

array([0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0,
       0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0,
       0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0,
       0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0,
       1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1,
       1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1,
       1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0,
       0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0,
       0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0,
       0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0,
       1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1,
       0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0,
       1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0,
       0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1,
       0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1,
       1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1,
       0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0,
       1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1,
       1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1,
       1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0,
       0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1,
       1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0,
       0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1,
       1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0,
       0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0,
       1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1,
       1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1],
      dtype=int64)

# get probabilities for each observation in the test set

pred_prob = model_rf.predict_proba(X_test)[:,1]

model_rf.classes_

array([0, 1], dtype=int64)

# Source for code below: https://stackoverflow.com/questions/25009284/how-to-plot-roc-curve-in-python
fpr, tpr, thresholds = metrics.roc_curve(y_test, pred_prob)
roc_auc = metrics.auc(fpr, tpr)
plt.figure(figsize = (9,8))
plt.title('Receiver Operating Characteristic')
plt.plot(fpr, tpr, 'b', label = 'AUC = %0.2f' % roc_auc)
plt.legend(loc = 'lower right')
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
for i, txt in enumerate(thresholds):
    if i in np.arange(1, len(thresholds), 10): # print every 10th point to prevent overplotting:
        plt.annotate(text = round(txt,3), xy = (fpr[i], tpr[i]),  
                     xytext=(-44, 0), textcoords='offset points',
                     arrowprops={'arrowstyle':"simple"}, color='green',fontsize=8)
plt.show()

threshold_dataframe = pd.DataFrame({'fpr':fpr, 'tpr': tpr, 'threshold':thresholds}).sort_values(by='threshold')
threshold_dataframe.head()

Random Forest for Regression

The Random Forest algorithm can also be used effectively for regression problems. Let us try a larger dataset this time.

We will try to predict diamond prices based on all the other attributes we know about the diamonds.

However, our data contains a number of categorical variables. We will need to convert these into numerical using one-hot encoding. Let us do that next!

from sklearn.ensemble import RandomForestRegressor

diamonds = sns.load_dataset("diamonds")

diamonds.head()

diamonds = pd.get_dummies(diamonds)

diamonds.head()

# Define X and y as arrays. y is the price column, X is everything else

X = diamonds.loc[:, diamonds.columns != 'price'].values
y = diamonds.price.values

# Train test split

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.20)

# Fit model

model_rf_regr = RandomForestRegressor(max_depth=2, random_state=0)
model_rf_regr.fit(X_train, y_train)
model_rf_regr.predict(X_test)

array([1054.29089419, 1054.29089419, 1054.29089419, ..., 6145.62603236,
       1054.29089419, 1054.29089419])

# Evaluate model

y_pred  =  model_rf_regr.predict(X_test)
from sklearn.metrics import mean_absolute_error, mean_squared_error
print('MSE = ', mean_squared_error(y_test,y_pred))
print('RMSE = ', np.sqrt(mean_squared_error(y_test,y_pred)))
print('MAE = ', mean_absolute_error(y_test,y_pred))

MSE =  2757832.1354701095
RMSE =  1660.6721938631083
MAE =  1036.4110791707412

# Evaluate residuals

plt.figure(figsize = (8,8))
plt.scatter(y_test, y_pred, alpha=0.5)
plt.title('Actual vs Predicted Diamond Value\n Closer to red line (identity) means more accurate prediction')
plt.plot( [0,19000],[0,19000], color='red' )
plt.xlabel("Actual")
plt.ylabel("Predicted")

Text(0, 0.5, 'Predicted')

diamonds.price.describe()

count    53940.000000
mean      3932.799722
std       3989.439738
min        326.000000
25%        950.000000
50%       2401.000000
75%       5324.250000
max      18823.000000
Name: price, dtype: float64

# R-squared calculation
pd.DataFrame({'actual':y_test, 'predicted':y_pred}).corr()**2

importance = model_rf_regr.feature_importances_
feature_names = diamonds.loc[:, diamonds.columns != 'price'].columns
pd.DataFrame({'Feature':feature_names, 'Importance':importance}).sort_values(by='Importance', ascending=False)

Random Forest Regression - Another Example

Let us look at our California Housing Dataset that we examined before to predict home prices.

# Load the data

from sklearn.datasets import fetch_california_housing
housing = fetch_california_housing()

X = housing['data']
y = housing['target']
features = housing['feature_names']
DESCR = housing['DESCR']

cali_df = pd.DataFrame(X, columns = features)
cali_df.insert(0,'medv', y)

cali_df

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.20)

from sklearn.ensemble import RandomForestRegressor
model = RandomForestRegressor(max_depth=2, random_state=0)
model.fit(X_train, y_train)

RandomForestRegressor(max_depth=2, random_state=0)

RandomForestRegressor(max_depth=2, random_state=0)

y_pred  =  model.predict(X_test)
from sklearn.metrics import mean_absolute_error, mean_squared_error
print('MSE = ', mean_squared_error(y_test,y_pred))
print('RMSE = ', np.sqrt(mean_squared_error(y_test,y_pred)))
print('MAE = ', mean_absolute_error(y_test,y_pred))

MSE =  0.7109737317347218
RMSE =  0.8431925828271509
MAE =  0.6387472402358885

print(cali_df.medv.describe())
cali_df.medv.plot.hist(bins=20)

count    20640.000000
mean         2.068558
std          1.153956
min          0.149990
25%          1.196000
50%          1.797000
75%          2.647250
max          5.000010
Name: medv, dtype: float64





<Axes: ylabel='Frequency'>

plt.figure(figsize = (8,8))
plt.scatter(y_test, y_pred, alpha=0.5)
plt.title('Actual vs Predicted Home Value in $000s \n Closer to red line (identity) means more accurate prediction')
plt.plot( [0,5],[0,5], color='red' )
plt.xlabel("Actual")
plt.ylabel("Predicted")

Text(0, 0.5, 'Predicted')

XGBoost

Like Random Forest, XGBoost is a tree based algorithm. In Random Forest, multiple trees are built in parallel, and averaged. In XGBoost, trees are built sequentially, with each tree correcting the errors of the previous one.

Trees are built in sequence, with each next tree in the sequence targeting the errors of the previous one. The trees are then added, with a multiplicative constant ‘learning rate’ between 0 and 1 applied to each tree.

XGBoost has by far exceeded the performance of other algorithms, and is one of the most used algorithms on Kaggle. In many cases, it outperforms Neural Nets.

Extensive documentation is available at https://xgboost.readthedocs.io/en/latest

Example
Let us consider our college placement dataset, and check if we are able to predict the ‘PlacedOrNot’ variable correctly.

We will convert the categorical variables (stream of study, gender, etc) into numerical using one-hot encoding.

We will keep 20% of the data as the test set, and fit a model using the XGBoost algorithm.

XGBoost - Classification

# load the data
college = pd.read_csv('collegePlace.csv')
college = pd.get_dummies(college)

college

# Test train split
X = college.loc[:, college.columns != 'PlacedOrNot']
y = college['PlacedOrNot']
feature_names = college.loc[:, college.columns != 'PlacedOrNot'].columns

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

# Fit the model
from xgboost import XGBClassifier

model_xgb = XGBClassifier(use_label_encoder=False, objective= 'binary:logistic')
model_xgb.fit(X_train, y_train)

XGBClassifier(base_score=None, booster=None, callbacks=None,
              colsample_bylevel=None, colsample_bynode=None,
              colsample_bytree=None, device=None, early_stopping_rounds=None,
              enable_categorical=False, eval_metric=None, feature_types=None,
              gamma=None, grow_policy=None, importance_type=None,
              interaction_constraints=None, learning_rate=None, max_bin=None,
              max_cat_threshold=None, max_cat_to_onehot=None,
              max_delta_step=None, max_depth=None, max_leaves=None,
              min_child_weight=None, missing=nan, monotone_constraints=None,
              multi_strategy=None, n_estimators=None, n_jobs=None,
              num_parallel_tree=None, random_state=None, ...)

XGBClassifier(base_score=None, booster=None, callbacks=None,
              colsample_bylevel=None, colsample_bynode=None,
              colsample_bytree=None, device=None, early_stopping_rounds=None,
              enable_categorical=False, eval_metric=None, feature_types=None,
              gamma=None, grow_policy=None, importance_type=None,
              interaction_constraints=None, learning_rate=None, max_bin=None,
              max_cat_threshold=None, max_cat_to_onehot=None,
              max_delta_step=None, max_depth=None, max_leaves=None,
              min_child_weight=None, missing=nan, monotone_constraints=None,
              multi_strategy=None, n_estimators=None, n_jobs=None,
              num_parallel_tree=None, random_state=None, ...)

# Perform predictions, and store the results in a variable called 'pred'
pred = model_xgb.predict(X_test)

# Check the classification report and the confusion matrix
print(classification_report(y_true = y_test, y_pred = pred))
ConfusionMatrixDisplay.from_estimator(model_xgb, X_test, y_test);

              precision    recall  f1-score   support

           0       0.82      0.92      0.87       269
           1       0.93      0.84      0.88       325

    accuracy                           0.88       594
   macro avg       0.88      0.88      0.88       594
weighted avg       0.88      0.88      0.88       594

Class Probabilities
We can obtain class probabilities from an XGBoost model. These can help us use different thresholds for cutoff and decide on the error rates we are comfortable with.

model_xgb.classes_

array([0, 1])

y_test

1435    0
1899    1
1475    1
1978    1
100     1
       ..
1614    1
1717    0
556     0
1773    0
1294    0
Name: PlacedOrNot, Length: 594, dtype: int64

pred_prob = model_xgb.predict_proba(X_test).round(3)
pred_prob

array([[0.314, 0.686],
       [0.004, 0.996],
       [0.002, 0.998],
       ...,
       [0.984, 0.016],
       [0.758, 0.242],
       [0.773, 0.227]], dtype=float32)

# Source for code below: https://stackoverflow.com/questions/25009284/how-to-plot-roc-curve-in-python
fpr, tpr, thresholds = metrics.roc_curve(y_test, pred_prob[:,1])
roc_auc = metrics.auc(fpr, tpr)
plt.figure(figsize = (9,8))
plt.title('Receiver Operating Characteristic')
plt.plot(fpr, tpr, 'b', label = 'AUC = %0.2f' % roc_auc)
plt.legend(loc = 'lower right')
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
for i, txt in enumerate(thresholds):
    if i in np.arange(1, len(thresholds), 10): # print every 10th point to prevent overplotting:
        plt.annotate(text = round(txt,3), xy = (fpr[i], tpr[i]),  
                     xytext=(-44, 0), textcoords='offset points',
                     arrowprops={'arrowstyle':"simple"}, color='green',fontsize=8)
plt.show()

threshold_dataframe = pd.DataFrame({'fpr':fpr, 'tpr': tpr, 'threshold':thresholds}).sort_values(by='threshold')
threshold_dataframe.head()

Change results by varying threshold

# Look at how the probabilities look for the first 10 observations
# The first column is class 0, and the second column is class 1

model_xgb.predict_proba(X_test)[:10]

array([[3.143e-01, 6.857e-01],
       [4.300e-03, 9.957e-01],
       [2.000e-03, 9.980e-01],
       [8.000e-04, 9.992e-01],
       [8.798e-01, 1.202e-01],
       [5.881e-01, 4.119e-01],
       [8.866e-01, 1.134e-01],
       [3.000e-04, 9.997e-01],
       [9.956e-01, 4.400e-03],
       [1.463e-01, 8.537e-01]], dtype=float32)

# Let us round the above as to make it a bit easier to read...
# same thing as prior cell, just presentation
np.round(model_xgb.predict_proba(X_test)[:10], 3)

array([[0.314, 0.686],
       [0.004, 0.996],
       [0.002, 0.998],
       [0.001, 0.999],
       [0.88 , 0.12 ],
       [0.588, 0.412],
       [0.887, 0.113],
       [0.   , 1.   ],
       [0.996, 0.004],
       [0.146, 0.854]], dtype=float32)

# Now see what the actual prediction is for the first 10 items
# You can see the model has picked the most probable item
# for identifying which category it should be assigned.
#
# We can vary the threshold to change the predictions.
# We do this next

model_xgb.predict(X_test)[:10]

array([1, 1, 1, 1, 0, 0, 0, 1, 0, 1])

# Set threshold for identifying class 1
threshold = 0.9

# Create predictions.  Note that predictions give us probabilities, not classes!
pred_prob = model_xgb.predict_proba(X_test)

# We drop the probabilities for class 0, and keep just the second column
pred_prob = pred_prob[:,1]

# Convert probabilities to 1s and 0s based on threshold
pred = (pred_prob>threshold).astype(int)

# confusion matrix
cm = confusion_matrix(y_test, pred)
print ("Confusion Matrix : \n", cm)
ConfusionMatrixDisplay(confusion_matrix=cm).plot();

# accuracy score of the model
print('Test accuracy = ', accuracy_score(y_test, pred))
print(classification_report(y_true = y_test, y_pred = pred,))

Confusion Matrix : 
 [[272   1]
 [ 59 262]]
Test accuracy =  0.898989898989899
              precision    recall  f1-score   support

           0       0.82      1.00      0.90       273
           1       1.00      0.82      0.90       321

    accuracy                           0.90       594
   macro avg       0.91      0.91      0.90       594
weighted avg       0.92      0.90      0.90       594

Feature Importance
Using the method feature_importances_, we can get a sense for what the model considers more important than others. However, feature importance identified in this way should be reviewed in the context of domain knowledge. Refer article at https://explained.ai/rf-importance/

# Check feature importance
# This can be misleading though - check out https://explained.ai/rf-importance/
importance = model_xgb.feature_importances_
pd.DataFrame({'Feature':feature_names, 'Importance':importance}).sort_values(by='Importance', ascending=False)

from xgboost import plot_importance

# plot feature importance
plot_importance(model_xgb)

plt.show()

XGBoost for Regression

Let us try to predict diamond prices again, this time using XGBoost. As we can see below, RMSE is half of what we had with Random Forest.

# Load data

diamonds = sns.load_dataset("diamonds")

diamonds.head()

# Get dummy variables

diamonds = pd.get_dummies(diamonds)

diamonds.head()

# Define X and y as arrays. y is the price column, X is everything else

X = diamonds.loc[:, diamonds.columns != 'price'].values
y = diamonds.price.values

# Define X and y as arrays. y is the price column, X is everything else

X = diamonds.loc[:, diamonds.columns != 'price']
y = diamonds.price

# Train test split

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.20)

# Fit model

from xgboost import XGBRegressor

model_xgb_regr = XGBRegressor()
model_xgb_regr.fit(X_train, y_train)
model_xgb_regr.predict(X_test)

array([ 7206.3213,  3110.482 ,  5646.054 , ..., 13976.481 ,  5555.7554,
       11428.439 ], dtype=float32)

# Evaluate model

y_pred  =  model_xgb_regr.predict(X_test)
from sklearn.metrics import mean_absolute_error, mean_squared_error
print('MSE = ', mean_squared_error(y_test,y_pred))
print('RMSE = ', np.sqrt(mean_squared_error(y_test,y_pred)))
print('MAE = ', mean_absolute_error(y_test,y_pred))

MSE =  280824.1563066477
RMSE =  529.9284445155287
MAE =  276.8015830181774

# Evaluate residuals

plt.figure(figsize = (8,8))
plt.scatter(y_test, y_pred, alpha=0.5)
plt.title('Actual vs Predicted Diamond Value\n Closer to red line (identity) means more accurate prediction')
plt.plot( [0,19000],[0,19000], color='red' )
plt.xlabel("Actual")
plt.ylabel("Predicted");

# R-squared calculation
pd.DataFrame({'actual':y_test, 'predicted':y_pred}).corr()**2

from xgboost import plot_importance

# plot feature importance
plot_importance(model_xgb_regr);

As we can see, XGBoost has vastly improved the prediction results. R-squared is 0.98, and the residual plot looks much better than with Random Forest.

diamonds.price.describe()

count    53940.000000
mean      3932.799722
std       3989.439738
min        326.000000
25%        950.000000
50%       2401.000000
75%       5324.250000
max      18823.000000
Name: price, dtype: float64

Linear Methods

Linear methods are different from tree based methods that we looked at earlier. They approach the problem from the perspective of plotting the points and drawing a line (or a plane) that separates the categories.

Let us consider a toy dataset that we create at random. The dataset has two features (Feature_1 and Feature_2), that help us distinguish between two classes - 0 and 1.

The data is graphed in the scatterplot below. The point to note here is that it is pretty easy to distinguish between the two classes by drawing a straight line between the two classes. The question though is which line is the best possible line for classification, given an infinite number of such lines can be drawn?

# Import libraries

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import make_blobs, make_classification


# Generate random data
X, y, centers = make_blobs(n_samples=30, centers=2,
                           n_features=2, random_state=14,
                           return_centers=True,
                           center_box=(0,20), cluster_std = 5)

# Round to one place of decimal

X = np.round_(X,1)
y = np.round_(y,1)

# Create a dataframe with the features and the y variable

df = pd.DataFrame(dict(Feature_1=X[:,0], Feature_2=X[:,1], Label_y=y))
df = round(df,ndigits=2)


# Plot the data

plt.figure(figsize=(9,9))
sns.scatterplot(data = df, x = 'Feature_1', y = 'Feature_2', style = 'Label_y', hue = 'Label_y', 
                alpha = .8, palette="deep",edgecolor = 'None')

# Plot possible lines to discriminate between classes

plt.plot([0,30],[2.5,9], 'k--')
plt.plot([0,30], [0,12], 'k--')
plt.plot([0,20], [-10,20], 'k--');

If we were to create a decision tree, the problem is solved as the decision tree draws two straight line boundaries - first at Feature_2 > 4.55, and the second at Feature_1 > 20. While this works for the current data, we can obviously see that a more robust and simpler solution would be to draw a straight line between the data that is at an angle separating the two classes.

import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn import tree
# iris = load_iris()

clf = tree.DecisionTreeClassifier()
clf = clf.fit(X, y)

import graphviz 
dot_data = tree.export_graphviz(clf, out_file=None) 
graph = graphviz.Source(dot_data) 
# graph.render("iris") 

dot_data = tree.export_graphviz(clf, out_file=None, 
                         feature_names=['Feature_1', 'Feature_2'],  
                         class_names=['Class_0', 'Class_1'],  
                         filled=True, rounded=True,  
                         special_characters=True)  
graph = graphviz.Source(dot_data)  
graph 

# Plot the data

plt.figure(figsize=(9,9))
sns.scatterplot(data = df, x = 'Feature_1', y = 'Feature_2', style = 'Label_y', hue = 'Label_y', 
                alpha = .8, palette="deep",edgecolor = 'None')

# Plot possible lines to discriminate between classes

plt.plot([0,20],[4.55,4.55], color='green')
plt.plot([20,20], [-12,22], color = 'purple');

However, we can see that a single linear split provides better results.

This is an example of a Linear Classifier. The decision boundary is essentially a line represented as the weighted sum of the two axes. This is called a linear discriminant because it discriminates between the two classes using a linear combination of the independent attributes.

A general linear model would look as follows:

$f(x) = w_0 + w_1 x_1 + w_2 x_2+ ...$

For our example, the linear classifier line is defined by the following example:

$[\mbox{Constant Intercept}] + [\mbox{Coefficient 1} * \mbox{Feature2}] + [\mbox{Coefficient 2} * \mbox{Feature2}] = 0$

The coefficients, or weights, are often loosely interpreted as the importance of the features, assuming all feature values have been normalized.

The question is: How do we identify the correct line as many different lines are possible.

There are many methods to determine the line that serves as our linear discriminant. Each method differs in the ‘objective function’ that is optimized to arrive at the solution.

Two of the common methods used are:
- Linear Discriminant Analysis, and - Support Vector Machines

# Fit linear model
from sklearn.svm import SVC
model_svc = SVC(kernel="linear")
model_svc.fit(X, y)

SVC(kernel='linear')

SVC(kernel='linear')

# Plot the data, and line dividing the classification

plt.figure(figsize=(9,9))
sns.scatterplot(data = df, x = 'Feature_1', y = 'Feature_2', style = 'Label_y', hue = 'Label_y', 
                alpha = .8, palette="deep",edgecolor = 'None');



# Plot the equation of the linear discriminant

w = model_svc.coef_[0]
a = -w[0] / w[1]
xx = np.linspace(df.Feature_1.min()-1, df.Feature_1.max()+1)
yy = a * xx - (model_svc.intercept_[0]) / w[1]
plt.plot(xx, yy, 'k-')

# Identify the support vectors, ie the points that decide the decision boundary

plt.scatter(
    model_svc.support_vectors_[:, 0],
    model_svc.support_vectors_[:, 1],
    s=80,
    facecolors="none",
    zorder=10,
    edgecolors="k",
   )

# Plot the margin lines

margin = 1 / np.sqrt(np.sum(model_svc.coef_**2))
yy_down = yy - np.sqrt(1 + a**2) * margin
yy_up = yy + np.sqrt(1 + a**2) * margin

plt.plot(xx, yy_down, "k--")
plt.plot(xx, yy_up, "k--");

# Another way to plot the decision boundary for SVM models
# Source: https://stackoverflow.com/questions/51297423/plot-scikit-learn-sklearn-svm-decision-boundary-surface

from sklearn.svm import SVC
import matplotlib.pyplot as plt
from mlxtend.plotting import plot_decision_regions

svm_graph = SVC(kernel='linear')
svm_graph.fit(X, y)
plot_decision_regions(X, y, clf=svm_graph, legend=2)
plt.show()

Linear Discriminant Analysis

LDA assumes a normal distribution for the data points for the different categories, and attempts to create a 1D projection in a way that separates classes well.

Fortunately, there are libraries available that do all the tough math for us.

LDA expects predictor variables to be continuous due to its distributional assumption of independent variables being multivariate normal. This limits its use in situations where the predictor variables are categorical.

You do not need to standardize the feature set prior to using linear discriminant analysis. You should rule out logistic regression as a better alternative before using linear discriminant analysis.

LDA in Action

We revisit the collegePlace.csv data.

About the data:
A University Announced Its On-Campus Placement Records For The Engineering Course. The Data Is From The Years 2013 And 2014.

Data Fields:
- Age: Age At The Time Of Final Year - Gender: Gender Of Candidate - Stream: Engineering Stream That The Candidate Belongs To - Internships: Number Of Internships Undertaken During The Course Of Studies, Not Necessarily Related To College Studies Or Stream - CGPA: CGPA Till 6th Semester - Hostel: Whether Student Lives In College Accomodation - HistoryOfBacklogs: Whether Student Ever Had Any Backlogs In Any Subjects - PlacedOrNot: Target Variable

# load the data

college = pd.read_csv('collegePlace.csv')

college

college.columns

Index(['Age', 'Gender', 'Stream', 'Internships', 'CGPA', 'Hostel', 'HistoryOfBacklogs', 'PlacedOrNot'], dtype='object')

# divide the dataset into train and test sets, separating the features and target variable
X = college[['Age', 'Internships', 'CGPA', 'Hostel', 'HistoryOfBacklogs']].values
y = college['PlacedOrNot'].values

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

# apply Linear Discriminant Analysis
LDA = LinearDiscriminantAnalysis()
model_lda = LDA.fit(X = X_train, y = y_train)
pred = model_lda.predict(X_test)

college.PlacedOrNot.value_counts()

PlacedOrNot
1    1639
0    1327
Name: count, dtype: int64

1639/(1639+1327)

0.552596089008766

# evaluate performance

ConfusionMatrixDisplay.from_estimator(model_lda, X_test, y_test);
print(classification_report(y_true = y_test, y_pred = pred))

              precision    recall  f1-score   support

           0       0.74      0.73      0.73       275
           1       0.77      0.78      0.77       319

    accuracy                           0.76       594
   macro avg       0.75      0.75      0.75       594
weighted avg       0.76      0.76      0.76       594

confusion_matrix(y_true = y_test, y_pred = pred)

array([[200,  75],
       [ 70, 249]], dtype=int64)

# Get predictions

model_lda.predict(X_test)

array([0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1,
       1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1,
       0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1,
       1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1,
       0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0,
       0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0,
       0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1,
       1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0,
       1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1,
       1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0,
       1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1,
       0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1,
       1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1,
       0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1,
       0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1,
       1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1,
       0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0,
       0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0,
       1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0,
       1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1,
       0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1,
       1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
       1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0,
       0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1,
       0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0,
       1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1,
       0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1],
      dtype=int64)

# Get probability of class membership 

pred_prob = model_lda.predict_proba(X_test)
pred_prob

array([[0.876743  , 0.123257  ],
       [0.0832324 , 0.9167676 ],
       [0.23516243, 0.76483757],
       ...,
       [0.05691089, 0.94308911],
       [0.11393595, 0.88606405],
       [0.06217806, 0.93782194]])

y_test

array([0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1,
       1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0,
       0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1,
       1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1,
       1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0,
       0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1,
       1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1,
       1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0,
       1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1,
       1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1,
       1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1,
       0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1,
       1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1,
       1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1,
       1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1,
       0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1,
       1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1,
       0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0,
       1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1,
       0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1,
       0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1,
       0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1,
       0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0,
       1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1,
       0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1,
       0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1,
       1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1],
      dtype=int64)

# Source for code below: https://stackoverflow.com/questions/25009284/how-to-plot-roc-curve-in-python
fpr, tpr, thresholds = metrics.roc_curve(y_test, pred_prob[:,1])
roc_auc = metrics.auc(fpr, tpr)
plt.figure(figsize = (9,8))
plt.title('Receiver Operating Characteristic')
plt.plot(fpr, tpr, 'b', label = 'AUC = %0.2f' % roc_auc)
plt.legend(loc = 'lower right')
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
for i, txt in enumerate(thresholds):
    if i in np.arange(1, len(thresholds), 10): # print every 10th point to prevent overplotting:
        plt.annotate(text = round(txt,3), xy = (fpr[i], tpr[i]),  
                     xytext=(-44, 0), textcoords='offset points',
                     arrowprops={'arrowstyle':"simple"}, color='green',fontsize=8)
plt.show()

threshold_dataframe = pd.DataFrame({'fpr':fpr, 'tpr': tpr, 'threshold':thresholds}).sort_values(by='threshold')
threshold_dataframe.head()

Closing remarks on LDA: - LDA can not be applied to regression problems, it is useful only for classification. - LDA does provide class membership probabilities, using the predict_proba() method. - There are additional variations to LDA, eg Quadratic Discriminant Analysis, and those may yield better results by allowing a non-linear decision boundary.

Support Vector Machines

Classification with SVM

SVMs use linear classification techniques, ie, they classify instances based on a linear function of the features. The idea behind SVMs is simple: instead of thinking about separating with a line, fit the fattest possible bar between the two classes.

The objective function for SVM incorporates the idea that a wider bar is better.

Once the widest bar is found, the linear discriminant will be the center line through the bar.

The distance between the dashed parallel lines is called the margin around the linear discriminant, and the objective function attempts to maximize the margin.

SVMs require data to be standardized for best results

SVM Example

We will use the same data as before – collegePlace.csv. However this time we will include all the variables, including the categorical variables. We convert the categorical variables to numerical using dummy variables with pd.get_dummies().

# load the data & convert categoricals into numerical variables
college = pd.read_csv('collegePlace.csv')

college = pd.get_dummies(college)
college

At this point, fitting a simple SVM SVC (Support Vector Classification) model is trivial. Refer code below.

SVM has several variations, including LinearSVC, SVC with Polynomial, etc, refer documentation at https://scikit-learn.org/stable/modules/svm.html.

Note that we have chosen to pre-process and standardize the input data first.

# divide the dataset into train and test sets, separating the features and target variable
X = college.drop(['PlacedOrNot'], axis=1).values
y = college['PlacedOrNot'].values

scale = preproc.StandardScaler().fit(X)
X = scale.transform(X)
# X = preproc.StandardScaler().fit_transform(X) 

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

# fit the model
clf = SVC(probability=True) # setting probability=True here can allow us to get probabilities later
model_svm = clf.fit(X_train, y_train)
pred = model_svm.predict(X_test)

# evaluate performance
ConfusionMatrixDisplay.from_estimator(model_svm, X_test, y_test);
print(classification_report(y_true = y_test, y_pred = pred))

              precision    recall  f1-score   support

           0       0.83      0.93      0.88       254
           1       0.94      0.86      0.90       340

    accuracy                           0.89       594
   macro avg       0.89      0.89      0.89       594
weighted avg       0.89      0.89      0.89       594

pred_prob = model_svm.predict_proba(X_test)
pred_prob

array([[3.28070341e-06, 9.99996719e-01],
       [6.65439974e-01, 3.34560026e-01],
       [8.46441402e-01, 1.53558598e-01],
       ...,
       [1.44786221e-02, 9.85521378e-01],
       [3.99281913e-01, 6.00718087e-01],
       [6.46269728e-01, 3.53730272e-01]])

# Source for code below: https://stackoverflow.com/questions/25009284/how-to-plot-roc-curve-in-python
fpr, tpr, thresholds = metrics.roc_curve(y_test, pred_prob[:,1])
roc_auc = metrics.auc(fpr, tpr)
plt.figure(figsize = (9,8))
plt.title('Receiver Operating Characteristic')
plt.plot(fpr, tpr, 'b', label = 'AUC = %0.2f' % roc_auc)
plt.legend(loc = 'lower right')
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
for i, txt in enumerate(thresholds):
    if i in np.arange(1, len(thresholds), 10): # print every 10th point to prevent overplotting:
        plt.annotate(text = round(txt,3), xy = (fpr[i], tpr[i]),  
                     xytext=(-44, 0), textcoords='offset points',
                     arrowprops={'arrowstyle':"simple"}, color='green',fontsize=8)
plt.show()

threshold_dataframe = pd.DataFrame({'fpr':fpr, 'tpr': tpr, 'threshold':thresholds}).sort_values(by='threshold')
threshold_dataframe.head()

SVMs can predict class probabilities, if probability calculations have been set to True as part of the model fitting process. However, these are not calculated by default by the sklearn algorithm.

SVMs can also be used for regression problems, using the model type SVR (‘R’ standing for regression), which we examine next.

Regression with SVM

We perform regression using SVR from sklearn.

# Load the data

from sklearn.datasets import fetch_california_housing
housing = fetch_california_housing()

X = housing['data']
y = housing['target']
features = housing['feature_names']
DESCR = housing['DESCR']

cali_df = pd.DataFrame(X, columns = features)
cali_df.insert(0,'medv', y)

cali_df

# Train test split

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.20)

# Fit model

from sklearn.svm import SVR
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
model_svr = make_pipeline(StandardScaler(), SVR())
model_svr.fit(X, y)
model_svr = model_svr.fit(X_train, y_train)
model_svr.predict(X_test)

array([2.89622439, 1.92606932, 1.55771122, ..., 1.60987874, 0.82130714,
       2.96297243])

# Evaluate model

y_pred  =  model_svr.predict(X_test)
from sklearn.metrics import mean_absolute_error, mean_squared_error
print('MSE = ', mean_squared_error(y_test,y_pred))
print('RMSE = ', np.sqrt(mean_squared_error(y_test,y_pred)))
print('MAE = ', mean_absolute_error(y_test,y_pred))

MSE =  0.3382225528180732
RMSE =  0.5815690438959704
MAE =  0.39084152978427034

# Look at residuals

plt.figure(figsize = (8,8))
plt.scatter(y_test, y_pred, alpha=0.5)
plt.title('Actual vs Predicted Home Value in $000s \n Closer to red \
line (identity) means more accurate prediction')
plt.plot( [y_test.min(),y_test.max()],[y_test.min(),y_test.max()], color='red' )
plt.xlabel("Actual")
plt.ylabel("Predicted");

print(cali_df.medv.describe())
cali_df.medv.plot.hist(bins=20)

count    20640.000000
mean         2.068558
std          1.153956
min          0.149990
25%          1.196000
50%          1.797000
75%          2.647250
max          5.000010
Name: medv, dtype: float64





<Axes: ylabel='Frequency'>

Naive Bayes

Essentially, the logic behind Naive Bayes is as follows: Instead of taking the absolute probability of something happening, we look at the probability of something happening given other things we know have already happened.

So the probability of a flood in the next 1 week may be say 0.1%, but this probability would be different if we already know that 6 inches of rain has already fallen in the past 24 hours.

For each of the categories to be predicted, Naive Bayes considers the conditional probability given the values of other independent variables.

Naïve Bayes uses categorical predictors. For continuous predictors, it assumes a distribution with a mean and standard deviation, which are used to calculate probabilities used in the algorithm. We do not need to standardize the feature set before using Naïve Bayes.

from sklearn import datasets

X = datasets.load_wine()['data']
y = datasets.load_wine()['target']
features = datasets.load_wine()['feature_names']
DESCR = datasets.load_wine()['DESCR']
classes = datasets.load_wine()['target_names']

wine_df = pd.DataFrame(X, columns = features)
wine_df.insert(0,'class', y)

wine_df

# Let us look at the distribution of the observations across classes
wine_df['class'].value_counts()

class
1    71
0    59
2    48
Name: count, dtype: int64

# Train test split
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.20, random_state=123)

# Fit a NB model, and predict
from sklearn.naive_bayes import GaussianNB
model_nb = GaussianNB()
model_nb.fit(X_train, y_train)
pred = model_nb.predict(X_test)

# Evaluate the model accuracy
y_pred  =  model_nb.predict(X_test)
from sklearn.metrics import confusion_matrix,accuracy_score
print('Accuracy = ', accuracy_score(y_test,y_pred))

Accuracy =  1.0

# Model evaluation using the classification report and the confusion matrix

ConfusionMatrixDisplay.from_estimator(model_nb, X_test, y_test);
print(classification_report(y_true = y_test, y_pred = pred))

              precision    recall  f1-score   support

           0       1.00      1.00      1.00         8
           1       1.00      1.00      1.00        11
           2       1.00      1.00      1.00        17

    accuracy                           1.00        36
   macro avg       1.00      1.00      1.00        36
weighted avg       1.00      1.00      1.00        36

wine_df['class'].value_counts()

class
1    71
0    59
2    48
Name: count, dtype: int64

model_nb.classes_

array([0, 1, 2])

model_nb.predict_proba(X_test)

array([[5.93957711e-24, 6.55944182e-07, 9.99999344e-01],
       [1.17999186e-18, 9.99998706e-01, 1.29386818e-06],
       [1.28963191e-29, 2.29833079e-06, 9.99997702e-01],
       [4.24204725e-14, 9.99995948e-01, 4.05248392e-06],
       [1.11923629e-15, 9.99999937e-01, 6.32987582e-08],
       [4.45783147e-18, 3.45360852e-18, 1.00000000e+00],
       [9.96369215e-01, 3.63078524e-03, 1.41680626e-18],
       [8.63917919e-31, 1.78055017e-06, 9.99998219e-01],
       [1.32009669e-20, 1.51909745e-18, 1.00000000e+00],
       [5.37112173e-07, 9.99999463e-01, 6.89392222e-24],
       [2.30365630e-25, 1.06216462e-08, 9.99999989e-01],
       [1.91497199e-16, 4.50596573e-04, 9.99549403e-01],
       [5.09908113e-24, 4.58569400e-16, 1.00000000e+00],
       [9.99999633e-01, 3.66916909e-07, 1.18786467e-29],
       [1.00000000e+00, 1.68595020e-15, 7.56482058e-39],
       [9.82033272e-15, 2.78899850e-05, 9.99972110e-01],
       [3.32691011e-13, 1.00000000e+00, 1.16501960e-15],
       [1.17131670e-10, 1.00000000e+00, 4.47996641e-15],
       [9.99999998e-01, 2.24766861e-09, 3.94855627e-43],
       [2.20594330e-14, 1.00000000e+00, 1.33543221e-21],
       [3.51550512e-16, 4.13954961e-02, 9.58604504e-01],
       [3.55974671e-27, 1.53210419e-11, 1.00000000e+00],
       [4.61982286e-23, 3.19318406e-17, 1.00000000e+00],
       [3.10631399e-21, 6.97558616e-05, 9.99930244e-01],
       [1.01889857e-05, 9.99989811e-01, 1.85439383e-22],
       [3.97166946e-18, 1.95609517e-11, 1.00000000e+00],
       [2.73836673e-25, 1.21551153e-08, 9.99999988e-01],
       [9.82596161e-05, 9.99901740e-01, 5.51859667e-22],
       [1.00000000e+00, 5.58263839e-13, 2.48529450e-33],
       [1.00000000e+00, 2.89393898e-15, 1.82467084e-47],
       [9.93380089e-01, 6.61991141e-03, 5.81320602e-27],
       [9.99179926e-01, 8.20073577e-04, 1.33912510e-17],
       [2.91531470e-22, 1.00486802e-03, 9.98995132e-01],
       [1.25246874e-09, 9.99999999e-01, 7.80563851e-24],
       [2.40587871e-24, 7.15436698e-17, 1.00000000e+00],
       [2.27484144e-06, 9.99997725e-01, 5.43313331e-25]])

k-Nearest Neighbors

Neighbors-based classification is a type of instance-based learning or non-generalizing learning: it does not attempt to construct a general internal model, but simply stores instances of the training data. Classification is computed from a simple majority vote of the nearest neighbors of each point: a query point is assigned the data class which has the most representatives within the nearest neighbors of the point.

The basic nearest neighbors classification uses uniform weights: that is, the value assigned to a query point is computed from a simple majority vote of the nearest neighbors.

Under some circumstances, it is better to weight the neighbors such that nearer neighbors contribute more to the fit. This can be accomplished through the weights parameter (weights = 'uniform’, or weights = 'distance’).

kNN is an easy to use, intuitive algorithm.

kNN requires variables to be normalized or scaled, else distance calculations can be skewed by numerically large features.

kNN can be used to predict categories as well as continuous variables.

kNN classifier

Example
Let us consider our college placement dataset again. We load the data, and perform a train-test split. We also standard-scale the data $((x - mean)/stdev)$.

# load the data
college = pd.read_csv('collegePlace.csv')
college = pd.get_dummies(college)

college

# Test train split

X = college.loc[:, college.columns != 'PlacedOrNot'].values
X = preproc.StandardScaler().fit_transform(X) 
y = college['PlacedOrNot'].values
feature_names = college.loc[:, college.columns != 'PlacedOrNot'].columns

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=123)

pd.DataFrame(X, columns = college.loc[:, college.columns != 'PlacedOrNot'].columns).describe()

# Fit the model using 5 neighbors

from sklearn.neighbors import KNeighborsClassifier
model_knn = KNeighborsClassifier(n_neighbors=5)
model_knn.fit(X_train, y_train)

KNeighborsClassifier()

KNeighborsClassifier()

# Perform predictions, and store the results in a variable called 'pred'

pred = model_knn.predict(X_test)

# Check the classification report and the confusion matrix

print(classification_report(y_true = y_test, y_pred = pred))
ConfusionMatrixDisplay.from_estimator(model_knn, X_test, y_test);

              precision    recall  f1-score   support

           0       0.82      0.85      0.83       265
           1       0.87      0.85      0.86       329

    accuracy                           0.85       594
   macro avg       0.85      0.85      0.85       594
weighted avg       0.85      0.85      0.85       594

pred_prob = model_knn.predict_proba(X_test).round(3)
pred_prob

array([[0.6, 0.4],
       [0. , 1. ],
       [0. , 1. ],
       ...,
       [0. , 1. ],
       [1. , 0. ],
       [0. , 1. ]])

# Source for code below: https://stackoverflow.com/questions/25009284/how-to-plot-roc-curve-in-python
fpr, tpr, thresholds = metrics.roc_curve(y_test, pred_prob[:,1])
roc_auc = metrics.auc(fpr, tpr)
plt.figure(figsize = (9,8))
plt.title('Receiver Operating Characteristic')
plt.plot(fpr, tpr, 'b', label = 'AUC = %0.2f' % roc_auc)
plt.legend(loc = 'lower right')
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
for i, txt in enumerate(thresholds):
    if i in np.arange(1, len(thresholds), 1): # print every n-th point to prevent overplotting:
        plt.annotate(text = round(txt,3), xy = (fpr[i], tpr[i]),  
                     xytext=(-44, 0), textcoords='offset points',
                     arrowprops={'arrowstyle':"simple"}, color='green',fontsize=8)
plt.show()

threshold_dataframe = pd.DataFrame({'fpr':fpr, 'tpr': tpr, 'threshold':thresholds}).sort_values(by='threshold')
threshold_dataframe.head()

What is the right number of neighbors to use?
We used 5 as the count of neighbors in the example above. But how do we know that 5 is the correct number of neighbors? The fact is, we don't know.

So we can try several counts of number of neighbors to see what gives us the best result for the metric we are interested in. Let us consider accuracy as the measure we are looking to improve. We will build the model several times, each time with a different count of the number of neighbors, and calculate the accuracy each time.

# Loop through n from 1 to 25 to find the best n  

acc = []
for n in range(1,25):
    model_knn = KNeighborsClassifier(n_neighbors=n)
    model_knn.fit(X_train, y_train)
    pred = model_knn.predict(X_test)
    acc.append([n, accuracy_score(y_test, pred)])
sns.lineplot(data = pd.DataFrame(acc, columns=['n','accuracy']), x = 'n', y = 'accuracy')
plt.show()
pd.DataFrame(acc, columns=['n','accuracy']).sort_values(by='accuracy', ascending=False).head()

kNN Regressor

# Load data
diamonds = sns.load_dataset("diamonds")

diamonds.head()

# Get dummy variables
diamonds = pd.get_dummies(diamonds)

diamonds.head()

# Define X and y as arrays. y is the price column, X is everything else
X = diamonds.loc[:, diamonds.columns != 'price'].values
X = preproc.StandardScaler().fit_transform(X) 
y = diamonds.price.values

X.shape

(53940, 26)

# Train test split
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.20)

# Fit model
from sklearn.neighbors import KNeighborsRegressor
model_knn_regress = KNeighborsRegressor(n_neighbors=1)
model_knn_regress.fit(X_train, y_train)

KNeighborsRegressor(n_neighbors=1)

KNeighborsRegressor(n_neighbors=1)

# Evaluate model
y_pred  =  model_knn_regress.predict(X_test)
from sklearn.metrics import mean_absolute_error, mean_squared_error
print('MSE = ', mean_squared_error(y_test,y_pred))
print('RMSE = ', np.sqrt(mean_squared_error(y_test,y_pred)))
print('MAE = ', mean_absolute_error(y_test,y_pred))

MSE =  819278.694475343
RMSE =  905.1401518413284
MAE =  448.56952169076754

# Evaluate residuals
plt.figure(figsize = (8,8))
plt.scatter(y_test, y_pred, alpha=0.5)
plt.title('Actual vs Predicted Diamond Value\n Closer to red line (identity) means more accurate prediction')
plt.plot( [0,19000],[0,19000], color='red' )
plt.xlabel("Actual")
plt.ylabel("Predicted");

# R-squared calculation
pd.DataFrame({'actual':y_test, 'predicted':y_pred}).corr()**2

diamonds.price.describe()

count    53940.000000
mean      3932.799722
std       3989.439738
min        326.000000
25%        950.000000
50%       2401.000000
75%       5324.250000
max      18823.000000
Name: price, dtype: float64

k-Means Clustering

k-means clustering is a non-hierarchical approach to clustering. The goal is to divide the observations into a number of non-overlapping groups, or clusters, in a way that the clusters are as homogenous as possible.

k stands for the number of clusters the observations are divided into. There is no natural number of clusters, it is a user provided parameter.

Homogeneity within the cluster is measured using some measure of dispersion, for example, the sum of Euclidean distances.

The algorithm is iterative, and roughly works as follows:
1. Select any k data points as cluster centers (the centroid).
2. Assign all observations to the cluster centroid closest to the observations.
3. Recompute the location of the centroids once all data points have been assigned.
4. Repeat steps 2 and 3.
5. Stop when the measure of dispersion stops improving, or a certain number of repetitions have been performed.

Limitations of k-Means clustering
- k is chosen manually and the correct value for k may be difficult to know. (Algorithms, such as the ‘elbow method’, are available to identify an appropriate value of k.)
- Clusters have no intuitive meaning.
- You may get different results each time due to dependence on initial values. You can overcome this by running k-means several times and picking the best result.
- As the count of dimensions increases (say 1000), PCA may need to be used to prevent similarity measures converging to a constant value.
- Outliers may impact k-means disproportionately.

In spite of the above, k-means clustering remains a preferred clustering technique given its simplicity, scalability to large data sets, and adaptability to different kinds of data.

k-means - example

Let us use the Iris dataset to create a classification model. We will try to cluster the Iris data into 3 clusters using the data in the first four columns. We would like to see the clusters correspond to the species as members of the same species have similar features.

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.cluster import KMeans

iris = sm.datasets.get_rdataset('iris').data

iris

x = iris.iloc[:, 0:4]
kmeans = KMeans(3, n_init='auto')
clusters = kmeans.fit_predict(x)

iris['clusters'] = clusters

df = iris.loc[:,['Species', 'clusters']]
pd.crosstab(index = df['Species'], 
            columns = df['clusters'], margins=True)

plt.figure(figsize = (8,6))
sns.scatterplot(x='Sepal.Width', y='Petal.Width', data=iris, 
                hue='Species', style='clusters', 
                markers= {0: "s", 1: "X", 2: "P"});

We run the KMeans algorithm (from scikit learn) on the data and obtain 3 clusters.

How do the clusters look?

• All setosa are neatly included in cluster 1.
• Versicolor are mostly in cluster 0, but 2 are in a different cluster.
• Verginica is spread across two clusters.

kmeans.fit_predict(x)

array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 1, 2, 2, 2,
       2, 2, 2, 1, 1, 2, 2, 2, 2, 1, 2, 1, 2, 1, 2, 2, 1, 1, 2, 2, 2, 2,
       2, 1, 2, 2, 2, 2, 1, 2, 2, 2, 1, 2, 2, 2, 1, 2, 2, 1])

len(kmeans.fit_transform(x))

150

kmeans.fit_transform(x)[:3]

array([[0.14135063, 5.03132789, 3.41251117],
       [0.44763825, 5.08750645, 3.38963991],
       [0.4171091 , 5.25229169, 3.56011415]])

x[:3]

Right number of clusters

kmeans score is a measure of how far the data points are from the cluster centroids, expressed as a negative number. The closer it is to zero, the better it is. Of course, if we have the number of clusters equal to the number of observations, the score will be zero as each point will be its own centroid, with a sum of zero. If we have only one cluster, we will have a large negative score.

The ideal number of clusters is somewhere when we start getting diminished returns to adding more clusters. We can run the kmeans algorithm for a range of cluster numbers, and compare the score.

KMeans works by minimizing the sum of squared distance of each observation to their respective cluster center. In an extreme situation, all observations would coincide with their centroid center, and the sum of squared distances will be zero.

With sklearn, we can get sum of squared distances of samples to their closest cluster center using _model_name.intertia__.

The negative of inertia_ is model_name.score(x), where x is the dataset kmeans was fitted on.

Elbow Method

The elbow method tracks the sum of squares against the number of clusters, and we can make a subjective judgement on the appropriate number of clusters based on graphing the sum of squares as below. The sum of squares is calculated using the distance between cluster centers and each observation in that cluster. As an extreme case, when the number of clusters is equal to the number of observations, the sum of squares will be zero.

num_clusters = []
score = []
for cluster_count in range(1,15):
    kmeans = KMeans(cluster_count, n_init='auto')
    kmeans.fit(x)
    kmeans.score(x)
    num_clusters.append(cluster_count)
    # score.append(kmeans.score(x)) # score is just the negative of inertia_
    score.append(kmeans.inertia_)

plt.plot(num_clusters, score)

[<matplotlib.lines.Line2D at 0x1e7a09c3ad0>]

print(kmeans.score(x))
kmeans.inertia_

-21.55109126984127





21.55109126984127

# Alternative way of listing labels for the training data
kmeans.labels_

array([ 1, 11, 11, 11,  1,  7, 11,  1, 11, 11,  1,  1, 11, 11,  7,  7,  7,
        1,  7,  1,  1,  1, 11,  1,  1, 11,  1,  1,  1, 11, 11,  1,  7,  7,
       11, 11,  1,  1, 11,  1,  1, 10, 11,  1,  1, 11,  1, 11,  1,  1,  5,
        5,  5, 12,  5,  6,  5,  9,  5, 12,  9,  6,  0,  5, 12,  5,  6, 12,
        0, 12,  8, 12,  3,  5,  5,  5,  5,  5,  5, 12, 12, 12, 12,  3,  6,
        5,  5,  0,  6, 12,  6,  5, 12,  9,  6,  6,  6,  5,  9,  6, 13,  8,
        2, 13, 13,  4,  6,  4, 13,  2, 13, 13,  2,  8,  8, 13, 13,  4,  4,
        3,  2,  8,  4,  3,  2,  2,  3,  8, 13,  2,  4,  4, 13,  3,  3,  4,
       13, 13,  8,  2,  2,  2,  8,  2,  2, 13,  3, 13, 13,  8])

Silhouette Plot

The silhouette plot is a measure of how close each point in one cluster is to points in the neighboring clusters. It provides a visual way to assess parameters such as the number of clusters visually. It does so using the silhouette coefficient.

Silhouette coefficient - This measure has a range of [-1, 1]. Higher the score the better, so +1 is the best result.

The silhouette coefficient is calculated individually for every observation in a cluster as follows: (b - a) / max(a, b). 'b' is the distance between a sample and the nearest cluster that the sample is not a part of. 'a' is the distance between the sample and the cluster it is a part of. One would expect b - a to be a positive number, but if it is not, then likely the point is misclassified.

sklearn.metrics.silhouette_samples(X) - gives the silhouette coefficient for every point in X.
sklearn.metrics.silhouette_score(X) - gives mean of the above.

The silhouette plot gives the mean (ie silhouette_score) as a red vertical line for the entire dataset for all clusters. Then each cluster is presented as a sideways histogram of the distances of each of the datapoints. The fatter the representation of a cluster, the more datapoints are included in that cluster.

Negative points on the histogram indicate misclassifications that may be difficult to correct as moving them changes the centroid center.

# Source: https://scikit-learn.org/stable/auto_examples/cluster/plot_kmeans_silhouette_analysis.html

from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_samples, silhouette_score

import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np

# Generating the sample data from make_blobs
# This particular setting has one distinct cluster and 3 clusters placed close
# together.
# X, y = make_blobs(
#     n_samples=500,
#     n_features=2,
#     centers=4,
#     cluster_std=1,
#     center_box=(-10.0, 10.0),
#     shuffle=True,
#     random_state=1,
# )  # For reproducibility

range_n_clusters = [2, 3, 4, 5, 6]

for n_clusters in range_n_clusters:
    # Create a subplot with 1 row and 2 columns
    fig, (ax1, ax2) = plt.subplots(1, 2)
    fig.set_size_inches(18, 7)

    # The 1st subplot is the silhouette plot
    # The silhouette coefficient can range from -1, 1 but in this example all
    # lie within [-0.1, 1]
    ax1.set_xlim([-.1, 1])
    # The (n_clusters+1)*10 is for inserting blank space between silhouette
    # plots of individual clusters, to demarcate them clearly.
    ax1.set_ylim([0, len(x) + (n_clusters + 1) * 10])

    # Initialize the clusterer with n_clusters value and a random generator
    # seed of 10 for reproducibility.
    clusterer = KMeans(n_clusters=n_clusters, n_init="auto", random_state=10)
    cluster_labels = clusterer.fit_predict(x)

    # The silhouette_score gives the average value for all the samples.
    # This gives a perspective into the density and separation of the formed
    # clusters
    silhouette_avg = silhouette_score(x, cluster_labels)
    print(
        "For n_clusters =",
        n_clusters,
        "The average silhouette_score is :",
        silhouette_avg,
    )

    # Compute the silhouette scores for each sample
    sample_silhouette_values = silhouette_samples(x, cluster_labels)

    y_lower = 2
    for i in range(n_clusters):
        # Aggregate the silhouette scores for samples belonging to
        # cluster i, and sort them
        ith_cluster_silhouette_values = sample_silhouette_values[cluster_labels == i]

        ith_cluster_silhouette_values.sort()

        size_cluster_i = ith_cluster_silhouette_values.shape[0]
        y_upper = y_lower + size_cluster_i

        color = cm.nipy_spectral(float(i) / n_clusters)
        ax1.fill_betweenx(
            np.arange(y_lower, y_upper),
            0,
            ith_cluster_silhouette_values,
            facecolor=color,
            edgecolor=color,
            alpha=0.7,
        )

        # Label the silhouette plots with their cluster numbers at the middle
        ax1.text(-0.05, y_lower + 0.5 * size_cluster_i, str(i))

        # Compute the new y_lower for next plot
        y_lower = y_upper + 10  # 10 for the 0 samples

    ax1.set_title("The silhouette plot for the various clusters.")
    ax1.set_xlabel("The silhouette coefficient values")
    ax1.set_ylabel("Cluster label")

    # The vertical line for average silhouette score of all the values
    ax1.axvline(x=silhouette_avg, color="red", linestyle="--")

    ax1.set_yticks([])  # Clear the yaxis labels / ticks
    ax1.set_xticks([ 0, 0.2, 0.4, 0.6, 0.8, 1])

    # 2nd Plot showing the actual clusters formed
    colors = cm.nipy_spectral(cluster_labels.astype(float) / n_clusters)
    ax2.scatter(
        np.array(x)[:, 0], np.array(x)[:, 1], marker=".", s=30, lw=0, alpha=0.7, c=colors, edgecolor="k"
    )

    # Labeling the clusters
    centers = clusterer.cluster_centers_
    # Draw white circles at cluster centers
    ax2.scatter(
        centers[:, 0],
        centers[:, 1],
        marker="o",
        c="white",
        alpha=1,
        s=200,
        edgecolor="k",
    )

    for i, c in enumerate(centers):
        ax2.scatter(c[0], c[1], marker="$%d$" % i, alpha=1, s=50, edgecolor="k")

    ax2.set_title("The visualization of the clustered data.")
    ax2.set_xlabel("Feature space for the 1st feature")
    ax2.set_ylabel("Feature space for the 2nd feature")

    plt.suptitle(
        "Silhouette analysis for KMeans clustering on sample data with n_clusters = %d"
        % n_clusters,
        fontsize=14,
        fontweight="bold",
    )

plt.show()

For n_clusters = 2 The average silhouette_score is : 0.6810461692117465
For n_clusters = 3 The average silhouette_score is : 0.5511916046195927
For n_clusters = 4 The average silhouette_score is : 0.49535632852885064
For n_clusters = 5 The average silhouette_score is : 0.48989824728439524
For n_clusters = 6 The average silhouette_score is : 0.47711750058213453

Hierarchical Clustering

Hierarchical clustering is a general family of clustering algorithms that build nested clusters by merging or splitting them successively. The point of hierarchical clustering is to organize observations that are close together, and separate them out into groups/clusters. Closeness is generally a measure of distance between observations, the primary measures being Euclidean, Manhattan or Cosine. You have to pick the one that makes sense for your situation.

For most uses, Euclidean distance (often the default) does a great job. Cosine distances are more useful when doing natural language analysis.

Agglomerative Clustering
Agglomerative Clustering is a bottom up approach: each observation starts in its own cluster, and closest clusters are successively merged together. The ‘linkage criteria’ is a parameter passed to the sklearn function for performing the clustering.
- Single linkage (default) minimizes the distance between the closest observations of pairs of clusters.
- Ward minimizes the sum of squared differences within all clusters.
- Maximum or complete linkage minimizes the maximum distance between observations of pairs of clusters.
- Average linkage minimizes the average of the distances between all observations of pairs of clusters.

Agglomerative cluster has a "rich get richer" behavior that leads to uneven cluster sizes. In this regard, single linkage is the worst strategy, and Ward gives the most regular sizes.

Example
Recall the mtcars dataset that has numerical information on 32 models of cars. Let us apply agglomerative clustering to it. We will first standardize, or rescale the data, to make sure no individual feature overwhelms the other due to its scale. Then we will run the clustering algorithm, and present the result as a dendrogram.

mtcars = sm.datasets.get_rdataset('mtcars').data
data = mtcars.iloc[:, :]
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
data_scaled = scaler.fit_transform(data) 
data_scaled = pd.DataFrame(data_scaled, columns=data.columns)
data_scaled.head()