Основной вопрос

Dataset "Airlines Customer satisfaction" -- содержит данные перелетов и клиентов, оценки различных фичей при перелёте и уровень удовлетворенности, характеристики клиентов и характеристики перелётов.

Необходимо предсказать, будет ли лояльный клиент, который летит по личным мотивам, бизнес - классом удовлетворен своим обслуживанием, учитывая детали значений других параметров.

Данные и библиотеки

# загружаем библиотеки и данные import numpy as np import pandas as pd dataAT = pd.read_csv("/work/Invistico_Airline2.csv")

# загружаем библиотеки import matplotlib.pyplot as plt import seaborn as sns; sns.set(style='white') %matplotlib inline import warnings warnings.filterwarnings("ignore") from sklearn import tree from sklearn import metrics from sklearn.tree import DecisionTreeClassifier from sklearn.tree import export_graphviz from sklearn.model_selection import train_test_split from sklearn.model_selection import GridSearchCV from sklearn.preprocessing import LabelEncoder from sklearn.model_selection import StratifiedKFold from sklearn.model_selection import cross_val_score from sklearn.linear_model import LogisticRegression from sklearn.ensemble import VotingClassifier from sklearn.ensemble import AdaBoostClassifier from sklearn.ensemble import BaggingClassifier from sklearn.ensemble import GradientBoostingClassifier from sklearn.ensemble import RandomForestClassifier from sklearn.datasets import load_breast_cancer from sklearn.metrics import accuracy_score from sklearn.metrics import confusion_matrix from sklearn.metrics import plot_confusion_matrix from sklearn.metrics import roc_auc_score , roc_curve from sklearn.model_selection import cross_val_score from scipy import stats

!pip install statsmodels import statsmodels.api as sm

!pip install pydotplus==2.0.2 from six import StringIO from IPython.display import Image import pydotplus

!pip install shap import shap from shap import TreeExplainer, Explanation from shap.plots import waterfall

!pip install pycebox from pycebox.ice import ice, ice_plot

!pip install eli5==0.13.0 import eli5 from eli5.sklearn import PermutationImportance

Обработка данных

Основной вопрос ещё раз: необходимо предсказать, будет ли лояльный клиент, который летит по личным мотивам, бизнес - классом удовлетворен своим обслуживанием, учитывая детали значений других параметров.

#переименую для удобства dataAT.rename(columns = { 'Flight Distance':'flightDis', 'Customer Type':'custType', 'Age':'age', 'Class':'class', 'Type of Travel':'travelType', 'Seat comfort':'seat', 'Departure/Arrival time convenient':'depArr', 'Food and drink': 'foodDrink', 'Gate location': 'gate', 'Inflight wifi service' : 'wifi', 'Inflight entertainment' : 'entertainment', 'Online support' : 'onlineSup', 'Ease of Online booking' : 'easyBooking', 'On-board service' : 'service', 'Leg room service' : 'LegRoom', 'Baggage handling' : 'baggageH', 'Checkin service' : 'checking', 'Online boarding' : 'onlineBoar', 'Inflight service' : 'infServ', 'Cleanliness' : 'clean', 'Departure Delay in Minutes' : 'depDelay', 'Arrival Delay in Minutes' : 'arrDelay' }, inplace = True)

#проверка dataAT.isnull().sum()

# заполним пропущенные значения задержки прибытия средним значением dataAT['arrDelay'] = dataAT['arrDelay'].fillna(dataAT['arrDelay'].mean())

dataCa = dataAT.copy() # переведем из object в int dataCa['class'] = dataCa['class'].replace({ 'Business': 2, 'Eco': 1, 'Eco Plus' : 0}) dataCa['travelType'] = dataCa['travelType'].replace({ 'Business travel': 0, 'Personal travel': 1}) #использую гет дамис для преобразования dataCa = pd.get_dummies(dataCa, columns = ['custType'], drop_first=True)

