Продолжая начатую ранее тему диагностики линейных регрессионных моделей, рассмотрим некоторые распространенные методы выявления необычных наблюдений.
Говоря о необычных наблюдениях в контексте регрессионного анализа, можно выделить следующие три ситуации:
Говоря о необычных наблюдениях в контексте регрессионного анализа, можно выделить следующие три ситуации:
- Наблюдение представлено необычным сочетанием значений предикторов (англ. leverage point).
- Наблюдение не согласуется с рассматриваемой моделью, т.е. является выбросом (англ. outlier).
- Наблюдение оказывает существенное влияние на оценки параметров модели (англ. influential point или influential observation). Другими словами, удаление такого влиятельного наблюдения из выборки приведет к значительному изменению предсказываемых моделью значений. Влиятельные наблюдения обладают как минимум одним из двух указанных выше свойств (т.е. являются либо "leverage point", либо "outlier"), но чаще всего сочетают их.
Потенциал воздействия отдельных наблюдений на параметры модели
Наблюдения, представленные необычным сочетанием предикторов и, как следствие, расположенные далеко от центра (многомерного) распределения наблюдаемых значений предикторов, потенциально могут оказывать существенное влияние на оценки параметров модели (рис. 1). О таких наблюдениях говорят, что они обладают "высоким потенциалом воздействия" (англ. high leverage) на параметры модели.
 |
| Рис. 1. Пример ядерной оценки плотности вероятности двумерного распределения. Примерные координаты центра этого распределения: (10, 2). Видно, что некоторые точки находятся достаточно далеко от центра - например, точка, выделенная желтым цветом. Это необычное для данного распределения наблюдение потенциально (но необязательно) может оказывать существенное влияние на оценки параметров линейной модели, включающей предикторы x1 и x2. |
Имеется возможность выразить этот потенциал воздействия количественно. Распространенным подходом является расчет диагональных элементов т.н. матрицы проекции на пространство предикторов (англ. hat matrix; другие названия - "матрица влияния" и "матрица воздействия" (influence matrix)), которые обозначаются как \(h_{ii}\). Английский термин "hat matrix" происходит от названия значка ^ (hat - "шапка"), который обычно используется в обозначении предсказываемых моделью значений зависимой переменной - \(\hat{y}\).
Чтобы понять, откуда берется матрица проекции и что она собой представляет, необходимо прибегнуть к записи линейных моделей с использованием обозначений линейной алгебры. В общем виде мы можем представить линейную модель следующим образом:
\[\boldsymbol{y = X\beta + \epsilon},\]
где \(\boldsymbol{X}\) - матрица модели (размером \(n\times p\), где \(p\) - число регрессионных коэффициентов + свободный член уравнения), \(\boldsymbol{\beta}\) - вектор подлежащих оценке регрессионных коэффициентов, а \(\boldsymbol{\epsilon}\) - вектор остатков с нулевым математическим ожиданием.
Вектор с оценками регрессионных коэффициентов \(\boldsymbol{\hat{\beta}}\) получают следующим образом:
\[\boldsymbol{\hat{\beta} = (X^{T}X)^{-1}X^{T}y,}\]
откуда предсказываемые моделью значения зависимой переменной можно записать как
\[\boldsymbol{\hat{y} = X\hat{\beta} = X(X^{T}X)^{-1}X^{T}y}\]
или, в более простой форме, как
или, в более простой форме, как
\[\boldsymbol{\hat{y} = Hy}.\]
Матрица \(\boldsymbol{H}\) (размером \(n \times n\)) как раз и является матрицей проекции, поскольку она позволяет выполнить линейное преобразование наблюдаемых значений \(\boldsymbol{y}\) так, чтобы получить предсказываемые значения (подробнее см., например, Hoaglin and Welsch 1978).
Рассмотрим на примере некоторые полезные свойства матрицы проекции. Предположим, у нас имеются два предиктора x1 и x2 (см. также рис. 1):
Предположим далее, что переменные x1 и x2 связаны с некоторой переменной y следующей зависимостью (подробнее о способах записи линейных моделей и генерации данных на их основе см. здесь):
По имеющимся данным мы можем восстановить регрессионную модель при помощи функции lm():
> M <- lm(y ~ x1 + x2, data = df) > summary(M) Call: lm(formula = y ~ x1 + x2, data = df) Residuals: Min 1Q Median 3Q Max -4.3855 -1.2553 0.3463 1.2528 2.8974 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -12.1994 6.0746 -2.008 0.0504 . x1 2.2655 0.1766 12.827 < 2e-16 *** x2 -20.2806 2.9861 -6.792 1.69e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.811 on 47 degrees of freedom Multiple R-squared: 0.8137, Adjusted R-squared: 0.8057 F-statistic: 102.6 on 2 and 47 DF, p-value: < 2.2e-16
Для расчета матрицы проекции сначала сформируем матрицу модели, для чего воспользуемся базовой R-функцией model.matrix():
X = model.matrix(y ~ x1 + x2, data = df) head(X) (Intercept) x1 x2 1 1 7.235169 1.961452 2 1 12.788041 1.971398 3 1 9.640045 1.995138 4 1 11.353070 1.983749 5 1 11.242469 1.963742 6 1 11.038517 2.060605
Получить \(\boldsymbol{H}\) теперь очень просто (в приведенной ниже команде %*% - это оператор умножения матриц, solve() - функция, выполняющая обращение матриц, а t() - функция, выполняющая транспонирование матриц):
Для расчета предсказываемых моделью значений y воспользуемся приведенной выше формулой:
Yhat <- H %*% df$y
Легко проверить, что полученный таким образом вектор предсказанных значений \(\boldsymbol{\hat{y}}\) идентичен вектору, который можно было бы получить более удобным способом при помощи стандартной функции fitted():
Диагональные элементы матрицы проекции изменяются от 0 до 1 и отражают потенциал воздействия отдельных наблюдений на оценки регрессионных коэффициентов. Чем дальше то или иное наблюдение находится от центра многомерного распределения значений предикторов, тем выше будет соответствующий диагональный элемент. Эти элементы можно получить несколькими эквивалентными способами:
Оказывается, что сумма диагональных элементов матрицы проекции равна числу коэффициентов регрессионного уравнения, включая свободный член:
Соответственно, среднее значение \(h_{ii}\) можно рассчитать как \(p/n\). Отсюда вытекает эмпирическое правило, позволяющее судить о том, оказывает ли некоторое наблюдение существенное влияние на параметры модели - значения \(h_{ii} > 2p/n\) являются достаточно большими, чтобы считать соответствующие наблюдения стоящими внимания (Fox 2010).
Значения \(h_{ii}\) удобнее всего анализировать графически. Так, для нашего примера имеем (рис. 2):
 |
