首页 文章

ggplot2时间序列的阴影信封

提问于
浏览
9

我正在绘制50-100个实验的结果 . 每个实验都会产生一个时间序列 . 我可以绘制所有时间序列的意大利面条图,但我想要的是时间序列羽流的密度图 . (类似于此图中下方面板中的灰色阴影:http://www.ipcc.ch/graphics/ar4-wg1/jpg/fig-6-14.jpg

enter image description here

我可以通过2d binning或binhex“做”这样做,但结果可能更漂亮(参见下面的示例) .

这是一个代码,用于再现模拟数据的羽流图(使用ggplot2和reshape2) .

# mock data: random walk plus a sinus curve.
# two envelopes for added contrast.
tt=10*sin(c(1:100)/(3*pi))
rr=apply(matrix(rnorm(5000),100,50),2,cumsum) +tt
rr2=apply(matrix(rnorm(5000),100,50),2,cumsum)/1.5 +tt

# stuff data into a dataframe and melt it.
df=data.frame(c(1:100),cbind(rr,rr2) )
names(df)=c("step",paste("ser",c(1:100),sep=""))
dfm=melt(df,id.vars = 1)

# ensemble average
ensemble_av=data.frame(step=df[,1],ensav=apply(df[,-1],1,mean))
ensemble_av$variable=as.factor("Mean")


ggplot(dfm,aes(step,value,group=variable))+
  stat_binhex(alpha=0.2) + geom_line(alpha=0.2) +
  geom_line(data=ensemble_av,aes(step,ensav,size=2))+
  theme(legend.position="none")

有没有人知道一个很好的方法来获得带渐变的阴影信封 . 我也试过geom_ribbon,但没有给出任何沿羽流密度变化的迹象 . binhex做到了这一点,但没有美学上令人愉悦的结果 .

3 回答

  • 8

    计算分位数:

    qs = data.frame(
       do.call(
           rbind,
           tapply(
             dfm$value, dfm$step, function(i){quantile(i)})),
       t=1:100)
    
    head(qs)
             X0.      X25.      X50.     X75.     X100. t
    1 -0.8514179 0.4197579 0.7681517 1.396382  2.883903 1
    2 -0.6506662 1.2019163 1.6889073 2.480807  5.614209 2
    3 -0.3182652 2.0480082 2.6206045 4.205954  6.485394 3
    4 -0.1357976 2.8956990 4.2082762 5.138747  8.860838 4
    5  0.8988975 3.5289219 5.0621513 6.075937 10.253379 5
    6  2.0027973 4.5398120 5.9713921 7.015491 11.494183 6
    

    情节色带:

    ggplot() + 
     geom_ribbon(data=qs, aes(x=t, ymin=X0., ymax=X100.),fill="gray30", alpha=0.2) +
     geom_ribbon(data=qs, aes(x=t, ymin=X25., ymax=X75.),fill="gray30", alpha=0.2)
    

    quantile intervals

    这是针对两个分位数间隔,(0-100)和(25-75) . 对于更多的分位数,你需要更多的args到 quantile 和更多的带状层,并且还需要调整颜色 .

  • 0

    根据Spacedman的想法,我找到了一种以自动方式添加更多间隔的方法:我首先计算每个 step 的分位数,用成对的对称值对它们进行分组,然后按正确的顺序使用 geom_ribbon ...

    library(tidyr)
    library(dplyr)
    condquant <- dfm %>% group_by(step) %>%
      do(quant = quantile(.$value, probs = seq(0,1,.05)), probs = seq(0,1,.05)) %>%
      unnest() %>%
      mutate(delta = 2*round(abs(.5-probs)*100)) %>% 
      group_by(step, delta) %>%
      summarize(quantmin = min(quant), quantmax= max(quant))
    
    ggplot() +
      geom_ribbon(data = condquant, aes(x = step, ymin = quantmin, ymax = quantmax,
                                        group = reorder(delta, -delta), fill = as.numeric(delta)),
                  alpha = .5) +
      scale_fill_gradient(low = "grey10", high = "grey95") + 
      geom_line(data = dfm, aes(x = step, y = value, group=variable), alpha=0.2) +
      geom_line(data=ensemble_av,aes(step,ensav),size=2)+
      theme(legend.position="none")
    
  • 1

    谢谢Erwan和Spacedman .

    避免'tidyr'('dplyr'和'magrittr')我的Erwans回答版本变为

    probs=c(0:10)/10  # use fewer quantiles than Erwan
    arr=t(apply(df[,-1],1,quantile,prob=probs))
    dfq=data.frame(step=df[,1],arr)
    names(dfq)=c("step",colnames(arr))
    dfqm=melt(dfq,id.vars=c(1))
    # add inter-quantile (per) range as delta 
     dfqm$delta=dfqm$variable
     levels(dfqm$delta)=abs(probs-rev(probs))*100
    
    
    dfplot=ddply(dfqm,.(step,delta),summarize,
      quantmin=min(value),
      quantmax=max(value) )
    
    ggplot() +
      geom_ribbon(data = dfplot, aes(x = step, ymin = quantmin, 
                                   ymax =quantmax,group=rev(delta),
                                   fill = as.numeric(delta)),
                 alpha = .5) +
      scale_fill_gradient(low = "grey25", high = "grey75") +
      geom_line(data=ensemble_av,aes(step,ensav),size=2) + 
      theme(legend.position="none")
    

    Result of code

相关问题