14 января 2015

Диагностика линейных регрессионных моделей. Часть 2



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

Говоря о необычных наблюдениях в контексте регрессионного анализа, можно выделить следующие три ситуации:
  • Наблюдение представлено необычным сочетанием значений предикторов (англ. 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):

set.seed(604)
df = data.frame(x1 = rnorm(50, 10, 1.5), x2 = rnorm(50, 2, 0.1))
 
head(df)
        x1       x2
1  7.235169 1.961452
2 12.788041 1.971398
3  9.640045 1.995138
4 11.353070 1.983749
5 11.242469 1.963742
6 11.038517 2.060605

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

set.seed(101)
df$y = 2*df$x1 - 25*df$x2 + rnorm(50, 0, 2)

По имеющимся данным мы можем восстановить регрессионную модель при помощи функции 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() - функция, выполняющая транспонирование матриц):

H <- X %*% solve( t(X) %*% X ) %*% t(X)

Для расчета предсказываемых моделью значений y воспользуемся приведенной выше формулой:

Yhat <- H %*% df$y

Легко проверить, что полученный таким образом вектор предсказанных значений \(\boldsymbol{\hat{y}}\) идентичен вектору, который можно было бы получить более удобным способом при помощи стандартной функции fitted():

identical(round(fitted(M), 2), round(H %*% df$y, 2)[, 1])
[1] TRUE

Диагональные элементы матрицы проекции изменяются от 0 до 1 и отражают потенциал воздействия отдельных наблюдений на оценки регрессионных коэффициентов. Чем дальше то или иное наблюдение находится от центра многомерного распределения значений предикторов, тем выше будет соответствующий диагональный элемент. Эти элементы можно получить несколькими эквивалентными способами:

hat(X# Способ 1
diag(H) # Способ 2
# Способ 3: Minf <- influence(M) Minf$hat

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

sum(hat(X))
[1] 3

Соответственно, среднее значение \(h_{ii}\) можно рассчитать как \(p/n\). Отсюда вытекает эмпирическое правило, позволяющее судить о том, оказывает ли некоторое наблюдение существенное влияние на параметры модели - значения \(h_{ii} > 2p/n\) являются достаточно большими, чтобы считать соответствующие наблюдения стоящими внимания (Fox 2010).

Значения \(h_{ii}\) удобнее всего анализировать графически. Так, для нашего примера имеем (рис. 2):

plot(hat(X), type = "h")
abline(h = 2*3/50, lty = 2, col = "blue")

Рис. 2. Показатели потенциала воздействия для приведенной выше модели M. Согласно пороговому значению \(2p/n\) (показано прерывистой горизонтальной линией), четыре из 50 имеющихся наблюдений имеют высокий потенциал воздействия на коэффициенты этой модели.

Согласно полученному результату, четыре из 50 имеющихся наблюдений (под номерами 24, 25, 40 и 44), имеют высокий потенциал воздействия на коэффициенты модели M и требуют дополнительного изучения.


Выбросы

В контексте регрессионного анализа выбросом называют наблюдение, которое на согласуется с рассматриваемой моделью. Другими словами, выброс - это наблюдение с большим остатком, возникающим из-за того, что соответствующее выборочное значение зависимой переменной \(y_i\) значительно отличается от предсказанного значения \(\hat{y}_{i}\).

Для выявления выбросов можно использовать уже рассмотренные нами ранее диагностические графики, на которых по оси X откладывают предсказываемые моделью значения, а по оси Y - остатки. Пример такого графика приведен на рис. 3.

plot(fitted(M), resid(M))
abline(h = 0, lty = 2)

Рис. 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()):

r <- residuals(M)/(summary(M)$sig*sqrt(1-hat(X)))
plot(fitted(M), r, ylim = c(-3, 3))
abline(h=0, lty = 2)
abline(h = c(-2, 2), lty = 2, col = "blue")

Рис. 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 1978Kutner et al. 2004Fox 2010).

Стьюдентизированные остатки имеют t-распределение с \(n - p - 1\) степенями свободы. Соответственно, мы можем использовать квантили этого распределения для проверки того, насколько статистически значимо определенное наблюдение является выбросом. Так, в случае с нашим примером, мы можем проверить, является ли статистически значимым выбросом наблюдение с наибольшим (абсолютным) значением стьюдентизированного остатка (в R для расчета стьюдентизированных остатков служит функция rstudent()):

max(abs(rstudent(M)))
[1] 2.671676
 
abs(qt(.05/(50*2), 50-3-1)) # с поправкой Бонферрони
[1] 3.514957

Как видим, максимальное наблюдаемое значение стьюдентизированного остатка (2.672) не превышает критическое t-значение = 3.515 (рассчитанное с применением поправки Бонферрони). Из этого следует вывод, что данное наблюдение (№20) не является статистически значимым выбросом. Аналогичный вывод справедлив и для других наблюдений (чьи стьюдентизированные остатки еще меньше, чем 3.515).

Таким образом, применяя стьюдентизированные остатки, мы не обнаружили каких-либо наблюдений, явно не согласующихся с моделью M. Для сравнения отметим, что в выполненном ранее анализе мы выявили три наблюдения с высоким потенциалом воздействия на параметры модели. Такая ситуация обычна - наблюдения с высоким потенциалом воздействия необязательно являются выбросами (и наоборот).

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


Комментариев нет :

Отправить комментарий