Session 2: Classification and Clustering

Getting started

First, we will load the packages we need for these exercises.

# This line is needed to force matplotlib to display inline in the notebook
%matplotlib inline

# The 3 lines below here suppress ConvergenceWarnings -- this is not necessarily something
# I would recommend doing in practice, but avoids spammed warnings in a particular cell
# towards the end of this file
from warnings import simplefilter
from sklearn.exceptions import ConvergenceWarning
simplefilter("ignore", category=ConvergenceWarning)

import os
import pickle
import tempfile

import polars as pl                        # Work with data frames
import numpy as np                         # Simple mathematics function and linear algebra
import matplotlib.pyplot as plt            # Charting functions
plt.rcParams['figure.figsize'] = [12, 8]   # Make plots larger by default
from sklearn import cluster                # For using kmeans
from sklearn import linear_model           # GLM-type models -- both econometric and ML
from sklearn import model_selection        # Training and testing splits
from sklearn import metrics                # Model performance evaluation
from sklearn import neighbors              # For KNN
from sklearn import preprocessing          # For standardizing data for LASSO and Elastic net
from sklearn import svm                    # For SVM-type models
import xgboost as xgb                      # For xgboost models

# Numba cannot cache beside packages installed on some mapped/network drives.
# Set a writable cache location before importing UMAP.
os.environ.setdefault('NUMBA_CACHE_DIR', os.path.join(tempfile.gettempdir(), 'numba_cache'))
import umap.umap_ as umap                  # Needed for the custom function `umap_compare_svm()`

Next, we need to import the dataset for the sessions exercises.

We will also split the data file like we did for Session 1 and apply the needed transformations to it using sklearn.

# Load data
# Scan the full file when inferring types because some numeric columns begin with integers
# but contain floating-point values later in the file.
df = pl.read_parquet('../../Data/S1_data.parquet')

# Define variable sets
vars_financial = ['logtotasset', 'rsst_acc', 'chg_recv', 'chg_inv', 'soft_assets', 'pct_chg_cashsales', 'chg_roa',
                  'issuance', 'oplease_dum', 'book_mkt', 'lag_sdvol', 'merger', 'bigNaudit', 'midNaudit', 'cffin',
                  'exfin', 'restruct']
vars_style = ['bullets', 'headerlen', 'newlines', 'alltags', 'processedsize', 'sentlen_u', 'wordlen_s', 'paralen_s',
              'repetitious_p', 'sentlen_s', 'typetoken', 'clindex', 'fog', 'active_p', 'passive_p', 'lm_negative_p',
              'lm_positive_p', 'allcaps', 'exclamationpoints', 'questionmarks']
vars_topic = ['Topic_' + str(i+1) + '_n_oI' for i in range(0,31)]

# Subset the final year to be the testing year
train = df.filter(pl.col('year') < 2004)
test = df.filter(pl.col('year') == 2004)

# Set up the data for the linear problem
vars_linear = vars_topic
scaler_X = preprocessing.StandardScaler()
scaler_X.fit(train.select(vars_linear).to_numpy())
train_X_linear = scaler_X.transform(train.select(vars_linear).to_numpy())
test_X_linear = scaler_X.transform(test.select(vars_linear).to_numpy())

scaler_Y = preprocessing.StandardScaler()
train_sdvol = train.get_column('sdvol1').to_numpy().reshape(-1, 1)
test_sdvol = test.get_column('sdvol1').to_numpy().reshape(-1, 1)
scaler_Y.fit(train_sdvol)
train_Y_linear = scaler_Y.transform(train_sdvol)
test_Y_linear = scaler_Y.transform(test_sdvol)

# Set up the data for the binary classification problem
vars_logistic = vars_topic + vars_financial + vars_style
scaler_X = preprocessing.StandardScaler()
scaler_X.fit(train.select(vars_logistic).to_numpy())
train_X_logistic = scaler_X.transform(train.select(vars_logistic).to_numpy())
test_X_logistic = scaler_X.transform(test.select(vars_logistic).to_numpy())

train_Y_logistic = train.get_column('Restate_Int').to_numpy()
test_Y_logistic = test.get_column('Restate_Int').to_numpy()

Below I have also defined some custom functions. Don't worry about these until we get to using them.

