Comment calculer le chevauchement entre les densités de probabilité empiriques?

14

Je cherche une méthode pour calculer la zone de chevauchement entre deux estimations de densité de noyau dans R, comme mesure de similitude entre deux échantillons. Pour clarifier, dans l'exemple suivant, il me faudrait quantifier l'aire de la région de chevauchement violacé:

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)

entrez la description de l'image ici

Une question similaire a été discutée ici , la différence étant que je dois le faire pour des données empiriques arbitraires plutôt que pour des distributions normales prédéfinies. Le overlappackage répond à cette question, mais apparemment uniquement pour les données d'horodatage, ce qui ne fonctionne pas pour moi. L'indice Bray-Curtis (tel qu'implémenté dans veganla vegdist(method="bray")fonction du package ) semble également pertinent mais là encore pour des données quelque peu différentes.

Je m'intéresse à la fois à l'approche théorique et aux fonctions R que je pourrais utiliser pour l'implémenter.

mmk
la source
2
«quantifier la zone violette» est un problème d'estimation, pas de test d'hypothèse, vous ne pouvez donc pas espérer «accomplir cela en utilisant un test statistique standard citable ». Vous vous contredisez. Veuillez clarifier ce que vous voulez réellement . Si tout ce que vous voulez est une estimation de la zone de chevauchement de deux KDE, c'est un calcul simple.
Glen_b -Reinstate Monica
@Glen_b merci pour le commentaire, a aidé à clarifier ma pensée non statisticienne. Je crois que le domaine de chevauchement entre les KDE est en effet ce que je recherche - j'ai édité la question pour refléter cela.
mmk
2
Je serais très préoccupé par le risque d'arbitraire dans cette méthode. En fonction de la bande passante du noyau, le chevauchement calculé entre les deux ensembles de données pourrait être faite pour égaler une valeur choisie dans l'intervalle . Les bandes passantes par défaut ne sont pas optimisées à cet effet et pourraient donc donner des résultats surprenants, arbitraires ou incohérents. Des ensembles de données avec des limites naturelles (telles que des données ou des proportions non négatives, etc.) introduiraient en outre des effets de bord indésirables. Que faire à la place? Commençons par la raison de ce calcul: que signifie cette "similitude"? (0,1)
whuber
La même question est apparue quelques mois plus tard mais faisait référence à des points d'intersection mais il y avait quelques notes valables qui pouvaient être prises en considération. Dans la question référée, il s'agit de deux distributions empiriques. J'ajoute le lien car cet article ne répond qu'à cela via l'estimation de la densité du noyau et pour les distributions normales. Je pense que le lien ci-dessous s'étend sur la question des paires de distributions empiriques. stats.stackexchange.com/questions/122857/… - Barnaby il y a 7 heures
Barnaby

Réponses:

9

La zone de chevauchement de deux estimations de densité de noyau peut être approximée à tout degré de précision souhaité.

min(K1(x),K2(x))

Si les deux sont sur des grilles différentes et ne peuvent pas être facilement recalculés sur la même grille, une interpolation pourrait être utilisée.

1hK(xxih)

Cependant , les commentaires de whuber ci-dessus doivent être clairement gardés à l'esprit - ce n'est pas nécessairement une chose très significative à faire.

Glen_b -Reinstate Monica
la source
Comment calculez-vous l'erreur associée à la méthode un et à la méthode 2?
olliepower
Dans des circonstances normales, les deux seront minuscules par rapport à l'erreur dans les estimations de densité du noyau, donc je ne m'inquiéterais pas trop. Les limites d'erreur peuvent être calculées sur des méthodes trapézoïdales et d'autres intégrations numériques bien sûr - de tels calculs sont assez standard - mais il est inutile de s'inquiéter étant donné que les KDE ont de grandes incertitudes. La méthode 2 sera précise à l'erreur d'arrondi cumulée des calculs.
Glen_b -Reinstate Monica
1
Ces suggestions de méthodologie ont du sens, merci beaucoup pour votre réponse. Je vais travailler sur l'implémentation de cela dans R, mais en tant que novice, je serais intéressé par des suggestions sur la façon de coder cela proprement.
mmk
10

Par souci d'exhaustivité, voici comment j'ai fini par le faire dans 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

Comme indiqué, il existe une incertitude et une subjectivité inhérentes à la génération de KDE ainsi qu'à l'intégration.

mmk
la source
2
Il existe maintenant un package sur CRAN appelé overlappingqui estime l'aire du chevauchement de 2 (ou plus) distributions empiriques. Consultez la documentation ici: rdocumentation.org/packages/overlapping/versions/1.5.0/topics/…
Stefan Avey
x,dx,dx,d
@mmk pouvez-vous faire cela pour les densités 2D?
No Lie
4

Tout d'abord, je peux me tromper, mais je pense que votre solution ne fonctionnerait pas dans le cas où il existe plusieurs points où les estimations de densité de noyau (KDE) se croisent. Deuxièmement, bien que le overlappackage ait été créé pour être utilisé avec des données d'horodatage, vous pouvez toujours l'utiliser pour estimer la zone de chevauchement de deux KDE. Vous devez simplement redimensionner vos données afin qu'elles s'étendent de 0 à 2π.
Par exemple :

# 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
la source