In [1]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import mglearn
from IPython.display import display
%matplotlib inline

Chapter 4. Representing Data and Engineering Features

So far, we’ve assumed that our data comes in as a two-dimensional array of floating-point numbers, where each column is a continuous feature that describes the data points.
For many applications, this is not how the data is collected.
A particularly common type of feature is the categorical features.
Also known as discrete features, these are usually not numeric.
The distinction between categorical features and continuous features is analogous to the distinction between classification and regression, only on the input side rather than the output side.
Examples of continuous features that we have seen are pixel brightnesses and size measurements of plant flowers.
Examples of categorical features are the brand of a product, the color of a product, or the department (books, clothing, hardware) it is sold in.
These are all properties that can describe a product, but they don’t vary in a continuous way.
A product belongs either in the clothing department or in the books department.
There is no middle ground between books and clothing, and no natural order for the different categories (books is not greater or less than clothing, hardware is not between books and clothing, etc.).

Regardless of the types of features your data consists of, how you represent them can have an enormous effect on the performance of machine learning models.
We saw in Chapters 2 and 3 that scaling of the data is important.
In other words, if you don’t rescale your data (say, to unit variance), then it makes a difference whether you represent a measurement in centimeters or inches.
We also saw in Chapter 2 that it can be helpful to augment your data with additional features, like adding interactions (products) of features or more general polynomials.

The question of how to represent your data best for a particular application is known as feature engineering, and it is one of the main tasks of data scientists and machine learning practitioners trying to solve real-world problems.
Representing your data in the right way can have a bigger influence on the performance of a supervised model than the exact parameters you choose.
In this chapter, we will first go over the important and very common case of categorical features, and then give some examples of helpful transformations for specific combinations of features and models.

Categorical Variables

As an example, we will use the dataset of adult incomes in the United States, derived from the 1994 census database.
The task of the adult dataset is to predict whether a worker has an income of over \$50,000 or under \$50,000.
The features in this dataset include the workers’ ages, how they are employed (self employed, private industry employee, government employee, etc.), their education, their gender, their working hours per week, occupation, and more.
Table 4-1 shows the first few entries in the dataset:

Table 4-1. The first few entries in the adult employment dataset

The task is phrased as a classification task with the two classes being income <=50k and >50k.
It would also be possible to predict the exact income, and make this a regression task.
However, that would be much more difficult, and the 50K division is interesting to understand on its own.
In this dataset, age and hours-per-week are continuous features, which we know how to treat.
The workclass, education, sex, and occupation features are categorical, however.
All of them come from a fixed list of possible values, as opposed to a range, and denote a qualitative property, as opposed to a quantity.
As a starting point, let’s say we want to learn a logistic regression classifier on this data.
We know from Chapter 2 that a logistic regression makes predictions, $ŷ$, using the following formula:

$ŷ = w[0] * x[0] +w[1] * x[1] + \ldots + w[p] * x[p] + b > 0$

where $w[i]$ and $b$ are coefficients learned from the training set and $x[i]$ are the input features.
This formula makes sense when $x[i]$ are numbers, but not when $x[2]$ is "Masters" or "Bachelors".
Clearly we need to represent our data in some different way when applying logistic regression.
The next section will explain how we can overcome this problem.

One-Hot-Encoding (Dummy Variables)

By far the most common way to represent categorical variables is using the one-hot-encoding or one-out-of-N encoding, also known as dummy variables.
The idea behind dummy variables is to replace a categorical variable with one or more new features that can have the values 0 and 1.
The values 0 and 1 make sense in the formula for linear binary classification (and for all other models in scikit-learn), and we can represent any number of categories by introducing one new feature per category, as described here.
Let’s say for the workclass feature we have possible values of "Government Employee", "Private Employee", "Self Employed", and "Self Employed Incorporated".
To encode these four possible values, we create four new features, called "Government Employee", "Private Employee", "Self Employed", and "Self Employed Incorporated".
A feature is 1 if workclass for this person has the corresponding value and 0 otherwise, so exactly one of the four new features will be 1 for each data point.
This is why this is called one-hot or one-out-of-N encoding.

The principle is illustrated in Table 4-2.
A single feature is encoded using four new features.
When using this data in a machine learning algorithm, we would drop the original workclass feature and only keep the 0-1 features.

Table 4-2. Encoding the workclass feature using one-hot encoding

