Phystech@DataScience

Линейная регрессия

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

Поставить библиотеку можно командой pip install scikit-learn.

import numpy as np
import pandas as pd
pd.options.mode.chained_assignment = None

import matplotlib.pyplot as plt
import seaborn as sns

# обратите внимание, что Scikit-Learn импортируется как sklearn
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn import metrics
from sklearn.metrics import r2_score
from typing import Tuple, Optional, Any, Dict, List
from datetime import datetime


import scipy.stats as sps

sns.set(font_scale=1.3, palette='Set2')
COLOR_1 = (0.5529411764705883, 0.6274509803921569, 0.796078431372549)  # цвета для графиков
COLOR_2 = (0.6509803921568628, 0.8470588235294118, 0.32941176470588235)
RANDOM_STATE = 42
CURR_DATE = str(datetime.now().date())

1. Линейная регрессия на искусственных данных

Данные

Рассмотрим задачу обучения модели на искусственных данных, которые генерируются с помощью нелинейной функции с аддитивным шумом. Данные описываются зависимостью $y = f(x) + \varepsilon$.

• $ f(x) = 3x + 4\sin{x} $ — истинная целевая функция.
• $\varepsilon$ — аддитивный шум, следующий нормальному распределению $\varepsilon \sim \mathcal{N}(0, 1)$.

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

def generate_linear_data(n_samples: int = 100,
                         theta1: float = 3.0,
                         theta2: float = 4.0,
                         X_max: float = 10.0) -> Tuple[np.ndarray, np.ndarray]:
    """
    Генерирует данные для линейной модели y = features * theta1 + sin(features) * theta2 + ε, где ε ~ N(0, 2).

    :param n_samples: Количество выборок.
    :param theta1: Коэффициент для линейной зависимости.
    :param theta2: Коэффициент для зависимости от синуса features.
    :param X_max: Максимальное значение для генерации features.
    :return: Кортеж из двух массивов numpy: features и target_values.
    """
    # Генерация случайных значений для features
    features = np.random.uniform(low=0.0, high=1.0, size=(n_samples, 1)) * X_max  # случайные числа от 0 до X_max
    # Генерация шума
    noise = np.random.normal(loc=0.0, scale=2.0, size=(n_samples, 1))

    # Вычисление target_values согласно функции
    target_values = features * theta1 + np.sin(features) * theta2 + noise

    return features, target_values

n_samples = 60  # количество элементов в выборке
theta1 = 3  # коэффициент для линейной зависимости
theta2 = 4  # коэффициент при sin(x)
X_max = 3   # максимальное значение X

X_grid = np.linspace(0, X_max, n_samples) # равномерная сетка для X от 0 до X_max
X, y = generate_linear_data(n_samples=n_samples, theta1=theta1, theta2=theta2, X_max=X_max)

Визуализируем данные и истинную зависимость $$f(x) = 3x + 4\sin{x}$$

def plot_regression_results(X: np.ndarray,
                            y: np.ndarray,
                            theta1: float,
                            theta2: float,
                            X_grid: np.ndarray,
                            models: Dict[str, Optional[Any]],
                            features_grid_list: List[np.ndarray]) -> None:
    """
    Отображает результаты линейной регрессии вместе с истинной зависимостью и данными.

    :param X: Входные данные, массив numpy.
    :param y: Целевые значения, массив numpy.
    :param theta1: Коэффициент для линейной зависимости.
    :param theta2: Коэффициент для зависимости от синуса.
    :param X_grid: Cетка для отображения исходной зависимости.
    :param models: Словарь, где ключи — подписи, а значения — модели, реализующие метод predict (могут быть None).
    :param features_grid_list: Список массивов сетки признаков для каждой модели.
    """

    plt.figure(figsize=(10, 6))

    # Отображение данных
    sns.scatterplot(x=X.reshape(-1), y=y.reshape(-1), label="Данные",
                    alpha=0.7, s=80, color=COLOR_2)

    # Истинная зависимость
    sns.lineplot(x=X_grid, y=X_grid * theta1 + np.sin(X_grid) * theta2, label="Истинная зависимость",
                 linewidth=3, linestyle='--')

    # Предсказания каждой модели
    for (label, model), features_grid in zip(models.items(), features_grid_list):
        if model is not None:
            sns.lineplot(x=X_grid, y=model.predict(features_grid).reshape(-1),
                         label=f"Предсказание: {label}", linewidth=3)

    plt.title("Сравнение предсказаний моделей")
    plt.xlabel("Признак x")
    plt.ylabel("Таргет")
    plt.ylim(-2, 17)
    plt.legend()
    plt.show()

plot_regression_results(X, y, theta1, theta2, X_grid, {}, [])

Модель

Попробуем предсказать значение $y$ по значению $x$, обучив для начала модель линейной регрессии по признаку $X_1 = x$. Множество модель линейной регрессии в данном случае имеет вид: $$ \widehat{y} = \theta_0 + \theta_1 x, $$ где:

• $ \theta_0 \in \mathbb{R} $ — свободный член (intercept),
• $ \theta_1 \in \mathbb{R} $ — коэффициент при признаке $ x $.

Модель предполагает, что зависимость между $ x $ и $ y $ линейна. В нашей задаче истинная зависимость нелинейная, поэтому модель линейной регрессии может не дать хорошего приближения.

---

Будем использовать LinearRegression.

• fit_intercept — нужно ли включать в модель свободный член. В случае fit_intercept=True модели не нужно передавать признак из всех единиц для того, чтобы она оценивала свободный член. По умолчанию fit_intercept=True.

Важные аргументы конструктора:

Основные методы класса:

• fit(X, y) — обучить линейную регрессию на основе данных X предсказывать целевой признак y. В данном случае под термином "обучить" поднимается вычисление оценки коэффициентов $\widehat{\theta}$.

• predict(X) — предсказать по данным X целевой признак y. В данном случае под термином "предсказать" поднимается вычисление оценки целевого признака $\widehat{y}$.

Используем fit_intercept=True для оценки свободного коэффициента, что позволяет не добавлять в матрицу признаков столбец из единиц.

Замечание. О нормализации данных пойдет речь в домашнем задании.

---

Вспомним, как работают методы fit, predict:

❓ Вопрос ❓

Как получаются оценки параметра $\theta$ и целевого признака в методах fit и predict?

Кликни для показа ответа

fit

Делает оценку параметра $\theta$, минимизируя среднеквадратичную ошибку между предсказанными значениями $ \hat{y} $ и истинными значениями $y$:

$$ \widehat{\theta} = (X^TX)^{-1}X^Ty $$

predict

Получает оценку целевого признака по формуле:

$$ \widehat{y} = X\widehat{\theta}$$

Обучаем модель линейной регрессии на данных, используя только один признак $ x $.

model = LinearRegression(fit_intercept=True) # объявляем модель
model.fit(X, y) # обучаем на признаке x

LinearRegression()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Оценки свободного коэффициента $\theta_0$ и коэффициента $\theta_1$ при $x$

model.intercept_, model.coef_

(array([1.49377888]), array([[3.92862336]]))

Сравним предсказания модели с истинной зависимостью.

models = {
    "Линейная регрессия": model,
}
features_grids = [X_grid.reshape(-1, 1)]

plot_regression_results(X, y, theta1, theta2, X_grid, models, features_grids)

❓ Вопрос ❓

Что можно сделать для увеличения точности предсказаний?

Кликни для показа ответа

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

2. Добавить новые признаки. Примеры преобразований:

• Из временных данных можно извлечь день недели, месяц, год или временные интервалы.

• Полиномиальные признаки. Добавление квадратов, кубов или произведений существующих признаков может помочь модели уловить нелинейные зависимости.

• Создание новых признаков как комбинаций существующих (например, умножение или деление двух признаков).

• Группировка числовых признаков в интервалы (например, возрастные группы).

• Кодирование категориальных признаков. Использование One-Hot Encoding для преобразования категориальных данных в числовые.

Далее рассмотрим второй вариант.

---

Модель с нелинейными признаками

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

Одно нелинейных преобразований - полиномиальное преобразование (в данном случае попробуем использовать признак $x^2$). Добавим $x^2$ в матрицу признаков.

model_2 = LinearRegression(fit_intercept=True) # объявляем модель

data = pd.DataFrame({
    'X': X.flatten(),
    'X^2': X.flatten() ** 2
}) # создаем новый датасет с использованием X^2

data.head()

	X	X^2
0	1.731087	2.996661
1	2.990116	8.940793
2	2.383224	5.679758
3	1.704514	2.905367
4	2.206452	4.868432

❓ Вопрос ❓

Какой вид будет иметь модель линейной регрессии при добавлении нового признака?

Кликни для показа ответа

$$ \widehat{y} = \theta_0 + \theta_1 x + \theta_2 x^2, $$

• $ \theta_0 \in \mathbb{R} $ — свободный член (intercept),
• $ \theta_1$, $\theta_2 \in \mathbb{R} $ — коэффициенты при $ x $ и $x^2$.

model_2.fit(data, y) # обучим модель на новых данных

LinearRegression()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

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

models = {
    "Линейная регрессия": model,
    "Линейная регрессия с x^2": model_2
}
features_grids = [X_grid.reshape(-1, 1),
                  pd.DataFrame({'X': X_grid.flatten(), 'X^2': X_grid.flatten() ** 2})]

plot_regression_results(X, y, theta1, theta2, X_grid, models, features_grids)

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

Результаты

1. Линейная регрессия с одним признаком $ x $ не способна хорошо аппроксимировать нелинейную зависимость.
2. Добавление полиномиального признака $ x^2 $ улучшает качество модели, но всё равно не позволяет точно описать истинную зависимость $ f(x) = 3x + 4\sin{x} $.
3. Для более точного описания данных можно рассмотреть другие нелинейные признаки (например, $ \sin{x} $) или использовать более сложные модели.

2. Постановка задачи с использованием реальных данных

Рассмотрим данные медицинского страхования. В каждой строке представлены признаки для каждого клиента страховой организации.

В данных содержатся следующие признаки:

• birthday — день рожденияж
• sex — пол, возможные значения: female, male;
• bmi — соотношение массы тела квадрату его высоты (body mass index), измеряется в кг/м$^2$, хорошие значения лежат в диапазоне от 18.5 до 24.9;
• children — количество детей;
• smoker — курит ли клиент;
• region — район в США, возможные значения: northeast, southeast, southwest, northwest;
• charges — индивидуальные медицинские расходы, оплачиваемые медицинским страхованием.

Задача: предсказать индивидуальные медицинские расходы по остальным признакам.

Это может быть полезно для таких целей:

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

Оригинальные данные можно посмотреть на Kaggle.

Загрузим данные

data = pd.read_csv('./insurance_miptstats.csv', parse_dates=[0])
data.head()

	birthday	sex	bmi	children	smoker	region	charges
0	2001-12-20	female	27.900	0	yes	southwest	16884.92400
1	2003-03-18	male	33.770	1	no	southeast	1725.55230
2	1992-11-02	male	33.000	3	no	southeast	4449.46200
3	1987-07-27	male	22.705	0	no	northwest	21984.47061
4	1988-11-04	male	28.880	0	no	northwest	3866.85520

Посмотрим на размер таблицы

data.shape

(1338, 7)

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

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

Применим эту функцию к нашим данным, отнеся в тестовую часть 20% данных, а в обучающую — остальные 80%.

train, test = train_test_split(data, test_size=0.2, random_state=RANDOM_STATE)
train.shape, test.shape

((1070, 7), (268, 7))

3. Обучение

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

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

train['age'] = (pd.Timestamp(CURR_DATE) - train['birthday']) / pd.Timedelta(days=365)

Выделим группы признаков

categorial_features = ['sex', 'smoker', 'region']  # категориальные признаки
real_features = ['age', 'bmi', 'children']  # вещественные признаки
target_feature = 'charges'  # целевой признак

Посмотрим на визуализацию совместных распределений вещественных признаков при помощи PairGrid, причем будем разбивать данные по одному признаку из числа категориальных. На графиках приведены:

• данные в виде точек для каждой пары вещественных признаков;
• ядерные оценки плотности для каждой пары вещественных признаков;
• ядерные оценки плотности для всех вещественных признаков по отдельности.

for hue in categorial_features:
    g = sns.PairGrid(train[['bmi', 'age', 'charges', hue]],
                     hue=hue, diag_sharey=False, height=3)

    g.fig.set_size_inches(6,6)
    g.map_lower(sns.kdeplot, alpha=0.6)
    g.map_upper(plt.scatter, alpha=0.3)
    g.map_diag(sns.kdeplot, lw=3, alpha=0.6,
               common_norm=False)  # каждая плотность по отдельности должна давать 1 при интегрировании

    g.add_legend()

По графикам сразу можно сделать следующие выводы:

• расходы растут с увеличением возраста клиента;
• величина расходов больше для курящих людей.

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

---

Напомним, как работают методы fit и predict:

fit

Делает оценку параметра $\theta$, минимизируя среднеквадратичную ошибку между предсказанными значениями $ \hat{y} $ и истинными значениями $y$:

$$ \widehat{\theta} = (X^TX)^{-1}X^Ty $$

predict

Получает оценку целевого признака по формуле:

$$ \widehat{y} = X\widehat{\theta}$$

---

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

Примеры с использованием одного или двух признаков

Признак age

model_age = LinearRegression(fit_intercept=True) # объявляем модель
model_age.fit(train[['age']], train[target_feature]) # обучаем на признаке age

LinearRegression()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

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

# Получаем возраст клиента по дате рождения
test['age'] = (pd.Timestamp(CURR_DATE) - test['birthday']) / pd.Timedelta(days=365)

Получим предсказание модели для тестовых и обучающих данных

y_pred = model_age.predict(test[['age']])

Выведем коэффициент, полученный при обучении модели регресии

print(model_age.coef_[0].round(2))

240.47

❓ Вопрос ❓

Как можно интерпретировать значение коэффициентов в модели линейной регрессии?

Кликни для показа ответа

Коэффициент в линейной регрессии равен изменению таргета при изменении признака на 1. Численно это значит, что в данной модели при увеличении возраста человека на 1 год предсказание затрат увеличивается на ~ 240. С каждым прожитым годом предполагаемые затраты увеличиваются

pred33 = model_age.predict(pd.DataFrame({'age': [33]}))[0]
pred34 = model_age.predict(pd.DataFrame({'age': [34]}))[0]

print('Предсказание расходов.')
print(f'33 года: {pred33:.0f}')
print(f'34 года: {pred34:.0f}')

Предсказание расходов.
33 года: 10764
34 года: 11004

Признак bmi

model_bmi = LinearRegression(fit_intercept=True)
model_bmi.fit(train[['bmi']], train[target_feature])
y_pred = model_bmi.predict(test[['bmi']])

print(f"Bmi coefficient: {model_bmi.coef_[0].round(2)}")

Bmi coefficient: 392.44

По предсказанию модели при увеличении индекса массы тела на 1 расходы увеличиваются на ~392

Признаки bmi и age

model_bmi_age = LinearRegression(fit_intercept=True)
model_bmi_age.fit(train[['bmi', 'age']], train[target_feature])
y_pred = model_bmi_age.predict(test[['bmi', 'age']])

Выведем коэффициенты

print(model_bmi_age.coef_)

[330.40058921 223.55689301]

Изменение индекса массы тела на 1 приводит к большему изменению предсказания медицинских затрат в данной модели, чем изменение возраста на 1 год. Численные значения изменений в зависимости от увеличения bmi/age на 1 будут в точности равны соответствующим коэффициентам в данной модели.

Рассмотрим предказания медицинских затрат для объекта:

• age = 20, bmi = 18 \

И со сдвигом на единицу по одному или обоим признакам:

• age = 21, bmi = 19
• age = 20, bmi = 19
• age = 21, bmi = 18

objects = pd.DataFrame({'bmi': [18, 19, 19, 18],
                        'age': [20, 21, 20, 21]})
preds =  model_bmi_age.predict(objects)
objects['charges_predicted'] = preds

objects

	bmi	age	charges_predicted
0	18	20	3889.012966
1	19	21	4442.970448
2	19	20	4219.413555
3	18	21	4112.569859

Посмотрим на визуализацию

# Грид для предсказаний
bmi_range = np.linspace(17.75, 19.35, 200)
age_range = np.linspace(19.75, 21.25, 200)
bmi_grid, age_grid = np.meshgrid(bmi_range, age_range)
pred_grid = model_bmi_age.predict(pd.DataFrame({'bmi': bmi_grid.ravel(),
                                        'age': age_grid.ravel()})).reshape(bmi_grid.shape)

# Визуализация
plt.figure(figsize=(8, 5))
plt.grid(False)

# Градиент
im = plt.imshow(pred_grid, extent=(bmi_range.min(), bmi_range.max(), age_range.min(), age_range.max()),
           origin='lower', cmap='viridis', alpha=0.8, aspect='auto')

# Линии контуров
contour = plt.contour(bmi_grid, age_grid, pred_grid, colors='grey', linewidths=0.5)

# Отображение признаков
plt.scatter(objects['bmi'], objects['age'], color='red', edgecolor='black', zorder=5, s=80)
for i, row in objects.iterrows():
    plt.text(row['bmi'] + 0.01, row['age'] + 0.03,
             f" {row['charges_predicted']:.0f} ",
             fontsize=15, ha='left', color='white', zorder=6)

# Добавление стрелок с подписями
plt.arrow(18, 20, 1 - 0.06, 0, head_width=0.02, head_length=0.02, fc='lightblue', ec='lightblue', linewidth=2)
plt.text(18.5 - 0.1, 20.05, rf'$\widehat{{\theta}}_1 = {int(model_bmi_age.coef_[0])}$', fontsize=15, color='white')

plt.arrow(18, 20, 0, 1 - 0.06, head_width=0.02, head_length=0.02, fc='lightblue', ec='lightblue', linewidth=2)
plt.text(18.05, 20.5, rf'$\widehat{{\theta}}_2 = {int(model_bmi_age.coef_[1])}$', fontsize=15, color='white')

plt.arrow(18, 21, 1 - 0.06, 0, head_width=0.02, head_length=0.02, fc='lightblue', ec='lightblue', linewidth=2)
plt.text(18.5 - 0.1, 21.05, rf'$\widehat{{\theta}}_1 = {int(model_bmi_age.coef_[0])}$', fontsize=15, color='white')

plt.arrow(19, 20, 0, 1 - 0.06, head_width=0.02, head_length=0.02, fc='lightblue', ec='lightblue', linewidth=2)
plt.text(19.05, 20.5, rf'$\widehat{{\theta}}_2 = {int(model_bmi_age.coef_[1])}$', fontsize=15, color='white')

# Подписи осей
plt.title('Предсказание затрат')
plt.xlabel('BMI')
plt.ylabel('Age')

# Цветовая шкала
cbar = plt.colorbar(im, label='Предсказанные затраты')


# Добавление подписей к линиям контуров
plt.clabel(contour, inline=True, fontsize=8, fmt='%1.0f')

plt.show()

На графике представлены точки, каждая из которых соответствует человеку в возрасте 20 или 21 года с индексом массы тела (BMI) 18 или 19. Стрелки на графике указывают направление изменения положения точки при увеличении возраста или BMI на единицу, и как это влияет на предполагаемые медицинские затраты.

По графику видно, что изменение индекса массы тела на 1 дает больший прирост (+330 к оценке затрат на медицину) к предскаазанию затрат, чем изменение возраста на 1 (+224 к оценке затрат)

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

---

Обработка категориальных признаков

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

• drop — указывает методику, используемую для удаления одной из категорий для каждого объекта. Это может быть полезно для некоторых моделей, например, для линейной регрессии. Возможные значения указаны далее.
  • None — оставляем все признаки.
  • 'first' — удаляет первую категорию для каждого признака. Если признак имеет одну категрию, то он будет удален полностью.
  • 'if_binary' — удаляет первую категорию только для бинарных признаков.
  • массивdrop, где drop[i] — категория в признаке feature X[:, i], которая должна быть удалена.
• sparse_output — возвращает sparse-матрицу, если установлено значение True, иначе — массив.

Важные аргументы конструктора:

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

Основные методы класса:

• fit(X) — обучить кодировщик кодировать признаки на основе данных X. В данном случае под термином "обучить" поднимается определение функций кодирования и декодирования признаков.
• transform(X) — закодировать признаки в данных X.
• fit_transform(X) — обучить кодировщик по данным X и сразу их закодировать.
• inverse_transform(X) — декодировать признаки в данных X, то есть перевести бинарные признаки в исходные категориальные.

---

При построении кодировщика для наших данных учтем ряд особенностей:

• указываем drop='first', то есть одну категорию нужно исключить;
• указываем sparse_output=False, то есть вернуть нужно неразреженную матрицу;
• нужно выполнить обучение, что в данном случае подразумевает построение и сохранение правила преобразования;
• сразу же кодируем признаки из обучающего множества;

encoder = OneHotEncoder(drop='first', sparse_output=False)  # объявляем модель
train_cat = encoder.fit_transform(train[categorial_features])  # обучаем и кодируем
train_cat

array([[0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 0., 0., 1., 0.],
       ...,
       [1., 0., 0., 0., 0.],
       [0., 1., 0., 0., 1.],
       [1., 0., 0., 0., 1.]])

Можем посмотреть на то, как у нас "обучились" категории. Для каждого категориального признака приведен список его категорий

encoder.categories_

[array(['female', 'male'], dtype=object),
 array(['no', 'yes'], dtype=object),
 array(['northeast', 'northwest', 'southeast', 'southwest'], dtype=object)]

Соединим вместе вещественные признаки и закодированные категориальные

X_train = np.hstack([train[real_features], train_cat])
X_train.shape

(1070, 8)

Обучим модель, используя категориальные и вещественные признаки

model_full = LinearRegression(fit_intercept=True)
model_full.fit(X_train, train[target_feature])

LinearRegression()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Посмотрим на результат обучения. Оценки коэффициентов перед признаками

model_full.coef_

array([  257.14354577,   336.56325568,   423.94099187,   -25.48434935,
       23656.64811639,  -370.88646373,  -659.67773002,  -818.2905385 ])

Оценка свободного коэффициента

model_full.intercept_

-13042.059944009723

4. Тестирование и оценка качества

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

# Получаем возраст клиента по дате рождения
test['age'] = (pd.Timestamp(CURR_DATE) - test['birthday']) / pd.Timedelta(days=365)

# Кодируем категориальные признаки с помощью метода transform обученного ранее кодировщика
test_cat = encoder.transform(test[categorial_features])

# Соединяем данные
X_test = np.hstack([test[real_features], test_cat])

Выполним предсказания построенными ранее моделями с помощью метода predict

train_data = [train[['age']], train[['bmi']], train[['bmi', 'age']], X_train]
test_data = [test[['age']], test[['bmi']], test[['bmi', 'age']], X_test]
models = [model_age, model_bmi, model_bmi_age, model_full]
model_names = ['age', 'bmi', 'age+bmi', 'все признаки']

train_preds = {}
test_preds = {}

for model, train_input, test_input, name in zip(models, train_data, test_data, model_names):
    train_pred = model.predict(train_input)
    test_pred = model.predict(test_input)

    train_preds[name] = train_pred
    test_preds[name] = test_pred

Посчитать ошибку предсказания можно разными способами. Наиболее популярный способ — метрика MSE (mean squared error). Пусть $Y_1, ..., Y_n$ — истинные значения, а $\widehat{Y}_1, ..., \widehat{Y}_n$ — предсказания. Тогда метрика MSE определяется как $$MSE = \frac{1}{n}\sum_{i=1}^n \left(Y_i - \widehat{Y}_i\right)^2.$$

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

Посчитаем ее, а точнее — корень из нее, который еще обозначается как RMSE (root MSE)

np.sqrt(((test[target_feature] - test_preds['все признаки']) ** 2).mean())

5793.112670753037

Это значение уже имеет конкретный смысл — насколько в среднем модель отклоняется от истинного значения. То есть в среднем отклонения предсказания величины страховых расходов имеют порядок 6 000 условных единиц. Напомним, что в данных страховые расходы в основной массе принимают значения до 50 000 у.е..

Готовая реализация также есть в sklearn:

metrics.mean_squared_error(test[target_feature], test_preds['все признаки']) ** 0.5

5793.112670753037

---

Другой вариант — метрика MAE (mean absolute error), определяемая как $$MAE = \frac{1}{n}\sum_{i=1}^n \left|Y_i - \widehat{Y}_i\right|.$$

Ее большое преимущество в том, что она не сильно подстраивается по выбросы по сравнению с MSE. Однако ее труднее оптимизировать.

metrics.mean_absolute_error(test[target_feature], test_preds['все признаки'])

4180.1207157432555

---

Еще один популярный вариант — метрика MAPE (mean absolute percentage error), определяемая как $$MAPE = 100\% \cdot \frac{1}{n}\sum_{i=1}^n \frac{\left|Y_i - \widehat{Y}_i\right|}{Y_i}.$$

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

Реализуем ее

def mean_absolute_percentage_error(y_true, y_pred):
    return 100 * (np.abs(y_true - y_pred) / y_true).mean()

Посчитаем для ее для наших данных

mean_absolute_percentage_error(test[target_feature], test_preds['все признаки'])

46.870457081881696

---

Еще для измерения качества модели можно использовать коэффициент детерминации $ R^2 $, который также часто используется для оценки качества предсказательной модели. Этот коэффициент измеряет, какую долю дисперсии зависимой переменной модель сумела объяснить.

Коэффициент детерминации $R^2$ рассчитывается как:

$$ R^2 = 1 - \frac{\sum_{i=1}^n (Y_i - \widehat{Y}_i)^2}{\sum_{i=1}^n (Y_i - \overline{Y})^2}, $$

где

• $\sum_{i=1}^n \left(Y_i - \widehat{Y}_i\right)^2$ — сумма квадратов ошибок (SSE)
• $\overline{Y} = \frac{1}{n} \sum_{i=1}^n Y_i$ — среднее из всех истинных значений
• $\sum_{i=1}^n \left(Y_i - \overline{Y}\right)^2$ — полная сумма квадратов (SST)

❓ Вопрос ❓

Чему равно значение метрики $R^2$ для "идеальной модели"?

Кликни для показа ответа

$R^2$ = 1

Интерпретация $R^2$

• $R^2 = 1$: Такая модель объясняет всю дисперсию данных (идеальный случай).
• $R^2 = 0$: Модель не объясняет дисперсию данных лучше, чем просто среднее значение $Y$.
• $R^2 < 0$: Модель объясняет дисперсию хуже, чем простая модель без регрессии (может случиться, если модель переобучается или плохо подходит для использования на конкретных данных).

Для сравнения посчитаем предсказания и ошибки на обучающем и тестовом множествах.

def get_regression_metrics_df(y_true: pd.Series, y_preds: Dict[str, List[float]]) -> pd.DataFrame:
    """
    Вычисляет метрики качества регрессионной модели для нескольких моделей:
    RMSE, MAE, MAPE, R2
    и выводит их в виде pandas DataFrame.

    Параметры:
    y_true -- истинные значения целевой переменной (pandas Series).
    y_preds -- словарь предсказанных значений целевой переменной для каждой модели (словарь, где ключи - имена моделей, а значения - списки предсказаний).

    Возвращает:
    DataFrame с метриками.
    """
    metric_names = ["RMSE",
                    "MAE",
                    "MAPE",
                    "R2"]

    # Подготовка данных для DataFrame
    data = []
    for model_name, y_pred in y_preds.items():
        rmse = metrics.mean_squared_error(y_true, y_pred) ** 0.5
        mae = metrics.mean_absolute_error(y_true, y_pred)
        mape = metrics.mean_absolute_percentage_error(y_true, y_pred)
        r2 = metrics.r2_score(y_true, y_pred)
        data.append([model_name, rmse, mae, mape, r2])

    # Создание DataFrame
    metrics_df = pd.DataFrame(data, columns=["Model"] + metric_names)

    return metrics_df

Метрики для тестовой выборки

get_regression_metrics_df(test[target_feature], test_preds)

	Model	RMSE	MAE	MAPE	R2
0	age	11659.381084	9172.410208	1.265044	0.124365
1	bmi	12210.039191	9784.652596	1.703504	0.039702
2	age+bmi	11464.256085	9221.883313	1.303051	0.153428
3	все признаки	5793.112671	4180.120716	0.468705	0.783830

Метрики для обучающей выборки

get_regression_metrics_df(train[target_feature], train_preds)

	Model	RMSE	MAE	MAPE	R2
0	age	11528.038976	9042.863931	1.163218	0.079247
1	bmi	11777.698282	9067.951789	1.480413	0.038934
2	age+bmi	11356.453102	9022.142299	1.171528	0.106452
3	все признаки	6105.021537	4208.139727	0.422071	0.741770

---

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

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

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

5. Преобразование признаков

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

Список преобразований, которые могут быть полезны при обучении регрессии:

1. Стандартизация данных. Пусть $x_{ij}$ — значение признака $j$ для объекта $i$. Обозначим $m_j$ — выборочное среднее значение признака $j$ (функция np.mean), а $s^2_j$ — выборочную дисперсию признака $j$ (функция np.var). Тогда стандартизацией является преобразование $$\widetilde{x}_{ij} = \frac{x_{ij} - m_j}{s_j}.$$

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

3. Логарифмирование

4. Обработка выбросов

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

Посмотрим еще раз на истинные и предсказанные медицинские расходы в зависимости от признака smoker.

fig, axes = plt.subplots(1, 2, figsize=(18, 8))

title_fontsize = 18
label_fontsize = 18

sns.scatterplot(ax=axes[0], x=train['bmi'], y=train['charges'],
                hue=train['smoker'], alpha=0.6, s=150)
axes[0].set_title('Зависимость медицинских расходов от BMI (истинные данные)', fontsize=title_fontsize)
axes[0].set_xlabel('BMI', fontsize=label_fontsize)
axes[0].set_ylabel('Charges', fontsize=label_fontsize)

sns.scatterplot(ax=axes[1], x=train['bmi'], y=model_full.predict(X_train),
                hue=train['smoker'], alpha=0.6, s=150)
axes[1].set_title('Зависимость предсказаний медицинских расходов от BMI', fontsize=title_fontsize)
axes[1].set_xlabel('BMI', fontsize=label_fontsize)
axes[1].set_ylabel('Charges', fontsize=label_fontsize)

plt.tight_layout()
plt.show()

Мы видим, что в предсказаниях charges (расходы) увеличиваются на один и тот же коэффициент при росте bmi (индекса массы тела) на 1.

Если посмотреть на график с реальными данными, то можно заметить, что зависимость расходов от bmi отличается в зависимости от признака smoker. Есть предположение, что после добавления признака [smoker * bmi] предсказания для объектов со значением smoker = "yes" улучшатся.

# для обучающей выборки
train['smoker_bmi'] = (train['smoker'] == 'yes') * train['bmi']


# для тестовой выборки
test['smoker_bmi'] = (test['smoker'] == 'yes') * test['bmi']

real_features.extend(['smoker_bmi'])

Обучим модель с новыми признаками

X_train = np.hstack([train[real_features], train_cat])
X_test = np.hstack([test[real_features], test_cat])

model_with_new_feature = LinearRegression(fit_intercept=True)
model_with_new_feature.fit(X_train, train[target_feature])

LinearRegression()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

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

# Создаем подграфики
fig, axes = plt.subplots(1, 3, figsize=(36, 12))

# Первый график
sns.scatterplot(ax=axes[0], x=train['bmi'], y=train['charges'],
                hue=train['smoker'], alpha=0.6, s=170)
axes[0].set_title('Зависимость медицинских расходов от BMI', fontsize=20)
axes[0].set_xlabel('BMI', fontsize=20)
axes[0].set_ylabel('Charges', fontsize=20)

# Второй график
sns.scatterplot(ax=axes[1], x=train['bmi'], y=model_full.predict(np.hstack([train[real_features[:-1]], train_cat])),
                hue=train['smoker'], alpha=0.6, s=170)
axes[1].set_title('Зависимость предсказаний медицинских расходов от BMI', fontsize=20)
axes[1].set_xlabel('BMI', fontsize=20)
axes[1].set_ylabel('Charges', fontsize=20)

