Машинное обучение с Python: регрессия (полное руководство)

Машинное обучение с Python: регрессия (полное руководство)

Анализ и визуализация данных, Разработка и выбор функций, Разработка и тестирование моделей, Оценка и объяснение

Резюме

В этой статье, используя Data Science и Python, я объясню основные шаги варианта использования регрессии, от анализа данных до понимания выходных данных модели.

Я представлю полезный код Python, который можно легко использовать в других подобных случаях (просто скопировать, вставить, запустить) и пройтись по каждой строке кода с комментариями, чтобы вы могли легко воспроизвести этот пример (ссылка на полный код ниже).

Я буду использовать «Набор данных цен на дома» (ссылка ниже), в котором вам предоставляется несколько объясняющих переменных, описывающих различные аспекты некоторых жилых домов, и задача состоит в том, чтобы предсказать окончательную цену каждого дома.

В частности, я перейду:

• Настройка среды: импорт библиотек и чтение данных
• Анализ данных: понять значение и предсказательную силу переменных
• Разработка функций: извлечение функций из необработанных данных
• Предварительная обработка: разделение данных, обработка отсутствующих значений, кодирование категориальных переменных, масштабирование
• Выбор функций: сохраняйте только самые релевантные переменные
• Дизайн модели: базовый уровень, поезд, проверка, тест
• Оценка производительности: ознакомьтесь с метриками
• Объяснимость: понять, как модель делает прогнозы

Настраивать

Прежде всего, мне нужно импортировать следующие библиотеки.

## for data
import pandas as pd
import numpy as np
## for plotting
import matplotlib.pyplot as plt
import seaborn as sns
## for statistical tests
import scipy
import statsmodels.formula.api as smf
import statsmodels.api as sm
## for machine learning
from sklearn import model_selection, preprocessing, feature_selection, ensemble, linear_model, metrics, decomposition
## for explainer
from lime import lime_tabular

Затем я считываю данные в фрейм данных pandas. Исходный набор данных содержит 81 столбец, но для целей этого руководства я буду работать с подмножеством из 12 столбцов.

dtf = pd.read_csv("data_houses.csv")
cols = ["OverallQual","GrLivArea","GarageCars", 
        "GarageArea","TotalBsmtSF","FullBath",
        "YearBuilt","YearRemodAdd",
        "LotFrontage","MSSubClass"]
dtf = dtf[["Id"]+cols+["SalePrice"]]
dtf.head()

Подробную информацию о столбцах можно найти в предоставленной ссылке на набор данных.

Обратите внимание, что каждая строка таблицы представляет собой конкретный дом (или наблюдение). Если вы работаете с другим набором данных, у которого нет такой структуры, в которой каждая строка представляет собой наблюдение, вам необходимо суммировать данные и преобразовать их.

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

Приступим, ладно?

Анализ данных

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

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

'''
Recognize whether a column is numerical or categorical.
:parameter
    :param dtf: dataframe - input data
    :param col: str - name of the column to analyze
    :param max_cat: num - max number of unique values to recognize a column as categorical
:return
    "cat" if the column is categorical or "num" otherwise
'''
def utils_recognize_type(dtf, col, max_cat=20):
    if (dtf[col].dtype == "O") | (dtf[col].nunique() < max_cat):
        return "cat"
    else:
        return "num"

Эта функция очень полезна и может использоваться несколько раз. В качестве иллюстрации я построю тепловую карту фрейма данных и визуализирую тип столбцов и отсутствующие данные.

dic_cols = {col:utils_recognize_type(dtf, col, max_cat=20) for col in dtf.columns}
heatmap = dtf.isnull()
for k,v in dic_cols.items():
 if v == "num":
   heatmap[k] = heatmap[k].apply(lambda x: 0.5 if x is False else 1)
 else:
   heatmap[k] = heatmap[k].apply(lambda x: 0 if x is False else 1)
sns.heatmap(heatmap, cbar=False).set_title('Dataset Overview')
plt.show()
print("\033[1;37;40m Categerocial ", "\033[1;30;41m Numeric ", "\033[1;30;47m NaN ")

Всего 1460 строк и 12 столбцов:

• каждая строка таблицы представляет собой конкретный дом (или наблюдение), идентифицированный Id, поэтому я установлю его как индекс (или первичный ключ таблицы для любителей SQL).
• SalePrice - это зависимая переменная, которую мы хотим понять и спрогнозировать, поэтому я переименую столбец «Y».
• TotalQuall, GarageCars, FullBath и MSSubClass - категориальные переменные, а остальные - числовые.
• Только LotFrontage содержит недостающие данные.

dtf = dtf.set_index("Id")
dtf = dtf.rename(columns={"SalePrice":"Y"})

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

Во-первых, давайте посмотрим на одномерные распределения (распределение вероятностей только одной переменной). Гистограмма идеально подходит, чтобы дать приблизительное представление о плотности основного распределения отдельных числовых данных. Я рекомендую использовать коробчатую диаграмму для графического изображения групп данных через их квартили. Например, изобразим целевую переменную:

x = "Y"
fig, ax = plt.subplots(nrows=1, ncols=2,  sharex=False, sharey=False)
fig.suptitle(x, fontsize=20)
### distribution
ax[0].title.set_text('distribution')
variable = dtf[x].fillna(dtf[x].mean())
breaks = np.quantile(variable, q=np.linspace(0, 1, 11))
variable = variable[ (variable > breaks[0]) & (variable < 
                    breaks[10]) ]
sns.distplot(variable, hist=True, kde=True, kde_kws={"shade": True}, ax=ax[0])
des = dtf[x].describe()
ax[0].axvline(des["25%"], ls='--')
ax[0].axvline(des["mean"], ls='--')
ax[0].axvline(des["75%"], ls='--')
ax[0].grid(True)
des = round(des, 2).apply(lambda x: str(x))
box = '\n'.join(("min: "+des["min"], "25%: "+des["25%"], "mean: "+des["mean"], "75%: "+des["75%"], "max: "+des["max"]))
ax[0].text(0.95, 0.95, box, transform=ax[0].transAxes, fontsize=10, va='top', ha="right", bbox=dict(boxstyle='round', facecolor='white', alpha=1))
### boxplot 
ax[1].title.set_text('outliers (log scale)')
tmp_dtf = pd.DataFrame(dtf[x])
tmp_dtf[x] = np.log(tmp_dtf[x])
tmp_dtf.boxplot(column=x, ax=ax[1])
plt.show()

Средняя цена дома для этой группы населения составляет 181 тыс. Долларов, распределение сильно неравномерно и с обеих сторон наблюдаются выбросы.

Кроме того, столбчатая диаграмма подходит для понимания частоты ярлыков для одной категориальной переменной. Возьмем, к примеру, переменную FullBath (количество ванных комнат): она имеет порядковый характер (2 ванные комнаты ›1 ванная), но не непрерывна (в доме не может быть 1,5 ванных комнат), поэтому ее можно проанализировать. как категоричный.

x = "Y"
ax = dtf[x].value_counts().sort_values().plot(kind="barh")
totals= []
for i in ax.patches:
    totals.append(i.get_width())
total = sum(totals)
for i in ax.patches:
     ax.text(i.get_width()+.3, i.get_y()+.20, 
     str(round((i.get_width()/total)*100, 2))+'%', 
     fontsize=10, color='black')
ax.grid(axis="x")
plt.suptitle(x, fontsize=20)
plt.show()

В большинстве домов есть 1 или 2 ванные комнаты, есть некоторые выбросы с 0 и 3 ванными комнатами.

Я перейду на следующий уровень анализа и изучу двумерное распределение, чтобы понять, обладает ли FullBath прогностической силой для прогнозирования Y. Это будет случай категориального (FullBath) против числового (Y), поэтому я буду действовать следующим образом:

