Ggplot2 time series envelope shading

Question

I draw the results of 50 - 100 experiments. Each experiment leads to a time series. I can plot the spaghetti chart of all time series, but what I would like to have is a density map for the time series. (something similar to gray shading in the bottom panel in this illustration: http://www.ipcc.ch/graphics/ar4-wg1/jpg/fig-6-14.jpg )

I can sort this using 2d binning or binhex, but the result may be more beautiful (see example below).

Here is the code that reproduces the plume graph for the layout data (uses ggplot2 and 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") 

Does anyone know a good way to get a shaded envelope with gradients. I also tried geom_ribbon, but this did not give any indication of density changes along the plume. binhex does this, but not with aesthetically pleasing results.

Answer 1

Calculation of quantiles:

 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 

Graphic tapes:

 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) 

This is for two quantile intervals (0-100) and (25-75). You will need more arguments for quantile and more tape layers for more quantiles, and you will also need to adjust the colors.

Answer 2

Based on the idea of ​​Spacedman, I found a way to add extra intervals automatically: first I calculated the quantiles for each step , grouped them into pairs of symmetric values, and then used geom_ribbon in the right order ...

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

Answer 3

Thanks Erwan and Spacedman.

Avoiding "tidyr" ("dplyr" and "magrittr"), my version of Erwans answer becomes

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