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

J'utilise fréquemment des graphiques de densité de noyau pour illustrer les distributions. Ceux-ci sont faciles et rapides à créer dans R comme ceci:

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:

Je voudrais ombrer la zone sous le PDF du 75e au 95e centile. Il est facile de calculer les points en utilisant la quantilefonction:

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

Mais comment ombrer la zone entre q75et q95?

Avec la polygon()fonction, consultez sa page d'aide et je pense que nous avons eu des questions similaires ici aussi.

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

Edit: Ici vous allez:

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ée par JDL)

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:

Une solution étendue:

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:

Cette question a besoin d'une latticeréponse. En voici une très basique, en adaptant simplement 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))

Voici une autre ggplot2variante basée sur une fonction qui se rapproche de la densité du noyau aux valeurs de données d'origine:

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

L'utilisation des données d'origine (plutôt que de produire une nouvelle base de données avec les valeurs x et y de l'estimation de densité) a l'avantage de travailler également dans des graphiques à facettes où les valeurs de quantile 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 13/07/2018 par le package reprex (v0.2.0).