Множественная регрессия: Расширяем горизонты прогнозирования

В предыдущей статье разобрали простую линейную регрессию ^[1], где целевая переменная зависела от одного фактора, но в реальной жизни всё сложнее. Представьте, что мы прогнозируем стоимость квартиры: она зависит не только от площади, но и от количества комнат, этажа, района, года постройки, наличия парковки и десятков других важных характеристик.

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

1. Учитывать комплекс факторов – строить прогнозы на основе множества признаков одновременно,

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

3. Повышать точность прогнозов по сравнению с одномерными моделями.

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

Математическая основа

Если простая линейная регрессия описывается уравнением прямой, то множественная регрессия стремится описать гиперплоскость в многомерном пространстве.

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

Y=β0+β1X1+β2X2+⋯+βnXn+ϵ

Где:
Y — зависимая переменная (то, что мы предсказываем)
X1, X2, X_n — независимые переменные (те самые предикторы)
β0 — интерсепт (свободный член), это, по сути, значение Y, когда все предикторы равны нулю
βi — коэффициенты регрессии, он нтерпретируются как ожидаемое изменение Y при увеличении Xi на одну единицу, но, при условии, что все остальные предикторы остаются неизменными
ϵ — случайная ошибка ^[2] (шум), которую модель не может объяснить

Как и в случае с простой регрессией, для нахождения оптимальных коэффициентов β используется метод наименьших квадратов, он минимизирует сумму квадратов вертикальных расстояний (остатков ϵ) между фактическими значениями Y и предсказанными моделью значениями. Решение этой задачи в матричном виде выглядит так:

β^{^}=(X^TX)⁻¹X^Tyβ^​

Где:
X — матрица признаков (с добавленным столбцом единиц для интерсепта)
y — вектор целевых значений

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

Расщепление данных

Критически важный этап в машинном обучении ^[3] – это разделение данных на обучающую и тестовую выборки. Модель обучается на одной части данных, а проверяется на другой, которую она “не видела”, это помогает оценить, насколько хорошо модель будет работать с новыми данными, и избежать переобучения (overfitting), когда модель “запоминает” шум в обучающих данных вместо выявления общей закономерности.

Категориальные переменные

Линейная регрессия работает с числами, если у вас есть категориальные признаки (например, город или тип недвижимости), их необходимо преобразовать, самый распространенный способ – это One-Hot Encoding (создание фиктивных переменных), когда каждая категория превращается в отдельную колонку из 0 и 1.

Масштабирование признаков

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

Практическая реализация на Python

Перейдем к практике, рассмотрим построение модели на классическом датасете Advertising, он содержит данные о продажах (Sales) в зависимости от бюджетов на рекламу в трёх каналах: TV, Radio и Newspaper.

Техническое задание: Предсказать объем продаж на основе затрат на рекламу в различных медиа.

Импорт библиотек и загрузка данных

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from sklearn.preprocessing import PolynomialFeatures

# Настройка визуализации
plt.style.use('seaborn-v0_8-whitegrid')
sns.set_palette("husl")

# Загрузка полного датасета Advertising (200 наблюдений)
url = "https://www.tu-chemnitz.de/mathematik/numa/lehre/ds-2018/exercises/Advertising.csv"
df = pd.read_csv(url)

# Проверяем и удаляем лишний столбец индекса, если он есть
if 'Unnamed: 0' in df.columns:
    df = df.drop('Unnamed: 0', axis=1)

# Проверяем названия столбцов и приводим к стандартному виду
df.columns = [col.strip().lower() for col in df.columns]
if 'tv' in df.columns and 'radio' in df.columns and 'newspaper' in df.columns and 'sales' in df.columns:
    df.columns = ['TV', 'Radio', 'Newspaper', 'Sales']

print("Первые 5 строк данных:")
print(df.head())
print(f"nРазмер данных: {df.shape}")
print(f"nНазвания столбцов: {df.columns.tolist()}")

Результат выполнения:

Первые 5 строк данных:
      TV  Radio  Newspaper  Sales
0  230.1   37.8       69.2   22.1
1   44.5   39.3       45.1   10.4
2   17.2   45.9       69.3    9.3
3  151.5   41.3       58.5   18.5
4  180.8   10.8       58.4   12.9

Размер данных: (200, 4)

Разведочный анализ данных (EDA)

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

# Статистика данных
print("Статистика данных:")
print(df.describe())

# Корреляционная матрица
correlation_matrix = df.corr()
print("nКорреляция с целевой переменной Sales:")
print(correlation_matrix['Sales'].sort_values(ascending=False))

# Визуализация корреляционной матрицы
plt.figure(figsize=(8, 6))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', center=0, fmt='.3f')
plt.title('Корреляционная матрица признаков', fontsize=14)
plt.tight_layout()
plt.show()

# Парные графики зависимостей
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

for idx, feature in enumerate(['TV', 'Radio', 'Newspaper']):
    axes[idx].scatter(df[feature], df['Sales'], alpha=0.6, color='steelblue')
    axes[idx].set_xlabel(feature, fontsize=12)
    axes[idx].set_ylabel('Sales', fontsize=12)
    axes[idx].set_title(f'{feature} vs Sales', fontsize=14)

plt.tight_layout()
plt.show()

Статистика данных:
              TV       Radio   Newspaper       Sales
count  200.00000  200.000000  200.000000  200.000000
mean   147.04250   23.264000   30.554000   14.022500
std     85.85426   14.846809   21.778621    5.217457
min      0.70000    0.000000    0.300000    1.600000
25%     74.37500    9.975000   12.750000   10.375000
50%    149.75000   22.900000   25.750000   12.900000
75%    218.82500   36.525000   45.100000   17.400000
max    296.40000   49.600000  114.000000   27.000000

Корреляция с целевой переменной Sales:
Sales        1.000000
TV           0.782224
Radio        0.576223
Newspaper    0.228299
Name: Sales, dtype: float64

Подготовка данных для обучения

# Определяем признаки (X) и целевую переменную (y)
X = df[['TV', 'Radio', 'Newspaper']]
y = df['Sales']

# Разделение данных на обучающую и тестовую выборки (80% / 20%)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

print(f"Размер обучающей выборки: {X_train.shape[0]} samples")
print(f"Размер тестовой выборки: {X_test.shape[0]} samples")

Результат выполнения:

Размер обучающей выборки: 160 samples
Размер тестовой выборки: 40 samples

Обучение модели множественной регрессии

# Создание и обучение модели
model = LinearRegression()
model.fit(X_train, y_train)

# Вывод коэффициентов модели
print("=" * 50)
print("Коэффициенты модели множественной регрессии:")
print("=" * 50)
print(f"Интерсепт (Intercept): {model.intercept_:.4f}")

print("nКоэффициенты при признаках:")
for name, coef in zip(X.columns, model.coef_):
    print(f"  {name:10} : {coef:.4f}")

Результат выполнения:

==================================================
Коэффициенты модели множественной регрессии:
==================================================
Интерсепт (Intercept): 2.9791

Коэффициенты при признаках:
  TV         : 0.0447
  Radio      : 0.1892
  Newspaper  : 0.0028

Коэффициенты модели дают важную информацию о влиянии каждого фактора:

TV (0.0447): При увеличении бюджета на TV-рекламу на 1 тысячу долларов, продажи вырастут в среднем на 44.7 единиц, при неизменных бюджетах Radio и Newspaper.

Radio (0.1892): При увеличении бюджета на радио-рекламу на 1 тысячу долларов, продажи вырастут в среднем на 189.2 единиц это наибольший эффект среди всех каналов.

Newspaper (0.0028): Коэффициент близок к нулю и статистически незначим. Это означает, что при наличии TV и Radio в модели, реклама в газетах практически не влияет на продажи.

Оценка качества модели

# Прогноз на тестовой выборке
y_pred = model.predict(X_test)

# Метрики качества
mse = mean_squared_error(y_test, y_pred)
rmse = np.sqrt(mse)
mae = mean_absolute_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)

print("n" + "=" * 50)
print("Оценка качества модели:")
print("=" * 50)
print(f"Среднеквадратичная ошибка (MSE): {mse:.2f}")
print(f"Корень из MSE (RMSE): {rmse:.2f}")
print(f"Средняя абсолютная ошибка (MAE): {mae:.2f}")
print(f"Коэффициент детерминации (R²): {r2:.4f}")

Результат выполнения:

==================================================
Оценка качества модели:
==================================================
Среднеквадратичная ошибка (MSE): 3.17
Корень из MSE (RMSE): 1.78
Средняя абсолютная ошибка (MAE): 1.46
Коэффициент детерминации (R²): 0.8994

Интерпретация метрик:

Что означают эти цифры на практике:

Если предсказанный объем продаж составляет, например, 15 тысяч единиц, то реальное значение будет находиться в интервале 15 ± 1.78 тысяч единиц (если ориентироваться на RMSE) или 15 ± 1.46 тысяч единиц (если ориентироваться на MAE).

Разница между RMSE (1.78) и MAE (1.46) указывает на наличие небольшого количества ошибок с большим отклонением, которые RMSE “штрафует” сильнее из-за возведения в квадрат.

Визуализация результатов

# Сравнение реальных и предсказанных значений
plt.figure(figsize=(8, 6))
plt.scatter(y_test, y_pred, alpha=0.7)
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--', lw=2)
plt.xlabel('Фактические значения Sales')
plt.ylabel('Предсказанные значения Sales')
plt.title('Сравнение фактических и предсказанных значений')
plt.tight_layout()
plt.show()

# Анализ остатков
residuals = y_test - y_pred

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Гистограмма остатков
axes[0].hist(residuals, bins=10, edgecolor='black', alpha=0.7)
axes[0].axvline(x=0, color='red', linestyle='--')
axes[0].set_xlabel('Остатки')
axes[0].set_ylabel('Частота')
axes[0].set_title('Распределение остатков')

# Остатки vs предсказанные значения
axes[1].scatter(y_pred, residuals, alpha=0.7)
axes[1].axhline(y=0, color='red', linestyle='--')
axes[1].set_xlabel('Предсказанные значения')
axes[1].set_ylabel('Остатки')
axes[1].set_title('Остатки vs предсказанные значения')

plt.tight_layout()
plt.show()

Результат выполнения:

График 1: Сравнение фактических и предсказанных значений

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

График 2: Распределение остатков (гистограмма)

Гистограмма остатков показывает распределение ошибок модели. Остатки приблизительно следуют нормальному распределению с центром в нуле, что подтверждает корректность модели.

График 3: Остатки vs предсказанные значения

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

Сравнение с простой линейной регрессией

Давайте теперь еще сравним, насколько множественная регрессия лучше, чем одномерные модели:

# Простая регрессия только на TV
model_tv = LinearRegression()
model_tv.fit(X_train[['TV']], y_train)
y_pred_tv = model_tv.predict(X_test[['TV']])
r2_tv = r2_score(y_test, y_pred_tv)

# Простая регрессия только на Radio
model_radio = LinearRegression()
model_radio.fit(X_train[['Radio']], y_train)
y_pred_radio = model_radio.predict(X_test[['Radio']])
r2_radio = r2_score(y_test, y_pred_radio)

# Простая регрессия только на Newspaper
model_newspaper = LinearRegression()
model_newspaper.fit(X_train[['Newspaper']], y_train)
y_pred_newspaper = model_newspaper.predict(X_test[['Newspaper']])
r2_newspaper = r2_score(y_test, y_pred_newspaper)

print("=" * 50)
print("Сравнение моделей:")
print("=" * 50)
print(f"R² только TV:           {r2_tv:.4f}")
print(f"R² только Radio:        {r2_radio:.4f}")
print(f"R² только Newspaper:    {r2_newspaper:.4f}")
print(f"R² множественная регрессия: {r2:.4f}")

Результат выполнения:

==================================================
Сравнение моделей:
==================================================
R² только TV:           0.6767
R² только Radio:        0.2634
R² только Newspaper:    0.0299
R² множественная регрессия: 0.8994

Анализ результатов:
Модель только с TV-рекламой объясняет 67.67% дисперсии продаж
Модель только с Radio-рекламой объясняет 26.34% дисперсии
Модель только с Newspaper-рекламой объясняет лишь 2.29% дисперсии
Множественная регрессия, учитывающая TV, Radio и Newspaper одновременно, объясняет 89.94% дисперсии

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

Расширение модели: Полиномиальная регрессия

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

from sklearn.preprocessing import PolynomialFeatures

# Создаем полиномиальные признаки степени 2
poly = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly.fit_transform(X)

# Разделяем данные с полиномиальными признаками
X_train_poly, X_test_poly, y_train, y_test = train_test_split(
    X_poly, y, test_size=0.2, random_state=42
)

# Обучаем модель на полиномиальных признаках
model_poly = LinearRegression()
model_poly.fit(X_train_poly, y_train)

# Оценка качества
y_pred_poly = model_poly.predict(X_test_poly)
r2_poly = r2_score(y_test, y_pred_poly)

print("=" * 50)
print("Полиномиальная регрессия (degree=2):")
print("=" * 50)
print(f"Количество признаков после преобразования: {X_poly.shape[1]}")
print(f"R² полиномиальной модели: {r2_poly:.4f}")
print(f"R² линейной модели: {r2:.4f}")

Результат выполнения:

==================================================
Полиномиальная регрессия (degree=2):
==================================================
Количество признаков после преобразования: 9
R² полиномиальной модели: 0.9869
R² линейной модели: 0.8994

Анализ результатов:
– Полиномиальная модель (со степенью 2) использует 9 признаков вместо исходных 3;
– R² увеличился с 0.8994 до 0.9869  – модель объясняет уже 98.69% дисперсии;
– Однако на малых выборках (всего 10 наблюдений) есть риск переобучения, так как на каждый признак приходится всего 1 наблюдение.

Важно: Добавление слишком большого количества степеней может привести к переобучению. Контролировать сложность модели помогают кросс-валидация и метрика Adjusted R².

Заключение

Множественная линейная регрессия – это мощный и интерпретируемый инструмент анализа данных и ее ключевые преимущества:

1. Интерпретируемость – благодаря ей мы можем сказать: “при увеличении X₁ на единицу, Y изменится на β₁, при неизменных остальных факторах”

2. Учет комплексных зависимостей – возможность моделировать влияние множества факторов одновременно

3. Фундамент для сложных методов – понимание принципов множественной регрессии необходимо для освоения Ridge, LASSO и нейронных сетей

Ключевые моменты, которые нужно запомнить:

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

Датасет Advertising на Kaggle ^[4]

✔️Больше про будни и задачи аналитика данных в моем тг канале 🌸Таня и Данные📊 ^[5]