¿Métodos estadísticos para trazar datos de manera más eficiente cuando hay millones de puntos presentes?

Creo que R puede tomar mucho tiempo generar gráficos cuando hay millones de puntos presentes, lo que no es sorprendente dado que los puntos se trazan individualmente. Además, tales tramas a menudo son demasiado desordenadas y densas para ser útiles. Muchos de los puntos se superponen y forman una masa negra y se dedica mucho tiempo a trazar más puntos en esa masa.

¿Hay alguna alternativa estadística para representar datos grandes de en un diagrama de dispersión estándar? He considerado una gráfica de densidad, pero ¿qué otras alternativas hay?$n$

Esta es una tarea difícil sin soluciones listas (esto es, por supuesto, porque la gráfica de densidad es una alternativa tan tentadora que a nadie le importa) ¿Entonces que puedes hacer?

Si realmente se superponen (es decir, tienen exactamente las mismas coordenadas X e Y) y no está usando alfa, la mejor idea sería simplemente reducir la superposición usando unique(con alfa, se puede sumar sobre dichos grupos).

De lo contrario, puede redondear manualmente las coordenadas a los píxeles más cercanos y utilizar el método anterior (sin embargo, esta es una solución sucia).

Finalmente, puede hacer un diagrama de densidad solo para usarlo para submuestrear los puntos en las áreas más densas. Esto, por otro lado, no hará exactamente la misma trama y puede introducir artefactos si no se ajusta con precisión.

Mire el paquete hexbin que implementa el papel / método de Dan Carr. La viñeta pdf tiene más detalles que cito a continuación:

1. Información general

La agrupación de hexágonos es una forma de histograma bivariado útil para visualizar la estructura en conjuntos de datos con n grande. El concepto subyacente de hexagon binning es extremadamente simple;

1. El plano xy sobre el conjunto (rango (x), rango (y)) está teselado por una cuadrícula regular de hexágonos.
2. El número de puntos que caen en cada hexágono se cuentan y se almacenan en una estructura de datos.
3. $n \ge 10^6$

Si el tamaño de la cuadrícula y los cortes en la rampa de color se eligen de una manera inteligente, la estructura inherente a los datos debería surgir en las parcelas agrupadas. Las mismas advertencias se aplican a la agrupación hexagonal que se aplican a los histogramas y se debe tener cuidado al elegir los parámetros de agrupación

Debo admitir que no entiendo completamente tu último párrafo:

"No estoy buscando un diagrama de densidad (aunque a menudo son útiles), me gustaría obtener el mismo resultado que una simple llamada de diagrama pero mucho más rápido que millones de diagramas excesivos si es posible".

Tampoco está claro qué tipo de trama (función) está buscando.

Dado que tiene variables métricas, puede encontrar útiles las parcelas agrupadas en hexágono o las parcelas de girasol. Para más referencias, ver

• Gráficos de grandes conjuntos de datos de Unwin / Theus / Hofmann
• Quick-R en " Diagramas de dispersión de alta densidad "
• stat_hexbin de ggplot2

Otra respuesta directa a la pregunta es el paquete rgl, que puede trazar millones de puntos usando OpenGL. Además, especifique un tamaño de punto (por ejemplo, 3) y aleje para ver estos centros de masas como bloques monolíticos, o acerque y vea la estructura de lo que solía ser monolítico: los tamaños de los puntos son constantes pero las distancias entre ellos en la pantalla Depende del zoom. Los niveles alfa también se pueden usar.

Aquí hay un archivo que llamo bigplotfix.R. Si lo obtiene, definirá un contenedor para el plot.xycual "comprime" los datos de la trama cuando es muy grande. El reiniciador no hace nada si la entrada es pequeña, pero si la entrada es grande, la divide en trozos y solo traza el valor máximo y mínimo x e y para cada fragmento. El abastecimiento bigplotfix.Rtambién se vuelve graphics::plot.xya unir para apuntar al contenedor (el abastecimiento varias veces está bien).

Tenga en cuenta que plot.xyes la función de "caballo de batalla" para los métodos de trazado estándar como plot(), lines(), y points(). Por lo tanto, puede continuar utilizando estas funciones en su código sin modificaciones, y sus grandes parcelas se comprimirán automáticamente.

Este es un ejemplo de salida. Es esencialmente plot(runif(1e5)), con puntos y líneas, y con y sin la "compresión" implementada aquí. El gráfico de "puntos comprimidos" pierde la región central debido a la naturaleza de la compresión, pero el gráfico de "líneas comprimidas" se parece mucho más al original sin comprimir. Los tiempos son para el png()dispositivo; por alguna razón, los puntos son mucho más rápidos en el pngdispositivo que en el X11dispositivo, pero las aceleraciones X11son comparables ( X11(type="cairo")fue más lento que X11(type="Xlib")en mis experimentos).

La razón por la que escribí esto es porque estaba cansado de correr plot()por accidente en un gran conjunto de datos (por ejemplo, un archivo WAV). En tales casos, tendría que elegir entre esperar varios minutos para que finalice la trama y finalizar mi sesión de R con una señal (perdiendo así mi historial de comandos y variables recientes). Ahora, si puedo recordar cargar este archivo antes de cada sesión, en realidad puedo obtener una trama útil en estos casos. Un pequeño mensaje de advertencia indica cuando los datos de la trama se han "comprimido".

# bigplotfix.R
# 28 Nov 2016

# This file defines a wrapper for plot.xy which checks if the input
# data is longer than a certain maximum limit. If it is, it is
# downsampled before plotting. For 3 million input points, I got
# speed-ups of 10-100x. Note that if you want the output to look the
# same as the "uncompressed" version, you should be drawing lines,
# because the compression involves taking maximum and minimum values
# of blocks of points (try running test_bigplotfix() for a visual
# explanation). Also, no sorting is done on the input points, so
# things could get weird if they are out of order.
test_bigplotfix = function() {
  oldpar=par();
  par(mfrow=c(2,2))
  n=1e5;
  r=runif(n)
  bigplotfix_verbose<<-T
  mytitle=function(t,m) { title(main=sprintf("%s; elapsed=%0.4f s",m,t["elapsed"])) }
  mytime=function(m,e) { t=system.time(e); mytitle(t,m); }

  oldbigplotfix_maxlen = bigplotfix_maxlen
  bigplotfix_maxlen <<- 1e3;

  mytime("Compressed, points",plot(r));
  mytime("Compressed, lines",plot(r,type="l"));
  bigplotfix_maxlen <<- n
  mytime("Uncompressed, points",plot(r));
  mytime("Uncompressed, lines",plot(r,type="l"));
  par(oldpar);
  bigplotfix_maxlen <<- oldbigplotfix_maxlen
  bigplotfix_verbose <<- F
}

bigplotfix_verbose=F

downsample_xy = function(xy, n, xlog=F) {
  msg=if(bigplotfix_verbose) { message } else { function(...) { NULL } }
  msg("Finding range");
  r=range(xy$x);
  msg("Finding breaks");
  if(xlog) {
    breaks=exp(seq(from=log(r[1]),to=log(r[2]),length.out=n))
  } else {
    breaks=seq(from=r[1],to=r[2],length.out=n)
  }
  msg("Calling findInterval");
  ## cuts=cut(xy$x,breaks);
  # findInterval is much faster than cuts!
  cuts = findInterval(xy$x,breaks);
  if(0) {
    msg("In aggregate 1");
    dmax = aggregate(list(x=xy$x, y=xy$y), by=list(cuts=cuts), max)
    dmax$cuts = NULL;
    msg("In aggregate 2");
    dmin = aggregate(list(x=xy$x, y=xy$y), by=list(cuts=cuts), min)
    dmin$cuts = NULL;
  } else { # use data.table for MUCH faster aggregates
    # (see http://stackoverflow.com/questions/7722493/how-does-one-aggregate-and-summarize-data-quickly)
    suppressMessages(library(data.table))
    msg("In data.table");
    dt = data.table(x=xy$x,y=xy$y,cuts=cuts)
    msg("In data.table aggregate 1");
    dmax = dt[,list(x=max(x),y=max(y)),keyby="cuts"]
    dmax$cuts=NULL;
    msg("In data.table aggregate 2");
    dmin = dt[,list(x=min(x),y=min(y)),keyby="cuts"]
    dmin$cuts=NULL;
    #  ans = data_t[,list(A = sum(count), B = mean(count)), by = 'PID,Time,Site']
  }
  msg("In rep, rbind");
  # interleave rows (copied from a SO answer)
  s <- rep(1:n, each = 2) + (0:1) * n
  xy = rbind(dmin,dmax)[s,];
  xy
}

library(graphics);
# make sure we don't create infinite recursion if someone sources
# this file twice
if(!exists("old_plot.xy")) {
  old_plot.xy = graphics::plot.xy
}

bigplotfix_maxlen = 1e4

# formals copied from graphics::plot.xy
my_plot.xy = function(xy, type, pch = par("pch"), lty = par("lty"),
  col = par("col"), bg = NA, cex = 1, lwd = par("lwd"),
  ...) {

  if(bigplotfix_verbose) {
    message("In bigplotfix's plot.xy\n");
  }

  mycall=match.call();
  len=length(xy$x)
  if(len>bigplotfix_maxlen) {
    warning("bigplotfix.R (plot.xy): too many points (",len,"), compressing to ",bigplotfix_maxlen,"\n");
    xy = downsample_xy(xy, bigplotfix_maxlen, xlog=par("xlog"));
    mycall$xy=xy
  }
  mycall[[1]]=as.symbol("old_plot.xy");

  eval(mycall,envir=parent.frame());
}

# new binding solution adapted from Henrik Bengtsson
# https://stat.ethz.ch/pipermail/r-help/2008-August/171217.html
rebindPackageVar = function(pkg, name, new) {
  # assignInNamespace() no longer works here, thanks nannies
  ns=asNamespace(pkg)
  unlockBinding(name,ns)
  assign(name,new,envir=asNamespace(pkg),inherits=F)
  assign(name,new,envir=globalenv())
  lockBinding(name,ns)
}
rebindPackageVar("graphics", "plot.xy", my_plot.xy);

Tal vez me rechacen por mi método, malos recuerdos de uno de mis profesores de investigación gritando a la gente por tirar buenos datos traduciéndolos en categorías (por supuesto, estoy de acuerdo ahora un día jajaja), no lo sé. De todos modos, si estás hablando de un diagrama de dispersión, entonces he tenido los mismos problemas. Ahora, cuando tengo datos numéricos, no tiene mucho sentido clasificarlos para su análisis. Pero visualizar es una historia diferente. Lo que he encontrado que funciona mejor para mí es primero (1) dividir su variable independiente en grupos usando la función de corte. Puede jugar con el número de grupos, y luego (2) simplemente trazando el DV contra la versión cortada del IV. R generará gráficos de caja en lugar de ese gráfico de dispersión desagradable. Recomiendo eliminar los valores atípicos del diagrama (use la opción del esquema = FALSO en el comando de diagrama). Nuevamente, NUNCA desperdiciaría datos numéricos perfectamente buenos categorizando y luego analizando. Demasiados problemas para hacer eso. Aunque sé que es un tema de debate delicado. Pero hacer eso específicamente con el objetivo de al menos darle un sentido visual a los datos, no hay mucho daño que haya visto. Tracé datos tan grandes como 10M y aún así logré darle sentido a partir de este método. ¡Espero que ayude! ¡Atentamente! Lo he visto. Tracé datos tan grandes como 10M y aún así logré darle sentido a partir de este método. ¡Espero que ayude! ¡Atentamente! Lo he visto. Tracé datos tan grandes como 10M y aún así logré darle sentido a partir de este método. ¡Espero que ayude! ¡Atentamente!

Para series de tiempo grandes, he llegado a amar smoothScatter (parte de la base R no menos). A menudo tengo que incluir algunos datos adicionales, y preservar la API de trama básica es realmente útil, por ejemplo:

set.seed(1)
ra <- rnorm(n = 100000, sd = 1, mean = 0)
smoothScatter(ra)
abline(v=25000, col=2)
text(25000, 0, "Event 1", col=2)

Lo que te da (si perdonas el diseño):

Siempre está disponible y funciona bien con enormes conjuntos de datos, por lo que es bueno al menos echar un vistazo a lo que tienes.