27 марта 2015

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



Это сообщение является последним по теме диагностики линейных моделей (см. первые две части здесь и здесь). Обсуждение этой темы мы завершим рассмотрением наиболее распространенных количественных показателей, позволяющих обнаружить влиятельные наблюдения. Как следует из их названия, влиятельные наблюдения оказывают существенное воздействие на оценки коэффициентов регрессионной модели. Такие наблюдения обладают как минимум одним из рассмотренных ранее свойств, т.е. являются либо выбросами, либо имеют высокий потенциал воздействия, но чаще всего в определенной степени сочетают в себе оба этих свойства.




DFBETA и DFBETAS

Говоря, что некоторое наблюдение является влиятельным, мы имеем в виду, что оно вносит существенный вклад в оценки параметров модели. Поэтому наиболее простой и очевидный способ измерения степени влиятельности заключается в удалении конкретного наблюдения из выборки и последующем расчете оценок параметров модели без него. Если удаление наблюдения приводит к значительному "скачку" в оценках тех или иных параметров модели, значит это наблюдение является влиятельным. Подобные изменения оценок параметров модели можно записать следующим образом:

\[d_{ij} = b_j - b_{j(-i)},\]

где \(b_{j(-i)}\) обозначает оценку \(j\)-го параметра модели (т.е. \(\beta_j\)), полученную по методу наименьших квадратов после удаления \(i\)-го наблюдения (\(i = 1, \dots, n\), а \(j = 1, \dots, k\)).

Для облегчения интерпретации полезно стандартизовать каждое значение \(d_{ij}\) путем его деления на стандартную ошибку соответствующего коэффициента:

\[d_{ij}^{*} = \frac{d_{ij}}{SE_{(-i)}(b_j)}.\]

Во многих источниках по регрессионному анализу величину \(d_{ij}\) называют еще  \({DFBETA}_{ij}\), а \(d_{ij}^*\) - \({DFBETAS}_{ij}\). Для расчета эти двух показателей в R имеются одноименные функции - dfbeta() и dfbetas(). В качестве примера, используем следующую простую модель (эта модель уже была исследована нами ранее):

set.seed(604)
df = data.frame(x1 = rnorm(50, 10, 1.5), x2 = rnorm(50, 2, 0.1))
 
set.seed(101)
df$y = 2*df$x1 - 25*df$x2 + rnorm(50, 0, 2)
 
M <- lm(y ~ x1 + x2, data = df)

Удобнее всего анализировать величины \({DFBETA}_{ij}\) и \({DFBETAS}_{ij}\) графически. Например, следующим образом можно изобразить стандартизованные показатели влиятельности каждого наблюдения из таблицы df на параметры модели M:

par(mfrow = c(1, 3))
plot(dfbetas(M)[, 1], ylab = "Изменение свободного члена", xlab = "Порядковый номер")
plot(dfbetas(M)[, 2], ylab = "Изменение коэффициента x1", xlab = "Порядковый номер")
plot(dfbetas(M)[, 3], ylab = "Изменение коэффициента x2", xlab = "Порядковый номер")


Из приведенного рисунка видно, что в большинстве случаев значения показателя \({DFBETAS}_{ij}\) близки к нулю. Тем не менее несколько наблюдений демонстрируют повышенную влиятельность. Например, есть два хорошо выделяющихся наблюдения с повышенной влиятельностью на оценки коэффициента предиктора x2.


Расстояние Кука

Cook (1977) предложил выражать влиятельность отдельных наблюдений при помощи следующего показателя:

\[D_{i} = \frac{r_{i}^2}{k+1} \times \frac{h_{ii}}{1-h_{ii}},\]

где \(r_{i}\) - стандартизованный остаток \(i\)-го наблюдения, а \(h_{ii}\) - диагональные элементы матрицы влияния. Как можно видеть из приведенной формулы, первая ее часть отражает степень того, насколько \(i\)-е наблюдение не согласуется с моделью, тогда как вторая часть отражает потенциал воздействия этого наблюдения на параметры модели. Как и в случае с \({DFBETA}_{ij}\) и \({DFBETAS}_{ij}\), большие значения расстояния Кука указывают на высокую влиятельность конкретного наблюдения. Расстояния Кука также принято анализировать графически.

Родственным расстоянию Кука показателем является показатель \({DFFITS}_{i}\) (Belsley et al. 1980):

\[ DFFITS_{i} = t_{i} \sqrt{\frac{h_{ii}}{1 - h_{ii}}},\]

где \(t_{i}\) - стьюдентизированный остаток \(i\)-го наблюдения.

Для расчета расстояний Кука в R служит функция cooks.distance(). Ниже приведены первые 6 расстояний Кука для наблюдений, использованных для подгонки модели M:

head(cooks.distance(M))
           1            2            3            4            5            6 
0.0017350568 0.0038310490 0.0035057901 0.0002404343 0.0010931090 0.0216748540


