Rimuovere valori anomali, dal calcolo coefficiente di correlazione

Question

Supponiamo di avere due numerico vettori x e y. Il coefficiente di correlazione di Pearson tra x e y è dato da

  

cor (x, y)

Come posso considerare automaticamente solo un sottoinsieme di x e y nel calcolo (diciamo il 90%) da massimizzare il coefficiente di correlazione?

Answer 1

Se davvero voglio fare questo (togliere i più grandi residui (assoluti)), allora possiamo utilizzare il modello lineare per stimare la soluzione meno quadrati e residui associati e quindi selezionare la centrale n% dei dati. Ecco un esempio:

In primo luogo, generare alcuni dati fittizi:

require(MASS) ## for mvrnorm()
set.seed(1)
dat <- mvrnorm(1000, mu = c(4,5), Sigma = matrix(c(1,0.8,1,0.8), ncol = 2))
dat <- data.frame(dat)
names(dat) <- c("X","Y")
plot(dat)

Successivamente, montare il modello lineare ed estrarre i residui:

res <- resid(mod <- lm(Y ~ X, data = dat))

La funzione quantile() ci può dare i quantili richieste dei residui. È suggerito conservando il 90% dei dati, in modo da vogliamo che il superiore e inferiore 0,05 quantili:

res.qt <- quantile(res, probs = c(0.05,0.95))

Seleziona tali osservazioni con residui in mezzo il 90% dei dati:

want <- which(res >= res.qt[1] & res <= res.qt[2])

Si può quindi visualizzare questo, con i punti rossi sono quelli tratterremo:

plot(dat, type = "n")
points(dat[-want,], col = "black", pch = 21, bg = "black", cex = 0.8)
points(dat[want,], col = "red", pch = 21, bg = "red", cex = 0.8)
abline(mod, col = "blue", lwd = 2)

Le correlazioni per la completa dei dati e il sottoinsieme selezionato sono i seguenti:

> cor(dat)
          X         Y
X 1.0000000 0.8935235
Y 0.8935235 1.0000000
> cor(dat[want,])
          X         Y
X 1.0000000 0.9272109
Y 0.9272109 1.0000000
> cor(dat[-want,])
         X        Y
X 1.000000 0.739972
Y 0.739972 1.000000

Si noti che qui potremmo buttare fuori perfettamente buoni dati, perché abbiamo appena scelto il 5%, con grandi residui positivi e 5% con il più grande negativo. Un'alternativa è quella di selezionare il 90% con piccole residui assoluti:

ares <- abs(res)
absres.qt <- quantile(ares, prob = c(.9))
abswant <- which(ares <= absres.qt)
## plot - virtually the same, but not quite
plot(dat, type = "n")
points(dat[-abswant,], col = "black", pch = 21, bg = "black", cex = 0.8)
points(dat[abswant,], col = "red", pch = 21, bg = "red", cex = 0.8)
abline(mod, col = "blue", lwd = 2)

Con questo leggermente diverso sottoinsieme, la correlazione è leggermente inferiore:

> cor(dat[abswant,])
          X         Y
X 1.0000000 0.9272032
Y 0.9272032 1.0000000

Un altro punto è che anche in questo caso stiamo gettando fuori buoni dati. Si potrebbe desiderare di guardare distanza di Cook come una misura della forza dei valori anomali, e scartare solo quei valori di sopra di una certa soglia di distanza di Cook. Wikipedia ha informazioni sulla distanza di Cook e soglie proposte. La funzione cooks.distance() può essere utilizzata per recuperare i valori da mod:

> head(cooks.distance(mod))
           1            2            3            4            5            6 
7.738789e-04 6.056810e-04 6.375505e-04 4.338566e-04 1.163721e-05 1.740565e-03

e se si calcola la soglia di (s) consigliato su Wikipedia e rimuovere solo quelli che superano la soglia. Per questi dati:

> any(cooks.distance(mod) > 1)
[1] FALSE
> any(cooks.distance(mod) > (4 * nrow(dat)))
[1] FALSE

nessuna delle distanze del cuoco supera le soglie proposte (non sorprende, dato il modo in cui ho generato i dati.)

Dopo aver detto tutto questo, perché si vuole fare questo? Se si sta solo cercando di sbarazzarsi di dati per migliorare la correlazione o generare una relazione significativa, che suona un po 'di pesce e un po' come i dati di dragaggio per me.

Answer 2

Uso method = "spearman" in cor sarà robusta alla contaminazione ed è facile da implementare in quanto coinvolge solo sostituendo cor(x, y) con cor(x, y, method = "spearman").

Ripetendo l'analisi di Prasad ma utilizzando correlazioni Spearman invece troviamo che la correlazione di Spearman è infatti robusta alla contaminazione qui, recuperando il sottostante correlazione nulla:

set.seed(1)

# x and y are uncorrelated
x <- rnorm(1000)
y <- rnorm(1000)
cor(x,y)
## [1] 0.006401211

# add contamination -- now cor says they are highly correlated
x <- c(x, 500)
y <- c(y, 500)
cor(x, y)
## [1] 0.995741

# but with method = "spearman" contamination is removed & they are shown to be uncorrelated
cor(x, y, method = "spearman")
## [1] -0.007270813

Answer 3

Questo può essere stato già evidente per l'OP, ma solo per assicurarsi che ... Bisogna stare attenti perché si cerca di maxmimize correlazione può in realtà tendono a includono valori anomali. (@Gavin toccato questo punto nella sua risposta / commenti.) Sarei prima rimozione valori anomali, poi calcolo di una correlazione. Più in generale, vogliamo essere calcolare una correlazione che è robusto per valori anomali (e ci sono molti tali metodi in R).

Proprio per illustrare questo drammaticamente, creiamo due vettori x e y che non sono correlate:

set.seed(1)
x <- rnorm(1000)
y <- rnorm(1000)
> cor(x,y)
[1] 0.006401211

Ora aggiungiamo un (500,500) punto outlier:

x <- c(x, 500)
y <- c(y, 500)

Ora la correlazione di qualsiasi sottoinsieme che include il punto anomalo sarà vicino al 100%, e la correlazione di qualsiasi sufficientemente grande sottoinsieme che esclude il valore anomalo sarà vicino a zero. In particolare,

> cor(x,y)
[1] 0.995741

Se si vuole stimare una "vera" la correlazione che non è sensibile a valori anomali, si potrebbe provare il pacchetto robust:

require(robust)
> covRob(cbind(x,y), corr = TRUE)
Call:
covRob(data = cbind(x, y), corr = TRUE)

Robust Estimate of Correlation: 
            x           y
x  1.00000000 -0.02594260
y -0.02594260  1.00000000

Si può giocare con i parametri di covRob per decidere come tagliare i dati. UPDATE:. C'è anche la rlm (robusta regressione lineare) nel pacchetto MASS

Answer 4

Ecco un'altra possibilità con i valori erratici catturati. Utilizzando uno schema simile a Prasad:

library(mvoutlier)    
set.seed(1)    
x <- rnorm(1000)    
y <- rnorm(1000)    
xy <- cbind(x, y)    
outliers <- aq.plot(xy, alpha=0.975) #The documentation/default says alpha=0.025.  I think the functions wants 0.975   
cor.plot(x, y)    
color.plot(xy)   
dd.plot(xy)   
uni.plot(xy)    

Nelle altre risposte, 500 è stato bloccato sulla fine del X e Y come un outlier. Questo può, o non può causare un problema di memoria con la vostra macchina, così ho lasciato cadere giù a 4 per evitare che.

x1 <- c(x, 4)     
y1 <- c(y, 4)    
xy1 <- cbind(x1, y1)    
outliers1 <- aq.plot(xy1, alpha=0.975) #The documentation/default says alpha=0.025.  I think the functions wants 0.975
cor.plot(x1, y1)    
color.plot(xy1)    
dd.plot(xy1)    
uni.plot(xy1)    

Ecco le immagini dai X1, Y1, dati xy1:

Answer 5

Si potrebbe provare bootstrapping tuoi dati per trovare il più alto coefficiente di correlazione, per esempio:.

x <- cars$dist
y <- cars$speed
percent <- 0.9         # given in the question above
n <- 1000              # number of resampling
boot.cor <- replicate(n, {tmp <- sample(round(length(x)*percent), replace=FALSE); cor(x[tmp], y[tmp])})

E dopo max(boot.cor) corsa. Non essere delusi se tutti i coefficienti di correlazione saranno tutti uguali:)