# From umap.plot source code on Github
def _get_embedding(umap_object):
    if hasattr(umap_object, "embedding_"):
        return umap_object.embedding_
    elif hasattr(umap_object, "embedding"):
        return umap_object.embedding
    else:
        raise ValueError("Could not find embedding attribute of umap_object")


def umap_color(data_map, data_color, cmap='viridis', subset=None, title=None):
    """Plot a UMAP embedding colored by numeric or categorical labels."""
    data_map = np.asarray(data_map)
    data_color = np.asarray(data_color)
    embed = _get_embedding(umap.UMAP(random_state=42, n_jobs=1).fit(data_map))

    if subset is not None:
        embed = embed[subset]
        data_color = data_color[subset]

    fig, ax = plt.subplots(figsize=(12, 8))
    point_size = 100.0 / np.sqrt(len(embed))
    if np.issubdtype(data_color.dtype, np.number):
        points = ax.scatter(embed[:, 0], embed[:, 1], s=point_size, c=data_color, cmap=cmap)
        fig.colorbar(points, ax=ax)
    else:
        categories = np.unique(data_color)
        colors = plt.get_cmap(cmap)(np.linspace(0, 1, len(categories)))
        for category, color in zip(categories, colors):
            mask = data_color == category
            ax.scatter(embed[mask, 0], embed[mask, 1], s=point_size, color=color, label=category)
        ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left')

    ax.set(xticks=[], yticks=[], title=title)
    return ax


# Cut down version of umap.plot.points to remove dependencies on  datashader, bokeh, holoviews, scikit-image, and colorcet
def umap_compare_svm(X, Yhat, Y, clip = None, cmap='viridis', subset=None, binary=False, title=None):
    reducer = umap.UMAP()
    umap_object = reducer.fit(X)
    embed = _get_embedding(umap_object)
    if clip is not None:
        Yhat = np.clip(Yhat, clip[0][0], clip[0][1])
        Y = np.clip(Y, clip[1][0], clip[1][1])

    fig, (ax1, ax2) = plt.subplots(1, 2)
    if subset is not None:
        embed_X = embed[subset,0]
        embed_Y = embed[subset,1]
        Y = np.array(Y[subset])
        X = np.array(X[subset])
        Yhat = np.array(Yhat[subset])
    else:
        embed_X = embed[:, 0]
        embed_Y = embed[:, 1]
    
    point_size = 100.0 / np.sqrt(len(embed_X))
    
    if binary:
        point_size = point_size * (1 + Y * binary)
    
    # color by values for Yhat
    points = ax1.scatter(embed_X, embed_Y, s=point_size, c=Yhat, cmap=cmap)
    fig.colorbar(points, ax=ax1, orientation='horizontal')

    ax1.set(xticks=[], yticks=[])
    ax1.set_title("Predicted values")
    
    # color by values for Y
    points = ax2.scatter(embed_X, embed_Y, s=point_size, c=Y, cmap=cmap)
    
    fig.colorbar(points, ax=ax2, orientation='horizontal')

    ax2.set(xticks=[], yticks=[])
    ax2.set_title("Actual values")
    
    if title is not None:
        fig.suptitle(title)
    
    if clip is not None:
        foot = 'Predicted values winsorized to [{}, {}]; Actual values winsorized to [{}, {}]'.format(clip[0][0], clip[0][1], clip[1][0], clip[1][1])
        plt.figtext(0.2, 0.3, foot, horizontalalignment='left')
    
    return (ax1, ax2)


def logistic(x):
    return 1 / (1 + np.exp(-1 * x))


def coefplot(names, coef, title=None):
    # Make sure coef is list, cast to list if needed.
    if isinstance(coef, np.ndarray):
        if len(coef.shape) > 1:
            coef = list(coef[0])
        else:
            coef = list(coef)
    
    # Drop unneeded vars
    data = []
    for i in range(0, len(coef)):
        if coef[i] != 0:
            data.append([names[i], coef[i]])
    data.sort(key=lambda x: x[1])
    
    # Add in a key for the plot axis
    data = [data[i] + [i+1] for i in range(0,len(data))]
    
    fig, ax = plt.subplots(figsize=(4,0.25*len(data)))

    ax.scatter([i[1] for i in data], [i[2] for i in data])
    
    ax.grid(axis='y')
    ax.set(xlabel="Fitted value", ylabel="Residual", title=(title if title is not None else "Coefficient Plot"))
    
    ax.axvline(x=0, linestyle='dotted')
    ax.set_yticks([i[2] for i in data])
    ax.set_yticklabels([i[0] for i in data])
    
    return ax

