Un programa en STATA para identificar y clasificar outliers multivariados

En el análisis de regresión, pocas observaciones con comportamiento diferente a la mayoría (outliers) pueden distorsionar completamente las estimaciones de los parámetros, e inducir a inferencias erróneas sobre las relaciones entre las variables. La identificación de outliers en el análisis multivariante no es sencillo por dos motivos:
(i) Observaciones con valores extremos en una o dos dimensiones pueden no serlo en más dimensiones; así, desacertadamente se podría eliminar observaciones que sí pertenecen a la distribución de los datos.
(ii) Las medidas tradicionales para identificar outliers multivariados en el análisis de regresión como los residuos estandarizados, medidas de apalancamiento o distancias de Cook, etc., sufren los llamados "masking y swamping effects": estas medidas utilizadas para identificar outliers, se encuentran por su construcción ya distorsionadas por los outliers; así las mismas erróneamente pueden sub categorizar (masking) o sobre categorizar (swamping) observaciones con el status de outliers.
Este post presenta un programa en STATA del tipo do-file para identificar y categorizar outliers multivariados para datos de corte transversal, implementando la herramienta sugerida por Rousseeauw y vanZomeren (1990), la misma que no sufre de los mencionados "masking y swamping effects". Espero les sea útil.

a) Tipos de outliers en datos de corte transversal

La literatura sobre Estadística Robusta sugiere la existencia de tres tipos de outliers, véase la siguiente figura: (i) los llamados "vertical outliers", son las observaciones con valores extremos en la dimensión-y (la variable dependiente), pero que tienen un comportamiento similar al de la muestra en la dimensión-x; las mismas pueden sesgar el valor del intercepto y la pendiente. (ii) Los "good leverage points" (puntos de influencia buenos), son observaciones con valores extremos en las dimensiones y y x (son outliers en la variable dependiente y en los regresores), pero que cumplen la característica que no se encuentran muy lejos de la línea de regresión, justamente por ello, las mismas no afectan a los coeficientes, pero si la inferencia (inducen a rechazar la hipótesis nula de no significatividad, Dehon et al., 2009). (iii) Los "bad leverage points" (puntos de influencia malos) son observaciones que se encuentran alejados de la muestra en la dimensión-x, y a su vez lejos de la línea de regresión; los mismos sesgan fácilmente ambos coeficientes, el intercepto y la pendiente.

Fuente: Verardi y Croux (2009)

b) El esquema de Rousseauw y vanZomeren (1990) y su código en STATA

El mismo consiste en graficar en el eje de las ordenadas los residuos estandarizados y en el eje de las abscisas las distancias de Mahalanobis. Para mayores detalles y evidencia sobre su "robustez", véase Rousseeauw y vanZomeren (1990), Maronna et. al. (2006).

Los residuos estandarizados son la medida utilizada para identificar valores extremos en la dimensión-y, éstas consisten en el cociente residuos sobre su desviación estándar, \({r_i}/\sigma{(r)_i}\). Para que ésta sea resistente a outliers (robusto), en el numerador se utilizan los residuos de la regresión S, y en el denominador, como medida de dispersión se utiliza la desviación absoluta mediana normalizada: \(MAD(r)={Med\left| {{r}_{i}}-Med(r)\right|}/{0.6745}\;\). Valores de estos residuos mayores a 2,25 en valor absoluto requieren atención ya que pueden ser o "good leverage points" (si a su vez estas observaciones son valores extremos en x), o "vertical outliers" (si no son outliers en x). 

Para medir valores extremos en la dimensión-x, para cada observación se calcula la distancia de Mahalanobis, \(M{{D}_{i}}=\sqrt{({{x}_{i}}-\mu){{\Sigma }^{-1}}{{({{x}_{i}}-\mu )}^{T}}}\), donde \(\mu \) representa el vector de medias y \(\Sigma \) la matriz de covarianzas; las \(M{{D}_{i}}\) pueden entenderse como la distancia estandarizada de cada observación al centro de los datos. Al tener una distribución chi-2, observaciones con valores mayor a \(M{{D}_{i}}>\sqrt{\chi _{NumdeVariables,0.975}^{2}}\) pueden definirse como valores extremos en la dimensión-x, y éstas pueden ser o "good leverage points" (si a su vez estas observaciones son valores extremos en y), o "bad leverage points" (si no son outliers en x). Asimismo, para que las \(M{{D}_{i}}\) sean robustas, \(\mu\) y \(\Sigma \) se calculan a través del estimador S-multivariado; éste último detalle define la contribución del presente programa.

El programa en STATA para obtener este gráfico es el siguiente. El programa hace uso de los comandos smultiv, sregress, msregress de Verardi y Croux (2010) es muy importante que el lector instale los mismos en su computadora. Puede hacerlo utilizando los comandos findit o ssc install.

*******************

