Ombrage d'un graphique de densité de noyau entre deux points.

Question

J'utilise fréquemment les diagrammes de densité de noyau pour illustrer les distributions. Il est facile et rapide de les créer dans R comme suit :

set.seed(1)
draws <- rnorm(100)^2
dens <- density(draws)
plot(dens)
#or in one line like this: plot(density(rnorm(100)^2))

Ce qui me donne ce joli petit PDF :

J'aimerais ombrer la zone sous le PDF du 75ème au 95ème percentiles. Il est facile de calculer les points à l'aide de la fonction quantile fonction :

q75 <- quantile(draws, .75)
q95 <- quantile(draws, .95)

Mais comment ombrager la zone entre q75 y q95 ?

Answer 1

Avec le polygon() voir sa page d'aide et je crois que nous avons eu des questions similaires ici aussi.

Vous devez trouver l'indice des valeurs du quantile pour obtenir les valeurs réelles. (x,y) paires.

Edit : Voilà :

x1 <- min(which(dens$x >= q75))  
x2 <- max(which(dens$x <  q95))
with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray"))

Sortie (ajouté par JDL)

Answer 2

Une autre solution :

dd <- with(dens,data.frame(x,y))

library(ggplot2)

qplot(x,y,data=dd,geom="line")+
  geom_ribbon(data=subset(dd,x>q75 & x<q95),aes(ymax=y),ymin=0,
              fill="red",colour=NA,alpha=0.5)

Résultat :

Answer 3

Une solution élargie :

Si vous vouliez ombrer les deux queues (copier-coller du code de Dirk) et utiliser des valeurs x connues :

set.seed(1)
draws <- rnorm(100)^2
dens <- density(draws)
plot(dens)

q2     <- 2
q65    <- 6.5
qn08   <- -0.8
qn02   <- -0.2

x1 <- min(which(dens$x >= q2))  
x2 <- max(which(dens$x <  q65))
x3 <- min(which(dens$x >= qn08))  
x4 <- max(which(dens$x <  qn02))

with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray"))
with(dens, polygon(x=c(x[c(x3,x3:x4,x4)]), y= c(0, y[x3:x4], 0), col="gray"))

Résultat :

Answer 4

Cette question nécessite une lattice réponse. En voici une très basique, qui ne fait qu'adapter la méthode employée par Dirk et d'autres :

#Set up the data
set.seed(1)
draws <- rnorm(100)^2
dens <- density(draws)

#Put in a simple data frame   
d <- data.frame(x = dens$x, y = dens$y)

#Define a custom panel function;
# Options like color don't need to be hard coded    
shadePanel <- function(x,y,shadeLims){
    panel.lines(x,y)
    m1 <- min(which(x >= shadeLims[1]))
    m2 <- max(which(x <= shadeLims[2]))
    tmp <- data.frame(x1 = x[c(m1,m1:m2,m2)], y1 = c(0,y[m1:m2],0))
    panel.polygon(tmp$x1,tmp$y1,col = "blue")
}

#Plot
xyplot(y~x,data = d, panel = shadePanel, shadeLims = c(1,3))

Answer 5

En voici un autre ggplot2 variante basée sur une fonction qui approxime la densité du noyau aux valeurs des données originales :

approxdens <- function(x) {
    dens <- density(x)
    f <- with(dens, approxfun(x, y))
    f(x)
}

L'utilisation des données d'origine (plutôt que la production d'un nouveau cadre de données avec les valeurs x et y de l'estimation de la densité) présente l'avantage de fonctionner également dans des graphiques à facettes où les valeurs des quantiles dépendent de la variable par laquelle les données sont regroupées :

Code utilisé

library(tidyverse)
library(RColorBrewer)

# dummy data
set.seed(1)
n <- 1e2
dt <- tibble(value = rnorm(n)^2)

# function that approximates the density at the provided values
approxdens <- function(x) {
    dens <- density(x)
    f <- with(dens, approxfun(x, y))
    f(x)
}

probs <- c(0.75, 0.95)

dt <- dt %>%
    mutate(dy = approxdens(value),                         # calculate density
           p = percent_rank(value),                        # percentile rank 
           pcat = as.factor(cut(p, breaks = probs,         # percentile category based on probs
                                include.lowest = TRUE)))

ggplot(dt, aes(value, dy)) +
    geom_ribbon(aes(ymin = 0, ymax = dy, fill = pcat)) +
    geom_line() +
    scale_fill_brewer(guide = "none") +
    theme_bw()

# dummy data with 2 groups
dt2 <- tibble(category = c(rep("A", n), rep("B", n)),
              value = c(rnorm(n)^2, rnorm(n, mean = 2)))

dt2 <- dt2 %>%
    group_by(category) %>% 
    mutate(dy = approxdens(value),    
           p = percent_rank(value),
           pcat = as.factor(cut(p, breaks = probs,
                                include.lowest = TRUE)))

# faceted plot
ggplot(dt2, aes(value, dy)) +
    geom_ribbon(aes(ymin = 0, ymax = dy, fill = pcat)) +
    geom_line() +
    facet_wrap(~ category, nrow = 2, scales = "fixed") +
    scale_fill_brewer(guide = "none") +
    theme_bw()

Créé le 2018-07-13 par le paquet reprex (v0.2.0).