Bons métodos para gráficos de densidade de variáveis ​​não negativas em R?

36
plot(density(rexp(100))

Obviamente, toda a densidade à esquerda de zero representa viés.

Estou procurando resumir alguns dados para não estatísticos e quero evitar perguntas sobre por que os dados não negativos têm densidade à esquerda de zero. Os gráficos são para verificação aleatória; Quero mostrar as distribuições de variáveis ​​por grupos de tratamento e controle. As distribuições geralmente são exponenciais. Os histogramas são complicados por várias razões.

Uma pesquisa rápida no google me fornece trabalho de estatísticos em núcleos não negativos, por exemplo: isto .

Mas alguma coisa foi implementada no R? Dos métodos implementados, algum deles é "melhor" de alguma forma para estatística descritiva?

EDIT: mesmo que o fromcomando possa resolver meu problema atual, seria bom saber se alguém implementou kernels com base na literatura sobre estimativa de densidade não negativa

generic_user
fonte
3
Não é o que você está perguntando, mas eu não aplicaria a estimativa de densidade do kernel a algo que deveria ser exponencial, especialmente para apresentação a públicos não estatísticos. Eu usaria um gráfico quantil-quantil e explicaria que o gráfico deveria ser reto se a distribuição fosse exponencial.
24413 Nick Cox
6
plot(density(rexp(100), from=0))?
Stéphane Laurent
4
Uma coisa que às vezes tenho feito com bastante êxito é colocar um kde nos logs e depois transformar a estimativa de densidade (sem esquecer o jacobiano). Outra possibilidade seria usar uma estimativa de densidade de log-spline configurada para que ele saiba sobre o limite.
Glen_b -Reinstar Monica
1
Discuti o método de transformação mencionado por @Glen_b em stata-journal.com/sjpdf.html?articlenum=gr0003 (consulte as pp.76-78). Os zeros podem ser acomodados usando log (x + 1) em vez de log e modificando o jacobiano.
Nick Cox

Respostas:

21

Uma solução, emprestada de abordagens para ponderar a borda de estatísticas espaciais, é truncar a densidade à esquerda em zero, mas aumentar a ponderação dos dados mais próximos de zero. A idéia é que cada valor seja "espalhado" em um núcleo da área total da unidade centralizada em x ; qualquer parte do núcleo que se espalhe para território negativo é removida e o núcleo é normalizado para a área da unidade.xx

Por exemplo, com um kernel gaussiano , o peso da renormalização éKh(y,x)=exp(-12((y-x)/h)2)/2π

W(x)=1/0 0K(y,x)dy=11-Φx,h(0 0)

onde é a função de distribuição cumulativa de uma variável normal de média xΦx desvio padrão . Fórmulas comparáveis ​​estão disponíveis para outros kernels.h

Isso é mais simples - e muito mais rápido na computação - do que tentar restringir as larguras de banda perto de . É difícil prescrever exatamente como as larguras de banda devem ser alteradas perto de 0 , de qualquer maneira. No entanto, este método também é0 00 00

Figura

0


Código R

A densityfunção in Rirá reclamar que a soma dos pesos não é unidade, porque deseja que a integral sobre todos os números reais seja unidade, enquanto essa abordagem torna a integral sobre números positivos igual à unidade. Como verificação, a última integral é estimada como uma soma de Riemann.

set.seed(17)
x <- rexp(1000)
#
# Compute a bandwidth.
#
h <- density(x, kernel="gaussian")$bw # $
#
# Compute edge weights.
#
w <- 1 / pnorm(0, mean=x, sd=h, lower.tail=FALSE)
#
# The truncated weighted density is what we want.
#
d <- density(x, bw=h, kernel="gaussian", weights=w / length(x))
d$y[d$x < 0] <- 0
#
# Check: the integral ought to be close to 1:
#
sum(d$y * diff(d$x)[1])
#
# Plot the two density estimates.
#
par(mfrow=c(1,1))
plot(d, type="n", main="Default and truncated densities", xlim=c(-1, 5))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
whuber
fonte
21

Uma alternativa é a abordagem de Kooperberg e colegas, com base na estimativa da densidade usando splines para aproximar a densidade de log dos dados. Vou mostrar um exemplo usando os dados da resposta do @ whuber, que permitirá uma comparação de abordagens.

set.seed(17)
x <- rexp(1000)

Você precisará do pacote de logspline instalado para isso; instale-o se não estiver:

install.packages("logspline")

Carregue o pacote e estime a densidade usando a logspline()função:

require("logspline")
m <- logspline(x)

A seguir, presumo que o objeto dda resposta do @ whuber esteja presente no espaço de trabalho.

plot(d, type="n", main="Default, truncated, and logspline densities", 
     xlim=c(-1, 5), ylim = c(0, 1))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
plot(m, add = TRUE, col = "red", lwd = 3, xlim = c(-0.001, max(x)))
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)

