Диагностика линейных регрессионных моделей. Часть 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\)