¿Cómo calcular el solapamiento entre densidades de probabilidad empíricas?

Question

Estoy buscando un método para calcular el área de solapamiento entre dos estimaciones de densidad de kernel en R, como medida de similitud entre dos muestras. Para aclarar, en el siguiente ejemplo, necesitaría cuantificar el área de la región superpuesta de color púrpura:

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)

Se debatió una cuestión similar h La diferencia es que necesito hacer esto para datos empíricos arbitrarios en lugar de distribuciones normales predefinidas. La dirección overlap aborda esta cuestión, pero aparentemente sólo para datos de fecha y hora, lo que no me funciona. El índice de Bray-Curtis (como se implementa en vegan del paquete vegdist(method="bray") ) también parece pertinente, pero de nuevo para datos algo diferentes.

Me interesa tanto el enfoque teórico como las funciones de R que podría emplear para aplicarlo.

Answer 1

El área de solapamiento de dos estimaciones de densidad de kernel puede aproximarse con el grado de precisión que se desee.

1) Dado que las KDE originales probablemente han sido evaluadas sobre alguna rejilla, si la rejilla es la misma para ambas (o puede hacerse fácilmente igual), el ejercicio podría ser tan fácil como simplemente tomar $\min(K_1(x),K_2(x))$ en cada punto y, a continuación, utilizar la regla trapezoidal, o incluso una regla de punto medio.

Si las dos están en cuadrículas diferentes y no se pueden recalcular fácilmente en la misma cuadrícula, se podría recurrir a la interpolación.

2) Podrías encontrar el punto (o puntos) de intersección e integrar la menor de las dos KDEs en cada intervalo donde cada una sea menor. En tu diagrama anterior, integrarías la curva azul a la izquierda de la intersección y la rosa a la derecha por cualquier medio que te guste/dispongas. Esto puede hacerse exactamente considerando el área bajo cada componente del núcleo $\frac{1}{h}K(\frac{x-x_i}{h})$ a la izquierda o a la derecha de ese punto de corte.

Sin embargo Hay que tener muy en cuenta los comentarios de Whuber más arriba.

Answer 2

En aras de la exhaustividad, así es como terminé haciendo esto en 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 ya se ha señalado, hay incertidumbre y subjetividad inherentes en la generación de KDE y también en la integración.

Answer 3

En primer lugar, puede que me equivoque, pero creo que su solución no funcionaría en caso de que haya varios puntos en los que se crucen las estimaciones de densidad del núcleo (KDE). En segundo lugar, aunque el overlap se creó para su uso con datos de marcas de tiempo, puede seguir utilizándolo para estimar el área de solapamiento de dos KDE cualesquiera. Simplemente tienes que reescalar tus datos para que oscilen entre 0 y 2π.
Por ejemplo:

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

Answer 4

Un método alternativo para la estimación empírica es utilizar la tecnología ROC (Receiver Operating Curve) para la estimación. El umbral de Youden nos proporciona una estimación empírica del principal punto de intersección (véase https://journals.lww.com/epidem/Fulltext/2005/01000/Optimal_Cut_point_and_Its_Corresponding_Youden.11.aspx y https://math.stackexchange.com/questions/2404750/intersection-normal-distributions-and-minimal-decision-error/2435957#2435957 ).

El umbral de Youden es el umbral en el que se maximiza la suma de la sensibilidad y la especificidad de la prueba y se minimiza la suma de las tasas de error (tasa de falsos positivos y tasa de falsos negativos). El solapamiento es igual a esta suma mínima de tasas de error.

library(UncertainInterval)
simple_roc2 <- function(ref, test){
  tab = table(test, ref) # head(tab)
  data.frame(threshold=paste('>=',rownames(tab)), 
             ref0 = tab[,1], 
             ref1 = tab[,2],  
             FPR = rev(cumsum(rev(tab[,1])/sum(tab[,1]))), # 1-Sp
             TPR = rev(cumsum(rev(tab[,2])/sum(tab[,2]))), # Se
             row.names=1:nrow(tab))
}
a <- rnorm(10000)
b <- rnorm(10000, 2)
test=c(a,b)
ref=c(rep(0, length(a)), rep(1, length(b)))
# table(test, ref)
res = simple_roc2(ref, test)
res$FNR = 1-res$TPR # 1-Se
pos.optimal.threshold = which.min(res$FPR+res$FNR)
optimal.threshold=row.names(table(test, ref))[pos.optimal.threshold] # Youden threshold
plotMD(ref, test) # library(UncertainInterval) # includes kernel intersection estimate
abline(v=optimal.threshold, col='red')
overlap1(a, b)
(overlap2 = min(res$FPR+res$FNR))

En este caso, esta estimación no paramétrica tiene una ligera tendencia a subestimar el valor real. Esta técnica roc sólo maneja un único punto de intersección (principal). No depende de ninguna distribución específica. Asegúrese de que la distribución b tiene los valores más altos (media(b) > media(a)).

Observando repetidamente los gráficos producidos por plotMD se observa que con 2 * 100 casos el solapamiento de la muestra difiere considerablemente. La mayoría de las diferencias se deben a diferencias de muestra, pero, dependiendo de las distribuciones, todos los métodos tienen condiciones en las que no funcionan correctamente. El uso de un kernel de densidad gaussiana es sensible a los picos en los datos, que pueden ser subestimados. Los métodos de densidad kernel dependen de los parámetros de ajuste que se den a la función de densidad. El método roc no tiene parámetros, pero asume un único punto de intersección. En consecuencia, puede sobrestimar el solapamiento cuando hay un punto de intersección adicional (el punto crítico es la presencia de más de un punto de intersección, no la varianza). Esta sobreestimación puede ser insignificante cuando este punto de intersección secundario se encuentra en las colas de ambas distribuciones.

¿Cómo dar sentido a los distintos métodos y sugerencias? Idear una prueba es más sencillo cuando conocemos el valor verdadero de dos distribuciones. El valor verdadero del solapamiento de las dos distribuciones normales es fácil de calcular. El punto de intersección es simplemente la media de las dos medias de las distribuciones, ya que tienen igual varianza: 1. El solapamiento verdadero es entonces 0,3173105:

(true.overlap = pnorm(1,2,1)+ 1-pnorm(1,0,1))

Véase https://stackoverflow.com/questions/16982146/point-of-intersection-2-normal-curves/45184024#45184024 para un método general de cálculo del punto de intersección de dos distribuciones normales.

En el problema original, hay una mezcla de distribución normal y uniforme. El valor verdadero está en ese caso:

    true.value=sum(pmin(diff(pnorm(0:3)),1/3)) 

Una simulación puede mostrarnos qué método de estimación se aproxima más al valor real:

library(sfsmisc)
overlap1 <- function(a,b){
  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
  total <- integrate.xy(d$x, d$a) + integrate.xy(d$x, d$b)
  intersection <- integrate.xy(d$x, d$w)

  # compute overlap coefficient
  2 * intersection / total
}

library(overlap)
library(scales)
# For explanation of the next function see the answer of S. Venne
overlapEstimates =function(a, b){

  a = data.frame( value = a, Source = "a" )
  b = data.frame( value = b, Source = "b" )
  d = rbind(a, b)

  d$value <- scales::rescale( d$value, to = c(0,2*pi) )

  a <- d[d$Source == "a", 1]
  b <- d[d$Source == "b", 1]

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

nsim=1000; nobs=100; m1=4; sd1=1; m2=6; sd2=1; poi=5
(true.overlap= 1-pnorm(poi, m1, sd1)+pnorm(poi,m2,sd2))
out=matrix(NA,nrow=nsim,ncol=4)
set.seed(0)
for (i in 1:nsim){
  x <- rnorm( nobs, m1, sd1 )
  y <- rnorm( nobs, m2, sd2 )

  out[i,1] = overlap1(x,y)
  out[i,2] = overlapping::overlap(list( x = x, y = y ))$OV
  out[i,3] = overlapEstimates(x,y)['Dhat4']
  out[i,4] = roc.overlap(x,y)
}
(true.overlap=pnorm(poi,m2,sd2)+1-pnorm(poi,m1,sd1))
colMeans(out-true.overlap) # estimation errors
apply(out, 2, sd) # # sd of the estimation errors
apply(out, 2, range)-true.overlap
par(mfrow=c(2,2))
br = seq(-.33,+.33,by=0.05)
hist(out[,1]-true.overlap, breaks=br, ylim=c(0,500), 
     xlim=c(-.33,.33), main='overlap1'); 
abline(v=0, col='red')
hist(out[,2]-true.overlap, breaks=br, ylim=c(0,500), 
     xlim=c(-.33,.33), main='overlapping::overlap')
abline(v=0, col='red')
hist(out[,3]-true.overlap, breaks=br, ylim=c(0,600), 
     xlim=c(-.33,.33), main='overlap::overlapEst')
abline(v=0, col='red')
hist(out[,4]-true.overlap, breaks=br, ylim=c(0,500), 
     xlim=c(-.33,.33), main="ROC estimate"); 
abline(v=0, col='red')

En este caso, especialmente la función overlapping::overlap tiene una tendencia a la (ligera) subestimación, mientras que overlap1 muestra el menor error de estimación. Las estimaciones que utilizan la función de densidad de un modo u otro, pueden producir mejores o peores resultados dependiendo de los parámetros dados a la función de densidad. El método roc no tiene parámetros, lo que puede ser una ventaja.

Siempre es aconsejable observar detenidamente un gráfico de las distribuciones solapadas e idear un método de prueba pertinente para comprobar si la técnica de estimación del solapamiento funciona como se espera para el tipo de datos que se tiene. Es preferible no utilizar técnicas que produzcan sistemáticamente estimaciones demasiado bajas o demasiado altas.