| Рис. 2. Показатели потенциала воздействия для приведенной выше модели M. Согласно пороговому значению \(2p/n\) (показано прерывистой горизонтальной линией), четыре из 50 имеющихся наблюдений имеют высокий потенциал воздействия на коэффициенты этой модели. |
Выбросы
В контексте регрессионного анализа выбросом называют наблюдение, которое на согласуется с рассматриваемой моделью. Другими словами, выброс - это наблюдение с большим остатком, возникающим из-за того, что соответствующее выборочное значение зависимой переменной \(y_i\) значительно отличается от предсказанного значения \(\hat{y}_{i}\).
Для выявления выбросов можно использовать уже рассмотренные нами ранее диагностические графики, на которых по оси X откладывают предсказываемые моделью значения, а по оси Y - остатки. Пример такого графика приведен на рис. 3.
 |
| Рис. 3. График остатков модели М. |
На практике, однако, при работе с такими "сырыми" значениями остатков (англ. raw residuals) бывает трудно сказать, какие из них достаточно велики. Кроме того, значения "сырых" остатков зависят от шкалы, по которой измеряется зависимая переменная. Для устранения этих недостатков применяют разного рода стандартизацию. Следует отметить, что терминология, используемая в разных источниках для обозначения таких преобразованных остатков, крайне запутанна. В R для диагностики линейных моделей, чьи параметры оцениваются по методу наименьших квадратов, используются следующие два типа остатков.
Standardized residuals - собственно стандартизованные остатки:
\[r_{i} = \frac{\epsilon_i}{S_{\epsilon}\sqrt{1-h_{ii}}},\]
где \(\epsilon_i\) - остаток \(i\)-го наблюдения, \(S_{\epsilon}\) - стандартное отклонение всех остатков модели, а \(h_{ii}\) - уже знакомый нам показатель потенциала воздействия \(i\)-го наблюдения на коэффициенты модели.
Такие остатки имеют приближенно стандартное нормальное распределение. Соответственно, наблюдения, чьи стандартизованные остатки выходят за пределы диапазона от -2 до 2, считают выбросами. Для нашего примера получаем (заметьте, что для расчета r ниже мы могли бы применить также базовую R-функцию rstandard()):
Как видно из рис. 4, наблюдения 20, 34 и (с натяжкой) 43 можно было бы считать выбросами.
Одним из важных недостатков стандартизованных остатков является тот факт, что при расчете \(i\)-го значения \(r\) используется оценка \(S_{\epsilon}\), которая включает значение соответствующего "сырого" остатка. Другими словами, любое значение \(r_i\) и \(S_{\epsilon}\) не являются независимыми, затрудняя формальную проверку статистической гипотезы о том, что некоторое \(i\)-е наблюдение не является выбросом.
Для устранения указанного недостатка используют Studentized residuals - стьюдентизированные остатки:
\[t_{i} = \frac{\epsilon_i}{S_{\epsilon_{(-i)}}\sqrt{1-h_{ii}}},\]
Главным отличием этой формулы от предыдущей является наличие индекса \((-i)\) в обозначении стандартного отклонения остатков \(S_{\epsilon_{(-i)}}\). Этот индекс указывает на то, что стандартное отклонение рассчитывается по остаткам модели, подогнанной после исключения из данных \(i\)-го наблюдения (т.е., по сути, используется техника jackknife). Поскольку мы рассчитываем как бы "внешнюю" (по отношению к \(i\)-му наблюдению) оценку стандартного отклонения (\(i\)-е наблюдение исключается из расчета и не влияет на эту оценку), такие остатки известны также под названием externally studentized residuals (внешне стьюдентизированные остатки). Соответственно, то, что мы выше определили как стандартизованные остатки, в ряде руководств называют internally studentized residuals (внутренне стьюдентизированные остатки) (см., например, Hoaglin and Welsch 1978; Kutner et al. 2004; Fox 2010).
Стьюдентизированные остатки имеют t-распределение с \(n - p - 1\) степенями свободы. Соответственно, мы можем использовать квантили этого распределения для проверки того, насколько статистически значимо определенное наблюдение является выбросом. Так, в случае с нашим примером, мы можем проверить, является ли статистически значимым выбросом наблюдение с наибольшим (абсолютным) значением стьюдентизированного остатка (в R для расчета стьюдентизированных остатков служит функция rstudent()):
Как видим, максимальное наблюдаемое значение стьюдентизированного остатка (2.672) не превышает критическое t-значение = 3.515 (рассчитанное с применением поправки Бонферрони). Из этого следует вывод, что данное наблюдение (№20) не является статистически значимым выбросом. Аналогичный вывод справедлив и для других наблюдений (чьи стьюдентизированные остатки еще меньше, чем 3.515).
Таким образом, применяя стьюдентизированные остатки, мы не обнаружили каких-либо наблюдений, явно не согласующихся с моделью M. Для сравнения отметим, что в выполненном ранее анализе мы выявили три наблюдения с высоким потенциалом воздействия на параметры модели. Такая ситуация обычна - наблюдения с высоким потенциалом воздействия необязательно являются выбросами (и наоборот).
Выполнив проверку на наличие наблюдений с необычными сочетаниями значений предикторов и наблюдений, которые не согласуются с рассматриваемой моделью, обычно приступают к обобщению полученной информации в виде количественных показателей, отражающих собственно степень влияния каждого наблюдения на параметры модели. Этим показателям будет посвящена следующая, 3-я часть этой статьи.
Standardized residuals - собственно стандартизованные остатки:
\[r_{i} = \frac{\epsilon_i}{S_{\epsilon}\sqrt{1-h_{ii}}},\]
где \(\epsilon_i\) - остаток \(i\)-го наблюдения, \(S_{\epsilon}\) - стандартное отклонение всех остатков модели, а \(h_{ii}\) - уже знакомый нам показатель потенциала воздействия \(i\)-го наблюдения на коэффициенты модели.
Такие остатки имеют приближенно стандартное нормальное распределение. Соответственно, наблюдения, чьи стандартизованные остатки выходят за пределы диапазона от -2 до 2, считают выбросами. Для нашего примера получаем (заметьте, что для расчета r ниже мы могли бы применить также базовую R-функцию rstandard()):
 |
