Como calcular a sobreposição entre densidades empíricas de probabilidade?

14

Estou procurando um método para calcular a área de sobreposição entre duas estimativas de densidade de kernel em R, como uma medida de similaridade entre duas amostras. Para esclarecer, no exemplo a seguir, eu precisaria quantificar a área da região sobreposta arroxeada:

library(ggplot2)
set.seed(1234)
d <- data.frame(variable=c(rep("a", 50), rep("b", 30)), value=c(rnorm(50), runif(30, 0, 3)))
ggplot(d, aes(value, fill=variable)) + geom_density(alpha=.4, color=NA)

insira a descrição da imagem aqui

Uma pergunta semelhante foi discutida aqui , a diferença é que eu preciso fazer isso para dados empíricos arbitrários, em vez de distribuições normais predefinidas. O overlappacote aborda essa questão, mas aparentemente apenas para dados de registro de data e hora, o que não funciona para mim. O índice Bray-Curtis (conforme implementado na função vegando pacote vegdist(method="bray")) também parece relevante, mas novamente para dados um pouco diferentes.

Estou interessado na abordagem teórica e nas funções R que eu poderia empregar para implementá-la.

mmk
fonte
2
"quantificar a área roxa" é um problema na estimativa, não no teste de hipóteses; portanto, você não pode esperar "realizar isso usando um teste estatístico citável padrão ". Você se contradiz. Por favor, esclareça o que você realmente deseja. Se tudo o que você deseja é uma estimativa da área de sobreposição de dois KDEs, esse é um cálculo simples.
Glen_b -Reinstate Monica
@Glen_b obrigado pelo comentário, ajudou a esclarecer meu pensamento não estatístico. Acredito que a área de sobreposição entre os KDEs é realmente o que estou procurando - editei a pergunta para refletir isso.
mmk 14/05
2
Eu ficaria muito preocupado com o risco de arbitrariedade neste método. Dependendo da largura de banda kernel, a sobreposição computadorizada entre quaisquer dois conjuntos de dados poderia ser feito para igualar qualquer valor escolhido no intervalo . As larguras de banda padrão não são otimizadas para essa finalidade e, portanto, concebivelmente, podem dar resultados surpreendentes, arbitrários ou inconsistentes. Conjuntos de dados com limites naturais (como dados ou proporções não negativos, etc.) introduziriam ainda efeitos de borda indesejados. O que fazer em vez disso? Comece com o motivo desse cálculo: o que essa "semelhança" pretende significar? (0 0,1)
whuber
A mesma pergunta apareceu alguns meses depois, mas se referia a pontos de interseção, no entanto, havia algumas notas válidas que poderiam ser levadas em consideração. Na questão referida, trata-se de duas distribuições empíricas. Eu adiciono o link, pois este post só responde isso via estimativa de densidade do kernel e para distribuições normais. Acho que o link abaixo se estende à questão de pares de distribuições empíricas. stats.stackexchange.com/questions/122857/… - Barnaby há 7 horas
Barnaby

Respostas:

9

A área de sobreposição de duas estimativas de densidade de kernel pode ser aproximada para qualquer grau de precisão desejado.

1) Como os KDEs originais provavelmente foram avaliados em alguma grade, se a grade é a mesma para ambos (ou pode ser facilmente igual), o exercício pode ser tão fácil quanto simplesmente tomar em cada ponto e, em seguida, usando a regra trapezoidal ou mesmo uma regra do ponto médio.min(K1(x),K2(x))

Se os dois estiverem em grades diferentes e não puderem ser recalculados facilmente na mesma grade, a interpolação poderá ser usada.

2) Você pode encontrar o (s) ponto (s) de interseção e integrar o mais baixo dos dois KDEs em cada intervalo em que cada um é mais baixo. No diagrama acima, você integraria a curva azul à esquerda da interseção e a rosa à direita por qualquer meio que desejar / disponível. Isso pode ser feito essencialmente exatamente considerando a área sob cada componente do kernel à esquerda ou à direita desse ponto de corte.1hK(x-xEuh)

No entanto , os comentários acima devem ser claramente lembrados - isso não é necessariamente uma coisa muito significativa a se fazer.

