import numpy as np
import seaborn as sns
from sklearn.datasets import load_diabetes

from sklearn.metrics import (
    mean_squared_error,
    mean_absolute_error,
    mean_absolute_percentage_error,
)

sns.set(font_scale=1.3)

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

Рассмотрим данные load_diabetes из библиотеки scikit-learn, которые были созданы для изучения прогрессирования диабета у пациентов на основании медицинских показателей. Идея заключается в использовании десяти базовых физиологических параметров, таких как возраст, индекс массы тела (BMI), среднее артериальное давление и показатели анализа крови, чтобы предсказать, как будет развиваться заболевание через год после начального обследования.

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

1. Загрузка данных

data = load_diabetes()
data.keys()

dict_keys(['data', 'target', 'frame', 'DESCR', 'feature_names', 'data_filename', 'target_filename', 'data_module'])

Функция sklearn.datasets.load_diabetes() возвращает словарь. В поле data записана матрица регрессоров, в которой данные предварительно центрированы и нормированы. В поле target записана мера прогрессирования заболевания в течение года. В поле DESCR можно прочитать подробнее о данных.

Посмотрим на описание датасета.

print(data["DESCR"])

.. _diabetes_dataset:

Diabetes dataset
----------------

Ten baseline variables, age, sex, body mass index, average blood
pressure, and six blood serum measurements were obtained for each of n =
442 diabetes patients, as well as the response of interest, a
quantitative measure of disease progression one year after baseline.

**Data Set Characteristics:**

  :Number of Instances: 442

  :Number of Attributes: First 10 columns are numeric predictive values

  :Target: Column 11 is a quantitative measure of disease progression one year after baseline

  :Attribute Information:
      - age     age in years
      - sex
      - bmi     body mass index
      - bp      average blood pressure
      - s1      tc, total serum cholesterol
      - s2      ldl, low-density lipoproteins
      - s3      hdl, high-density lipoproteins
      - s4      tch, total cholesterol / HDL
      - s5      ltg, possibly log of serum triglycerides level
      - s6      glu, blood sugar level

Note: Each of these 10 feature variables have been mean centered and scaled by the standard deviation times the square root of `n_samples` (i.e. the sum of squares of each column totals 1).

Source URL:
https://www4.stat.ncsu.edu/~boos/var.select/diabetes.html

For more information see:
Bradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani (2004) "Least Angle Regression," Annals of Statistics (with discussion), 407-499.
(https://web.stanford.edu/~hastie/Papers/LARS/LeastAngle_2002.pdf)

Поле data содержит матрицу размера 442 $\times$ 10, где 442 — количество пациентов, а 10 — количество признаков (возраст, пол, и т.д.).

data["data"].shape

(442, 10)

Целевая метка — мера прогрессирования заболевания в течение года.

data["target"].shape

(442,)

Создайте матрицу регрессоров $X$ (data) и столбец наблюдений $Y$ (целевая переменная).

X, Y = <...>

2. Обучение моделей

Разбейте данные случайно на две части — обучающую и тестовую в соотношении 80:20. Для этого используйте функцию train_test_split.

При разбиении датасета стоит зафиксировать случайность для воспроизводимости результатов, зафиксировав random_state=42 в функции разбиения (глобально можно ставить любое число, необязательно именно 42).

Для этого создайте переменную RANDOM_STATE = 42 в начале ноутбука и присвойте ее значение аргументу random_state.

X_train, X_test, y_train, y_test = <...>

print(X_train.shape, X_test.shape, y_train.shape, y_test.shape)

Создайте модель линейной регрессии из sklearn и обучите ее на обучающей части данных.

model = <...>
<...>

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

<...>

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

y_pred = <...>

Реализуйте метрики MSE, MAE, MAPE без использования sklearn и других готовых реализаций.

Пусть $Y_1, ..., Y_n$ — истинные значения, а $\widehat{Y}_1, ..., \widehat{Y}_n$ — предсказания.

Метрика MSE (mean squared error) определяется как $$MSE = \frac{1}{n}\sum_{i=1}^n \left(Y_i - \widehat{Y}_i\right)^2.$$
Метрика MAE (mean absolute error) определяется как $$MAE = \frac{1}{n}\sum_{i=1}^n \left|Y_i - \widehat{Y}_i\right|.$$
Метрика MAPE (mean absolute percentage error) определяется как $$MAPE = \frac{1}{n}\sum_{i=1}^n \left|\frac{Y_i - \widehat{Y}_i}{Y_i}\right|.$$

def MSE(y_true, y_pred):
    <...>

def MAE(y_true, y_pred):
    <...>

def MAPE(y_true, y_pred):
    <...>

Сравните MSE, MAE, MAPE на тренировочной и тестовой выборках. Что вы можете сказать о полученных значениях?

y_pred_train = model.predict(X_train)

print(f"MSE_train = {round(MSE(y_train, y_pred_train), 3)}")
print(f"MSE_test = {round(MSE(y_test, y_pred), 3)}")

print(f"MAE_train = {round(MAE(y_train, y_pred_train), 3)}")
print(f"MAE_test = {round(MAE(y_test, y_pred), 3)}")

print(f"MAPE_train = {round(MAPE(y_train, y_pred_train), 3)}")
print(f"MAPE_test = {round(MAPE(y_test, y_pred), 3)}")

Ответ:

Также ниже вы можете увидеть эти же метрики, реализованные в sklearn.metrics, и сравнить результаты вычислений с ними.

metrics_to_check = [
    (MSE, mean_squared_error, "MSE"), # ваша реализация MSE, реализация из sklearn.metrics, название метрики
    (MAE, mean_absolute_error, "MAE"), # ваша реализация MAE, реализация из sklearn.metrics, название метрики
    (MAPE, mean_absolute_percentage_error, "MAPE"), # ваша реализация MAPE, реализация из sklearn.metrics, название метрики
]

for your_metrics, sklearn_metrics, name in metrics_to_check:
    assert (
        np.abs(your_metrics(y_test, y_pred) - sklearn_metrics(y_test, y_pred))
        < 1e-4
    ), f"Ошибка в реализации {name}"

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

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

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

Проверьте, что в нашем датасете все признаки уже имеют одинаковое среднее и дисперсию.

<...>

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

feature_names = <...>
coefficients = <...>

fig, ax = plt.subplots(figsize=(15, 5))
bars = plt.bar(feature_names, coefficients)

norm = mcolors.Normalize(vmin=0, vmax=max(abs(coefficients)))
cmap = plt.cm.viridis

for bar, value in zip(bars, coefficients):
    bar.set_color(cmap(norm(abs(value))))

sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
cbar = plt.colorbar(sm)
cbar.set_label('Абсолютное значение')

ax.set_title('Коэффициенты модели');

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

print(data["DESCR"].split("baseline\n\n")[-1].split("\n\nNote")[0])

  :Attribute Information:
      - age     age in years
      - sex
      - bmi     body mass index
      - bp      average blood pressure
      - s1      tc, total serum cholesterol
      - s2      ldl, low-density lipoproteins
      - s3      hdl, high-density lipoproteins
      - s4      tch, total cholesterol / HDL
      - s5      ltg, possibly log of serum triglycerides level
      - s6      glu, blood sugar level

Вывод: