LSD_BEAST_Plot <- function(num=c(4,5,6)){
par(mfrow=c(length(num),1))
for (j in 1:length(num)){
load(paste("C:/Users/Oliver/Google Drive/Academic Work/Imperial/O15-12/BEAST/Data Collection/LSD",num[j],"/mu_1e-04_R0_2.4_Tg_50/casetree.RData",sep=""))
dir=paste("C:/Users/Oliver/Google Drive/Academic Work/Imperial/O15-12/BEAST/Data Collection/LSD",num[j],"/",sep="")
#dir.create(path=dir)
network <- paste(network.path,tail(strsplit(dir,"/")[[1]],1),"\\",sep="")
#dir.create(path=network)
network.dir <- paste(network.path,tail(strsplit(dir,"/")[[1]],1),"\\","mu_",mu,"_R0_",R0,"_Tg_",Tg,"\\",sep="")
#dir.create(path=network.dir)
param.dir <- paste(dir,"mu_",mu,"_R0_",R0,"_Tg_",Tg,"/",sep="")
#dir.create(path=param.dir)
res <- matrix(nrow = n.bs, ncol = length(sigma.reps))
m.res <- matrix(nrow = n.bs, ncol = length(sigma.reps))
## import beast
network.root.dir=paste("LSD",num[j],"\\mu_1e-04_R0_2.4_Tg_50\\",sep="")
local.root.dir=paste("LSD",num[j],"/mu_1e-04_R0_2.4_Tg_50/",sep="")
beast.res <- Network_Pull_Tracer(network.root.dir = network.root.dir, local.root.dir = local.root.dir)
root <- unique(beast.res$E_TMRCA)
## set up parameter loops
for (i in 1:length(sigma.reps)){
## create dirs for results
sub.dir <- paste(param.dir,"sigma",sigma.reps[i],"/",sep="")
#dir.create(path=sub.dir)
network.sub.dir <- paste(network.dir,"sigma",sigma.reps[i],"\\",sep="")
#dir.create(path=network.sub.dir)
## write and execute LSD related script
#sh.sub.dir <- paste("\"",gsub(pattern = "/",replacement = "\\",fixed=TRUE,sub.dir),sep="")
#cmd <- paste(lsd.dir," -i ",sh.sub.dir,"bs_tree.txt\" -d ",sh.sub.dir,"bs_date.date\" -r a -n ",n.bs,sep="")
#system
## colate and plotbootstrap error against sigma
lsd.tb <- read.table(file=paste(sub.dir,"bs_tree.txt.result",sep=""),skip = 20,nrows = n.bs)
res[,i] <- lsd.tb[,6]
m.res[,i] <- as.numeric(substr(lsd.tb[,4], 1, nchar(as.character(lsd.tb[,4]))-1))
}
boxplot(res,xlab = "sigma",ylab = "TMRCA Error",ylim=c(-100, 100))
par(new=T)
for (i in 1:length(sigma.reps)){
points(x=c(i,i),y=c(root-beast.res$O_TMRCA_95L[i],root-beast.res$O_TMRCA_95H[i]), type = "l", col = "red",lwd = 2)
points(x=i,y=root-beast.res$O_TMRCA_Med[i], col = "green",lwd = 2)
}
abline(h=case.tree$Time[1],col="blue")
par(new=F)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.