スカイラインプロットのエラーバーですか？

Question

私はコンセンサスツリー（newick形式）を作成し、apeパッケージの関数を使用して再バインドしました。

スカイラインのプロットにエラーバーを追加することができますか？

乾杯！ 入力：

>write.tree(ccc) 
[1] "(13.1,43.1,28.1,21.1,20.1,4.1,37.1,23.1,29.1,15.1,36.1,40.1,26.1,16.1,33.1,10.1,18.1,6.1,11.1,34.1,2.1,38.1,35.1,8.1,42.1,22.1,5.1,27.1,39.1,17.1,30.1,31.1,41.1,7.1,25.1,3.1,12.1,14.1,1.1,24.1,9.1,32.1,19.1);" 

分岐ツリーと最終入力：

> write.tree(new_tree) 
[1] " (13.1:1.244090638,43.1:0.1755118382,28.1:4.442071925,21.1:4.866247957,20.1:3.863557777,4.1:0.08206463186,37.1:1.180065627,23.1:4.778251834,29.1:4.65377018,15.1:3.726316827,36.1:3.799621998,40.1:1.707243007,26.1:3.655532246,16.1:2.436848865,33.1:1.315581813,10.1:0.5973169219,18.1:2.561718575,6.1:1.409437305,11.1:4.000977813,34.1:2.769218051,2.1:1.143412833,38.1:2.636477078,35.1:3.435094416,8.1:0.3714309109,42.1:4.429772993,22.1:0.2805716533,5.1:2.991251546,27.1:3.269689628,39.1:1.192241808,17.1:1.866855259,30.1:1.363407132,31.1:3.150236446,41.1:4.079005712,7.1:2.695786237,25.1:2.867184281,3.1:3.351797838,12.1:1.958669939,14.1:3.119812493,1.1:4.864444738,24.1:1.734493783,9.1:3.283970281,32.1:3.055829309,19.1:0.8193949074);" 

require(ape) 
new_tree<-compute.brlen(ccc,runif,min=0,max=5) 
sk1 <- skyline(new_tree) 
sk2 <- skyline(new_tree, 0.0119) 
plot(sk1, show.years=TRUE, subst.rate=0.0023, present.year = 2011, col="red") 
lines(sk2, show.years=TRUE, subst.rate=0.0023, present.year = 2011) 
legend(.15,500, c("classic", "generalized"), col="red",lty=1) 

Answer 1

あなたがパッケージplotrixからplotCIを使用することができます。 plot.skylineツリーを追加してエラーバーを実装するのは楽しいですか？

とにかく、これは私が下の画像を生成するために使用したコードです。私はapeパッケージに付属の別のデータセットを使用しました。あなたは詳細の世話をする必要があります。

library(ape) 
library(plotrix) 
data(bird.orders) 

new_tree <- compute.brlen(bird.orders, runif, min = 0, max = 5) 
sk1 <- skyline(new_tree) 

subst.rate <- 0.0023 
year <- 2011 
m <- sk1$population.size 
lm <- length(m) 

plot(sk1, show.years=TRUE, subst.rate = subst.rate, present.year = year, col="red") 

plotCI(x = (-c(0, sk1$time))/subst.rate + year, y = c(m, m[lm]), 
     uiw = runif(length(sk1$time)+1, min = 1, max = 30), # some random deviation 
     #liw = runif(length(sk1$time)+1, min = 1, max = 30), # this one can be missing, see ?plotCI 
     add = TRUE, 
     pt.bg = NA 
) 

legend(.15,500, c("classic", "generalized"), col="red",lty=1)