Eliminar puntos extraños cerca del centro de un gráfico QQ

Estoy tratando de trazar un diagrama QQ con dos conjuntos de datos de aproximadamente 1,2 millones de puntos, en R (usando qqplot y alimentando los datos en ggplot2). El cálculo es bastante fácil, pero el gráfico resultante es muy lento de cargar, porque hay muchos puntos. He intentado una aproximación lineal para reducir el número de puntos a 10000 (esto es lo que hace la función qqplot de todos modos, si uno de sus conjuntos de datos es más grande que el otro), pero luego pierde muchos detalles en las colas.

La mayoría de los puntos de datos hacia el centro son básicamente inútiles: se superponen tanto que probablemente haya aproximadamente 100 por píxel. ¿Hay alguna forma simple de eliminar datos que estén demasiado juntos, sin perder los datos más dispersos hacia las colas?

Las gráficas QQ están increíblemente autocorrelacionadas, excepto en las colas. Al revisarlos, uno se enfoca en la forma general de la trama y en el comportamiento de la cola. Ergo , lo harás bien submuestreando groseramente en los centros de las distribuciones e incluyendo una cantidad suficiente de las colas.

Aquí hay un código que ilustra cómo muestrear un conjunto de datos completo y cómo tomar valores extremos.

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
}

Para ilustrar, este conjunto de datos simulado muestra una diferencia estructural entre dos conjuntos de datos de aproximadamente 1,2 millones de valores, así como una muy pequeña cantidad de "contaminación" en uno de ellos. Además, para que esta prueba sea estricta, se excluye un intervalo de valores de uno de los conjuntos de datos por completo: el gráfico QQ debe mostrar un salto para esos valores.

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)

Podemos submuestrear el 0.1% de cada conjunto de datos e incluir otro 0.1% de sus extremos, dando 2420 puntos para trazar. El tiempo total transcurrido es inferior a 0,5 segundos:

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

No se pierde información alguna:

En otra parte de este hilo propuse una solución simple pero algo ad hoc de submuestreo de puntos. Es rápido, pero requiere algo de experimentación para producir grandes parcelas. La solución a punto de describirse es un orden de magnitud más lento (que demora hasta 10 segundos por 1.2 millones de puntos) pero es adaptativo y automático. Para grandes conjuntos de datos, debería dar buenos resultados la primera vez y hacerlo razonablemente rápido.

$D_n$ , la desviación vertical máxima de una línea ajustada. En consecuencia, el algoritmo es este:

$(x,y)$$t$$y$ , reemplace la gráfica con esta línea. De lo contrario, divida los datos en los que preceden al punto de máxima desviación vertical y los que siguen y aplique el algoritmo recursivamente a las dos piezas.

Hay algunos detalles a tener en cuenta, especialmente para hacer frente a conjuntos de datos de diferente longitud. Hago esto reemplazando el más corto por los cuantiles correspondientes al más largo: en efecto, se usa una aproximación lineal por partes del FED del más corto en lugar de sus valores de datos reales. ("Más corto" y "más largo" se pueden revertir configurandouse.shortest=TRUE ).

Aquí hay una Rimplementación.

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
}

Como ejemplo, uso datos simulados como en mi respuesta anterior (con un valor atípico extremadamente alto yy bastante más contaminación en xeste momento):

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)

Tracemos varias versiones, utilizando valores cada vez más pequeños del umbral. Con un valor de .0005 y que se muestra en un monitor de 1000 píxeles de altura, estaríamos garantizando un error de no más de la mitad de un píxel vertical en todas partes del gráfico. Esto se muestra en gris (solo 522 puntos, unidos por segmentos de línea); las aproximaciones más gruesas se trazan en la parte superior: primero en negro, luego en rojo (los puntos rojos serán un subconjunto de los negros y los trazarán en exceso), luego en azul (que nuevamente son un subconjunto y una sobreparcela). Los tiempos varían de 6.5 (azul) a 10 segundos (gris). Dado que escalan tan bien, uno podría usar aproximadamente medio píxel como valor predeterminado universal para el umbral ( por ejemplo , 1/2000 para un monitor de 1000 píxeles de altura) y terminar con él.

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

Editar

He modificado el código original para qqdevolver una tercera columna de índices en la más larga (o más corta, como se especifica) de las dos matrices originales, xyy , correspondiente a los puntos que se seleccionan. Estos índices apuntan a valores "interesantes" de los datos y, por lo tanto, podrían ser útiles para su posterior análisis.

También eliminé un error que ocurría con valores repetidos de x(que causaban betaser indefinidos).

Eliminar algunos de los puntos de datos en el medio cambiaría la distribución empírica y, por lo tanto, el qqplot. Dicho esto, puede hacer lo siguiente y trazar directamente los cuantiles de la distribución empírica frente a los cuantiles de la distribución teórica:

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) 

Tendrás que ajustar la secuencia dependiendo de qué tan profundo quieras llegar a las colas. Si quieres ser inteligente, también puedes reducir esa secuencia en el medio para acelerar la trama. Por ejemplo usando

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

Es una posibilidad.

Podrías hacer 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)

Otra alternativa es un diagrama de caja paralelo; dijiste que tenías dos conjuntos de datos, así que algo como:

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

y puede ajustar las diversas opciones para mejorarlo con sus datos.