2015-10-12 12 views
5

Mam wykres punktowy, chcę wiedzieć, w jaki sposób mogę znaleźć geny powyżej i poniżej linii przedziału ufności?Znajdź punkty powyżej i poniżej przedziału ufności podczas korzystania z geom_stat/geom_smooth w ggplot2

enter image description here


EDIT: Powtarzalne przykład:

library(ggplot2) 
#dummy data 
df <- mtcars[,c("mpg","cyl")] 

#plot 
ggplot(df,aes(mpg,cyl)) + 
    geom_point() + 
    geom_smooth() 

enter image description here

+7

można uruchomić poprzez włączenie kodu i danych. – nrussell

+0

'ident (x, y ...)' ale część twoich danych jest potrzebna – Mateusz1981

+0

Należy zauważyć, że linie przedziału ufności są przedziałem ufności dla średniej danych, a nie dla samych danych. A ponieważ masz tak dużo danych, spodziewam się, że większość wartości będzie poza przedziałem. – bramtayl

Odpowiedz

7

Musiałam wziąć głęboki nur do github repo ale w końcu dostał. Aby to zrobić, musisz wiedzieć, jak działa stat_smooth. W tym konkretnym przypadku funkcja loess nazywa wykonać wygładzanie (różne funkcje wygładzania może być wykonana przy użyciu tego samego procesu, jak poniżej):

Więc korzystając loess przy tej okazji chcielibyśmy zrobić:

#data 
df <- mtcars[,c("mpg","cyl"), with=FALSE] 
#run loess model 
cars.lo <- loess(cyl ~ mpg, df) 

Następnie musiałem przeczytać this, aby zobaczyć, jak przewidywania są dokonywane wewnętrznie w stat_smooth. Widocznie Hadley wykorzystuje funkcję predictdf (co nie jest eksportowany do przestrzeni nazw) w następujący sposób dla naszego przypadku:

predictdf.loess <- function(model, xseq, se, level) { 
    pred <- stats::predict(model, newdata = data.frame(x = xseq), se = se) 

    if (se) { 
    y = pred$fit 
    ci <- pred$se.fit * stats::qt(level/2 + .5, pred$df) 
    ymin = y - ci 
    ymax = y + ci 
    data.frame(x = xseq, y, ymin, ymax, se = pred$se.fit) 
    } else { 
    data.frame(x = xseq, y = as.vector(pred)) 
    } 
} 

Po przeczytaniu powyższego udało mi się stworzyć własną data.frame z przewidywaniami, używając:

#get the predictions i.e. the fit and se.fit vectors 
pred <- predict(cars.lo, se=TRUE) 
#create a data.frame from those 
df2 <- data.frame(mpg=df$mpg, fit=pred$fit, se.fit=pred$se.fit * qt(0.95/2 + .5, pred$df)) 

Patrząc na predictdf.loess widzimy, że górna granica przedziału ufności jest tworzony jako pred$fit + pred$se.fit * qt(0.95/2 + .5, pred$df) a dolna granica jako pred$fit - pred$se.fit * qt(0.95/2 + .5, pred$df).

Używanie tych możemy stworzyć flagę punktów powyżej lub poniżej tych granic:

#make the flag 
outerpoints <- +(df$cyl > df2$fit + df2$se.fit | df$cyl < df2$fit - df2$se.fit) 
#add flag to original data frame 
df$outer <- outerpoints 

Kolumna df$outer jest chyba to, co PO szuka (to przyjmuje wartość 1, jeśli jest poza granice lub 0 inaczej), ale tylko ze względu na to spiskuję go poniżej.

Zauważ, że powyższa funkcja + służy tylko do przekształcenia flagi logicznej w numeryczną.

Teraz jeśli wykreślić jako to:

ggplot(df,aes(mpg,cyl)) + 
    geom_point(aes(colour=factor(outer))) + 
    geom_smooth() 

Możemy faktycznie patrz punkty wewnątrz i na zewnątrz przedziału ufności.

wyjściowa:

enter image description here

PS:Dla każdego, kto jest zainteresowany w górnych i dolnych granic, są one tworzone tak (spekulacji: chociaż zacienione obszary są prawdopodobnie utworzone z geom_ribbon - lub coś podobnego - co czyni je bardziej okrągłe i ładny):

#upper boundary 
ggplot(df,aes(mpg,cyl)) + 
    geom_point(aes(colour=factor(outer))) + 
    geom_smooth() + 
    geom_line(data=df2, aes(mpg , fit + se.fit , group=1), colour='red') 

#lower boundary 
ggplot(df,aes(mpg,cyl)) + 
    geom_point(aes(colour=factor(outer))) + 
    geom_smooth() + 
    geom_line(data=df2, aes(mpg , fit - se.fit , group=1), colour='red') 
+1

miło, miał zamiar opublikować porównywalną odpowiedź ;-) – Jaap

+0

Dzięki @Jaap :). Przepraszam za to, wiem jak to jest z doświadczenia :). Publikuj to, jeśli uważasz, że dodaje dodatkowe informacje. – LyzandeR

+1

nie ma potrzeby, nie mam nic do poprawienia na twoją odpowiedź :-) (oprócz drobnych zmian) – Jaap

8

rozwiązanie to wykorzystuje dysku ggplot2 pracy robi dla Ciebie:

library(sp) 

# we have to build the plot first so ggplot can do the calculations 
ggplot(df,aes(mpg,cyl)) + 
    geom_point() + 
    geom_smooth() -> gg 

# do the calculations 
gb <- ggplot_build(gg) 

# get the CI data 
p <- gb$data[[2]] 

# make a polygon out of it 
poly <- data.frame(
    x=c(p$x[1], p$x, p$x[length(p$x)], rev(p$x)), 
    y=c(p$ymax[1], p$ymin, p$ymax[length(p$x)], rev(p$ymax)) 
) 

# test for original values in said polygon and add that to orig data 
# so we can color by it 
df$in_ci <- point.in.polygon(df$mpg, df$cyl, poly$x, poly$y) 

# re-do the plot with the new data 
ggplot(df,aes(mpg,cyl)) + 
    geom_point(aes(color=factor(in_ci))) + 
    geom_smooth() 

enter image description here

to wymaga trochę szczypanie (czyli tej ostatniej kwestii uzyskiwanie wartości 2) ale jestem ograniczony na czas. Zauważ, że wartości zwracane point.in.polygon są:

  • 0: punkt jest ściśle zewnętrzne Pol
  • 1: punkt jest ściśle wnętrze Pol
  • 2: punkt leży na względnej wewnętrznej krawędzi pol
  • 3: punkt jest wierzchołkiem pol

więc powinno być łatwe wystarczy zmienić kod do TRUE/FALSE czy wartość jest 0 czy nie.

6

Korzystając z niezłego rozwiązania @ hrbrmstr, można to zrobić, po prostu przekazując sekwencję wartości x do geom_smooth określając, gdzie powinny zostać obliczone granice błędów, i sprawić, by były równe wartościom x punktów. Następnie zobaczysz, czy wartości y są w tym zakresie.

library(ggplot2) 

## dummy data 
df <- mtcars[,c("mpg","cyl")] 

ggplot(df, aes(mpg, cyl)) + 
    geom_smooth(params=list(xseq=df$mpg)) -> gg 

## Find the points within bounds 
bounds <- ggplot_build(gg)[[1]][[1]] 
df$inside <- with(df, bounds$ymin < cyl & bounds$ymax > cyl) 

## Add the points 
gg + geom_point(data=df, aes(color=inside)) + theme_bw() 

enter image description here

Powiązane problemy