Distanza di Mahalanobis e outliers multivariati

1 Metrica di MahalaNobis e ricerca di outlier multivariati

Si introduce la distanza di Mahalanobis, che ci dà anche uno strumento per definire degli outliers multivariati

Preferisco in questo caso il termine anglosassone outlier ossia che sta lontano invece della pessima traduzione italiana Dato anomalo, invalsa qualche decina di anni fa, che identifica gli outlier con dei dati sbagliati . Semplicemente, per il mio modo di vedere la statistica, più che anomali si tratta di dati che stanno lontano, forse sono lontani, o forse sono un po’ differenti, forse sono indizio di un’eterogeneità insita nei dati.

In effetti nell’epoca attuale di grande disponibilità di dati, se non si ha a che fare con indagini pianificate, molto spesso dal punto di vista esplorativo si cercano tecniche che a dare una buona descrizione della maggior parte dei dati (non necessariamente tutti) ## Inizializzazione

Vengono caricati i packages necessari per realizzare questo documento

1.1 Variabili osservate.

matrice di grafici per le sole variabili Altezza, Span, Torace rilevate su 1427 righe

source

round(      cor(dati),3   )

output

        Altezza  Span Torace
Altezza   1.000 0.892  0.587
Span      0.892 1.000  0.578
Torace    0.587 0.578  1.000

source

MLA.explor.pairs(dati)

plot

1.2 Definizione di distanza Mahalanobis.

E’ una forma di distanza standardizzata, in cui si tiene conto non solo della diversa dispersione delle variabili, ma anche della loro correlazione.

1.3 esempio con tre variabili:

output

        Altezza  Span Torace
Altezza   1.000 0.892  0.587
Span      0.892 1.000  0.578
Torace    0.587 0.578  1.000

You must enable Javascript to view this page properly.

1.4 La metrica euclidea è sempre adeguata?

La presenza di forti correlazioni fra le variabili può rendere inappropriato l’uso della distanza euclidea.

Osserviamo la nuvola di punti riportati nella figura.

Con una sfera rossa è indicato il centroide \(C\) dei punti, di coordinate \(C=\{0, 0, 0\}\) perchè si sono standardizzati i dati.

Poi vengono evidenziati due punti, ossia quelli di indice 894, 178, che approssimativamente sono equidistanti dal centroide. Le loro distanze euclidee dal centroide sono infatti, rispettivamente,

\[d(P_{894},C)=\sum_{j=1}^3 (x_{ij}-C_j)^2 = 2.517\] e \[d(P_{178},C)=\sum_{j=1}^3 (x_{ij}-C_j)^2 =2.44\].

Eppure palesemente il punto \(P_{894}\) sembra più lontano dal centro, rispetto al punto \(P_{178}\), che invece sembra più interno rispetto alla massa dei dati. Il punto \(P_{894}\) invece sembra più esterno, sta in una zona più in cui non vi sono altri punti prossimi.

Questo accade perchè la distanza eucclidea non tiene conto delle correlazioni fra le variabili

You must enable Javascript to view this page properly.

D’altra parte, se ora rappresentiamo una sfera di raggio 2.45 centrata nel centroide, sulla sua superficie abbiamo i punti che distano 2.45 dal centro, secondo una metrica euclidea, vediamo che in effetti questa sfera ha una forma che non somiglia alla nostra nuvola di punti 3D.

1.5 Ellissi di concentrazione

Consideriamo la forma della nuvola di punti e costruiamo un ellissoide di concentrazione, ossia un ellissoide centarto sul centroide dei dati e che ha come assi gli assi principali della nuvola dei punti.

Sarebbe l’ellissoide di equiprobabilità di una normale multivariata con gli stessi primi due momenti multivariati dei nostri dati (vettore medie e matrice di varianze e covarianze).

(onestamente non è necessario ricorrrere alla normale multivariata per giustificare il nostro ragionamento: basta dire che prendiamo un ellisse centrato nel centroide dei dati e con gli stessi assi principali della nostra nuvola di punti)

output

        Altezza  Span Torace
Altezza   1.000 0.892  0.587
Span      0.892 1.000  0.578
Torace    0.587 0.578  1.000

output

glX 
  4 

You must enable Javascript to view this page properly.

Nella figura è rappresentato l’ellissoide ad un livello di probabilità di 0.95. E’ importante ruotare la figura interattivamente per rendersi conto dell’effetiva distanza del punto \(P_{894}\) dal centro

