Оценка параметров линейной регрессионной модели вида \(y_i = \beta_0 + \beta_1x_{i1} \dots + \beta_px_{ip} + \epsilon_i\), равно как и выводы в отношении статистической значимости этих параметров, базируются на выполнении ряда математических допущений. Диагностика выполнения этих допущений является составной частью процесса построения регрессионной модели и сводится к следующим составляющим (Faraway 2004):
- проверка допущений в отношении остатков модели;
- проверка адекватности структуры систематической части модели;
- обнаружение необычных наблюдений.
Существуют как графические, так и формальные методы диагностики линейных моделей. Хотя формальные методы используются реже, они доступны в нескольких R-пакетах (см., в частности, car и lmtest). Настоящее сообщение посвящено более распространенным графическим методам. Более того, здесь рассмотрены только первые два из указанных выше пунктов. Обнаружение необычных наблюдений - большая тема, которой будет посвящено отдельное сообщение.
Проверка допущений в отношении остатков модели
Одно из наиболее важных допущений при работе с линейными моделями, параметры которых оцениваются методом наименьших квадратов, состоит в том, что остатки модели независимы (т.е. не коррелируют) и имеют нормальное распределение со средним значением 0 и некоторым фиксированным стандартным отклонением \(\sigma_{\epsilon}\), т.е. \(\epsilon_i \sim N(0, \sigma_{\epsilon}) \). Мы уже сталкивались с этими условиями при рассмотрении моделей дисперсионного анализа (что неудивительно, поскольку дисперсионный анализ является частным случаем линейной модели), и поэтому детально останавливаться на них не будем. Вспомним лишь, что для проверки нормальности распределения количественных переменных можно использовать визуальный анализ гистограмм и квантильных графиков, а также применять формальные методы, вроде теста Шапиро-Уилка (подробнее см. здесь и здесь). Для проверки наличия корреляции в остатках можно использовать графики автокорреляционной функции, а также соответствующие формальные тесты (в частности, тест Дарбина-Уотсона, реализованный, например, в пакете lmtest).
Остановимся несколько подробнее лишь на условии однородности дисперсии остатков модели. Это условие предполагает, что остатки имеют постоянную дисперсию во всем диапазоне предсказанных моделью значений зависимой переменной, а также во всем пространстве значений предикторов. Чаще всего данное условие проверяется графическим способом. Прежде всего, строится график, на котором по оси X откладывают предсказанные моделью значения зависимой переменной, а по оси Y - соответствующие значения остатков. Если условие однородности дисперсии выполняется, то на таком графике точки будут располагаться совершенно случайно относительно горизонтальной линии Y = 0, без каких-либо закономерностей (см. пример для дисперсионного анализа, где это условие не выполняется).
К сожалению, иногда интерпретация подобных диаграмм может вызывать затруднения, особенно при небольших объемах наблюдений. Поэтому для "тренировки глаза" рекомендуется построить несколько калибровочных графиков со случайно сгенерированными данными, свойства которых заведомо известны (Faraway 2004). Примеры таких графиков (по 9 для каждого случая) и код для их построения приведены ниже.
Случай однородной дисперсии остатков:
Случай явно выраженной непостоянной дисперсии остатков (увеличивается при возрастании предсказанных значений зависимой переменной):
Умеренный случай непостоянной дисперсии остатков (увеличивается при возрастании предсказанных значений зависимой переменной):
Помимо построения графиков с предсказанными моделью значениями зависимой переменной, рекомендуется также анализировать графики, где по оси X откладываются значения каждого из включенных и не включенных в модель предикторов. При наличии проблем с неоднородностью дисперсии остатков этот подход позволит выявить "проблемные" предикторы. В зависимости от свойств данных и стоящей перед исследователем задачи, способы устранения неоднородности дисперсии остатков могут включать преобразование исходных значений зависимой переменной и/или предикторов, добавление в модель недостающих предикторов, а также выбора другого, более подходящего ситуации, типа модели (Faraway 2004; Harrell 2002).
Как было отмечено выше, существуют не только графические, но и формальные методы проверки условия однородности дисперсии остатков. В частности, в пакете car имеется функция ncvTest() (от non-constant variance test - "тест на непостоянную дисперсию"), который позволяет проверить нулевую гипотезу о том, что дисперсия остатков никак не связана с предсказанными моделью значениями. Пример использования этой функции здесь не приводится - синтаксис соответствующей команды очень прост (например, ncvTest(model)), равно как и интерпретация результатов теста.
Как было отмечено выше, существуют не только графические, но и формальные методы проверки условия однородности дисперсии остатков. В частности, в пакете car имеется функция ncvTest() (от non-constant variance test - "тест на непостоянную дисперсию"), который позволяет проверить нулевую гипотезу о том, что дисперсия остатков никак не связана с предсказанными моделью значениями. Пример использования этой функции здесь не приводится - синтаксис соответствующей команды очень прост (например, ncvTest(model)), равно как и интерпретация результатов теста.
Проверка адекватности структуры систематической части модели
Как было отмечено ранее, все параметрические статистические модели имеют две части: систематическую и случайную. Так, например, в модели \(y_i = \beta_0 + \beta_1x_{i1} + \epsilon_i\) систематическая часть представлена выражением \([\beta_0 + \beta_1x_{i1}]\), которое определяет математическое ожидание (=среднее значение) зависимой перемененной. В свою очередь, остатки \(\epsilon_i\) отражают наличие случайной вариации наблюдаемых значений зависимой переменной относительно ожидаемого среднего значения.
Предположим, что истинная модель, описывающая процесс генерации значений некоторой переменной \(y\), выглядит следующим образом: \(y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i1}^2 + \epsilon_i\), где \(\beta_0 = 0.1\), \(\beta_1 = 0.2\), \(\beta_2 = 0.035\) и \(\epsilon_i \sim N(0, 10) \). Ниже приведен код, позволяющий сгенерировать 50 значений \(y\) на основе этой модели для \(x = 1, 2, \dots 50\) и изобразить соответствующую зависимость графически:
На практике исследователю никогда не известны ни структура истинной модели, ни значения ее параметров. Все, что он видит - это данные вроде тех, что изображены на рисунке выше. По имеющимся данным бывает очень сложно (а порой и просто невозможно) сделать верное предположение касательно структуры систематической части модели. Так, в рассматриваемом примере, мы могли бы начать с подгонки простой регрессионной модели \(y_i = \beta_0 + \beta_1x_{i1} + \epsilon_i\):
M1 <- lm(y ~ x) summary(M1) Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -21.927 -5.853 1.328 7.870 15.264 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -14.754 2.959 -4.986 8.45e-06 *** x 1.929 0.101 19.103 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 10.31 on 48 degrees of freedom Multiple R-squared: 0.8838, Adjusted R-squared: 0.8813 F-statistic: 364.9 on 1 and 48 DF, p-value: < 2.2e-16
Казалось бы, все указывает на то, что модель M1 хорошо описывает данные: сама модель в целом (Р << 0.05, F-тест) и все ее коэффициенты статистически значимы (Р << 0.05, t-тесты), а скорректированный коэффициент детерминации очень высок (0.881). Графически эта модель представлена на рисунке ниже в виде прямой линии:
Тем не менее, глядя на графическое представление модели M1, можно заметить, что она все же не так хорошо подходит для описания имеющихся данных: видно, что в области низких и высоких значений \(x\) модель "недооценивает" (англ. underfitting) значения \(y\), тогда как в середине диапазона значений \(x\) модель имеет тенденцию к "переоцениванию" (overfitting) значений \(y\). Другими словами, есть основания полагать, что зависимость между этими двумя переменными не является прямолинейной. Это становится особенно очевидным при анализе графика зависимости остатков модели от подогнанных значений - на этом графике остатки образуют U-образную фигуру, что прямо указывает на нелинейный характер зависимости \(y\) от \(x\):
Таким образом, анализ остатков приводит нас к заключению о необходимости изменения структуры систематической части модели. Учесть нелинейных характер зависимости \(y\) от \(x\) это можно несколькими способами. Например, в модель M1 мы можем добавить еще один предиктор, значения которого получены путем возведения в квадрат значений исходного предиктора \(x\):
Модель M2 воспроизводит структуру истинной модели, на основе которой были сгенерированы данные. Как следствие, полученные точечные оценки параметров модели M2 близки к соответствующим истинным значениям: 0.58 против \(\beta_0 = 0.1\), 0.16 против \(\beta_1 = 0.2\), 0.035 против \(\beta_2 = 0.035\) и 8.0 против \(\sigma_{\epsilon} = 10\). Интересно, что оценки коэффициентов \(\beta_0\) и \(\beta_1\) оказались незначимы (Р = 0.871 и P = 0.619 соответственно), тогда как оценка параметра \(\beta_3\) значима. Означает ли этот результат, что предиктор \(x\) можно исключить из модели, оставив только \(x^2\)? Нет. Наблюдаемая ситуация (т.е. незначимый регрессионный коэффициент предиктора \(x\) при значимом коэффициенте производного предиктора - \(x^2\)) вполне может быть обусловлена недостаточно большим числом наблюдений и свойствами данной случайно сгенерированной выборки. В таких случаях рекомендуется оставлять в модели оба предиктора (Hurrell 2002). Кроме того, в целом модель M2 статистически значима (Р << 0.001, F-тест) и имеет более высокий скорректированный коэффициент детерминации, чем M1 (0.928 против 0.881). Наконец, анализ остатков модели M2 показывает, что структура ее систематической части адекватна имеющимся данным (сравните с предыдущим рисунком):
Подчеркнем еще раз, что структура модели M2 совпадает со структурой истинной модели, на основе которой были сгенерированы анализируемые данные. Однако в этом примере мы знали, с чем сравнивать. В реальной жизни мы бы так и не узнали о том, что нам повезло подобрать правильную структуру. Например, для описания этих данных вполне можно было попробовать и другие типы моделей, которые тоже показали бы хорошие результаты (например, модель экспоненциального роста и сплайн-модель). На практике выбор окончательной модели должен подкрепляться хорошим пониманием изучаемого явления, а также использованием соответствующих статистических критериев (в частности, информационным критерием Акаике).
Таким образом, анализ остатков приводит нас к заключению о необходимости изменения структуры систематической части модели. Учесть нелинейных характер зависимости \(y\) от \(x\) это можно несколькими способами. Например, в модель M1 мы можем добавить еще один предиктор, значения которого получены путем возведения в квадрат значений исходного предиктора \(x\):
M2 <- lm(y ~ x + I(x^2)) summary(M2) Call: lm(formula = y ~ x + I(x^2)) Residuals: Min 1Q Median 3Q Max -18.5271 -4.2275 0.4925 4.8585 17.2978 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.579515 3.537067 0.164 0.871 x 0.159959 0.319949 0.500 0.619 I(x^2) 0.034692 0.006082 5.704 7.53e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 8.006 on 47 degrees of freedom Multiple R-squared: 0.9313, Adjusted R-squared: 0.9284 F-statistic: 318.6 on 2 and 47 DF, p-value: < 2.2e-16
Модель M2 воспроизводит структуру истинной модели, на основе которой были сгенерированы данные. Как следствие, полученные точечные оценки параметров модели M2 близки к соответствующим истинным значениям: 0.58 против \(\beta_0 = 0.1\), 0.16 против \(\beta_1 = 0.2\), 0.035 против \(\beta_2 = 0.035\) и 8.0 против \(\sigma_{\epsilon} = 10\). Интересно, что оценки коэффициентов \(\beta_0\) и \(\beta_1\) оказались незначимы (Р = 0.871 и P = 0.619 соответственно), тогда как оценка параметра \(\beta_3\) значима. Означает ли этот результат, что предиктор \(x\) можно исключить из модели, оставив только \(x^2\)? Нет. Наблюдаемая ситуация (т.е. незначимый регрессионный коэффициент предиктора \(x\) при значимом коэффициенте производного предиктора - \(x^2\)) вполне может быть обусловлена недостаточно большим числом наблюдений и свойствами данной случайно сгенерированной выборки. В таких случаях рекомендуется оставлять в модели оба предиктора (Hurrell 2002). Кроме того, в целом модель M2 статистически значима (Р << 0.001, F-тест) и имеет более высокий скорректированный коэффициент детерминации, чем M1 (0.928 против 0.881). Наконец, анализ остатков модели M2 показывает, что структура ее систематической части адекватна имеющимся данным (сравните с предыдущим рисунком):
Подчеркнем еще раз, что структура модели M2 совпадает со структурой истинной модели, на основе которой были сгенерированы анализируемые данные. Однако в этом примере мы знали, с чем сравнивать. В реальной жизни мы бы так и не узнали о том, что нам повезло подобрать правильную структуру. Например, для описания этих данных вполне можно было попробовать и другие типы моделей, которые тоже показали бы хорошие результаты (например, модель экспоненциального роста и сплайн-модель). На практике выбор окончательной модели должен подкрепляться хорошим пониманием изучаемого явления, а также использованием соответствующих статистических критериев (в частности, информационным критерием Акаике).
Встроенные диагностические графики
В рассмотренных выше примерах диагностические графики были построены "вручную". Однако в R имеется и встроенный метод, позволяющий создавать соответствующие графики автоматически. Для реализации этого метода достаточно просто подать модельный объект на функцию plot(). Например:
Слева вверху представлен уже знакомый нам тип графиков: он отражает разброс остатков в зависимости от предсказанных моделью значений и служит для диагностики условия однородности дисперсии. К графику добавлена сглаживающая линия, облегчающая выявление закономерностей в расположении остатков. Примерно горизонтальное расположение этой линии (как в примере выше) указывает на отсутствие каких-либо проблем.
Слева внизу находится очень похожий график - единственное отличие здесь состоит в том, что вместо исходных остатков по оси Y откладываются значения квадратного корня из абсолютных значений стандартизованных остатков (стандартизация выполняется путем деления каждого остатка на стандартное отклонение, рассчитанное по всей совокупности остатков). Стандартизация остатков облегчает выявление неоднородности их дисперии. Сглаживающая линия в данном случае не является горизонтальной. Однако это скорее обусловлено неодинаковой плотностью точек в разных участках диаграммы, чем нарушением условия однородности дисперсии.
Справа вверху представлен график квантилей, позволяющий проверить условие нормальности распределения остатков (подробнее об интерпретации таких графиков см. здесь и здесь).
Справа внизу находится график, предназначенный для выявления необычных наблюдений. Его интерпретация будет описана в отдельном сообщении.
Обратите внимание: на всех четырех графиках есть несколько точек, отмеченных их порядковыми номерами. Этот прием используется для того, чтобы сообщить исследователю о "необычных" свойствах этих конкретных точек. Интерпретация термина "необычный" зависит от контекста (т.е. от того, о каком из приведенных графиков идет речь). Отмеченные точки не обязательно "проблемные" - просто исследователю следует обратить на них внимание.
Слева вверху представлен уже знакомый нам тип графиков: он отражает разброс остатков в зависимости от предсказанных моделью значений и служит для диагностики условия однородности дисперсии. К графику добавлена сглаживающая линия, облегчающая выявление закономерностей в расположении остатков. Примерно горизонтальное расположение этой линии (как в примере выше) указывает на отсутствие каких-либо проблем.
Слева внизу находится очень похожий график - единственное отличие здесь состоит в том, что вместо исходных остатков по оси Y откладываются значения квадратного корня из абсолютных значений стандартизованных остатков (стандартизация выполняется путем деления каждого остатка на стандартное отклонение, рассчитанное по всей совокупности остатков). Стандартизация остатков облегчает выявление неоднородности их дисперии. Сглаживающая линия в данном случае не является горизонтальной. Однако это скорее обусловлено неодинаковой плотностью точек в разных участках диаграммы, чем нарушением условия однородности дисперсии.
Справа вверху представлен график квантилей, позволяющий проверить условие нормальности распределения остатков (подробнее об интерпретации таких графиков см. здесь и здесь).
Справа внизу находится график, предназначенный для выявления необычных наблюдений. Его интерпретация будет описана в отдельном сообщении.
Обратите внимание: на всех четырех графиках есть несколько точек, отмеченных их порядковыми номерами. Этот прием используется для того, чтобы сообщить исследователю о "необычных" свойствах этих конкретных точек. Интерпретация термина "необычный" зависит от контекста (т.е. от того, о каком из приведенных графиков идет речь). Отмеченные точки не обязательно "проблемные" - просто исследователю следует обратить на них внимание.
А что делать с выбросами, которые мы нашли? Избавляться от этих наблюдений и делать расчет без них?
Отправить комментарий