<jupyter_start><jupyter_text>Telco Customer Churn: Exploration**Author:** <NAME> -- **Description:** in this notebook, we attempt to model the [Telco Customer Churn](https://www.kaggle.com/blastchar/telco-customer-churn) problem using a Machine Learning algorithm called Random Forest.**Results:** We were able to achieve an average of 0.828 RoC AUC score over 7 weeks of test data, without loosing significant performance as time passed.The rest of this notebook is organized in the following way:1. The *Reading* section contains the code that read the dataset from my google drive2. *Preprocessing* describes all pre-training procedure applied to the data, such as conversions and empty string handling3. *Visualizing* describes the data acording to its features and provides some insights which might help when modeling the problem4. *Learning* contains all training-validation strategy5. Finally, Section *Evaluating* tests thetrained model over previously-separated, unseen data  Setup<jupyter_code>! pip -q install shap
from math import ceil

import pandas as pd
import numpy as np
import seaborn as sns

import matplotlib.pyplot as plt

import warnings
sns.set()

np.random.seed(4012808)<jupyter_output><empty_output><jupyter_text>ReadingThe data was read from a virtual directory contained my personal google drive.This part can be executed by simply downloading the [Telco Churn Dataset](https://www.kaggle.com/blastchar/telco-customer-churn) dataset and adjusting the path described below.<jupyter_code>DATASET = '/content/drive/My Drive/Colab Notebooks/churn/WA_Fn-UseC_-Telco-Customer-Churn.csv'
x = pd.read_csv(DATASET)
x.head(2)<jupyter_output><empty_output><jupyter_text>Preprocessing<jupyter_code>CATEGORICAL_FIELDS = ('gender', 'partner', 'dependents', 'phoneservice',
                      'multiplelines', 'internetservice', 'contract',
                      'paymentmethod', 'onlinesecurity', 'onlinebackup',
                      'deviceprotection', 'techsupport', 'streamingtv',
                      'streamingmovies', 'paperlessbilling', 'churn',
                      'seniorcitizen')

DISCRETE_FIELDS = ('tenure',)
CONTINUOUS_FIELDS = ('monthlycharges', 'totalcharges')


def preprocess(x):
    print(f'Preprocessing dataframe of shape {x.shape}')

    x = x.copy()

    x.columns = list(map(str.lower, x.columns))  # lower col names

    x = boolean_field_as_text(x, 'seniorcitizen') 
    x = apply_to_cols(x, CATEGORICAL_FIELDS, str.lower)
    x = convert_empty_fields_to_nan(x,
                                    CATEGORICAL_FIELDS
                                    + DISCRETE_FIELDS
                                    + CONTINUOUS_FIELDS)

    x = cast_fields(x, DISCRETE_FIELDS, 'int32')
    x = cast_fields(x, CONTINUOUS_FIELDS, 'float32')

    print('Done')

    return x

def boolean_field_as_text(x, field):
    m = x[field] == 1

    x.loc[m, field] = 'yes'
    x.loc[~m, field] = 'no'

    return x

def apply_to_cols(x, fields, fn):
    for f in fields:
        x[f] = x[f].apply(fn)

    return x

def convert_empty_fields_to_nan(x, fields):
    for f in fields:
        print(f'Converting ` ` in {f} to nan...')
        x[f] = x[f].replace(r'^\s*$', np.nan, regex=True)
    
    return x

def cast_fields(x, fields, dtype):
    for f in fields:
        print(f'Converting {f} to {dtype}...')
        x[f] = x[f].astype(dtype)

    return x

x = preprocess(x)
def validate(x):
    print(f'Validating dataframe of shape {x.shape}...')

    assert_fields_are_numeric(x, DISCRETE_FIELDS + CONTINUOUS_FIELDS)

    assert_contains_expected_values(x, ('gender', ('female', 'male')),
                                       ('seniorcitizen', ('yes', 'no')),
                                       ('partner', ('yes', 'no')),
                                       ('dependents', ('yes', 'no')),
                                       ('phoneservice', ('yes', 'no')),
                                       ('onlinesecurity', ('yes', 'no', 'no internet service')),
                                       ('onlinebackup', ('yes', 'no', 'no internet service')),
                                       ('deviceprotection', ('yes', 'no', 'no internet service')),
                                       ('techsupport', ('yes', 'no', 'no internet service')),
                                       ('streamingtv', ('yes', 'no', 'no internet service')),
                                       ('streamingmovies', ('yes', 'no', 'no internet service')),
                                       ('paperlessbilling', ('yes', 'no')),
                                       ('churn', ('yes', 'no')))
    
    assert_does_not_contain_nan(x, CATEGORICAL_FIELDS + DISCRETE_FIELDS + ('monthlycharges',))

    print('Done.')
    

def assert_contains_expected_values(x, *fields):
    for field, expected in fields:
        expected = set(expected)
        actual = set(x[field].unique())

        if actual - expected:
            raise ValueError(f'Field `{field}` contains non-expected values:',
                             actual - expected)


def assert_fields_are_numeric(x, fields):
    malformed = []

    for f in fields:
        if x[f].dtype not in ('int32', 'float32'):
            malformed.append((f, x[f].dtype))
    
    if malformed:
        raise ValueError(f'Fields `{malformed}` is not of a type numeric. '
                         f'Check the input data.')

def assert_does_not_contain_nan(x, fields):
    malformed = []

    for f in fields:
        nans = x[f].isna()
        if nans.any():
            malformed.append((f, nans.sum()))
    
    if malformed:
        raise ValueError(f'Unexpected NaN values were found in the following fields:\n  '
                         + '\n  '.join(map(str, malformed)))
    

validate(x)
pd.DataFrame(x.dtypes).T
x.head().round(1)<jupyter_output><empty_output><jupyter_text>Exploring  Checking the distribution of the dataset's features  Categorical ColumnsCategorical features seem relatively balanced, with the exception of `dependents`, `phoneservice`, `contract`.Important note: the most undersampled value is `no`, in the `phoneservice` column, present in at most 800 samples. This is might be of some concern, as we will further split the data for testing (30%, hence approx. 240 test samples will present `phoneservice=no`). However, this number is still above the magic number $30$ from the t-student distribution-rule. For the time being, we will keep exploring as it is and maybe handle this attention point in the future.<jupyter_code>#@title

def fields_histograms(x, fields, cols=3):
    rows = ceil(len(fields) / cols)

    plt.figure(figsize=(12, 3*rows))

    for ix, f in enumerate(fields):
        plt.subplot(rows, cols, ix + 1)
        sns.countplot(x=f, hue='churn', data=x, palette='Accent')

    plt.tight_layout();
fields_histograms(x, ['gender', 'phoneservice', 'multiplelines', 'churn'])
    fields_histograms(x, ['contract', 'streamingtv', 'streamingmovies',
                        'onlinesecurity', 'deviceprotection', 'partner',
                        'paperlessbilling', 'onlinebackup', 'seniorcitizen',
                        'internetservice', 'paymentmethod', 'dependents',
                        'techsupport'])<jupyter_output><empty_output><jupyter_text>Discrete ColumnsWe have a large amount of "novice" users and a reasonable amount of long-lifetime users.<jupyter_code>#@title

def fields_kdes(x, fields):
    for ix, f in enumerate(fields):
        sns.displot(data=x, x=f, hue='churn', kde=True)
        plt.tight_layout();
fields_kdes(x, DISCRETE_FIELDS)<jupyter_output><empty_output><jupyter_text>Continuous ColumnsColumn `TotalCharges` is a long-right tail distribution.This is not ideal for parametric models, which tend to saturate over the median line.  We will handle this in the following sections.<jupyter_code>fields_kdes(x, CONTINUOUS_FIELDS);<jupyter_output><empty_output><jupyter_text>Relationship between The Variables DefinitionsWe account for data unbalance by undersampling the most frenquent class.This is important for the distribution charts, which are highly sensible to un-normalized data volumes. Results* Churning clients seem to have concentrated in samples with low tenure, indicating users with high tenure are less likely to evade* Churning clients seem to correlate more strongly with samples with high monthly-charges, indicating users with lower bills are less prone to evade<jupyter_code>from imblearn.under_sampling import RandomUnderSampler

xn = x.dropna()
xr, yr = RandomUnderSampler().fit_resample(xn, xn.churn)
xr = pd.DataFrame(xr, columns=xn.columns)
xr = xr.astype(dict(zip(x.columns, x.dtypes)))
sns.pairplot(xr.sample(frac=.3), hue='churn');<jupyter_output><empty_output><jupyter_text>Fixing Irregular Features  Normalizing `totalcharges`We can try to normalize the column `totalcharges` using the `log` function,which will crunch the higher values that are utterly sparse, while maintaining the lowest values relatively separated.The result is added to the `totalcharges_log` column in the frame. Now, the distribution of users that have churned seems closer to the uniform one.<jupyter_code>x['totalcharges_log'] = np.log(x.totalcharges)

fields_kdes(x, ['totalcharges_log']);<jupyter_output><empty_output><jupyter_text>Estimating the missing values for `totalcharges`Considering only 11 samples are missing, a cleanest solution would be to just drop these lines with `DataFramedropna()`.However, assuming new samples would eventually become available to us as time pass, *linear regressing* these features might generate better results at the time.  Notice we are not using regression in the *machine-learning* context (we do not concern ourselves with generalization/reproducibility), but as a mere statistical tool to infer missing data through correlation.<jupyter_code>from sklearn.pipeline import make_pipeline
from sklearn.compose import make_column_transformer, make_column_selector
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.linear_model import LinearRegression

def estimate_missing_values(x, cats, numerics, target):
    x = x.copy()

    training = ~x[target].isna()
    print(f'{training.sum()} `{target}` samples used as training')

    r = make_pipeline(
            make_column_transformer(
                (StandardScaler(), list(numerics)),
                (OneHotEncoder(), list(cats)),
            ),
            LinearRegression())

    r.fit(x[training], x[target][training])
    x.loc[~training, target] = r.predict(x.loc[~training])

    print(f'{(~training).sum()} `{target}` examples were estimated')

    sns.displot(data=x, x=target, hue=training, kde=True)
    return x
CATEGORICAL_FEATURES = set(CATEGORICAL_FIELDS) - {'churn'}
NUMERIC_FEATURES = set(DISCRETE_FIELDS + CONTINUOUS_FIELDS) - {'totalcharges', 'totalcharges_log'}
TARGET = 'totalcharges_log'

z = estimate_missing_values(x, CATEGORICAL_FEATURES, NUMERIC_FEATURES, TARGET)<jupyter_output><empty_output><jupyter_text>The estimated samples are barely visible in the distribution curve above, considering how few samples there are to begin with.  Checking for Data Separability<jupyter_code>TARGET = 'churn'

CATEGORICAL_FEATURES = set(CATEGORICAL_FIELDS) - {TARGET}
NUMERIC_FEATURES = set(DISCRETE_FIELDS + CONTINUOUS_FIELDS) - {'totalcharges'}
#@title

print('Using the following columns as features:')
print('target variable:', TARGET)
print('Available features:',
      *CATEGORICAL_FEATURES,
      *NUMERIC_FEATURES,
      sep='\n  ')
class NamedStandardScaler(StandardScaler):
    def __init__(self, copy=True, with_mean=True, with_std=True, columns=None):
        super().__init__(copy=copy, with_mean=with_mean, with_std=with_std)
        self.columns = columns

    def get_feature_names(self):
        return self.columns


t = make_column_transformer(
        (NamedStandardScaler(columns=list(NUMERIC_FEATURES)), list(NUMERIC_FEATURES)),
        (OneHotEncoder(), list(CATEGORICAL_FEATURES)))

z = t.fit_transform(x)
print('Using features:', list(NUMERIC_FEATURES))
from sklearn.manifold import TSNE
from sklearn.cluster import MiniBatchKMeans


def display_dataset_distrib(x, true_labels):
    a, b = TSNE(perplexity=50).fit_transform(x).T
    c = MiniBatchKMeans(n_clusters=3)
    y = c.fit_predict(x)

    sns.scatterplot(x=a, y=b, hue=y, style=true_labels, palette='Set1')

    return c.cluster_centers_


plt.figure(figsize=(12, 6))
extracted_clusters = display_dataset_distrib(x[NUMERIC_FEATURES], x[TARGET]);

print('clusters:')
print(extracted_clusters.round(2))<jupyter_output><empty_output><jupyter_text>The method [T-SNE](https://scikit-learn.org/stable/modules/generated/sklearn.manifold.TSNE.html) will attempt to create a low-dimension visualization for the data while preserving its local and global distributions, where "locality" is described by its `perplexity` parameter.The KMeans, on the other hand, uses $L_2$ the distance to find data cluster.From inspection, cluster `1` seem to contain fewer churning samples than other ones. Although it is also the smallest. Further investigation must be performed in order to determine if these results are of some importance, but so far it seems the variable `churn` does not match the most prominant projections of the dataset's distribution, with respect to its features.  LearningIn many problems, samples are distributed according to a temporal dimension.This usually induces sub-groups of samples that are farthest from the training set, which might result in worse general interpretation by models.To account for this time-shifting effect and make this problem more interesting, I added an artificial column named `created_at`, which contains a random date for each sample.  Additionally, we generate noise that grows linearly as time passes, and apply it to the numeric variables.<jupyter_code>#@title

def random_dates(start, end, n, unit='D'):
    """Generate random dates.

    Reference: https://stackoverflow.com/questions/50559078/generating-random-dates-within-a-given-range-in-pandas
    """
    start = pd.to_datetime(start)
    end = pd.to_datetime(end)
    
    ndays = (end - start).days + 1
    return pd.to_timedelta(np.random.rand(n) * ndays, unit=unit) + start

def add_linear_noise(x, col, rate):
    dt_to_first = (x.created_at - x.created_at.min()).dt.days
    dt_to_first /= dt_to_first.max()

    a, b = np.quantile(x.tenure, [.05, .95])
    s = dt_to_first * np.random.rand(len(x)) * rate * (b - a)

    x[col] += s
x['created_at'] = random_dates('2020-01-01', '2020-10-30', len(x))
NOISE_RATE = {
    'tenure': .05,
    'monthlycharges': .05,
    'totalcharges_log': .05
}

for col, rate in NOISE_RATE.items():
    add_linear_noise(x, col, rate)<jupyter_output><empty_output><jupyter_text>TrainingThis procedure will occur as follows:1. Data is split according to its temporal dimension, so we can train using the past and evaluating considering present/future samples.2. Most models are sensitive to data unbalancing and will overfit over the most common class. We use the [SMOTTENN](https://imbalanced-learn.readthedocs.io/en/stable/generated/imblearn.combine.SMOTEENN.htmlimblearn.combine.SMOTEENN) method to rebalance the dataset by creating samples that respect the classes original distributions and elimiate the most unsure/diverging outliers.3. Grid-search is used to find for the model's most appropriate hyper-parameters.    a. The validation strategy used here is `RepeatedStratifiedKFold`. `2-Fold` was employed, as it yields a large chunk of evaluation data without undersampling training data given to the model.    The `10 repetitions` add variability to the procedure, fortifying it    against evaluation noise.    b. The metric [roc-auc](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_auc_score.html) is used as evaluating function to the optimization process, which represents more granuar information with respect to the models overall generability when compared to accuracy or class-balanced accuracy.<jupyter_code>TEST_SPLIT = .3
TRAIN_SIZE = ceil((1 - TEST_SPLIT) * len(x))
t = make_column_transformer(
        (NamedStandardScaler(columns=list(NUMERIC_FEATURES)), list(NUMERIC_FEATURES)),
        (OneHotEncoder(), list(CATEGORICAL_FEATURES)))

x = x.sort_values('created_at')

train, test = x[:TRAIN_SIZE], x[TRAIN_SIZE:]
z_train, z_test = t.fit_transform(train), t.transform(test)
y_train, y_test = train[TARGET], test[TARGET]
from imblearn.combine import SMOTEENN

zr_train, yr_train = SMOTEENN().fit_resample(z_train, y_train)
from sklearn.model_selection import GridSearchCV, RepeatedStratifiedKFold
from sklearn import metrics
from sklearn.ensemble import RandomForestClassifier


def train(model, grid, x, y, scoring='roc_auc', cv=3):
    return (GridSearchCV(model, grid, scoring=scoring, cv=cv,
                         n_jobs=-1, verbose=1)
            .fit(x, y))

g = train(
    model=RandomForestClassifier(),
    grid=dict(
        n_estimators=[20, 50, 100, 200],
        class_weight=[None, 'balanced'],
        max_depth=[None, 5, 10]
    ),
    x=zr_train,
    y=yr_train,
    cv=RepeatedStratifiedKFold(n_splits=2, n_repeats=10))<jupyter_output><empty_output><jupyter_text>Describing the trained modelThe most impacting features' values over the classification of a sample as a probable churning users are:1. Contract2. Tenure3. Tech-support4. Payment method5. Monthly charges<jupyter_code>#@title

def describe(grid):
    r = pd.DataFrame(grid.cv_results_)
    r = r[r.rank_test_score == 1].round(3).reset_index()

    print('best parameters:', grid.best_params_)
    print(f'Avg test score: {r.mean_test_score[0]}. Std: {r.std_fit_time[0]}')
    print(f'Avg fit time: {r.mean_fit_time[0]}. Std: {r.std_fit_time[0]}')

describe(g)
#@title

import shap

explainer = shap.TreeExplainer(g.best_estimator_)


def decode_feature(f):
    for i, n in reversed(list(enumerate(CATEGORICAL_FEATURES))):
        # We reverse this listing so ocurrences like "x17" are
        # not replaced by sub-contained indices, such as "x1".
        f = f.replace(f'x{i}', n)
    return f

features = [decode_feature(f)
            for f in t.get_feature_names()]

_z = pd.DataFrame(z_train, columns=features)
shap_values = explainer.shap_values(_z)

shap.initjs()
shap.summary_plot(shap_values, _z)<jupyter_output><empty_output><jupyter_text>Evaluating<jupyter_code>def evaluate(model, x, y):
    p = model.predict(x)
    pa = model.predict_proba(x)[:, 1]

    print('RoC AUC:', metrics.roc_auc_score(y, pa).round(3))
    print('Accuracy:', metrics.accuracy_score(y, p).round(3))
    print('Class-balanced accuracy:', metrics.balanced_accuracy_score(y, p).round(3))
    print()

    print('Classification Report:')
    print(metrics.classification_report(y, p))


evaluate(g, z_test, y_test)<jupyter_output><empty_output><jupyter_text>ROC AuC Curve for Churn:yes Class<jupyter_code>#@title

from sklearn.preprocessing import label_binarize


def calculate_roc_auc_curve(x, y):
    """Calculate the roc auc curve for an estimator.

    Ref: https://scikit-learn.org/stable/auto_examples/model_selection/plot_roc.html#sphx-glr-auto-examples-model-selection-plot-roc-py
    """
    y = label_binarize(y, classes=['no', 'yes'])
    n_classes = y.shape[1]

    y_score = g.predict_proba(x)

    fpr = dict()
    tpr = dict()
    roc_auc = dict()
    
    i = 1
    fpr[i], tpr[i], _ = metrics.roc_curve(y, y_score[:, i])
    roc_auc[i] = metrics.auc(fpr[i], tpr[i])

    # Compute micro-average ROC curve and ROC area
    fpr["micro"], tpr["micro"], _ = metrics.roc_curve(y.ravel(), y_score[:, i].ravel())
    roc_auc["micro"] = metrics.auc(fpr["micro"], tpr["micro"])

    return fpr, tpr, roc_auc

fpr, tpr, roc_auc = calculate_roc_auc_curve(z_test, y_test)

plt.figure()
lw = 2
plt.plot(fpr[1], tpr[1], color='darkorange',
         lw=lw, label='ROC curve (area = %0.2f)' % roc_auc[1])
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic example')
plt.legend(loc="lower right")
plt.show()<jupyter_output><empty_output><jupyter_text>Evaluating model degradation over timeWe do not observe a degradation in the model's performance over the course of almost two months after training, hinting this model is either strongly resilient to temporal changes or that this domain does not shift in an aggressive fashion.Further experimentation must be conducted after this two month period in order to assert if this behavior can be further generalized to longer periods, although I cannot recommend in good conscience maintaining a production model without re-training for more than a few weeks/month.<jupyter_code>def backtest_one_step(model, x, y, created_at):
    p = model.predict(x)

    report = pd.DataFrame({
        'p': model.predict(x),
        'pa': model.predict_proba(x)[:, 1],
        'y': y,
        'created_at': created_at,
    })

    results = [(
        week,
        metrics.accuracy_score(g.y, g.p),
        metrics.balanced_accuracy_score(g.y, g.p),
        metrics.roc_auc_score(g.y, g.pa),
        len(g),
        (g.y != g.p).sum(),
    )
    for week, g in report.groupby(report.created_at.dt.isocalendar().week)]

    return pd.DataFrame(results,
                        columns=['week', 'acc', 'balanced_acc',
                                 'roc_auc', 'samples', 'misses'])

r = backtest_one_step(g, z_test, y_test, test.created_at)

r.round(2)
#@title

m = (r.melt(id_vars=['week', 'samples'],
            value_vars=['acc', 'balanced_acc', 'roc_auc'])
      .rename(columns={'variable': 'metric'}))

sns.lmplot(x='week', y='value', hue='metric', data=m, scatter=False, height=5,
           aspect=16/9)
sns.scatterplot(x='week', y='value', hue='metric', size='samples', data=m,
                legend=False)

plt.ylim(.3, 1);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11408815903197926, 'fraction_numerical': 0.007505799936314425, 'mean_word_length': 3.8052459016393443, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 66, 'https://': 7, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '467822', 'license': 'Apache-2.0', 'n_tokens_mistral': 7167, 'n_tokens_neox': 6841, 'n_words': 1918, 'path': 'notebooks/supervised/classification/churn/exploration.ipynb', 'repo_name': 'lucasdavid/algorithms-in-tensorflow', 'size': '21997'}