Entfernen von Randpunkten in der Nähe der Mitte eines QQ-Diagramms

14

Ich versuche, ein QQ-Diagramm mit zwei Datensätzen von ungefähr 1,2 Millionen Punkten in R zu zeichnen (unter Verwendung von qqplot und Eingabe der Daten in ggplot2). Die Berechnung ist einfach genug, aber das resultierende Diagramm ist schmerzhaft langsam zu laden, weil es so viele Punkte gibt. Ich habe versucht, die Anzahl der Punkte durch lineare Näherung auf 10000 zu reduzieren (das macht die Funktion qqplot ohnehin, wenn einer Ihrer Datensätze größer als der andere ist), aber dann geht viel Detail in den Endpunkten verloren.

Die meisten Datenpunkte in Richtung des Zentrums sind im Grunde genommen unbrauchbar - sie überlappen sich so sehr, dass es wahrscheinlich ungefähr 100 pro Pixel gibt. Gibt es eine einfache Möglichkeit, Daten zu entfernen, die zu nahe beieinander liegen, ohne dass die spärlicheren Daten zu den Endpunkten hin verloren gehen?

naught101
quelle
Ich hätte erwähnen sollen, dass ich tatsächlich einen Datensatz (Klimabeobachtungen) mit einem Ensemble vergleichbarer Datensätze (Modellläufe) vergleiche. Ich vergleiche also tatsächlich 1,2 Millionen Beobachtungspunkte mit 87 Millionen Modellpunkten, daher kommt die approx()Funktion in der Funktion ins Spiel qqplot().
Naught101

Antworten:

12

QQ-Diagramme sind mit Ausnahme der Schwänze unglaublich automatisch korreliert. Bei der Überprüfung werden die Gesamtform der Handlung und das Schwanzverhalten im Vordergrund stehen. Ergo , Sie werden es gut machen, indem Sie in den Zentren der Verteilungen eine grobe Unterabtastung durchführen und eine ausreichende Menge der Schwänze einbeziehen.

Der folgende Code zeigt, wie ein Sample für einen gesamten Datensatz erstellt und wie Extremwerte ermittelt werden.

quant.subsample <- function(y, m=100, e=1) {
  # m: size of a systematic sample
  # e: number of extreme values at either end to use
  x <- sort(y)
  n <- length(x)
  quants <- (1 + sin(1:m / (m+1) * pi - pi/2))/2
  sort(c(x[1:e], quantile(x, probs=quants), x[(n+1-e):n]))
  # Returns m + 2*e sorted values from the EDF of y
}

Zur Veranschaulichung zeigt dieser simulierte Datensatz einen strukturellen Unterschied zwischen zwei Datensätzen von ungefähr 1,2 Millionen Werten sowie eine sehr geringe Menge an "Kontamination" in einem von ihnen. Um diesen Test stringent zu machen, wird ein Intervall von Werten aus einem der Datasets ausgeschlossen: Der QQ-Plot muss für diese Werte einen Umbruch anzeigen.

set.seed(17)
n.x <- 1.21 * 10^6
n.y <- 1.20 * 10^6
k <- floor(0.0001*n.x)
x <- c(rnorm(n.x-k), rnorm(k, mean=2, sd=2))
x <- x[x <= -3 | x >= -2.5]
y <- rbeta(n.y, 10,13)

Wir können 0,1% jedes Datensatzes subsamplen und weitere 0,1% ihrer Extreme einbeziehen, was 2420 Punkte zum Plotten ergibt. Die verstrichene Gesamtzeit beträgt weniger als 0,5 Sekunden:

m <- .001 * max(n.x, n.y)
e <- floor(0.0005 * max(n.x, n.y))

system.time(
  plot(quant.subsample(x, m, e), 
       quant.subsample(y, m, e), 
       pch=".", cex=4,
       xlab="x", ylab="y", main="QQ Plot")
  )

Es gehen keinerlei Informationen verloren:

QQ Handlung

whuber
quelle
Sollten Sie Ihre Antworten nicht zusammenführen?
Michael R. Chernick
2
@Michael Ja, normalerweise hätte ich die erste Antwort (die vorliegende) bearbeitet. Jede Antwort ist jedoch lang und verwendet ganz unterschiedliche Ansätze mit unterschiedlichen Leistungsmerkmalen. Daher schien es am besten, die zweite Antwort als separate Antwort zu veröffentlichen. Tatsächlich war ich versucht, das erste zu löschen, nachdem mir das zweite (adaptive) einfiel, aber seine relative Geschwindigkeit mag einigen Leuten gefallen, so dass es ungerecht wäre, es insgesamt zu entfernen.
whuber
Dies ist im Grunde, was ich wollte, aber was ist das Grundprinzip hinter der Verwendung von sin? Habe ich recht, dass eine normale CDF eine bessere Funktion wäre, wenn Sie annehmen würden, dass das x normal verteilt ist? Hast du dich gerade für Sünde entschieden, weil es einfacher zu berechnen ist?
Naught101
Sollen dies die gleichen Daten sein wie Ihre andere Antwort? Wenn ja, warum sind die Handlungen so unterschiedlich? Was ist mit allen Daten für x> 6 passiert?
Naught101
(3-2x)x2
11

