Cox Modeling Using Lifelines

Extending from our notebook on the math and intuition behind the Cox Model let’s do a practical example using real data.

The Data

We’ll use the Telco Customer Churn dataset on Kaggle, which is basically a bunch of client records for a telecom company, where the goal is to predict churn (Churn) and the duration it takes for churn to happen (tenure).

%pylab inline

import pandas as pd

Populating the interactive namespace from numpy and matplotlib

Lot of potentially-useful information here, this notebook does a particularly good job exploring the data.

df = pd.read_csv('https://raw.githubusercontent.com/'
                 'treselle-systems/customer_churn_analysis/'
                 'master/WA_Fn-UseC_-Telco-Customer-Churn.csv')

df.head().T

There’s about a 3-1 split between churn and not churn

df['Churn'].value_counts() / len(df)

No     0.73463
Yes    0.26537
Name: Churn, dtype: float64

And the difference in tenure distribution between the two is pretty stark– lot of right-censored data in our No group.

fig, ax = plt.subplots(figsize=(12, 10))

for idx, group in df.groupby('Churn'):
    group['tenure'].hist(ax=ax, bins=20, label=idx)

ax.legend();

There’s a lot more here than is illustrative to use, so let’s pare down this dataset a bit.

First, we’ll move customerID into the index as well as drop the TotalCharges variables, as there’s an obvious, uninformative, colinearity between that and tenure.

simple = df.set_index('customerID', drop=True)
simple = simple.drop('TotalCharges', axis=1)

We’ll keep the numeric tenure and MonthlyCharges columns, but here, I’m narrowing down to categorical features that I think might be useful, then dummy-ing them for model consumption.

simple = pd.get_dummies(simple,
                        columns=['gender', 'SeniorCitizen', 'InternetService',
                                 'Contract', 'PaymentMethod', 'Churn'],
                        drop_first=True).select_dtypes('number')

simple.head().T

Printing some simple statistics, a few things stand out:

• The gender breakup is about ⁵⁰⁄₅₀
• About a third of the pop has DSL internet
• Over half of the userbase is going month-to-month

simple.describe().T

Fitting a Model

So if we take what we’ve got and fit a simple model to it, we can get an easy glimpse at the significance of our features in the p column.

from lifelines import CoxPHFitter

cph = CoxPHFitter().fit(simple, 'tenure', 'Churn_Yes')

cph.print_summary()

I was negative-surprised to see that the Electronic Check and Mailed Check payment methods led to a significant dip in the Surival Curve– me and my Planet Fitness membership are all too familiar with the “out of sight, out of mind” business model. Credit card (automatic), on the other hand, is not significant, because the 4th payment method is Bank transfer (automatic), which basically achieves the same thing

df['PaymentMethod'].value_counts()

Electronic check             2365
Mailed check                 1612
Bank transfer (automatic)    1544
Credit card (automatic)      1522
Name: PaymentMethod, dtype: int64

Plotting

Of course, we can use the model and records from our dataset to predict the Survival Curve

sample = simple.sample(50, random_state=12)

cph.predict_survival_function(sample).plot(figsize=(12, 10),
                                           legend=False);

Additionally, we can use plot_covariate_groups() to hold everything equal, save for one attribute, then examine the relative effects of different values on the population.

Here, we can see that having fiber internet tends to make users churn faster.

cph.plot_covariate_groups('InternetService_Fiber optic',
                          values=[0, 1],
                          figsize=(10, 5));

But why do you suppose that is?

Comparing Two Populations

My immediate hypothesis is that it’s got something to do with how the client has to pay. Indeed, at all levels of contract, someone with fiber internet pays more. So it could be a fiber thing, or could be a money thing.

simple.groupby(['Contract_Two year', 'Contract_One year',
                'InternetService_Fiber optic'])['MonthlyCharges'].agg(['mean', 'std'])

Let’s pare down to just Fiber users to do some investigation.

fiber = simple[simple['InternetService_Fiber optic'] == 1]
fiber = fiber[['tenure', 'Contract_One year',
               'Contract_Two year', 'MonthlyCharges',
               'gender_Male', 'Churn_Yes']]

cph = CoxPHFitter().fit(fiber, 'tenure', 'Churn_Yes')

cph.print_summary()

Actually, it looks like having a higher monthly charge means you’re more likely to stick around.

seed = 42

