Revista Colombiana de Estad´ıstica Volumen 27 No 2. P´ ags. 109 a 121. Diciembre 2004
Un criterio para identificar datos at´ıpicos ´ Alfredo Jime ´nez Moscoso* Jose
Resumen En este art´ıculo se presenta un m´etodo para determinar las observaciones que son at´ıpicas en un modelo de regresi´ on lineal m´ ultiple; estos datos se estableceran de acuerdo al cambio que ejercen sobre la suma de los cuadrados de residuales del modelo. Palabras Claves: Modelos lineales, m´ınimos cuadrados, formas cuadr´ aticas, observaciones at´ıpicas, estad´ıstica Qk .
Abstract This paper present a method to determine the observations that are outliers in a model of multiple linear regression; these data will be established according to the change that is presented on the sum of the squares of residual of the model. Key words: Linear models, Least squares, Quadratic forms, Outliers, Qk Statistics.
1.
Introducci´ on
Draper & John (1981) proponen una metodolog´ıa para detectar un grupo de k observaciones at´ıpicas, an´aloga a la propuesta de Bartlett (1937), citada en Little & Rubin (1987), para estimar los par´ametros del modelo de regresi´on lineal cuando existen observaciones faltantes en la variable respuesta. En el * Profesor asistente, Universidad Nacional de Colombia, Departamento de Matem´ aticas. E-mail:
[email protected]
109
110
Jos´e A. Jim´enez M.
planteamiento de Draper & John (1981) se considera el modelo de regresi´on lineal m´ ultiple: Y = X
n×1
β + ,
(1)
β + 1 , γ 2
(2)
n×r r×1
n×1
particionado de la siguiente manera: Y1 X1 = Y2 X2
I 0
donde Y1 es el bloque conformado por las observaciones consideradas at´ıpicas. Para el modelo (2) establecen las estimaciones de β y γ mediante: βˆ =(X20 X2 )−1 X20 Y2 ,
γˆ = (I − H11 )−1 ˆ1 ,
donde Hij = Xi (X 0 X)−1 Xj0 es una submatriz de la matriz
H = X(X 0 X)−1 X 0 ,
para X =
X1 . X2
La notaci´ on de H y el nombre de matriz hat fue introducido por Tukey (1977); por otra parte, el cambio en la suma de cuadrados de residuales lo calculan usando la estad´ıstica: Qk =ˆ 01 (I − H11 )−1 ˆ1 ,
con k = dim (Y1 ) .
(3)
En resumen, el m´etodo descrito permite detectar el grupo de observaciones at´ıpicas en base al cambio en la suma de cuadrados de residuales, lo cual se cuantifica con la estad´ıstica Qk , es decir, mediante este procedimiento se selecciona el bloque Y1 que posee el Qk m´as alto, como el bloque m´as at´ıpico, y en muchos casos quedan datos at´ıpicos dentro de un bloque y el m´etodo no los identifica. En este art´ıculo se muestra un criterio para identificar el bloque Y1 que contiene el grupo m´ as grande de observaciones at´ıpicas.
111
Un criterio para identificar datos at´ıpicos
2.
Resultados b´ asicos del ajuste del modelo de regresi´ on lineal m´ ultiple
Mediante el m´etodo de estimaci´on m´ınimos cuadrados ordinarios (MCO) se obtiene para el modelo dado en (1) los siguientes estimadores: 0 βb = (X X)−1 X 0 Y, Yb = X βˆ = X(X 0 X)−1 X 0 Y = HY, b = Y − Yˆ = Y − HY = (I − H)Y,
(4)
0
SCE = b 0b = [(I − H)Y ] (I − H)Y = Y 0 (I − H)Y. Obs´ervese que la matriz H determina muchos de los resultados de las estimaciones por MCO; por ejemplo, cuando premultiplica al vector de respuestas Y se obtienen los valores predichos de la variable dependiente, por eso en la literatura estad´ıstica en algunos casos la denominan matriz de predicci´ on, y a la matriz I − H la llaman matriz residual, puesto que al antepon´ersele a la variable dependiente Y se obtienen los respectivos residuales.
2.1.
Propiedades de las componentes de la matriz H
En Hoaglin & Welsch (1978) se establece para la matriz H = [hij ] de tama˜ no n × n, las siguientes propiedades: n P P 2 (a) hii = h2ij = h2ii + hij ya que Hes sim´etrica e idempotente. j=1
j6=i
(b) 0 < hii ≤ 1, si i = 1, 2, . . . , n. (c) −0,5 ≤ hij ≤ 0,5, para i 6= j. (d) (1 − hii )(1 − hjj ) − h2ij ≥ 0. (e) hii hjj − h2ij ≥ 0. (f ) Si hii = 1, entonces hij = 0, para todo j 6= i. Si la matriz X de tama˜ no n × r es de rango r, entonces n n n P P P (g) hii = h2ij = r = tr(H), (h)
i=1 n P
i=1
hij =
i=1 j=1 n P
hij = 1,
j=1
donde tr(H) denota la traza de la matriz H. Dado que hij = xi (X 0 X)−1 x0j , entonces hii est´a determinado por la locali-
112
Jos´e A. Jim´enez M.
zaci´ on de xi en el espacio X, es decir, un valor peque˜ no (grande) de hii indica que xi se encuentra cerca (lejos) de la masa de los otros puntos. Adem´as, sugieren que xi es un punto influyente si hii > 2r/n.
3.
C´ alculo de la estad´ıstica Qk
En Jim´enez (2001b) se establece para la estad´ıstica dada en (3), la siguiente expresi´ on: Qk = SCE − SCE ∗ = −2γ 0 ˆ − γ 0 (I − H)γ , (5) donde SCE es obtenida en t´erminos algebraicos como en (4) y SCE ∗ , representa la estimaci´ on v´ıa m´ınimos cuadrados (EM C) de SCE sin el bloque Y1 de observaciones. Adem´ as, muestra que si el inter´es es minimizar la SCE ∗ , esto se logra haciendo: ∂Qk = 0, ∂γ lo cual equivalente a hacer: b + (I − H)b γ = 0,
(6)
donde b es la estimaci´ on v´ıa m´ınimos cuadrados (EM C) de del modelo (1). Al remplazar (6) en (5) se tiene: Qk = γ b0 (I − H)b γ=γ b0 γ b−γ b0 Hb γ.
(7)
Esta nueva expresi´ on de Qk tiene la ventaja de que est´a en t´erminos de la estimaci´ on del γ arbitrario, la cual para los objetivos de este trabajo es m´as atractiva, ya que se podr´ a establecer su distribuci´on de probabilidad correspondiente.
4.
Distribuci´ on de probabilidad de Qk γˆ En Jim´enez (2001a) al asumir la restricci´on γˆ = 1 , se llega a: 0 0 0 γ b −Ik X1 (X2 X2 )−1 X2 Y1 γ b= 1 = , 0 Y2 0 0
(8)
donde Ik es la matriz identidad de tama˜ no k × k, con k igual a la dimensi´on del bloque Y1 y Mij = Xi (X20 X2 )−1 Xj0 .
113
Un criterio para identificar datos at´ıpicos
Si se reemplaza (8) en el primer t´ermino de la expresi´on (7) se obtiene 0 0 0 0 −Ik 0 −Ik X1 (X2 X2 )−1 X2 0 0 γ bγ b=Y Y X2 (X2 X2 )−1 X1 0 0 0 0 Ik −M12 =Y Y. (9) −M21 M21 M12 Por otra parte, si se sustituye (8) en el segundo t´ermino de la expresi´on (7) y se emplean los resultados dados en Jim´enez (2001a), se tiene que: 0 0 H11 H12 − M12 γ b Hb γ =Y Y H21 − M21 H22 + M21 M12 − M22 0 0 H11 H12 0 M12 =Y Y −Y Y. (10) H21 H22 M21 M22 − M21 M12 Finalmente, al sustituir (9) y (10) en la ecuaci´on (7), se obtiene que: 0
Qk =b γ (I − H)b γ 0 0 Ik −M12 H11 H12 − M12 =Y Y −Y Y −M21 M21 M12 H21 − M21 H22 + M21 M12 − M22 0 0 Ik 0 H11 H12 =Y Y −Y Y 0 M22 H21 H22 0
0
0
=Y M Y − Y HY = Y (M − H) Y.
(11)
N´ otese que la matriz (M − H) es sim´etrica; adem´as, es idempotente. Esto se puede verificar de la siguiente manera: (M − H) (M − H) =M 2 − M H − HM + H 2 , pero M 2 = M , ya que: Ik 0 Ik 0 M22 0
0 I = k M22 0
0 I = k M22 M22 0
0 . M22
Esto se tiene, ya que para i, j = 1, 2: 0
0
0
0
0
0
Mi2 M2j = [Xi (X2 X2 )−1 X2 ][X2 (X2 X2 )−1 Xj ] = Xi (X2 X2 )−1 Xj = Mij ; por otra parte, HM = H lo cual se puede verificar como sigue: H11 H12 Ik 0 H11 H12 M22 H11 H12 = = . H21 H22 0 M22 H21 H22 M22 H21 H22
114
Jos´e A. Jim´enez M.
Aqu´ı cabe notar que cuando X = 0
0
X1 es de rango completo, entonces: X2 0
0
0
0
Hi2 M2j = [Xi (X X)−1 X2 ][X2 (X2 X2 )−1 Xj ] = Xi (X X)−1 Xj = Hij , para i, j = 1, 2; adem´ as, como las matrices H y M son sim´etricas se tiene que H = (M H)t = HM . En consecuencia, (M − H) (M − H) = M − H. Para establecer la distribuci´on de Qk , se presentan, sin demostraci´on, los teoremas 1 y 2, mencionados en Searle (1971). Teorema 1. Si Y es un vector aleatorio de tama˜ no n × 1, distribuido N (µ, V ), donde µ es en si mismo un vector entonces: E [Y 0 AY ] = tr(AV ) + µ0 Aµ
y
Var [Y 0 AY ] =2 tr(AV )2 + 4µ0 AV Aµ. 0
Teorema 2. Si Y ∼ N (µ, V ), entonces Y 0 AY ∼ χ2(ν,λ) , con grados de libertad ν = ρ(A) y par´ ametro de no centralidad λ = 21 µ0 Aµ, si y s´olo si AV es idempotente. Puesto que, bajo el supuesto de normalidad en los residuales se tiene que Y ∼ N (Xβ, σ 2 In ).
(12)
Como la expresi´ on dada en (11) es una forma cuadr´atica se establecer´a a continuaci´ on la respectiva distribuci´on asociada. Por el teorema 1, se tiene que " 0 # n h 0 io 0 Y (M − H)Y −1 E = k − r + tr (X X ) (X X ) , 2 2 2 2 σ2 # " 0 n h 0 io 0 Y (M − H)Y −1 Var =2 k − r + tr (X X ) (X X ) , 2 2 2 2 σ2 donde r es el rango de la matriz X definida el modelo (1). h en i Cuando esta 0 0 −1 matriz es de rango completo se tiene que tr (X2 X2 ) (X2 X2 ) = r. Utilizando el teorema 2, tambi´en se concluye que Qk /σ 2 tiene distribuci´on ji-cuadrado central: Qk ∼ χ2(ν) , (13) σ2 h 0 i 0 donde ν = k − r + tr (X2 X2 )−1 (X2 X2 ) . Aqu´ı el teorema 2 es aplicable ya 1 que 2 (M − H)σ 2 In es una matriz idempotente. σ
115
Un criterio para identificar datos at´ıpicos
5.
Metodolog´ıa para establecer datos at´ıpicos Dado que la estad´ıstica Qk se puede obtener de la forma cuadr´atica: 0
Qk = γ b (I − H)b γ,
(14)
al expresarla en t´erminos del vector de respuestas Y , queda como: 0 0 Ik 0 H11 H12 Qk =Y Y −Y Y. 0 M22 H21 H22
(15)
Y1 , el bloque Y1 est´a conformado Y2 por las observaciones at´ıpicas, dicho bloque afectar´a todas las EMC del modelo dado en (1). Por otra parte, si se reescribe la expresi´on (5), se tiene que: Si se considera que en la partici´on Y =
SCE = SCE ∗ + Qk , y dado que SCE ∗ puede expresarse en forma matricial como sigue 0 0 0 0 ∗ SCE = Y Y = Y [In − M ] Y ; 0 In−k − M22
(16)
usando (12), se puede establecer que las expresiones, SCE σ2
y
SCE ∗ , σ2
(17)
tienen distribuci´ on ji-cuadrado central. Luego, si se divide la ecuaci´on (13) por cualquiera de las expresiones dadas en (17), se elimina el t´ermino σ 2 y queda el cociente entre dos formas cuadr´aticas que se distribuyen ji-cuadrado. Por la teor´ıa estad´ıstica se sabe que cuando se realiza el cociente entre dos variables aleatorias independientes con distribuci´on ji-cuadrado y cada una se divide por sus respectivos grados de libertad, se obtiene una nueva variable con distribuci´ on F . Para llevar a cabo el cociente mencionado anteriormente se debe verificar con cu´ al de las distribuciones asociadas a las expresiones dadas en (17) la distribuci´ on de probabilidad expresada en (13) es independiente; para ello, se enuncia sin demostraci´ on el teorema 3, citado en Searle (1971). Teorema 3. Cuando Y ∼ N (µ, V ), las formas cuadr´aticas Y 0 AY y Y 0 BY , est´ an distribuidas independientemente si y s´olo si AV B = 0.
116
Jos´e A. Jim´enez M.
Veamos si las distribuciones asociadas a Qk y SCE son independientes. Si se retoman las ecuaciones dadas en (11) y (4), se tiene por el teorema 3 que Qk y SCE no son independientes, pues, (M − H)(σ 2 In )(In − H) = σ 2 (M − H)(In − H) = σ 2 [M − M H − H + H 2 ] = σ 2 (M − H) 6= 0; en la u ´ltima ecuaci´ on se tuvo en cuenta que H es idempotente y que M H = H. De manera an´ aloga, se verifica si son independientes las distribuciones de probabilidad de Qk y SCE ∗ ; de las ecuaciones (11) y (16) utilizando el teorema 3, se concluye que son independientes, ya que: (M − H)(σ 2 In ) (In − M ) = σ 2 (M − H) (In − M ) = σ 2 M − M 2 − H + HM = 0. En esta u ´ltima expresi´ on se utilizaron los resultados: M H = H y M 2 = M . La media y varianza de la SCE ∗ se obtienen por el teorema 1, como sigue: # 0 n h 0 io 0 Y (In − M ) Y −1 = n − k − tr (X X ) (X X ) , E 2 2 2 2 σ2 " 0 # n h 0 io 0 Y (In − M ) Y −1 V ar =2 n − k − tr (X X ) (X X ) . 2 2 2 2 σ2 "
Como la media y la varianza de la distribuci´on χ2η son η y 2η respectivamente, se 0 deduce que Y (In − M ) Y /σ 2 tiene distribuci´on ji-cuadrado central. Se llega 1 a la misma conclusi´ on, ya que 2 (In − M ) σ 2 In es idempotente, utilizando el σ teorema 2. As´ı pues, 0
Y (In − M ) Y ∼ χ2η , σ2
(18)
0 0 con η = n−k−tr (X2 X2 )−1 (X2 X2 ) . Cuando la matriz X es de rango completo 0 0 se tiene que tr (X2 X2 )−1 (X2 X2 ) = r. Como las distribuciones de probabilidad asociadas a las expresiones (15) y (16) son independientes, al hacer el cociente entre las relaciones (13) y (18),
117
Un criterio para identificar datos at´ıpicos
dividiendo cada una por sus correspondientes grados de libertad, se llega a: Qk n−r−kγ b(I − H)b γ kσ 2 = , ∗ SCE ∗ k SCE (n − r − k)σ 2 n−r−k Qk ∼ F(k,n−r−k) . k SCE ∗ Estos resultados se pueden resumir en los siguientes teoremas. Teorema 4. Si en un modelo de regresi´on lineal m´ ultiple particionado como: Y1 X1 = β+ 1 , Y2 X2 2 se elimina el bloque Y1 de dimensi´on k, entonces el cambio que se presenta en la SCE se calcula mediante la expresi´on: ∆(Y1 ) =
b0 [In − H] γ b 1 γ , 2 k S(Y1 )
(19)
SCE ∗ es la estimaci´on usual de σ 2 , despu´es de eliminar n−k−r 0 0 γ b las observaciones del bloque Y1 , y γ b = 1 , con γ b1 = −Y1 +X1 (X2 X2 )−1 X2 Y2 . 0 2 donde S(Y =σ b2 = 1)
Teorema 5. En un modelo de regresi´on lineal m´ ultiple Y = Xβ + , bajo el supuesto de que ∼ N (0, σ 2 In ), se tiene que: ∆(Y1 ) ∼ F(k,n−r−k) ,
con
k = dimensi´on del bloque Y1 , r = rango de la matriz X.
En este caso, se clasifica como at´ıpico al bloque Y1 de observaciones, si con un nivel de significancia α se satisface que: ∆(Y1 ) > F(k,n−r−k,α/2) .
6.
(20)
Ejemplo
En la Tabla 1, se considera el conjunto de 21 observaciones (x, y), dado por Mickey, Dunn & Clark (1967). Para este conjunto de datos, se presentan los siguientes resultados:
118
Jos´e A. Jim´enez M.
Tabla 1: Datos de Mickey, Dunn, and Clark (1967) Obs.
x
y
1 2 3 4 5 6 7
15 26 10 9 15 20 18
95 71 83 91 102 87 93
Obs.
x
y
8 9 10 11 12 13 14
11 8 20 7 9 10 11
100 104 94 113 96 83 84
Obs.
x
y
15 16 17 18 19 20 21
11 10 12 42 17 11 10
102 100 105 57 121 86 100
1. La estimaci´ on del modelo de regresi´on lineal, con las 21 observaciones. 2. Los elementos de la diagonal de la matriz H, las estimaciones de los γi y al eliminar el i-´esimo dato se establecen la estad´ıstica Q1 , la distancia de Cook y la estad´ıstica ∆(i) con su p-valor correspondiente. 3. La estimaci´ on del modelo de regresi´on lineal, despu´es de eliminar la observaci´ on influyente determinada mediante distancia de Cook. 4. La estimaci´ on del modelo de regresi´on lineal, sin la observaci´on que se considera influyente por la estad´ıstica ∆(i) . 1. An´ alisis de varianza para el conjunto completo de datos: Fuente de variaci´ on
Grados libertad
Suma de cuadrados
Cuadrados Medios
F
Valor cr´ıtico de F
Regresi´ on Residuos Total
1 19 20
1604,0809 2308,5858 3912,6667
1604,0809 121,5045
13,2018
0,00177
Coeficiente de determinaci´on R2 = 0,409971261:
Intercepto Variable X
Coeficientes
Error t´ıpico
Estad´ıstico t
109,8738 -1,1270
5,0678 0,3102
21,6808 -3,6334
119
Un criterio para identificar datos at´ıpicos
2. Compendio de estad´ısticas: Obs. Elim.
hii
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
0,0479 0,1545 0,0628 0,0705 0,0479 0,0726 0,0580 0,0567 0,0799 0,0726 0,0908 0,0705 0,0628 0,0567 0,0567 0,0628 0,0521 0,6516 0,0531 0,0567 0,0628
γbi
Qk
D∗i
(k=1)
(100∗Di )
-2,1332 11,3214 16,6498 9,3936 -9,4856 0,3602 -3,6220 -2,6746 -3,4148 -7,1879 -12,1145 4,0141 16,6498 14,2866 -4,7948 -1,4896 -9,1255 15,9026 -31,9816 12,1664 -1,4896
4,333 108,370 259,803 82,015 85,664 0,120 12,358 6,748 10,729 47,914 133,443 14,976 259,803 192,540 21,687 2,080 78,936 88,105 968,562 139,634 2,080
0,09 8,15 7,17 2,56 1,77 0,00 0,31 0,17 0,38 1,54 5,48 0,47 7,17 4,76 0,54 0,06 1,79 67,81 22,33 3,45 0,06
Valor hii > 4/21 |γi | ≥ γj para todoj Inusual
el m´ as grande
Di > 0, 5
∆(i)
p−valor
0,0338 0,8866 2,2826 0,6630 0,6937 0,0009 0,0969 0,0528 0,0840 0,3815 1,1043 0,1175 2,2826 1,6378 0,1707 0,0162 0,6373 0,7142 13,0103 1,1588 0,0162
0,8561 0,3589 0,1482 0,4261 0,4158 0,9759 0,7592 0,8209 0,7752 0,5445 0,3072 0,7357 0,1482 0,2169 0,6844 0,9000 0,4351 0,4091 0,0020 0,2959 0,9000 p