План исследования
import os
import numpy as np
import pandas as pd
import scipy.stats as stats
from matplotlib import pyplot as plt
import seaborn as sns
from sklearn.preprocessing import OneHotEncoder, LabelEncoder
from sklearn.linear_model import Ridge
from sklearn.model_selection import train_test_split, KFold, GridSearchCV, validation_curve, learning_curve
from catboost import Pool, CatBoostRegressor
from sklearn.metrics import r2_score
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
Задача данного проекта - предсказать количество посетителей в национальных парках и других достопримечательностях США для последнего доступного года в датасете. Будем рещать задачу регрессии.
Данные взяты с этого сайта. Датасет включает в себя исторические данные о посещаемости национальных парков и достопримечательной США с 1904 по 2016 год.
Загрузим наш основной датасет.
PATH_TO_DATA = ('data')
df = pd.read_csv(os.path.join(PATH_TO_DATA, 'national-parks-visitation.csv'))
df.shape
df.head(2)
Полный датасет состоит из 21560 объектов и 18 признаков.
Загрузим дополнительный датасет, чтобы по GNIS ID вытащить широту и долготу объектов. Данный датасет взят с этого сайта.
gnis_df = pd.read_csv(os.path.join(PATH_TO_DATA, 'gnis-april2018.txt'), sep='|').rename(columns={'FEATURE_ID':'Gnis Id'})
gnis_df.head(3)
Удалим ненужные колонки, преобразуем GNIS ID в текстовый формат и объединим два датасета по GNIS ID.
gnis_df.drop(['FEATURE_NAME', 'FEATURE_CLASS', 'STATE_ALPHA','STATE_NUMERIC',
'COUNTY_NAME', 'COUNTY_NUMERIC', 'PRIMARY_LAT_DMS','PRIM_LONG_DMS',
'SOURCE_LAT_DMS', 'SOURCE_LONG_DMS', 'SOURCE_LAT_DEC','SOURCE_LONG_DEC',
'ELEV_IN_M','ELEV_IN_FT', 'MAP_NAME', 'DATE_CREATED','DATE_EDITED'], axis=1, inplace=True)
gnis_df['Gnis Id'] = gnis_df['Gnis Id'].astype('str')
df = pd.merge(df, gnis_df, how='left', on='Gnis Id')
Посмотрим на тип признаков в главном датасете.
df.info()
Из 18 признаков только 4 вещественных, включая целевой признак Visitors, все остальные категориальные или текстовые. Посмотрим еще раз на данные.
df.head(3)
Первым делом необходимо удалить некоторые признаки из датасета, так как большинство из них являются метаданными датасета и не несут большой информативности для анализа. Кроме того, удалим признак Year в формате даты, так как нам будет достатотчно признака YearRaw.
df.drop(['Created By', 'Measure Selector', 'Year','Date Edit', 'ScrapeURL',
'GIS Notes', 'Geometry', 'Metadata', 'Number of Records'], axis=1, inplace=True)
Проверим данные на пропуски.
df.isnull().sum()
4 признака имеют пропуски. Для каждого из признаков разные решения. Например, для признака Parkname можно заменить пропуски на Unknown. На самом деле, этот признак можно удалить (Сделаем это чуть позднее), т.к. есть другой более точной признак, который содержит в себе полное название достопримечательни Unit Name и не содержит пропусков.
df['Parkname'].fillna('Unknown', inplace=True)
Посмотрим на другие пропуски в датасете.
df[df.isnull().any(axis=1)]
У двух объектов отсутсвует широта и долгота. Проблема в том, что эти два объекта не имеют GNIS ID, по которому мы находили географические данные. В таком случае можно поискать долготу и широту в интернете и заменить их в датасете вручную.
df['PRIM_LAT_DEC'] = df['PRIM_LAT_DEC'].replace(float('nan'), 38.892235)
df['PRIM_LONG_DEC'] = df['PRIM_LONG_DEC'].replace(float('nan'), -77.003689)
Еще есть пропуски в целевой переменной Visitors. YearRaw у данных объектов тоже весьма странный - Total. Посмотрим на YearRaw.
df['YearRaw'].value_counts().head()
df.loc[df['YearRaw'] == 'Total'].head()
386 объектов имеют Total вместо года. Предполагается, что данные объекты были созданы в результате ошибки во время сбора данных. Предлагается удалить эти объекты.
df = df.loc[df['YearRaw'] != 'Total']
Изменим формат целевой переменой Visitors и признака YearRaw.
df['YearRaw'] = df['YearRaw'].astype('int64')
df['Visitors'] = df['Visitors'].astype('int64')
Кроме того, для удобства переименуем YearRaw в Year.
df.rename(columns={'YearRaw': 'Year'}, inplace=True)
Т.к. наша задача предсказать количество посетителей по историческим данных, большое количество исторических данных может затруднить задачу, в связи с чем облегчим себе задачу и оставим данные, начиная с 2000 года.
df = df.loc[df['Year'] >= 2000]
Теперь можно разбить данных не обучающую и тестовую выборки. Посмотрим, какой год последний в датасете.
np.sort(pd.unique(df['Year']))
Будем предсказывать посетителей для 2016 года по данным с 2000 по 2015.
train = df.loc[df['Year'] != 2016]
train.shape
test = df.loc[df['Year'] == 2016]
test.shape
Обучающая выборка включает 5795 объектов, тестовая - 379. Всего 11 признаков.
Сначала будем изучать целевую переменную Visitors. Проведем тесты на нормальность распределения и рассчитаем скошенность и эксцесс распределения.
sns.set(rc={"figure.figsize": (10, 6)})
fig = plt.figure()
sns.distplot(train['Visitors'])
Из графика выше видно, что распределение совсем не похоже на нормальное. Чтобы убедиться в этом, добавим нормальное распределение (из scipy) к графику и сравним его с распределением целевой переменной.
sns.set(rc={"figure.figsize": (10, 6)})
fig = plt.figure()
sns.distplot(train['Visitors'], fit=stats.norm)
Распределение целевой переменной скошенно и сильно отличается от нормального. Проведем стат.тесты.
print("Скошенность: %f" % stats.skew(train['Visitors']))
print("Эксцесс: %f" % stats.kurtosis(train['Visitors']))
stats.normaltest(df['Visitors'].values)
Как видим, скошенность и эксцесс распределения очень большие. В дальнейшем нам будет необходимо преобразовать целевую переменную и попытаться сделать распределение близким к нормальному.
Рассмотрим вещественный признаки.
train.describe()
В датасете всего три вещественных признака, построим графики, чтобы найти какие-нибудь закономерности или зависимости.
visitors_df = train[['Visitors'] + ['Year']]
visitors_df.groupby('Year').sum().plot(figsize=(12,6))
Шкала графика выше вводит заблуждение, на самом деле количество посетителей меняется не так сильно от года в год. Чтобы убедиться в этом, посмотрим на bar-chart.
visitors_df.groupby('Year').sum().plot(kind='bar', rot=0,figsize=(12,6))
Что и требовалось доказать. Тем не менее, число посителей все же выросло в 2015 году. Теперь посмотрим на корреляции между веществеными признаками и целевой переменной.
sns.set(rc={"figure.figsize": (12, 6)})
sns.heatmap(train.corr(), vmin=0, vmax=1,annot=True, cmap="Blues")
plt.show()
Как видно из графика, корреляция между вещественными признаками и целевой переменной достаточно слабая.
Изучим датасет на выбросы. Для этого построим пару графиков и попытаемся найти и интерпретировать выбросы.
sns.set(rc={"figure.figsize": (12, 6)})
sns.violinplot(x='Year', y="Visitors", data=train, palette="muted", split=True)
sns.set(rc={"figure.figsize": (12, 6)})
plt.plot(train['Year'], train['Visitors'], 'r.')
plt.show()
Графики выше показывают, что у датасета есть выбросы, но необходимо изучить природу таких больших значений.
train.sort_values(by='Visitors', ascending=False).head(5)
Самое большое количество посещений на протяжении всего периода 2000-2015 у Parkway, что вполне логично. Таким образом, большие значения объясняются природой данных.
Посмотрим на другие достопримечательности, которые собирают большое количество посетителей.
pd.unique(train['Unit Name'].loc[(train['Visitors'] >= 5000000) & (train['Unit Name'] != 'Blue Ridge Parkway')])
Названия объектов говорят сами за себя, они действительно набирают миллионы посетителей каждый год, поэтому причина больших значений вполне объяснима - природа данных.
В датасете 6 категориальных признаков. Один Parkname из них имеет пропуски, и кроме того, полные названия объектов представлены в признаке Unit Name. Таким образом, предлагается удалить этот признак.
train.drop(['Parkname'], axis=1, inplace=True)
test.drop(['Parkname'], axis=1, inplace=True)
Теперь у нас 5 категориальных признаком. Посмотрим на количество уникальных категорий в каждом из них.
print(len(pd.unique(df['Region'])))
print(len(pd.unique(df['State'])))
print(len(pd.unique(df['Unit Code'])))
print(len(pd.unique(df['Unit Name'])))
print(len(pd.unique(df['Unit Type'])))
Предлагается, каждый такой признак преобразовать с помощью One-Hot кодирования.
В качестве метрики будем использовать r2_score (r_squared) - коэффициент детерминации. Под данной метрикой подразумевают долю дисперсии зависимой переменной, которая объясняется рассматриваемой моделью зависимости, то есть объясняющими переменными (признаками). Выбор обусловлен:
В качестве модели будем использовать:
Некоторые признаки мы уже предобработали ранее. Далее предлагается удалить ненужный признак GNIS ID, т.к. сам по себе он не несет никакой информации, необходимой для модели, но с помощью него мы уже достали широту и долготу.
train.drop(['Gnis Id'], axis=1, inplace=True)
test.drop(['Gnis Id'], axis=1, inplace=True)
Как было отмечено ранее, распределение целевой переменной не прошло тест на нормальность. Поэтому необходимо попытаться сделать его близким к нормальному. Рассмотрим два метотда: логарифмирование и метод Бокса-Кокса
Сначала применим первый метод - логарифмирование.
sns.set(rc={"figure.figsize": (10, 6)})
fig = plt.figure()
sns.distplot(np.log1p((df['Visitors'].values)), fit=stats.norm)
print("Скошенность: %f" % stats.skew(np.log1p((df['Visitors'].values))))
print("Эксцесс: %f" % stats.kurtosis(np.log1p((df['Visitors'].values))))
stats.normaltest(np.log1p((df['Visitors'].values)))
Как видно из результатов теста на нормальность и графика, нам удалось сделать распределение более нормальным, но все еще присутствует скошенность и эксцесс. Применим второй метод - меттод Бокса-Кокса.
# добавим единицу к посетителям, т.к. для данного метода значения должны быть строго больше нуля
stats.boxcox(df['Visitors'].values+1)
sns.set(rc={"figure.figsize": (10, 6)})
fig = plt.figure()
sns.distplot(stats.boxcox(train['Visitors'].values+1)[0], fit=stats.norm)
stats.normaltest(stats.boxcox(train['Visitors'].values+1)[0])
print("Скошенность: %f" % stats.skew(stats.boxcox(train['Visitors'].values+1)[0]))
print("Эксцесс: %f" % stats.kurtosis(stats.boxcox(train['Visitors'].values+1)[0]))
В данном случае, метод Бокса-Кокса был более эффективен, поэтому преобразуем нашу целевую переменную с помощью этого метода. Также отделим целевую переменную от других признаков.
X_train, y_train = train.drop('Visitors', axis=1).values, stats.boxcox(train['Visitors'].values+1)[0]
X_test, y_test = test.drop('Visitors', axis=1).values, stats.boxcox(test['Visitors'].values+1)[0]
Будем использовать готовые выборки (X_train, X_test) для CatBoostRegressor, т.к. он сам преобразует категориальные признаки. Для Ridge преобразуем категориальные признаки с помощью One-Hot Encoding. Напишем функцию для преобразования, а затем применим ее ко всем категоримальным признакам.
def ohe_encoder(column,dftype):
enc = OneHotEncoder(sparse=False)
label_encoder = LabelEncoder()
values = np.array(df['{}'.format(column)].values)
values_labeled = label_encoder.fit_transform(values)
enc.fit(values_labeled.reshape(len(values_labeled), 1))
if dftype == 'train':
df_labeled = label_encoder.transform(train['{}'.format(column)])
else:
df_labeled = label_encoder.transform(test['{}'.format(column)])
df_encoded = enc.transform(df_labeled.reshape(len(df_labeled), 1))
return df_encoded
%%time
train_region_ohe = ohe_encoder(column='Region',dftype = 'train')
train_state_ohe = ohe_encoder(column='State',dftype = 'train')
train_ucode_ohe = ohe_encoder(column='Unit Code',dftype = 'train')
train_uname_ohe = ohe_encoder(column='Unit Name',dftype = 'train')
train_utype_ohe = ohe_encoder(column='Unit Type',dftype = 'train')
%%time
test_region_ohe = ohe_encoder(column='Region',dftype = 'test')
test_state_ohe = ohe_encoder(column='State',dftype = 'test')
test_ucode_ohe = ohe_encoder(column='Unit Code',dftype = 'test')
test_uname_ohe = ohe_encoder(column='Unit Name',dftype = 'test')
test_utype_ohe = ohe_encoder(column='Unit Type',dftype = 'test')
Соединим все полученные признаки в два датасета train_stacked и test_stacked.
train_stacked = np.hstack([train_region_ohe,train_state_ohe,
train_ucode_ohe,train_uname_ohe,
train_utype_ohe,
train[['Year', 'PRIM_LAT_DEC','PRIM_LONG_DEC']].values])
test_stacked = np.hstack([test_region_ohe,test_state_ohe,
test_ucode_ohe,test_uname_ohe,
test_utype_ohe,
test[['Year', 'PRIM_LAT_DEC','PRIM_LONG_DEC']].values])
Для Ridge укажем различные значения для основного параметра alpha - порога регуляризации. Чем выше значение, тем сильнее регуляризация. Далее обучим GridSearchCV для выявления лучших параметров.
params = {'alpha' : [1,0.1,1e2,1e3,1e4]}
%%time
linreg = Ridge()
kf = KFold(n_splits=3, shuffle=True, random_state=17)
linreg_cv = GridSearchCV(linreg, param_grid=params, scoring='r2', cv=kf)
linreg_cv.fit(train_stacked, y_train)
print('Best parameters for Ridge: ', linreg_cv.best_params_)
print('R2 score for Ridge: ', round(linreg_cv.best_score_, 5))
GridSearchCV показал лучшее значение порога регуляризации 0.1, его и будем использовать на тестовой выборке.
Для CatBoostRegressor явно укажем какие из признаков являются категориальными (features). Кроме того, зададим различные значения для глубины дерева (depths) и количества итераций (iterations).
features = [0,1,2,3,4]
depths = [7,8,9,10]
iterations = [250,300]
Напишем собственную кастомную версию кросс-валидации для CatBoost и сравним резльтуты различных комбинаций параметров.
counter = 1
for x in depths:
for y in iterations:
print('Starting '+ str(counter) + ' training: depth - ' + str(x) + ', iterations - ' + str(y))
kf = KFold(n_splits=3,shuffle=True,random_state=17)
results = 0
fold = 1
for train_index, test_index in kf.split(X_train):
cat_train, cat_test = X_train[train_index], X_train[test_index]
train_labels, test_labels = y_train[train_index], y_train[test_index]
cat = CatBoostRegressor(depth=x,iterations=y,random_seed=17,logging_level='Silent')
cat.fit(cat_train, train_labels, cat_features=features)
cat_pred = cat.predict(cat_test)
results += r2_score(test_labels, cat_pred)
print(fold, 'fold score is: ', r2_score(test_labels, cat_pred))
fold += 1
print('Average score for',str(counter),'training is:',results/3)
print('Done...' + '\n')
counter += 1
Лучший резлуьтат (0.8170822248344393) на кросс-валидации показала модель с параметрами depth - 10, iterations - 300. Именно их будем использовать на тестовой выборке.
Построим кривые валидации для Ridge.
plt.figure(figsize=(10,6))
param_range = [1,0.1,1e2,1e3,1e4]
train_scores, test_scores = validation_curve(
Ridge(), train_stacked, y_train, param_name='alpha', param_range=param_range,
cv=3, scoring="r2")
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.title("Validation Curve with Ridge")
plt.xlabel("alpha")
plt.ylabel("Score")
plt.ylim(0.0, 1.1)
lw = 2
plt.semilogx(param_range, train_scores_mean, label="Training score",
color="darkorange", lw=lw)
plt.fill_between(param_range, train_scores_mean - train_scores_std,
train_scores_mean + train_scores_std, alpha=0.2,
color="darkorange", lw=lw)
plt.semilogx(param_range, test_scores_mean, label="Cross-validation score",
color="navy", lw=lw)
plt.fill_between(param_range, test_scores_mean - test_scores_std,
test_scores_mean + test_scores_std, alpha=0.2,
color="navy", lw=lw)
plt.legend(loc="best")
plt.show()
Как мы видим из графика, кривые валидации с увеличением порога регуляризации alpha немного расходятся (справа-налево), это означает, что с увеличением alpha может наблюдаться переобучение. Теперь построим кривые обучения для Ridge.
plt.figure(figsize=(10,6))
train_sizes, train_scores, test_scores = learning_curve(
Ridge(), train_stacked, y_train, cv=3)
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.title("Validation Curve with Ridge")
plt.xlabel("Training examples")
plt.ylabel("Score")
plt.legend(loc="best")
plt.show()
На графике видно, что с увеличением количества примеров. качества на валидации улучшается. Т.к. кривые обучения не сблизились в самом конце, необходимо добавить больше примеров обуччающей выборки.
Теперь, когда у нас есть лучшие параметры, мы готовы обучить модели на тестовой выборке, т.е. предсказать количество посетителей национальных достопримечательностей США.
%%time
ridge = Ridge(alpha=0.1)
ridge.fit(train_stacked, y_train)
ridge_pred = ridge.predict(test_stacked)
r2_score(y_test, ridge_pred)
Коэффициент детерминации достаточно высокий но сильно отличается от результата, показанного на кросс-валидации (0.95302).
%%time
cat = CatBoostRegressor(depth=10, iterations=300,eval_metric='R2',logging_level='Silent')
cat.fit(X_train, y_train, cat_features=features)
cat_pred = cat.predict(X_test)
r2_score(y_test, cat_pred)
CatBoostRegressor показал себя хуже на тестовой выборке (0.7529929633300172), чем лучшая модель на валидации (0.8170822248344393). Вполне вероятно, что виной такой разницы может быть переобучение.
В результате исследования, были построены две модели Ridge из sklearn и CatBoostRegressor из CatBoost. На тестовой выборке CatBoostRegressor (0.7529929633300172) показал результат хуже, чем Ridge (0.8711647097580936). Обе модели переобучаются. Если Ridge переобучается из-за высокого порога alpha, то для выяснения причин переобучения CatBoostRegressor потребуется более детальное исследование. Кривые обучения также показали, что необходимо добавить больше примеров для обучающей выборки.
Предполагаемые улучшения: