Данный материал описывает участие в соревновании "DengAI: Predicting Disease Spread" на портале Drivedata.com. Для исследования предоставляется набор данных о выявленных случаях заражения в двух городах: Сан-Хуан, Пуэрто-Рико и Икитос, Перу в течение 18 и 10 лет соответственно. Население городов составляет 395 и 406 тыс. человек, т.е. количества случаев заболевания можно не нормализовывать.
Целью исследования является исследование зависимостей всплесков частоты заболевания от погодных условий. Для проверки требуется предсказать количество выявленных заболеваний в двух городах за период в 5 и 3 года соответственно.
Информация о заболевании из Википедии:
Данный вирус переносится комарами рода Aedes, чаще всего вида Aedes aegypti. По этой причине динамика переноса лихорадки связана с климатическими условиями, такими как температура и количество осадков.
Информация о переносчике заболевания из Википедии:
Данные представляют собой объединение информации из различных источников: Центр по Контролю и Предотврщению Заболеваний США(CDC USA) и Национальное управление океанических и атмосферных исследований США(NOAA).
Глубокое понимание зависимости болезни от погодных условий должно помочь в предотвращении эпидемий.
total_cases - количество выявленных за соответствующую неделю заболеваний
Тестовая выборка - последовательная по времени выборка из будущего по отношению к доступным для обучения данным.
city – название города: sj - Сан-Хуан и iq - Икитос
year – год
weekofyear – номер недели в текущем году
week_start_date – дата в формате гггг-мм-дд
station_max_temp_c – максимальная температура воздуха в течение суток в градусах Цельсия
station_min_temp_c – минимальная температура воздуха в течение суток в градусах Цельсия
station_avg_temp_c – средняя температура воздуха в течение суток в градусах Цельсия
station_precip_mm – размер осадков, измеренный датчиками на поверхности
station_diur_temp_rng_c – дневные колебания температуры, т.е. изменение температуры в период от восхода до заката, в градусах Цельсия
precipitation_amt_mm – замер осадков, но произведенный с помощью спутника
reanalysis_sat_precip_amt_mm – количество осадков в мм
reanalysis_dew_point_temp_k – средняя температура росы в градусах Кельвина
reanalysis_air_temp_k – средняя температура воздуха в течение суток на уровне 2м от земли в градусах Кельвина
reanalysis_relative_humidity_percent – относительная влажность
reanalysis_specific_humidity_g_per_kg – абсолютная влажность в г воды на кг воздуха
reanalysis_precip_amt_kg_per_m2 – количество осадков в кг на м2
reanalysis_max_air_temp_k – максимальная температура воздуха в течение суток в градусах Кельвина
reanalysis_min_air_temp_k – минимальная температура воздуха в течение суток в градусах Кельвина
reanalysis_avg_temp_k – средняя температура воздуха в течение суток на уровне земл в градусах Кельвина
reanalysis_tdtr_k – дневные колебания температуры в градусах Кельвина
ndvi_se – юго-восточной части города
ndvi_sw – юго-западной части города
ndvi_ne – северо-восточной части города
ndvi_nw – северо-западной части города
# Импортируем необходимые библиотеки для обработки и визуализации данных
from hyperopt import fmin, tpe, rand, hp, STATUS_OK, Trials
from mlxtend.feature_selection import SequentialFeatureSelector as SFS
from sklearn.linear_model import Lasso
from sklearn.metrics import mean_absolute_error
from sklearn.model_selection import TimeSeriesSplit, cross_val_score, GridSearchCV, learning_curve
from sklearn.preprocessing import StandardScaler
from warnings import filterwarnings
filterwarnings('ignore')
import math
import matplotlib.pylab as pylab
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import statsmodels.api as sm
import statsmodels.formula.api as smf
import xgboost as xgb
plt.style.use('ggplot')
%matplotlib inline
Данные предоставлены в виде 3 файлов: признаки для обучения, целевой признак и признаки для предсказания.
features = pd.read_csv('../../data/dengue_features_train.csv', infer_datetime_format=True)
labels = pd.read_csv('../../data/dengue_labels_train.csv', infer_datetime_format=True)
test = pd.read_csv('../../data/dengue_features_test.csv', infer_datetime_format=True)
train_df = pd.merge(features, labels, how='outer', on=labels.columns.tolist()[:-1])
test_df = test.copy()
train_df.iloc[:, :15].head()
train_df.iloc[:, 15:].head()
Два города сильно различаются по климату и качеству сбора информации, поэтому будем рассматривать признаки отдельно по городам
Проверим, все ли данные в наличии:
train_df[train_df.city=='sj'].isnull().sum()
train_df[train_df.city=='iq'].isnull().sum()
test_df[test_df.city=='sj'].isnull().sum()
test_df[test_df.city=='iq'].isnull().sum()
Взглянем на переходы из года в год:
train_df[(train_df.weekofyear==1)|(train_df.weekofyear==52)|(train_df.weekofyear==53)].iloc[:, :15]
train_df.year = np.where((train_df.week_start_date.apply(lambda x: x[5:]) == '01-01')&(train_df.weekofyear != 1),
train_df.year-1, train_df.year)
test_df.year = np.where((test_df.week_start_date.apply(lambda x: x[5:]) == '01-01')&(test_df.weekofyear != 1),
test_df.year-1, test_df.year)
train_df = train_df.sort_values(by='week_start_date')
train_df.week_start_date = train_df.week_start_date.apply(pd.to_datetime)
train_df.set_index('week_start_date', inplace=True)
test_df.set_index('week_start_date', inplace=True)
Теперь посмотрим на распределение признаков:
cols = train_df.columns[1:3].tolist() + train_df.columns[3:25].tolist()
train_df[train_df.city=='sj'][cols[:13]].describe()
train_df[train_df.city=='iq'][cols[:13]].describe()
train_df[train_df.city=='sj'][cols[13:]].describe()
train_df[train_df.city=='iq'][cols[13:]].describe()
Найдем пары признаков с максимальной корреляцией:
def mostcorr(dataframe, numtoreport):
cormatrix = dataframe.corr()
cormatrix *= np.tri(*cormatrix.values.shape, k=-1).T
cormatrix = cormatrix.stack()
cormatrix = cormatrix.reindex(cormatrix.abs().sort_values(ascending=False).index).reset_index()
cormatrix.columns = ["Признак №1", "Признак №2", "Корреляция"]
return cormatrix.head(numtoreport)
mostcorr(train_df[train_df.city=='sj'],10)
mostcorr(train_df[train_df.city=='iq'],10)
del train_df['reanalysis_sat_precip_amt_mm']
del train_df['reanalysis_dew_point_temp_k']
del train_df['reanalysis_avg_temp_k']
del test_df['reanalysis_sat_precip_amt_mm']
del test_df['reanalysis_dew_point_temp_k']
del test_df['reanalysis_avg_temp_k']
Посмотрим на связь признаков с целевой переменной:
sns.set(font_scale = 1.5)
(abs(train_df[(train_df.city=='sj')].corr())
.total_cases
.drop('total_cases')
.sort_values()
.plot
.barh())
sns.set(font_scale = 1.5)
(abs(train_df[(train_df.city=='iq')].corr())
.total_cases
.drop('total_cases')
.sort_values()
.plot
.barh())
Часть признаков присутствует в обоих группах. Это означает, что абсолютная влажность и оба замера минимальной температуры воздуха могут быть наиболее полезны в предсказании целевого признака. Интересно также присутствие признака 'year'.
Рассмотрим график относительной частотности количества случаев заболевания для двух городов. Оба распределения явно относятся к отрицательным биномиальным.
fig = sns.FacetGrid(train_df, hue='city', aspect=4, size=4)
fig.map(sns.kdeplot,'total_cases',shade=True)
max_x = train_df.total_cases.max()
min_x = train_df.total_cases.min()
fig.set(xlim=(min_x,max_x))
fig.set(ylim=(0, 0.09))
fig.add_legend()
fig.fig.suptitle("Количество случаев")
Несмотря на визуальную разницу распределения отличаются только масштабом:
train_df['total_cases_temp'] = train_df.total_cases.where(train_df.city=='sj', train_df.total_cases*3.3)
fig = sns.FacetGrid(train_df, hue='city', aspect=4, size=4)
fig.map(sns.kdeplot,'total_cases_temp',shade=True)
max_x = train_df.total_cases.max()
min_x = train_df.total_cases.min()
fig.set(xlim=(min_x,max_x))
fig.add_legend()
fig.fig.suptitle("Количество случаев масштабированное")
del train_df['total_cases_temp']
Вычислим параметры распределений для целевой переменной: скошенность, а также p и r.
sj_std = train_df[train_df.city=='sj'].total_cases.std()
sj_mean = train_df[train_df.city=='sj'].total_cases.mean()
sj_median = train_df[train_df.city=='sj'].total_cases.median()
iq_std = train_df[train_df.city=='iq'].total_cases.std()
iq_mean = train_df[train_df.city=='iq'].total_cases.mean()
iq_median = train_df[train_df.city=='iq'].total_cases.median()
#Скошенность
sj_skew = (sj_mean - sj_median) / sj_std
iq_skew = (iq_mean - iq_median) / iq_std
#p
sj_p = 1 - sj_mean / sj_std**2
iq_p = 1 - iq_mean / iq_std**2
#r
sj_r = sj_mean**2 / sj_std**2 / sj_p
iq_r = iq_mean**2 / iq_std**2 / iq_p
print('Скошенность распределения целевой переменной для Сан-Хуана равна:%0.3f' % sj_skew)
print('Скошенность распределения целевой переменной для Икитос равна:%0.3f' % iq_skew)
print('')
print('Параметры p и r распределения целевой переменной для Сан-Хуана равны:%0.3f' % sj_p + ' и %0.3f' % sj_r)
print('Параметры p и r распределения целевой переменной для Икитос равны:%0.3f' % iq_p + ' и %0.3f' % iq_r)
Это еще одно доказательство схожести распределений целевой переменной для обоих городов
Посмотрим на историческое распределение случаев заболевания в обоих городах:
train_df[train_df.city=='sj'].plot(kind='line', x=['year', 'weekofyear'], y='total_cases',
figsize=[15,7],title='Распределение случаев болезни в Сан-Хуане по времени')
train_df[train_df.city=='iq'].plot(kind='line', x=['year', 'weekofyear'], y='total_cases',
figsize=[15,7],title='Распределение случаев болезни в Икитос по времени')
Данные по Икитос шумнее, но в целом видно, что всплески заболевания имеют схожую продолжительность.
Теперь посмотрим на распределение случаев заболевания понедельно:
fig = sns.FacetGrid(train_df, hue='city', aspect=4, size=4)
fig.map(sns.pointplot,'weekofyear','total_cases')
max_x = train_df.weekofyear.max()
min_x = train_df.weekofyear.min()
fig.set(xlim=(min_x,max_x))
fig.set(ylim=(0, 130))
fig.add_legend()
fig.fig.suptitle("Распределение случаев болезни в течение года")
Наблюдается интересная сезонность с минимумом в районе 13-19 недель и максимумом - 40-46. Более точные выводы невозможно сделать ввиду малого размера выборки.
Для приближения распределения целевого признака к нормальной форме используем операцию логарифмирования:
train_df['total_cases_adj'] = train_df.total_cases.apply(math.log1p)
Принимая во внимание срок инкубации и жизненный цикл распространителя, попробуем рассмотреть влияние текущих погодных условий на уровень заболевания в будущем на 1, 2, 3 и 4 недели. Для этого создадим новые переменные из целевой сдвигом на соответствующее число недель.
sj_df = train_df[train_df.city=='sj'].copy()
iq_df = train_df[train_df.city=='iq'].copy()
for i in range(1, 5):
sj_df['total_cases'+str(i)] = sj_df.total_cases.shift(i)
iq_df['total_cases'+str(i)] = iq_df.total_cases.shift(i)
train_df = pd.concat((sj_df, iq_df))
Так как города имеют разный климат, взглянем на взаимодействие признаков для каждого из них отдельно:
f, ax = plt.subplots(figsize=(30, 30))
corr = train_df[(train_df.city=='sj')].corr()
sns.heatmap(corr, mask=np.zeros_like(corr, dtype=np.bool),
cmap=sns.diverging_palette(220, 10, as_cmap=True), annot=True, square=True, center=0, ax=ax)
f, ax = plt.subplots(figsize=(30, 30))
corr = train_df[train_df.city=='iq'].corr()
sns.heatmap(corr, mask=np.zeros_like(corr, dtype=np.bool),
cmap=sns.diverging_palette(220, 10, as_cmap=True), annot=True, center=0, square=True, ax=ax)
Удалим созданные сдвиговые переменные:
del train_df['total_cases1']
del train_df['total_cases2']
del train_df['total_cases3']
del train_df['total_cases4']
Посмотрим теперь на распределение признаков:
plt.rc("font", size=13)
plt.figure(figsize=(20,50))
alpha=0.6
ax1 = plt.subplot2grid((6,3), (0,0))
train_df.ndvi_ne.plot(kind="kde",alpha=alpha)
ax1.set_xlabel("ndvi_ne")
ax2 = plt.subplot2grid((6,3),(0,1))
train_df.ndvi_nw.plot(kind="kde",alpha=alpha)
ax2.set_xlabel("ndvi_nw")
ax3 = plt.subplot2grid((6,3),(0,2))
train_df.ndvi_se.plot(kind="kde",alpha=alpha)
ax3.set_xlabel("ndvi_se")
ax4 = plt.subplot2grid((6,3),(1,0))
train_df.ndvi_sw.plot(kind="kde",alpha=alpha)
ax4.set_xlabel("ndvi_sw")
ax5 = plt.subplot2grid((6,3),(1,1))
train_df.reanalysis_air_temp_k.plot(kind="kde",alpha=alpha)
ax5.set_xlabel("reanalysis_air_temp_k")
ax6 = plt.subplot2grid((6,3),(1,2))
train_df.reanalysis_max_air_temp_k.plot(kind="kde",alpha=alpha)
ax6.set_xlabel("reanalysis_max_air_temp_k")
ax7 = plt.subplot2grid((6,3),(2,0))
train_df.reanalysis_min_air_temp_k.plot(kind="kde",alpha=alpha)
ax7.set_xlabel("reanalysis_min_air_temp_k")
ax8 = plt.subplot2grid((6,3),(2,1))
train_df.reanalysis_precip_amt_kg_per_m2.plot(kind="kde",alpha=alpha)
ax8.set_xlabel("reanalysis_precip_amt_kg_per_m2")
ax9 = plt.subplot2grid((6,3),(2,2))
train_df.reanalysis_relative_humidity_percent.plot(kind="kde",alpha=alpha)
ax9.set_xlabel("reanalysis_relative_humidity_percent")
ax10 = plt.subplot2grid((6,3),(3,0))
train_df.precipitation_amt_mm.plot(kind="kde",alpha=alpha)
ax10.set_xlabel("precipitation_amt_mm")
ax11 = plt.subplot2grid((6,3),(3,1))
train_df.reanalysis_specific_humidity_g_per_kg.plot(kind="kde",alpha=alpha)
ax11.set_xlabel("reanalysis_specific_humidity_g_per_kg")
ax12 = plt.subplot2grid((6,3),(3,2))
train_df.reanalysis_tdtr_k.plot(kind="kde",alpha=alpha)
ax12.set_xlabel("reanalysis_tdtr_k")
ax13 = plt.subplot2grid((6,3),(4,0))
train_df.station_avg_temp_c.plot(kind="kde",alpha=alpha)
ax13.set_xlabel("station_avg_temp_c")
ax14 = plt.subplot2grid((6,3),(4,1))
train_df.station_diur_temp_rng_c.plot(kind="kde",alpha=alpha)
ax14.set_xlabel("station_diur_temp_rng_c")
ax15 = plt.subplot2grid((6,3),(4,2))
train_df.station_max_temp_c.plot(kind="kde",alpha=alpha)
ax15.set_xlabel("station_max_temp_c")
ax16 = plt.subplot2grid((6,3),(5,0))
train_df.station_min_temp_c.plot(kind="kde",alpha=alpha)
ax16.set_xlabel("station_min_temp_c")
ax17 = plt.subplot2grid((6,3),(5,1))
train_df.station_precip_mm.plot(kind="kde",alpha=alpha)
ax17.set_xlabel("station_precip_mm")
По своей природе признаки делятся на три части: уровень вегетации, температура, влажность и количество осадков. Для первых трех наблюдается распределение близкое к нормальному, а вот последняя группа имеет биномиальное распределение.
Воспользуемся логарифмированием для нормализации их распределения:
train_df.precipitation_amt_mm = train_df.precipitation_amt_mm.apply(math.log1p)
test_df.precipitation_amt_mm = test_df.precipitation_amt_mm.apply(math.log1p)
train_df.station_precip_mm = train_df.station_precip_mm.apply(math.log1p)
test_df.station_precip_mm = test_df.station_precip_mm.apply(math.log1p)
Визуальный анализ:
del train_df['ndvi_ne']
del train_df['ndvi_nw']
del train_df['ndvi_sw']
del train_df['ndvi_se']
del train_df['station_diur_temp_rng_c']
del test_df['ndvi_ne']
del test_df['ndvi_nw']
del test_df['ndvi_sw']
del test_df['ndvi_se']
del test_df['station_diur_temp_rng_c']
Выбор метрики задан организаторами: средняя абсолютная ошибка(mean absolute error, MAE). Она вычисляется по формуле:
$$\large \begin{array}{rcl}\mathcal{MAE} &=& \frac{1}{n} \sum_{i=1}^n \left |f_i - y_i \right| \end{array}$$Обоснование выбора именно этой метрики организаторами дано не было, однако можно предположить, что они сделали выбор, основываясь на главном преимуществе MAE - простоте интерпретации.
Учитывая заявленное желание предсказывать эпидемии, использование средней квадратичной ошибки, которая гораздо лучше отслеживает выбросы, которыми являются вспышки эпидемий, было бы предпочтительнее.
В рамках данного проекта нам необходимо решить задачу регрессии.
Выбор ограничен сложными начальными условиями: небольшие набор данных и количество признаков.
Как потенциальные варианты решения рассмотрим:
Посмотрим, какие результаты дадут модели на кросс-валидации и выберем лучшую.
cols = ['year', 'weekofyear', 'precipitation_amt_mm', 'reanalysis_air_temp_k', 'reanalysis_max_air_temp_k',
'reanalysis_min_air_temp_k', 'reanalysis_precip_amt_kg_per_m2', 'reanalysis_relative_humidity_percent',
'reanalysis_specific_humidity_g_per_kg', 'reanalysis_tdtr_k', 'station_avg_temp_c', 'station_max_temp_c',
'station_min_temp_c', 'station_precip_mm', 'total_cases', 'total_cases_adj']
features_to_scale = ['year', 'weekofyear', 'precipitation_amt_mm', 'reanalysis_air_temp_k', 'reanalysis_max_air_temp_k',
'reanalysis_min_air_temp_k', 'reanalysis_precip_amt_kg_per_m2',
'reanalysis_relative_humidity_percent', 'reanalysis_specific_humidity_g_per_kg',
'reanalysis_tdtr_k', 'station_avg_temp_c', 'station_max_temp_c',
'station_min_temp_c', 'station_precip_mm']
Разделим данные по городам:
sj_df = pd.concat([train_df[train_df.city=='sj'].copy(), test_df[test_df.city=='sj'].copy()])
iq_df = pd.concat([train_df[train_df.city=='iq'].copy(), test_df[test_df.city=='iq'].copy()])
sj_df = sj_df.drop(['city'], axis=1)
iq_df = iq_df.drop(['city'], axis=1)
id_sj = train_df[train_df.city=='sj'].shape[0]
id_iq = train_df[train_df.city=='iq'].shape[0]
vid_sj = int(0.7*id_sj)
vid_iq = int(0.7*id_iq)
Сдвинем данные для Икитос на 3 недели назад. Образовавшийся пробел ввиду природы данных лучше заполнить последними значениями. Оставим также изначальный вариант, чтобы проверить гипотезу о сдвиге.
iq_base_df = iq_df.copy()
iq_df[cols] = iq_df[cols].shift(-3)
iq_df.iloc[iq_df.shape[0]-3] = iq_df.iloc[iq_df.shape[0]-4].values
iq_df.iloc[iq_df.shape[0]-2] = iq_df.iloc[iq_df.shape[0]-4].values
iq_df.iloc[iq_df.shape[0]-1] = iq_df.iloc[iq_df.shape[0]-4].values
Заполним пропущенные данные методом заполнения вперед:
sj_df.fillna(method='ffill', inplace=True)
iq_df.fillna(method='ffill', inplace=True)
iq_base_df.fillna(method='ffill', inplace=True)
Отмасштабируем признаки:
scaler = StandardScaler()
features_scaled_sj = sj_df[features_to_scale].values
sjx_train = scaler.fit_transform(features_scaled_sj[:id_sj])
sjx_test = scaler.transform(features_scaled_sj[id_sj:])
sj_df[features_to_scale] = np.concatenate([sjx_train, sjx_test])
features_scaled_iq = iq_df[features_to_scale].values
iqx_train = scaler.fit_transform(features_scaled_iq[:id_iq])
iqx_test = scaler.transform(features_scaled_iq[id_iq:])
iq_df[features_to_scale] = np.concatenate([iqx_train, iqx_test])
features_scaled_iq_base = iq_base_df[features_to_scale].values
iq_base_train = scaler.fit_transform(features_scaled_iq_base[:id_iq])
iq_base_test = scaler.transform(features_scaled_iq_base[id_iq:])
iq_base_df[features_to_scale] = np.concatenate([iq_base_train, iq_base_test])
Учитывая, что мы имеем дело с временным рядом, корректным является только выделение валидационной выборки из конца обучающей:
sjx_valid_df = pd.DataFrame(sjx_train[vid_sj:id_sj], columns=features_to_scale)
sjx_valid_df['y'] = sj_df.total_cases[vid_sj:id_sj].values
sjx_valid_df['y_adj'] = sj_df.total_cases_adj[vid_sj:id_sj].values
sjx_train_df = pd.DataFrame(sjx_train[:vid_sj], columns=features_to_scale)
sjx_train_df['y'] = sj_df.total_cases[:vid_sj].values
sjx_train_df['y_adj'] = sj_df.total_cases_adj[:vid_sj].values
sjx_test_df = pd.DataFrame(sjx_test, columns=features_to_scale)
iqx_valid_df = pd.DataFrame(iqx_train[vid_iq:id_iq], columns=features_to_scale)
iqx_valid_df['y'] = iq_df.total_cases[vid_iq:id_iq].values
iqx_valid_df['y_adj'] = iq_df.total_cases_adj[vid_iq:id_iq].values
iqx_train_df = pd.DataFrame(iqx_train[:vid_iq], columns=features_to_scale)
iqx_train_df['y'] = iq_df.total_cases[:vid_iq].values
iqx_train_df['y_adj'] = iq_df.total_cases_adj[:vid_iq].values
iqx_test_df = pd.DataFrame(iqx_test, columns=features_to_scale)
iqx_base_valid_df = pd.DataFrame(iq_base_train[vid_iq:id_iq], columns=features_to_scale)
iqx_base_valid_df['y'] = iq_base_df.total_cases[vid_iq:id_iq].values
iqx_base_valid_df['y_adj'] = iq_base_df.total_cases_adj[vid_iq:id_iq].values
iqx_base_train_df = pd.DataFrame(iq_base_train[:vid_iq], columns=features_to_scale)
iqx_base_train_df['y'] = iq_base_df.total_cases[:vid_iq].values
iqx_base_train_df['y_adj'] = iq_base_df.total_cases_adj[:vid_iq].values
iqx_base_test_df = pd.DataFrame(iq_base_test, columns=features_to_scale)
(+) Кросс-валидация выполнена технически верно, нет утечек данных. Разумно выбрано количество фолдов и разбиение (Random/Stratified или иное), зафиксирован seed. Присутствует объяснение. Объяснены гиперпараметры модели и способ их выбора. Выбор основан на некотором исследовании гипрепараметров модели для данной задачи;
(+/-) Присутствуют незначительные ошибки (например, не зафиксирован seed) или отсутствует объяснение кросс-валидации. Но гиперпараметры модели и способ их выбора должны быть объяснены;
(-/+) Кросс-валидация выполнена со значительными ошибками (например, преобразования данных проводится на всей выборке, таким образом возникает утечка данных из тестовой части выборки и, соответственно, результат кросс-валидации может иметь слишком оптимистичное значение). Гиперпараметры модели и способ их выбора не объяснены;
(-) Отсутствуют.
Выделим обучающую и тестовую выборки, а также целевую переменную для трех наборов данных:
X_sj = sj_df.iloc[:id_sj].drop(['total_cases', 'total_cases_adj'], axis=1)
y_sj = sj_df.iloc[:id_sj].total_cases_adj
y_true_sj = sj_df['total_cases'].iloc[:id_sj]
test_sj = sj_df.iloc[id_sj:].drop(['total_cases', 'total_cases_adj'], axis=1)
X_iq = iq_df.iloc[:id_iq].drop(['total_cases', 'total_cases_adj'], axis=1)
y_iq = iq_df.iloc[:id_iq].total_cases_adj
y_true_iq = iq_df['total_cases'].iloc[:id_iq]
test_iq = iq_df.iloc[id_iq:].drop(['total_cases', 'total_cases_adj'], axis=1)
X_iq_base = iq_base_df.iloc[:id_iq].drop(['total_cases', 'total_cases_adj'], axis=1)
y_iq_base = iq_base_df.iloc[:id_iq].total_cases_adj
y_true_iq_base = iq_base_df['total_cases'].iloc[:id_iq]
test_iq_base = iq_base_df.iloc[id_iq:].drop(['total_cases', 'total_cases_adj'], axis=1)
Будем производить поиск гиперпараметров XGBRegressor с помощью GridSearchCV, используя методику TimeSeriesSplit для осуществления кросс-валидации.
Также будем осуществлять поиск лучшего набора признаков при помощи метода SFS библиотеки mlxtend.
tss = TimeSeriesSplit(n_splits=10)
set_params = {'n_estimators': [100, 200, 400],
'learning_rate': [0.1, 0.3, 0.5, 1],
'max_depth': [4, 5, 6]
}
clf_xgb = xgb.XGBRegressor(objective='reg:linear', seed=13)
grid = GridSearchCV(clf_xgb,
set_params,
cv=tss,
scoring='neg_mean_absolute_error',
n_jobs=-1
)
grid.fit(X_sj, y_sj)
sj_best = grid.best_estimator_
pred_sj = sj_best.predict(test_sj)
SFS1 = SFS(sj_best, k_features=(1, 14), forward=True, floating=True,
verbose=1, scoring='neg_mean_absolute_error', cv=tss, n_jobs=-1)
sfs1 = SFS1.fit(X_sj.as_matrix(), y_sj)
print([X_sg.columns[i-1] for i in sfs1.k_features])
print(sfs1.k_score_)
grid.fit(X_iq, y_iq)
iq_best = grid.best_estimator_
pred_iq = iq_best.predict(test_iq)
print(grid.best_score_)
SFS2 = SFS(iq_best, k_features=(1, 14), forward=True, floating=True,
verbose=1, scoring='neg_mean_absolute_error', cv=tss, n_jobs=-1)
sfs2 = SFS2.fit(X_iq.as_matrix(), y_iq)
print([X_iq.columns[i-1] for i in sfs2.k_features])
print(sfs2.k_score_)
grid.fit(X_iq_base, y_iq_base)
iq_base_best = grid.best_estimator_
pred_iq_base = iq_base_best.predict(test_iq_base)
print(grid.best_score_)
SFS3 = SFS(iq_base_best, k_features=(1, 14), forward=True, floating=True,
verbose=1, scoring='neg_mean_absolute_error', cv=tss, n_jobs=-1)
sfs3 = SFS3.fit(X_iq_base.as_matrix(), y_iq_base)
print([X_iq_base.columns[i-1] for i in sfs3.k_features])
print(sfs3.k_score_)
Мы получили интересный результат: кросс-валидация показала, что для обоих городов лучший набор признаков - это ['precipitation_amt_mm', 'year'].
Гипотеза о том, что несмещенные признаки хуже для Икитос не подтверждается.
Построим модель с новыми признаками и узнаем значение метрики для нее:
grid.fit(X_sj[['precipitation_amt_mm', 'year']], y_sj)
pred_sj = grid.best_estimator_.predict(test_sj[['precipitation_amt_mm', 'year']])
grid.fit(X_iq_base[['precipitation_amt_mm', 'year']], y_iq_base)
pred_iq_base = grid.best_estimator_.predict(test_iq_base[['precipitation_amt_mm', 'year']])
base1 = pd.read_csv('/media/sadworker/DataRed/data/dengai/submission.csv', index_col=[0, 1, 2])
base1.total_cases = np.concatenate([np.round(np.exp(pred_sj) - 1).astype(int),
np.round(np.exp(pred_iq_base) - 1).astype(int)])
base1.to_csv('/media/sadworker/DataRed/data/dengai/subm1imp.csv')
Попробуем использовать в качестве обучающей выборку из 4 наиболее коррелирующих с целевым признаков:
sj_best.fit(X_sj[['year', 'reanalysis_min_air_temp_k',
'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']], y_sj)
pred2_sj = sj_best.predict(test_sj[['year', 'reanalysis_min_air_temp_k',
'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']])
iq_base_best.fit(X_iq_base[['year', 'reanalysis_min_air_temp_k',
'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']], y_iq_base)
pred2_iq_base = iq_base_best.predict(test_iq_base[['year', 'reanalysis_min_air_temp_k',
'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']])
base2 = pd.read_csv('/media/sadworker/DataRed/data/dengai/submission.csv', index_col=[0, 1, 2])
base2.total_cases = np.concatenate([np.round(np.exp(pred2_sj) - 1).astype(int),
np.round(np.exp(pred2_iq_base) - 1).astype(int)])
base2.to_csv('/media/sadworker/DataRed/data/dengai/subm1corr.csv')
Повторим те же операции, но уже с моделью Лассо:
set_params = {'alpha': np.logspace(-10, 4, 50)}
clf_lasso = Lasso(random_state = 13)
gridl = GridSearchCV(clf_lasso,
set_params,
cv=tss,
scoring='neg_mean_absolute_error',
n_jobs=-1
)
gridl.fit(X_sj, y_sj)
predl_sj = gridl.best_estimator_.predict(test_sj)
SFSl1 = SFS(gridl.best_estimator_, k_features=(1, 14), forward=True, floating=True,
verbose=1, scoring='neg_mean_absolute_error', cv=tss, n_jobs=-1)
sfsl1 = SFSl1.fit(X_sj.as_matrix(), y_sj)
print([X_sg.columns[i-1] for i in sfsl1.k_features])
print(sfsl1.k_score_)
gridl.fit(X_iq, y_iq)
predl_iq = gridl.best_estimator_.predict(test_iq)
print(grid.best_score_)
SFSl2 = SFS(gridl.best_estimator_, k_features=(1, 14), forward=True, floating=True,
verbose=1, scoring='neg_mean_absolute_error', cv=tss, n_jobs=-1)
sfsl2 = SFSl2.fit(X_iq.as_matrix(), y_iq)
print([X_iq.columns[i-1] for i in sfsl2.k_features])
print(sfsl2.k_score_)
gridl.fit(X_iq_base, y_iq_base)
predl_iq_base = gridl.best_estimator_.predict(test_iq_base)
print(gridl.best_score_)
SFSl3 = SFS(gridl.best_estimator_, k_features=(1, 14), forward=True, floating=True,
verbose=1, scoring='neg_mean_absolute_error', cv=tss, n_jobs=-1)
sfsl3 = SFSl3.fit(X_iq_base.as_matrix(), y_iq_base)
print([X_iq_base.columns[i-1] for i in sfsl3.k_features])
print(sfsl3.k_score_)
Лассо уступает xgboost на кросс-валидации.
(+) Указаны результаты на тестовой выборке или LB score. Результаты на тестовой выборке сравнимы с результатами на кросс-валидации. Если тестовая выборка создавалась автором проекта, то механизм создания должен быть непредвзят и объяснен (применен разумный механизм выборки, в простейшем случае – рандомизация);
(+/-) Значения метрик на тестовой выборке не сильно отличаются от значений метрик на кросс-валидации и/или тестовая выборка создана предвзято, но есть разумное обоснование этому (пример: заказчик взял тестовую выборку из другого распределения);
(-/+) Значения метрик на тестовой выборке сильно отличаются от значений метрик на кросс-валидации и/или тестовая выборка создана предвзято;
(-) Прогноз для тестовой или отложенной выборки отсутствует.
Модель xgbregressor со всеми признаками показала наилучший результат: MAE=27.6514.
Модель с двумя признаками, отобранными с помощью mlxtend не дала улучшения: MAE=34.0433
Модель с коррелирующими признаками дала: MAE=30.1058
Модель Лассо отстает от xgb на 2-4 единицы по метрике по всем наборам признаков
Посмотрим на результаты по отложенной выборке:
sjx_valid_df = pd.DataFrame(sjx_train[vid_sj:id_sj], columns=features_to_scale)
sjx_valid_df['y'] = sj_df.total_cases[vid_sj:id_sj].values
sjx_valid_df['y_adj'] = sj_df.total_cases_adj[vid_sj:id_sj].values
sjx_train_df = pd.DataFrame(sjx_train[:vid_sj], columns=features_to_scale)
sjx_train_df['y'] = sj_df.total_cases[:vid_sj].values
sjx_train_df['y_adj'] = sj_df.total_cases_adj[:vid_sj].values
sjx_test_df = pd.DataFrame(sjx_test, columns=features_to_scale)
iqx_valid_df = pd.DataFrame(iqx_train[vid_iq:id_iq], columns=features_to_scale)
iqx_valid_df['y'] = iq_df.total_cases[vid_iq:id_iq].values
iqx_valid_df['y_adj'] = iq_df.total_cases_adj[vid_iq:id_iq].values
iqx_train_df = pd.DataFrame(iqx_train[:vid_iq], columns=features_to_scale)
iqx_train_df['y'] = iq_df.total_cases[:vid_iq].values
iqx_train_df['y_adj'] = iq_df.total_cases_adj[:vid_iq].values
iqx_test_df = pd.DataFrame(iqx_test, columns=features_to_scale)
iqx_base_valid_df = pd.DataFrame(iq_base_train[vid_iq:id_iq], columns=features_to_scale)
iqx_base_valid_df['y'] = iq_base_df.total_cases[vid_iq:id_iq].values
iqx_base_valid_df['y_adj'] = iq_base_df.total_cases_adj[vid_iq:id_iq].values
iqx_base_train_df = pd.DataFrame(iq_base_train[:vid_iq], columns=features_to_scale)
iqx_base_train_df['y'] = iq_base_df.total_cases[:vid_iq].values
iqx_base_train_df['y_adj'] = iq_base_df.total_cases_adj[:vid_iq].values
iqx_base_test_df = pd.DataFrame(iq_base_test, columns=features_to_scale)
sj_best.fit(sjx_train_df.drop(['y', 'y_adj'], axis=1), sjx_train_df['y_adj'])
valf_sj = sj_best.predict(sjx_valid_df.drop(['y', 'y_adj'], axis=1))
print('МАЕ на отложенной выборке для всех признаков и Сан-Хуана равно:%0.3f' % mean_absolute_error(sjx_valid_df['y'],
np.exp(valf_sj)-1))
iq_base_best.fit(iqx_base_train_df.drop(['y', 'y_adj'], axis=1), iqx_base_train_df['y_adj'])
valf_iq_base = iq_base_best.predict(iqx_base_valid_df.drop(['y', 'y_adj'], axis=1))
print('МАЕ на отложенной выборке для всех признаков и Икитос равно:%0.3f' % mean_absolute_error(iqx_base_valid_df['y'],
np.exp(valf_iq_base)-1))
sj_best.fit(sjx_train_df[['precipitation_amt_mm', 'year']], sjx_train_df['y_adj'])
val2_sj = sj_best.predict(sjx_valid_df[['precipitation_amt_mm', 'year']])
print('МАЕ на отложенной выборке для 2 признаков и Сан-Хуана равно:%0.3f' % mean_absolute_error(sjx_valid_df['y'],
np.exp(val2_sj)-1))
iq_base_best.fit(iqx_base_train_df[['precipitation_amt_mm', 'year']], iqx_base_train_df['y_adj'])
val2_iq_base = iq_base_best.predict(iqx_base_valid_df[['precipitation_amt_mm', 'year']])
print('МАЕ на отложенной выборке для 2 признаков и Икитос равно:%0.3f' % mean_absolute_error(iqx_base_valid_df['y'],
np.exp(val2_iq_base)-1))
sj_best.fit(sjx_train_df[['year', 'reanalysis_min_air_temp_k',
'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']], sjx_train_df['y_adj'])
val4_sj = sj_best.predict(sjx_valid_df[['year', 'reanalysis_min_air_temp_k',
'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']])
print('МАЕ на отложенной выборке для 4 признаков и Сан-Хуана равно:%0.3f' % mean_absolute_error(sjx_valid_df['y'],
np.exp(val4_sj)-1))
iq_base_best.fit(iqx_base_train_df[['year', 'reanalysis_min_air_temp_k',
'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']], iqx_base_train_df['y_adj'])
val4_iq_base = iq_base_best.predict(iqx_base_valid_df[['year', 'reanalysis_min_air_temp_k',
'reanalysis_precip_amt_kg_per_m2', 'station_min_temp_c']])
print('МАЕ на отложенной выборке для 4 признаков и Икитос равно:%0.3f' % mean_absolute_error(iqx_base_valid_df['y'],
np.exp(val4_iq_base)-1))
Таким образом, результаты на отложенной выборке показывают плохую корреляцию с тестовыми данными.
В условиях малой выборки и шумных данных единственным способом повышать результат видится валидирование по Leaderboard.
В рамках проекта исследованы предоставленные данные, выявлены отдельные тренды и использованы 2 модели с подбором гиперпараметров.
На данный момент результат на Leaderboard очень посредственный - только топ-75%. В процессе кросс-валидации сделан вывод о ее сомнительной ценности для решения данной задачи.
В дальнейшем стоит: