Lab 6 - Feature Engineering and Ensemble Models

Up until this point we've applied machine learning models to "tabular data", or 2-dimensional data objects that have a spreadsheet-like format, where each row represents a single observation, and each column represents a predictive feature. However, many modern data collection protocols do not innately produce tabular data, for example:

• Social media and news publishers produce textual data, where each observation is a sequence of words or characters.
• Wearable devices and other sensors produce time-series data, where each observation consists of many time-ordered measurements, often spanning multiple channels.
• Medical technologies and satellites produce image data, where each observation consists of spatially structured pixel intensities.
• Audio recording devices produce audio wave data, which can often be viewed as a specialized form of time-series data.

None of these data types can be used directly by the types of models we've encountered so far in the course (we'll talk about deep learning models soon), which raises an important question:

How can we create effective representations of raw data so that it can be used by machine learning models?

The answer to this question is feature engineering, the process of using domain expertise to create meaningful predictive features that can be used by models designed for tabular data.

## Load the standard libraries for later
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import sklearn
import math

Part 1 - Feature Engineering and Textual Data

The "clickbait" data set given below contains $n=3000$ news article headlines from "clickbait publishers", including Upworthy, Buzzfeed, Thatscoop, ViralNova, Scoopwhoop, and ViralStories, and $n=3000$ headlines from websites deemed not to publish clickbait, including The New York Times, WikiNews, The Guardian, and The Hindu.

clickbait = pd.read_csv("https://remiller1450.github.io/data/clickbait_headlines.csv", encoding='latin-1')

Our goal is to develop a model that can predict whether an article comes from a clickbait publisher (Label = 1) based upon the article's headline.

There are three broad approaches to engineering useful features:

1. Domain expertise - we might know from experience that clickbait headlines frequently start with a number (ie: "15 Things You'll Only Understand If..."), so we might create a binary feature indicating whether or not a headline begins with a number.
2. Pattern discovery - we might graph the frequencies of the top $k$ words that appear in clickbait and non-clickbait headlines, then based upon these differences we could create a clickbait keyword score for each headline.
3. Automated extraction - we might generate a Term Frequency-Inverse Document Frequency (TF-IDF) for each word appearing in the data, then use these scores as features.

This lab will focus on the first two approaches, but the third is something I'd encourage you to explore on the course project.

To begin, we'll split off training and testing sets:

from sklearn.model_selection import train_test_split

X = clickbait['Headline']
y = clickbait['Label']

train_X, test_X, train_y, test_y = train_test_split(X, y, test_size=0.2, random_state=7)

Next, we'll define our own function that creates a feature by applying the domain expertise example from earlier:

def first_word_q(input_text):
    first_word = input_text.str.split(' ').str[0]
    return 1*((first_word == "What") | (first_word == "Which") | (first_word == "Would") | (first_word == "How") | (first_word == "Can") | (first_word == "Are"))

We should always check to make sure our function is vectorized, meaning it gets applied to each headline when given train_X:

first_word_q(train_X)

3033    0
2016    0
4025    1
4356    0
2801    0
       ..
5699    0
2550    0
537     0
1220    0
4271    0
Name: Headline, Length: 4800, dtype: int64

This confirms that first_word_q() producing output that will be easy to assemble into a data frame of engineered features.

