Suppression de points étrangers près du centre d'un tracé QQ

14

J'essaie de tracer un QQ-plot avec deux ensembles de données d'environ 1,2 million de points, dans R (en utilisant qqplot et en introduisant les données dans ggplot2). Le calcul est assez facile, mais le graphique résultant est douloureusement lent à charger, car il y a tellement de points. J'ai essayé l'approximation linéaire pour réduire le nombre de points à 10000 (c'est ce que fait de toute façon la fonction qqplot, si l'un de vos ensembles de données est plus grand que l'autre), mais vous perdez ensuite beaucoup de détails dans les queues.

La plupart des points de données vers le centre sont essentiellement inutiles - ils se chevauchent tellement qu'il y en a probablement environ 100 par pixel. Existe-t-il un moyen simple de supprimer des données trop rapprochées, sans perdre les données les plus clairsemées vers les queues?

rien101
la source
J'aurais dû mentionner que je compare en fait un ensemble de données (observations climatiques) avec un ensemble d'ensembles de données comparables (séries de modèles). Je compare donc en fait 1,2 m de points obs, avec 87 m de points de modèle, d'où la approx()fonction entre en jeu dans la qqplot()fonction.
naught101

Réponses:

12

Les tracés QQ sont incroyablement autocorrélés sauf dans les queues. En les examinant, on se concentre sur la forme globale de l'intrigue et sur le comportement de la queue. Ergo , vous ferez bien en sous-échantillonnant grossièrement au centre des distributions et en incluant une quantité suffisante de queues.

Voici du code illustrant comment échantillonner dans un ensemble de données entier ainsi que comment prendre des valeurs extrêmes.

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
}

Pour illustrer, cet ensemble de données simulé montre une différence structurelle entre deux ensembles de données d'environ 1,2 million de valeurs ainsi qu'une très petite quantité de «contamination» dans l'un d'eux. De plus, pour rendre ce test rigoureux, un intervalle de valeurs est exclu de l'un des ensembles de données: le tracé QQ doit montrer une rupture pour ces valeurs.

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)

Nous pouvons sous-échantillonner 0,1% de chaque ensemble de données et inclure un autre 0,1% de leurs extrêmes, ce qui donne 2420 points à tracer. Le temps total écoulé est inférieur à 0,5 seconde:

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")
  )

Aucune information n'est perdue:

Graphique QQ

whuber
la source
Ne devriez-vous pas fusionner vos réponses?
Michael R. Chernick
2
@Michael Oui, normalement j'aurais édité la première réponse (la présente). Mais chaque réponse est longue et utilise des approches sensiblement différentes, avec des caractéristiques de performance différentes, il a donc semblé préférable de publier la seconde comme réponse distincte. En fait, j'ai été tenté de supprimer le premier après le second (adaptatif) qui m'est venu à l'esprit, mais sa vitesse relative peut plaire à certaines personnes, il serait donc injuste de le supprimer complètement.
whuber
C'est essentiellement ce que je voulais, mais quelle est la justification de l'utilisation de sin? Ai-je raison de dire qu'un CDF normal serait une meilleure fonction, si vous supposiez que le x était normalement distribué? Vous venez de choisir le péché parce qu'il est plus facile à calculer?
naught101
Est-ce censé être les mêmes données que votre autre réponse? Si oui, pourquoi les parcelles sont-elles si différentes? qu'est-il arrivé à toutes les données pour x> 6?
naught101
(3-2X)X2
11

Ailleurs dans ce fil, j'ai proposé une solution simple mais quelque peu ad hoc de sous-échantillonnage des points. Il est rapide, mais nécessite une certaine expérimentation pour produire de superbes parcelles. La solution sur le point d'être décrite est d'un ordre de grandeur plus lent (prenant jusqu'à 10 secondes pour 1,2 million de points) mais est adaptative et automatique. Pour les grands ensembles de données, il devrait donner de bons résultats la première fois et le faire assez rapidement.

n

(X,y)ty