Rispetto a questa rappresentazione \(P_{178}\) è senz’altro interno a tale ellissoide, mentre \(P_{894}\) è palesemente esterno.

E' come dire che in una distribuzione normale multivariata questo punto più esterno ha una densità molto bassa

1.6 Definizione della distanza di Mahalanobis

Praticamente utilizziamo la formula: \[ d^2(P_i,C)=[\mathbf{x}_i-C]^T\Sigma^{-1}[\mathbf{x}_i-C] \] che sostanzialmente costituisce una sorta di standardizzazione della differenza fra \(\mathbf{x}_i\) e \(C\), che tiene conto non solo degli elementi diagonali di \(\Sigma\), ma anche delle correlazioni.

In effetti l’equazione del generico ellisse centrato nel centroide e con assi principali coincidenti con gli assi principali dei dati, è data da:

\[ [\mathbf{x}_i-C]^T\Sigma^{-1}[\mathbf{x}_i-C]=k^2; \ k\ge 0 \] Al variare di \(k\) variano le distanze dei punti. Fissato \(k\) l’equazione descrive il luogo dei punti equidistanti dal centroidi (equidistanti in questo nuovo senso, secondo la distanza di Mahalanobis)

Vediamo anche che è una quantità proporzionale all’esponente della normale multivariata di parametri \(C, \Sigma\).

source

x1=rid[ind[1],]
x2=rid[ind[2],]
d1=sqrt(crossprod((x1-med),solve(sigma))%*%(x1-med))
d2=sqrt(crossprod((x2-med),solve(sigma))%*%(x2-med))
print(c(d1,d2))

output

[1] 5.249053 1.774988

Adesso vediamo che la diversa distanza dal centro dei dati è ben rispecchiata, e che con la distanza di Mahalanobis il punto \(P_{894}\) risulta quasi tre volte più lontano dal centroide rispetto a \(P_{178}\)!

output

glX 
  6 

output

Importance of components:
                          Comp.1    Comp.2     Comp.3
Standard deviation     1.5430776 0.7134693 0.32828463
Proportion of Variance 0.7942527 0.1697985 0.03594879
Cumulative Proportion  0.7942527 0.9640512 1.00000000

You must enable Javascript to view this page properly.

Ancora ci possiamo rifare alla rappresentazione degli ellissoidi centrati nel centroide, e adesso rappresentiamo proprio quelli che passano dai due punti, così da apprezzare ancora di più il concetto.

Perchè \(P_{178}\) è più vicino a C secondo la metrica di Mahalanobis rispetto a \(P_{894}\) ? Perchè sebbene le distanze euclidee siano simili, \(P_{178}\) è su un ellisse di equiprobabilità molto più interno rispetto a \(P_{894}\)

1.7 Luogo dei punti equidistanti

Ciascun ellissoide è il luogo dei punti equidistanti dal centroide secondo la metrica di Mahalanobis.

1.8 Outliers multivariati

Possiamo anche utilizzare questa distanza per definire in qualche modo il concetto di outlier multivariato.

E’ logico secondo quanto visto finora pensare a \(P_{894}\) come a un outlier in senso multivariato.

Infatti se guardiamo le coordinate (standardizzate) univariate di \(P_{894}\), ossia 2.28, 0.13, -0.85, in effetti rispetto al primo asse forse il punto è un po’ lontano, ma poco più di due unità standardizzate.

In sostanza il punto \(P_{894}\) da un punto di vista univariato, ossia marginalmente, non verrebbe considerato un outlier rispetto a nessuno dei tre assi come si vede dai grafici, come si vede anche dagli istogrammi delle distribuzioni delle tre variabili, ove ho rappresentato con un punto rosso la coordinata del punto \(P_{894}\)

source

for (i in 1:3){
  hist(rid[,i],main=paste("Posizione del punto rispetto alla distribuzione della var.",i))
points(rid[ind[1],i],0,col=2,cex=1.5,pch=19)
}

plot

plot

plot

1.9 Ordinamento dei punti secondo la distanza di Mahalanobis

Se invece rappresentiamo l’istogramma di tutte le distanze di mahalanobis, e rappresentiamo ancora con un punto rosso il punto \(P_{894}\), otteniamo una situazione ben diversa.

source

d.tot   =sqrt(rowSums(crossprod(t(rid),solve(sigma))*rid))
hist(d.tot)
d1=d.tot[ind[1]]

points(d1,0,col=2,cex=1.5,pch=19)

plot

source

print(d1)

output

