1.0 MiB
Постановка задачи¶
Крупная авиакомпания хочет повысить эффективность своей работы за счет сокращения времени задержки прибытия внутренних рейсов в будущем году.
Сокращение времени задержки позволит снизить расходы компании и повысить ее привлекательность для клиентов.
Для анализа доступны данные о полетах за текущий год.
Задачи:
- определение факторов, влияющих на задержку рейсов в текущем году;
- прогнозировать задержку рейса для принятия соответствующих мер по сокращению задержки (при возможности).
Эксперты сообщили, что на задержку рейсов влияют разные типы факторов:
- Некоторые факторы могут быть устранены или возможно снижение влияния данных факторов на задержку рейса (загруженность аэропорта, влияние типа аэропорта на задержу, проблемы с техобслуживанием или экипажем, чистка воздушного судна, загрузка багажа, заправка топливом).
- Некоторые факторы выходят из под контроля авиакомпании и не могут быть устранены (погода, проблемы с судном, проблемы обеспечения безопасности).
- Задержка судна в предыдущем рейсе не учитывается.
- Опозданием считается задержка более 15 минут.
Загрузка и распаковка данных¶
from urllib.request import urlretrieve
from zipfile import ZipFile
ds_url = "https://github.com/PacktPublishing/Interpretable-Machine-Learning-with-Python/raw/master/datasets/aa-domestic-delays-2018.csv.zip"
ds_zip_filename = "data/aa-domestic-delays-2018.csv.zip"
urlretrieve(ds_url, ds_zip_filename)
with ZipFile(ds_zip_filename, "r") as zObject:
zObject.extractall(path="data")
Загрузка данных в Dataframe¶
import pandas as pd
orig_df = pd.read_csv("data/aa-domestic-delays-2018.csv")
orig_df.info()
Признаки¶
Общие признаки:
- FL_NUM – номер рейса.
- ORIGIN – код аэропорта отправления.
- DEST – код аэропорта назначения.
Признаки вылета:
- PLANNED_DEP_DATETIME – запланированные дата и время рейса.
- CRS_DEP_TIME – запланированное время отправления.
- DEP_TIME – фактическое время отправления.
- DEP_AFPH – количество фактических вылетов в час. Данный признак показывает степень загруженности аэропорта отправления во время взлета.
- DEP_RFPH – относительное количество вылетов в час. Данный признак показывает степень загруженности аэропорта отправления во время взлета.
- TAXI_OUT – время, прошедшее между выездом самолета на взлетную полосу аэропорта отправления и отрывом колес самолета от земли.
- WHEELS_OFF – момент времени, когда колеса самолета отрываются от земли.
Полетные признаки:
- CRS_ELAPSED_TIME – запланированное время, необходимое для полета.
- PCT_ELAPSED_TIME – отношение фактического времени полета к запланированному времени полета для измерения относительной скорости самолета.
- DISTANCE – расстояние между двумя аэропортами.
Признаки прибытия:
- CRS_ARR_TIME – запланированное время прибытия.
- ARR_AFPH – число фактических прибытия в час. Этот признак показывает степень загруженности аэропорта назначения во время посадки самолета.
- ARR_RFPH – относительное число прибытий в час. Этот признак показывает степень загруженности аэропорта назначения во время посадки самолета.
- DEP_DELAY – суммарная задержка отправления в минутах.
- ARR_DELAY – суммарная задержка прибытия в минутах, которая включает следующие параметры.
- CARRIER_DELAY (целевая переменная) – задержка в минутах, вызванная обстоятельствами, находящимися в пределах контроля авиакомпании.
- WEATHER_DELAY – задержка в минутах, вызванная метеорологическими условиями.
- NAS_DELAY – задержка в минутах, предусмотренная национальной авиационной системой.
- SECURITY_DELAY – задержка в минутах, вызванная процедурами обеспечения безопасности.
- LATE_AIRCRAFT_DELAY – задержка в минутах, вызванная предыдущим рейсом того же самолета, который прибыл с опозданием.
Подготовка данных и конструирование признаков¶
- Преобразование даты (PLANNED_DEP_DATETIME) из строки в datetime.
- Получение месяца и дня недели вылета из даты для учета сезонности и особенностей дня недели (DEP_MONTH и DEP_DOW).
- Удаление столбца с датой (PLANNED_DEP_DATETIME).
- Определение признака хаба для аэропортов вылета и назначения (ORIGIN_HUB и DEST_HUB).
- Удаление лишних столбцов (FL_NUM, ORIGIN, DEST).
- Удаление столбца с общим временем задержки прибытия, так как данные значения будут иметь сильное влияние на результат (ARR_DELAY).
df = orig_df.copy()
# Преобразование даты из строки в datetime
df["PLANNED_DEP_DATETIME"] = pd.to_datetime(df["PLANNED_DEP_DATETIME"])
# Получение месяца и дня недели вылета из даты для учета сезонности и особенностей дня недели
df["DEP_MONTH"] = df["PLANNED_DEP_DATETIME"].dt.month
df["DEP_DOW"] = df["PLANNED_DEP_DATETIME"].dt.dayofweek
# Удаление столбца с датой
df = df.drop(["PLANNED_DEP_DATETIME"], axis=1)
# Список аэропортов-хабов
hubs = ["CLT", "ORD", "DFW", "LAX", "MIA", "JFK", "LGA", "PHL", "PHX", "DCA"]
# Определение признака хаба для аэропортов вылета и назначения
is_origin_hub = df["ORIGIN"].isin(hubs)
is_dest_hub = df["DEST"].isin(hubs)
# Установка признака хаба для данных
df["ORIGIN_HUB"] = 0
df.loc[is_origin_hub, "ORIGIN_HUB"] = 1
df["DEST_HUB"] = 0
df.loc[is_dest_hub, "DEST_HUB"] = 1
# Удаление лишних столбцов
df = df.drop(["FL_NUM", "ORIGIN", "DEST"], axis=1)
# Удаление столбца с общим временем задержки прибытия, так как данные значения будут иметь сильное влияние на результат
df = df.drop(["ARR_DELAY"], axis=1)
df
Формирование тестовой и обучающей выборок данных¶
- Задание фиксированного случайного состояния для воспроизводимости результатов (rand = 9).
- Выделение признака, который модель должна предсказать (CARRIER_DELAY).
- Формирование множества признаков, на основе которых модель будет обучаться (удаление столбца CARRIER_DELAY).
- Формирование выборок (85 % обучающая, 15 % тестовая).
- Создание классов для классификаторов в виде двоичных меток (опоздание свыше 15 минут – 1, иначе – 0).
from sklearn.model_selection import train_test_split
# Задание фиксированного случайного состояния для воспроизводимости результатов
rand = 9
# Выделение признака, который модель должна предсказать
y = df["CARRIER_DELAY"]
# Формирование множества признаков, на основе которых модель будет обучаться (удаление столбца с y)
X = df.drop(["CARRIER_DELAY"], axis=1).copy()
X_train, X_test, y_train_reg, y_test_reg = train_test_split(
X, y, test_size=0.15, random_state=rand
)
# Создание классов для классификаторов в виде двоичных меток (опоздание свыше 15 минут - 1, иначе - 0)
y_train_class = y_train_reg.apply(lambda x: 1 if x > 15 else 0)
y_test_class = y_test_reg.apply(lambda x: 1 if x > 15 else 0)
display(X_train)
display(y_train_reg)
display(y_train_class)
X_train.info()
Определение линейной корреляции признаков с целевым признаком с помощью корреляции Пирсона¶
Существует сильная прямая корреляция между целевым признаком CARRIER_DELAY и признаком DEP_DELAY (70 %).
Но это не значит, что в данных отсутствуют нелинейные зависимости, или что совокупность признаков не оказывает влияние на значение целевого признака.
corr = df.corr()
abs(corr["CARRIER_DELAY"]).sort_values(ascending=False)
Использование регрессионных моделей для предсказания задержки рейса¶
- Линейная регрессия (linear).
- Полиномиальная регрессия (linear_poly).
- Полиномиальная регрессия без квадратичных членов (linear_interact).
- Гребневая регрессия (ridge).
- Дерево решений (decision_tree).
- k ближайших соседей (kNN).
- Случайный лес (random_forest).
- Многослойный персептрон (mlp).
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn import linear_model, tree, neighbors, ensemble, neural_network
reg_models = {
# Линейные модели
"linear": {"model": linear_model.LinearRegression(n_jobs=-1)},
"linear_poly": {
"model": make_pipeline(
PolynomialFeatures(degree=2, interaction_only=False),
linear_model.LinearRegression(fit_intercept=False, n_jobs=-1),
memory=None
)
},
"linear_interact": {
"model": make_pipeline(
PolynomialFeatures(degree=2, interaction_only=True),
linear_model.LinearRegression(fit_intercept=False, n_jobs=-1),
memory=None
)
},
"ridge": {"model": linear_model.RidgeCV()},
# Деревья
"decision_tree": {
"model": tree.DecisionTreeRegressor(max_depth=7, random_state=rand)
},
# Ближайшие соседи
"knn": {"model": neighbors.KNeighborsRegressor(n_neighbors=7, n_jobs=-1)},
# Ансамблевые методы
"random_forest": {
"model": ensemble.RandomForestRegressor(
max_depth=7, random_state=rand, n_jobs=-1
)
},
# Нейронные сети
"mlp": {
"model": neural_network.MLPRegressor(
hidden_layer_sizes=(21,),
max_iter=500,
early_stopping=True,
random_state=rand,
)
},
}
Обучение и оценка регрессионных моделей¶
Метрики:
- RMSE – корень из средней квадратичной ошибки.
- MAE – средняя абсолютная ошибка.
- R-квадрат – коэффициент детерминации.
import math
from sklearn import metrics
for model_name in reg_models.keys():
print(f'Model: {model_name}')
fitted_model = reg_models[model_name]["model"].fit(
X_train.values, y_train_reg.to_numpy().ravel()
)
y_train_pred = fitted_model.predict(X_train.values)
y_test_pred = fitted_model.predict(X_test.values)
reg_models[model_name]["fitted"] = fitted_model
reg_models[model_name]["preds"] = y_test_pred
reg_models[model_name]["RMSE_train"] = math.sqrt(
metrics.mean_squared_error(y_train_reg, y_train_pred)
)
reg_models[model_name]["RMSE_test"] = math.sqrt(
metrics.mean_squared_error(y_test_reg, y_test_pred)
)
reg_models[model_name]["MAE_train"] = metrics.mean_absolute_error(y_train_reg, y_train_pred)
reg_models[model_name]["MAE_test"] = metrics.mean_absolute_error(
y_test_reg, y_test_pred
)
reg_models[model_name]["R2_test"] = metrics.r2_score(y_test_reg, y_test_pred)
Вывод оценки в виде таблицы¶
reg_metrics = pd.DataFrame.from_dict(reg_models, "index")[
["RMSE_train", "RMSE_test", "MAE_train", "MAE_test", "R2_test"]
]
reg_metrics.sort_values(by="RMSE_test").style.background_gradient(
cmap="viridis", low=1, high=0.3, subset=["RMSE_train", "RMSE_test"]
).background_gradient(
cmap="plasma", low=0.3, high=1, subset=["MAE_train", "MAE_test", "R2_test"]
)
Наилучший результат показали модели с типом черный ящик и аквариумные модели.
Определение сбалансированности выборки для классификации¶
y_train_class[y_train_class == 1].shape[0] / y_train_class.shape[0]
Выборка данных имеет сильный дисбаланс в сторону отрицательного класса (y=0).
Для решения данной проблемы используется гиперпараметр class_weight="balanced" во многих моделях классификации.
Использование классификаторов для предсказания задержки рейса¶
- Логистическая регрессия (logistic).
- Гребневая классификация (ridge).
- Дерево решений (decision_tree).
- k ближайших соседей (knn).
- Наивный Байес (naive_bayes).
- Градиентно-бустированные деревья (gradient_boosting).
- Случайный лес (random_forest).
- Многослойный персептрон (mlp).
Сбалансированность выборки: 0,061. Использование гиперпараметра class_weight="balanced" позволяет адаптировать модели к данной особенности данных.
from sklearn import naive_bayes
class_models = {
# Линейные модели
"logistic": {"model": linear_model.LogisticRegression()},
"ridge": {
"model": linear_model.LogisticRegression(penalty="l2", class_weight="balanced")
},
# Дерево
"decision_tree": {
"model": tree.DecisionTreeClassifier(max_depth=7, random_state=rand)
},
# Ближайшие соседи
"knn": {"model": neighbors.KNeighborsClassifier(n_neighbors=7)},
# Наивный Байес
"naive_bayes": {"model": naive_bayes.GaussianNB()},
# Ансамблевые методы
"gradient_boosting": {
"model": ensemble.GradientBoostingClassifier(n_estimators=210)
},
"random_forest": {
"model": ensemble.RandomForestClassifier(
max_depth=11, class_weight="balanced", random_state=rand
)
},
# Нейронные сети
"mlp": {
"model": make_pipeline(
StandardScaler(),
neural_network.MLPClassifier(
hidden_layer_sizes=(7,),
max_iter=500,
early_stopping=True,
random_state=rand,
),
memory=None
)
},
}
Обучение и оценка классификаторов¶
Метрики:
- Accuracy – верность/аккуратность.
- Recall – полнота.
- ROC AUC – площадь под ROC-кривой.
- F1-мера.
- MCC – коэффициент корреляции Мэтьюса.
import numpy as np
for model_name in class_models.keys():
print(f"Model: {model_name}")
fitted_model = class_models[model_name]["model"].fit(
X_train.values,
y_train_class.to_numpy().ravel(),
)
y_train_pred = fitted_model.predict(X_train.values)
y_test_prob = fitted_model.predict_proba(X_test.values)[:, 1]
y_test_pred = fitted_model.predict(X_test.values)
class_models[model_name]["fitted"] = fitted_model
class_models[model_name]["probs"] = y_test_prob
class_models[model_name]["preds"] = y_test_pred
class_models[model_name]["Accuracy_train"] = metrics.accuracy_score(
y_train_class, y_train_pred
)
class_models[model_name]["Accuracy_test"] = metrics.accuracy_score(
y_test_class, y_test_pred
)
class_models[model_name]["Recall_train"] = metrics.recall_score(
y_train_class, y_train_pred
)
class_models[model_name]["Recall_test"] = metrics.recall_score(
y_test_class, y_test_pred
)
class_models[model_name]["ROC_AUC_test"] = metrics.roc_auc_score(
y_test_class, y_test_prob
)
class_models[model_name]["F1_test"] = metrics.f1_score(y_test_class, y_test_pred)
class_models[model_name]["MCC_test"] = metrics.matthews_corrcoef(
y_test_class, y_test_pred
)
Вывод оценки в виде таблицы¶
class_metrics = pd.DataFrame.from_dict(class_models, "index")[
[
"Accuracy_train",
"Accuracy_test",
"Recall_train",
"Recall_test",
"ROC_AUC_test",
"F1_test",
"MCC_test",
]
]
class_metrics.sort_values(by="ROC_AUC_test", ascending=False).style.background_gradient(
cmap="plasma", low=0.3, high=1, subset=["Accuracy_train", "Accuracy_test"]
).background_gradient(
cmap="viridis",
low=1,
high=0.3,
subset=["Recall_train", "Recall_test", "ROC_AUC_test", "F1_test", "MCC_test"],
)
Использование различных меток позволяют оценить качество моделей классификации с точки зрения различных свойств (точность и полнота), что особенно важно при несбалансированной выборке.
Интерпретация результатов для моделей типа белый ящик¶
К моделям типа белого ящика в нашем случае относятся:
- Линейные модели: разные виды регрессии.
- Деревья решений.
- Метод ближайших соседей.
- Наивный Байес.
К аквариумным моделям (серый ящик) относятся ансамблевые модели: случайный лес, градиентный бустинг.
К моделям типа черного ящика относятся искусственные нейронные сети (MLP).
Линейные модели¶
Линейная регрессия
Коэффициенты и пересечение модели
coefs_lm = reg_models["linear"]["fitted"].coef_
intercept_lm = reg_models["linear"]["fitted"].intercept_
print("coefficients:\t%s" % coefs_lm)
print("intercept:\t%s" % intercept_lm)
print(
"y = %0.2f + %0.4fX1 + %0.4fX2 + %0.3fX3 + ..."
% (intercept_lm, coefs_lm[0], coefs_lm[1], coefs_lm[2])
)
Вывод признаков и соответствующих им коэффициентам модели регрессии
coef_df = pd.DataFrame({"feature": X_train.columns.values.tolist(), "coef": coefs_lm})
display(coef_df)
Интерпретация признаков:
- Числовой признак (DEP_DELAY). Увеличение значения признака на 1 единицу увеличивает целевой признак на значение коэффициента в соответствующих единицах.
- Бинарный признак (ORIGIN_HUB). Если значение отрицательное, то истинное значение признака уменьшает целевой признак на значение коэффициента.
- Категориальные признаки (DEP_MONTH). При условии кодирования категориальные признаки позволяют отслеживать влияние отдельных признаков на целевой признак.
y_train_class[y_train_class == 1].shape[0] / y_train_class.shape[0]
Получение сведений о значимости признаков средствами библиотеки statsmodels
import statsmodels.api as sm
linreg_mdl = sm.OLS(y_train_reg, sm.add_constant(X_train))
linreg_mdl = linreg_mdl.fit()
summary_df = linreg_mdl.summary2().tables[1]
summary_df = (
summary_df.drop(["const"]).reset_index().rename(columns={"index": "feature"})
)
summary_df["t_abs"] = abs(summary_df["t"])
summary_df.sort_values(by="t_abs", ascending=False).style.background_gradient(
cmap="plasma_r", low=0, high=0.1, subset=["P>|t|"]
).background_gradient(cmap="plasma_r", low=0, high=0.1, subset=["t_abs"])
Важными являются p-значение (P>|t|) и t-статистика (t_abs).
Из результатов видно, что самыми значимыми признаками модели являются:
- DEP_DELAY;
- LATE_AIRCRAFT_DELAY;
- WEATHER_DELAY;
- NAS_DELAY.
Гребневая регрессия
Интерпретация гребневой регрессии производится аналогично интерпретации линейной регрессии
coefs_ridge = reg_models["ridge"]["fitted"].coef_
coef_ridge_df = pd.DataFrame(
{
"feature": X_train.columns.values.tolist(),
"coef_linear": coefs_lm,
"coef_ridge": coefs_ridge,
}
)
coef_ridge_df.style.background_gradient(cmap="viridis_r", low=0.3, high=0.2, axis=1)
Полиномиальная регрессия
Для интерпретации и оценки важности признаков полиномиальной регрессии рекомендуется использовать возможности библиотеки statsmodels, так как модель содержит 253 линейных члена.
display(
reg_models["linear_poly"]["fitted"].get_params()["linearregression"].coef_.shape[0]
)
display(
reg_models["linear_interact"]["fitted"]
.get_params()["linearregression"]
.coef_.shape[0]
)
Логистическая регрессия
Интерпретация логистической регрессии выполняется аналогично интерпретации линейной регрессии
Способ оценки важности признаков модели на основе логистической регрессии рассмотрен в предыдущей лекции
coefs_log = class_models["logistic"]["fitted"].coef_
intercept_log = class_models["logistic"]["fitted"].intercept_
print("coefficients:\t%s" % coefs_log)
print("intercept:\t%s" % intercept_log)
stdv = np.std(X_train, 0)
abs(
coefs_log.reshape(
21,
)
* stdv
).sort_values(ascending=False)
Дерево решений
В библиотеке scikit-learn для построения дерева решений используется алгоритм CART https://ru.wikipedia.org/wiki/CART_(алгоритм)
Оценочная функция, используемая алгоритмом CART, базируется на идее уменьшения неопределенности (неоднородности) в узле
Степень неопределенности вычисляется на основе индекса Gini-Simpson (Gini impurity)
Значение данного индекса позволяет рекурсивно расщеплять узел дерева в местах, где две ветви отличаются максимально
Визуализация дерева решений
Каждый узел содержит:
- Условие
- Значение индекса Джини-Симпсона
- Количество наблюдений, которые поступили на вход
- Количество наблюдения для левого и правого выходов
Если в узле отсутствует условие, то такой узел является листом дерева и содержит конечное решение
Путь от корня дерева до листа описывает некоторое правило дерева решений
import matplotlib.pyplot as plt
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(16, 8), dpi=600)
tree.plot_tree(
class_models["decision_tree"]["fitted"],
feature_names=X_train.columns.values.tolist(),
filled=True,
max_depth=2,
)
fig.show()
В левом узле 629167 наблюдений отнесены к классу 0 (без опозданий)
Для данного листа можно сформулировать следующее правило: Если DEP_DELAY <= 15.5 Тогда 0
Аналогичное правило можно наблюдать в множестве правил, который можно извлечь из дерева решений:\ IF DEP_DELAY <= 20.5 AND DEP_DELAY <= 15.5 THEN class: 0
Вывод правил из модели дерева решений
text_tree = tree.export_text(
class_models["decision_tree"]["fitted"],
feature_names=X_train.columns.values.tolist(),
)
print(text_tree)
Библиотека позволяет получить важность признаков модели дерева решений
dt_imp_df = pd.DataFrame(
{
"feature": X_train.columns.values.tolist(),
"importance": class_models["decision_tree"]["fitted"].feature_importances_,
}
).sort_values(by="importance", ascending=False)
dt_imp_df
k ближайших соседей
Интерпретировать модель на основе алгоритма k ближайших соседей можно только локально: на уровне отдельных наблюдений
Например, используем наблюдение с id = 721043 из тестовой выборки
Значения целевого признака для случая 721043:
- Фактическое ($y$): 1
- Предсказанное ($\hat{y}$): 0
display(X_test.loc[721043, :])
display(y_test_class[721043])
display(class_models["knn"]["preds"][X_test.index.get_loc(721043)])
Библиотека scikit-learn позволяет получить информацию об указанном количестве соседей (7) для некоторого элемента из тестовой выборки (id = 721043)
Первый элемент кортежа – расстояние от соседей в обучающей выборке до текущей точки в тестовой выборке
Второй элемент – идентификаторы соседей в тестовой выборке
class_models["knn"]["fitted"].kneighbors(
X_test.loc[721043, :].values.reshape(1, 21), 7
)
Данные о целевом признаке CARRIER_DELAY соседей для случая 721043 из обучающей выборки
y_train_class.iloc[[105172, 571912, 73409, 89450, 77474, 705972, 706911]]
Можно сделать вывод, что значение $\hat{y}$ для случая 721043 было определено как мода от значений целевой переменной соседей.
Мода – одно или несколько значений во множестве наблюдений, которое встречается наиболее часто.
Наивный байесовский классификатор
Алгоритм обучения модели на основе наивного байесовского классификатора вычисляет априорные вероятности принадлежности наблюдения некоторому классу
$\hat{y} = P(y=1|X) = P(y=1) \prod^{n}_{i=1} P(x_{i} | y=1)$
$P(x_{i} | y=1) = \frac{ 1 }{ \sqrt{2 \pi \sigma_{i}^{2}}} e^{\frac{ (x_{i} \theta_{i})^2 }{ 2 \sigma^{2}_{i} }}$
Библиотека scikit-learn позволяет получить вероятности $P(y=0|X)$ и $P(y=1|X)$
class_models["naive_bayes"]["fitted"].class_prior_
Вероятность принадлежности классу 1 соответствует показателю сбалансированности выборки
Для получения значения дисперсии $\sigma_{i}$ и среднего значения $\theta_{i}$ признаков в одном классе по сравнению с другим используются атрибуты var и theta
Дисперсия
class_models["naive_bayes"]["fitted"].var_
Средние значения
class_models["naive_bayes"]["fitted"].theta_
Эти значения можно использовать для глобальной и локальной интерпретации модели.
Результат¶
Все алгоритмы обучения показали высокие значения метрик качества, но выявили некорректные закономерности.
Как видно из результатов интерпретации моделей, значение целевого признака CARRIER_DELAY сильно зависит от признака DEP_DELAY.
Зная суммарную задержку отправления в минутах (DEP_DELAY) можно с высокой степенью качества спрогнозировать задержку в минутах, вызванную обстоятельствами, находящимися в пределах контроля авиакомпании (CARRIER_DELAY).
Если в производственной среде значение признака DEP_DELAY недоступно или станет доступным после применения модели, то модель будет иметь низкое качество.