# Третий график
sns.scatterplot(ax=axes[2], x=train['bmi'], y=model_with_new_feature.predict(X_train),
                hue=train['smoker'], alpha=0.5, s=170)
axes[2].set_title('Зависимость предсказаний медицинских расходов от BMI с новым признаком', fontsize=20)
axes[2].set_xlabel('BMI', fontsize=20)
axes[2].set_ylabel('Charges', fontsize=20)

plt.tight_layout()
plt.show()

Теперь можно увидеть, что при увеличении bmi на 1, charges меняется на значение, которое зависит от признака smoker, т.е. модель учла эту зависимость.

Получим метрики на тестовой выборке

test_preds_with_new_feature = model_with_new_feature.predict(X_test)
test_preds['модель с новым признаком'] = test_preds_with_new_feature
get_regression_metrics_df(test[target_feature], test_preds)

	Model	RMSE	MAE	MAPE	R2
0	age	11659.381084	9172.410208	1.265044	0.124365
1	bmi	12210.039191	9784.652596	1.703504	0.039702
2	age+bmi	11464.256085	9221.883313	1.303051	0.153428
3	все признаки	5793.112671	4180.120716	0.468705	0.783830
4	модель с новым признаком	4572.470758	2750.922752	0.290165	0.865329

Значение ошибок значительно уменьшилось, с добавлением нового признака получилось повысить качество предсказаний

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

6. Интерпретация результатов

Существует область машинного обучения Explainable AI (XAI), которая занимается тем, что пытается разными способами понять, как работает модель, на какие признаки обращает внимание, насколько является устойчивой... То есть пытается интерпретировать работу модели. Это очень актуальная задача, поскольку таким образом мы можем:

• узнавать новую информацию о наших данных;
• понимать, умеет ли модель обращать внимание на действительно важные признаки;
• находить слабости модели;
• валидировать ее работу и т. д.

Сам подход имеет широкое распространение в сфере бизнеса, поскольку модели, используемые для принятия важных решений, таких как постановка медицинского диагноза, решение о выдаче кредита, беспилотное управление транспортом и т.п., должны соответствовать ГОСТ 59276-2020, регулирующему доверие к ИИ и объяснимость их результатов. Однако не меньший интерес область вызывает и в науке, например, чтобы понимать важность разных генов для формирования фенотипа и формулировать новые гипотезы.

---

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

fit(X)
Метод fit используется для вычисления основных статистик, таких как среднее значение и стандартное отклонение, для каждого признака на предоставленных данных X. Эти статистики затем используются для стандартизации признаков.

transform(X):
После вычисления статистик с помощью fit, метод transform применяется для преобразования данных X. Он вычитает среднее и делит на стандартное отклонение для каждого признака. fit_transform(X): Назначение: Этот метод сочетает в себе fit и transform в одной операции.

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

---

Обучим модель после стандратизации данных

# Инициализация и применение StandardScaler
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Обучение модели линейной регрессии
model_scaled = LinearRegression()
model_scaled.fit(X_train_scaled, train[target_feature])

LinearRegression()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Посмотрим на значения ошибок и коэффициенты в последней линейной модели и сделаем выводы.

# Коэффициенты модели
coefficients = model_scaled.coef_

# Признаки
encoded_features = encoder.get_feature_names_out()
feature_names = real_features + list(encoded_features)

df = pd.DataFrame({
    'Признак': feature_names,
    'Значение коэффициента': coefficients
})

df

	Признак	Значение коэффициента
0	age	3707.567529
1	bmi	118.517407
2	children	561.827116
3	smoker_bmi	18612.936237
4	sex_male	-266.120775
5	smoker_yes	-8572.434745
6	region_northwest	-269.509658
7	region_southeast	-427.542403
8	region_southwest	-535.544448

Визуализируем полученные полученные коэффициенты

fig, ax = plt.subplots(nrows=1, ncols=1, sharey=True, figsize=(8, 6), tight_layout=True)  # Один график


df_sorted = df.sort_values(by='Значение коэффициента', ascending=False)

ax.barh(df_sorted['Признак'], df_sorted['Значение коэффициента'])
ax.set_xlabel('Коэф.')
ax.set_title('Коэффициенты модели')

plt.show()

❓ Вопрос ❓

Какие выводы можно сделать на основе полученных коэффициентов?

Кликни для показа ответа

Коэффициенты регрессии

1. Возраст (age): Коэффициент 3707.57 предполагает, что с увеличением возраста на 1 год, при прочих равных, ожидается увеличение предполагаемых затрат на 3707.57. Это указывает на то, что возраст является существенным фактором.

2. Произведение курильщика и индекса массы тела (smoker_bmi): Высокий коэффициент 18612.94 указывает на сильный положительный эффект. Это может означать, что взаимодействие между курением и bmi значительно увеличивает целевую переменную.

---

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