Лучший способ построить доверительные полосы вокруг среднего/медианы наблюдаемой выборки с помощью ggplot2

Итак, у меня есть фрейм данных с тремя столбцами, в котором есть испытания, Ind. Переменная, Наблюдение. Что-то вроде:

df1<- data.frame(Trial=rep(1:10,5), Variable=rep(1:5, each=10), Observation=rnorm(1:50))

Я пытаюсь построить 95% conf. Интервал вокруг среднего значения для каждого испытания с использованием довольно неэффективного метода выглядит следующим образом: 

    b<-NULL
    b$mean<- aggregate(Observation~Variable, data=df1,mean)[,2]
    b$sd  <- aggregate(Observation~Variable, data=df1,sd)[,2]
    b$Variable<- df1$Variable
    b$Observation <- df1$Observation 
    b$ucl <- rep(qnorm(.975, mean=b$mean, sd=b$sd), each=10)
    b$lcl <- rep(qnorm(.025, mean=b$mean, sd=b$sd), each=10)
    b<- as.data.frame(b)
    c <- ggplot(b, aes(Variable, Observation))  
    c + geom_point(color="red") + 
    geom_smooth(aes(ymin = lcl, ymax = ucl), data=b, stat="summary", fun.y="mean")

Это неэффективно, так как дублирует значения для ymin, ymax. Я видел методы geom_ribbon, но мне все равно нужно будет дублировать. Однако, если бы я использовал какой-либо вид сглаживания, такой как glm, код намного проще без дублирование. Есть ли лучший способ сделать это?

Ссылки: 1. R построение доверительных полос с помощью ggplot 2. затенение доверительных интервалов вручную с помощью ggplot2 3. http://docs.ggplot2.org/current/geom_smooth.html

2 8
ggplot2 confidence-interval

2 ответа:

С помощью этого метода я получаю тот же результат, что и с вашим методом. Это было вдохновлено документами для ggplot. Опять же, это будет иметь смысл до тех пор, пока каждое значение x имеет несколько точек.

set.seed(1)
df1 <- data.frame(Trial=rep(1:10,5), Variable=rep(1:5, each=10), Observation=rnorm(1:50))    my_ci <- function(x) data.frame(y=mean(x), ymin=mean(x)-2*sd(x), ymax=mean(x)+2*sd(x))

my_ci <- function(x) data.frame(
  y=mean(x), 
  ymin=mean(x) - 2 * sd(x), 
  ymax=mean(x) + 2 * sd(x)
)
ggplot(df1, aes(Variable, Observation)) + geom_point(color="red") +
  stat_summary(fun.data="my_ci", geom="smooth")

Введите описание изображения здесь

7

Пакет ggplot поставляется с оболочками для ряда суммирующих функций в пакете Hmisc, включая

  1. mean_cl_normal, который вычисляет доверительные пределы на основе t-распределения,
  2. mean_cl_boot, который использует метод bootstrap, не предполагающий распределения среднего,
  3. mean_sdl, который использует кратное стандартное отклонение (по умолчанию=2).

Этот последний метод такой же, как и в ответе выше, но являетсяне 95% CL. Доверительные пределы, основанные на T-распределении, задаются следующим образом:

CL = t × s / √n

Где t-соответствующий квантиль t-распределения, а s-стандартное отклонение выборки. Сравните доверительные полосы: 
ggplot(df1, aes(x=Variable, y=Observation)) + 
  stat_summary(fun.data="mean_sdl", geom="line", colour="blue")+
  stat_summary(fun.data="mean_sdl", mult=2, geom="errorbar", 
               width=0.1, linetype=2, colour="blue")+
  geom_point(color="red") +
  labs(title=expression(paste(bar(x)," \u00B1 ","2 * sd")))


ggplot(df1, aes(x=Variable, y=Observation)) + 
  geom_point(color="red") +
  stat_summary(fun.data="mean_cl_normal", geom="line", colour="blue")+
  stat_summary(fun.data="mean_cl_normal", conf.int=0.95, geom="errorbar", 
               width=0.1, linetype=2, colour="blue")+
  stat_summary(fun.data="mean_cl_normal", geom="point", size=3, 
               shape=1, colour="blue")+
  labs(title=expression(paste(bar(x)," \u00B1 ","t * sd / sqrt(n)")))

Наконец, поворот этого последнего графика с помощью coord_flip() создает что-то очень близкое к Forest Plot, что является стандартным методом для суммирования данных, подобных вашему. 
ggplot(df1, aes(x=Variable, y=Observation)) + 
  geom_point(color="red") +
  stat_summary(fun.data="mean_cl_normal", conf.int=0.95, geom="errorbar", 
               width=0.2, colour="blue")+
  stat_summary(fun.data="mean_cl_normal", geom="point", size=3, 
               shape=1, colour="blue")+
  geom_hline(aes(yintercept=mean(Observation)), linetype=2)+
  labs(title="Forest Plot")+
  coord_flip()

4