Cieniowanie wykresu gęstości jądra między dwoma punktami.

97

Często używam wykresów gęstości jądra, aby zilustrować rozkłady. Są łatwe i szybkie do utworzenia w R, jak na przykład:

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

Co daje mi ten ładny mały plik PDF:

wprowadź opis obrazu tutaj

Chciałbym zacienić obszar pod plikiem PDF od 75 do 95 centyla. Punkty można łatwo obliczyć za pomocą quantilefunkcji:

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

Ale jak zacienić obszar między q75a q95?

JD Long
źródło
Czy możesz podać przykład cieniowania poza zasięgiem w porównaniu z jego wnętrzem? Dzięki.
Milktrader

Odpowiedzi:

76

Jeśli chodzi o tę polygon()funkcję, zobacz jej stronę pomocy i wydaje mi się, że mieliśmy tutaj podobne pytania.

Musisz znaleźć indeks wartości kwantylowych, aby uzyskać rzeczywiste (x,y)pary.

Edycja: proszę bardzo:

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

Wyjście (dodane przez JDL)

wprowadź opis obrazu tutaj

Dirk Eddelbuettel
źródło
3
Nigdy bym tego nie zrobił, gdybyś nie dostarczył struktury. Dzięki!
JD Long
2
To jedna z tych rzeczy ... które miały miejsce demo(graphics)od zarania dziejów, więc od czasu do czasu można je spotkać. Ten sam pomysł na cieniowanie regresji NBER itp.
Dirk Eddelbuettel
1
ohhhh. WIEDZIAŁEM, że gdzieś to widziałem, ale nie mogłem wyciągnąć z mojego indeksu mentalnego miejsca, w którym to widziałem. Cieszę się, że twój indeks mentalny jest lepszy niż mój.
JD Long
72

Inne rozwiązanie:

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)

Wynik:

tekst alternatywny

Ben Bolker
źródło
22

Rozbudowane rozwiązanie:

Jeśli chcesz przyciemnić oba ogony (skopiuj i wklej kod Dirka) i użyj znanych wartości x:

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

Wynik:

2-ogonowy poli

Milktrader
źródło
Mam plik png i umieściłem go na serwerze freeimagehosting i może się nie ładować, ponieważ ... nie jestem pewien.
Milktrader
Bardzo rozmyty plik. Czy możesz go odtworzyć i przesłać bezpośrednio tutaj SO ma do tego własną usługę serwerową?
Dirk Eddelbuettel
Przepraszam, ale nie widzę, jak przesłać go bezpośrednio do SO.
Milktrader
19

To pytanie wymaga latticeodpowiedzi. Oto bardzo podstawowy, po prostu dostosowujący metodę zastosowaną przez Dirka i innych:

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

wprowadź opis obrazu tutaj

joran
źródło
3

Oto inny ggplot2wariant oparty na funkcji przybliżającej gęstość jądra przy oryginalnych wartościach danych:

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

Korzystanie z oryginalnych danych (zamiast tworzenia nowej ramki danych z wartościami x i y oszacowania gęstości) ma tę zaletę, że działa również na wykresach aspektowych, w których wartości kwantylowe zależą od zmiennej, według której grupowane są dane:

Używany kod

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

Utworzony 2018-07-13 przez pakiet reprex (v0.2.0).

Matt Flor
źródło