SVM and SVR

In computer science and many other discplines, a simple type of classifier used in training models is SVM. This is particularly common for algorithms that are being used as supervised methods for predicting some variable in a model.

SVM: Support Vector Machine for Classification

This is the most common approach for these algorithms. The algorithm will allow us to approach a classification problem, and it can either report the most likely class for each observation, or, like with logistic regression, it can report the probability of belonging to a given class.

Multiple implementations are available in sklearn for this algorithm:

• sklearn.svm.LinearSVC(): Fast and memory efficient; assumes a linear kernel. Does not output probabilities!
• sklearn.svm.SVC(): More flexible as it allows other kernel functions, but less efficient when using a linear kernel
• sklearn.linear.SGDClassifier(): Can emulate sklearn.svm.LinearSVC() depending on specified parameters, allows some additional flexibility in parameters, is more memory efficient, and can do online (batch) training (good for very large data sets)

Since we are using a linear kernel in this example, we can also make sense of the coefficients assigned to each input in our model. We can actually use the exact same coefplot() as in Session 1 to visualize the algorithm. However, if we change the kernel, such as using the radial basis function kernel, then we cannot use coefplot() meaningfully. There are other ways to visualize this, however, as we will see later.

model_svc = svm.LinearSVC(C=1, dual=False)
model_svc.fit(train_X_logistic, train_Y_logistic)

LinearSVC(C=1, dual=False)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

display = metrics.RocCurveDisplay.from_estimator(model_svc, test_X_logistic, test_Y_logistic)
display.plot()

<sklearn.metrics._plot.roc_curve.RocCurveDisplay at 0x18893b823c0>

coefplot(vars_logistic, model_svc.coef_)

<Axes: title={'center': 'Coefficient Plot'}, xlabel='Fitted value', ylabel='Residual'>

Optimizing C

The above classifier just took the default of C=1. However, like with LASSO and related models, C is a hyperparameter. Unlike LASSO and related models, sklearn does not provide us with a built-in CV function. Instead, we will need to use sklearn's built-in cross validation constructors to customize a cross validation for this. This is built into the same sklearn.model_selection library that we used for splitting data in Session 1.

Note that the model_selection.GridSearchCV() function allows for parallel processing -- you can specify the number of processes to use using the n_jobs= parameter.

Also note that model_selection.StratifiedShuffleSplit() is a bit more complex than the standard cross validation we did in session 1, in that it is making stratified random samples across classes. This is rather important in our case, since there is pretty severe class imbalance (less than 3% of data is an irregularity). If you want to directly replicate the method from before, model_selection.KFold() will do the trick

C_range = np.logspace(-2, 6, 9)
param_grid = dict(C=C_range)
cv = model_selection.StratifiedShuffleSplit(n_splits=5, test_size=0.2, random_state=1)
grid_svc = model_selection.GridSearchCV(svm.LinearSVC(dual=False), param_grid=param_grid, cv=cv, n_jobs=20)
grid_svc.fit(train_X_logistic, train_Y_logistic)
print("The best parameters are %s with a score of %0.2f"
      % (grid_svc.best_params_, grid_svc.best_score_))

The best parameters are {'C': np.float64(0.01)} with a score of 0.99

display = metrics.RocCurveDisplay.from_estimator(grid_svc, test_X_logistic, test_Y_logistic)
display.plot()

<sklearn.metrics._plot.roc_curve.RocCurveDisplay at 0x188939c9bd0>

coefplot(vars_logistic, grid_svc.best_estimator_.coef_)

<Axes: title={'center': 'Coefficient Plot'}, xlabel='Fitted value', ylabel='Residual'>

Visualizing SVC

Given the high dimensionality of the SVC algorithm, it is difficult to visualize. To aid in this, the below visualizations implement a dimensionality reduction technique (UMAP) and show the resulting probabilities embedded into the 2-dimensional projection from UMAP. This will allow us to see if there are any clear patterns in our data as it relates to our dependent variable.

Note that LinearSVC does not provide probabilities. However, it reports a score that roughly approximates the output of a logistic regression via .decision_function(). Thus, applying a logistic (i.e., sigmoid) function of $f(x)=\frac{1}{1+e^{-x}}$ to the output will convert to probabilities. Note that there is some weighting internally applied within LinearSVC which we are unable to recover, introducing some potential error in the below charts.

Probabilities are fully recoverable under SVC(), however, SVC() uses a "hinge" loss function whereas LinearSVC() uses a "squared hinge" loss function. As the problem we are focused on is quite sensative to the choice of loss function (i.e., performs poorly under hinge loss), we will use the output from our linearSVC() above.

train_Yhat_logistic = logistic(grid_svc.decision_function(train_X_logistic))

umap_compare_svm(train_X_logistic, train_Yhat_logistic, train_Y_logistic,
                 clip=[[0.25, 0.3], [0, 1]], binary=5,
                title="Full sample")

(<Axes: title={'center': 'Predicted values'}>,
 <Axes: title={'center': 'Actual values'}>)

umap_compare_svm(train_X_logistic, train_Yhat_logistic, train_Y_logistic, clip=[[0.25, 0.3], [0, 1]], cmap='coolwarm', binary=5,
                 subset=((train_Y_logistic==1) | (np.random.rand(len(train_Y_logistic))<0.05)),
                title="Performance on actual irregularities (Large) and random sample of non-irregularities")

(<Axes: title={'center': 'Predicted values'}>,
 <Axes: title={'center': 'Actual values'}>)

umap_compare_svm(train_X_logistic, train_Yhat_logistic, train_Y_logistic, clip=[[0.25, 0.3], [0, 1]], cmap='coolwarm', binary=5,
                 subset=((train_Y_logistic==0) & (np.random.rand(len(train_Y_logistic))<0.05)),
                title="Performance on a random sample of non-irregularities")

(<Axes: title={'center': 'Predicted values'}>,
 <Axes: title={'center': 'Actual values'}>)

SVR: Support Vector Regression

Notes:

• I am specifying dual=False here because we have more observations than regressors. If you have more regressors than datapoints, set dual=True.

model_svr = svm.LinearSVR(C=1, dual=False, loss='squared_epsilon_insensitive')
model_svr.fit(train_X_linear, np.ravel(train_Y_linear))

LinearSVR(C=1, dual=False, loss='squared_epsilon_insensitive')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

coefplot(vars_linear, model_svr.coef_)

<Axes: title={'center': 'Coefficient Plot'}, xlabel='Fitted value', ylabel='Residual'>

C_range = np.logspace(-4, 6, 11)
param_grid = dict(C=C_range)
cv = model_selection.KFold(n_splits=5)
grid_svr = model_selection.GridSearchCV(svm.LinearSVR(dual=False, loss="squared_epsilon_insensitive"), param_grid=param_grid, cv=cv, n_jobs=8)
grid_svr.fit(train_X_linear, np.ravel(train_Y_linear))
print("The best parameters are %s with a score of %0.2f"
      % (grid_svr.best_params_, grid_svr.best_score_))

The best parameters are {'C': np.float64(0.0001)} with a score of 0.06

coefplot(vars_linear, grid_svr.best_estimator_.coef_)

<Axes: title={'center': 'Coefficient Plot'}, xlabel='Fitted value', ylabel='Residual'>

train_Yhat_linear = model_svr.predict(train_X_linear)

umap_compare_svm(train_X_linear, train_Yhat_linear, train_Y_linear, clip=[[0, 2], [0, 2]])

(<Axes: title={'center': 'Predicted values'}>,
 <Axes: title={'center': 'Actual values'}>)

XGBoost

This algorithm differs from the previous algorithms we have looked at primarily by being non-linear (unless you request it to be linear). Furthermore, it is tree-based and leverages a technique called "boosting" -- a type of ensembling.

XGBoost is a flexible algorithm supporting manyt functional forms (both regression structures and loss functions). It supports a variety of linear, binary, count, survival, multiclass, ranking, gamma, and tweedie forms.

Implementing XGBoost in python

Python has a very good library for XGBoost, aptly called xgboost. We will use this library to both fit and visualize our models.

Note: The below does not exhaust all possible parameters. There are also $L_2$ and $L_1$ regularization options as well as sampling parameters, among others. See a full list here.

The algorithm requires that data is passed to it in its own format. It can take our current datasets as input, and will create its desired format using the xgb.DMatrix() function.

