2013-11-29 3 views
6

Я пытаюсь нарисовать ECDF некоторых данных с «доверительным интервалом», представленным через затененную область, используя ggplot2. У меня возникли проблемы с объединением geom_ribbon() с stat_ecdf() для достижения эффекта, которым я пользуюсь.Как совместить stat_ecdf с geom_ribbon?

Рассмотрим следующий пример данных:

set.seed(1) 
dat <- data.frame(variable = rlnorm(100) + 2) 
dat <- transform(dat, lower = variable - 2, upper = variable + 2) 

> head(dat) 
    variable  lower upper 
1 2.534484 0.5344838 4.534484 
2 3.201587 1.2015872 5.201587 
3 2.433602 0.4336018 4.433602 
4 6.929713 4.9297132 8.929713 
5 3.390284 1.3902836 5.390284 
6 2.440225 0.4402254 4.440225 

Я могу производить ECDF из variable использованием

library("ggplot2") 
ggplot(dat, aes(x = variable)) + 
    geom_step(stat = "ecdf") 

Однако я не могу использовать lower и upper как в ymin и ymax эстетики geom_ribbon(), чтобы наложить доверительный интервал на график как другой слой. Я пробовал:

ggplot(dat, aes(x = variable)) + 
    geom_ribbon(aes(ymin = lower, ymax = upper), stat = "ecdf") + 
    geom_step(stat = "ecdf") 

, но при этом возникает следующую ошибку

Error: geom_ribbon requires the following missing aesthetics: ymin, ymax 

Есть ли способ, чтобы задобрить geom_ribbon() в работе с stat_ecdf() производить тенистый доверительный интервал? Или может ли кто-нибудь предложить альтернативное средство добавления затененного многоугольника, определенного lower и upper в качестве слоя на участок ECDF?

ответ

3

Попробуйте это (немного выстрел в темноте):

ggplot(dat, aes(x = variable)) + 
    geom_ribbon(aes(x = variable,ymin = ..y..-2,ymax = ..y..+2), stat = "ecdf",alpha=0.2) + 
    geom_step(stat = "ecdf") 

Ok, так что это не то же самое, что вы пытаетесь сделать, но он должен объяснить, что происходит. stat возвращает фрейм данных только с исходным x и вычисленным y, поэтому я думаю, что это все, с чем вам нужно работать. то есть stat_ecdf вычисляет только кумулятивную функцию распределения для одного x за раз.

Единственная вещь, которую я могу думать очевидное, вычислением нижнее и верхнее отдельно, что-то вроде этого:

l <- ecdf(dat$lower) 
u <- ecdf(dat$upper) 
v <- ecdf(dat$variable) 
dat$lower1 <- l(dat$variable) 
dat$upper1 <- u(dat$variable) 
dat$variable1 <- v(dat$variable) 

ggplot(dat,aes(x = variable)) + 
    geom_step(aes(y = variable1)) + 
    geom_ribbon(aes(ymin = upper1,ymax = lower1),alpha = 0.2) 
+0

Спасибо, Джоран. Не могли бы вы расширить свое последнее предложение? Не уверен, что я полностью следую этому, но насколько я могу судить по вашему ответу, я не могу сделать это через 'stat_ecdf', если' lower' и 'upper' уже существуют? +/- 2 бит - это просто фиктивные данные; информация CI, которую я имею, является результатом последующего моделирования производной статистики, вычисленной по установленной модели. –

+1

@GavinSimpson Да, я думаю, что это невозможно сразу в ggplot (хотя это будет хорошая возможность добавить, я думаю). Все, что я имел в виду с последним, было то, что вам, возможно, придется вычислять все значения ECDF вручную, а затем строить их. – joran

+0

Спасибо, я понимаю, что вы имеете в виду, вычислите совокупную долю напрямую. Я дам это. +1 –

2

Не уверен, как именно вы хотите, чтобы отразить CI, но ggplot_build() позволяет получить сгенерированные данные обратно из сюжета, вы можете переопределить то, что вам нравится.

Эта диаграмма показывает:

  • красный = оригинальная лента
  • синий = принимает исходные векторы CI и применяет к кривой ecdf
  • зеленый = вычисляет ecdf верхнего и нижнего ряда и участков

enter image description here

g<-ggplot(dat, aes(x = variable)) + 
     geom_step(stat = "ecdf") + 
     geom_ribbon(aes(ymin = lower, ymax = upper), alpha=0.5, fill="red") 

    inside<-ggplot_build(g) 
    matched<-merge(inside$data[[1]],data.frame(x=dat$variable,dat$lower,dat$upper),by=("x")) 

    g + 
     geom_ribbon(data=matched, aes(x = x, 
             ymin = y + dat.upper-x, 
             ymax = y - x + dat.lower), 
        alpha=0.5, fill="blue") + 
     geom_ribbon(data=matched, aes(x = x, 
             ymin = ecdf(dat.lower)(x), 
             ymax = ecdf(dat.upper)(x)), 
        alpha=0.5, fill="green") 
+0

Thanks Troy; ваша последняя идея, как только я понял, что показал сюжет, похож на идею @ joran, а именно, что можно вычислить соответствующие данные 'y' для нижнего и верхнего CI, используя' ecdf() '. Зеленая лента - это то, что я хочу изобразить. –

Смежные вопросы