[1] 5.249053

Adesso il punto si trova all’estremità dell’istogramma delle distanze. Questo succede perchè le tre coordinate del punto \(P_{894}\), ossia 2.28, 0.13, -0.85, non rispecchiano la correlazione generale fra le tre variabili nell’intero insieme di dati.

Possiamo addirittura ordinare i punti secondo questa distanza di Mahalanobis ed evidenziare i più lontani

source

ind.mahala=order(d.tot, decreasing = TRUE)
print(cbind(rid[ind.mahala[1:10],],d.tot[ind.mahala[1:10]]))

output

           Altezza        Span     Torace         
 [1,]  0.005137188  2.62796607 -0.3368174 5.945929
 [2,]  2.283619083  0.12533213 -0.8489520 5.249053
 [3,] -2.174280276  0.03595234 -1.3610866 4.940315
 [4,]  0.302330479 -1.75164333  0.3033509 4.557688
 [5,] -1.381764835  0.57223105 -1.6171539 4.518618
 [6,]  0.500459339  2.18106715  2.6079568 4.513095
 [7,]  0.005137188  0.57223105  3.5041924 4.371025
 [8,] -1.480829265  0.30409170 -2.0012549 4.349332
 [9,]  0.599523770  1.46602888  3.5041924 4.217928
[10,]  0.896717060 -1.03660506  0.4313846 4.189782

1.10 Un test di multinormalità

Andiamo ancora avanti! Se adesso prendiamo i quadrati delle distanze di Mahalanobis, si può dimostrare che in una normale multivariata si distribuiscono come una \(\chi^2\) con \(p\) gradi di libertà:

\[ Y \sim \mathcal{N}(\mu,\Sigma) \ \ \Rightarrow \ \ (Y-\mu)^T\Sigma^{-1}(Y-\mu)\sim \chi^2_p \]

Quindi dato un campione casuale semplice di \(n\) unità realmente provenienti da una normale bivariata, i quadrati delle distanze di mahalanobis dal centroide dovrebbero costituire un campione di ampiezza \(n\) proveniente da una distribuzione \(\chi^2\) (per valori grandi n, per i quali è poco influente l’effetto di utilizzare le stime di \(\mu,\Sigma\), invece dei veri valori incogniti )

source

isort =order(d.tot)
p     =pchisq((d.tot[isort])^2,3)  
plot(p, ((1:n)-0.5)/n)
abline(a=0,b=1,col=3)
d1    =d.tot[ind[1]]
print(d.tot[n])

output

[1] 1.466693

source

points(p[n],(n-0.5)/n,col=2,cex=1.5,pch=19)

plot

source

rgl.close()

1.11 Collegamento con l’analisi delle componenti principali

Con un po’ di pratica di algebra delle matrici non è difficile spiegare la distanza di Mahalanobis sfruttando le componenti principali. notebook introduttivo su ACP

Che un collegamento ci sia è ovvio per il modo in cui abbiamo costruito la distanza, considerando equidistanti i punti che si trovano sulla superficie (o iper-superficie) di un ellisoide orientato come gli assi principali definiti dalla nuvola dei punti dei dati.

Sappiamo che \(\Sigma\), definita positiva, si può scomporre in: \[\Sigma=\Gamma \Lambda \Gamma^T\] e quindi \[\Sigma^{-1}=\Gamma\Lambda^{-1} \Gamma^T \] perchè l’inversa ha gli stessi autovettori e autovalori inversi.

Sostituendo questa espressione nella distanza di Mahalanobis: \[ d^2(P_i,C)=[\mathbf{x}_i-C]^T\Sigma^{-1}[\mathbf{x}_i-C]= [\mathbf{x}_i-C]^T\Gamma \Lambda^{-1} \Gamma^T[\mathbf{x}_i-C] \] ma \([\mathbf{x}_i-C]^T\Gamma= \mathbf{y}_i^T\) è il vettore delle coordinate dell\(i-\)esimo punto nel nuovo sistema di riferimento e quindi: \[ d^2(P_i,C)= \mathbf{y}_i^T \Lambda^{-1} \mathbf{y}_i=\mathbf{y}_i^T \Lambda^{-1/2} \mathbf{y}_i \Lambda^{-1/2} \] Gli elementi di \(\Lambda^{-1/2}\) sono gli iversi delle deviazioni standard degli assi principali! e quindi la distanza di Mahalanobis corrisponde a una distanza euclidea standardizzata nel sistema di riferimento delle componenti principali