Intersection entre les tracés de densité de plusieurs groupes

Question

J'utilise ggplot/easyGgplot2 pour créer des graphiques de densité de deux groupes. J'aimerais avoir une mesure ou une indication de la distance entre les deux courbes. Je pourrais même utiliser n'importe quelle autre solution sans les courbes, tant que cela me permet d'avoir une mesure de quels groupes sont plus distincts (de plusieurs groupes de données différents).

Existe-t-il un moyen simple de faire cela dans R?

Par exemple en utilisant cet échantillon, qui génère cette parcelle

Comment puis-je estimer le pourcentage de la superficie qui est commune à la fois?

ggplot2.density(data=weight, xName='weight', groupName='sex', 
    legendPosition="top", 
    alpha=0.5, fillGroupDensity=TRUE) 

Answer 1

J'aime la réponse précédente, mais cela peut être un peu plus intuitive, aussi je me suis assuré d'utiliser une bande passante commune:

library ("caTools") 

# Extract common bandwidth 
Bw <- (density (iris$Petal.Width))$bw 

# Get iris data 
Sample <- with (iris, split (Petal.Width, Species))[ 2:3 ] 

# Estimate kernel densities using common bandwidth 
Densities <- lapply (Sample, density, 
         bw = bw, 
         n = 512, 
         from = -1, 
         to = 3) 

# Plot 
plot(Densities [[ 1 ]], xlim = c (-1, 3), 
     col = "steelblue", 
     main = "") 
lines (Densities [[ 2 ]], col = "orange") 

# Overlap 
X <- Densities [[ 1 ]]$x 
Y1 <- Densities [[ 1 ]]$y 
Y2 <- Densities [[ 2 ]]$y 

Overlap <- pmin (Y1, Y2) 
polygon (c (X, X [ 1 ]), c (Overlap, Overlap [ 1 ]), 
    lwd = 2, col = "hotpink", border = "n", density = 20) 

# Integrate 
Total <- trapz (X, Y1) + trapz (X, Y2) 
(Surface <- trapz (X, Overlap)/Total) 
SText <- paste (sprintf ("%.3f", 100*Surface), "%") 
text (X [ which.max (Overlap)], 1.2 * max (Overlap), SText) 

Answer 2

D'abord, faites quelques données à utiliser. Ici, nous allons regarder la largeur de pétale de deux espèces végétales de l'ensemble de données intégré iris.

## Some sample data from iris 
dat <- droplevels(with(iris, iris[Species %in% c("versicolor", "virginica"), ])) 

## make a similar graph 
library(ggplot2) 
ggplot(dat, aes(Petal.Width, fill=Species)) + 
    geom_density(alpha=0.5) 

Pour trouver la zone de l'intersection, vous pouvez utiliser approxfun pour approcher la fonction décrivant le chevauchement. Ensuite, intégrez-le pour obtenir la zone. Puisque ce sont des courbes de densité, leur aire est 1 (ish), donc l'intégrale sera le pourcentage de chevauchement.

## Get density curves for each species 
ps <- lapply(split(dat, dat$Species), function(x) { 
    dens <- density(x$Petal.Width) 
    data.frame(x=dens$x, y=dens$y) 
}) 

## Approximate the functions and find intersection 
fs <- sapply(ps, function(x) approxfun(x$x, x$y, yleft=0, yright=0)) 
f <- function(x) fs[[1]](x) - fs[[2]](x) # function to minimize (difference b/w curves) 
meet <- uniroot(f, interval=c(1, 2))$root # intersection of the two curves 

## Find overlapping x, y values 
ps1 <- is.na(cut(ps[[1]]$x, c(-Inf, meet))) 
ps2 <- is.na(cut(ps[[2]]$x, c(Inf, meet))) 
shared <- rbind(ps[[1]][ps1,], ps[[2]][ps2,]) 

## Approximate function of intersection 
f <- with(shared, approxfun(x, y, yleft=0, yright=0)) 

## have a look 
xs <- seq(0, 3, len=1000) 
plot(xs, f(xs), type="l", col="blue", ylim=c(0, 2)) 

points(ps[[1]], col="red", type="l", lty=2, lwd=2) 
points(ps[[2]], col="blue", type="l", lty=2, lwd=2) 

polygon(c(xs, rev(xs)), y=c(f(xs), rep(0, length(xs))), col="orange", density=40) 

## Integrate it to get the value 
integrate(f, lower=0, upper=3)$value 
# [1] 0.1548127