§6. Парадоксы регрессии

Существуют три вида лжи: ложь, наглая ложь и статистика.

— Марк Твен

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

1. Неоднородность данных

Существенно исказить вид регрессионной зависимости могут не только выделяющиеся наблюдения («выбросы») по оси отклика $Y$, но и аномальные значения предиктора $X$.

Пример 9. У Лагутина приведены данные о числе телевизионных точек $Y$ (в десятках тысяч), установленных в 1953 г. в девяти городах США (Денвере, Сан-Антонио, Канзас-Сити, Сиэтле, Цинциннати, Буффало, Нью-Орлеане, Милуоки, Хьюстоне), и о численности населения $X$ (в десятках тысяч) этих городов.

Выборочный коэффициент корреляции между наборами $x_1, \dots, x_9$ и $y_1, \dots, y_9$ $\hat\rho = 0{,}403$, что при $n = 9$ свидетельствует о весьма малой степени линейной связи между $X$ и $Y$.

Если же к этим данным добавить соответствующие сведения о Нью-Йорке ($x_{10} = 802$, $y_{10} = 345$), то пересчитанный для $n = 10$ коэффициент $\hat\rho = 0{,}995$. Регрессионная прямая, по сути, соединяет центр масс первых девяти точек $(\bar x, \bar y)$ с $(x_{10}, y_{10})$.

Сами данные по девяти городам в книге приведены только на рисунке, поэтому в виджете — модельные значения, подобранные так, чтобы получились те же 0,403 и 0,995.

Смысл по-пацански: девять мелких городов — облачко без особой связи. Добавили одного гиганта — и «связь» стала почти идеальной. Это не открытие, а одна точка, которая тянет на себя всю прямую. Всегда смотри на картинку, прежде чем верить корреляции.

2. Коррелированность предикторов

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

При подгонке полинома на отрезке $[0, 1]$ (см. пример 2) уже для предикторов $X_3 = U^2$ и $X_4 = U^3$, измеряемых в «узлах» $u_i = i/20$, $i = 0, 1, \dots, 20$, выборочный коэффициент корреляции $\hat\rho$ между соответствующими столбцами матрицы $\boldsymbol X$ составляет 0,986.

Сильная коррелированность предикторов опасна тем, что приводит к неустойчивости МНК-оценок, вычисляемых по формуле (7), к малым возмущениям наблюдений $\eta_i$. Дело в том, что в этом случае столбцы матрицы $\boldsymbol X$ оказываются практически линейно зависимыми, вследствие чего матрица $\boldsymbol B = \boldsymbol X^T\boldsymbol X$ становится почти вырожденной, а задача поиска решения линейной системы (6) — плохо обусловленной.

Рассмотрим подробнее влияние возмущений на интерполяционный полином степени $n - 1$, проходящий через точки плоскости с координатами $(u_i, \eta_i)$, $i = 1, \dots, n$. (Интерполяция — частный случай регрессии при $n = m$.) В форме Лагранжа он выглядит так:

То, что $p(u)$ интерполирует заданные точки, вытекает из равенств $L_i(u_i) = 1$ и $L_i(u_j) = 0$ при $j \ne i$.

У Лагутина на рис. 13 построен полином степени 8, а затем продемонстрированы последствия добавления новой (десятой) точки: график даже не поместился в окне рисунка. Другой пример: даже для полинома степени 6 удаление одной точки приводит к значительным изменениям — возникновению осцилляций.

В регрессионных задачах полиномы высоких степеней имеют тенденцию сглаживать ошибки наблюдений вместо того, чтобы выражать истинный вид зависимости отклика от предиктора. На рис. 15 у Лагутина десять точек, разбросанных вблизи прямой, сглажены МНК-полиномом степени 6. Если зависимость на самом деле достаточно сложная, то для её обнаружения могут оказаться полезными методы непараметрической регрессии, обсуждаемые в гл. 22.

В виджете — своя «галочка» из 9 точек. Подвинь десятую точку так, чтобы она легла рядом с узлом, но не на «галочку»: многочлен начнёт раскачиваться по всему отрезку.

Степень 9 проходит через все точки (RSS ≈ 0), но в $x = 11$ даёт полную чушь: подгонка «выучила» шум.

Интерпретация коэффициентов

Кроме вычислительной неустойчивости, коррелированность предикторов приводит к затруднениям в интерпретации результатов расчётов. Например, при исследовании зависимости веса $Z$ студентов двух групп от их роста $X$ и размера обуви $Y$ методом наименьших квадратов (после предварительного выравнивания масштабов данных) в первой группе было получено регрессионное уравнение

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

Как объяснить существенное различие коэффициентов этих двух моделей?