• разделите совокупность (весь набор наблюдений) на 4 выборки: часть домов с 0 ванной (FullBath = 0), 1 ванной (FullBath = 1 ) и так далее…
• Постройте и сравните плотности 4 выборок, если распределения отличаются, тогда переменная является прогнозной, потому что 4 группы имеют разные шаблоны.
• Сгруппируйте числовую переменную (Y) по ячейкам (подвыборкам) и постройте состав каждой ячейки; если пропорция категорий одинакова во всех из них, то переменная не является прогнозирующей.
• Постройте и сравните ящичные диаграммы для 4 выборок, чтобы определить различное поведение выбросов.

cat, num = "FullBath", "Y"
fig, ax = plt.subplots(nrows=1, ncols=3,  sharex=False, sharey=False)
fig.suptitle(x+"   vs   "+y, fontsize=20)
            
### distribution
ax[0].title.set_text('density')
for i in dtf[cat].unique():
    sns.distplot(dtf[dtf[cat]==i][num], hist=False, label=i, ax=ax[0])
ax[0].grid(True)
### stacked
ax[1].title.set_text('bins')
breaks = np.quantile(dtf[num], q=np.linspace(0,1,11))
tmp = dtf.groupby([cat, pd.cut(dtf[num], breaks, duplicates='drop')]).size().unstack().T
tmp = tmp[dtf[cat].unique()]
tmp["tot"] = tmp.sum(axis=1)
for col in tmp.drop("tot", axis=1).columns:
     tmp[col] = tmp[col] / tmp["tot"]
tmp.drop("tot", axis=1).plot(kind='bar', stacked=True, ax=ax[1], legend=False, grid=True)
### boxplot   
ax[2].title.set_text('outliers')
sns.catplot(x=cat, y=num, data=dtf, kind="box", ax=ax[2])
ax[2].grid(True)
plt.show()

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

Если вас не убеждает интуиция глаза, вы всегда можете прибегнуть к старой доброй статистике и провести тест. В этом случае категориального (FullBath) против числового (Y) я бы использовал односторонний ANOVA-тест ». По сути, он проверяет, существенно ли различаются средние значения двух или более независимых выборок, поэтому, если значение p достаточно мало (

cat, num = "FullBath", "Y"
model = smf.ols(num+' ~ '+cat, data=dtf).fit()
table = sm.stats.anova_lm(model)
p = table["PR(>F)"][0]
coeff, p = None, round(p, 3)
conclusion = "Correlated" if p < 0.05 else "Non-Correlated"
print("Anova F: the variables are", conclusion, "(p-value: "+str(p)+")")

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

Чтобы проверить справедливость этого первого вывода, мне нужно будет проанализировать поведение целевой переменной по отношению к GrLivArea (жилая площадь над землей в квадратных футах). Это случай числового (GrLivArea) против числового (Y), поэтому я построю 2 графика:

• сначала я сгруппирую значения GrLivArea в интервалы и сравню среднее значение (и медиану) Y в каждом интервале. Если кривая не плоская, то переменная прогнозирующий, потому что у бункеров разные шаблоны.
• Во-вторых, я буду использовать диаграмму рассеяния с распределениями двух переменных по сторонам.

x, y = "GrLivArea", "Y"
### bin plot
dtf_noNan = dtf[dtf[x].notnull()]
breaks = np.quantile(dtf_noNan[x], q=np.linspace(0, 1, 11))
groups = dtf_noNan.groupby([pd.cut(dtf_noNan[x], bins=breaks, 
           duplicates='drop')])[y].agg(['mean','median','size'])
fig, ax = plt.subplots(figsize=figsize)
fig.suptitle(x+"   vs   "+y, fontsize=20)
groups[["mean", "median"]].plot(kind="line", ax=ax)
groups["size"].plot(kind="bar", ax=ax, rot=45, secondary_y=True,
                    color="grey", alpha=0.3, grid=True)
ax.set(ylabel=y)
ax.right_ax.set_ylabel("Observazions in each bin")
plt.show()
### scatter plot
sns.jointplot(x=x, y=y, data=dtf, dropna=True, kind='reg', 
              height=int((figsize[0]+figsize[1])/2) )
plt.show()

GrLivArea является прогнозным, здесь прослеживается четкая закономерность: в среднем, чем больше дом, тем выше цена, хотя есть некоторые выбросы с размером выше среднего и относительно низкой ценой.

Как и раньше, мы можем проверить корреляцию между этими двумя переменными. Поскольку обе они числовые, я бы не оценил коэффициент корреляции Пирсона: предполагая, что две переменные независимы (нулевая гипотеза), он проверяет, имеют ли две выборки линейную связь. Если p-значение достаточно мало (

x, y = "GrLivArea", "Y"
dtf_noNan = dtf[dtf[x].notnull()]
coeff, p = scipy.stats.pearsonr(dtf_noNan[x], dtf_noNan[y])
coeff, p = round(coeff, 3), round(p, 3)
conclusion = "Significant" if p < 0.05 else "Non-Significant"
print("Pearson Correlation:", coeff, conclusion, "(p-value: "+str(p)+")")

FullBath и GrLivArea - это примеры функций прогнозирования, поэтому я оставлю их для моделирования.

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

Функциональная инженерия

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

sns.catplot(x="MSSubClass", y="Y", data=dtf, kind="box")

Есть много категорий, и трудно понять, каково распределение внутри каждой из них. Поэтому я сгруппирую эти категории в кластеры: классы с более высоким значением Y (например, MSSubClass 60 и 120) войдут в кластер «max», классы с более низкими ценами. (например, MSSubClass 30, 45, 180) перейдет в «минимальный» кластер, остальные будут сгруппированы в «средний» кластер.

## define clusters
MSSubClass_clusters = {"min":[30,45,180], "max":[60,120], "mean":[]}
## create new columns
dic_flat = {v:k for k,lst in MSSubClass_clusters.items() for v in lst}
for k,v in MSSubClass_clusters.items():
    if len(v)==0:
        residual_class = k 
dtf[x+"_cluster"] = dtf[x].apply(lambda x: dic_flat[x] if x in 
                          dic_flat.keys() else residual_class)
## print
dtf[["MSSubClass","MSSubClass_cluster","Y"]].head()

Таким образом я уменьшил количество категорий с 15 до 3, что намного лучше для анализа:

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

Предварительная обработка

Предварительная обработка данных - это этап подготовки необработанных данных, чтобы сделать их пригодными для модели машинного обучения. Особенно:

1. каждое наблюдение должно быть представлено одной строкой, другими словами, у вас не может быть двух строк, описывающих одного и того же пассажира, потому что они будут обрабатываться моделью отдельно (набор данных уже в такой форме, поэтому ✅). Более того, каждый столбец должен быть функцией, поэтому вы не должны использовать Id в качестве предиктора, поэтому такая таблица называется «матрица характеристик».
2. Набор данных должен быть разделен как минимум на два набора: модель должна быть обучена на значительной части вашего набора данных (так называемый «набор поездов») и протестирована на меньшем наборе («тестовый набор» ).
3. Отсутствующие значения следует заменить на что-нибудь, иначе ваша модель может взбеситься.
4. Категориальные данные должны быть закодированы, что означает преобразование меток в целые числа, поскольку машинное обучение ожидает числа, а не строки.
5. Рекомендуется масштабировать данные, это помогает нормализовать данные в определенном диапазоне и ускорять вычисления в алгоритме.

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

## split data
dtf_train, dtf_test = model_selection.train_test_split(dtf, 
                      test_size=0.3)
## print info
print("X_train shape:", dtf_train.drop("Y",axis=1).shape, "| X_test shape:", dtf_test.drop("Y",axis=1).shape)
print("y_train mean:", round(np.mean(dtf_train["Y"]),2), "| y_test mean:", round(np.mean(dtf_test["Y"]),2))
print(dtf_train.shape[1], "features:", dtf_train.drop("Y",axis=1).columns.to_list())

Следующий шаг: столбец LotFrontage содержит некоторые отсутствующие данные (17%), которые необходимо обработать. С точки зрения машинного обучения правильно сначала разделить на обучение и тестирование, а затем заменить NA на среднее значение обучающей выборки.

dtf_train["LotFrontage"] = dtf_train["LotFrontage"].fillna(dtf_train["LotFrontage"].mean())

Созданный мной новый столбец MSSubClass_cluster содержит категориальные данные, которые следует закодировать. Я буду использовать метод One-Hot-Encoding, преобразовывая 1 категориальный столбец с n уникальными значениями в n-1 фиктивные.

## create dummy
dummy = pd.get_dummies(dtf_train["MSSubClass_cluster"], 
                       prefix="MSSubClass_cluster",drop_first=True)
dtf_train= pd.concat([dtf_train, dummy], axis=1)
print( dtf_train.filter(like="MSSubClass_cluster",axis=1).head() )
## drop the original categorical column
dtf_train = dtf_train.drop("MSSubClass_cluster", axis=1)

И последнее, но не менее важное: я собираюсь масштабировать функции. Для задач регрессии часто желательно преобразовать как входные, так и целевые переменные. Я буду использовать RobustScaler, который преобразует функцию путем вычитания медианы и последующего деления на межквартильный диапазон (значение 75% - значение 25%). Преимущество этого средства масштабирования в том, что на него меньше влияют выбросы.

## scale X
scalerX = preprocessing.RobustScaler(quantile_range=(25.0, 75.0))
X = scaler.fit_transform(dtf_train.drop("Y", axis=1))
dtf_scaled= pd.DataFrame(X, columns=dtf_train.drop("Y", 
                        axis=1).columns, index=dtf_train.index)
## scale Y
scalerY = preprocessing.RobustScaler(quantile_range=(25.0, 75.0))
dtf_scaled[y] = scalerY.fit_transform(
                    dtf_train[y].values.reshape(-1,1))
dtf_scaled.head()

Выбор функции

Выбор функций - это процесс выбора подмножества соответствующих переменных для построения модели машинного обучения. Это упрощает интерпретацию модели и снижает переобучение (когда модель слишком сильно адаптируется к обучающим данным и плохо работает за пределами набора поездов).

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

Я объясню на примере: GarageCars сильно коррелирует с GarageArea, потому что они оба предоставляют одинаковую информацию (насколько велик гараж, один с точки зрения того, сколько автомобилей может подходит, другой - в квадратных футах). Давайте вычислим корреляционную матрицу, чтобы увидеть это:

corr_matrix = dtf_train.corr(method="pearson")
sns.heatmap(corr_matrix, vmin=-1., vmax=1., annot=True, fmt='.2f', cmap="YlGnBu", cbar=True, linewidths=0.5)
plt.title("pearson correlation")

Один из GarageCars и GarageArea может быть ненужным, и мы можем отказаться от него и оставить наиболее полезный (т. Е. С наименьшим p -значение или то, которое больше всего снижает энтропию).

Линейная регрессия - это линейный подход к моделированию взаимосвязи между скалярным откликом и одной или несколькими независимыми переменными. Одномерные линейные регрессионные тесты широко используются для тестирования индивидуального эффекта каждого из множества регрессоров: сначала вычисляется корреляция между каждым регрессором и целью, затем выполняется F-тест ANOVA.

РИДЖЕВАЯ регуляризация особенно полезна для смягчения проблемы мультиколлинеарности в линейной регрессии, которая обычно возникает в моделях с большим количеством параметров.

X = dtf_train.drop("Y", axis=1).values
y = dtf_train["Y"].values
feature_names = dtf_train.drop("Y", axis=1).columns
## p-value
selector = feature_selection.SelectKBest(score_func=  
               feature_selection.f_regression, k=10).fit(X,y)
pvalue_selected_features = feature_names[selector.get_support()]

## regularization
selector = feature_selection.SelectFromModel(estimator= 
              linear_model.Ridge(alpha=1.0, fit_intercept=True), 
                                 max_features=10).fit(X,y)
regularization_selected_features = feature_names[selector.get_support()]
 
## plot
dtf_features = pd.DataFrame({"features":feature_names})
dtf_features["p_value"] = dtf_features["features"].apply(lambda x: "p_value" if x in pvalue_selected_features else "")
dtf_features["num1"] = dtf_features["features"].apply(lambda x: 1 if x in pvalue_selected_features else 0)
dtf_features["regularization"] = dtf_features["features"].apply(lambda x: "regularization" if x in regularization_selected_features else "")
dtf_features["num2"] = dtf_features["features"].apply(lambda x: 1 if x in regularization_selected_features else 0)
dtf_features["method"] = dtf_features[["p_value","regularization"]].apply(lambda x: (x[0]+" "+x[1]).strip(), axis=1)
dtf_features["selection"] = dtf_features["num1"] + dtf_features["num2"]
dtf_features["method"] = dtf_features["method"].apply(lambda x: "both" if len(x.split()) == 2 else x)
sns.barplot(y="features", x="selection", hue="method", data=dtf_features.sort_values("selection", ascending=False), dodge=False)

Синие объекты выбираются как ANOVA, так и RIDGE, остальные выбираются только первым статистическим методом.

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

X = dtf_train.drop("Y", axis=1).values
y = dtf_train["Y"].values
feature_names = dtf_train.drop("Y", axis=1).columns.tolist()
## call model
model = ensemble.GradientBoostingRegressor()
## Importance
model.fit(X,y)
importances = model.feature_importances_
## Put in a pandas dtf
dtf_importances = pd.DataFrame({"IMPORTANCE":importances, 
            "VARIABLE":feature_names}).sort_values("IMPORTANCE", 
            ascending=False)
dtf_importances['cumsum'] =  
            dtf_importances['IMPORTANCE'].cumsum(axis=0)
dtf_importances = dtf_importances.set_index("VARIABLE")
    
## Plot
fig, ax = plt.subplots(nrows=1, ncols=2, sharex=False, sharey=False)
fig.suptitle("Features Importance", fontsize=20)
ax[0].title.set_text('variables')
    dtf_importances[["IMPORTANCE"]].sort_values(by="IMPORTANCE").plot(
                kind="barh", legend=False, ax=ax[0]).grid(axis="x")
ax[0].set(ylabel="")
ax[1].title.set_text('cumulative')
dtf_importances[["cumsum"]].plot(kind="line", linewidth=4, 
                                 legend=False, ax=ax[1])
ax[1].set(xlabel="", xticks=np.arange(len(dtf_importances)), 
          xticklabels=dtf_importances.index)
plt.xticks(rotation=70)
plt.grid(axis='both')
plt.show()

Действительно интересно, что во всех представленных методах преобладают TotalQual, GrLivArea и TotalBsmtSf.

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

X_names = ['OverallQual', 'GrLivArea', 'TotalBsmtSF', "GarageCars"]
X_train = dtf_train[X_names].values
y_train = dtf_train["Y"].values
X_test = dtf_test[X_names].values
y_test = dtf_test["Y"].values

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

Модельный дизайн

Наконец, пришло время построить модель машинного обучения. Сначала я запущу простую линейную регрессию и буду использовать ее в качестве основы для более сложной модели, такой как алгоритм повышения градиента.

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

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

Я визуализирую результаты проверки, построив график прогнозируемых значений против фактического Y. В идеале все точки должны располагаться близко к диагональной линии, где прогнозируемое = фактическое.

## call model
model = linear_model.LinearRegression()
## K fold validation
scores = []
cv = model_selection.KFold(n_splits=5, shuffle=True)
fig = plt.figure()
i = 1
for train, test in cv.split(X_train, y_train):
    prediction = model.fit(X_train[train],
                 y_train[train]).predict(X_train[test])
    true = y_train[test]
    score = metrics.r2_score(true, prediction)
    scores.append(score)
    plt.scatter(prediction, true, lw=2, alpha=0.3, 
                label='Fold %d (R2 = %0.2f)' % (i,score))
    i = i+1
plt.plot([min(y_train),max(y_train)], [min(y_train),max(y_train)], 
         linestyle='--', lw=2, color='black')
plt.xlabel('Predicted')
plt.ylabel('True')
plt.title('K-Fold Validation')
plt.legend()
plt.show()

Среднее значение R в квадрате линейной регрессии составляет 0,77. Давайте посмотрим, как проходит проверка повышения градиента:

Модель повышения градиента обеспечивает лучшую производительность (средний R в квадрате 0,83), поэтому я буду использовать ее для прогнозирования тестовых данных:

## train
model.fit(X_train, y_train)
## test
predicted = model.predict(X_test)

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

predicted = scalerY.inverse_transform( 
                  predicted.reshape(-1,1) ).reshape(-1)

Оценка

Момент истины, мы собираемся увидеть, стоит ли вся эта тяжелая работа того. Все дело в том, чтобы изучить, насколько вариативность Y может объяснить модель и как распределяются ошибки.

Я буду оценивать модель, используя следующие общие показатели: R в квадрате, средняя абсолютная ошибка (MAE) и среднеквадратичная ошибка (RMSD). Последние два являются мерой ошибки между парными наблюдениями, выражающими одно и то же явление. Поскольку ошибки могут быть как положительными (фактическое ›прогноз), так и отрицательными (фактическое‹ прогнозирование), вы можете измерить абсолютное значение и квадрат значения каждой ошибки.

## Kpi
print("R2 (explained variance):", round(metrics.r2_score(y_test, predicted), 2))
print("Mean Absolute Perc Error (Σ(|y-pred|/y)/n):", round(np.mean(np.abs((y_test-predicted)/predicted)), 2))
print("Mean Absolute Error (Σ|y-pred|/n):", "{:,.0f}".format(metrics.mean_absolute_error(y_test, predicted)))
print("Root Mean Squared Error (sqrt(Σ(y-pred)^2/n)):", "{:,.0f}".format(np.sqrt(metrics.mean_squared_error(y_test, predicted))))
## residuals
residuals = y_test - predicted
max_error = max(residuals) if abs(max(residuals)) > abs(min(residuals)) else min(residuals)
max_idx = list(residuals).index(max(residuals)) if abs(max(residuals)) > abs(min(residuals)) else list(residuals).index(min(residuals))
max_true, max_pred = y_test[max_idx], predicted[max_idx]
print("Max Error:", "{:,.0f}".format(max_error))

Модель объясняет 86% дисперсии целевой переменной. В среднем прогнозы имеют ошибку в 20 тысяч долларов, или они ошибаются на 11%. Самая большая ошибка на тестовом наборе составила более 170 тысяч долларов. Мы можем визуализировать ошибки, сопоставляя прогнозируемые значения с фактическими данными и остатком (ошибкой) каждого прогноза.

## Plot predicted vs true
fig, ax = plt.subplots(nrows=1, ncols=2)
from statsmodels.graphics.api import abline_plot
ax[0].scatter(predicted, y_test, color="black")
abline_plot(intercept=0, slope=1, color="red", ax=ax[0])
ax[0].vlines(x=max_pred, ymin=max_true, ymax=max_true-max_error, color='red', linestyle='--', alpha=0.7, label="max error")
ax[0].grid(True)
ax[0].set(xlabel="Predicted", ylabel="True", title="Predicted vs True")
ax[0].legend()
    
## Plot predicted vs residuals
ax[1].scatter(predicted, residuals, color="red")
ax[1].vlines(x=max_pred, ymin=0, ymax=max_error, color='black', linestyle='--', alpha=0.7, label="max error")
ax[1].grid(True)
ax[1].set(xlabel="Predicted", ylabel="Residuals", title="Predicted vs Residuals")
ax[1].hlines(y=0, xmin=np.min(predicted), xmax=np.max(predicted))
ax[1].legend()
plt.show()

Вот она, самая большая ошибка -170k: модель предсказала около 320k, в то время как истинное значение этого наблюдения составляет около 150k. Похоже, что большинство ошибок лежат между 50k и -50k, давайте лучше посмотрим на распределение остатков и посмотрим, выглядит ли оно примерно нормально:

fig, ax = plt.subplots()
sns.distplot(residuals, color="red", hist=True, kde=True, kde_kws={"shade":True}, ax=ax)
ax.grid(True)
ax.set(yticks=[], yticklabels=[], title="Residuals distribution")
plt.show()

Объяснимость

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

Пакет Lime может помочь нам создать объяснитель. Для иллюстрации я возьму случайное наблюдение из набора тестов и посмотрю, что предсказывает модель:

print("True:", "{:,.0f}".format(y_test[1]), "--> Pred:", "{:,.0f}".format(predicted[1]))

Модель предсказывала цену этого дома в размере 194 870 долларов. Почему? Давайте воспользуемся объяснителем:

explainer = lime_tabular.LimeTabularExplainer(training_data=X_train, feature_names=X_names, class_names="Y", mode="regression")
explained = explainer.explain_instance(X_test[1], model.predict, num_features=10)
explained.as_pyplot_figure()

Основными факторами для этого конкретного прогноза являются то, что у дома есть большой подвал (TotalBsmft ›1,3k), он был построен из высококачественных материалов (CompleteQual› 6) и был построен недавно (Год постройки ›2001).

График Прогнозируемое против фактических - отличный инструмент, чтобы показать, как прошло тестирование, но я также построил плоскость регрессии, чтобы наглядно показать выбросы, которые модель не спрогнозировала правильно. Поскольку он лучше работает для линейных моделей, я буду использовать линейную регрессию для подбора двумерных данных. Чтобы отобразить данные в двух измерениях, требуется некоторое уменьшение размерности (процесс уменьшения количества функций путем получения набора основных переменных). Я приведу пример использования алгоритма PCA для суммирования данных в 2 переменные, полученные с помощью линейных комбинаций признаков.

## PCA
pca = decomposition.PCA(n_components=2)
X_train_2d = pca.fit_transform(X_train)
X_test_2d = pca.transform(X_test)
## train 2d model
model_2d = linear_model.LinearRegression()
model_2d.fit(X_train, y_train)
## plot regression plane
from mpl_toolkits.mplot3d import Axes3D
ax = Axes3D(plt.figure())
ax.scatter(X_test[:,0], X_test[:,1], y_test, color="black")
X1 = np.array([[X_test.min(), X_test.min()], [X_test.max(), 
               X_test.max()]])
X2 = np.array([[X_test.min(), X_test.max()], [X_test.min(), 
               X_test.max()]])
Y = model_2d.predict(np.array([[X_test.min(), X_test.min(), 
                     X_test.max(), X_test.max()], 
                    [X_test.min(), X_test.max(), X_test.min(), 
                     X_test.max()]]).T).reshape((2,2))
Y = scalerY.inverse_transform(Y)
ax.plot_surface(X1, X2, Y, alpha=0.5)
ax.set(zlabel="Y", title="Regression plane", xticklabels=[], 
       yticklabels=[])
plt.show()

Заключение

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

В исследовательском разделе я проанализировал случай одной категориальной переменной, одной числовой переменной и то, как они взаимодействуют друг с другом. Я привел пример разработки функций, извлекающей функцию из необработанных данных. Что касается предварительной обработки, я объяснил, как обрабатывать отсутствующие значения и категориальные данные. Я показал разные способы выбора правильных функций, как использовать их для построения регрессионной модели и как оценить производительность. В последнем разделе я дал несколько предложений о том, как улучшить объяснимость вашей модели машинного обучения.

Важное примечание: я не рассказал, что происходит после утверждения вашей модели для развертывания. Просто имейте в виду, что вам нужно построить конвейер для автоматической обработки новых данных, которые вы будете периодически получать.

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

Надеюсь, вам понравилось! Не стесняйтесь обращаться ко мне с вопросами и отзывами или просто для того, чтобы поделиться своими интересными проектами.

LinkedIn | Instagram | Твиттер | GitHub

Эта статья является частью серии Машинное обучение с помощью Python, см. Также: