Множественная регрессия: матричная форма, R² и VIF с примером

Две квартиры в одном районе: 56 м² в восьми минутах ходьбы от метро и 64 м² в двадцати. Первая стоит 8,2 млн, вторая — 8,6. Что двигает цену сильнее — лишние восемь метров или лишние двенадцать минут? Парная регрессия из урока о корреляции и регрессии такой вопрос не возьмёт: она умеет один фактор за раз, и её коэффициент при площади впитает в себя всё, что с площадью связано, — включая удалённость от центра. Нужна модель, которая держит в руках несколько причин сразу, — множественная регрессия.

Что значит «при прочих равных»?#

Модель линейна по каждому фактору одновременно: $\hat{y}=b_0+b_1{x}_1+b_2{x}_2$, где $x_1$ — площадь в м², $x_2$ — минуты пешком до метро, $\hat{y}$ — предсказанная цена в млн руб. Коэффициент $b_1$ читается так: на сколько изменится предсказанная цена при увеличении площади на один квадратный метр, если удалённость при этом не менялась. Коэффициент $b_2$ — наоборот, про минуты при зафиксированной площади. Фраза «при прочих равных» — не украшение, а вся суть метода: каждый коэффициент отвечает за свой фактор при замороженных остальных.

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

Возьмём пять квартир — минимальная выборка, на которой уже видно всё:

Как матрицы решают систему разом?#

Метод наименьших квадратов из парной регрессии остаётся тем же: подбираем коэффициенты, чтобы сумма квадратов остатков $\sum (y_i-{\hat{y}}_i{)}^2$ была минимальной. Минимум ищут частными производными по каждому коэффициенту, приравненными к нулю. В матричной записи весь вывод сжимается в одну строку: столбец единиц и столбцы факторов складываются в матрицу $X$ размера $5\times 3$, цены — в столбец $y$, коэффициенты — в вектор $\beta$:

Расшифровка: $X^{T}X$ накапливает суммы, квадраты и произведения факторов, $X^{T}y$ — произведения факторов с откликом, а $(X^{T}X{)}^{-1}$ переворачивает систему. Эта формула — компактная запись системы нормальных уравнений: по одному уравнению на каждый коэффициент. Для двух факторов после центрирования остаётся система 2×2 — решим её руками, по формулам Крамера, знакомым по определителям; с тремя и больше факторами включается метод Гаусса.

В сыром виде, без центрирования, та же система — три уравнения с тремя неизвестными: $5b_0+250b_1+60b_2=36,4$; $250b_0+13500b_1+2700b_2=1938$; $60b_0+2700b_1+850b_2=399$. Когда дорешаете до ответа $(3;0,1;-0,06)$, подставьте его сюда — все три равенства сойдутся копейка в копейку. Центрирование ничего не меняет по сути: оно лишь аккуратнее перегруппирует те же числа и выносит свободный член из борьбы.

Считаем суммы по таблице: $\sum x_1=250$, $\sum x_2=60$, $\sum y=36,4$; средние ${\bar{x}}_1=50$, ${\bar{x}}_2=12$, $\bar{y}=7,28$. Дальше работают центрированные суммы квадратов и произведений — при таком ходе система становится чище, а свободный член возвращается в самом конце:

1. 1
   Система для наклонов: $S_{11}{b}_1+S_{12}{b}_2=S_{1y}$ и $S_{12}{b}_1+S_{22}{b}_2=S_{2y}$, то есть $1000b_1-300b_2=118$ и $-300b_1+130b_2=-37,8$.
2. 2
   Определитель системы: $\Delta =1000\cdot 130-(-300{)}^2=130000-90000=40000$ — не ноль, решение существует.
3. 3
   Первый столбец меняем на правые части: ${\Delta }_1=118\cdot 130-(-37,8)\cdot (-300)=15340-11340=4000$, откуда $b_1=4000/40000=0,1$.
4. 4
   Второй столбец меняем на правые части: ${\Delta }_2=1000\cdot (-37,8)-(-300)\cdot 118=-37800+35400=-2400$, откуда $b_2=-2400/40000=-0,06$.
5. 5
   Свободный член через средние: $b_0=\bar{y}-b_1{\bar{x}}_1-b_2{\bar{x}}_2=7,28-0,1\cdot 50+0,06\cdot 12=3,0$.
6. 6
   Итог: $\hat{y}=3+0,1x_1-0,06x_2$. Контроль: для первой квартиры $3+3-1,2=4,8$ против фактических $4,7$ — остаток $-0,1$.

Читаем коэффициенты. Каждый квадратный метр добавляет к предсказанной цене $0,1$ млн — сто тысяч рублей — при той же удалённости от метро. Каждая минута ходьбы отнимает $0,06$ млн — шестьдесят тысяч — при той же площади. Свободный член $3$ млн — «базовая» цена при нулевых факторах; за пределами данных (квартира меньше 30 м²) его трактовать не стоит. Прогноз для квартиры 55 м² в восьми минутах: $3+0,1\cdot 55-0,06\cdot 8=8,02$ млн.

Почему R² всегда растёт — и чем это плохо?#

Качество подгонки измеряет коэффициент детерминации $R^2$ — доля разброса отклика, объяснённая моделью. Считаем по нашей пятёрке: полный разброс цен вокруг среднего $SST=14,428$, модель объяснила $SSR=b_1{S}_{1y}+b_2{S}_{2y}=0,1\cdot 118+0,06\cdot 37,8=14,068$, в остатках $SSE=0,36$:

Для парной модели только по площади наклон $b_1=S_{1y}/S_{11}=0,118$ и $R^2=13,924/14,428\approx 0,965$. Добавили метро — R² вырос: с 0,965 до 0,975. Так бывает всегда, и это не достижение: новый регрессор физически не может ухудшить подгонку — худший случай, когда МНК присвоит ему нулевой коэффициент и R² не изменится. Модель с лишними переменными выглядит лучше чисто механически, как отчёт, где к успехам команды записали и погоду.

Лекарство — скорректированный R², который штрафует за число факторов:

Считаем: для парной модели ${\bar{R}}^2=1-0,035\cdot 4/3\approx 0,953$, для множественной ${\bar{R}}^2=1-0,025\cdot 4/2\approx 0,950$. Вот и развязка: обычный R² вырос, а скорректированный упал — второй фактор не окупил себя на этой пятёрке данных. Причина видна в следующем разделе: площадь и удалённость в выборке дублируют друг друга.

Что такое мультиколлинеарность и как её поймать?#

В нашей пятёрке большие квартиры стоят ближе к метро — типичная картина новых кварталов. Корреляция факторов: $r_{12}=S_{12}/\sqrt{S_{11}{S}_{22}}=-300/\sqrt{130000}\approx -0,83$. Это и есть мультиколлинеарность — регрессоры несут пересекающуюся информацию. Модель не может решить, кому начислять общий эффект: коэффициенты становятся неустойчивыми — уберите одну квартиру, и $b_2$ перепрыгнет с $-0,06$ на $-0,15$ или сменит знак. Стандартные ошибки растут, t-статистики падают, хотя предсказания страдают меньше всего — модель хорошо описывает данные, но беспомощно разделяет причины.

Стандартный измеритель — фактор инфляции дисперсии. Для каждого регрессора строят вспомогательную регрессию его на все остальные и берут её R²:

Для $x_2$: $R_2^2=(-0,83{)}^2\approx 0,69$, значит $VIF=1/(1-0,69)=3,25$. Классический порог тревоги — 10, в жёстких рекомендациях — 5: до них далеко, но стандартная ошибка $b_2$ уже раздута в $\sqrt{3,25}\approx 1,8$ раза по сравнению с бесколлинеарным миром. Лечение: оставить один из дублирующих факторов, собрать данные, где факторы независимы, или перейти к гребневой регрессии, которая чуть смещает оценку ради устойчивости.

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

Следующий шаг: возьмите реальные данные — десяток квартир и два фактора — и пройдите весь маршрут руками: суммы, система 2×2, коэффициенты, оба R², VIF. Пары «площадь — цена» достаточно, чтобы каждую формулу почувствовать пальцами. Матрицы, если забылись, поджидают в уроке о матрицах и операциях, формулы по статистике — в шпаргалке, тренировка задач — в тренажёре по теории вероятностей. А из сравнения групп начните, если пропустили: дисперсионный анализ — родной брат регрессии, тот же МНК в другой упаковке.

• Модель $\hat{y}=b_0+b_1{x}_1+b_2{x}_2$: каждый коэффициент читается «при прочих равных»
• МНК в матричной форме: $\hat{\beta }=(X^{T}X{)}^{-1}{X}^{T}y$; для двух факторов — система 2×2, решается по Крамеру
• R² не убывает от новых регрессоров; скорректированный R² штрафует за их число и умеет падать
• Мультиколлинеарность ловят корреляцией регрессоров и $VI{F}_j=1/(1-R_j^2)$; лечат выбором факторов или регуляризацией
• Плоскость МНК проходит через точку средних, сумма остатков равна нулю — две бесплатные проверки