На практике, подсчитав МНК-оценки $\hat\theta_1, \dots, \hat\theta_m$ в линейной регрессионной модели, исследователь в первую очередь обращает внимание на предикторы с самыми большими (по абсолютной величине) $\hat\theta_j$, так как именно их изменение сильнее всего сказывается на отклике. Исследователь нередко хочет не только точно предсказывать отклик для произвольных значений предикторов, но и желает, задавая эти значения, управлять откликом, надеясь, что регрессия отражает причинно-следственную связь между откликом и предикторами. Понятно, что «нажимать» надо на самые эффективные «рычаги».

В приведённом выше примере для первой группы важнейшим предиктором оказался $X$, а для второй — $Y$. Дело здесь в том, что $X$ и $Y$ сильно коррелируют друг с другом, вследствие чего общий «весовой» коэффициент при $(X - \bar X) + (Y - \bar Y)$ случайным образом распределился между слагаемыми.

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

Смысл по-пацански: рост и размер обуви — почти одно и то же. МНК честно говорит «вместе они весят 1», но как разделить этот вес между ними — ему всё равно, и делит он как попало. Не спрашивай модель о том, чего данные не знают.

3. Неадекватность модели

Когда простейшая линейная зависимость $Y = \theta_1X_1 + \dots + \theta_mX_m$ неадекватно описывает данные, для построения более сложной модели обычно пытаются отдельно изучить влияние каждого предиктора $X_j$ на отклик $Y$. Для этого сглаживают двумерное «облако» точек при помощи некоторой нелинейной функции. Список наиболее часто используемых монотонных функций содержится в следующей таблице.

(В книге [51] у Лагутина содержится любопытная классификация углов из книги по альпинизму, изданной около 1900 г.: «Перпендикулярно — 60°, мой дорогой сэр, абсолютно перпендикулярно — 65°, нависающе — 70°».)

Последняя функция называется логистической кривой. При $a > 0$, $b > 0$ она возрастает, имеет две горизонтальные асимптоты $y = 0$, $y = 1/a$ и перегиб в точке с координатами $(\ln(b/a), 1/(2a))$.

Переход к новым переменным $x'$ и $y'$ сводит задачу к подгонке прямой $y' = a + bx'$ из §1.

Пример 10. В таблице для каждого из годов с 1865 по 1910 ($T$ — номер года) указано количество чугуна $Y$ (в млн тонн), которое выплавлялось за год во всём мире. Постараемся подогнать регрессионную кривую к этим данным.

Положим $X = T - 1864$. Ван дер Варден пишет: «Эта кривая поднимается вверх значительно быстрее, чем прямая линия или квадратная парабола». Не соглашаясь с этим мнением, попытаемся сгладить кривую с помощью параболы. Для подогнанной с помощью МНК параболы остаточная сумма квадратов $RSS = 364{,}7$ и (согласно теореме 1) оценка стандартного отклонения ошибок $\hat\sigma = \sqrt{RSS/(n - m)} = 2{,}91$, где $m = 3$.

Поскольку заметна некоторая тенденция усиления колебаний кривой относительно параболы с ростом $X$, применим, следуя Ван дер Вардену, преобразование $Y' = \ln Y$. Кривая, соединяющая точки с координатами $(X_i, Y_i')$, хорошо согласуется с подогнанной МНК-прямой, имеющей уравнение $y = \hat a + \hat bx$, где $\hat a = 2{,}203$, $\hat b = 0{,}0413$.

Однако, если попытаться использовать последнюю подгонку для предсказания $Y_i$ на основе формулы $\hat Y_i = \exp(\hat a + \hat bX_i)$, то получим значение $RSS = 381{,}8$, которое несколько больше, чем вычисленное ранее при подгонке параболы. Справедливости ради отметим, что сумма модулей остатков для параболы, наоборот, немного больше: $95{,}6 > 93{,}6$. В целом можно считать подгонки примерно одинаковыми по точности.

(Наш пересчёт: с округлёнными $\hat a = 2{,}203$ и $\hat b = 0{,}0413$ получаются ровно книжные 381,8 и 93,6; с неокруглёнными коэффициентами — 379,7 и 93,5.)

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

Экстраполяция

Следует иметь в виду, что в случае ошибки при выборе типа сглаживающей кривой результаты экстраполяции (то есть восстановления значений отклика по значениям предиктора, расположенным вне обследованного диапазона) могут оказаться совершенно неудовлетворительными. Это наглядно демонстрирует рис. 18, на котором зависимость $y = \log_2x$ аппроксимируется полиномом степени 3, интерполирующим точки с координатами $(1, 0)$, $(2, 1)$, $(4, 2)$ и $(16, 4)$.

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

Пример 11. Вес коровы

Рассмотрим в качестве отклика $Z$ вес коровы, а в качестве предикторов — окружность её туловища $X$ и расстояние от хвоста до холки $Y$. Сравнительному анализу были подвергнуты три регрессионные модели:

• (a) линейная: $Z = \theta_1 + \theta_2X + \theta_3Y$;
• (b) степенная: $Z = \theta_1'X^{\theta_2'}Y^{\theta_3'}$;
• (c) учитывающая содержательный смысл задачи: $Z = \theta_0X^2Y$.

Происхождение последней модели легко объяснить. Для этого следует представить себе приближённо тушу коровы в форме цилиндра с длиной образующей, равной $Y$, и радиусом основания, равным $X/(2\pi)$. Если $\rho$ — средняя плотность, то вес $Z = \rho\pi[X/(2\pi)]^2Y = \theta_0X^2Y$ с точностью до головы и ног («рогов и копыт»).

Вначале по всем ($n = 20$) имеющимся наблюдениям для каждой из моделей была вычислена МНК-оценка $\hat{\boldsymbol\theta}$ векторного параметра $\boldsymbol\theta$ и оценка $\hat\sigma = \sqrt{RSS/(n - m)}$ стандартного отклонения ошибок. Результаты расчётов:

Из них как будто следует, что формальные модели 1 и 2 оказались несколько точнее содержательной модели 3. Однако это лишь кажущееся благополучие, что сразу выявляется при перепроверке путём разбиения данных на обучающую и контрольную подвыборки.

В качестве обучающей была взята подвыборка из 10 самых тяжёлых коров, а в качестве контрольной — из 10 оставшихся лёгких коров. На основе обучающей подвыборки для каждой из моделей была заново подсчитана МНК-оценка $\hat{\boldsymbol\theta}_{\text{тяж}}$. Видим, что модели 1 и 2 не выдержали испытание на устойчивость коэффициентов ($\hat\theta_3$ даже поменял знак).

Кроме того, при попытке предсказания веса лёгких коров с помощью модели с такими коэффициентами оценка стандартного отклонения ошибок $\hat\sigma_{\text{лег}}$ увеличилась для формальных моделей более чем в 3 раза!

Этот пример убедительно демонстрирует, что не следует переусложнять модель, ориентируясь на минимизацию $\hat\sigma$: за счёт трёх управляемых параметров («рычагов») удалось «подогнать» формальные модели к данным лучше, чем однопараметрическую содержательную.

Исходных измерений в книге нет, поэтому в виджете — модельное стадо, в котором вес действительно равен $\theta_0X^2Y$ с шумом:

Смысл по-пацански: три «рычага» подгонят что угодно, но не объяснят ничего. Корова — цилиндр, и одна честная формула с одним параметром предсказывает лучше, чем три параметра, подобранные «на глаз» по тяжёлым коровам. Проверяй модель на данных, которые она не видела.

4. Скрытый фактор

Желание истолковывать регрессионную связь как причинно-следственную может приводить к парадоксам, подобным приведённым в двух следующих примерах (см. также близкое к этой теме понятие частной или «очищенной» корреляции из §8 гл. 20).

Пример 12. Во время Второй мировой войны англичане исследовали зависимость точности бомбометания $Z$ от ряда факторов, в число которых входили высота бомбардировщика $H$, скорость ветра $V$, количество истребителей противника $X$. Как и ожидалось, $Z$ увеличивалась при уменьшении $H$ и $V$. Однако (что поначалу представлялось необъяснимым) точность бомбометания возрастала также и при увеличении $X$.

Дальнейший анализ позволил понять причину этого парадокса. Дело оказалось в том, что первоначально в модель не был включён такой важный фактор, как $Y$ — облачность. Он сильно влияет и на $Z$ (уменьшая точность), и на $X$ (бессмысленно высылать истребители, если ничего не видно). Сильные отрицательные причинно-следственные связи в парах $(Y, Z)$ и $(X, Y)$ привели к появлению положительного коэффициента при $X$ в линейной регрессионной модели для $Z$.

Пример 13. Если найти корреляцию между ежегодным количеством родившихся в Голландии детей $Z$ и количеством прилетевших аистов $X$, то она окажется довольно значительной. Можно ли на основе этого статистического результата заключить, что детей приносят аисты?

Рассмотрим проблему на содержательном уровне. Аисты появляются там, где им удобно вить гнёзда; излюбленным же местом их гнездовья являются высокие дымовые трубы, какие строят в голландских сельских домах. По традиции новая семья строит себе новый дом — появляются новые трубы и, естественно, рождаются дети. Таким образом, и увеличение числа гнёзд аистов, и увеличение числа детей являются следствиями одной причины $Y$ — образования новых семей.

Смысл по-пацански: регрессия говорит «вместе меняются», а не «одно вызывает другое». Если забыл про третьего — получишь либо аистов, которые приносят детей, либо истребителей, которые помогают бомбить. Причинность — это не из формулы, это из головы.

Проверь себя