ggplot2時系列の陰影エンベロープ

Question

私は50〜100回の実験結果をプロットしています。 各実験では時系列が得られます。 私はすべての時系列のスパゲッティプロットをプロットすることができますが、 は、時系列プルームの密度マップの一種です。 （この図の下のパネル の灰色の影のようなもの：http://www.ipcc.ch/graphics/ar4-wg1/jpg/fig-6-14.jpg）

私は2Dビニングまたは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はそれを行いますが、審美的に満足な結果は得られません。

Answer 1

計算変位値：

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) 

これは、2つの変位値の間隔、（0-100）および（25〜75）のためのものです。あなたはquantileとより多くのquantilesのためのより多くのリボン層にもっと多くのargsを必要とし、また色を調整する必要があります。

Spacedmanの考えに基づいて

Answer 2

、私は自動的な方法でより多くの間隔を追加する方法を発見した：私は最初対称値の対により、各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") 

Answer 3

Thanks Erwan and Spacedman。

Erwansの答えの「tidyr」（「dplyr」と「magrittr」）私のバージョンは、私はあなたが、各時点で、あなたのシリーズの分位数を計算したい疑う

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