An anderer Stelle in diesem Thread habe ich eine einfache, aber etwas spontane Lösung für die Unterabtastung der Punkte vorgeschlagen. Es ist schnell, erfordert jedoch einige Experimente, um großartige Diagramme zu erstellen. Die zu beschreibende Lösung ist um eine Größenordnung langsamer (bis zu 10 Sekunden für 1,2 Millionen Punkte), ist jedoch adaptiv und automatisch. Bei großen Datenmengen sollte dies beim ersten Mal zu guten Ergebnissen führen und dies relativ schnell.

Dn

Ermitteln Sie die maximale vertikale Abweichung zwischen der Verbindungslinie der Extrema von (x,y)ty

Es sind einige Details zu beachten, insbesondere, um mit Datensätzen unterschiedlicher Länge fertig zu werden. Ich tue dies, indem ich die kürzere durch die Quantile ersetze, die der längeren entsprechen: In der Tat wird anstelle der tatsächlichen Datenwerte eine stückweise lineare Approximation der EDF der kürzeren verwendet. ("Kürzer" und "länger" können durch Einstellen umgekehrt werden use.shortest=TRUE.)

Hier ist eine RImplementierung.

qq <- function(x0, y0, t.y=0.0005, use.shortest=FALSE) {
  qq.int <- function(x,y, i.min,i.max) {
    # x, y are sorted and of equal length
    n <-length(y)
    if (n==1) return(c(x=x, y=y, i=i.max))
    if (n==2) return(cbind(x=x, y=y, i=c(i.min,i.max)))
    beta <- ifelse( x[1]==x[n], 0, (y[n] - y[1]) / (x[n] - x[1]))
    alpha <- y[1] - beta*x[1]
    fit <- alpha + x * beta
    i <- median(c(2, n-1, which.max(abs(y-fit))))
    if (abs(y[i]-fit[i]) > thresh) {
      assemble(qq.int(x[1:i], y[1:i], i.min, i.min+i-1), 
               qq.int(x[i:n], y[i:n], i.min+i-1, i.max))
    } else {
      cbind(x=c(x[1],x[n]), y=c(y[1], y[n]), i=c(i.min, i.max))
    }
  }
  assemble <- function(xy1, xy2) {
    rbind(xy1, xy2[-1,])
  }
  #
  # Pre-process the input so that sorting is done once
  # and the most detail is extracted from the data.
  #
  is.reversed <- length(y0) < length(x0)
  if (use.shortest) is.reversed <- !is.reversed
  if (is.reversed) {
    y <- sort(x0)
    n <- length(y)
    x <- quantile(y0, prob=(1:n-1)/(n-1))    
  } else {
    y <- sort(y0)
    n <- length(y)
    x <- quantile(x0, prob=(1:n-1)/(n-1))    
  }
  #
  # Convert the relative threshold t.y into an absolute.
  #
  thresh <- t.y * diff(range(y))
  #
  # Recursively obtain points on the QQ plot.
  #
  xy <- qq.int(x, y, 1, n)
  if (is.reversed) cbind(x=xy[,2], y=xy[,1], i=xy[,3]) else xy
}

Als Beispiel verwende ich Daten, die wie in meiner früheren Antwort simuliert wurden (mit einem extrem hohen Ausreißer, der in dieser Zeit stark yverschmutzt ist x):

set.seed(17)
n.x <- 1.21 * 10^6
n.y <- 1.20 * 10^6
k <- floor(0.01*n.x)
x <- c(rnorm(n.x-k), rnorm(k, mean=2, sd=2))
x <- x[x <= -3 | x >= -2.5]
y <- c(rbeta(n.y, 10,13), 1)

Zeichnen wir mehrere Versionen mit immer kleineren Schwellenwerten. Bei einem Wert von .0005 und einer Bildschirmdiagonale von 1000 Pixeln würde ein Fehler von höchstens einem halben vertikalen Pixel überall auf dem Plot garantiert . Dies ist grau dargestellt (nur 522 Punkte, verbunden durch Liniensegmente); Darüber werden die gröberen Näherungen aufgetragen: zuerst in Schwarz, dann in Rot (die roten Punkte sind eine Teilmenge der schwarzen und überzeichnen sie), dann in Blau (die wiederum eine Teilmenge und eine Überzeichnung sind). Die Timings reichen von 6,5 (blau) bis 10 Sekunden (grau). Da sie so gut skalieren, kann man genauso gut ungefähr ein halbes Pixel als universellen Standardwert für den Schwellenwert verwenden ( z. B. 1/2000 für einen 1000 Pixel hohen Monitor) und damit fertig werden.

