<jupyter_start><jupyter_text>Summary*Author: <NAME>*In this notebook we will finally start getting to the exciting stuff: machine learning! We did the necessary preprocessing, splitting, and feature engineering and now it's time to start the iterative building process. As mentioned in the previous notebook, rather than creating all possible useful features, we will begin with the few we constructed and build a very basic classifier. Depending on how well that does we may decide we need to go back to the drawing board and either collect more data, design better features, or choose a different model.This notebook will introduce a few topics which you may have seen before but we feel are particularly useful in an applied setting. These topics include:- Cross validation- The notion of a test set score and a training set score as they relates to bias and variance- Various performance metrics that are utilized both in machine learning circles as well as clinical circles  Learning objectivesBy the time you work through this notebook you should be able to:- Use cross validation for the purposes of hyperparameter tuning on the development (training) set- Understand the implications of disparities in training and testing set performance for model choice and hyperparameter tuning- Conceptually understand why regularization is often necesassary for prediction tasks- Implement a basic machine learning pipeline for a clasification task  InitializeDependencies. . .<jupyter_code>import pandas as pd
import numpy as np
import os
import matplotlib.pyplot as plt
import utils # Utility functions containained in ./utils/
import scipy
import seaborn as sns

from IPython.display import Image
from IPython.core.display import HTML 

%matplotlib inline
%load_ext autoreload
%autoreload 2

pd.set_option('max.rows', 100)
pd.set_option('max.columns', 100)
sns.set_style('darkgrid')

# Change this location to the path where you would like your data saved to
data_dir = '/Users/bbradshaw/nhanes/'

# Path to hdf store we will create later
hdf_path = 'nhanes_data.h5'
# Read in only training set features for now
train_features = pd.read_hdf(os.path.join(data_dir, hdf_path), 'train_features')
train_features.head()<jupyter_output><empty_output><jupyter_text>Building a machine learning pipelineNow it's time to get our hands dirty. In order to build a functional pipeline there are several design choices we need to consider:- How do we handle missing values?- How do we choose a model?- How do we choose model hyperparameters?We'll go over each of these, but first let's talk about cross validation.  Cross validationCross validation is a powerful tool for both building a model that generalizes to unseen examples, as well as tuning hyperparameters and slecting a model. The idea is rather than splitting our dev set into a separate training and validation set (like we did for the development-testing set) we will use the entire development set and split it into k disjoint folds. We will then train the model on (k-1) of the folds, and use the kth fold for validation. We will repeat this process k times, each time using a different fold for validation. The benefits of this method are twofold. First, this process is much more efficient in term sof data usage. Essentially we are using the entire dataset (at one point or another) to formulate generalization estimates. Second, this method yields a *distribution* over metrics of interest. For example, suppose we are interested in comparing the AUROC of several models. Since many models contain stochastic components, a comparison of a single train-validation split may give us misleading information. However, computing k AUROC scores for each model allows us to build bootstrapped confidence intervals around the true value of interest.The image below gives a visual explanation of how cross validation works.<jupyter_code>Image(filename='../supplementary_materials/cv_image.jpg')<jupyter_output><empty_output><jupyter_text>Constructing a preprocessign pipelinePretty much all machine learning models require all values for all features to be specified. For that to be the case, we need to impute missing values before they are input to the model. In order to do this we need to talk a bit about feature types. Take a look at the figure blow. It gives a break down of the two main types of data we will need to address for input into a model.<jupyter_code>Image(filename='../supplementary_materials/data_types.jpg')
# Let's look at our variables and see what dtypes
# they map to
train_features.info()
# Map variables to types

# Features to exclude from the model
excl_feats = ['seqn', 'day']

# Numeric features
numeric_feats = [
    'first_step_minute',
    'last_step_minute',
    'max_steps',
    'sum_steps',
    'pct_nonzero_steps',
    'max_rolling_30',
    'first_morning_step_minute',
    'age_in_years'
]

# Categorical features
categorical_feats = [
    'gender',
    'currently_in_school',
    'household_size'
]

# Dummy variable names
dummy_feats = list(
    pd.get_dummies(
        train_features[categorical_feats]
    ).columns
)
# Check how many null values there are for each feature
train_features.isnull().sum()
# Turn categorical variables into strings for one-hot encoding
for feat in categorical_feats:
    train_features[feat] = train_features[feat].astype('category')<jupyter_output><empty_output><jupyter_text>You can see most features have no missing values, but there are a few. We need to come up with a strategy for filling these values. There are two items to note here:- Typically the strategy you use for imputing numeric variables is not the same as for computing categorical variables- When imputing values we almost always impute missing values in the test set with statistics computed in the training set. This prevents us from a very nefarious situation known as *label leaking* which we will talk about a bit later.For simplicity, we will use the median value of each numeric feature to fill missing values, and the mode of each feature to fill categorical features. In practice there are much more complicated schemes, and you can select the right strategy by using cross validation. However, this choice is usually a reasonable baseline and tends to work just fine.<jupyter_code># These are utilities for preprocessing- you should 