Below is an example of a function that requires some extra effort due to it not being vectorized (you'll get an error if you try running it on train_X):

def count_US_Europe(input_text):
    return input_text.count('US') + input_text.count('U.S.') + input_text.count('UK') + input_text.count('U.K.') +  input_text.count('EU') + input_text.count('E.U.') 

For situations like this, the apply() method of pandas DataFrames can be useful:

train_X.apply(count_US_Europe)

3033    0
2016    0
4025    0
4356    0
2801    0
       ..
5699    0
2550    0
537     0
1220    0
4271    0
Name: Headline, Length: 4800, dtype: int64

Once we've engineered enough features to properly represent the predictive patterns in the data, we can assemble a data frame of engineered features to use in the subsequent steps of a machine learning pipeline:

train_X_FE = pd.DataFrame({
    'First_Word_Q': first_word_q(train_X),
    'count_US_Europe': train_X.apply(count_US_Europe)
})

train_X_FE.head()

	First_Word_Q	count_US_Europe
3033	0	0
2016	0	0
4025	1	0
4356	0	0
2801	0	0

You should note that I don't expect you to be an expert on string processing in Python, and you're encouraged to look up useful string methods using this reference or other online resources. However, for the purposes of this lab and Homework 3, you should limit yourself to using only functions/methods from pandas, numpy, and sklearn.

Question #1:

• Part A: Write a function that creates a new feature that takes on a value of 1 if the first word in a headline is a number, and 0 otherwise. Check to make sure your function is vectorized and produces the intended output.
• Part B: Use the .sample() method to view a few different random samples clickbait and non-clickbait headlines from the training data. Based upon your observations, write a function that creates a new feature to represent one of the patterns you observed.
• Part C: Assemble a pandas DataFrame using the two features you created in Parts A and B, as well as the two features created using this sections examples. Fit a decision tree model with a maximum depth of 4 to the training data and display the fitted modeling using plot_tree().
• Part D: Briefly explain why classification accuracy is a reasonable performance metric for this analysis. You should consider the characteristics of the data and the interpretability of the metric.
• Part E: Suppose the decision tree from Part C is deemed sufficiently accurate on the training data (ie: assume that cross-validation showed good enough performance for it to be useful). Apply this classifier to the test set in order to obtain a final evaluation of the model's classification accuracy.

Part 2 - Feature Engineering and Time-Series Data

This part of the lab provides a second example of feature engineering, this time using time-series data collected by an accelerometer worn by $n=15$ subjects as they performed 7 different types of activities.

The data were recorded at 52 Hz, meaning that 52 measurements are present for every second of time recorded by the accelerometer recorded. Our analysis will focus on predicting a subject's activity using 10-second samples of accelerometer data, which are delineated by the variable Sample_Num.

act_ts = pd.read_csv("https://remiller1450.github.io/data/activity.csv")

For simplicity, we will build a model that uses data from Subject #1 for training, and Subject #7 for validation, and we will restrict our analysis to only two activities: "Work", which records when the subject is working at their desk, and "Talk", which records when the subject is talking while standing.

## Training Data
train_act_X = act_ts[act_ts['Subject'] == 1].drop(['Label'], axis = 1)
train_act_y = act_ts[act_ts['Subject'] == 1].groupby('Sample_Num').nth(0)['Label']

## Validation Data
test_act_X = act_ts[act_ts['Subject'] == 7].drop(['Label'], axis = 1)
test_act_y = act_ts[act_ts['Subject'] == 7].groupby('Sample_Num').nth(0)['Label']

## Distribution of classes (10-sec samples of each activity)
train_act_y.value_counts()

Label
Talk    160
Work     64
Name: count, dtype: int64

In this application we might start with a pattern discovery approach to help us identify useful predictive features. The code below graphs the raw time-series for Subject #1 (the training data) color-coding the two activities we're trying to predict:

## Create plots of select variables over time and color by activity (label)
my_vars = ['X_accel', 'Y_accel', 'Z_accel']
for cur_var in my_vars:
    plt.figure().set_figheight(2)
    plt.scatter(x = range(len(train_act_X[cur_var])), 
                y=train_act_X[cur_var], 
                c = act_ts[act_ts['Subject'] == 1]['Label'].map({'Work': 1, 'Talk': 0}))
    plt.title(cur_var)
    plt.show()

When working with time-series data, we're generally looking to engineer features that capture one of the following characteristics:

1. The typical values of the signal (mean, median, percentiles, etc.)
2. The amount of variation in the signal (standard deviation, range, etc.)
3. The shape or behavior of the signal over time (number of peaks, rates of change, etc.)

For these data, we can clearly observe:

1. Working tends to have higher average values, particularly along the X and Z axes.
2. Working tends to have lower amounts of variability in the accelerometer's measurements (aside from some unusual behavior at the very beginning of the series)

Representing these patterns in a tabular dataset can be done using grouped summarization (first introduced in Lab 1!):

## Get the mean and std dev of each selected variable for each sample
train_X_FE = pd.concat([train_act_X.groupby('Sample_Num')[my_vars].mean(),
                        train_act_X.groupby('Sample_Num')[my_vars].std()], axis=1)
train_X_FE.columns = ['Mean_X', 'Mean_Y', 'Mean_Z', 'SD_X', 'SD_Y', 'SD_Z']
train_X_FE.head(5)

	Mean_X	Mean_Y	Mean_Z	SD_X	SD_Y	SD_Z
Sample_Num
1	1897.030769	2292.657692	2064.600000	170.098461	143.059786	125.358252
2	1957.184615	2379.226923	2108.601923	11.067517	7.863383	16.358731
3	1962.921154	2378.396154	2112.584615	5.970345	3.332676	5.728319
4	1969.034615	2375.501923	2120.165385	14.418689	7.394212	15.416865
5	1976.882692	2373.157692	2130.376923	7.419165	6.416962	7.635497

You should notice that we grouped by Sample_Num here, which is because all of our data are from a single subject and we want to make predictions for 10-second samples. More generally, we should group by the variables needed to uniquely define an individual sample, which might be a combination of a sample ID and subject ID in other applications.

Features relating to the shape or behavior of the signal over time are often calculated using rolling windows. The example below creates a feature representing the maximum standard deviation occurring over any 1-second period in a sample:

## Function
def max_one_sec_sd(x):
    return x.rolling(52).std().max()

## Applying it to each sample (for each input channel)
train_act_X.groupby('Sample_Num')[my_vars].apply(max_one_sec_sd)

	X_accel	Y_accel	Z_accel
Sample_Num
1	182.212334	111.997069	181.242432
2	19.522225	17.247467	32.027250
3	7.132764	4.726294	8.229458
4	17.936582	11.411192	21.629717
5	11.748159	11.781364	10.603651
...	...	...	...
221	26.512327	9.133747	34.016077
222	22.813320	10.985439	19.616256
223	33.497456	36.880509	26.136491
224	33.187024	19.268524	28.775408
225	32.478279	16.053173	30.759124

224 rows × 3 columns

Question #2:

• Part A: Engineer an additional feature of your choosing that is created to capture the shape or behavior of the signal over time. You may choose to look at a rolling measure, one that captures peaks or valleys, or one that captures any other type of temporal trend.
• Part B: Assemble a data frame that includes the features created in this section's examples as well as the one you created in Part A.
• Part C: Fit a decision tree model with a maximum depth of 4 to the training data and display the fitted model using plot_tree().
• Part D: Apply the fitted model from Part C to the test set (the data from Subject 7) using classification accuracy as the scoring metric. Briefly comment on the performance of the model. Note: You might plot the Subject 7 data to better understand the model's performance.

Part 3 - Random Forests

Recall that a random forest is an ensemble model that combines predictions from many decision trees. By training each tree on a different bootstrapped sample of the data and randomly restricting which features can be considered at each split, random forests often achieve better predictive performance and generalize better than a single decision tree.

A unique aspect of random forests are that they naturally exclude some observations from the training set used by each tree. Bootstrapping will include multiple of some observations, while others will not be sampled. The data-points that were excluded during bootstrapping can be used to calculate "out-of-bag performance" (sometimes abbreviated OOB), which can be interpreted similarly to cross-validated measures of performance.

The existence of OOB performance estimates, combined with the fact that adding more trees will not harm performance, alters how we should approach hyperparameter tuning. In general our workflow will be:

1. Track OOB performance by n_estimators using default values for all other hyperparameters - the goal is to determine how many trees appear to be enough for performance to stabilize. While there's no guarantee that this will be optimal for all other hyperparameter configurations, it will most likely provide a reasonable value.
2. Tune the other main hyperparameters: max_features, max_depth, min_samples_leaf using a cross-validated grid search search (or random search) using the n_estimators chosen in step 1.
3. Plot the OOB learning curve for the final configuration of max_features, max_depth, min_samples_leaf to confirm enough trees were used.

Shown below is a template that plots the OOB performance for n_estimators going from 20 to 600 by increments of 20. Note that the argument warm_start=True greatly improves computational efficiency, as new trees are added to the existing forest (rather than recreating the entire forest from scratch each time we call rf.fit()).

from sklearn.ensemble import RandomForestClassifier
rf = RandomForestClassifier(n_estimators=1, warm_start=True, oob_score=True, random_state=0)

n_trees = []
oob_acc = []

## Increment the number of base learners, tracking OOB accuracy 
for n in range(20, 600, 20):
    rf.n_estimators = n
    rf.fit(train_X_FE, train_act_y)

    n_trees.append(n)
    oob_acc.append(rf.oob_score_)

plt.plot(n_trees, oob_acc, marker='o')
plt.xlabel("Number of Trees")
plt.ylabel("OOB Accuracy")
plt.title("Random Forest OOB Learning Curve")
plt.show()

Question #3: For this question you should use the data frame you created in Question 1 Part C, which contains 4 different predictive features.

• Part A: Using the default values for all hyperparameters except oob_score=True, create an OOB learning curve for a random forest classifier with n_estimators ranging from 20 to 400 by increments of 20. Based upon your learning curve, approximately how many trees appear necessary for OOB performance to stabilize?
• Part B: Using your answer from Part A for n_estimators, conduct a cross-validated grid search to tune each of the following hyperparameters:
  • max_depth - either 3 or None.
  • max_features - either 1 or 2
  • min_samples_leaf - either 1 or 50
• Part C: Create an OOB learning curve for the best random forest configuration from Part B. Does your initially chosen value of n_estimators (from Part A) appear to be sufficient, or is there evidence that additional trees might improve performance?
• Part D: Decide upon a final random forest model and evaluate it on the test data from Question 1 (reporting the classification accuracy). How does this performance compare to the decision tree model from Question 1?
• Part E: Create and display a confusion matrix showing the classification errors made by your final model from Part D on the test data. Briefly comment on how well the model performed within each class (clickbait and non-clickbait). If you decide that further feature engineering is needed, how might this influence the types of features you engineer?