dtrain = xgb.DMatrix(train_X_logistic, label=train_Y_logistic, feature_names=vars_logistic)
dtest = xgb.DMatrix(test_X_logistic, label=test_Y_logistic, feature_names=vars_logistic)

There are many parameters to set for the model. Every parameter has a default value, so you can focus on just those that are needed for your usage. Below is a sample starting point.

param = {
    'booster': 'gbtree',             # default -- tree based
    'nthread': 8,                    # number of threads to use for parallel processing
    'objective': 'binary:logistic',  # binary, output probabilities
    'eval_metric': 'auc',            # maximize ROC AUC
    'eta': 0.3,                      # shrinkage; [0, 1], default 0.3
    'max_depth': 6,                  # maximum depth of each tree; default 6
    'gamma': 0.1,                    # set above 0 to prune trees, [0, inf], default 0
    'min_child_weight': 1,           # higher leads to more pruning of tress, [0, inf], default 1
    'subsample': 0.8,                # Randomly subsample rows if in (0, 1), default 1
    'colsample_bytree': 0.8,         # Randomly subsample variables if in (0, 1), default 1
    'random_state': 70
}
num_round = 30

Next we can train the model using xgb.train().

model_xgb_logistic = xgb.train(param, dtrain, num_round)

We can check the ROC AUC just like we did with prior models.

test_Yhat_xgb_logistic = model_xgb_logistic.predict(dtest)
auc = metrics.roc_auc_score(test_Y_logistic, test_Yhat_xgb_logistic)
auc

0.5940013491775207

fpr, tpr, thresholds = metrics.roc_curve(test_Y_logistic, test_Yhat_xgb_logistic)
display = metrics.RocCurveDisplay(fpr=fpr, tpr=tpr, roc_auc=auc)
display.plot()

<sklearn.metrics._plot.roc_curve.RocCurveDisplay at 0x188adc32490>

To aid with interpretation, XGBoost also can output a chart showing the relative importance of each variable in the model. This can help us to see which data is important for the model. Note that since the model is non-linear, there are no signed coefficients. The chart simply tells us which variables have the greatest impact on its classification.

fig, ax = plt.subplots(figsize=(8,16))
xgb.plot_importance(model_xgb_logistic, ax=ax)

<Axes: title={'center': 'Feature importance'}, xlabel='Importance score', ylabel='Features'>

Seeing the trees

There are two ways to see the tress in the model:

1. Dump the model to a text file
2. Use graphviz to draw the trees (which you can output to the notebook or save to a file.

Both of these methods are shown below.

# If you want to see the final tree model, this will dump it into a file in the same folder as this notebook.
model_xgb_logistic.dump_model('dump.raw.txt')

# To run this, you need to:
#     1. Install the python graphviz package
#     2. Install graphviz from https://graphviz.org/download/
#         - MAKE SURE TO SELECT one of the "add to path" options during installation if on Windows!

xgb.to_graphviz(model_xgb_logistic, num_trees=0)

M:\Python_environments\miniconda\envs\MLSS2\Lib\site-packages\xgboost\plotting.py:268: FutureWarning: The `num_trees` parameter is deprecated, use `tree_idx` instead. 
  warnings.warn(

xgb.to_graphviz(model_xgb_logistic, num_trees=1)

M:\Python_environments\miniconda\envs\MLSS2\Lib\site-packages\xgboost\plotting.py:268: FutureWarning: The `num_trees` parameter is deprecated, use `tree_idx` instead. 
  warnings.warn(

xgb.to_graphviz(model_xgb_logistic, num_trees=2)

M:\Python_environments\miniconda\envs\MLSS2\Lib\site-packages\xgboost\plotting.py:268: FutureWarning: The `num_trees` parameter is deprecated, use `tree_idx` instead. 
  warnings.warn(

# If you want to view every tree contained in your final model, the below code will dump a PNG file of each tree
# into a "trees/" directory in the same folder as this file.
num_trees = len(model_xgb_logistic.get_dump())
for tree_index in range(0, num_trees):
    dot = xgb.to_graphviz(model_xgb_logistic, num_trees=tree_index)
    dot.format = 'png'
    dot.render("xgb_trees/tree{}".format(tree_index))

M:\Python_environments\miniconda\envs\MLSS2\Lib\site-packages\xgboost\plotting.py:268: FutureWarning: The `num_trees` parameter is deprecated, use `tree_idx` instead. 
  warnings.warn(

Optimizing parameters

Here we will start with the same parameters as above. We will then train a model based on iteratively optimizing various parameters:

1. max_depth and min_child_weight
2. eta
3. gamma
4. subsample and colsample_bytree
5. Number of rounds

The below example is reasonably rigorous for a properly trained model, but with some intuition applied to keep the search space smaller and to to keep execution time reasonable.

To make it easier to follow, we will use the same cross-validation methods as before, leaning on Scikit-learn. The XGBoost package comes with an interface to Scikit-learn, which is accessed via the xgb.XGBClassifier() function rather than the xgb.train() function.

param = {
    'booster': 'gbtree',             # default -- tree based
    'nthread': 8,                    # number of threads to use for parallel processing
    'objective': 'binary:logistic',  # binary, output probabilities
    'eval_metric': 'auc',            # maximize ROC AUC
    'eta': 0.3,                      # shrinkage; [0, 1], default 0.3
    'max_depth': 6,                  # maximum depth of each tree; default 6
    'gamma': 0.1,                    # set above 0 to prune trees, [0, inf], default 0
    'min_child_weight': 1,           # higher leads to more pruning of tress, [0, inf], default 1
    'subsample': 0.8,                # Randomly subsample rows if in (0, 1), default 1
    'colsample_bytree': 0.8,         # Randomly subsample variables if in (0, 1), default 1
    'random_state': 70
}
n_rounds = 20
random_states= [231534]

param_test = {
 'max_depth': [3, 4, 5],
 'min_child_weight': [1, 5, 10],
 'subsample': [0.7, 0.8, 0.9],
 'colsample_bytree': [0.8, 1.0],
 'eta': [0.05, 0.08, 0.1],
 'gamma': [0, 0.1, 0.3],
 }

del param['max_depth']
del param['min_child_weight']
del param['subsample']
del param['colsample_bytree']
del param['eta']
del param['gamma']

cv = model_selection.StratifiedShuffleSplit(n_splits=5, test_size=0.2, random_state=random_states[0])
search1 = model_selection.GridSearchCV(xgb.XGBClassifier(**param, n_estimators=n_rounds), 
             param_grid = param_test, scoring='roc_auc', n_jobs=12, cv=cv)
search1.fit(train_X_logistic,train_Y_logistic)
print(search1.best_params_, search1.best_score_)

param.update(search1.best_params_)

{'colsample_bytree': 1.0, 'eta': 0.05, 'gamma': 0, 'max_depth': 4, 'min_child_weight': 1, 'subsample': 0.8} 0.6834784818705524

num_round = 200
cv_results = xgb.cv(
    param, dtrain, num_round, nfold=10, stratified=False,
    early_stopping_rounds=50, as_pandas=False,
)
n_rounds = len(next(iter(cv_results.values())))
print('Stopping after {} rounds'.format(n_rounds))

Stopping after 28 rounds

final = xgb.XGBClassifier(**param, n_estimators=n_rounds)
final.fit(train_X_logistic,train_Y_logistic)

XGBClassifier(base_score=None, booster='gbtree', callbacks=None,
              colsample_bylevel=None, colsample_bynode=None,
              colsample_bytree=1.0, device=None, early_stopping_rounds=None,
              enable_categorical=True, eta=0.05, eval_metric='auc',
              feature_types=None, feature_weights=None, gamma=0,
              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=4, max_leaves=None,
              min_child_weight=1, missing=nan, monotone_constraints=None,
              multi_strategy=None, n_estimators=28, n_jobs=None, ...)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

metrics.RocCurveDisplay.from_estimator(final, test_X_logistic, test_Y_logistic)

<sklearn.metrics._plot.roc_curve.RocCurveDisplay at 0x188ccedae90>

fig, ax = plt.subplots(figsize=(8,16))
xgb.plot_importance(final, ax=ax)
locs, labels = plt.yticks()
new_labels = []
for i in locs:
    new_labels.append(vars_logistic[int(labels[i].get_text()[1:])])
_ = plt.yticks(locs, new_labels)  # The `_ = ` here captures the excessive matplotlib output and deletes it

with open('../../Data/S2_models.pkl', 'wb') as f:
    pickle.dump({'SVC': grid_svc, 'XGBoost': final}, f)

Clustering

Kmeans

topic_names = ['Topic_' + str(i) + '_n_oI' for i in range(1, 32)]
industry_from_sic = (
    pl.when(pl.col('sic').is_between(0, 999)).then(pl.lit('Agriculture'))
    .when(pl.col('sic').is_between(1000, 1499)).then(pl.lit('Mining'))
    .when(pl.col('sic').is_between(1500, 1799)).then(pl.lit('Construction'))
    .when(pl.col('sic').is_between(2000, 3999)).then(pl.lit('Manufacturing'))
    .when(pl.col('sic').is_between(4000, 4999)).then(pl.lit('Utilities'))
    .when(pl.col('sic').is_between(5000, 5199)).then(pl.lit('Wholesale Trade'))
    .when(pl.col('sic').is_between(5200, 5999)).then(pl.lit('Retail Trade'))
    .when(pl.col('sic').is_between(6000, 6799)).then(pl.lit('Finance'))
    .when(pl.col('sic').is_between(7000, 8999)).then(pl.lit('Services'))
    .when(pl.col('sic').is_between(9100, 9999)).then(pl.lit('Public Admin'))
    .otherwise(pl.lit('Unknown'))
    .alias('industry')
)
train = train.with_columns(industry_from_sic)

train.select('industry', 'sic')

shape: (11_478, 2)

industry	sic
str	i64
"Wholesale Trade"	5080
"Wholesale Trade"	5080
"Manufacturing"	3661
"Manufacturing"	3661
"Manufacturing"	2834
…	…
"Retail Trade"	5531
"Wholesale Trade"	5040
"Manufacturing"	2835
"Services"	8051
"Services"	7370

umap_color(
    train.select(topic_names).to_numpy(),
    train.get_column('industry').to_numpy(),
    title='Industries',
)

<Axes: title={'center': 'Industries'}>

model = cluster.KMeans(n_clusters=9, n_init='auto')
kmeans = model.fit(train.select(topic_names).to_numpy())
train = train.with_columns(pl.Series('cluster', kmeans.labels_))

umap_color(
    train.select(topic_names).to_numpy(),
    train.get_column('cluster').cast(pl.String).to_numpy(),
    title='K-means clusters',
)

<Axes: title={'center': 'K-means clusters'}>

Optimizing kmeans

The Gap statistic provides us a way to optimize the number of clusters for kmeans in a way that is unsupervised and statistically grounded. It is simulation based though (bootstrap-ish), which means it will take a while to run.

iterations = 10
ks = []
gaps = []
sks = []
inertias = []

k = 2
init_k = k
optimizing = True
while optimizing:
    model = cluster.KMeans(n_clusters=k, n_init='auto')
    kmeans = model.fit(train.select(topic_names).to_numpy())
    inertia_at_k = kmeans.inertia_
    # run 50 iterations to determine s.d.
    sim_inertias = []
    for i in range(0,iterations):
        model = cluster.KMeans(n_clusters=k, n_init='auto')
        # Simulate on random data
        kmeans = model.fit(np.random.rand(train.height, len(topic_names)))
        sim_inertias.append(kmeans.inertia_)
    l = np.mean(np.log(sim_inertias))
    gap = l - np.log(inertia_at_k)
    sk = (sum((1/iterations) * (np.log(sim_inertias) - l)**2)**(1/2)) * (1+1/iterations)**(1/2)
    ks.append(k)
    sks.append(sk)
    gaps.append(gap)
    inertias.append(inertia_at_k)
    if k > init_k:
        thresh = gap - gaps[-2] - sks[-2]
        if thresh < 0:
            print('Optimal k found: ' + str(k-1))
            print(thresh)
            break
        else:
            print('k={} Not optimal.  To optimize: {:.2f}'.format(k-1, thresh))
            k += 1
    else:
        k += 1

k=2 Not optimal.  To optimize: 0.13
k=3 Not optimal.  To optimize: 0.09
k=4 Not optimal.  To optimize: 0.09
k=5 Not optimal.  To optimize: 0.08
k=6 Not optimal.  To optimize: 0.04
k=7 Not optimal.  To optimize: 0.07
k=8 Not optimal.  To optimize: 0.07
k=9 Not optimal.  To optimize: 0.05
k=10 Not optimal.  To optimize: 0.05
k=11 Not optimal.  To optimize: 0.04
k=12 Not optimal.  To optimize: 0.03
k=13 Not optimal.  To optimize: 0.04
k=14 Not optimal.  To optimize: 0.00
k=15 Not optimal.  To optimize: 0.06
k=16 Not optimal.  To optimize: 0.03
Optimal k found: 17
-0.011332538640951958

print('Using optimal k = ' + str(k-1))
model = cluster.KMeans(n_clusters=k-1, n_init='auto')
kmeans = model.fit(train.select(topic_names).to_numpy())
train = train.with_columns(pl.Series('cluster_opt', kmeans.labels_))

Using optimal k = 17

tables = [
    train.filter(pl.col('cluster_opt') == cluster_id).select('industry').sample(n=10, seed=42)
    for cluster_id in (1, 2, 3)
]
with open('../Results/S2_cluster_tables.pkl', 'wb') as f:
    pickle.dump(tables, f, pickle.HIGHEST_PROTOCOL)

umap_color(
    train.select(topic_names).to_numpy(),
    train.get_column('cluster_opt').cast(pl.String).to_numpy(),
    title='Optimized K-means clusters',
)

<Axes: title={'center': 'Optimized K-means clusters'}>

KNN for multiclass prediction

To set up this problem, we need to apply our industry classification to the testing data as well.

test = test.with_columns(industry_from_sic)

To set up a KNN model, we will use the KNeighborsClassifier() function from Scikit Learn. The most important parameter is n_neighbors, which corresponds to the k in KNN -- how many neighbors to use for classification. There are other parameters for the Scikit Learn implementation as well, which can be seen in the documentation.

knn = neighbors.KNeighborsClassifier(n_neighbors=5)
knn.fit(
    train.select(topic_names).to_numpy(),
    train.get_column('industry').to_numpy(),
)

KNeighborsClassifier()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

To see how well the algorithm works, we can make multiclass predictions using .predict(). We can then pass the predictions to sklearn.metrics.accuracy_score() to get accuracy as a percentage.

in_pred = knn.predict(train.select(topic_names).to_numpy())
out_pred = knn.predict(test.select(topic_names).to_numpy())

print('In sample: {},\nOut of sample: {}'.format(
    metrics.accuracy_score(train.get_column('industry').to_numpy(), in_pred),
    metrics.accuracy_score(test.get_column('industry').to_numpy(), out_pred),
))

In sample: 0.9123540686530754,
Out of sample: 0.8597236981934112

To optimize parameters for KNN, we can use a Grid Search. It is fairly efficient to run, so we can do a full grid search.

knn_opt = neighbors.KNeighborsClassifier(algorithm='auto')
knn_param = {
    'n_neighbors': [1,2,3,4,5,6,7,8,9,10],
    'leaf_size': [10, 20, 30, 40],
    'p': [1, 2],
    'weights': ['uniform', 'distance'],
}
                   
# with GridSearch
grid_knn = model_selection.GridSearchCV(
    estimator=knn_opt,
    param_grid=knn_param,
    scoring = 'accuracy',
    n_jobs = -1,
    cv = 5
)

grid_knn.fit(
    train.select(topic_names).to_numpy(),
    train.get_column('industry').to_numpy(),
)

GridSearchCV(cv=5, estimator=KNeighborsClassifier(), n_jobs=-1,
             param_grid={'leaf_size': [10, 20, 30, 40],
                         'n_neighbors': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
                         'p': [1, 2], 'weights': ['uniform', 'distance']},
             scoring='accuracy')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

params = grid_knn.best_params_
print(params) 

{'leaf_size': 10, 'n_neighbors': 10, 'p': 1, 'weights': 'distance'}

knn2 = neighbors.KNeighborsClassifier(**params)
knn2.fit(
    train.select(topic_names).to_numpy(),
    train.get_column('industry').to_numpy(),
)

KNeighborsClassifier(leaf_size=10, n_neighbors=10, p=1, weights='distance')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

in_pred = knn2.predict(train.select(topic_names).to_numpy())
out_pred = knn2.predict(test.select(topic_names).to_numpy())

print('In sample: {},\nOut of sample: {}'.format(
    metrics.accuracy_score(train.get_column('industry').to_numpy(), in_pred),
    metrics.accuracy_score(test.get_column('industry').to_numpy(), out_pred),
))

In sample: 1.0,
Out of sample: 0.8852284803400637