# укорачиваю датасет по требованиям dataCa.drop(dataCa[dataCa["travelType"] == 0].index, inplace = True) #только персональная поездка dataCa.drop(dataCa[dataCa["class"] == 1].index, inplace = True) #НЕ eco dataCa.drop(dataCa[dataCa["class"] == 0].index, inplace = True) #НЕ eco plus dataCa.drop(dataCa[dataCa["custType_disloyal Customer"] == 1].index, inplace = True) #только лояльные dataCa.shape

dataCa = dataCa.drop(["id", "Unnamed: 0", "travelType", "class", "custType_disloyal Customer"], axis=1) dataCa.head(2)

#визуализирую данные warnings.filterwarnings('ignore') dataCa.hist(bins=50, figsize=(20,15)) plt.show()

Также я решила обработать данные статистическими вопросами, чтобы попробовать получить какую-то дополнительную информацию из данных, которая не очевидна при простом взгляде на визуализацию данных.

Рассматривая визуализацию данных меня заинтересовали вопросы касающиеся расстояния:

1)Зависит ли удовлетворенность полётом от расстояния, который пролетели пассажиры?

Нулевая гипотеза: удовлетворенность полётом не зависит от расстояния.

Для начала я рассмотрела распределение расстояния, далее был рассмотрен размер эффекта. При сравнении моего размера эффекта с разбросом данных, выяснилось, что размер эффекта небольшой, следовательно ассоциация между независимой и зависимой переменной имеет небольшую силу.

2)А есть ли разница в расстоянии у довольных и недовольных пассажиров?

Нулевая гипотеза: разницы в расстоянии у довольных и недовольных пассажиров нет.

Так как с помощью теста Шапиро-Уилка выяснилось, что данные распределены не нормально, для ответа на этот вопрос был проведён непараметрический тест (Манна-Уитни) вместо стандартного t-теста. P-value оказался большим -- pvalue=0.9014101643714322 - больше 0.05, поэтому у меня не достаточно данных, чтобы отвергнуть нулевую гипотезу. И я НЕ могу сделать вывод, что эффект, скорее всего, существует и в популяции, и в выборке.

Следующее, что мне захотелось рассмотреть - это есть ли какая-нибудь взаимосвязь между полом пассажира и его удовлетворенность полётом.

Вопрос: Интересно, зависит ли удовлетворенность полётом от пола?

Нулевая гипотеза: удовлетворенность полётом у мужчин и женщин одинаковая.

Так как этот статистический вопрос касается двух категориальных переменных - пол (м/ж) и удовлетворенность (да/нет) я использовала хи-квадрат. pvalue оказался больше 0.05, таким образом у меня также недостаточно оснований отвергнуть нулевую гипотезу.

3)Также мне стало интересно, а что вообще влияет на удовлетворенность пассажиров, точнее какая из фич делает пассажиров довольными.

Для этого я построила логистическую регрессию. Выделяя переменные главной была переменная "satisfaction". Получившаяся мной модель оказалась статистически значима, так как LLR p-value == 3.596e-08, что <0,05.

Также оказалось, что для всех параметров, кроме "wifi", "onlineBoar" p-значение получилось больше определенного уровня (0.05), поэтому переменная-предиктора, кроме этих двух не имеет статистически значимую связь с переменной ответа в модели.

У получившейся модели Log-Likelihood: -125.54, а L-Null (нулевой модели): -161.76, поэтому она лучше, чем нулевая, т.е. модель без предикторов, также она имела показатель качества модели псевдо-R^2=0.2239, значит показатель ближе к 0, а не к 1, поэтому у меня недостаточно данных для определения, что делает пассажиров удовлетворенными.

Построение моделей для предсказания

Следующим этапом было построение моделей для предсказания, чтобы ответить на основной вопрос. Зависимой переменной является "satisfaction" (y), а остальные независимые (x). Выборку поделю в соотношении 1/3, так как у меня не так много данных - всего 514 наблюдений.

В процессе ведения дневников были выполнены clustering + subgroup discovery, где задача состояла в том, чтобы разделить клиентов на группы, чтобы понять, какие услуги на борту больше влияют на удовлетворенность, чтобы улучшить какие-то "отстающие" услуги и сделать в будущем персонализированную рекламу, но мне они не подходят для ответа на мой основной вопрос, то есть бесполезны, поэтому я не включаю их в итоговый отчёт.

Однако были получены следующие результаты, на которые можно обратить внимание: кластер, который был более удовлетворен -- по средним оценкам выставили более плохие оценки, нежели, чем кластер наиболее удовлетворенный. С помощью subgroup discovery я заметила тот факт, что модель поделила данные на 3 подгруппы и из полученных данных выяснилось, что параметры повторяются (низкая оценка еды и напитков). Возможно авиакомпании следует пересмотреть предлагаемое меню, напитки.

Также были построены модели регрессии и регрессии с регуляризацией -- они представлены кратко при построении ансамблей.

Дерево

Задача типа классификации - метрика "accuracy".

dataTr = dataCa.copy() # переведем из object в int dataTr['Gender'] = dataTr['Gender'].replace({ 'Female': 1, 'Male' : 0}) dataTr['satisfaction'] = dataTr['satisfaction'].replace({ 'neutral or dissatisfied': 0, 'satisfied': 1})

y = dataTr.satisfaction X = dataTr.drop(["satisfaction"], axis=1) X.head()

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=1/3, random_state = 0)

treeTA = DecisionTreeClassifier(max_depth = 3) treeTA.fit(X_train,y_train)

train_score = treeTA.score(X_train, y_train) test_score = treeTA.score(X_test, y_test) print('На тренировочном массиве модель работает с точностью {0}, \ а на тестовом массиве с точностью {1:6.2f}'.format \ (train_score, test_score) )

На тестовом качество получилось даже лучше, чем на тренировочном

tree_pred = treeTA.predict(X_test)

print('Метрики:', metrics.accuracy_score(y_test, tree_pred), ';', metrics.f1_score(y_test, tree_pred, pos_label = 0 ))

plot_confusion_matrix(treeTA, X_test, y_test, display_labels = [0, 1])

Модель ошибается 0 раз, когда он на самом деле недоволен. И ошибается 4 раза когда он на самом деле доволен, а модель думает, что он не доволен.

Строю дерево :

dot_data = StringIO() export_graphviz(treeTA, out_file=dot_data, filled=True, rounded=True, special_characters=True, feature_names = X.columns) graph = pydotplus.graph_from_dot_data(dot_data.getvalue()) Image(graph.create_png())

Для самого левого листа голубого цвета value = [0, 8], что значит, что в этот лист попало 0 представитель из числа недовольных и 8 из числа довольных -- предсказываем, что будет недоволен полётом. Или самый правый, нижний лист - [35, 13]. То есть 35 недовольных и 13 довольных - предсказываем, что будет недовольный. Качество дерева получилось хорошим, на тестовой выборе данных даже лучше, чем на тренировочной.

Случайный лес

Поиск параметров для случайного леса:

def search_param(model, param, X_train, y_train, X_test, y_test, area=range(1, 11), msg=True, plot=True, seed=None): import matplotlib.pyplot as plt import time score_list = [] if msg: print('# точность время') for abc in area: start = time.time() rfc = eval(model + '(' + param + '=' + str(abc) + ', random_state=' + str(seed) + ')') rfc.fit(X_train, y_train) s = rfc.score(X_test, y_test) end = time.time() score_list.append(s) if msg: print("%-3d %10f %7f" % (abc, s, end - start)) if plot: plt.plot(list(area), score_list) return list(area)[score_list.index(max(score_list))]

Поиск n_estimators -- количество деревьев в лесу

Поиск max_depth -- максимальной глубины

Поиск min_samples_split -- минимальное кол-во выборок, необходимых для разделения

Поиск min_samples_leaf -- минимальное кол-во примеров, требуемое для нахождения в листовом узле

Кросс-валидация. Я решила использовать перекрестную проверку, чтобы выявить наилучшие параметры.

randFCl = RandomForestClassifier(random_state=15) param = {'n_estimators': [i for i in range(38, 51)], 'max_depth': [i for i in range(20, 25)]} gscv = GridSearchCV(randFCl, param, cv=3, n_jobs=-1, verbose=1) gscv.fit(X_train, y_train)

gscv.best_params_

Создаю модель по полученным данным:

best_c = gscv.best_estimator_ imp = pd.DataFrame(best_c.feature_importances_, index=X_train.columns, columns=['importance']) imp.sort_values('importance').plot(kind='barh', figsize=(12, 8))

print('На тренировочном массиве с точностью', best_c.score(X_train, y_train), 'На тестовом массиве с точностью', best_c.score(X_test, y_test))

Разница не сильно большая, но кажется есть небольшое переобучение.

dtc_proba=best_c.predict_proba(X_test) dtc_proba=dtc_proba[:,1] auc=roc_auc_score(y_test, dtc_proba) print('Случайный лес: ROC AUC=%.3f' % (auc))

Ансамбли

Основа для ансамблей:

# лог. регрессия logReg = LogisticRegression(solver='lbfgs', multi_class='multinomial',random_state=1) logReg.fit(X_train, y_train) train_score_log = logReg.score(X_train, y_train) test_score_log = logReg.score(X_test, y_test) print('Регрессия : на тренировочном массиве модель работает с точностью {0}, \ а на тестовом массиве с точностью {1:6.2f}'.format \ (train_score_log, test_score_log) ) # с регуляризацией logReg2 = LogisticRegression(C=0.0002, solver='lbfgs', multi_class='multinomial',random_state=1) logReg2.fit(X_train, y_train) train_score_logred = logReg2.score(X_train, y_train) test_score_logreg = logReg2.score(X_test, y_test) print('Регрессия с регуляризацией: на тренировочном массиве модель работает с точностью {0}, \ а на тестовом массиве с точностью {1:6.2f}'.format \ (train_score_logred, test_score_logreg) ) # дерево print('Дерево: на тренировочном массиве модель работает с точностью {0}, \ а на тестовом массиве с точностью {1:6.2f}'.format \ (train_score, test_score) ) # лес train_rf_score = best_c.score(X_train, y_train) test_rf_score = best_c.score(X_test, y_test) print('Лес: на тренировочном массиве модель работает с точностью {0}, \ а на тестовом массиве с точностью {1:6.2f}'.format \ (train_rf_score, test_rf_score) )

Строю ансамбли:

# голосование modelVote = VotingClassifier(estimators=[('lr', logReg), ('dt', treeTA), ('lrreg', logReg2), ('rf', best_c)], voting='hard') modelVote.fit(X_train,y_train) train_mv = modelVote.score(X_train,y_train) test_mv = modelVote.score(X_test,y_test) print('Голосование на тренировочном массиве модель работает с точностью {0}, \ а на тестовом массиве с точностью {1:6.2f}'.format \ (train_mv, test_mv) ) # бэггинг modelBag = BaggingClassifier(base_estimator=logReg, n_estimators=100, random_state=5) modelBag.fit(X_train, y_train) train_mb = modelBag.score(X_train, y_train) test_mb = modelBag.score(X_test, y_test) print('Бэггинг на тренировочном массиве модель работает с точностью {0}, \ а на тестовом массиве с точностью {1:6.2f}'.format \ (train_mb, test_mb) ) # бустинг - идентифицирую слабые модели modelAB = AdaBoostClassifier(base_estimator=logReg, n_estimators=100) #на основании высоких значений modelAB.fit(X_train, y_train) train_AB = modelAB.score(X_train, y_train) test_AB = modelAB.score(X_test, y_test) print('Бустинг (выс.знач.) на тренировочном массиве модель работает с точностью {0}, \ а на тестовом массиве с точностью {1:6.2f}'.format \ (train_AB, test_AB) ) modelGB = GradientBoostingClassifier(n_estimators=100) #на основании градиентов функции потерь modelGB.fit(X_train, y_train) train_GB = modelGB.score(X_train, y_train) test_GB = modelGB.score(X_test, y_test) print('Бустинг (град.) на тренировочном массиве модель работает с точностью {0}, \ а на тестовом массиве с точностью {1:6.2f}'.format \ (train_GB, test_GB) )

Выбираю модель:

for clf, label in zip([treeTA, logReg,logReg, best_c, modelVote, modelBag, modelAB, modelGB], ['Дерево', 'Регрессия', 'Регрессия с регуляризацией', 'Лес', 'Ансамбль голосование', 'Ансамбль bagging', 'Ансамбль ada boosting', 'Ансамбль gradient boosting']): scores = cross_val_score(clf, X_test, y_test, cv=5, scoring='accuracy') print("Accuracy: %0.2f (+/- %0.2f) [%s]" % (scores.mean(), scores.std(), label))

print('Случайный лес: ROC AUC=%.3f' % (auc))

В целом из всех наилучшие показатели у Леса и Ансамбля gradient boosting. Лес: при сравнении с Accuracy тестового и тренировочного были признаки переобучения, также у леса ROC 0.95 - хороший показатель, modelGB также имела признаки переобучения, однако у леса признаки менее выражены (1.0 -> 0.97 vs. 1-> 0.959) поэтому выбираю лес.

Исследование модели

Как влияют конкретные переменные:

Так как случайный лес -- использую feature_importances_

best_c.feature_importances_

rf = eli5.explain_weights_df(best_c, feature_names=names[0:19]) rf

# Permutation Importance permRF = PermutationImportance(best_c).fit(X_test, y_test) eli5.explain_weights_df(permRF, feature_names=names[0:19])

eli5.show_weights(permRF, feature_names=names[0:19])

Выбираю переменную wifi, LegRoom так как для данной модели они наиболее важные, чаще всего встречаются, разделение более значимо и больше разделений по feature_importances_. Также по feature_importances_ важной переменной был возраст(age), но по Permutation Importance он не показал значимости. ( + также infServ показал значимость в обеих) И мне захотелось рассмотреть переменную easyBooking.

"wifi" :

cmap2 = plt.get_cmap('Pastel1') def predict_new(X): return best_c.predict_proba(X)[:,1] wifiRF = ice(data=X_train, column='wifi', predict=predict_new) ice_plot(wifiRF, linewidth=0.5, color_by='infServ', cmap=cmap2, plot_pdp=True, pdp_kwargs={'c': 'blue', 'linewidth': 3}) plt.ylabel('Predicted') plt.xlabel('wifi');

Согласно модели, предсказание резко увеличивается при оценке wifi меньше "1" или больше "4". Также заметно небольшое увеличение качества предсказания больше "3".

"easyBooking" :

def predict_new(X): return best_c.predict_proba(X)[:,1] ebRF = ice(data=X_train, column='easyBooking', predict=predict_new) ice_plot(ebRF, linewidth=0.5, color_by='infServ', cmap=cmap2, plot_pdp=True, pdp_kwargs={'c': 'blue', 'linewidth': 3}) plt.ylabel('Predicted') plt.xlabel('easyBooking');

Согласно графику предсказание модели увеличивается ясно при оценке больше 4, то есть при "5". Небольшое увеличение заметно при оценка "4".

"LegRoom" :

def predict_new(X): return best_c.predict_proba(X)[:,1] lrRF = ice(data=X_train, column='LegRoom', predict=predict_new) ice_plot(lrRF, linewidth=0.5, color_by='infServ', cmap=cmap2, plot_pdp=True, pdp_kwargs={'c': 'blue', 'linewidth': 3}) plt.ylabel('Predicted') plt.xlabel('LegRoom');

Согласно графику у модели явное улучшение предсказания при оценке меньше "1" ( то есть "0"), далее улучшения предсказания практически не наблюдается.

Что влияет на результат конкретного примера:

explainer = shap.Explainer(best_c) shap_values = explainer.shap_values(X_test)

График важности признаков для каждого класса:

Наиболее значимым признаком - параметром оказался вайфай, легкость бронировани. Оценка вайфая имеет среднее влияние на величину выходных данных модели более 0.12.

График зависимости как изменяется wifi при изменении infServ:

# довольные shap.dependence_plot("infServ", shap_values[1], X_test, feature_names=dataTr.columns, interaction_index="wifi")

График показывает влияние, которое переменные оказывают на прогнозируемый результат. По графику видно, что связь между целью среди довольных и переменными не является линейной, а скорее она более сложная. При оценках после 3 (включительно), то есть 3, 4, 5 - связь похожа на линейную и положительную, а от 1 до 3 она более монотонна с возможными выбросами(?).

График зависимости как изменяется wifi при изменении age:

# довольные shap.dependence_plot("age", shap_values[1], X_test, feature_names=dataTr.columns, interaction_index="wifi")

На графике можно заметить, что в данной модели связь между целью среди довольных и переменными не является линейной, а скорее она более сложная. После 67 лет наблюдается положительная тенденция, а с 20 до 60 лет связь сложная и более монотонная.

Другие графики зависимости:

for feature_name in X_test.columns: shap.dependence_plot(feature_name, shap_values[1], X_test, display_features=X_test)

Анализируя графики можно сделать вывод, что при взаимосвязи foodDrink и easyBooking наиболее высокое значение shap-value наблюдается прим foodDrink = "3", при foodDrink = "3","4" или "5" очевидно увеличение shap-value с увеличением значения easyBooking, что в свою очередь, означает, что удаление этого признака означает меньшую важность признака. Обратную картину можно наблюдать при анализе графика LegRoom и onlineBoar - небольшое значение shap-value(отрицательное) для onlineBoar при её наименьших оценках - "1","2" или "3". Также практически на всех графиках наблюдаются выбросы. Делать выводы о важности тех или иных признаков в модели сложно, так как shap-values зависят не только от модели, но и от распределения данных, при этом даже признаки, которые никак не используются моделью, могут иметь ненулевые shap-values.

Для случая "0" :

explainer1 = TreeExplainer(best_c) sv = explainer1(X) exp = Explanation(sv.values[:,:,1], sv.base_values[:,1], data=dataTr.values, feature_names=dataTr.columns) idx = 0 waterfall(exp[idx], max_display=20)

Можно заметить, что вайфай = "1" снижает вероятность. Такой же тренд можно заметить при оценке easyBooking или onlineBoar или LegRoom "0". Увеличение можно заметить на 7 параметрах. Но не большое, меньше 0.01. Среднее прогнозируемое значение 0.108, а минимальное = 0.

SHAP-значения для предсказания 1 (довольный) для примера i

shap_values[1][idx]

Для случая "15" :

idx = 15 waterfall(exp[idx], max_display=20)

При значении оценка вайфая в "3" снижает прогнозируемую вероятность того, что этот человек будет доволен. Такой же тренд можно заметить если бы оценка easyBooking или LegRoom также была равна "3". Более явное улучшение предсказания можно заметить при расстоянии полёта в 236 км. Среднее прогнозируемое значение 0.108, минимальное 0.026

shap_values[1][idx]

Для случая "120" :

idx = 120 waterfall(exp[idx], max_display=20)

Тут уже интереснее. Среднее прогнозируемое значение также 0.108(как и будет везде), минимальное 0.105. При значении оценка вайфая в "3" снижает прогнозируемую вероятность того, что этот человек будет доволен на -0.06. В отличие от прошлого примера для №15 при easyBooking в "3" предсказание увеличивается +0.03, такое же увеличение предсказания можно заметить при оценке onlineBoar "0". В целом 8 значений уменьшают предсказание, а 11 увеличивают.

explainer1.expected_value[1]

best_c.predict_proba(X_train)[:,1].mean()

Основной вопрос

Оглавление

Данные и библиотеки

Обработка данных

Построение моделей для предсказания

Дерево

Случайный лес

Ансамбли

Выбираю модель:

Исследование модели

Как влияют конкретные переменные:

"LegRoom" :

Что влияет на результат конкретного примера:

Для случая "15" :

Для случая "120" :