Множественная регрессия: расчёт коэффициентов методом МНК

Парная регрессия описывает отклик одним фактором, но реальные процессы почти всегда зависят от нескольких причин сразу: продажи определяются и рекламой, и ценой, и сезоном. Множественная регрессия связывает зависимую переменную $Y$ с набором факторов $X_1, X_2, \dots, X_k$ линейным уравнением, а её коэффициенты показывают, как меняется отклик при росте каждого фактора на единицу при фиксированных остальных. Главная вычислительная задача - найти эти коэффициенты так, чтобы модель как можно ближе ложилась на данные. Разберём, как строится нормальная система уравнений, откуда берётся матричная формула и как проверить расчёт на конкретном примере. Ниже стоит калькулятор, в котором можно подвигать наклоны и сразу увидеть, при каких значениях достигается минимум ошибки.

Что такое множественная регрессия

Линейная множественная регрессия - модель вида

$Y = b_0 + b_1 X_1 + b_2 X_2 + \dots + b_k X_k + \varepsilon,$

где $Y$ - зависимая переменная (отклик), $X_1, \dots, X_k$ - объясняющие факторы, $b_0$ - свободный член, $b_1, \dots, b_k$ - коэффициенты регрессии, а $\varepsilon$ - случайная ошибка. Коэффициент $b_j$ - это частная производная отклика по фактору $X_j$: на сколько единиц изменится $Y$, если $X_j$ вырастет на единицу, а все прочие факторы останутся неизменными («при прочих равных», ceteris paribus).

В отличие от парной регрессии, где можно просто провести прямую через облако точек, здесь модель геометрически задаёт гиперплоскость в $(k+1)$-мерном пространстве. Для двух факторов это наклонная плоскость над плоскостью $(X_1, X_2)$, и расчёт коэффициентов сводится к тому, чтобы поставить эту плоскость как можно ближе к наблюдаемым точкам.

Метод наименьших квадратов: что минимизируем

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

$S(b_0, b_1, \dots, b_k) = \sum_{i=1}^{n} \left( y_i - b_0 - b_1 x_{i1} - \dots - b_k x_{ik} \right)^2 \to \min.$

Это и есть метод наименьших квадратов (МНК). Квадрат, а не модуль, выбран не случайно: функция $S$ становится гладкой и выпуклой, у неё единственный минимум, который находится приравниванием частных производных к нулю. Идея минимизации именно квадратов остатков общая для всей регрессии - она же лежит в основе нормальной системы парного МНК, только там неизвестных всего два.

Рассчитать для своей задачи

Вертикальные отрезки от точек до плоскости - это остатки $e_i = y_i - \hat y_i$. МНК наклоняет и сдвигает плоскость так, чтобы суммарная длина этих отрезков (в квадрате) была наименьшей. В калькуляторе выше видно то же самое в двумерной проекции: чем дальше ваши наклоны от оптимальных, тем сильнее точки «вашей модели» уходят от биссектрисы, и тем больше растёт сумма квадратов остатков.

Нормальная система уравнений

Чтобы найти минимум $S$, берут частные производные по каждому коэффициенту и приравнивают их к нулю. Для модели с двумя факторами получается система трёх нормальных уравнений:

$$\begin{aligned} b_0 n + b_1 \sum x_1 + b_2 \sum x_2 &= \sum y, \\ b_0 \sum x_1 + b_1 \sum x_1^2 + b_2 \sum x_1 x_2 &= \sum x_1 y, \\ b_0 \sum x_2 + b_1 \sum x_1 x_2 + b_2 \sum x_2^2 &= \sum x_2 y. \end{aligned}$$

Здесь все суммы берутся по $n$ наблюдениям. Это обычная линейная система: три уравнения, три неизвестных $b_0, b_1, b_2$. Решают её любым методом - Крамера, Гаусса или обращением матрицы. Для $k$ факторов система разрастается до $(k+1)$ уравнений, но структура та же: в каждой строке стоят суммы попарных произведений факторов.

Матричная форма расчёта

Когда факторов много, систему компактнее записать через матрицы. Соберём матрицу наблюдений (дизайн-матрицу) $X$ размера $n \times (k+1)$, где первый столбец - единицы (под свободный член), и вектор откликов $y$:

$X = \begin{pmatrix} 1 & x_{11} & x_{12} \\ 1 & x_{21} & x_{22} \\ \vdots & \vdots & \vdots \\ 1 & x_{n1} & x_{n2} \end{pmatrix}, \quad y = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix}.$

Тогда вектор коэффициентов $b = (b_0, b_1, b_2)^\mathsf{T}$ находится по знаменитой формуле:

$b = (X^\mathsf{T} X)^{-1} X^\mathsf{T} y.$

Рассчитать для своей задачи

Произведение $X^\mathsf{T} X$ - это как раз матрица коэффициентов нормальной системы (с теми же суммами квадратов и произведений), а $X^\mathsf{T} y$ - её правая часть. Формула $b = (X^\mathsf{T} X)^{-1} X^\mathsf{T} y$ - это просто решение нормальной системы $X^\mathsf{T} X b = X^\mathsf{T} y$, записанное через обратную матрицу. Именно её применяют все пакеты (Excel, Python statsmodels, R) под капотом, хотя численно вместо явного обращения чаще используют разложения.

Пример расчёта по шагам

Возьмём учебный набор из калькулятора: 10 наблюдений, где $X_1$ - расходы на рекламу, $X_2$ - относительная цена, $Y$ - продажи. Порядок расчёта коэффициентов множественной регрессии такой:

1. Считаем суммы. Находим $n=10$, $\sum x_1$, $\sum x_2$, $\sum y$, суммы квадратов $\sum x_1^2$, $\sum x_2^2$ и попарные произведения $\sum x_1 x_2$, $\sum x_1 y$, $\sum x_2 y$.
2. Подставляем в нормальную систему. Получаем три уравнения с тремя неизвестными.
3. Решаем систему методом Крамера: считаем главный определитель и три вспомогательных, делим один на другой. Так и работает калькулятор внутри - он держит МНК-решение этой системы и сравнивает с тем, что вы подобрали ползунками.
4. Записываем уравнение $\hat Y = b_0 + b_1 X_1 + b_2 X_2$ и проверяем знаки: положительный $b_1$ означает, что реклама повышает продажи, отрицательный $b_2$ - что рост цены их снижает.

Подвигайте в калькуляторе наклоны $b_1$ и $b_2$: сумма квадратов остатков минимальна ровно в одной точке - это и есть единственное МНК-решение. Любое отклонение от него увеличивает ошибку, что наглядно подтверждает выпуклость функции $S$.

Коэффициент детерминации R²

После расчёта коэффициентов модель оценивают по качеству подгонки. Главная мера - коэффициент детерминации:

$R^2 = 1 - \frac{\sum (y_i - \hat y_i)^2}{\sum (y_i - \bar y)^2} = 1 - \frac{\text{RSS}}{\text{TSS}},$

где RSS - остаточная сумма квадратов, TSS - полная. $R^2$ показывает долю разброса отклика, объяснённую факторами: $R^2 = 0{,}95$ значит, что модель объясняет 95 % вариации $Y$. Калькулятор пересчитывает $R^2$ на лету - в МНК-минимуме он максимален.

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

Интерпретация и проверка значимости

Рассчитанные коэффициенты бессмысленно толковать без проверки значимости. Стандартный набор тестов:

• t-критерий для каждого коэффициента: проверяет гипотезу $b_j = 0$. Если $|t_j|$ больше критического, фактор значим. $t_j = b_j / S_{b_j}$, где $S_{b_j}$ - стандартная ошибка коэффициента.
• F-критерий для уравнения в целом: значимо ли уравнение, то есть объясняют ли факторы вместе хоть что-то. $F = \frac{R^2 / k}{(1 - R^2)/(n - k - 1)}$.
• Доверительные интервалы: $b_j \pm t_{\text{кр}} \cdot S_{b_j}$ показывают диапазон, в котором лежит истинный коэффициент.

Экономический смысл наклонов: коэффициент $b_1 = 3{,}2$ означает, что при росте рекламы на единицу продажи растут в среднем на 3,2 единицы - но только при неизменной цене. Эта оговорка «при прочих равных» - главное отличие частного эффекта в множественной регрессии от простой парной зависимости.

Частые ошибки

• Забывают столбец единиц в матрице $X$. Без него формула $b = (X^\mathsf{T} X)^{-1} X^\mathsf{T} y$ даёт модель без свободного члена - прямую через начало координат, что почти всегда неверно.
• Путают определитель в методе Крамера. В числитель идёт определитель с заменённым столбцом, в знаменатель - главный. Перестановка даёт обратный по смыслу результат.
• Толкуют коэффициент без оговорки «при прочих равных». В множественной регрессии $b_j$ - это частный эффект; сравнивать его с парным коэффициентом по тому же фактору некорректно, они отвечают на разные вопросы.
• Гонятся за высоким $R^2$. Добавление любого фактора поднимает $R^2$; ориентироваться надо на скорректированный $R^2$ и значимость коэффициентов, иначе модель переобучится.
• Игнорируют мультиколлинеарность. Если факторы сильно коррелируют между собой, матрица $X^\mathsf{T} X$ почти вырождена, обратная матрица «взрывается», и коэффициенты получаются неустойчивыми с огромными стандартными ошибками.

FAQ

Чем расчёт коэффициентов множественной регрессии отличается от парной? Принцип тот же - МНК и минимизация суммы квадратов остатков. Но вместо двух нормальных уравнений их становится $(k+1)$, и удобнее перейти к матричной форме. Главное смысловое отличие: каждый коэффициент теперь показывает частный эффект фактора при фиксированных остальных, а не общую зависимость.

Можно ли посчитать коэффициенты вручную или нужен компьютер? Для двух факторов и небольшой выборки реально решить нормальную систему вручную методом Крамера или Гаусса - это типичная учебная задача. Для трёх и более факторов или больших данных без матричных пакетов (Excel ЛИНЕЙН, Python, R) уже не обойтись: ручное обращение матрицы $4\times4$ слишком трудоёмко.

Что делать, если матрица $X^\mathsf{T} X$ необратима? Это значит, что факторы линейно зависимы (полная мультиколлинеарность) - например, один фактор является суммой двух других. Решение: убрать дублирующий фактор, либо применить регуляризацию (гребневая регрессия добавляет к диагонали $X^\mathsf{T} X$ малое число и делает её обратимой).

Коротко

Расчёт коэффициентов множественной регрессии - это минимизация суммы квадратов остатков методом наименьших квадратов. Приравнивание частных производных к нулю даёт нормальную систему из $(k+1)$ уравнений, коэффициенты которой - суммы квадратов и произведений факторов. Компактно решение записывается матричной формулой $b = (X^\mathsf{T} X)^{-1} X^\mathsf{T} y$. После расчёта модель оценивают по $R^2$ и скорректированному $R^2$, а коэффициенты проверяют t- и F-критериями. Каждый коэффициент $b_j$ толкуется как изменение отклика при росте фактора $X_j$ на единицу при прочих равных.