Detección de valores atípicos multivariantes con Mahalanobis robusto

Question

Estoy buscando algunos documentos y ejemplos relacionados con la detección multivariante de valores atípicos con distancia mahalanobis robusta (estimación de covarianza mínima). Tengo 6 variables y quiero trazarlas para mostrar también los valores atípicos. ¿Tienen alguna fuente?

Aquí están los códigos, pero creo que algo va mal. Porque tengo más de 2 millones de casos que ha tomado sólo n = 500.

> CovMcd(new)

Call:
CovMcd(x = new)
-> Method:  Fast MCD(alpha=0.5 ==> h=1342195); nsamp = 500; (n,k)mini = (300,5) 

Robust Estimate of Location: 
logdgr  logtr       lph       lpm       lpr  
   4.391     2.956    -2.722    -4.802    -4.802  

Robust Estimate of Covariance: 
          logdgr  logtr  lph     lpm     lpr   
logdeger  1.0183    0.8981    0.6427  0.7112  0.7113
loghacim  0.8981    1.0173    0.9613  0.9539  0.9541
lph       0.6427    0.9613    1.6921  1.3770  1.3772
lpd       0.7112    0.9539    1.3770  1.3085  1.3087
lpr       0.7113    0.9541    1.3772  1.3087  1.3089
> summary(mcd)

Call:
CovMcd(x = data)

Robust Estimate of Location: 
[1]  3.5  2.0

Robust Estimate of Covariance: 
      [,1]  [,2]
[1,]  3.5   0.8 
[2,]  0.8   0.8 

Eigenvalues of covariance matrix: 
[1]  3.7192  0.5808

Robust Distances: 
[1]  2.0833  0.8333  2.0833  2.0833  0.8333  2.0833
> mest <- CovMest(new)
> show(mcd)

Call:
CovMcd(x = data)
-> Method:  Fast MCD(alpha=0.5 ==> h=4); nsamp = 500; (n,k)mini = (300,5) 

Robust Estimate of Location: 
[1]  3.5  2.0

Robust Estimate of Covariance: 
      [,1]  [,2]
[1,]  3.5   0.8 
[2,]  0.8   0.8

Answer 1

Lo que se intenta hacer parte de la siguiente idea básica:

1. Calcular las distancias robustas (utilizando la estimación del Determinante de Covarianza Mínimo, MCD) de cada observación $$rd(x_{i}),\ i = 1 \ldots n,\ x_i \in \mathbb{R}^p, $$ donde $n$ y $p$ son el número de observaciones (filas) y de variables (columnas) respectivamente.

2. Compara cada uno $rd(x_{i})$ a $\sqrt{\chi^2_{p,.975}}$ . Declarar $x_i$ un valor atípico si $$rd(x_{i}) > \sqrt{\chi^2_{p,.975}}$$

3. Las distancias robustas vienen dadas por: $$rd(x_i) = \sqrt{(x_i - \mu_{mcd} )'S_{mcd}^{-1} (x_i - \mu_{mcd} )}$$ donde $ \mu_{mcd}$ y $S_{mcd}$ son las estimaciones MCD robustas de la localización (vector medio) y la dispersión (matriz de covarianza) respectivamente.

Por lo tanto, lo que tienes que hacer es:

A. Obtenga la estimación MCD robusta de la ubicación ( $ \mu_{mcd}$ ) y la dispersión ( $S_{mcd}$ ) para sus datos. Está en el camino con el CovMCD(x) en R. Tenga en cuenta que esta función implementa por defecto el algoritmo FastMCD, que toma repetidamente submuestras de sus datos (cada una de ellas de tamaño denotado por $h$ ) para hacer estimaciones. El número de submuestras que se toman por defecto en esta función es de 500, lo que explica la razón por la que ves 500 ahí. La función no está tomando 500 observaciones de sus datos. En su lugar, toma submuestras de tamaño $n_{subsample} = h$ repetidamente para hacer estimaciones. Se tomarán 500 de estas submuestras para hacer estimaciones. Compruebe Peter J. Rousseeuw y Katrien Van Driessen (1999) A Fast Algorithm for the Minimum Covariance Determinant Estimator, Technometrics, 41:2, 212-223 para más detalles.

B. Obtenga las distancias robustas para cada fila u observación utilizando la fórmula del punto 3 anterior.

C. Comparar cada distancia $rd(x_i)$ a $\sqrt{\chi^2_{p,.975}}$ y declarar $x_i$ un valor atípico si $rd(x_{i}) > \sqrt{\chi^2_{p,.975}}$ .

Y a continuación se proporciona un ejemplo de código R:

require(rrcov)
data(hbk)
mcd <- rrcov::CovMcd(hbk[,1:3]) # use only first three columns  
# get mcd estimate of location
mean_mcd <- mcd@raw.center
# get mcd estimate scatter
cov_mcd <- mcd@raw.cov

# get inverse of scatter
cov_mcd_inv <- solve(cov_mcd)

# compute distances

# compute the robust distance
robust_dist <- apply(hbk[,1:3], 1, function(x){
  x <- (x - mean_mcd)
  dist <- sqrt((t(x)  %*% cov_mcd_inv %*% x))
  return(dist)
})

# set cutoff using chi square distribution
threshold <- sqrt(qchisq(p = 0.975, df = ncol(hbk[,1:3]))) # df = no of columns

# find outliers
outliers <-  which(robust_dist >= threshold) # gives the row numbers of outliers