program define outid
syntax varlist , [dummies(varlist)]
local dv: word 1 of `varlist'
local expl: list varlist - dv
local ndum: word count `dummies'
local nvar: word count `varlist'
local p=`ndum'+`nvar'
local b=sqrt(invchi2(`p'),0.975)
capture drop outS rdS id
capture qui smultiv `expl', gen(outS rdS) dummies(`dummies')
label var rdS "Robust_distance_S"
gen id=_n
if `ndum'==0 {
capture drop S_outlier S_stdres
capture qui mmregress `dv' `expl', outlier
capture drop stdres
rename S_stdres stdres
}
else {
capture drop MS_outlier  MS_stdres
capture qui mmregress `dv' `expl', outlier dummies(`dummies')
capture drop stdres
rename MS_stdres stdres
}
label var stdres "Robust standardized residuals"
twoway (scatter stdres rdS if abs(stdres)<4&rdS<sqrt(2)*`b') (scatter stdres rdS if abs(stdres)>=4|rdS>=2*`b', mlabel(id) msymbol(circle_hollow)), xline(`b') yline(2.25) yline(-2.25) legend(off)
capture drop vo glp blp
gen vo=(abs(stdres)>2 & rdS<`b')
gen glp=(abs(stdres)<2 & rdS>`b')
gen blp=(abs(stdres)>2 & rdS>`b')
label var vo "Vertical outliers"
label var glp "Good leverage points"
label var blp "Bad leverage points"
edit id vo glp blp stdres rdS `varlist' `dummies' if vo==1 | glp==1 | blp==1
end
*******************

Para implementar estos códigos en STATA, el usuario debe copiar los mismos en un do-file, correr el programa y luego aplicar el mismo a su modelo de regresión de acuerdo a la siguiente sintaxis:

outid   Variable_dependiente    Variables_explicativas _continuas, dummies(Variables_explicativas_categóricas)

Un ejemplo de aplicación se presenta a continuación. En un do-file se puede escribir:

clear
set obs 300
set seed 1010
drawnorm x1-x5 e
gen i= _n
gen d1=(x4>0.7)
gen d2=(x5<-0.9)
gen y=x1+x2+x3+d1+d2+e
replace x1=invnorm(uniform())+10 in 1/20
replace y=invnorm(uniform())+10 in 15/30
scatter y x1, mlabel(i)

En este esquema las observaciones 1 al 14 son outliers del tipo bad leverage points, las observaciones 15 a 20 son good leverage points, mientras que las observaciones 20 a 30 son vertical outliers. En este escenario, para identificar y clasificar los outliers se puede implementar el programa descrito mediante la siguiente sintaxis en STATA:

outid y x1 x2 x3, dummies(d1 d2)

Los resultados que se obtienen al implementarse el programa outid son los siguientes: el gráfico de Rousseeauw y vanZomeren (1990) y la representación en la base de datos de las observaciones identificadas como outliers.

El siguiente diagrama es útil para fines de interpretación.

Una vez identificados los outliers se debería analizar cada caso, si corresponden a la muestra, o si son typos, o tal vez información errónea del proceso de levantamiento de los datos. Asumiendo que todos corresponden a la muestra pero no representan a la misma, una estrategia racional es correr las regresiones con y sin outliers, y comparar los coeficientes para ver el nivel de importancia de estas observaciones sobre las estimaciones (se podría usar un test de Haussman para este propósito).  Es decir:

Nótese que la primera regresión es con la muestra contaminada por los outliers, mientras que la segunda es con la muestra sin los outliers identificados con el programa outid.

Como se expresó anteriormente, la contribución de este programa radica en el uso del estimador S-multivariado para calcular el centro y dispersión de los datos de manera robusta. En STATA, los comandos de Verardi y Croux (2010) permiten obtener la representación de Rousseeauw y vanZomeren (1990), utilizando la opción graph, pero utilizando el estimador Determinante de Covarianza Mínima (Minimum Covariance Determinant, MCD) para identificar los outliers en x. Sin embargo, como demuestran Verardi y McCathie (2013), el MCD es inestable, presenta baja eficiencia, es computacionalmente pesado, y es superado en propiedades por el estimador S-multivariado, el mismo que se implementa en outid. Para mayores detalles y evidencia sobre las propiedades en muestras finitas del presente programa véase Ugarte y Morales (2018).

Referencias bibliográficas
- Dehon, C., Gassner, M., y Verardi, V. (2009). Beware of good leverage points, Oxford Bulletin of economics and statistics, 71,3.
-Maronna, R., Martin, D., y Yohai, V. 2006. Robust Statistics: Theory and Methods. Book Series: Wiley Series in Probability and Statistics.

- Rousseeuw, P., y van Zomeren, B. (1990). Unmasking multivariate outliers and leverage points. Journal of the American Statistical Association, 85,141.
- Ugarte y Morales (2018). Técnicas robustas y no robustas para identificar outliers multivariados en el análisis de regresión, mimeo.
- Verardi, V., y Croux, C. (2009). Robust Regresion in Stata. Stata Journal, 9,3.

- Verardi, V., y Dehon, C. (2010). Multivariate outlier detection in Stata. Stata Journal, 10,2.

- Verardi, V., y McCathie, A. (2012). The S-estimator of multivariate location and scatter in Stata. Stata Journal, 12,2.

°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°