Rimuovere i valori anomali dal calcolo del coefficiente di correlazione

Question

Supponiamo di disporre di due vettori numerici x e . 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 (ad esempio 90%) come massimizzare il coefficiente di correlazione?

Answer 1

Se davvero vuole fare questo (togliere i residui più grandi (assoluti)), allora possiamo utilizzare il modello lineare per stimare il minimo soluzione quadrata e residui associati e quindi selezionare il mezzo% dei dati. Ecco un esempio:

In primo luogo, generare alcuni dati dummy:

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() noi la richiesta può dare quantili dei residui. È suggerito conservando il 90% dei dati, quindi vogliamo le superiori e inferiori 0,05 quantili:

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

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

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

Possiamo 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 i dati completi e il sottoinsieme selezionato sono:

> 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 

essere consapevoli del fatto che qui ci potrebbe essere buttare fuori perfettamente buoni dati, perché abbiamo appena scelto il 5%, con grandi residui positivi e 5% con il più grande negativo. Un'alternativa consiste nel selezionare il 90% con piccole assoluti residui:

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 

altro punto è che anche allora stiamo gettando buoni dati. Potresti voler considerare la distanza di Cook come una misura della forza dei valori anomali e scartare solo quei valori al di sopra di una certa soglia della distanza di Cook.Wikipedia ha informazioni sulla distanza di Cook e le 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 inferiore alle soglie proposte (. Non sorprende, dato il modo in cui ho generato i dati)

Detto tutto questo, perché vuoi fare questo? Se stai solo cercando di sbarazzarti dei dati per migliorare una correlazione o generare una relazione significativa, sembra un po 'strano e un po' come i dati che mi draghano.

Answer 2

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

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 l'esecuzione max(boot.cor). Non siate delusi se tutti i coefficienti di correlazione saranno tutti uguali :)

Answer 3

Questo potrebbe essere stato già ovvio per l'OP, ma solo per essere sicuri ... Bisogna fare attenzione perché provare a massimizzare la correlazione può effettivamente tendere a includere valori anomali di. (@Gavin ha toccato questo punto nella sua risposta/commenti.) Sarei primo rimuovere i valori anomali, quindi calcolare una correlazione. Più in generale, vogliamo calcolare una correlazione che sia robusta rispetto ai valori anomali (e ci sono molti di questi 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 punto outlier (500,500):

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 sottoinsieme sufficientemente grande 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 a decidere come tagliare i dati. UPDATE: C'è anche il rlm (regressione lineare robusta) nel pacchetto MASS.

Answer 4

utilizzando method = "spearman" in cor sarà robusto 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 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 5

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)  

In altre risposte, 500 era bloccato sull'estremità di x ed y come valore aberrante. Ciò potrebbe o non potrebbe causare un problema di memoria con il tuo computer, quindi l'ho lasciato cadere a 4 per evitarlo.

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 dalla X1, Y1, dati xy1: