Remoção de pontos estranhos perto do centro de um gráfico QQ

14

Estou tentando plotar um QQ-plot com dois conjuntos de dados de cerca de 1,2 milhão de pontos, em R (usando qqplot e alimentando os dados no ggplot2). O cálculo é fácil, mas o gráfico resultante é muito lento para carregar, porque há muitos pontos. Eu tentei a aproximação linear para reduzir o número de pontos para 10000 (isso é o que a função qqplot faz de qualquer maneira, se um dos seus conjuntos de dados é maior que o outro), mas você perde muitos detalhes nos detalhes.

A maioria dos pontos de dados em direção ao centro é basicamente inútil - eles se sobrepõem tanto que provavelmente há cerca de 100 por pixel. Existe alguma maneira simples de remover dados muito próximos, sem perder os dados mais esparsos em direção às caudas?

naught101
fonte
Eu deveria ter mencionado, na verdade, estou comparando um conjunto de dados (observações climáticas) com um conjunto de conjuntos de dados comparáveis ​​(execuções de modelo). Então, na verdade, estou comparando 1,2 milhão de pontos obs, com 87 milhões de pontos do modelo; portanto, a approx()função entra em jogo na qqplot()função.
precisa saber é o seguinte

Respostas:

12

Os gráficos de QQ são incrivelmente autocorrelacionados, exceto nas caudas. Ao analisá-los, concentra-se na forma geral da trama e no comportamento da cauda. Portanto , você se sairá bem subamostrando grosseiramente nos centros das distribuições e incluindo uma quantidade suficiente de caudas.

Aqui está o código que ilustra como fazer a amostra em um conjunto de dados inteiro e como obter 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, esse conjunto de dados simulado mostra uma diferença estrutural entre dois conjuntos de dados de aproximadamente 1,2 milhão de valores, bem como uma quantidade muito pequena de "contaminação" em um deles. Além disso, para tornar esse teste rigoroso, um intervalo de valores é excluído de um dos conjuntos de dados completamente: o gráfico QQ precisa mostrar uma quebra para esses 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 subamostrar 0,1% de cada conjunto de dados e incluir outros 0,1% de seus extremos, fornecendo 2420 pontos para plotagem. O tempo total decorrido é 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")
  )

Nenhuma informação é perdida:

Gráfico de QQ

whuber
fonte
Você não deve mesclar suas respostas?
Michael R. Chernick 28/08/2012
2
@ Michael Sim, normalmente eu teria editado a primeira resposta (a atual). Mas cada resposta é longa e elas usam abordagens substancialmente diferentes, com diferentes características de desempenho, por isso parecia melhor publicar a segunda como uma resposta separada. De fato, fiquei tentado a excluir o primeiro após o segundo (adaptativo) que me ocorreu, mas sua velocidade relativa pode agradar a algumas pessoas, por isso seria injusto removê-lo completamente.
whuber
Isso é basicamente o que eu queria, mas qual é a lógica por trás do uso de sin? Estou certo de que um CDF normal seria uma função melhor se você assumisse que o x era normalmente distribuído? Você acabou de escolher o pecado porque é mais fácil de calcular?
precisa saber é o seguinte
Esses dados devem ser os mesmos da sua outra resposta? Se sim, por que as parcelas são tão diferentes? o que aconteceu com todos os dados para x> 6?
precisa saber é o seguinte
(3-2x)x2
11

Em outras partes deste tópico, propus uma solução simples, mas um tanto ad hoc , de subamostragem dos pontos. É rápido, mas requer algumas experiências para produzir grandes parcelas. A solução a ser descrita é uma ordem de magnitude mais lenta (levando até 10 segundos para 1,2 milhão de pontos), mas é adaptável e automática. Para conjuntos de dados grandes, deve fornecer bons resultados na primeira vez e fazê-lo razoavelmente rapidamente.

Dn

(x,y)ty

Há alguns detalhes a serem resolvidos, especialmente para lidar com conjuntos de dados de diferentes comprimentos. Eu faço isso substituindo o menor pelos quantis correspondentes ao maior: com efeito, uma aproximação linear por partes do EDF do menor é usada em vez de seus valores reais de dados. ("Mais curto" e "mais longo" podem ser revertidos pela configuração use.shortest=TRUE.)

Aqui está uma Rimplementação.

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 exemplo, eu uso dados simulados, como na minha resposta anterior (com um outlier extremamente alto jogado ye um pouco mais de contaminação xnesse período):

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)

Vamos plotar várias versões, usando valores cada vez menores do limite. Com um valor de .0005 e exibindo em um monitor com 1000 pixels de altura, estaríamos garantindo um erro não superior a metade de um pixel vertical em todo o gráfico. Isso é mostrado em cinza (apenas 522 pontos, unidos por segmentos de linha); as aproximações mais grosseiras são plotadas sobre ela: primeiro em preto, depois em vermelho (os pontos vermelhos serão um subconjunto dos pretos e os plotam em excesso), depois em azul (que novamente é um subconjunto e overplot). Os intervalos variam de 6,5 (azul) a 10 segundos (cinza). Dado que eles têm uma escala tão boa, pode-se usar da mesma maneira meio pixel como padrão universal para o limite ( por exemplo , 1/2000 para um monitor com 1000 pixels de altura) e terminar com ele.

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

Gráfico de QQ

Editar

Modifiquei o código original para qqretornar uma terceira coluna de índices na mais longa (ou mais curta, conforme especificado) das duas matrizes originais xe y, correspondendo aos pontos selecionados. Esses índices apontam para valores "interessantes" dos dados e, portanto, podem ser úteis para análises adicionais.

Também removi um erro que ocorria com valores repetidos de x(que betaeram indefinidos).

whuber
fonte
Como calculo qqos argumentos de um determinado vetor? Além disso, você poderia aconselhar sobre como usar sua qqfunção com o ggplot2pacote? Eu estava pensando em usar ggplot2's stat_functionpara isso.
Aleksandr Blekh
10

A remoção de alguns dos pontos de dados no meio alteraria a distribuição empírica e, portanto, o qqplot. Dito isto, você pode fazer o seguinte e plotar diretamente os quantis da distribuição empírica versus os quantis da distribuição 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) 

Você terá que ajustar o seq, dependendo da profundidade que deseja entrar na cauda. Se você quiser ser esperto, também pode afinar essa sequência no meio para acelerar o enredo. Por exemplo, usando

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

é uma possibilidade.

Erik
fonte
Desculpe, não pretendo remover os pontos dos conjuntos de dados, apenas os gráficos.
precisa saber é o seguinte
Mesmo removê-los da trama é uma má idéia. Mas você já tentou alterações de transparência e / ou amostragem aleatória do seu conjunto de dados?
Peter Flom - Restabelece Monica
2
Qual é o problema com a remoção de tinta redundante de pontos sobrepostos na plotagem, @Peter?
whuber
1

Você poderia fazer um hexbinenredo.

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
fonte
Não sei se isso é realmente aplicável aos dados plotados para qq (veja também meu comentário sobre minha pergunta sobre por que isso não funciona no meu caso específico). Ponto interessante embora. Talvez eu veja se consigo trabalhar com modelos individuais versus obs.
precisa saber é o seguinte
1

Outra alternativa é um boxplot paralelo; você disse que tinha dois conjuntos de dados, algo como:

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

e você pode ajustar as várias opções para melhorar seus dados.

Peter Flom - Restabelece Monica
fonte
Nunca fui muito fã de discretizar dados contínuos, mas é uma ideia interessante.
precisa saber é o seguinte