# review them for your own code edification
from utils import DataFrameSelector, DummyGenerator
from sklearn.pipeline import Pipeline, FeatureUnion
from sklearn.preprocessing import Imputer, StandardScaler
# The sklearn Pipeline class allows us to string together
# transformations into a single flow,
# Here we construct two pipelines: one for numeric features
# and one for categorical features. We then union them together
# into a single pipeline.

num_pipeline = Pipeline([
    ('selector', DataFrameSelector(numeric_feats)),
    ('imputer', Imputer(strategy='median')),
    ('std_scaler', StandardScaler())
])

cat_pipeline = Pipeline([
    ('dummy_generator', DummyGenerator(categorical_feats))
    ])

full_pipeline = FeatureUnion(transformer_list=[
    
        ("num_pipeline", num_pipeline),
        ("cat_pipeline", cat_pipeline),
    ])<jupyter_output><empty_output><jupyter_text>The above code isn't long but it's procedurally dense. It abstracts A LOT of the process of selecting the correct datatypes, applying an appropriate transform, and then glueing everything back together. Since this is meant to illustrate more the process than the details, I won't review what's going on under the hood in gory detail. With that being said: I highly recommend you take the time to review both the sklearn docs to better understand what each of the pieces of functionality do individually. Ideally, you should code these up yourself for practice.  Defining a performance metricPhew! We have finally arrived! We downloaded the data, did QC checks, constructed features, and built a preprocessing pipeline. It's time to test a model! Remember: building a machine learning system is an iterative process. We will implement a very basic model, and we will measure its performance. If it does "good enough" for the use case at hand then great, we are done! If not, we have decisions to make: do we need better features? Do we have enough/the right data? Is our model too simple? Too complex? These are some of the questions we will need to address as we build our system.How does one go about defining "good enough"? A lot of this depends on the specific use case. Suppose you were building a classifier that detected a rare disease. Would you care more about limiting false positives or limiting false negatives? Listed below are several important definitions you should be familiar with. They take some time to master but are foundational for reasoning about the performance of a machine learning system. Most of these metrics are functions of the *confusion matrix* which designates the distribution of a model's predictions across the true class labels. The figure below gives a nice overview.- **Accuracy:** The proportion of cases your algorithm accurately classified (TP + TN) / (TP + FP + TN + FN)- **Recall (sensitivity):** The proportion of positive cases your algorithm correctly identified as such (TP) / (TP + FN)- **Precision (positive predictive value):** The proportion of cases your algorithm identified as the positive case, that were correct (TP) / (TP + FP)- **Specificity**: The proportion of negative classes your algorithm correctly identified as such (TN) / (TN + FP)- **Negative predictive value:** The proportion of cases your algorithm identified as the negative class, that were correct (TN) / (TN + FN)<jupyter_code>Image('../supplementary_materials/confusion_matrix.png')
# First step - separate out y (our target) and X (our feature matrix)

# is_weekend is our target
y = train_features.is_weekend.values

# Extract out the features we care about as defined above
X = train_features[numeric_feats+categorical_feats].values

# Extract groups so we keep indivual peron data in the same fold
groups = train_features.seqn.values
# Classification metrics for validation
from sklearn.metrics import accuracy_score, precision_score, recall_score, roc_auc_score

# This is a helper function that abstracts the cross validation
# implementation
from sklearn.model_selection import cross_val_score
# Models to try out
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
# Random forest is a stochastic algorithm - let's set a seed for 
# reproduceability

# Nothing special about this number
random_seed = 42
# For now use pretty much the set defaults, except
# use class_weight-'balanced' since our classes are imbalanced 4:1
rf = RandomForestClassifier(random_state=42, class_weight='balanced', n_jobs=-1)
# Fit and transform features with our preprocessing pipeline
X_clean = full_pipeline.fit_transform(train_features)
# Compute a few of the metrics mentioned above
metrics = {
    'accuracy': None,
    'recall': None,
    'precision': None,
    'roc_auc': None
}

for key in metrics.keys():
    # Note this is inefficient: I am training the model (5 folds) x (4 metrics) = 20 times
    metrics[key] = cross_val_score(rf, X_clean, y, scoring=key, cv=5)
for m in metrics:
    print('metric:', m, '\n')
    print('mean:', metrics[m].mean(), '\n')
    print('std:', metrics[m].std(), '\n')
    print('---------------------------------------------------')<jupyter_output>metric: accuracy 

mean: 0.7780432432838047 

std: 0.004852076491476345 

---------------------------------------------------
metric: recall 

mean: 0.14444642063139174 

std: 0.018782187886612723 

---------------------------------------------------
metric: precision 

mean: 0.44631213869503084 

std: 0.034499310133308606 

---------------------------------------------------
metric: roc_auc 

mean: 0.6834217628288897 

std: 0.016201474979277937 

---------------------------------------------------<jupyter_text>Are the above metrics good? Let's think about it. Across folds my model has a mean accuracy of almost 78%: which at first glance seems pretty good. But on second though what would the accuracy of an extremely dumb model that always guess not weekend be? If I always guessed 0, then my sensitivity would be 100% and I would correctly identify all non-weekend days. Since the classes are unbalanced by definition (weekends vs weekdays) our accuracy will always look pretty good. Looking at the positive class proportion:<jupyter_code># Accuracy looks pretty good! Recall not so much. . . .
1 - y.sum() / len(y)<jupyter_output><empty_output><jupyter_text>The AUROC MetricThe above analysis tells us that in some ways our model is better than a naive model, in some ways it really isn't. How do we compare models using a metric that shows the trade offs in performance along different levels of operation? Perhaps the most wel known metric is the *Receiver Operating CHaracteristic Curve* which plots the true postive rate (recall or sensitivity) against the false postive rate (1-specificity) for various decision threshold boundaries. The area under this curve is known as AUROC (not suprisingly). A value of 0.5 implies that the classifier is no better than random, while  value of 1 indicates the classifier is able to perfectly discriminate between classes.<jupyter_code># This utility builds the ROC curve with k folds
# computing the AUROC for each held out set
from utils import plot_roc_with_cv
# Random forest AUROC is significantly better than random but not excellent.
plot_roc_with_cv(X=X_clean, y=y, classifier=rf, groups=groups, cv=10)<jupyter_output><empty_output><jupyter_text>Hyperparameter tuningNow that we have two candidate models we want to tune our models to get more performance out of them. There is a powerful tool for determining which hyperparameters are best and it again involves cross validation. This tool is called randomized search. The basic idea is to specify a distribution over the parameters that you want to try. The randomized search will then sample values over this multidimensional space and keep track of the combination that performs best using cross validation. Let's give it a shot.<jupyter_code>from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import randint

rf_param_distribs = {
    'n_estimators': randint(low=1, high=500),
    'max_features': randint(low=1, high=11),
    'max_depth': randint(low=1, high=15),
    'min_samples_split': randint(low=2, high=100)
    }

rf = RandomForestClassifier(
    random_state=random_seed,
    n_jobs=-1,
    class_weight='balanced'
)

rnd_search = RandomizedSearchCV(
    rf,
    param_distributions=rf_param_distribs,
    n_iter=50,
    cv=5,
    scoring='roc_auc',
    random_state=random_seed
)

rnd_search.fit(X_clean, y)
# Save the best estimator in a new object
best_rf = rnd_search.best_estimator_
# How much better is our AUROC with the tuned model?
plot_roc_with_cv(X=X_clean, y=y, classifier=best_rf, groups=groups, cv=10)
# Other metrics of interest
metrics = {
    'accuracy': None,
    'recall': None,
    'precision': None,
    'roc_auc': None
}

for key in metrics.keys():
    # Note this is inefficient: I am training the model (5 folds) x (4 metrics) = 20 times
    metrics[key] = cross_val_score(best_rf, X_clean, y, scoring=key, cv=5)
for m in metrics:
    print('metric:', m, '\n')
    print('mean:', metrics[m].mean(), '\n')
    print('std:', metrics[m].std(), '\n')
    print('---------------------------------------------------')<jupyter_output>metric: accuracy 

mean: 0.6939518341950816 

std: 0.01843571675808342 

---------------------------------------------------
metric: recall 

mean: 0.6321327997628575 

std: 0.023711093610297877 

---------------------------------------------------
metric: precision 

mean: 0.37452572081510777 

std: 0.019506193531910975 

---------------------------------------------------
metric: roc_auc 

mean: 0.7464932640545777 

std: 0.0077050108411386535 