| Рис. 4. График стандартизованных (=внутренне стьюдентизированных) остатков модели M. Согласно описанному в тексте эмпирическому правилу, наблюдения 20, 34 и 43 можно считать выбросами. С такими остатками мы уже сталкивались в первой части этой статьи. |
Как видно из рис. 4, наблюдения 20, 34 и (с натяжкой) 43 можно было бы считать выбросами.
Одним из важных недостатков стандартизованных остатков является тот факт, что при расчете \(i\)-го значения \(r\) используется оценка \(S_{\epsilon}\), которая включает значение соответствующего "сырого" остатка. Другими словами, любое значение \(r_i\) и \(S_{\epsilon}\) не являются независимыми, затрудняя формальную проверку статистической гипотезы о том, что некоторое \(i\)-е наблюдение не является выбросом.
Для устранения указанного недостатка используют Studentized residuals - стьюдентизированные остатки:
\[t_{i} = \frac{\epsilon_i}{S_{\epsilon_{(-i)}}\sqrt{1-h_{ii}}},\]
Главным отличием этой формулы от предыдущей является наличие индекса \((-i)\) в обозначении стандартного отклонения остатков \(S_{\epsilon_{(-i)}}\). Этот индекс указывает на то, что стандартное отклонение рассчитывается по остаткам модели, подогнанной после исключения из данных \(i\)-го наблюдения (т.е., по сути, используется техника jackknife). Поскольку мы рассчитываем как бы "внешнюю" (по отношению к \(i\)-му наблюдению) оценку стандартного отклонения (\(i\)-е наблюдение исключается из расчета и не влияет на эту оценку), такие остатки известны также под названием externally studentized residuals (внешне стьюдентизированные остатки). Соответственно, то, что мы выше определили как стандартизованные остатки, в ряде руководств называют internally studentized residuals (внутренне стьюдентизированные остатки) (см., например, Hoaglin and Welsch 1978; Kutner et al. 2004; Fox 2010).
Стьюдентизированные остатки имеют t-распределение с \(n - p - 1\) степенями свободы. Соответственно, мы можем использовать квантили этого распределения для проверки того, насколько статистически значимо определенное наблюдение является выбросом. Так, в случае с нашим примером, мы можем проверить, является ли статистически значимым выбросом наблюдение с наибольшим (абсолютным) значением стьюдентизированного остатка (в R для расчета стьюдентизированных остатков служит функция rstudent()):
Как видим, максимальное наблюдаемое значение стьюдентизированного остатка (2.672) не превышает критическое t-значение = 3.515 (рассчитанное с применением поправки Бонферрони). Из этого следует вывод, что данное наблюдение (№20) не является статистически значимым выбросом. Аналогичный вывод справедлив и для других наблюдений (чьи стьюдентизированные остатки еще меньше, чем 3.515).
Таким образом, применяя стьюдентизированные остатки, мы не обнаружили каких-либо наблюдений, явно не согласующихся с моделью M. Для сравнения отметим, что в выполненном ранее анализе мы выявили три наблюдения с высоким потенциалом воздействия на параметры модели. Такая ситуация обычна - наблюдения с высоким потенциалом воздействия необязательно являются выбросами (и наоборот).
Выполнив проверку на наличие наблюдений с необычными сочетаниями значений предикторов и наблюдений, которые не согласуются с рассматриваемой моделью, обычно приступают к обобщению полученной информации в виде количественных показателей, отражающих собственно степень влияния каждого наблюдения на параметры модели. Этим показателям будет посвящена следующая, 3-я часть этой статьи.
Отправить комментарий