qq.1 <- qq(x,y)
plot(qq.1, type="l", lwd=1, col="Gray",
     xlab="x", ylab="y", main="Adaptive QQ Plot")
points(qq.1, pch=".", cex=6, col="Gray")
points(qq(x,y, .01), pch=23, col="Black")
points(qq(x,y, .03), pch=22, col="Red")
points(qq(x,y, .1), pch=19, col="Blue")

QQ Handlung

Bearbeiten

Ich habe den ursprünglichen Code geändert qq, um eine dritte Spalte von Indizes in die längste (oder kürzeste, wie angegeben) der ursprünglichen zwei Arrays xund yentsprechend den ausgewählten Punkten zurückzugeben. Diese Indizes verweisen auf "interessante" Werte der Daten und könnten daher für die weitere Analyse nützlich sein.

Ich habe auch einen Fehler behoben, der bei wiederholten Werten von x( betaundefiniert) auftrat.

whuber
quelle
Wie berechne ich qqdie Argumente für einen bestimmten Vektor? Könnten Sie auch raten, Ihre qqFunktion mit ggplot2package zu verwenden? Ich habe darüber nachgedacht, ggplot2's stat_functiondafür zu verwenden.
Aleksandr Blekh
10

Das Entfernen einiger Datenpunkte in der Mitte würde die empirische Verteilung und damit den qqplot verändern. Vor diesem Hintergrund können Sie Folgendes tun und die Quantile der empirischen Verteilung direkt gegen die Quantile der theoretischen Verteilung zeichnen:

x <- rnorm(1200000)
mean.x <- mean(x)
sd.x <- sd(x)
quantiles.x <- quantile(x, probs = seq(0,1,b=0.000001))
quantiles.empirical <- qnorm(seq(0,1,by=0.000001),mean.x,sd.x)
plot(quantiles.x~quantiles.empirical) 

Sie müssen die Sequenz anpassen, je nachdem, wie tief Sie in die Schwänze gelangen möchten. Wenn Sie schlau werden möchten, können Sie diese Sequenz auch in der Mitte verdünnen, um die Handlung zu beschleunigen. Zum Beispiel mit

plogis(seq(-17,17,by=.1))

ist eine Möglichkeit.

Erik
quelle
Entschuldigung, ich meine nicht das Entfernen der Punkte aus den Datensätzen, nur aus den Plots.
Naught101
Sogar das Entfernen von ihnen aus der Handlung ist eine schlechte Idee. Aber haben Sie versucht, Transparenzänderungen und / oder Stichproben aus Ihrem Datensatz zu ermitteln?
Peter Flom - Wiedereinsetzung von Monica
2
Was ist los mit dem Entfernen überflüssiger Tinte von überlappenden Punkten im Plot, @Peter?
Whuber
1

Sie könnten eine hexbinHandlung machen.

x <- rnorm(1200000)
mean.x <- mean(x)
sd.x <- sd(x)
quantiles.x <- quantile(x, probs = seq(0,1,b=0.000001))
quantiles.empirical <- qnorm(seq(0,1,by=0.000001),mean.x,sd.x)

library(hexbin)
bin <- hexbin(quantiles.empirical[-c(1,length(quantiles.empirical))],quantiles.x[-c(1,length(quantiles.x))],xbins=100)
plot(bin)
Roland
quelle
Ich weiß nicht, ob dies wirklich für QQ-geplottete Daten gilt (siehe auch meinen Kommentar zu meiner Frage, warum dies für meinen speziellen Fall nicht funktioniert). Interessanter Punkt. Ich könnte sehen, ob ich es schaffen kann, an einzelnen Modellen gegen obs zu arbeiten.
Naught101
1

Eine andere Alternative ist ein paralleler Boxplot. Sie sagten, Sie hätten zwei Datensätze, also so etwas wie:

y <- rnorm(1200000)
x <- rnorm(1200000)
grpx <- cut(y,20)
boxplot(y~grpx)

und Sie können die verschiedenen Optionen anpassen, um es mit Ihren Daten besser zu machen.

Peter Flom - Setzen Sie Monica wieder ein
quelle
Ich war noch nie ein großer Fan von Diskretisierung kontinuierlicher Daten, aber das ist eine interessante Idee.
Naught101