Il y a quelques détails à prendre en compte, en particulier pour gérer des ensembles de données de différentes longueurs. Je fais cela en remplaçant le plus court par les quantiles correspondant au plus long: en effet, une approximation linéaire par morceaux de l'EDF du plus court est utilisée à la place de ses valeurs de données réelles. ("Plus court" et "plus long" peuvent être inversés en réglant use.shortest=TRUE.)

Voici une Rimplémentation.

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
}

À titre d'exemple, j'utilise des données simulées comme dans ma réponse précédente (avec une valeur aberrante extrêmement élevée yet beaucoup plus de contamination à xcette époque):

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)

Tracons plusieurs versions, en utilisant des valeurs de plus en plus petites du seuil. À une valeur de .0005 et affichant sur un moniteur de 1000 pixels de hauteur, nous garantirions une erreur ne dépassant pas un demi-pixel vertical partout sur le tracé. Ceci est indiqué en gris (seulement 522 points, joints par des segments de ligne); les approximations les plus grossières sont tracées dessus: d'abord en noir, puis en rouge (les points rouges seront un sous-ensemble des noirs et les superposeront), puis en bleu (qui sont à nouveau un sous-ensemble et un surplot). Les durées varient de 6,5 (bleu) à 10 secondes (gris). Étant donné qu'ils évoluent si bien, on pourrait tout aussi bien utiliser environ un demi-pixel comme valeur par défaut universelle pour le seuil ( par exemple , 1/2000 pour un moniteur de 1000 pixels de haut) et en finir avec.

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")

Graphique QQ

Éditer

J'ai modifié le code d'origine pour qqrenvoyer une troisième colonne d'index dans le plus long (ou le plus court, comme spécifié) des deux tableaux d'origine, xet y, correspondant aux points sélectionnés. Ces indices pointent vers des valeurs "intéressantes" des données et pourraient donc être utiles pour une analyse plus approfondie.

J'ai également supprimé un bogue se produisant avec des valeurs répétées de x(ce qui a causé betad'être indéfini).

whuber
la source
Comment calculer qqles arguments de pour un vecteur donné? Pouvez-vous également conseiller l'utilisation de votre qqfonction avec le ggplot2package? Je pensais à l' aide ggplot2« est stat_functionpour cela.
Aleksandr Blekh
10

La suppression de certains points de données au milieu modifierait la distribution empirique et donc le qqplot. Cela étant dit, vous pouvez faire ce qui suit et tracer directement les quantiles de la distribution empirique par rapport aux quantiles de la distribution théorique:

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) 

Vous devrez ajuster la séquence en fonction de la profondeur à laquelle vous souhaitez pénétrer dans les queues. Si vous voulez devenir intelligent, vous pouvez également affiner cette séquence au milieu pour accélérer l'intrigue. Par exemple en utilisant

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

est une possibilité.

Erik
la source
Désolé, je ne veux pas dire supprimer les points des ensembles de données, juste des tracés.
naught101
Même les retirer de l'intrigue est une mauvaise idée. Mais avez-vous essayé des modifications de transparence et / ou un échantillonnage aléatoire de votre ensemble de données?
Peter Flom - Réintègre Monica
2
Quel est le problème avec la suppression de l'encre redondante des points qui se chevauchent dans le tracé, @Peter?
whuber
1

Vous pourriez faire un hexbincomplot.

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
la source
Je ne sais pas si c'est vraiment applicable aux données tracées qq (voir aussi mon commentaire sur ma question pour savoir pourquoi cela ne fonctionnera pas pour mon cas spécifique). Point intéressant cependant. Je pourrais voir si je peux le faire fonctionner sur des modèles individuels vs obs.
naught101
1

Une autre alternative est un boxplot parallèle; vous avez dit que vous aviez deux ensembles de données, donc quelque chose comme:

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

et vous pouvez ajuster les différentes options pour améliorer vos données.

Peter Flom - Réintégrer Monica
la source
Je n'ai jamais été un grand fan de la discrétisation de données continues, mais c'est une idée intéressante.
naught101