Interaction Terms in Python

Vanilla OLS

Say we’ve got a dataset

from warnings import filterwarnings
filterwarnings('ignore')

from yellowbrick.datasets import load_bikeshare

X, y = load_bikeshare()

Where we measure the number of DC-area bikes rented

y.head()

0    16
1    40
2    32
3    13
4     1
Name: riders, dtype: int64

Based on a number of features

X.head()

And just shoving everything into a Logistic Regression seems to work… alright

from statsmodels.api import OLS
import statsmodels.api as sm

X_const = sm.add_constant(X)
model = OLS(y, X_const).fit()
model.summary()

But at the same time, have good reason to believe that there’s some colinearity/interaction at play with our features.

from yellowbrick.features import Rank2D

visualizer = Rank2D(algorithm='pearson', size=(600, 500))
visualizer.fit_transform(X)
visualizer.poof();

C:\Users\Nick\Anaconda3\lib\site-packages\yellowbrick\features\rankd.py:215: YellowbrickWarning: RankD plots may be clipped when using matplotlib v3.1.1, upgrade to matplotlib v3.1.2 or later to fix the plots.
  warnings.warn(msg, YellowbrickWarning)

Interaction Terms

From here, a good data scientist will take the time to do exploratory analysis and thoughtful feature engineering– this is the “More Art than Science” adage you hear so often.

But we’re trying to be home by 5, so how do we cram everything in and see what shakes out?

Getting Values

Thankfully, the PolynomialFeatures object in sklearn has us mostly-covered.

It’s originally used to generate sequences of (b_i1 * x_i) + (b_i2 * x_i^2) + ... for each feature in X, taking us from n features to 2^n features (in the case of PolynomialFeatures(degree=2), anyhow).

We’re not interested in polynomials, per se, but if you squint, the same itertools magic™ that powers the backend of this can also be used to provide all pairwise feature combinations, with minimal rewriting. They provide this, ez pz, with the interaction_only flag.

And so we go from

X.shape

(17379, 12)

To an expected $\frac{p * (p - 1)}{2}$ pairs, plus our original p features, plus a bias term

p = len(X.columns)
(p * (p-1)) / 2 + p + 1

79.0

Coolio

from sklearn.preprocessing import PolynomialFeatures

interaction = PolynomialFeatures(degree=2, interaction_only=True, include_bias=False)
X_inter = interaction.fit_transform(X)
X_inter_const = sm.add_constant(X_inter)

X_inter_const.shape

(17379, 79)

Now we brazenly throw it into a new model, and… oh. A lot of features called x

model = OLS(y, X_inter_const).fit()
model.summary()

Feature Names?

The interaction object keeps track of the names of our feature combinations… sort of

print(interaction.get_feature_names())

['x0', 'x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x9', 'x10', 'x11', 'x0 x1', 'x0 x2', 'x0 x3', 'x0 x4', 'x0 x5', 'x0 x6', 'x0 x7', 'x0 x8', 'x0 x9', 'x0 x10', 'x0 x11', 'x1 x2', 'x1 x3', 'x1 x4', 'x1 x5', 'x1 x6', 'x1 x7', 'x1 x8', 'x1 x9', 'x1 x10', 'x1 x11', 'x2 x3', 'x2 x4', 'x2 x5', 'x2 x6', 'x2 x7', 'x2 x8', 'x2 x9', 'x2 x10', 'x2 x11', 'x3 x4', 'x3 x5', 'x3 x6', 'x3 x7', 'x3 x8', 'x3 x9', 'x3 x10', 'x3 x11', 'x4 x5', 'x4 x6', 'x4 x7', 'x4 x8', 'x4 x9', 'x4 x10', 'x4 x11', 'x5 x6', 'x5 x7', 'x5 x8', 'x5 x9', 'x5 x10', 'x5 x11', 'x6 x7', 'x6 x8', 'x6 x9', 'x6 x10', 'x6 x11', 'x7 x8', 'x7 x9', 'x7 x10', 'x7 x11', 'x8 x9', 'x8 x10', 'x8 x11', 'x9 x10', 'x9 x11', 'x10 x11']

Assuming that a bevy of xi xjs aren’t much more useful than the output of model.summary() we can do some hacky, Python nonsense to decode a bit.

For starters, we want to create a dictionary that maps xi to its corresponding feature name in our dataset.

We’ll use the itertools.count() function, as it’s basically enumerate, but plays better with generator expressions.

from itertools import count

x_to_feature = dict(zip(('x{}'.format(i) for i in count()), X.columns))
x_to_feature

{'x0': 'season',
 'x1': 'year',
 'x2': 'month',
 'x3': 'hour',
 'x4': 'holiday',
 'x5': 'weekday',
 'x6': 'workingday',
 'x7': 'weather',
 'x8': 'temp',
 'x9': 'feelslike',
 'x10': 'humidity',
 'x11': 'windspeed'}

Next, you know it’s Sound Data Science™ when we break out

import re

Finally, this little diddy goes through and makes the appropriate substitutions for xi to their respective feature names, then replaces spaces with underscores so pandas references isn’t such a chore

# necessary so `x11` gets a chance to swap before 
# before `x1` leaves us with "season1" where we wanted
# "windspeed"
feature_keys = list(x_to_feature.keys())[::-1]

features = []

for feature in interaction.get_feature_names():
    for key in feature_keys:
        feature_name = x_to_feature[key]
        feature = re.sub(key, feature_name, feature)
    feature = re.sub(' ', '_', feature)
    features.append(feature)

Then we’ll slap these feature names onto our big ol’ dataset

import pandas as pd
import statsmodels.api as sm

X_all_inter = pd.DataFrame(X_inter, columns=features)
X_all_inter_const = sm.add_constant(X_all_inter)

And refit

model = OLS(y, X_all_inter_const).fit()
model.summary()

And that’s Data Science!

from IPython.display import Image
Image('images/we_did_it.jpg')