Glen_b -Reinstate Monica
fonte
Como você calcula o erro associado ao método um e ao método 2?
Olliepower 14/05
Em circunstâncias normais, ambos serão minúsculos em comparação com o erro nas estimativas de densidade do kernel, portanto, não me preocuparia muito. Os limites de erro podem ser calculados em métodos trapezoidais e outras integrações numéricas, é claro - esses cálculos são bastante padrão - mas é inútil se preocupar, pois os KDEs têm grandes incertezas. O método 2 será preciso no erro de arredondamento acumulado dos cálculos.
Glen_b -instala Monica
1
Essas sugestões de metodologia fazem sentido, muito obrigado pela sua resposta. Vou trabalhar para implementar isso em R, mas como iniciante, eu estaria interessado em sugestões sobre como codificar isso de maneira limpa.
mmk 14/05
10

Por uma questão de integridade, eis como eu acabei fazendo isso no R:

# simulate two samples
a <- rnorm(100)
b <- rnorm(100, 2)

# define limits of a common grid, adding a buffer so that tails aren't cut off
lower <- min(c(a, b)) - 1 
upper <- max(c(a, b)) + 1

# generate kernel densities
da <- density(a, from=lower, to=upper)
db <- density(b, from=lower, to=upper)
d <- data.frame(x=da$x, a=da$y, b=db$y)

# calculate intersection densities
d$w <- pmin(d$a, d$b)

# integrate areas under curves
library(sfsmisc)
total <- integrate.xy(d$x, d$a) + integrate.xy(d$x, d$b)
intersection <- integrate.xy(d$x, d$w)

# compute overlap coefficient
overlap <- 2 * intersection / total

Como observado, há incerteza e subjetividade inerentes envolvidas na geração do KDE e também na integração.

mmk
fonte
2
Agora existe um pacote chamado CRAN overlappingque estima a área de sobreposição de 2 (ou mais) distribuições empíricas. Confira a documentação aqui: rdocumentation.org/packages/overlapping/versions/1.5.0/topics/...
Stefan Avey
x,dx,dx,d
@mmk você pode fazer isso para densidades 2D?
No Lie
4

Primeiro, posso estar errado, mas acho que sua solução não funcionaria caso houvesse múltiplos pontos em que as estimativas de densidade de kernel (KDE) se cruzam. Segundo, embora o overlappacote tenha sido criado para ser usado com dados de registro de data e hora, você ainda pode usá-lo para estimar a área de sobreposição de dois KDEs. Você simplesmente precisa redimensionar seus dados para que eles variam de 0 a 2π.
Por exemplo :

# simulate two sample    
 a <- rnorm(100)
 b <- rnorm(100, 2)

# To use overplapTrue(){overlap} the scale must be in radian (i.e. 0 to 2pi)
# To keep the *relative* value of a and b the same, combine a and b in the
# same dataframe before rescaling. You'll need to load the ‘scales‘ library.
# But first add a "Source" column to be able to distinguish between a and b
# after they are combined.
 a = data.frame( value = a, Source = "a" )
 b = data.frame( value = b, Source = "b" )
 d = rbind(a, b)
 library(scales) 
 d$value <- rescale( d$value, to = c(0,2*pi) )

# Now you can created the rescaled a and b vectors
 a <- d[d$Source == "a", 1]
 b <- d[d$Source == "b", 1]

# You can then calculate the area of overlap as you did previously.
# It should give almost exactly the same answers.
# Or you can use either the overlapTrue() and overlapEst() function 
# provided with the overlap packages. 
# Note that with these function the KDE are fitted using von Mises kernel.
 library(overlap)
  # Using overlapTrue():
   # define limits of a common grid, adding a buffer so that tails aren't cut off
     lower <- min(d$value)-1 
     upper <- max(d$value)+1
   # generate kernel densities
     da <- density(a, from=lower, to=upper, adjust = 1)
     db <- density(b, from=lower, to=upper, adjust = 1)
   # Compute overlap coefficient
     overlapTrue(da$y,db$y)


  # Using overlapEst():            
    overlapEst(a, b, kmax = 3, adjust=c(0.8, 1, 4), n.grid = 500)

# You can also plot the two KDEs and the region of overlap using overlapPlot()
# but sadly I haven't found a way of changing the x scale so that the scale 
# range correspond to the initial x value and not the rescaled value.
# You can only change the maximum value of the scale using the xscale argument 
# (i.e. it always range from 0 to n, where n is set with xscale = n).
# So if some of your data take negative value, you're probably better off with
# a different plotting method. You can change the x label with the xlab
# argument.  
  overlapPlot(a, b, xscale = 10, xlab= "x metrics", rug=T)
S. Venne
fonte