---------------------------------------------------<jupyter_text>Tuning the model hyperparameters gave us a solid jump in AUC. We went from a mean AUC of 0.68 to a mean AUC of 0.74. That's a significant lift for a minimal amount of work. Interestingly, accuracy actually fell by quite a bit. This is because the model is getting more aggressive in predicting the postive class: We went from a recall of 0.12 to a recall of 0.59: an enormous jump. Unfortunately this came at the cost of a reduced precision. This is typical when building a classifier: there is a tradeoff between precision and recall.  Feature ImportancesOften times and especially in clinical contexts the interpretability of the model is extremely important since we will present our results to clinicians. One of the nice things about random forests and logistic regression is they feature importances/weights respectively that give insight into how the model is making its decisions.<jupyter_code>feature_set = numeric_feats + dummy_feats
# Write a little function that extracts feature importances (if
# they exist) and sort them by the absolute magnitude value
def get_feature_importances(trained_model, feature_list):
    if hasattr(trained_model, 'feature_importances_'):
        importances = pd.DataFrame(
            {'importance': trained_model.feature_importances_},
            index=feature_list
        )
        return importances.sort_values(by='importance', ascending=False)
    elif hasattr(trained_model, 'coef_'):
        importances = pd.DataFrame(
            {
                'coefficient': trained_model.coef_[0],
                'abs_weight': abs(trained_model.coef_[0])
            },
            index=feature_list
        )
        return importances.sort_values(by='abs_weight', ascending=False).drop('abs_weight', axis=1)
    else:
        raise AttributeError('Model does not have feature importance or coefficient attribute!')
# Extract random forest feature importances
rf_feature_imp = get_feature_importances(best_rf, feature_set)
rf_feature_imp.plot(kind='barh', figsize=(10, 6));<jupyter_output><empty_output><jupyter_text>The feature importances are reassuring, we see that when you wake up in the morning is the most important factor the model decides when it determines whether or not the target is a weekend. A cool result!  System evaluation and diagnosticsNow that we have a few working prototypes that are promising, if we decide our performance is not good enough we need to think about where the best place is to invest our efforts in improving the system. We have a few options:- **Collect better data:** Maybe the data we have isn't sufficient for the task at hand- **Collect more data:** Maybe we simply haven't exposed the algorithm to enough data in order for it to learn the optimal decision rule- **Construct better features:** Perhaps the amount of information in the underlying signal is sufficient for the problem at hand, but our features don't capture that signal. Note that this is the allure of neural networks: given enough data they can learn the features on their own. Unfortunately, often times we don't have the kind of data necessary to learn these features.- **Choose a different model:** It's possible our model either is too simple to capture the complexity of the data, or it is too complex and is overfitting the training data.We won't go into those next steps in this notebook, but if we have time we may cover them at the tail end of the module.  Testing our modelAt the very end of our machine learning pipeline, after we are satisfied with the results of our iterative evaluation process we will finally make use of our held out test set. This is the moment of truth! If we correctly constructed our test set, conducted appropriate hyperparameter tuning with cross validation, we should expect that the AUC of the classifier will be similar on the held out test set to the AUC we observed during cross validation on the training set. With that being said, it's normal for our AUC on the training set is *slightly* lower on the test set due to some marginal model overfitting on the training set due to the hyperparameter tuning process.<jupyter_code># Load our previously saved test features from notebook 2
test_features = pd.read_hdf(os.path.join(data_dir, hdf_path), 'test_features')
test_features.head()
# Reapply our pipeline to the testing set
for feat in categorical_feats:
    test_features[feat] = test_features[feat].astype('category')

X_test_clean = full_pipeline.fit_transform(test_features)
# Extract "true labels"
y_test_true = test_features.is_weekend.values
# Retrain our model on the full training set
best_rf.fit(X=X_clean, y=y)
# Let's compute some predictions!
# Note that we need to pass the roc_auc metric function a
# probability, so we extract positive class probs using the
# "predict_proba" method
test_set_probs = [x[1] for x in best_rf.predict_proba(X_test_clean)]
from sklearn.metrics import roc_auc_score
# Deep breathe. . . 
roc_auc_score(y_true=y_test_true, y_score=test_set_probs)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08215740126905957, 'fraction_numerical': 0.01761530447959087, 'mean_word_length': 4.188943488943489, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 40, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '152554', 'license': 'MIT', 'n_tokens_mistral': 5748, 'n_tokens_neox': 5324, 'n_words': 2685, 'path': 'notebooks/2-Classifier-Construction.ipynb', 'repo_name': 'evidation-health/EvidationDataScienceModule', 'size': '21129'}