В данной работе предлагаю рассмотреть достаточно известный датасет (он был симулирован, но достаточно мал и показателен):
Human Resources Analytics (https://www.kaggle.com/ludobenistant/hr-analytics/downloads/HR_comma_sep.csv)
import numpy as np
import pandas as pd
import seaborn as sns
import xgboost as xgb
import catboost as cat
import itertools
from matplotlib import pyplot as plt
from sklearn import preprocessing
from sklearn.model_selection import train_test_split
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_val_score, cross_val_predict
from sklearn.model_selection import learning_curve
from sklearn.metrics import roc_auc_score, confusion_matrix, auc, roc_curve
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import learning_curve
from sklearn.model_selection import ShuffleSplit
from mlxtend.feature_selection import SequentialFeatureSelector as SFS
from hyperopt import fmin, hp, STATUS_OK, Trials, tpe
%matplotlib inline
Процесс сбора данных:
Датасет был сгенерирован.
Описание решаемой задачи и ее ценность:
Предотвратить уход сотрудников из компании, определить причины.
Определить какие признаки влияют на то, что сотрудник уволится (профессиональное выгорание, отсутствие повышения и прочее). Датасет актуален для превентивного анализа ситуации в компании и сохранения ценных кадров. В нем содержаться основные количественные характеристики сотрудников и их работы в компании.
Описанная задача является задачей бинарной классификации.
Описание признаков и целевой переменной:
- left - целевая переменная - ушел ли сотрудник из компании;
- satisfaction_level - уровень удовлетворенности (0-1);
- last_evaluation - время, прошедшее с последнего повышения (в годах);
- number_project - количество законченных проектов за все время работы;
- average_montly_hours - среднее количество часов работы в месяц;
- time_spend_company - количество лет, которые сотрудник проработал в компании;
- work_accident - были ли особые случае (несчастный случай на работе и прочее);
- promotion_last_5years- было ли повышение за последние 5 лет;
- sales - департамент сотрудника;
- salary - оценочный уровень зарплаты (градации).
dataset = pd.read_csv('HR_comma_sep.csv')
dataset.shape
dataset.head()
Проверим наличие пропущенных значений в датасете:
dataset.isnull().sum()
Посмотрим основную статистическую информацию:
dataset.describe()
Наша целевая переменная left, посмотрим на ее распределение:
plt.figure(figsize=(6, 4))
plt.scatter(range(dataset.shape[0]), np.sort(dataset.left.values))
plt.xlabel('index', fontsize=12)
plt.ylabel('left', fontsize=12)
plt.show()
Проверим, какие типы признаков у нас есть в датасете:
dtype_df = dataset.dtypes.reset_index()
dtype_df.columns = ["Count", "Column Type"]
dtype_df.groupby("Column Type").aggregate('count').reset_index()
dtype_df.loc[:10,:]
В наличии у нас 2 категориальные переменные (sales и salary), остальные являются количественными.
Посмотрим на матрицу корреляции для датасета:
corr = dataset.corr()
plt.figure(figsize=(12, 8))
sns.heatmap(corr, annot=True, cmap="YlGnBu")
Как видим, сильно корреляции между признаками нет (< 0.7). Считаем, что они линейнонезависимы.
Пропущенные значения и выбросы отсутствуют (согласно гистрограммам)
Распределения для признаков имеют вид:
for column in ['satisfaction_level', 'last_evaluation', 'number_project',
'average_montly_hours', 'time_spend_company']:
plt.figure(figsize=(8, 6))
plt.title(column)
sns.distplot(dataset[column])
plt.show()
plt.figure(figsize=(10, 8))
ax = sns.countplot(dataset['sales'])
plt.title('sales')
plt.show()
plt.figure(figsize=(10, 8))
ax = sns.countplot(dataset['salary'])
plt.title('salary')
plt.show()
Проведем анализ влияния признаков на целевую переменную:
def plot_distribution(df, var, target, yl=4, **kwargs):
row = kwargs.get('row', None)
col = kwargs.get('col', None)
facet = sns.FacetGrid(df, hue=target, aspect=4, row=row, col=col)
facet.map(sns.kdeplot, var, shade=True)
facet.set(xlim=(0, df[var].max()), ylim=(0, yl))
facet.add_legend()
plt.show()
plot_distribution(dataset, 'satisfaction_level', 'left')
Из распределения satisfaction_level видно, что люди часто увольняются, когда уровень удовлетворение мал, в тоже время небольшая часть людей увольняется и при высоком уровне - возможно находят лучше или их переманивают.
plot_distribution(dataset, 'last_evaluation', 'left')
Из распределения last_evaluation видно люди уходят на другое место если их долго не повышают (есть 2 горба - первый это возле максимального времени пересмотра, второй примерно в середине)
plot_distribution(dataset, 'number_project', 'left', yl=1.5)
Четких закономерностей не видно.
plot_distribution(dataset, 'average_montly_hours', 'left', yl=0.015)
Четких закономерностей не видно: увольняются обычные сотрудники с ~140 часов/мес и в такой же степени трудоголики с 250+ часов.
plot_distribution(dataset, 'time_spend_company', 'left', yl=1)
Из распределения time_spend_company видно, что вероятность, того человек уволится растет в течении времени его работы в компании.
Построим scatter plot, предварительно преобразуя категориальные признаки с помощью LabelEncoder
sns.set(style="ticks")
labeled_dataset = dataset.copy()
salary_le = preprocessing.LabelEncoder()
sales_le = preprocessing.LabelEncoder()
labeled_dataset['salary'] = salary_le.fit_transform(labeled_dataset['salary'])
labeled_dataset['sales'] = sales_le.fit_transform(labeled_dataset['sales'])
sns.pairplot(labeled_dataset, hue="left")
Рассмотрим более подробно признаки:
plt.figure(figsize=(12, 10))
sns.factorplot(y="satisfaction_level",x="left",data=dataset,kind="box")
sns.factorplot(y="last_evaluation", x="left", data=dataset, kind="box")
sns.factorplot(y="number_project", x="left", data=dataset, kind="box")
sns.factorplot(y="average_montly_hours", x="left", data=dataset, kind="box")
sns.factorplot(y="time_spend_company", x="left", data=dataset, kind="box")
sns.factorplot(y="sales", x="left", data=dataset, kind="box")
sns.factorplot(y="salary", x="left", data=dataset, kind="box")
plot_distribution(labeled_dataset, 'sales', 'left', yl=0.4)
plot_distribution(labeled_dataset, 'salary', 'left', yl=3)
Инсайты:
- более склонны к увольнениям люди с большим количеством рабочих часов в месяц и уровнем удовлетворенности
ниже среднего;
- люди, работающие больше 7 лет, практически не увольняются, вне зависимости от уровня удовлетворенности;
- среди людей с низким уровнем компенсации высокая текучка;
- высокое количество увольнений в областях: HR, Accounting, Technical (20% - 30% от всех сотрудников).
Общие мысли:
Важными факторами являются:
- низкая зарплата;
- определенный департамент;
- низкий уровень удовлетворенности;
- большие переработки;
- большое количество проектов.
В качестве кодирование категориальных переменных будет использовать LabelEncoder.
Для нормализации датасета используем StandardScaler.
salary_le = preprocessing.LabelEncoder()
sales_le = preprocessing.LabelEncoder()
dataset['salary'] = salary_le.fit_transform(dataset['salary'])
dataset['sales'] = sales_le.fit_transform(dataset['sales'])
# может пригодиться
Norma = preprocessing.StandardScaler()
dataset_norm = Norma.fit_transform(dataset)
dataset.head()
Поделим датасет на тренировочную и тестовые выборки (используем ненормированный датасет).<
dataset_v1 = dataset.copy()
y = dataset_v1['left'].values
dataset_v1.drop(['left'], axis=1, inplace=True)
X = dataset_v1.copy()
В качестве метрики выберем основной ROC_AUC, дополнительно будем смотреть confusion_matrix.
Задача классификации - предположим, что нас интересуют вероятности (позволяет более точно определять степень уверенности модели).
Наша целевая переменная не очень сбалансирована, поэтому использовать accuracy выглядит не очень хорошей затеей.
В качестве пробных моделей проведем отбор с дефолтными параметрами среди следующих моделей:
- LogisticRegression;
- RandomForest;
- XGBoost;
- CATboost;
Все вышеописанные модели вполне применимы для задач бинарной классификации и можно использовать ненормализованный датасет.
clf_lm = LogisticRegression(random_state=42)
preds_lm = cross_val_predict(clf_lm, X, y, cv=5)
print('LogisticRegression: %s ROC AUC' % round(roc_auc_score(y, preds_lm), 4))
clf_rf = RandomForestClassifier(random_state=42)
preds_rf = cross_val_predict(clf_rf, X, y, cv=5)
print('RandomForest : %s ROC AUC' % round(roc_auc_score(y, preds_rf), 4))
clf_xgb = xgb.XGBClassifier(random_state=42)
preds_xgb = cross_val_predict(clf_xgb, X, y, cv=5)
print('XGBoost : %s ROC AUC' % round(roc_auc_score(y, preds_xgb), 4))
clf_cat = cat.CatBoostClassifier(random_seed=42)
preds_cat = cross_val_predict(clf_cat, X, y, cv=5)
print('CATBoost : %s ROC AUC' % round(roc_auc_score(y, preds_cat), 4))
RandomForest в этой задаче показывает впечатляющие результаты. Построим ROC кривую для полученных результатов:
fpr_rf, tpr_rf, threshold = roc_curve(y, preds_rf)
roc_auc_rf = auc(fpr_rf, tpr_rf)
fpr_xgb, tpr_xgb, threshold = roc_curve(y, preds_xgb)
roc_auc_xgb = auc(fpr_xgb, tpr_xgb)
fpr_lr, tpr_lr, threshold = roc_curve(y, preds_lm)
roc_auc_lr = auc(fpr_lr, tpr_lr)
plt.figure(figsize=(12, 8))
plt.title('Receiver Operating Characteristic')
plt.plot(fpr_rf, tpr_rf, 'b', label = 'AUC RF = %0.2f' % roc_auc_rf)
plt.plot(fpr_xgb, tpr_xgb, 'g', label = 'AUC XGB = %0.2f' % roc_auc_xgb)
plt.plot(fpr_lr, tpr_lr, 'r', label = 'AUC LR = %0.2f' % roc_auc_lr)
plt.legend(loc = 'lower right')
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.show()
fig_size = [8, 6]
plt.rcParams["figure.figsize"] = fig_size
def plt_matrix(c_matrix, names, title, normalize=False):
np.set_printoptions(precision=2)
plt.figure()
plt.imshow(c_matrix, interpolation='nearest', cmap=plt.cm.Blues)
plt.title(title)
plt.colorbar()
tick_marks = np.arange(len(names))
plt.xticks(tick_marks, names)
plt.yticks(tick_marks, names)
fmt = '.2f' if normalize else 'd'
thresh = c_matrix.max() / 2.
for i, j in itertools.product(range(c_matrix.shape[0]), range(c_matrix.shape[1])):
plt.text(j, i, format(c_matrix[i, j], fmt),
horizontalalignment="center",
color="white" if c_matrix[i, j] > thresh else "black")
plt.tight_layout()
plt.ylabel('True label')
plt.xlabel('Predicted label')
class_names = ['0', 'left']
baseline_matrix = confusion_matrix(y, preds_rf)
plt_matrix(baseline_matrix, class_names, 'Random Forest baseline')
Проведем поиск оптимальных параметров, посмотрим, получится ли улучшить модель.
Сделаем это следующим образом:
1. Зафиксируем random_seed (для воспроизведения и равных условий разных моделей);
2. Будем использовать Stratified разбиения (target не сбалансирован);
- у нас 15000 записей, будешь разбивать на 5 фолдов.
3. Подбирать параметры моделей будем с помощью модуля hyperopt.
4. Описание параметров для тюнинга:
- n_estimators - количество деревьев в лесу;
- max_features - количество признаков, используемых для разделения;
- max_depth - глубина дерева;
- class_weight - вес классов.
Возможно увеличение количество параметров для тюнинга, но вышеописанные параметры являются наиболее
важными для модели
Выберем отложенную выборку, на ней мы проверим финальный результат модели:
X_tuning, X_test, y_tuning, y_test = train_test_split(X, y, test_size=0.2, random_state=42, stratify=y)
print(X_tuning.shape, X_test.shape, y_tuning.shape, y_test.shape)
Попробуем сначала class_weight с дефолтными параметрами:
for weight in range(1, 20):
clf_rf = RandomForestClassifier(class_weight={0: 1, 1: weight}, random_state=42)
preds_rf = cross_val_predict(clf_rf, X_tuning, y_tuning, cv=5)
print('RandomForest %s : %s ROC AUC' % (weight, round(roc_auc_score(y_tuning, preds_rf), 4)))
Вес для значений таргета 1 равный 1 является лучшим. Пока остановимся на нем.
Перейдем к затюниванию модели:
# наиболее простой вариант использования модуля hyperopt
# стоит обратить внимание на факт, что модуль делает 20 случайных подборов параметров до начала настройки.
def hyperopt_train_test(hpparams):
# глобальный счетчик итераций для удобства и вывода информации.
global counter, result
counter += 1
# словарь подбираемых параметров
params_est = {
'n_estimators': int(hpparams['n_estimators']),
'max_features': hpparams['max_features'],
'max_depth': int(hpparams['max_depth']),
'random_state': 42,
'n_jobs': 8,
'class_weight': {0: 1, 1: 1},
}
# код обучения и получения результата метрики:
clf_rf = RandomForestClassifier(**params_est)
preds_rf = cross_val_predict(clf_rf, X_tuning, y_tuning, cv=5)
result_cur = round(roc_auc_score(y_tuning, preds_rf), 6)
if result_cur > result:
print('Best iteration: %s, roc_auc: %s, params: %s' % (counter, result_cur, hpparams))
result = result_cur
return result_cur
# пространство для поиска параметров
# более детально можно ознакомится в документации.
space4dt = {
'n_estimators': hp.quniform('n_estimators', 10, 100, 10),
'max_features': hp.quniform('max_features', 0.1, 1, 0.1),
'max_depth': hp.quniform('max_depth', 6, 20, 1),
}
def f(params):
metrics = hyperopt_train_test(params)
return {'loss': -metrics, 'status': STATUS_OK}
trials = Trials()
counter, result = 0, 0
best = fmin(f, space4dt, algo=tpe.suggest, max_evals=100, trials=None)
print('best: ', best)
Лучшие параметры (мы уперлись в 2 границы, значит улучшать еще есть куда, даже для этих параметров):
clf_rf = RandomForestClassifier(max_depth=20, max_features=0.5, n_estimators=100, random_state=42)
preds_rf = cross_val_predict(clf_rf, X, y, cv=5)
print('RandomForest : %s ROC AUC' % round(roc_auc_score(y, preds_rf), 4))
clf_rf = RandomForestClassifier(max_depth=20, max_features=0.5, n_estimators=100, random_state=42)
preds_rf = clf_rf.fit(X_tuning, y_tuning)
# отложенная выборка
print('RandomForest test : %s ROC AUC' % round(roc_auc_score(y_test, clf_rf.predict(X_test)), 4))
Для желающих - можно расширить пространство признаков и их значений. Будете служить богу вычислений :)
Итак, логически (поэтому и не совпадают номера глав) мы подошли генерации фичей (признаков).
Предлагаю попробовать автоматическую генерацию линейных комбинаций признаков. Можно сгенерировать на первый вгляд логичные признаки, но все они попадут в наше поле зрения при генерации.
Что мы имеем сейчас:
dataset.head()
Сгенерируем много много признаков:
dataset_generated = dataset.copy()
print('\nBefore transformation: ', dataset_generated.shape)
columns = [
'satisfaction_level', 'last_evaluation', 'number_project', 'average_montly_hours',
'time_spend_company', 'salary'
]
for i1, col1 in enumerate(columns):
for i2, col2 in enumerate(columns):
if col1 == col2:
dataset_generated['%s_%s_0' % (col1, col2)] = np.log(X[col1] + 1)
dataset_generated['%s_%s_1' % (col1, col2)] = X[col1] / (X[col2] + 1)
dataset_generated['%s_%s_2' % (col1, col2)] = X[col1] * X[col2]
print('\nAfter transformation: ', dataset_generated.shape)
Теперь проведем отбор признаков. Для этого воспользуемся или самописным жадным отбором или уже написанным для нас в пакете mlxtend. Используем вот такой вариант отбора - SequentialFeatureSelector (есть несколько вариантов, можно почитать документацию по пакету). Наш вариант добавляет по одному признаку в датасет, и выбирает лучшую комбинацию по достижению максимального количества признаков. Для нашей задачи давайте установим максимальное значение признаков равное 12 и минимальное 1.
y_sfs = dataset_generated['left'].values
dataset_generated.drop(['left'], axis=1, inplace=True)
X_sfs = dataset_generated.copy()
model = RandomForestClassifier(random_state=42)
sfs1 = SFS(model, k_features=(1, 12), forward=True, floating=False,
verbose=2, scoring='roc_auc', cv=3, n_jobs=-1)
sfs1 = sfs1.fit(X_sfs.as_matrix(), y_sfs)
print('Results:')
print(sfs1.k_feature_idx_)
print(sfs1.k_score_)
Собирем датасет только из отобранных фич:
dataset_best_features = pd.DataFrame()
columns = []
for elem in sfs1.k_feature_idx_:
dataset_best_features[dataset_generated.columns[elem]] = dataset_generated[dataset_generated.columns[elem]]
dataset_best_features['left'] = y_sfs
dataset_best_features.head()
Попробуем еще раз затюнить параметры модели для полученного датасета:
y_best = dataset_best_features['left'].values
dataset_best_features.drop(['left'], axis=1, inplace=True)
X_best = dataset_best_features.copy()
X_tuning_b, X_test_b, y_tuning_b, y_test_b = train_test_split(
X_best, y_best, test_size=0.2, random_state=42, stratify=y)
print(X_tuning.shape, X_test.shape, y_tuning.shape, y_test.shape)
clf_rf = RandomForestClassifier(random_state=42)
preds_rf = cross_val_predict(clf_rf, X_best, y_best, cv=5)
print('RandomForest : %s ROC AUC' % round(roc_auc_score(y_best, preds_rf), 4))
Как видим, добавление и отбор признаков в нашем случае не помог, значение целевой метрики уходшилось, попробуем исправить ситуацию тюнингом параметров модели:
# наиболее простой вариант использования модуля hyperopt
# стоит обратить внимание на факт, что модуль делает 20 случайных подборов параметров до начала настройки.
def hyperopt_train_test(hpparams):
# глобальный счетчик итераций для удобства и вывода информации.
global counter, result
counter += 1
# словарь подбираемых параметров
params_est = {
'n_estimators': int(hpparams['n_estimators']),
'max_features': hpparams['max_features'],
'max_depth': int(hpparams['max_depth']),
'random_state': 42,
'n_jobs': 8,
'class_weight': {0: 1, 1: 1},
}
# код обучения и получения результата метрики:
clf_rf = RandomForestClassifier(**params_est)
preds_rf = cross_val_predict(clf_rf, X_best, y_best, cv=5)
result_cur = round(roc_auc_score(y_best, preds_rf), 6)
if result_cur > result:
print('Best iteration: %s, roc_auc: %s, params: %s' % (counter, result_cur, hpparams))
result = result_cur
return result_cur
# пространство для поиска параметров
# более детально можно ознакомится в документации.
space4dt = {
'n_estimators': hp.quniform('n_estimators', 10, 200, 10),
'max_features': hp.quniform('max_features', 0.1, 1, 0.1),
'max_depth': hp.quniform('max_depth', 6, 24, 1),
}
def f(params):
metrics = hyperopt_train_test(params)
return {'loss': -metrics, 'status': STATUS_OK}
trials = Trials()
counter, result = 0, 0
best = fmin(f, space4dt, algo=tpe.suggest, max_evals=200, trials=None)
print('best: ', best)
clf_rf = RandomForestClassifier(max_depth=22, max_features=0.8, n_estimators=110, random_state=42)
clf_rf.fit(X_tuning_b, y_tuning_b)
# отложенная выборка
print('RandomForest test : %s ROC AUC' % round(roc_auc_score(y_test_b, clf_rf.predict(X_test_b)), 4))
clf_rf = RandomForestClassifier(max_depth=22, max_features=0.8, n_estimators=110, random_state=42)
preds_rf = cross_val_predict(clf_rf, X_best, y_best, cv=5)
print('RandomForest : %s ROC AUC' % round(roc_auc_score(y_best, preds_rf), 4))
clf_rf = RandomForestClassifier(max_depth=22, max_features=0.8, n_estimators=110, random_state=42)
preds_rf = cross_val_predict(clf_rf, X_best, y_best, cv=5)
print('RandomForest : %s ROC AUC' % round(roc_auc_score(y_best, preds_rf), 4))
Тюнинг модели улучшил результат, возможно нужно большее пространство параметров и больше степеней свободы (добавить параметры модели). Результаты получились на исходном датасете с улучшенными параметрами hyperopt:
RandomForest : 0.9859 ROC AUC (исходный датасет);
RandomForest : 0.986 ROC AUC (переработанные признаки в датасете);
def plot_learning_curve(estimator, title, X, y, ylim=None, cv=None,
n_jobs=1, train_sizes=np.linspace(.1, 1.0, 5)):
plt.figure()
plt.title(title)
if ylim is not None:
plt.ylim(*ylim)
plt.xlabel("Training examples")
plt.ylabel("Score")
train_sizes, train_scores, test_scores = learning_curve(
estimator, X, y, cv=cv, n_jobs=n_jobs, train_sizes=train_sizes)
train_scores_mean = np.mean(train_scores, axis=1)
train_scores_std = np.std(train_scores, axis=1)
test_scores_mean = np.mean(test_scores, axis=1)
test_scores_std = np.std(test_scores, axis=1)
plt.grid()
plt.fill_between(train_sizes, train_scores_mean - train_scores_std,
train_scores_mean + train_scores_std, alpha=0.1,
color="r")
plt.fill_between(train_sizes, test_scores_mean - test_scores_std,
test_scores_mean + test_scores_std, alpha=0.1, color="g")
plt.plot(train_sizes, train_scores_mean, 'o-', color="r",
label="Training score")
plt.plot(train_sizes, test_scores_mean, 'o-', color="g",
label="Cross-validation score")
plt.legend(loc="best")
return plt
title = "RandomForest default"
cv = ShuffleSplit(n_splits=10, test_size=0.2, random_state=0)
estimator = RandomForestClassifier(random_state=42)
plot_learning_curve(estimator, title, X_best, y, ylim=(0.7, 1.01), cv=cv, n_jobs=8)
title = "RandomForest tuned"
cv = ShuffleSplit(n_splits=10, test_size=0.2, random_state=0)
estimator = RandomForestClassifier(max_depth=20, max_features=0.5, n_estimators=100, random_state=42)
plot_learning_curve(estimator, title, X_best, y, (0.7, 1.01), cv=cv, n_jobs=8)
plt.show()
Модель достаточно хорошо сходится, хотя еще данные бы не помешали.
Проверим confusion_matrix для новых параметров:
clf_rf = RandomForestClassifier(max_depth=20, max_features=0.5, n_estimators=100, random_state=42)
preds_rf_tuned = cross_val_predict(clf_rf, X_best, y, cv=5)
tuned_matrix = confusion_matrix(y, preds_rf_tuned)
clf_rf = RandomForestClassifier(random_state=42)
preds_rf_baseline = cross_val_predict(clf_rf, X, y, cv=5)
baseline_matrix = confusion_matrix(y, preds_rf_baseline)
plt_matrix(baseline_matrix, class_names, 'Random Forest baseline')
plt_matrix(tuned_matrix, class_names, 'Random Forest tuned')
plt.show()
Как видно выше по проекту, отложенной выборкой мы несколько раз пользовались для оценки качества модели.
Выбранная метрика ROC_AUC близка к 1. Наша модель практически не ошибается и уверена в своих вероятностях. Как видно из confusion matrix в пункте 9. Модель неправильно идентифицирует всего 123 человек. Accuracy модели > 99%.