O gráfico resultante é mostrado abaixo, com a densidade da linha de logs mostrada pela linha vermelha

Densidades padrão, truncadas e de linha de logs

Além disso, o suporte para a densidade pode ser especificado através de argumentos lbounde ubound. Se quisermos supor que a densidade é 0 à esquerda de 0 e há uma descontinuidade em 0, poderíamos usar lbound = 0a chamada paralogspline() , por exemplo

m2 <- logspline(x, lbound = 0)

Rendimento da seguinte estimativa de densidade (mostrada aqui com o majuste original da linha de logs, pois a figura anterior já estava ficando ocupada).

plot.new()
plot.window(xlim = c(-1, max(x)), ylim = c(0, 1.2))
title(main = "Logspline densities with & without a lower bound",
      ylab = "Density", xlab = "x")
plot(m,  col = "red",  xlim = c(0, max(x)), lwd = 3, add = TRUE)
plot(m2, col = "blue", xlim = c(0, max(x)), lwd = 2, add = TRUE)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)
axis(1)
axis(2)
box()

O gráfico resultante é mostrado abaixo

Comparação das estimativas de densidade da linha de produção com e sem um limite inferior no suporte

xx=0 0x

Restabelecer Monica - G. Simpson
fonte
1
0 01
@whuber Boa pergunta. Só me deparei com essa abordagem recentemente. Eu suspeito que uma boa pergunta a ser feita aqui é que, como os métodos truncados e os logspline são apenas estimativas da densidade real, as diferenças de ajuste são significativas, estatisticamente? Não sei exatamente por que ele se sai tão bem em zero, no entanto. Gostaria muito de saber por que também.
Reinstate Monica - G. Simpson
@ GavinSimpson, Obrigado por esta resposta agradável. Você pode reproduzir o último enredo com a versão mais recente de logspline? Para mim, a densidade de ambos, a versão limitada e a ilimitada chega a zero em x = 0.
cel
4

Para comparar distribuições por grupos (que você diz ser o objetivo em um de seus comentários), por que não algo mais simples? Os gráficos de caixas paralelas funcionam bem se N for grande; os gráficos de faixas paralelas funcionam se N for pequeno (e ambos mostrarem valores discrepantes bem, o que você diz ser um problema em seus dados).

Peter Flom - Restabelece Monica
fonte
1
Sim, obrigado, isso funciona. Mas eu gosto de gráficos de densidade. Eles mostram mais sobre os dados do que os boxplots. Acho que estou meio surpreso que nada parece já ter sido implementado. Talvez eu mesmo implemente uma dessas coisas um dia. As pessoas provavelmente achariam útil.
generic_user
1
Também gosto de gráficos de densidade; mas você deve considerar seu público.
Peter Flom - Restabelece Monica
1
Tem que concordar com @PeterFlom neste. Não fique muito complicado se o seu público não tiver conhecimento estatístico. Você também pode fazer gráficos de caixa comparativos / paralelos com uma sobreposição de gráficos de borboleta na parte superior. Dessa forma, o resumo do gráfico de caixa é visível, assim como todos os dados.
Doug.numbers
A sugestão de que pessoas diferentes compreendem parcelas agregadas de maneira diferente é certamente correta. Apesar de entender o que é um gráfico de densidade (e entender que não é uma probabilidade), não compreendo o que possa ser um "boxplot paralelo". Ele sugere um gráfico de coordenadas paralelas, mas suspeito que isso não esteja correto.
Dwin
2

Como Stéphane comenta, você pode usar from = 0e, além disso, pode representar seus valores sob a curva de densidade comrug (x)

Aghila
fonte
4
Corrija-me se estiver errado, mas from=0parece que apenas suprime a plotagem para valores abaixo de 0; não corrige o cálculo pelo fato de que parte da distribuição foi manchada abaixo de 0.
Nick Cox
1
Está correto. O uso do fromcomando gera um gráfico que parece ter o pico exatamente à direita de zero. Mas se você observar histogramas com compartimentos continuamente menores, muitos dados mostrarão o pico em zero. O fromé apenas um truque gráfico.
generic_user
@ NickCox Não tenho certeza, mas acho que não from=0suprime nada. Apenas inicia a "grade" em zero.
Stéphane Laurent
A diferença é se a densidade estimada é diferente de zero para valores negativos, não se é ou não plotada. Os pesquisadores podem decidir não se preocupar com isso se tudo o que eles querem é uma visualização.
27613 Nick Cox
@NickCox O comando density(rexp(100), from=0)não tem nada a ver com o gráfico
Stéphane Laurent