2011-01-12 11 views
8

Angenommen, wir haben zwei numerische Vektoren x und y. Der Pearson-Korrelationskoeffizient zwischen x und y wird durchAusreißer aus der Korrelationskoeffizientenberechnung entfernen

gegeben

cor (x, y)

Wie kann ich nur eine Teilmenge von x und y bei der Berechnung automatisch berücksichtigen (zB 90%) als um den Korrelationskoeffizienten zu maximieren?

+0

Was halten Sie einen Ausreißer hier in Betracht ziehen? Abweichung von der Fit-Linie der kleinsten Quadrate (d. H. Größte Residuen) oder Werte an den Extremen der bivariaten Verteilung von "x" und "y"? –

+0

@Gavin Hier betrachte ich die größten Residuen als Ausreißer. – Leo

Antwort

22

Wenn Sie wirklich wollen dies tun (entfernen Sie die größten (absoluten) Residuen), dann können wir das lineare Modell verwenden, um die geringste Schätzung Quadrate Lösung und zugehörige Residuen und wählen Sie dann die mittleren n% der Daten. Hier ein Beispiel:

Zunächst einige Dummy-Daten erzeugen:

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) 

Als nächst wir das lineare Modell passen und die Residuen extrahieren:

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

Die quantile() Funktion kann uns die erforderlichen Quantile der Residuen. Sie schlugen vor 90% der Daten beibehalten wird, so wollen wir die oberen und unteren 0,05 quantiles:

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

Wählen Sie diese Beobachtungen mit Residuen in der Mitte 90% der Daten:

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

Wir können dann visualisieren diese, mit den roten Punkten diejenigen sind wir behalten:

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) 

The plot produced from the dummy data showing the selected points with the smallest residuals

Die Korrelationen für die vollständigen Daten und der ausgewählten Teilmenge sind:

> 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 

Beachten Sie, dass hier könnten wir ganz gute Daten werfen werden, weil wir die 5% mit der größten positiven Residuen und 5% mit der gerade wählen größtes Negativ. Eine Alternative ist es, das 90% mit kleinsten absoluten Residuen wählen:

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) 

Mit dieser etwas anderen Teilmenge ist die Korrelation etwas niedriger:

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

Ein weiterer Punkt ist, dass selbst dann wir werfen gute Daten aus. Sie sollten sich Cooks Entfernung als Maß für die Stärke der Ausreißer ansehen und nur die Werte oberhalb einer bestimmten Cook-Distanz verwerfen.Wikipedia hat Informationen über Cooks Entfernung und vorgeschlagene Schwellenwerte.

> 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 

und wenn Sie die Schwelle (n) vorgeschlagen diejenigen, die die Schwelle überschreiten nur auf Wikipedia, und entfernen Sie berechnen: Die cooks.distance() Funktion können die Werte von mod abzurufen verwendet werden. Für diese Daten:

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

keiner der Abstände Cook übersteigen die vorgeschlagenen Schwellenwerte (. Nicht überraschend angesichts der Art, wie ich die Daten generiert)

all dies gesagt ist, warum tun Sie dies tun wollen? Wenn Sie nur versuchen, Daten loszuwerden, um eine Korrelation zu verbessern oder eine signifikante Beziehung zu generieren, klingt das ein bisschen fischig und etwas wie Daten, die mich ausgraben.

+0

Vielen Dank für solch eine ausgezeichnete Antwort! Der Grund, warum ich das tun möchte, ist folgender. Ich benchmarkiere verschiedene Methoden zur Vorhersage experimenteller Beobachtungen (Änderungen der Bindungsenergie bei Mutation eines Proteinkomplexes) basierend auf experimentellen Strukturen der Komplexe. Die Zielwerte stammen aus verschiedenen Quellen mit unterschiedlicher Qualität. Und Fehler in den Strukturen können die Vorhersagen stark beeinflussen. Ich habe also einige Ausreißer, aber wenn ich eine "bereinigte" Korrelation für verschiedene Methoden ansehe, kann ich leichter die Methode auswählen, die am besten für die günstigen Fälle funktioniert. – Leo

2

Sie könnten versuchen, Ihre Daten Bootstrapping die höchsten Korrelationskoeffizienten zu finden, z.B .:

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])}) 

Und nach Laufe max(boot.cor). Seien Sie nicht enttäuscht, wenn alle Korrelationskoeffizienten alle gleich sind :)

9

Dies könnte schon offensichtlich für das OP gewesen sein, aber nur um sicherzustellen, dass ... Sie müssen vorsichtig sein, weil versuchen zu maxmimize Korrelation tatsächlich dazu neigen, gehören Ausreißer. (@ Gavin berührt diesen Punkt in seiner Antwort/Kommentare.) Ich wäre zuerst Entfernen von Ausreißern, dann Berechnung einer Korrelation. Allgemeiner wollen wir eine Korrelation zu Ausreißern berechnen (und es gibt viele solcher Methoden in R).

Nur diese dramatisch zu illustrieren, lassen Sie uns zwei Vektoren erzeugen x und y, die unkorreliert sind:

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

Lassen Sie uns jetzt (500,500) einen Ausreißer Punkt hinzufügen:

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

Nun ist die Korrelation von jeder Die Teilmenge, die den Ausreißerpunkt enthält, liegt nahe bei 100% und die Korrelation einer ausreichend großen Teilmenge, die den Ausreißer ausschließt, wird sein nahe bei Null. Insbesondere

> cor(x,y) 
[1] 0.995741 

Wenn Sie ein „true“ Korrelation schätzen wollen, die nicht empfindlich auf Ausreißer ist, können Sie versuchen, das robust Paket:

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 

Sie mit den Parametern von covRob rumspielen können Entscheiden Sie, wie die Daten zu trimmen sind. UPDATE: Es gibt auch die rlm (robuste lineare Regression) im MASS Paket.

+0

+1 Schöne Antwort Prasad. –

15

Mit method = "spearman" in cor wird Kontamination robust sein und ist einfach zu implementieren, da es sich um nur cor(x, y) mit cor(x, y, method = "spearman") ersetzen.

Wiederholung Prasad Analyse jedoch unter Verwendung von Spearman Korrelationen stattdessen finden wir, dass die Spearman Korrelation hier auf die Kontamination in der Tat robust ist, die zugrunde liegende Nullkorrelation erholt:

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 
+1

+1 für den Hinweis auf 'spearman' –

+0

' spearman' auf einige Arten von Verunreinigungen robust sein, nämlich einzelne hohen Wert Punkte korreliert perfekt ist in einem aufgeblasenen 'pearson' Korrelation ergibt. Es ist jedoch nicht vollständig robust gegenüber Verunreinigungen durch Ausreißer am unteren Ende der Skala. – cashoes

4

Hier ist eine andere Möglichkeit, mit den erfassten Ausreißer.Mit einem ähnlichen Schema wie 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 den anderen Antworten, 500 wurde am Ende von x und y als Ausreißer stecken. Das kann oder kann kein Speicherproblem mit Ihrem Computer verursachen, also habe ich es auf 4 reduziert, um das zu vermeiden.

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)  

Hier sind die Bilder aus dem x1, y1, xy1 Daten:

alt text

alt text

alt text

+3

ich per E-Mail den Maintainer für mvoutlier über das Problem, das ich in der oben aq.plot() Aussagen mit alpha hat. Er hat festgelegt, da das Problem und aktualisiert mvoutlier auf Version 1.6 (aktualisiert 14. Januar 2011) http://cran.r-project.org/web/packages/mvoutlier/index.html –