Ковариационное отношение (COVRATIO)

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

\[COVRATIO_{i} = \frac{1}{(1-h_{ii}) \left( \frac{n - k - 2 - t_{i}^2}{n - k - 1} \right)^{k+1}} \]

Как не сложно догадаться, в R для расчета ковариационных отношений служит функция covratio():

head(covratio(M))
       1        2        3        4        5        6 
1.183041 1.157888 1.067450 1.103658 1.093922 1.045941

Влиятельными считаются наблюдения, чьи значения ковариационного отношения значительно превышают 1.


Функция influence.measures()

Базовая R-функция influence.measures() позволяет одновременно рассчитать все перечисленные выше показатели влиятельности. Так, для нашего примера с моделью M получаем:

influence.measures(M)
Influence measures of
  lm(formula = y ~ x1 + x2, data = df) :
 
     dfb.1_    dfb.x1    dfb.x2    dffit cov.r   cook.d    hat inf
1   0.01813 -0.063954  0.002171  0.07141 1.183 1.74e-03 0.1010    
2  -0.02825  0.092827  0.003114  0.10620 1.158 3.83e-03 0.0860    
3   0.02172  0.032939 -0.036537 -0.10189 1.067 3.51e-03 0.0258    
4  -0.00862  0.016489  0.004771  0.02658 1.104 2.40e-04 0.0348    
5  -0.00838  0.033999  0.000268  0.05671 1.094 1.09e-03 0.0312    
6  -0.19156  0.089453  0.178146  0.25542 1.046 2.17e-02 0.0534    
7  -0.03018  0.031223  0.026096  0.10772 1.059 3.91e-03 0.0235    
8   0.09911  0.014119 -0.109647 -0.12430 1.165 5.25e-03 0.0932    
9   0.04351  0.184071 -0.094762  0.25601 1.050 2.18e-02 0.0552    
10  0.10655  0.030397 -0.122755 -0.14224 1.158 6.86e-03 0.0907    
11 -0.01049  0.060284 -0.003977  0.10205 1.077 3.52e-03 0.0307    
12  0.04534 -0.040584 -0.040190 -0.14030 1.038 6.59e-03 0.0241    
13 -0.13749 -0.055091  0.169883  0.29641 0.926 2.82e-02 0.0311    
14 -0.06980 -0.175889  0.116252 -0.30914 0.943 3.09e-02 0.0365    
15 -0.00290 -0.000487  0.002820 -0.00814 1.091 2.26e-05 0.0228    
16 -0.01789 -0.036047  0.028655 -0.05023 1.194 8.59e-04 0.1080   *
17 -0.11707 -0.036494  0.129081 -0.15562 1.130 8.19e-03 0.0744    
18  0.18590 -0.056689 -0.170842  0.21033 1.120 1.49e-02 0.0791    
19 -0.00872  0.008001  0.001914 -0.10741 1.050 3.88e-03 0.0201    
20  0.54367  0.124829 -0.619317 -0.74046 0.745 1.62e-01 0.0713   *
21  0.02583 -0.018782 -0.022792 -0.05013 1.095 8.54e-04 0.0312    
22  0.41428 -0.318250 -0.321729  0.51923 0.995 8.69e-02 0.0930    
23  0.01982  0.012162 -0.025866 -0.04588 1.098 7.16e-04 0.0322    
24 -0.24542  0.209544  0.185135 -0.30858 1.180 3.20e-02 0.1318    
25 -0.14183 -0.333578  0.256764  0.44853 1.140 6.67e-02 0.1372    
26  0.07147  0.040907 -0.096179 -0.24098 0.944 1.88e-02 0.0245    
27 -0.02902  0.052461  0.015710  0.06724 1.132 1.54e-03 0.0617    
28  0.02142  0.004683 -0.024778 -0.03594 1.108 4.40e-04 0.0390    
29  0.05729  0.042153 -0.067800  0.13023 1.065 5.71e-03 0.0315    
30 -0.06019 -0.001737  0.065973  0.09526 1.093 3.08e-03 0.0384    
31 -0.15968 -0.169586  0.224534  0.31968 1.076 3.39e-02 0.0805    
32  0.17321 -0.043257 -0.161027  0.20847 1.080 1.46e-02 0.0573    
33  0.12868 -0.136507 -0.082130  0.27150 0.951 2.39e-02 0.0310    
34  0.18868 -0.022618 -0.204170 -0.40772 0.758 5.01e-02 0.0269   *
35  0.07380 -0.028323 -0.058006  0.23506 0.931 1.79e-02 0.0217    
36  0.08553  0.010526 -0.097228 -0.15532 1.054 8.09e-03 0.0330    
37 -0.01878 -0.024655  0.028898  0.05667 1.100 1.09e-03 0.0354    
38 -0.06177 -0.241186  0.146278  0.33724 1.020 3.74e-02 0.0626    
39 -0.28845  0.069363  0.269177 -0.34481 1.004 3.89e-02 0.0589    
40  0.01576  0.127138 -0.052750  0.14792 1.216 7.43e-03 0.1308   *
41 -0.01879  0.068312  0.001283  0.09228 1.103 2.89e-03 0.0444    
42 -0.05374  0.077573  0.036675  0.13795 1.064 6.40e-03 0.0332    
43 -0.67300 -0.011346  0.687573 -0.75845 0.906 1.78e-01 0.1124   *
44  0.22057 -0.143563 -0.181125  0.24959 1.267 2.11e-02 0.1744   *
45 -0.02216  0.168941 -0.034342 -0.21376 1.075 1.53e-02 0.0560    
46 -0.00901  0.059482 -0.007638  0.06602 1.197 1.48e-03 0.1106   *
47 -0.04181  0.070595  0.023694  0.08739 1.149 2.60e-03 0.0772    
48  0.05957  0.050523 -0.081842 -0.14125 1.070 6.72e-03 0.0364    
49  0.02891 -0.017702 -0.023437  0.03949 1.116 5.31e-04 0.0466    
50 -0.20345  0.148449  0.158515 -0.27342 1.046 2.48e-02 0.0574    

В первых трех столбцах приведены DFBETAS для свободного члена и двух коэффициентов регрессионного уравнения. Далее следуют столбцы с DFFIT, COVRATIO, расстояниями Кука и показателями потенциала воздействия. Замыкает таблицу столбец inf - он содержит звездочки * напротив наблюдений, которые по совокупности всех показателей следует считать влиятельными.


Пороговые значения показателей влиятельности

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

  • \( |d_{ij}^{*}|  \geq 1\) или \( |d_{ij}^{*}|  \geq 2\)
  • \( D_i  > 4/(n-k-1)\)
  • \( DFFITS_{i} > 2 \sqrt{(k+1)(n-k-1)}\)
  • \( COVRATIO_{i} > 1\) или \( |COVRATIO_{i} -1 | > 3(k+1)/n\)

6 комментариев :

Alexandr комментирует...

Здравствуйте, спасибо за интересную статью. Непонятна фраза "т.е. являются либо выбросами, либо имеют высокий потенциал воздействия, но чаще всего в определенной степени сочетают в себе оба этих свойства". Как можно идентифицировать с каким случаем мы имеем дело и каковы дальнейшие действия? Если это выброс, то наверное его нужно исключить из набора данных (?). А если это наблюлюдение с высоким потенциалом действия, или наблюдение "сочетает в себе оба этих свойства" (кстати, как идентифицировать этот случай?), то что делать в этом случае? Оставить модель как есть, не переучивая её на данных без этих наблюдений?

Sergey Mastitsky комментирует...

Фраза, которая вызвала у Вас затруднение, означает только то, что означает: не каждое потенциально влиятельное наблюдение является выборосом, и не каждый выброс обязательно оказывает существенно воздействие на оценки параметров модели. Но когда эти два свойства сочетаются, тогда мы имеем дело с влиятельными наблюдениями. Советую перечитать еще раз предыдущие два сообщения по этой теме, особенно Часть 2, где приводится подробное описание понятий "выброс" и "потенциально влиятельное наблюдение". Кроме того, хорошие примеры можно найти в лит. источниках, упомянутых во второй части (в частности, попытайтесь раздобыть книгу Fox 2010).

Sergey Mastitsky комментирует...

В Части 3 описаны показатели, которые помогают учесть "сочетание обоих свойств". Насчет того, что делать с такими влиятельными наблюдениями: "сколько людей - столько мнений". Общий совет, я бы сказал таков: опишите эти необычные наблюдения в своей работе и приведите альтернативные модели - содержащие и не содержащие такие наблюдения. Так читатель будет видеть, наскольк велико влияние этих наблюдений на ту или иную модель и на получаемые с ее помощью выводы.

Alexandr комментирует...

После прочтения дискуссии на примерно такой же вопрос как возник у меня на stackoverflow:
http://stats.stackexchange.com/questions/46826/influential-residual-vs-outlier

нашёл ясное определение этих понятий тут:
http://stattrek.com/regression/influential-points.aspx?Tutorial=AP

Правда определение которое даётся там, отличается от того что Вы написали в первом комментарии.
Там написано что "An influential point is an _outlier_ that greatly affects the slope of the regression line.".

Т.е. не каждый outlier оказывает существенное влияние на оценку параметров, но каждая influential point является outlier-ом.

Sergey Mastitsky комментирует...

Alexandr, спасибо за ссылки. Да, мнения можно встретить разные. Свои сообщения я почти полностью строил по работам Fox'а, который дядька авторитетный в этих делах.

Alexandr комментирует...

Да, похоже определения в разных книгах разные, и иногда противоречащие друг другу в деталях. Я хоть и читал книги про регрессию, но такого понятия не встречал, а оно оказывается есть. Спасибо Вам за интересную статью!

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