#analysis with ntaxa held constant, generated to have ~100 taxa at every r #three levels of sampling: complete (rInf) and incomplete sampling with r=0.01 and 0.1 per Ltu #with 1000 datasets per mode and 100 taxa each #will do #PropRes - resolvable proportion of clades #measured now as Colless's Consensus Fork Index (Colless, 1980) #MaxNodeSize - max polytomy size #PolyDepth - how evenly is polytomy depth distributed? #p-value of chisq test of a contingency table of polytomies to dichotomous nodes #across the four depth-quartiles of nodes #non-significant values reflect equal distribution of unresolvable clades across depths #for use on geomips library(paleotree) nodeInfo<-function(tree){ #depth cannot be based on number of descendents or will most 'shallow' tips cannot be polytomies #instead make depth based on number of nodes away from root require(ape) nDesc<-sapply(Ntip(tree)+1:Nnode(tree),function(x) sum(x==tree$edge[,1])) tree$edge.length<-rep(1,Nedge(tree)) depths<-dist.nodes(tree)[Ntip(tree)+1,-(1:Ntip(tree))] return(cbind(nDesc=nDesc,Depths=depths)) } polyEvenness<-function(nodeData){ #uses four quantiles by default sorter<-sapply(nodeData[,2],function(x) sum(x>=quantile(nodeData[,2])[1:4])) mat<-t(sapply(1:4,function(x) c(sum(nodeData[sorter==x,1]>2),sum(nodeData[sorter==x,1]==2)))) return(chisq.test(mat)$p.value) } #set number of runs nruns<-100 output_propres_r0.01<-"r0.01_propRes_output.txt" output_propres_r0.1<-"r0.1_propRes_output.txt" output_propres_rInf<-"rInf_propRes_output.txt" output_maxnode_r0.01<-"r0.01_maxNodeSize_output.txt" output_maxnode_r0.1<-"r0.1_maxNodeSize_output.txt" output_maxnode_rInf<-"rInf_maxNodeSize_output.txt" output_polydepth_r0.01<-"r0.01_polyDepth_output.txt" output_polydepth_r0.1<-"r0.1_polyDepth_output.txt" output_polydepth_rInf<-"rInf_polyDepth_output.txt" output_ntaxa_r0.01<-"r0.01_ntaxa_output.txt" output_ntaxa_r0.1<-"r0.1_ntaxa_output.txt" output_ntaxa_rInf<-"rInf_ntaxa_output.txt" ########################################### set.seed(444) propRes_r0.01<-replicate(13,numeric(),simplify=FALSE) propRes_r0.1<-replicate(13,numeric(),simplify=FALSE) propRes_rInf<-replicate(13,numeric(),simplify=FALSE) maxNodeSize_r0.01<-replicate(13,numeric(),simplify=FALSE) maxNodeSize_r0.1<-replicate(13,numeric(),simplify=FALSE) maxNodeSize_rInf<-replicate(13,numeric(),simplify=FALSE) polyDepth_r0.01<-replicate(13,numeric(),simplify=FALSE) polyDepth_r0.1<-replicate(13,numeric(),simplify=FALSE) polyDepth_rInf<-replicate(13,numeric(),simplify=FALSE) ntaxa_r0.01<-replicate(13,numeric(),simplify=FALSE) ntaxa_r0.1<-replicate(13,numeric(),simplify=FALSE) ntaxa_rInf<-replicate(13,numeric(),simplify=FALSE) #just budding k<-1 for(i in 1:nruns){ #generate taxon datasets and pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0,prop.bifurc=0, prop.cryptic=0,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0,prop.bifurc=0, prop.cryptic=0,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0,prop.bifurc=0, prop.cryptic=0,nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa1_",Sys.Date(),".Rdata",sep="")) #just bifurcation k<-2 for(i in 1:nruns){ #generate pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0,prop.bifurc=1, prop.cryptic=0,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0,prop.bifurc=1, prop.cryptic=0,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0,prop.bifurc=1, prop.cryptic=0,nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa2_",Sys.Date(),".Rdata",sep="")) #cryptic + anagenesis k<-3 for(i in 1:nruns){ #generate pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=0, prop.cryptic=1,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=0, prop.cryptic=1,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=0, prop.cryptic=1,nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa3_",Sys.Date(),".Rdata",sep="")) #budding+bifurcation k<-4 for(i in 1:nruns){ #generate pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0,prop.bifurc=0.5, prop.cryptic=0,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0,prop.bifurc=0.5, prop.cryptic=0,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0,prop.bifurc=0.5, prop.cryptic=0,nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa4_",Sys.Date(),".Rdata",sep="")) #budding with anagenesis k<-5 for(i in 1:nruns){ #generate pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=0, prop.cryptic=0,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=0, prop.cryptic=0,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=0, prop.cryptic=0,nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa5_",Sys.Date(),".Rdata",sep="")) #bifurcation with anagenesis k<-6 for(i in 1:nruns){ #generate pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=1, prop.cryptic=0,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=1, prop.cryptic=0,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=1, prop.cryptic=0,nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa6_",Sys.Date(),".Rdata",sep="")) #budding+bifurcation with anagenesis k<-7 for(i in 1:nruns){ #generate pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=0.5, prop.cryptic=0,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=0.5, prop.cryptic=0,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=0.5, prop.cryptic=0,nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa7_",Sys.Date(),".Rdata",sep="")) #budding+cryptic k<-8 for(i in 1:nruns){ #generate pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0,prop.bifurc=0, prop.cryptic=0.5,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0,prop.bifurc=0, prop.cryptic=0.5,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0,prop.bifurc=0, prop.cryptic=0.5,nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa8_",Sys.Date(),".Rdata",sep="")) #bifurcation+cryptic k<-9 for(i in 1:nruns){ #generate pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0,prop.bifurc=1, prop.cryptic=0.5,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0,prop.bifurc=1, prop.cryptic=0.5,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0,prop.bifurc=1, prop.cryptic=0.5,nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa9_",Sys.Date(),".Rdata",sep="")) #budding+bifurcation+cryptic k<-10 for(i in 1:nruns){ #generate pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0,prop.bifurc=(1/3), prop.cryptic=(1/3),nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0,prop.bifurc=(1/3), prop.cryptic=(1/3),nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0,prop.bifurc=(1/3), prop.cryptic=(1/3),nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa10_",Sys.Date(),".Rdata",sep="")) #budding+cryptic with anagenesis k<-11 for(i in 1:nruns){ #generate pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=0, prop.cryptic=0.5,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=0, prop.cryptic=0.5,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=0, prop.cryptic=0.5,nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa11_",Sys.Date(),".Rdata",sep="")) #bifurcation+cryptic with anagenesis k<-12 for(i in 1:nruns){ #generate pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=1, prop.cryptic=0.5,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=1, prop.cryptic=0.5,nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=1, prop.cryptic=0.5,nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa12_",Sys.Date(),".Rdata",sep="")) #budding+bifurcation+cryptic with anagenesis k<-13 for(i in 1:nruns){ #generate pruned cladograms for three sampling regimes #since I'm sampling, I must use drop.cryptic=F #when I'm not sampling, I can ignore and use drop.cryptic=T #first r=0.01 taxa<-simFossilTaxa_SRCond(r=0.01,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=(1/3), prop.cryptic=(1/3),nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.01)[,1])))) propRes_r0.01[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.01[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.01[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.01[[k]][i]<-max(nodeData[,1]) #next r=0.1 taxa<-simFossilTaxa_SRCond(r=0.1,avgtaxa=100, p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=(1/3), prop.cryptic=(1/3),nruns=1, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-drop.tip(taxa2cladogram(taxa,drop.crypt=FALSE), names(which(is.na(sampleRanges(taxa,r=0.1)[,1])))) propRes_r0.1[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_r0.1[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_r0.1[[k]][i]<-polyEvenness(nodeData) maxNodeSize_r0.1[[k]][i]<-max(nodeData[,1]) #next r=Inf taxa<-simFossilTaxa(p=0.1,q=0.1,anag.rate=0.1,prop.bifurc=(1/3), prop.cryptic=(1/3),nruns=1,mintaxa=100, maxtime=10^6,maxExtant=0,plot=F,print.runs=T) cladogram<-taxa2cladogram(taxa,drop.crypt=TRUE) propRes_rInf[[k]][i]<-(Nnode(cladogram)-1)/(Ntip(cladogram)-2) ntaxa_rInf[[k]][i]<-Ntip(cladogram) nodeData<-nodeInfo(cladogram) polyDepth_rInf[[k]][i]<-polyEvenness(nodeData) maxNodeSize_rInf[[k]][i]<-max(nodeData[,1]) } dput(propRes_r0.01,output_propres_r0.01) dput(propRes_r0.1,output_propres_r0.1) dput(propRes_rInf,output_propres_rInf) dput(polyDepth_r0.01,output_polydepth_r0.01) dput(polyDepth_r0.1,output_polydepth_r0.1) dput(polyDepth_rInf,output_polydepth_rInf) dput(maxNodeSize_r0.01,output_maxnode_r0.01) dput(maxNodeSize_r0.1,output_maxnode_r0.1) dput(maxNodeSize_rInf,output_maxnode_rInf) dput(ntaxa_r0.01,output_ntaxa_r0.01) dput(ntaxa_r0.1,output_ntaxa_r0.1) dput(ntaxa_rInf,output_ntaxa_rInf) save.image(paste("clado_taxa13_",Sys.Date(),".Rdata",sep=""))