The one-hot encoding we use is quite similar, but not identical, to the dummy coding used in statistics.
For simplicity, we encode each category with a different binary feature.
In statistics, it is common to encode a categorical feature with k different possible values into k–1 features (the last one is represented as all zeros).
This is done to simplify the analysis (more technically, this will avoid making the data matrix rank-deficient.

There are two ways to convert your data to a one-hot encoding of categorical variables, using either pandas or scikit-learn.
At the time of writing, using pandas is slightly easier, so let’s go this route.
First we load the data using pandas from a comma-separated values (CSV) file:

In [2]:
import os
# The file has no headers naking the columns, so we pass 
# header=None and provide the column names explicitly in "names".
adult_path = os.path.join(mglearn.datasets.DATA_PATH, "")


In [3]:
data = pd.read_csv(adult_path, header=None, index_col=False,
       names=['age', 'workclass', 'fnlwgt', 'education',
              'education-num', 'marital-status', 'occupation',
              'relationship', 'race', 'gender', 'capital-gain',
              'capital-loss', 'hours-per-week', 'native-country',
# For illustrative purposes, we'll only select some of the columns:
data = data[['age', 'workclass', 'education', 'gender',
             'hours-per-week', 'occupation', 'income']]
# IPython.display allows nice output formatting within the 
# Jupyter notebook:

age workclass education gender hours-per-week occupation income
0 39 State-gov Bachelors Male 40 Adm-clerical <=50K
1 50 Self-emp-not-inc Bachelors Male 13 Exec-managerial <=50K
2 38 Private HS-grad Male 40 Handlers-cleaners <=50K
3 53 Private 11th Male 40 Handlers-cleaners <=50K
4 28 Private Bachelors Female 40 Prof-specialty <=50K
Table 4-3. The first five rows of the adult employment dataset

Checking string-encoded categorical data

After reading a dataset like this, it is often good to first check if a column actually contains meaningful categorical data.
When working with data that was input by humans (say, users on a website), there might not be a fixed set of categories, and differences in spelling and capitalization might require preprocessing.
For example, it might be that some people specified gender as “male” and some as “man,” and we might want to represent these two inputs using the same category.
A good way to check the contents of a column is using the value_counts method of a pandas Series (the type of a single column in a DataFrame), to show us what the unique values are and how often they appear:

In [4]:

 Male      21790
 Female    10771
Name: gender, dtype: int64

We can see that there are exactly two values for gender in this dataset, Male and Female, meaning the data is already in a good format to be represented using one-hot-encoding.
In a real application, you should look at all columns and check their values.
We will skip this here for brevity’s sake.
There is a very simple way to encode the data in pandas, using the get_dummies() function.
The get_dummies() function automatically transforms all columns that have object type (like strings) or are categorical (which is a special pandas concept that we haven’t talked about yet):

In [5]:
print("Original Features: \n", list(data.columns))

Original Features: 
 ['age', 'workclass', 'education', 'gender', 'hours-per-week', 'occupation', 'income']

In [6]:
data_dummies = pd.get_dummies(data)
print("Features After Applying get_dummies(): \n",

Features After Applying get_dummies(): 
 ['age', 'hours-per-week', 'workclass_ ?', 'workclass_ Federal-gov', 'workclass_ Local-gov', 'workclass_ Never-worked', 'workclass_ Private', 'workclass_ Self-emp-inc', 'workclass_ Self-emp-not-inc', 'workclass_ State-gov', 'workclass_ Without-pay', 'education_ 10th', 'education_ 11th', 'education_ 12th', 'education_ 1st-4th', 'education_ 5th-6th', 'education_ 7th-8th', 'education_ 9th', 'education_ Assoc-acdm', 'education_ Assoc-voc', 'education_ Bachelors', 'education_ Doctorate', 'education_ HS-grad', 'education_ Masters', 'education_ Preschool', 'education_ Prof-school', 'education_ Some-college', 'gender_ Female', 'gender_ Male', 'occupation_ ?', 'occupation_ Adm-clerical', 'occupation_ Armed-Forces', 'occupation_ Craft-repair', 'occupation_ Exec-managerial', 'occupation_ Farming-fishing', 'occupation_ Handlers-cleaners', 'occupation_ Machine-op-inspct', 'occupation_ Other-service', 'occupation_ Priv-house-serv', 'occupation_ Prof-specialty', 'occupation_ Protective-serv', 'occupation_ Sales', 'occupation_ Tech-support', 'occupation_ Transport-moving', 'income_ <=50K', 'income_ >50K']

You can see that the continuous features age and hours-per-week were not touched, while the categorical features were expanded into one new feature for each possible value:

In [7]:

age hours-per-week workclass_ ? workclass_ Federal-gov workclass_ Local-gov workclass_ Never-worked workclass_ Private workclass_ Self-emp-inc workclass_ Self-emp-not-inc workclass_ State-gov ... occupation_ Machine-op-inspct occupation_ Other-service occupation_ Priv-house-serv occupation_ Prof-specialty occupation_ Protective-serv occupation_ Sales occupation_ Tech-support occupation_ Transport-moving income_ <=50K income_ >50K
0 39 40 0 0 0 0 0 0 0 1 ... 0 0 0 0 0 0 0 0 1 0
1 50 13 0 0 0 0 0 0 1 0 ... 0 0 0 0 0 0 0 0 1 0
2 38 40 0 0 0 0 1 0 0 0 ... 0 0 0 0 0 0 0 0 1 0
3 53 40 0 0 0 0 1 0 0 0 ... 0 0 0 0 0 0 0 0 1 0
4 28 40 0 0 0 0 1 0 0 0 ... 0 0 0 1 0 0 0 0 1 0

5 rows × 46 columns

We can now use the values attribute to convert the data_dummies DataFrame into a NumPy array, and then train a machine learning model on it.
Be careful to separate the target variable (which is now encoded in two income columns) from the data before training a model.
Including the output variable, or some derived property of the output variable, into the feature representation is a very common mistake in building supervised machine learning models.

Be careful: column indexing in pandas includes the end of the range, so 'age':'occupation_ Transport-moving' is inclusive of occupation_ Transport-moving.
This is different from slicing a NumPy array, where the end of a range is not included: for example, np.arange(11)[0:10] doesn’t include the entry with index 10.

In this case, we extract only the columns containing features -- that is, all columns from age to occupation_ Transport-moving.
This range contains all of the features but not the target:

In [8]:
features = data_dummies.loc[:, 'age':'occupation_ Transport-moving']
# Extract NumPy arrays:
X = features.values
print("X[:5]: \n", X[:5], "\n")
y = data_dummies['income_ >50K'].values
print("y[:5]: \n", y[:5], "\n")
print("X.shape: {}  y.shape: {}".format(X.shape, y.shape))

 [[39 40  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  1  0  0  0
   0  0  0  0  1  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0]
 [50 13  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0
   0  0  0  0  1  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0]
 [38 40  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0
   0  0  0  0  1  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0]
 [53 40  0  0  0  0  1  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0
   0  0  0  0  1  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0]
 [28 40  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0
   0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0]] 

 [0 0 0 0 0] 

X.shape: (32561, 44)  y.shape: (32561,)

Now the data is represented in a way that scikit-learn can work with, and we can proceed as usual:

In [9]:
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y,
logreg = LogisticRegression(), y_train)
print("Test Score: {:.2f}".format(logreg.score(X_test, y_test)))

Test Score: 0.81

In this example, we called get_dummies() on a DataFrame containing both the training and the test data.
This is important to ensure categorical values are represented in the same way in the training set and the test set.
Imagine we have the training and test sets in two different DataFrames.
If the "Private Employee" value for the workclass feature does not appear in the test set, pandas will assume there are only three possible values for this feature and will create only three new dummy features.
Now our training and test sets have different numbers of features, and we can’t apply the model we learned on the training set to the test set anymore.
Even worse, imagine the workclass feature has the values "Government Employee" and "Private Employee" in the training set, and "Self Employed" and "Self Employed Incorporated" in the test set.
In both cases, pandas will create two new dummy features, so the encoded DataFrames will have the same number of features.
However, the two dummy features have entirely different meanings in the training and test sets.
The column that means "Government Employee" for the training set would encode "Self Employed" for the test set.
If we built a machine learning model on this data it would work very badly, because it would assume the columns mean the same things (because they are in the same position) when in fact they mean very different things.
To fix this, either call get_dummies() on a DataFrame that contains both the training and the test data points, or make sure that the column names are the same for the training and test sets after calling get_dummies(), to ensure they have the same semantics.

Numbers Can Encode Categoricals

In the example of the adult dataset, the categorical variables were encoded as strings.
On the one hand, that opens up the possibility of spelling errors, but on the other hand, it clearly marks a variable as categorical.
Often, whether for ease of storage or because of the way the data is collected, categorical variables are encoded as integers.
For example, imagine the census data in the adult dataset was collected using a questionnaire, and the answers for workclass were recorded as 0 (first box ticked), 1 (second box ticked), 2 (third box ticked), and so on.
Now the column will contain numbers from 0 to 8, instead of strings like "Private", and it won’t be immediately obvious to someone looking at the table representing the dataset whether they should treat this variable as continuous or categorical.
Knowing that the numbers indicate employment status, however, it is clear that these are very distinct states and should not be modeled by a single continuous variable.

Categorical features are often encoded using integers.
That they are numbers doesn’t mean that they should necessarily be treated as continuous features.
It is not always clear whether an integer feature should be treated as continuous or discrete (and one-hot-encoded).
If there is no ordering between the semantics that are encoded (like in the workclass example), the feature must be treated as discrete.
For other cases, like five-star ratings, the better encoding depends on the particular task and data and which machine learning algorithm is used.

The get_dummies() function in pandas treats all numbers as continuous and will not create dummy variables for them.
To get around this, you can either use scikit-learn’s OneHotEncoder, for which you can specify which variables are continuous and which are discrete, or convert numeric columns in the DataFrame to strings.
To illustrate, let’s create a DataFrame object with two columns, one containing strings and one containing integers:

In [10]:
# Create a DataFrame with an integer feature and a categorical
# string feature:
demo_df = pd.DataFrame({'Integer Feature': [0, 1, 2, 1],
                        'Categorical Feature': ['socks',

Categorical Feature Integer Feature
0 socks 0
1 fox 1
2 socks 2
3 box 1
Table 4-4. DataFrame containig categorical string features and integer features

Using get_dummies() will only encode the string feature and will not change the integer feature, as you can see in Table 4-5:

In [11]:

Integer Feature Categorical Feature_box Categorical Feature_fox Categorical Feature_socks
0 0 0 0 1
1 1 0 1 0
2 2 0 0 1
3 1 1 0 0
Table 4-5. One-hot-encoded version of the data from Table 4-4, leaving the integer feature unchanged

If you want dummy variables to be created for the "Integer Feature" column, you can explicitly list the columns you want to encode using the columns parameter.
Then, both features will be treated as categorical (see Table 4-6):

In [12]:
demo_df['Integer Feature'] = demo_df['Integer Feature'].astype(str)
                       columns=['Integer Feature',
                                'Categorical Feature']))

Integer Feature_0 Integer Feature_1 Integer Feature_2 Categorical Feature_box Categorical Feature_fox Categorical Feature_socks
0 1 0 0 0 0 1
1 0 1 0 0 1 0
2 0 0 1 0 0 1
3 0 1 0 1 0 0
Table 4-6. One-hot-encoding of the data shown in Table 4-4, encoding the integer and string features.

Binning, Discretization, Linear Models, and Trees

The best way to represent data depends not only on the semantics of the data, but also on the kind of model you are using.
Linear models and tree-based models (such as decision trees, gradient boosted trees, and random forests), two large and very commonly used families, have very different properties when it comes to how they work with different feature representations.
Let’s go back to the wave regression dataset that we used in Chapter 2.
It has only a single input feature.
Here is a comparison of a linear regression model and a decision tree regressor on this dataset (see Figure 4-1):

In [13]:
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor

X, y = mglearn.datasets.make_wave(n_samples=100)
line = np.linspace(-3, 3, 1000, endpoint=False).reshape(-1, 1)

reg = DecisionTreeRegressor(min_samples_split=3).fit(X, y)
plt.plot(line, reg.predict(line), label='decision tree')

reg = LinearRegression().fit(X, y)
plt.plot(line, reg.predict(line), label='linear regression')

plt.plot(X[:, 0], y, 'o', c='k')
plt.ylabel("Regression Output")
plt.xlabel("Input Feature")

/Users/benjamingrove/.pyenv/versions/3.6.1/lib/python3.6/site-packages/scipy/linalg/ RuntimeWarning: internal gelsd driver lwork query error, required iwork dimension not returned. This is likely the result of LAPACK bug 0038, fixed in LAPACK 3.2.2 (released July 21, 2010). Falling back to 'gelss' driver.
  warnings.warn(mesg, RuntimeWarning)
<matplotlib.legend.Legend at 0x10b760080>
Figure 4-1. Comparing linear regression and a decision tree on the wave dataset

As you know, linear models can only model linear relationships, which are lines in the case of a single feature.
The decision tree can build a much more complex model of the data.
However, this is strongly dependent on the representation of the data.
One way to make linear models more powerful on continuous data is to use binning (also known as discretization) of the feature to split it up into multiple features, as described here.

We imagine a partition of the input range for the feature (in this case, the numbers from –3 to 3) into a fixed number of bins -- say, 10.
A data point will then be represented by which bin it falls into.
To determine this, we first have to define the bins.
In this case, we’ll define 10 bins equally spaced between –3 and 3.
We use the np.linspace function for this, creating 11 entries, which will create 10 bins -- they are the spaces in between two consecutive boundaries:

In [14]:
bins = np.linspace(-3, 3, 11)
print("bins: {}".format(bins))

bins: [-3.  -2.4 -1.8 -1.2 -0.6  0.   0.6  1.2  1.8  2.4  3. ]

Here, the first bin contains all data points with feature values –3 to –2.4, the second bin contains all points with feature values from –2.4 to –1.8, and so on.
Next, we record for each data point which bin it falls into.
This can be easily computed using the np.digitize() function:

In [15]:
which_bin = np.digitize(X, bins=bins)
print("\nData Points:\n", X[:5])
print("\nBin Membership for Data Points:\n", which_bin[:5])

Data Points:
 [ 2.70428584]
 [ 1.39196365]
 [ 0.59195091]

Bin Membership for Data Points:
 [[ 4]
 [ 8]
 [ 6]
 [ 2]]

What we did here is transform the single continuous input feature in the wave dataset into a categorical feature that encodes which bin a data point is in.
To use a scikit-learn model on this data, we transform this discrete feature to a one-hot encoding using the OneHotEncoder from the preprocessing module.
The OneHotEncoder does the same encoding as pandas.get_dummies(), though it currently only works on categorical variables that are integers:

In [16]:
from sklearn.preprocessing import OneHotEncoder

# Transform data using the OneHotEncoder:
encoder = OneHotEncoder(sparse=False)
# Use to find the unique values that appear
# in which_bin:
# Transform creates the one-hot encoding:
X_binned = encoder.transform(which_bin)

[[ 0.  0.  0.  1.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  1.]
 [ 0.  0.  0.  0.  0.  0.  0.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  1.  0.  0.  0.  0.]
 [ 0.  1.  0.  0.  0.  0.  0.  0.  0.  0.]]

Because we specified 10 bins, the transformed dataset X_binned now is made up of 10 features:

In [17]:
print("X_binned.shape: {}".format(X_binned.shape))

X_binned.shape: (100, 10)

Now we can build a new linear regression model and a new decision tree model on the one-hot-encoded data.
The result is visualized in Figure 4-2, together with the bin boundaries, shown as dotted black lines:

In [18]:
line_binned = encoder.transform(np.digitize(line, bins=bins))

reg = LinearRegression().fit(X_binned, y)
plt.plot(line, reg.predict(line_binned),
         label="Linear Regression Binned")

reg = DecisionTreeRegressor(min_samples_split=3).fit(X_binned, y)
plt.plot(line, reg.predict(line_binned),
         label="Decision Tree Binned")

plt.plot(X[:, 0], y, 'o', c='k')
plt.vlines(bins, -3, 3, linewidth=1, alpha=.2)
plt.ylabel("Regression Output")
plt.xlabel("Input Feature")

<matplotlib.text.Text at 0x10c3bb748>
Figure 4-2. Comparing linear regression and decision tree regression on binned features. The lines for linear regression and decision tree regression are exactly on top of one another.

The dashed line and solid line are exactly on top of each other, meaning the linear regression model and the decision tree make exactly the same predictions.
For each bin, they predict a constant value.
As features are constant within each bin, any model must predict the same value for all points within a bin.
Comparing what the models learned before binning the features and after, we see that the linear model became much more flexible, because it now has a different value for each bin, while the decision tree model got much less flexible.

Binning features generally has no beneficial effect for tree-based models, as these models can learn to split up the data anywhere.
In a sense, that means decision trees can learn whatever binning is most useful for predicting on this data.
Additionally, decision trees look at multiple features at once, while binning is usually done on a per-feature basis.
However, the linear model benefited greatly in expressiveness from the transformation of the data.
If there are good reasons to use a linear model for a particular dataset -- say, because it is very large and high-dimensional, but some features have nonlinear relations with the output -- binning can be a great way to increase modeling power.

Interactions and Polynomials

Another way to enrich a feature representation, particularly for linear models, is adding interaction features and polynomial features of the original data.
This kind of feature engineering is often used in statistical modeling, but it’s also common in many practical machine learning applications.

As a first example, look again at Figure 4-2.
The linear model learned a constant value for each bin in the wave dataset.
We know, however, that linear models can learn not only offsets, but also slopes.
One way to add a slope to the linear model on the binned data is to add the original feature (the x-axis in the plot) back in.
This leads to an 11-dimensional dataset, as seen in Figure 4-3:

In [19]:
X_combined = np.hstack([X, X_binned])

(100, 11)

In [20]:
reg = LinearRegression().fit(X_combined, y)

line_combined = np.hstack([line, line_binned])
plt.plot(line, reg.predict(line_combined),
         label="Linear Regression Combined")

for bin in bins:
    plt.plot([bin, bin], [-3, 3], ':', c='k', linewidth=1)
plt.ylabel("Regression Output")
plt.xlabel("Input Feature")
plt.plot(X[:, 0], y, 'o', c='k')

[<matplotlib.lines.Line2D at 0x11204fc88>]
Figure 4-3. Linear Regression using binned features and a single global scope

In this example, the model learned an offset for each bin, together with a slope.
The learned slope is downward, and shared across all the bins --there is a single x-axis feature, which has a single slope.
Because the slope is shared across all bins, it doesn’t seem to be very helpful.
We would rather have a separate slope for each bin!
We can achieve this by adding an interaction or product feature that indicates which bin a data point is in and where it lies on the x-axis.
This feature is a product of the bin indicator and the original feature.
Let’s create this dataset:

In [21]:
X_product = np.hstack([X_binned, X * X_binned])

(100, 20)

The dataset now has 20 features: the indicators for which bin a data point is in, and a product of the original feature and the bin indicator.
You can think of the product feature as a separate copy of the x-axis feature for each bin.
It is the original feature within the bin, and zero everywhere else.
Figure 4-4 shows the result of the linear model on this new representation:

In [22]:
reg = LinearRegression().fit(X_product, y)

line_product = np.hstack([line_binned, line * line_binned])
plt.plot(line, reg.predict(line_product),
         label="Linear Regression Product")

for bin in bins:
    plt.plot([bin, bin], [-3, 3], ':', c='k', linewidth=1)
plt.plot(X[:, 0], y, 'o', c='k')
plt.ylabel("Regression Output")
plt.xlabel("Input Features")

<matplotlib.legend.Legend at 0x11208c5c0>
Figure 4-4. Linear regression with a separate slope per bin

As you can see, now each bin has its own offset and slope in this model.
Using binning is one way to expand a continuous feature.
Another one is to use polynomials of the original features.
For a given feature x, we might want to consider x ** 2, x ** 3, x ** 4, and so on.
This is implemented in PolynomialFeatures in the preprocessing module:

In [23]:
from sklearn.preprocessing import PolynomialFeatures
# Include polynomials up to x ** 10.
# The default "include_bias=True" adds a feature that's
# constantly 1.
poly = PolynomialFeatures(degree=10, include_bias=False)
X_poly = poly.transform(X)

Using a degree of 10 yields 10 features:

In [24]:
print("X_poly.shape: {}".format(X_poly.shape))

X_poly.shape: (100, 10)

Let's compare the entries of X_poly to those of X:

In [25]:
print("Entries of X: \n{}".format(X[:5]))

Entries of X: 
 [ 2.70428584]
 [ 1.39196365]
 [ 0.59195091]

In [26]:
print("Entries of X_poly: \n{}".format(X_poly[:5]))

Entries of X_poly: 
[[ -7.52759287e-01   5.66646544e-01  -4.26548448e-01   3.21088306e-01
   -2.41702204e-01   1.81943579e-01  -1.36959719e-01   1.03097700e-01
   -7.76077513e-02   5.84199555e-02]
 [  2.70428584e+00   7.31316190e+00   1.97768801e+01   5.34823369e+01
    1.44631526e+02   3.91124988e+02   1.05771377e+03   2.86036036e+03
    7.73523202e+03   2.09182784e+04]
 [  1.39196365e+00   1.93756281e+00   2.69701700e+00   3.75414962e+00
    5.22563982e+00   7.27390068e+00   1.01250053e+01   1.40936394e+01
    1.96178338e+01   2.73073115e+01]
 [  5.91950905e-01   3.50405874e-01   2.07423074e-01   1.22784277e-01
    7.26822637e-02   4.30243318e-02   2.54682921e-02   1.50759786e-02
    8.92423917e-03   5.28271146e-03]
 [ -2.06388816e+00   4.25963433e+00  -8.79140884e+00   1.81444846e+01
   -3.74481869e+01   7.72888694e+01  -1.59515582e+02   3.29222321e+02
   -6.79478050e+02   1.40236670e+03]]

You can obtain the semantics of the features by calling the get_feature_names() method, which provides the exponent for each feature:

In [27]:
print("Polynomial feature names: \n{}".format(

Polynomial feature names: 
['x0', 'x0^2', 'x0^3', 'x0^4', 'x0^5', 'x0^6', 'x0^7', 'x0^8', 'x0^9', 'x0^10']

You can see that the first column of X_poly corresponds exactly to X, while the other columns are the powers of the first entry.
It’s interesting to see how large some of the values can get.
The second row has entries above 20,000, orders of magnitude different from the rest.
Using polynomial features together with a linear regression model yields the classical model of polynomial regression (see Figure 4-5):

In [28]:
reg = LinearRegression().fit(X_poly, y)
line_poly = poly.transform(line)

plt.plot(line, reg.predict(line_poly),
         label="Polynomial Linear Regression")
plt.plot(X[:, 0], y, 'o', c='k')
plt.ylabel("Regression Output")
plt.xlabel("Input Feature")

<matplotlib.legend.Legend at 0x114960d68>
Figure 4-5. Linear regression with tenth-degree polynomial features

As you can see, polynomial features yield a very smooth fit on this one-dimensional data.
However, polynomials of high degree tend to behave in extreme ways on the boundaries or in regions with little data.
As a comparison, here is a kernel SVM model learned on the original data, without any transformation (see Figure 4-6):

In [29]:
from sklearn.svm import SVR

for gamma in [1, 10]:
    svr = SVR(gamma=gamma).fit(X, y)
    plt.plot(line, svr.predict(line),
             label="SVR gamma={}".format(gamma))
plt.plot(X[:, 0], y, 'o', c='k')
plt.ylabel("Regression Output")
plt.xlabel("Input Feature")

<matplotlib.legend.Legend at 0x11196b668>
Figure 4-6. Comparison of different gamma parameters for an SVM with RBF kernel

Using a more complex model, a kernel SVM, we are able to learn a similarly complex prediction to the polynomial regression without an explicit transformation of the features.

As a more realistic application of interactions and polynomials, let’s look again at the Boston Housing dataset.
We already used polynomial features on this dataset in Chapter 2.
Now let’s have a look at how these features were constructed, and at how much the polynomial features help.
First we load the data, and rescale it to between 0 and 1 using MinMaxScaler:

In [30]:
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import MinMaxScaler

boston = load_boston()
X_train, X_test, y_train, y_test = train_test_split(,, random_state=0)
# Rescale the data:
scaler = MinMaxScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

Now we extract polynomial features and interactions up to a degree of 2:

In [31]:
poly = PolynomialFeatures(degree=2).fit(X_train_scaled)
X_train_poly = poly.transform(X_train_scaled)
X_test_poly = poly.transform(X_test_scaled)
print("X_train.shape: {}".format(X_train.shape))
print("X_train_poly.shape: {}".format(X_train_poly.shape))

X_train.shape: (379, 13)
X_train_poly.shape: (379, 105)

The data originally had 13 features, which were expanded into 105 interaction features.
These new features represent all possible interactions between two different original features, as well as the square of each original feature.
degree=2 here means that we look at all features that are the product of up to two original features.
The exact correspondence between input and output features can be found using the get_feature_names() method:

In [32]:
print("Polynomial Feature Names: \n{}".format(

Polynomial Feature Names: 
['1', 'x0', 'x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x9', 'x10', 'x11', 'x12', 'x0^2', 'x0 x1', 'x0 x2', 'x0 x3', 'x0 x4', 'x0 x5', 'x0 x6', 'x0 x7', 'x0 x8', 'x0 x9', 'x0 x10', 'x0 x11', 'x0 x12', 'x1^2', 'x1 x2', 'x1 x3', 'x1 x4', 'x1 x5', 'x1 x6', 'x1 x7', 'x1 x8', 'x1 x9', 'x1 x10', 'x1 x11', 'x1 x12', 'x2^2', 'x2 x3', 'x2 x4', 'x2 x5', 'x2 x6', 'x2 x7', 'x2 x8', 'x2 x9', 'x2 x10', 'x2 x11', 'x2 x12', 'x3^2', 'x3 x4', 'x3 x5', 'x3 x6', 'x3 x7', 'x3 x8', 'x3 x9', 'x3 x10', 'x3 x11', 'x3 x12', 'x4^2', 'x4 x5', 'x4 x6', 'x4 x7', 'x4 x8', 'x4 x9', 'x4 x10', 'x4 x11', 'x4 x12', 'x5^2', 'x5 x6', 'x5 x7', 'x5 x8', 'x5 x9', 'x5 x10', 'x5 x11', 'x5 x12', 'x6^2', 'x6 x7', 'x6 x8', 'x6 x9', 'x6 x10', 'x6 x11', 'x6 x12', 'x7^2', 'x7 x8', 'x7 x9', 'x7 x10', 'x7 x11', 'x7 x12', 'x8^2', 'x8 x9', 'x8 x10', 'x8 x11', 'x8 x12', 'x9^2', 'x9 x10', 'x9 x11', 'x9 x12', 'x10^2', 'x10 x11', 'x10 x12', 'x11^2', 'x11 x12', 'x12^2']

The first new feature is a constant feature, called "1" here. The next 13 features are the original features (called "x0" to "x12").
Then follows the first feature squared ("x0^2") and combinations of the first and the other features.

Let’s compare the performance using Ridge on the data with and without interactions:

In [33]:
from sklearn.linear_model import Ridge

ridge = Ridge().fit(X_train_scaled, y_train)
print("Score without interactions: {:.3f}".format(
    ridge.score(X_test_scaled, y_test)))

ridge = Ridge().fit(X_train_poly, y_train)
print("Score with interactions: {:.3f}".format(
    ridge.score(X_test_poly, y_test)))

Score without interactions: 0.621
Score with interactions: 0.753

Clearly, the interactions and polynomial features gave us a good boost inperformance when using Ridge.
When using a more complex model like a random forest, the story is a bit different:

In [34]:
from sklearn.ensemble import RandomForestRegressor

rf = RandomForestRegressor(n_estimators=100).fit(
    X_train_scaled, y_train)
print("Score without interactions: {:.3f}".format(
    rf.score(X_test_scaled, y_test)))

rf = RandomForestRegressor(n_estimators=100).fit(
    X_train_poly, y_train)
print("Score with interactions: {:.3f}".format(
    rf.score(X_test_poly, y_test)))

Score without interactions: 0.805
Score with interactions: 0.775

You can see that even without additional features, the random forest beats the performance of Ridge.
Adding interactions and polynomials actually decreases performance slightly.

Univariate Nonlinear Transformations

We just saw that adding squared or cubed features can help linear models for regression.
There are other transformations that often prove useful for transforming certain features: in particular, applying mathematical functions like log(), exp(), or sin().
While tree-based models only care about the ordering of the features, linear models and neural networks are very tied to the scale and distribution of each feature, and if there is a nonlinear relation between the feature and the target, that becomes hard to model -- particularly in regression.
The functions log() and exp() can help by adjusting the relative scales in the data so that they can be captured better by a linear model or neural network.
We saw an application of that in Chapter 2 with the memory price data.
The sin() and cos() functions can come in handy when dealing with data that encodes periodic patterns.

Most models work best when each feature (and in regression also the target) is loosely Gaussian distributed -- that is, a histogram of each feature should have something resembling the familiar "bell curve" shape.
Using transformations like log() and exp() is a hacky but simple and efficient way to achieve this.
A particularly common case when such a transformation can be helpful is when dealing with integer count data.
By count data, we mean features like "how often did user A log in?"
Counts are never negative, and often follow particular statistical patterns.
We are using a synthetic dataset of counts here that has properties similar to those you can find in the wild.
The features are all integer-valued, while the response is continuous:

In [35]:
rnd = np.random.RandomState(0)
X_org = rnd.normal(size=(1000, 3))
w = rnd.normal(size=3)

# The Poisson distribution is a discrete probability distribution
# that expresses the probability of a given number of events 
# occurring in a fixed interval of time and/or space if these 
# events occur with a known constant rate and independently of
# the time since the last event.
X = rnd.poisson(10 * np.exp(X_org))
y =, w)

Let’s look at the first 10 entries of the first feature.
All are integer values and positive, but apart from that it’s hard to make out a particular pattern.
If we count the appearance of each value, the distribution of values becomes clearer:

In [36]:
print("Number of Feature Appearances: \n{}".format(
    np.bincount(X[:, 0])))

Number of Feature Appearances: 
[28 38 68 48 61 59 45 56 37 40 35 34 36 26 23 26 27 21 23 23 18 21 10  9 17
  9  7 14 12  7  3  8  4  5  5  3  4  2  4  1  1  3  2  5  3  8  2  5  2  1
  2  3  3  2  2  3  3  0  1  2  1  0  0  3  1  0  0  0  1  3  0  1  0  2  0
  1  1  0  0  0  0  1  0  0  2  2  0  1  1  0  0  0  0  1  1  0  0  0  0  0
  0  0  1  0  0  0  0  0  1  1  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0
  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1]

The value 2 seems to be the most common, with 68 appearances (bincount always starts at 0), and the counts for higher values fall quickly.
However, there are some very high values, like 84 and 85, that are appearing twice.
We visualize the counts in Figure 4-7:

In [37]:
bins = np.bincount(X[:, 0]), bins, color='r')
plt.ylabel("Number of Appearances")

<matplotlib.text.Text at 0x115bfa080>
Figure 4-7. Histogram of feature values for X[0]

Features X[:, 1] and X[:, 2] have similar properties.
This kind of distribution of values (many small ones and a few very large ones) is very common in practice.
However, it is something most linear models can’t handle very well.
Let’s try to fit a ridge regression to this model:

In [38]:
from sklearn.linear_model import Ridge

X_train, X_test, y_train, y_test = train_test_split(
                                X, y, random_state=0)
score = Ridge().fit(X_train, y_train).score(X_test, y_test)
print("Test score: {:.3f}".format(score))

Test score: 0.622

As you can see from the relatively low $R^2$ score, Ridge was not able to really capture the relationship between X and y.
Applying a logarithmic transformation can help, though.
Because the value 0 appears in the data (and the logarithm is not defined at 0), we can’t actually just apply log(), but we have to compute log(X + 1):

In [39]:
X_train_log = np.log(X_train + 1)
X_test_log = np.log(X_test + 1)
print("First entries of X_test: \n{}".format(X_test[:5]))
print("First entries of X_test_log: \n{}".format(X_test_log[:5]))

First entries of X_test: 
[[10  7 22]
 [16  5  4]
 [ 2 51  7]
 [11 24 14]
 [16 16  0]]

First entries of X_test_log: 
[[ 2.39789527  2.07944154  3.13549422]
 [ 2.83321334  1.79175947  1.60943791]
 [ 1.09861229  3.95124372  2.07944154]
 [ 2.48490665  3.21887582  2.7080502 ]
 [ 2.83321334  2.83321334  0.        ]]

After the transformation, the distibution of the data is less asymmetrical and doesn't have any large outliers anymore (see Figure 4-8 and Figure 4-8A):

In [40]:
plt.hist(X_train_log[:, 0], bins=25, color='gray')
plt.ylabel("Number of Appearances")

<matplotlib.text.Text at 0x115d377b8>
Figure 4-8. Histogram of feature values for X[0] after logarithmic transformation

Let's show a comparison with the un-transformed X_train data:

In [41]:
plt.hist(X_train[:, 0], bins=25, color='gray')
plt.ylabel("Number of Appearances")

<matplotlib.text.Text at 0x10b80bba8>
Figure 4-8A. Histogram of feature values for X[0] before logarithmic transformation

As you will see, building a ridge model on the new data provides a much better fit:

In [42]:
score = Ridge().fit(X_train_log, y_train).score(
                             X_test_log, y_test)
print("Test Score: {:.3f}".format(score))

Test Score: 0.875

Finding the transformation that works best for each combination of dataset and model is somewhat of an art.
In this example, all the features had the same properties.
This is rarely the case in practice, and usually only a subset of the features should be transformed, or sometimes each feature needs to be transformed in a different way.
As we mentioned earlier, these kinds of transformations are irrelevant for tree-based models but might be essential for linear models.
Sometimes it is also a good idea to transform the target variable y in regression.
Trying to predict counts (say, number of orders) is a fairly common task, and using the log(y + 1) transformation often helps.

As you saw in the previous examples, binning, polynomials, and interactions can have a huge influence on how models perform on a given dataset.
This is particularly true for less complex models like linear models and naive Bayes models.
Tree-based models, on the other hand, are often able to discover important interactions themselves, and don’t require transforming the data explicitly most of the time.
Other models, like SVMs, nearest neighbors, and neural networks, might sometimes benefit from using binning, interactions, or polynomials, but the implications there are usually much less clear than in the case of linear models.

Automatic Feature Selection

With so many ways to create new features, you might get tempted to increase the dimensionality of the data way beyond the number of original features.
However, adding more features makes all models more complex, and so increases the chance of overfitting.
When adding new features, or with high-dimensional datasets in general, it can be a good idea to reduce the number of features to only the most useful ones, and discard the rest.
This can lead to simpler models that generalize better.

But how can you know how good each feature is?
There are three basic strategies: univariate statistics, model-based selection, and iterative selection.
We will discuss all three of them in detail.
All of these methods are supervised methods, meaning they need the target for fitting the model.
This means we need to split the data into training and test sets, and fit the feature selection only on the training part of the data.

Univariate Statistics

In univariate statistics, we compute whether there is a statistically significant relationship between each feature and the target.
Then the features that are related with the highest confidence are selected.
In the case of classification, this is also known as analysis of variance (ANOVA).
A key property of these tests is that they are univariate, meaning that they only consider each feature individually.
Consequently, a feature will be discarded if it is only informative when combined with another feature.
Univariate tests are often very fast to compute, and don’t require building a model.
On the other hand, they are completely independent of the model that you might want to apply after the feature selection.

To use univariate feature selection in scikit-learn, you need to choose a test, usually either f_classif (the default) for classification or f_regression for regression, and a method to discard features based on the $p$-values determined in the test.
All methods for discarding parameters use a threshold to discard all features with too high a $p$-value (which means they are unlikely to be related to the target).
The methods differ in how they compute this threshold, with the simplest ones being SelectKBest, which selects a fixed number $k$ of features, and SelectPercentile, which selects a fixed percentage of features.
Let’s apply the feature selection for classification on the cancer dataset.
To make the task a bit harder, we’ll add some noninformative noise features to the data.
We expect the feature selection to be able to identify the features that are noninformative and remove them:

In [43]:
from sklearn.datasets import load_breast_cancer
from sklearn.feature_selection import SelectPercentile
from sklearn.model_selection import train_test_split

cancer = load_breast_cancer()

# Get deterministic random numbers
rng = np.random.RandomState(42)
noise = rng.normal(size=(len(, 50))
# Add noise features to the data.
# The first 30 features are from the dataset.
# The next 50 features are randomly generated noise.
X_w_noise = np.hstack([, noise])

array([  1.79900000e+01,   1.03800000e+01,   1.22800000e+02,
         1.00100000e+03,   1.18400000e-01,   2.77600000e-01,
         3.00100000e-01,   1.47100000e-01,   2.41900000e-01,
         7.87100000e-02,   1.09500000e+00,   9.05300000e-01,
         8.58900000e+00,   1.53400000e+02,   6.39900000e-03,
         4.90400000e-02,   5.37300000e-02,   1.58700000e-02,
         3.00300000e-02,   6.19300000e-03,   2.53800000e+01,
         1.73300000e+01,   1.84600000e+02,   2.01900000e+03,
         1.62200000e-01,   6.65600000e-01,   7.11900000e-01,
         2.65400000e-01,   4.60100000e-01,   1.18900000e-01,
         4.96714153e-01,  -1.38264301e-01,   6.47688538e-01,
         1.52302986e+00,  -2.34153375e-01,  -2.34136957e-01,
         1.57921282e+00,   7.67434729e-01,  -4.69474386e-01,
         5.42560044e-01,  -4.63417693e-01,  -4.65729754e-01,
         2.41962272e-01,  -1.91328024e+00,  -1.72491783e+00,
        -5.62287529e-01,  -1.01283112e+00,   3.14247333e-01,
        -9.08024076e-01,  -1.41230370e+00,   1.46564877e+00,
        -2.25776300e-01,   6.75282047e-02,  -1.42474819e+00,
        -5.44382725e-01,   1.10922590e-01,  -1.15099358e+00,
         3.75698018e-01,  -6.00638690e-01,  -2.91693750e-01,
        -6.01706612e-01,   1.85227818e+00,  -1.34972247e-02,
        -1.05771093e+00,   8.22544912e-01,  -1.22084365e+00,
         2.08863595e-01,  -1.95967012e+00,  -1.32818605e+00,
         1.96861236e-01,   7.38466580e-01,   1.71368281e-01,
        -1.15648282e-01,  -3.01103696e-01,  -1.47852199e+00,
        -7.19844208e-01,  -4.60638771e-01,   1.05712223e+00,
         3.43618290e-01,  -1.76304016e+00])

In [44]:
X_train, X_test, y_train, y_test = train_test_split(
    X_w_noise,, random_state=0, test_size=0.5)
# Use f_classif (the default) and SelectPercentile to
# select 50% of features:
select = SelectPercentile(percentile=50), y_train)
# Transform training set:
X_train_selected = select.transform(X_train)

print("X_train.shape: {}".format(X_train.shape))
print("X_train_selected.shape: {}".format(

X_train.shape: (284, 80)
X_train_selected.shape: (284, 40)

As you can see, the number of features was reduced from 80 to 40 (50 percent of the original features).
We can find out which features have been selected using the get_support() method, which returns a Boolean mask of the slected features (visualized in Figure 4-9):

In [45]:
mask = select.get_support()

[ True  True  True  True  True  True  True  True  True False  True False
  True  True  True  True  True  True False False  True  True  True  True
  True  True  True  True  True  True False False False  True False  True
 False False  True False False False False  True False False  True False
 False  True False  True False False False False False False  True False
  True False False False False  True False  True False False False False
  True  True False  True False False False False]

Nice, but visualizing the mask may be a bit more helpful.
Black is true, white is False:

In [46]:
plt.matshow(mask.reshape(1, -1), cmap="gray_r")
plt.xlabel("Sample Index")

([], <a list of 0 Text yticklabel objects>)
Figure 4-9. Features selected by SelectPercentile

As you can see from the visualization of the mask, most of the selected features are the original features, and most of the noise features were removed.
However, the recovery of the original features is not perfect.
Let’s compare the performance of logistic regression on all features against the performance using only the selected features:

In [47]:
from sklearn.linear_model import LogisticRegression

# Transform test data:
X_test_selected = select.transform(X_test)

lr = LogisticRegression(), y_train)
print("Score with all features: \n{:.3f}".format(
                        lr.score(X_test, y_test))), y_train)
print("Score with only the selected features: \n{:.3f}".format(
                            lr.score(X_test_selected, y_test)))

Score with all features: 
Score with only the selected features: 

In this case, removing the noise features improved performance, even though some of the original features were lost.
This was a very simple synthetic example, and outcomes on real data are usually mixed.
Univariate feature selection can still be very helpful, though, if there is such a large number of features that building a model on them is infeasible, or if you suspect that many features are completely uninformative.

Model-Based Feature Selection

Model-based feature selection uses a supervised machine learning model to judge the importance of each feature, and keeps only the most important ones.
The supervised model that is used for feature selection doesn’t need to be the same model that is used for the final supervised modeling.
The feature selection model needs to provide some measure of importance for each feature, so that they can be ranked by this measure.
Decision trees and decision tree–based models provide a feature_importances_ attribute, which directly encodes the importance of each feature.
Linear models have coefficients, which can also be used to capture feature importances by considering the absolute values.
As we saw in Chapter 2, linear models with L1 penalty learn sparse coefficients, which only use a small subset of features.
This can be viewed as a form of feature selection for the model itself, but can also be used as a preprocessing step to select features for another model.
In contrast to univariate selection, model-based selection considers all features at once, and so can capture interactions (if the model can capture them).
To use model-based feature selection, we need to use the SelectFromModel transformer:

In [48]:
from sklearn.feature_selection import SelectFromModel
from sklearn.ensemble import RandomForestClassifier

select = SelectFromModel(
    RandomForestClassifier(n_estimators=100, random_state=42),

SelectFromModel(estimator=RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',
            max_depth=None, max_features='auto', max_leaf_nodes=None,
            min_impurity_split=1e-07, min_samples_leaf=1,
            min_samples_split=2, min_weight_fraction_leaf=0.0,
            n_estimators=100, n_jobs=1, oob_score=False, random_state=42,
            verbose=0, warm_start=False),
        prefit=False, threshold='median')

The SelectFromModel class selects all features that have an importance measure of the feature (as provided by the supervised model) greater than the provided threshold.
To get a comparable result to what we got with univariate feature selection, we used the median as a threshold, so that half of the features will be selected.
We use a random forest classifier with 100 trees to compute the feature importances.
This is a quite complex model and much more powerful than using univariate tests.
Now let’s actually fit the model:

In [49]:, y_train)
X_train_l1 = select.transform(X_train)
print("X_train.shape: {}".format(X_train.shape))
print("X_train_l1.shape: {}".format(X_train_l1.shape))

X_train.shape: (284, 80)
X_train_l1.shape: (284, 40)

Again, we can have a look at the features that were selected (Figure 4-10):

In [50]:
mask = select.get_support()
# Visualize the mask; black is True, white is False
plt.matshow(mask.reshape(1, -1), cmap='gray_r')
plt.xlabel("Sample Index")

([], <a list of 0 Text yticklabel objects>)
Figure 4-10. Features selected by SelectFromModel using the RandomForestClassifier

This time, all but two of the original features were selected.
Because we told the machine learning classifier to select 40 features, some of the noise features are also selected.
Let's take a look at the performance:

In [51]:
X_test_l1 = select.transform(X_test)
score = LogisticRegression().fit(
    X_train_l1, y_train).score(X_test_l1, y_test)
print("Test score: {:.3f}".format(score))

Test score: 0.951

With the better feature selection, we also gained some improvements here.

Iterative Feature Selection

In univariate testing we used no model, while in model-based selection we used a single model to select features.
In iterative feature selection, a series of models are built, with varying numbers of features.
There are two basic methods: starting with no features and adding features one by one until some stopping criterion is reached, or starting with all features and removing features one by one until some stopping criterion is reached.
Because a series of models are built, these methods are much more computationally expensive than the methods we discussed previously.
One particular method of this kind is recursive feature elimination (RFE), which starts with all features, builds a model, and discards the least important feature according to the model.
Then a new model is built using all but the discarded feature, and so on until only a prespecified number of features are left.
For this to work, the model used for selection needs to provide some way to determine feature importance, as was the case for the model-based selection.
Here, we use the same random forest model that we used earlier, and get the results shown in Figure 4-11:

In [52]:
from sklearn.feature_selection import RFE

select = RFE(
    RandomForestClassifier(n_estimators=100, random_state=42),
                           n_features_to_select=40), y_train)
# Visualize the selected features:
mask = select.get_support()
plt.matshow(mask.reshape(1, -1), cmap='gray_r')
plt.xlabel("Sample Index")

([], <a list of 0 Text yticklabel objects>)
Figure 4-11. Features selected by recursive feature elimination with the random forest classifier model.

The feature selection got better compared to the univariate and model-based selection, but one feature was still missed.
Running this code also takes significantly longer than that for the model-based selection, because a random forest model is trained 40 times, once for each feature that is dropped.
Let’s test the accuracy of the logistic regression model when using RFE for feature selection:

In [53]:
X_train_rfe = select.transform(X_train)
X_test_rfe = select.transform(X_test)

score = LogisticRegression().fit(
    X_train_rfe, y_train).score(X_test_rfe, y_test)
print("Test Score: {:.3f}".format(score))

Test Score: 0.951

We can also use the model used inside the RFE to make predictions.
This uses only the feature set that was selected:

In [54]:
print("Test Score: {:.3f}".format(select.score(X_test, y_test)))

Test Score: 0.951

Here, the performance of the random forest used inside the RFE is the same as that achieved by training a logistic regression model on top of the selected features.
In other words, once we’ve selected the right features, the linear model performs as well as the random forest.

If you are unsure when selecting what to use as input to your machine learning algorithms, automatic feature selection can be quite helpful.
It is also great for reducing the amount of features needed—for example, to speed up prediction or to allow for more interpretable models.
In most real-world cases, applying feature selection is unlikely to provide large gains in performance.
However, it is still a valuable tool in the toolbox of the feature engineer.

Utilizing Expert Knowledge

Feature engineering is often an important place to use expert knowledge for a particular application.
While the purpose of machine learning in many cases is to avoid having to create a set of expert-designed rules, that doesn’t mean that prior knowledge of the application or domain should be discarded.
Often, domain experts can help in identifying useful features that are much more informative than the initial representation of the data.
Imagine you work for a travel agency and want to predict flight prices.
Let’s say you have a record of prices together with dates, airlines, start locations, and destinations.
A machine learning model might be able to build a decent model from that.
Some important factors in flight prices, however, cannot be learned.
For example, flights are usually more expensive during peak vacation months and around holidays.
While the dates of some holidays (like Christmas) are fixed, and their effect can therefore be learned from the date, others might depend on the phases of the moon (like Hanukkah and Easter) or be set by authorities (like school holidays).
These events cannot be learned from the data if each flight is only recorded using the (Gregorian) date.
However, it is easy to add a feature that encodes whether a flight was on, preceding, or following a public or school holiday.
In this way, prior knowledge about the nature of the task can be encoded in the features to aid a machine learning algorithm.
Adding a feature does not force a machine learning algorithm to use it, and even if the holiday information turns out to be noninformative for flight prices, augmenting the data with this information doesn’t hurt.

We'll now look at one particular case of using expert knowledge -- though in this case it might be more rightfully called "common sense".
The task is predicting bicycle rentals in front of Andreas’s house.
In New York, Citi Bike operates a network of bicycle rental stations with a subscription system.
The stations are all over the city and provide a convenient way to get around.
Bike rental data is made public in an anonymized form and has been analyzed in various ways.
The task we want to solve is to predict for a given time and day how many people will rent a bike in front of Andreas’s house—so he knows if any bikes will be left for him.
We first load the data for August 2015 for this particular station as a pandas DataFrame.
We resample the data into three-hour intervals to obtain the main trends for each day:

In [55]:
citibike = mglearn.datasets.load_citibike()
print("Citi Bike Data: \n{}".format(citibike.head()))

Citi Bike Data: 
2015-08-01 00:00:00     3.0
2015-08-01 03:00:00     0.0
2015-08-01 06:00:00     9.0
2015-08-01 09:00:00    41.0
2015-08-01 12:00:00    39.0
Freq: 3H, Name: one, dtype: float64

The following example shows a visualization of the rental frequencies for the whole month (Figure 4-12):

In [56]:
plt.figure(figsize=(10, 3))
xticks = pd.date_range(start=citibike.index.min(),
                       end=citibike.index.max(), freq='D')
plt.xticks(xticks, xticks.strftime("%a %m-%d"),
           rotation=90, ha="left")
plt.plot(citibike, linewidth=1)

<matplotlib.text.Text at 0x1164c21d0>
Figure 4-12. Number of bike rentals over time for a selected Citi Bike station.

Looking at the data, we can clearly distinguish day and night for each 24-hour interval.
The patterns for weekdays and weekends also seem to be quite different.
When evaluating a prediction task on a time series like this, we usually want to learn from the past and predict for the future.
This means when doing a split into a training and a test set, we want to use all the data up to a certain date as the training set and all the data past that date as the test set.
This is how we would usually use time series prediction: given everything that we know about rentals in the past, what do we think will happen tomorrow?
We will use the first 184 data points, corresponding to the first 23 days, as our training set, and the remaining 64 data points, corresponding to the remaining 8 days, as our test set.

The only feature that we are using in our prediction task is the date and time when a particular number of rentals occurred.
So, the input feature is the date and time—say, 2015-08-01 00:00:00 -- and the output is the number of rentals in the following three hours (three in this case, according to our DataFrame).
A (surprisingly) common way that dates are stored on computers is using POSIX time, which is the number of seconds since January 1970 00:00:00 (aka the beginning of Unix time).
As a first try, we can use this single integer feature as our data representation:

In [63]:
# Extract the target values (number of rentals):
y = citibike.values
print("y[:5]: {}".format(y[:5]))
# Convert to POSIX time by dividing by 10**9:
X = citibike.index.astype("int64").values.reshape(-1, 1) // 10**9
print("X[:5]: \n{}".format(X[:5]))

y[:5]: [  3.   0.   9.  41.  39.]

We first define a function to split the data into training and test sets, build the model, and visualize the result:

In [65]:
# Use the first 184 data points for training,
# and the rest for testing.
n_train = 184

# Function to evaluate and plot a regressor on
# a given test set.
def eval_on_features(features, target, regressor):
    # Split the given features into train and test sets.
    X_train, X_test = features[:n_train], features[n_train:]
    # Split the target array as well.
    y_train, y_test = target[:n_train], target[N-train:], y_train)
    print("Test-set R^2: {:.2f}".format(
        regressor.score(X_test, y_test)))
    y_pred = regressor.predict(X_test)
    y_pred_train = regressor.predict(X_train)
    plt.figure(figsize=(10, 3))
    plt.xticks(range(0, len(X), 8), xticks.strftime(
        "%a %m-%d"), rotation=90, ha="left")
    plt.plot(range(n_train), y_train, label="train")
    plt.plot(range(n_train, len(y_test) + n_train), y_test,
            '-', label="test")
    plt.plot(range(n_train), y_pred_train, '--',
            label="Training Prediction:")
    plt.plot(range(n_train, len(y_test) + n_train), y_pred,
            '--', label="Testing Prediction:")
    plt.legend(loc=(1.01, 0))

We saw earlier that random forests require very little preprocessing of the data, which makes this seem like a good model to start with.
We use the POSIX time feature X and pass a random forest regressor to our eval_on_features function.
Figure 4-13 shows the result:

In [ ]: