Error bars for skyline plot?

I have made a consensus tree (in newick format) and rebranched it using functions from the ape package.

Does anyone know if it is possible to add error bars to the skyline plot?

Cheers! Input:

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

Branched tree and final input:

> write.tree(new_tree)
[1] "   (13.1:1.244090638,43.1:0.175开发者_开发问答5118382,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)

---

You can use plotCI from package plotrix. Would be fun to append the plot.skyline tree to implement error bars...

Anyway, this is the code I used to produce the below picture. I used a different data set that comes with the ape package. You'll have to take care of the details.

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)