low = fiber[fiber['MonthlyCharges'] < 80].sample(50, random_state=seed)
hi = fiber[fiber['MonthlyCharges'] > 100].sample(50, random_state=seed)

Well that’s odd.

fig, ax = plt.subplots(figsize=(12, 10))

cph.predict_survival_function(low).plot(ax=ax, c='b', alpha=.5, legend=False)
cph.predict_survival_function(hi).plot(ax=ax, c='r', alpha=.5, legend=False);

from matplotlib.lines import Line2D
custom_lines = [Line2D([0], [0], color='b'),
                Line2D([0], [0], color='r')]

ax.legend(custom_lines, ['Low', 'High']);

It’s obvious visually, but if we wanted to say with statistical confidence that the high and low paying populations had different survival curve distributions, we could do so via the built-in logrank_test()

from lifelines.statistics import logrank_test

results = logrank_test(low['tenure'], hi['tenure'],
                       event_observed_A=low['Churn_Yes'],
                       event_observed_B=hi['Churn_Yes'])
results.print_summary()

Here’s a screengrab of the code. I’m tired of fighting with LaTex.

Image('images/i_hate_latex.PNG')

Re-fit and…

cph = CoxPHFitter()

cph.fit(m2m, duration_col='tenure', event_col='Churn_Yes')

cph.print_summary()

Damn it.

cph.plot_covariate_groups('MonthlyCharges', values=[80, 100])

<matplotlib.axes._subplots.AxesSubplot at 0x1d3611169b0>

Starting Over

So there’s obviously more at play here than just MonthlyCharges.

Let’s re-pull our original dataset with all of the features, but only for records that are going month-to-month on Fiber internet

m2m_index = m2m.index

a = df.set_index('customerID')
b = a.loc[m2m_index].drop('TotalCharges', axis=1)
c = pd.get_dummies(b, drop_first=True)

m2m_wide = c

Of course, fitting the model with everything but the kitchen sink yields us more than we need.

cph = CoxPHFitter()
cph.fit(m2m_wide, duration_col='tenure', event_col='Churn_Yes')
cph.print_summary()

Okay, so it looks like it’s less to do with contracts and money, and more to do with additional features. Perhaps a convincing-enough interpretation here is “if the user gets upsold on everything, they’re not probably not shopping around.”

Visually inspecting that p column again, let’s hold out all of the features that look significant and re-fit.

sig_cols = ['gender_Male', 'Partner_Yes', 'MultipleLines_Yes',
            'OnlineSecurity_Yes', 'OnlineBackup_Yes', 'DeviceProtection_Yes',
            'TechSupport_Yes', 'PaperlessBilling_Yes', 'PaymentMethod_Electronic check',
            'PaymentMethod_Mailed check']

cph = CoxPHFitter()
cph.fit(m2m_wide[sig_cols + ['tenure', 'Churn_Yes']], duration_col='tenure', event_col='Churn_Yes')
cph.print_summary()

Check Assumptions

This looks promising. But did we make a procedural blunder in all of our data transformation?

lifelines gives us an awesome tool that we can use to simply check the Cox Model assumptions

cph.check_assumptions(training_df=m2m_wide[sig_cols + ['tenure', 'Churn_Yes']])

The ``p_value_threshold`` is set at 0.01. Even under the null hypothesis of no violations, some
covariates will be below the threshold by chance. This is compounded when there are many covariates.
Similarly, when there are lots of observations, even minor deviances from the proportional hazard
assumption will be flagged.

With that in mind, it's best to use a combination of statistical tests and visual tests to determine
the most serious violations. Produce visual plots using ``check_assumptions(..., show_plots=True)``
and looking for non-constant lines. See link [A] below for a full example.

Re-checking, we can see that lifelines doesn’t hate the work we did, this time~

cph.check_assumptions(m2m_wide[sig_cols + ['tenure', 'Churn_Yes']])

Proportional hazard assumption looks okay.

Of course, this now means that our visual interpretability goes down

m2m_sample = m2m_wide[sig_cols + ['tenure', 'Churn_Yes']].sample(50)

cph.predict_cumulative_hazard(m2m_sample).plot(legend=False,
                                               figsize=(15, 10));

and while lifelines allows us to inspect the baseline survival curves, for the various strata

cph.baseline_survival_.plot(figsize=(15, 8))
plt.gcf().suptitle('{}'.format(cph.strata));

This gets us to the fundamental tension of “model accuracy” and “model interpretability” and is one of the many corners of “Data Science is more of an art than a science”