setwd("C:/R2/projects/RPANDA on all primates") library(ape); library(geiger); library(phytools);library(ggplot2);library(plotrix); library(picante); library(diversitree); library(utils);library(reshape); library(paleotree); library(expm);library(devtools); memory.limit(size=1000000) load("C:/R2/projects/RPANDA on all primates/RPANDA on all primates springer trees_w_quasse2.RData") ##read in springer trees trees<-read.tree(file="springer_trees.tre") trees plot(ladderize(trees[[1]])) axisPhylo() branching.times(trees[[4]])[1] ##save a backup trees_backup<-trees ##trees are scaled so that 1 = 100 million years - recalibrating to be in millions of years trees_timescaled<-vector("list") class(trees_timescaled)<-"multiPhylo" for(i in 1:length(trees)){ trees_timescaled[[i]]<-rescale(trees[[i]],model="depth",branching.times(trees[[i]])[1]*100) } trees_timescaled branching.times(trees_timescaled[[1]]) plot(trees_timescaled[[1]],cex=0.5) edgelabels(cex=0.5,frame="n") axisPhylo() is.extinct(trees_timescaled[[1]]) is.ultrametric(trees_timescaled[[1]]) trees<-trees_timescaled plot(trees[[1]]) edgelabels(cex=0.5,frame="n") nodelabels(cex=0.5,frame="n") trees[[1]] branching.times(trees[[1]]) par(mfrow=c(2,2)) plot(trees[[1]],show.tip.label = F) axisPhylo() plot(trees[[2]],show.tip.label = F) axisPhylo() plot(trees[[3]],show.tip.label = F) axisPhylo() plot(trees[[4]],show.tip.label = F) axisPhylo() tree<-trees[[1]] ############################### ###doing a for loop to run RPANDA on each taxonomic group through each tree in the file ############################### ##create the results objects to be populated by the loop mod.avg.params.results<-matrix(ncol=5) mod.avg.params.results<-as.data.frame(mod.avg.params.results) colnames(mod.avg.params.results)<-c("lambda","mu","time-var.par","clade","treenum") mod.fits.results<-matrix(ncol=7) mod.fits.results<-as.data.frame(mod.fits.results) colnames(mod.fits.results)<-c("lambda","mu","time.var.par","fit","delta","w","taxon") ######### begin for loop for(i in 1:length(trees)){ #################################################################################### ### fitting the RPANDA models to each taxonomic group ########################################## #################################################################################### ##isolate streps library(phytools) tree<-trees[[i]] ##isolate primates primates_node<-findMRCA(tree,tips=c("Daubentonia_madagascariensis","Homo_sapiens")) tree<-extract.clade(tree, primates_node, root.edge=1) ##isolate strepsirrhines streps_node<-findMRCA(tree,tips=c("Daubentonia_madagascariensis","Loris_tardigradus")) streps<-extract.clade(tree, streps_node, root.edge=1) ##isolate New World Monkeys (NWM) NWM_node<-findMRCA(tree,tips=c("Callicebus_moloch","Cebus_apella")) NWM<-extract.clade(tree, NWM_node, root.edge=1) ##isolate anthropoids anthropoids_node<-findMRCA(tree,tips=c("Papio_anubis","Pan_troglodytes")) anthropoids<-extract.clade(tree, anthropoids_node, root.edge=1) ##isolate Old World Monkeys (OWM) OWM_node<-findMRCA(tree,tips=c("Papio_anubis","Nasalis_larvatus")) OWM<-extract.clade(tree, OWM_node, root.edge=1) ##isolate tarsiers tars_node<-findMRCA(tree,tips=c("Tarsius_syrichta","Tarsius_tarsier")) tars<-extract.clade(tree, tars_node,root.edge = 1) ##isolate apes apes_node<-findMRCA(tree,tips=c("Pan_paniscus","Hylobates_lar")) apes<-extract.clade(tree, apes_node,root.edge = 1) ################################################ ## first for all primates tot_time<-max(node.age(tree)$ages)+tree$root.edge ######################################################################### # Fit the pure birth model (no extinction) with a constant speciation rate ######################################################################### library(RPANDA) f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){0} lamb_par<-c(0.15) mu_par<-c() result_cst_primates <- fit_bd(tree,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(tree$tip.label)/422,cst.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_cst_primates$model <- "pure birth with constant speciation rate" #plot_fit_bd(result_cst_streps, tot_time) ################################################################# # Fit the pure birth model (no extinction) with exponential variation (decrease) # of the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] * exp(y[2] * t)} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.05) mu_par<-c() result_exp_tree <- fit_bd(tree,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(tree$tip.label)/422,expo.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_exp_tree$model <- "pure birth with exponential variation in speciation rate" #windows() #par(mfrow=c(1,3)) #plot_fit_bd(result_exp_streps, tot_time) # plot_dtt(result_exp_tree,tot_time,N0=422) ################################################################# # Fit the a pure birth model (no extinction) with linear variation (decrease) of # the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] + y[2] * t} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.001) mu_par<-c() result_lin_tree <- fit_bd(tree,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(tree$tip.label)/422,fix.mu=TRUE,dt=1e-3) result_lin_tree$model <- "pure birth with linear variation in speciation rate" #plot_fit_bd(result_lin_streps) ################################################################# # Fit a birth-death model with constant speciation and constant extinction ################################################################# f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){y[1]} lamb_par <- c(0.15) mu_par <-c(0.05) result_bcst_dcst_tree <- fit_bd(tree,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(tree$tip.label)/422,cst.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bcst_dcst_tree$model <- "birth-death with Constant speciation rate and constant extinction" #par(mfrow=c(1,4)) #plot_fit_bd(result_bcst_dcst_tree) # plot_dtt(result_bcst_dcst_tree,tot_time,N0=422) ################################################################# # Fit a birth-death model with exponential variation of the speciation # rate with time and constant extinction ################################################################# f.lamb<-function(t,y){y[1] * exp(y[2] * t)} f.mu <-function(t,y){y[1]} lamb_par <- c(0.15, 0.05) mu_par <-c(0.05) result_bexp_dcst_tree <- fit_bd(tree,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(tree$tip.label)/422,expo.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bexp_dcst_tree$model <- "birth-death with exponential variation in speciation rate and constant extinction" #par(mfrow=c(1,3)) #plot_fit_bd(result_bexp_dcst_tree) # plot_dtt(result_bexp_dcst_tree,tot_time,N0=422) ################################################################# # Fit a birth-death model with exponential variation (decrease) of the extinction # rate with time and constant speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1]} mu_par <- c(0.05, 0.01) lambda_par <-c(0.15) result_bcst_dexp_tree <- fit_bd(tree,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(tree$tip.label)/422,cst.lamb=TRUE,expo.mu=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bcst_dexp_tree$model <- "birth-death with exponential decrease in extinction rate and constant speciation" #par(mfrow=c(4,1)) #plot_fit_bd(result_bcst_dexp_streps) #plot_dtt(result_bcst_dexp_tree,tot_time,N0=422) ################################################################# # Fit a birth-death model with exponential variation of the extinction # rate with time and exponential speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1] * exp(y[2] * t)} mu_par <- c(0.05, 0.01) lambda_par <-c(0.15, 0.01) result_bexp_dexp_tree <- fit_bd(tree,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(tree$tip.label)/422,expo.lamb=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bexp_dexp_tree$model <- "birth-death with exponential extinction rate and exponential speciation" #par(mfrow=c(4,1)) #plot_fit_bd(result_bexp_dexp_tree) # plot_dtt(result_bexp_dexp_tree,tot_time,N0=422) ################################################################# # Find the best model ################################################################# aiccs_tree <- c(result_cst_primates$aicc, result_exp_tree$aicc, result_lin_tree$aicc,result_bcst_dcst_tree$aicc,result_bexp_dcst_tree$aicc, result_bcst_dexp_tree$aicc, result_bexp_dexp_tree$aicc) aicws_tree<-aicw(aiccs_tree) #write.table(aicws_2, file="tree diversification model comp rpanda.csv", sep=",") ###model averaging #lambda lambdas_tree_tree <- c(result_cst_primates$lamb_par, result_exp_tree$lamb_par[1], result_lin_tree$lamb_par[1],result_bcst_dcst_tree$lamb_par[1], result_bexp_dcst_tree$lamb_par[1], result_bcst_dexp_tree$lamb_par[1], result_bexp_dexp_tree$lamb_par[1]) lambdas_tree_tree<-as.data.frame(lambdas_tree_tree) mod.names <- c(result_cst_primates$model, result_exp_tree$model, result_lin_tree$model,result_bcst_dcst_tree$model, result_bexp_dcst_tree$model, result_bcst_dexp_tree$model, result_bexp_dexp_tree$model) rownames(lambdas_tree_tree)<-mod.names lambdas_aicws_tree_tree<-cbind(lambdas_tree_tree, aicws_tree) lambdas_aicws_tree_tree str(lambdas_aicws_tree_tree) d.aics<-cbind(lambdas_aicws_tree_tree[1,2]-lambdas_aicws_tree_tree[2,2],lambdas_aicws_tree_tree[4,2]-lambdas_aicws_tree_tree[5,2]) colnames(d.aics)<-c("d.aic.pb","d.aic.bd") mod.av.lambda<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.lambda)<-"weighted.lambda" for (i in 1:length(lambdas_aicws_tree_tree$w)){ mod.av.lambda[i,1]<-lambdas_aicws_tree_tree[i,1]*lambdas_aicws_tree_tree[i,4] } mod.av.lambda[,1] sum(mod.av.lambda[,1]) #mu mu_tree_tree <- c(0, 0, 0,result_bcst_dcst_tree$mu_par[1], result_bexp_dcst_tree$mu_par[1], result_bcst_dexp_tree$mu_par[1], result_bexp_dexp_tree$mu_par[1]) mu_tree_tree<-as.data.frame(mu_tree_tree) mu_tree_tree colnames(mu_tree_tree)<-"mu_tree_tree" mu_aicws_tree_tree<-cbind(mu_tree_tree, aicws_tree) mu_aicws_tree_tree mod.av.mu<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.mu)<-"weighted.mu" for (i in 1:length(mu_aicws_tree_tree$w)){ mod.av.mu[i,1]<-mu_aicws_tree_tree[i,1]*mu_aicws_tree_tree[i,4] } mod.av.mu[,1] sum(mod.av.mu[,1]) #time varying parameter tim.var.par_tree_tree <- c(result_cst_primates$lamb_par[2], result_exp_tree$lamb_par[2], result_lin_tree$lamb_par[2],result_bcst_dcst_tree$lamb_par[2], result_bexp_dcst_tree$lamb_par[2], result_bcst_dexp_tree$lamb_par[2], result_bexp_dexp_tree$lamb_par[2]) tim.var.par_tree_tree<-as.data.frame(tim.var.par_tree_tree) tim.var.par_tree_tree tim.var.par_aicws_tree_tree<-cbind(tim.var.par_tree_tree, aicws_tree) tim.var.par_aicws_tree_tree mod.av.tim.var.par<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.tim.var.par)<-"weighted.tim.var.par" for (i in 1:length(tim.var.par_aicws_tree_tree$w)){ mod.av.tim.var.par[i,1]<-tim.var.par_aicws_tree_tree[i,1]*tim.var.par_aicws_tree_tree[i,4] } mod.av.tim.var.par[,1] sum(mod.av.tim.var.par[,1],na.rm=T) lam_mu_timvar_aicws_tree_tree<-cbind(lambdas_tree_tree,mu_tree_tree, tim.var.par_tree_tree, aicws_tree) str(lam_mu_timvar_aicws_tree_tree) lam_mu_timvar_aicws_tree_tree$taxon<-rep("tree",length(lam_mu_timvar_aicws_tree_tree$lambdas_tree_tree)) rownames(lam_mu_timvar_aicws_tree_tree)<-mod.names colnames(lam_mu_timvar_aicws_tree_tree)<-c("lambda","mu","time.var.par","fit","delta","w","taxon") mod_av_res_tree<-cbind(sum(mod.av.lambda[,1]),sum(mod.av.mu[,1]),sum(mod.av.tim.var.par[,1],na.rm=T)) #print(mod_av_res_tree) ################################################ ## now for streps tot_time<-max(node.age(streps)$ages)+streps$root.edge ######################################################################### # Fit the pure birth model (no extinction) with a constant speciation rate ######################################################################### library(RPANDA) f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){0} lamb_par<-c(0.15) mu_par<-c() result_cst_streps <- fit_bd(streps,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(streps$tip.label)/128,cst.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_cst_streps$model <- "pure birth with constant speciation rate" #plot_fit_bd(result_cst_streps, tot_time) ################################################################# # Fit the pure birth model (no extinction) with exponential variation (decrease) # of the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] * exp(y[2] * t)} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.05) mu_par<-c() result_exp_streps <- fit_bd(streps,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(streps$tip.label)/128,expo.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_exp_streps$model <- "pure birth with exponential variation in speciation rate" #windows() #par(mfrow=c(1,3)) # plot_fit_bd(result_exp_streps, tot_time) ################################################################# # Fit the a pure birth model (no extinction) with linear variation (decrease) of # the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] + y[2] * t} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.001) mu_par<-c() result_lin_streps <- fit_bd(streps,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(streps$tip.label)/128,fix.mu=TRUE,dt=1e-3) result_lin_streps$model <- "pure birth with linear variation in speciation rate" #plot_fit_bd(result_lin_streps) ################################################################# # Fit a birth-death model with constant speciation and constant extinction ################################################################# f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){y[1]} lamb_par <- c(0.15) mu_par <-c(0.05) result_bcst_dcst_streps <- fit_bd(streps,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(streps$tip.label)/128,cst.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bcst_dcst_streps$model <- "birth-death with Constant speciation rate and constant extinction" #par(mfrow=c(1,4)) #plot_fit_bd(result_bcst_dcst_streps) ################################################################# # Fit a birth-death model with exponential variation of the speciation # rate with time and constant extinction ################################################################# f.lamb<-function(t,y){y[1] * exp(y[2] * t)} f.mu <-function(t,y){y[1]} lamb_par <- c(0.15, 0.05) mu_par <-c(0.05) result_bexp_dcst_streps <- fit_bd(streps,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(streps$tip.label)/128,expo.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bexp_dcst_streps$model <- "birth-death with exponential variation in speciation rate and constant extinction" #par(mfrow=c(1,3)) # plot_fit_bd(result_bexp_dcst_streps) ################################################################# # Fit a birth-death model with exponential variation (decrease) of the extinction # rate with time and constant speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1]} mu_par <- c(0.05, 0.01) lambda_par <-c(0.15) result_bcst_dexp_streps <- fit_bd(streps,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(streps$tip.label)/128,cst.lamb=TRUE,expo.mu=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bcst_dexp_streps$model <- "birth-death with exponential decrease in extinction rate and constant speciation" #par(mfrow=c(4,1)) #plot_fit_bd(result_bcst_dexp_streps) #plot_dtt(result_bcst_dexp_streps, tot_time, 128) ################################################################# # Fit a birth-death model with exponential variation of the extinction # rate with time and exponential speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1] * exp(y[2] * t)} mu_par <- c(0.05, 0.01) lambda_par <-c(0.15, 0.01) result_bexp_dexp_streps <- fit_bd(streps,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(streps$tip.label)/128,expo.lamb=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bexp_dexp_streps$model <- "birth-death with exponential extinction rate and exponential speciation" #par(mfrow=c(4,1)) # plot_fit_bd(result_bexp_dexp_streps) # plot_dtt(result_bexp_dexp_streps, tot_time, 128) ################################################################# # Find the best model ################################################################# aiccs_streps <- c(result_cst_streps$aicc, result_exp_streps$aicc, result_lin_streps$aicc,result_bcst_dcst_streps$aicc,result_bexp_dcst_streps$aicc, result_bcst_dexp_streps$aicc, result_bexp_dexp_streps$aicc) aicws_streps<-aicw(aiccs_streps) #write.table(aicws_2, file="streps diversification model comp rpanda.csv", sep=",") ###model averaging #lambda lambdas_streps_tree <- c(result_cst_streps$lamb_par, result_exp_streps$lamb_par[1], result_lin_streps$lamb_par[1],result_bcst_dcst_streps$lamb_par[1], result_bexp_dcst_streps$lamb_par[1], result_bcst_dexp_streps$lamb_par[1], result_bexp_dexp_streps$lamb_par[1]) lambdas_streps_tree<-as.data.frame(lambdas_streps_tree) mod.names <- c(result_cst_streps$model, result_exp_streps$model, result_lin_streps$model,result_bcst_dcst_streps$model, result_bexp_dcst_streps$model, result_bcst_dexp_streps$model, result_bexp_dexp_streps$model) rownames(lambdas_streps_tree)<-mod.names lambdas_aicws_streps_tree<-cbind(lambdas_streps_tree, aicws_streps) lambdas_aicws_streps_tree str(lambdas_aicws_streps_tree) d.aics<-cbind(lambdas_aicws_streps_tree[1,2]-lambdas_aicws_streps_tree[2,2],lambdas_aicws_streps_tree[4,2]-lambdas_aicws_streps_tree[5,2]) colnames(d.aics)<-c("d.aic.pb","d.aic.bd") mod.av.lambda<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.lambda)<-"weighted.lambda" for (i in 1:length(lambdas_aicws_streps_tree$w)){ mod.av.lambda[i,1]<-lambdas_aicws_streps_tree[i,1]*lambdas_aicws_streps_tree[i,4] } mod.av.lambda[,1] sum(mod.av.lambda[,1]) #mu mu_streps_tree <- c(0, 0, 0,result_bcst_dcst_streps$mu_par[1], result_bexp_dcst_streps$mu_par[1], result_bcst_dexp_streps$mu_par[1], result_bexp_dexp_streps$mu_par[1]) mu_streps_tree<-as.data.frame(mu_streps_tree) mu_streps_tree colnames(mu_streps_tree)<-"mu_streps_tree" mu_aicws_streps_tree<-cbind(mu_streps_tree, aicws_streps) mu_aicws_streps_tree mod.av.mu<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.mu)<-"weighted.mu" for (i in 1:length(mu_aicws_streps_tree$w)){ mod.av.mu[i,1]<-mu_aicws_streps_tree[i,1]*mu_aicws_streps_tree[i,4] } mod.av.mu[,1] sum(mod.av.mu[,1]) #time varying parameter tim.var.par_streps_tree <- c(result_cst_streps$lamb_par[2], result_exp_streps$lamb_par[2], result_lin_streps$lamb_par[2],result_bcst_dcst_streps$lamb_par[2], result_bexp_dcst_streps$lamb_par[2], result_bcst_dexp_streps$lamb_par[2], result_bexp_dexp_streps$lamb_par[2]) tim.var.par_streps_tree<-as.data.frame(tim.var.par_streps_tree) tim.var.par_streps_tree tim.var.par_aicws_streps_tree<-cbind(tim.var.par_streps_tree, aicws_streps) tim.var.par_aicws_streps_tree mod.av.tim.var.par<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.tim.var.par)<-"weighted.tim.var.par" for (i in 1:length(tim.var.par_aicws_streps_tree$w)){ mod.av.tim.var.par[i,1]<-tim.var.par_aicws_streps_tree[i,1]*tim.var.par_aicws_streps_tree[i,4] } mod.av.tim.var.par[,1] sum(mod.av.tim.var.par[,1],na.rm=T) lam_mu_timvar_aicws_streps_tree<-cbind(lambdas_streps_tree,mu_streps_tree, tim.var.par_streps_tree, aicws_streps) str(lam_mu_timvar_aicws_streps_tree) lam_mu_timvar_aicws_streps_tree$taxon<-rep("streps",length(lam_mu_timvar_aicws_streps_tree$lambdas_streps_tree)) rownames(lam_mu_timvar_aicws_streps_tree)<-mod.names colnames(lam_mu_timvar_aicws_streps_tree)<-c("lambda","mu","time.var.par","fit","delta","w","taxon") mod_av_res_streps<-cbind(sum(mod.av.lambda[,1]),sum(mod.av.mu[,1]),sum(mod.av.tim.var.par[,1],na.rm=T)) #print(mod_av_res_streps) #################################################### ##RPANDA models for tarsiers #################################################### tot_time<-max(node.age(tars)$ages)+tars$root.edge ######################################################################### # Fit the pure birth model (no extinction) with a constant speciation rate ######################################################################### f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){0} lamb_par<-c(0.15) mu_par<-c() result_cst_tars <- fit_bd(tars,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(tars$tip.label)/10,cst.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_cst_tars$model <- "pure birth with constant speciation rate" #plot_fit_bd(result_cst_tars, tot_time) ################################################################# # Fit the pure birth model (no extinction) with exponential variation (decrease) # of the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] * exp(y[2] * t)} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.05) mu_par<-c() result_exp_tars <- fit_bd(tars,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(tars$tip.label)/10,expo.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_exp_tars$model <- "pure birth with exponential variation in speciation rate" #par(mfrow=c(1,3)) #plot_fit_bd(result_exp_tars, tot_time) ################################################################# # Fit the a pure birth model (no extinction) with linear variation (decrease) of # the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] + y[2] * t} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.001) mu_par<-c() result_lin_tars <- fit_bd(tars,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(tars$tip.label)/10,fix.mu=TRUE,dt=1e-3,cond="stem") result_lin_tars$model <- "pure birth with linear variation in speciation rate" #plot_fit_bd(result_lin_tars) ################################################################# # Fit a birth-death model with constant speciation and constant extinction ################################################################# f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){y[1]} lamb_par <- c(0.15) mu_par <-c(0.05) result_bcst_dcst_tars <- fit_bd(tars,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(tars$tip.label)/10,cst.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bcst_dcst_tars$model <- "birth-death with Constant speciation rate and constant extinction" #par(mfrow=c(1,4)) #plot_fit_bd(result_bcst_dcst_tars) ################################################################# # Fit a birth-death model with exponential variation of the speciation # rate with time and constant extinction ################################################################# f.lamb<-function(t,y){y[1] * exp(y[2] * t)} f.mu <-function(t,y){y[1]} lamb_par <- c(0.15, 0.05) mu_par <-c(0.05) result_bexp_dcst_tars <- fit_bd(tars,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(tars$tip.label)/10,expo.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bexp_dcst_tars$model <- "birth-death with exponential variation in speciation rate and constant extinction" #par(mfrow=c(1,3)) #plot_fit_bd(result_bexp_dcst_tars) ################################################################# # Fit a birth-death model with exponential variation (decrease) of the extinction # rate with time and constant speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1]} mu_par <- c(0.03, 0.01) lambda_par <-c(0.07) result_bcst_dexp_tars <- fit_bd(tars,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(tars$tip.label)/10,cst.lamb=TRUE,expo.mu=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bcst_dexp_tars$model <- "birth-death with exponential decrease in extinction rate and constant speciation" #par(mfrow=c(3,3)) #plot_fit_bd(result_bcst_dexp_tars) #plot_dtt(result_bcst_dexp_tars, tot_time, 10) ################################################################# # Fit a birth-death model with exponential variation of the extinction # rate with time and exponential speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1] * exp(y[2] * t)} mu_par <- c(0.05, 0.01) lambda_par <-c(0.15, 0.01) result_bexp_dexp_tars <- fit_bd(tars,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(tars$tip.label)/10,expo.lamb=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bexp_dexp_tars$model <- "birth-death with exponential extinction rate and exponential speciation" #par(mfrow=c(3,1)) #plot_fit_bd(result_bexp_dexp_tars) ################################################################# # Find the best model ################################################################# aiccs_tars <- c(result_cst_tars$aicc, result_exp_tars$aicc, result_lin_tars$aicc,result_bcst_dcst_tars$aicc,result_bexp_dcst_tars$aicc, result_bcst_dexp_tars$aicc, result_bexp_dexp_tars$aicc) aicws_tars<-aicw(aiccs_tars) #write.table(aicws_3, file="tars model comp rpanda.csv", sep=",") ###model averaging #lambda lambdas_tars_tree <- c(result_cst_tars$lamb_par, result_exp_tars$lamb_par[1], result_lin_tars$lamb_par[1],result_bcst_dcst_tars$lamb_par[1], result_bexp_dcst_tars$lamb_par[1], result_bcst_dexp_tars$lamb_par[1], result_bexp_dexp_tars$lamb_par[1]) lambdas_tars_tree<-as.data.frame(lambdas_tars_tree) lambdas_tars_tree lambdas_aicws_tars_tree<-cbind(lambdas_tars_tree, aicws_tars) lambdas_aicws_tars_tree mod.names <- c(result_cst_tars$model, result_exp_tars$model, result_lin_tars$model,result_bcst_dcst_tars$model, result_bexp_dcst_tars$model, result_bcst_dexp_tars$model, result_bexp_dexp_tars$model) rownames(lambdas_aicws_tars_tree)<-mod.names d.aics<-cbind(lambdas_aicws_tars_tree[1,2]-lambdas_aicws_tars_tree[2,2],lambdas_aicws_tars_tree[4,2]-lambdas_aicws_tars_tree[5,2]) colnames(d.aics)<-c("d.aic.pb","d.aic.bd") mod.av.lambda<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.lambda)<-"weighted.lambda" for (i in 1:length(lambdas_aicws_tars_tree$w)){ mod.av.lambda[i,1]<-lambdas_aicws_tars_tree[i,1]*lambdas_aicws_tars_tree[i,4] } mod.av.lambda[,1] sum(mod.av.lambda[,1]) #mu mu_tars_tree <- c(0, 0, 0,result_bcst_dcst_tars$mu_par[1], result_bexp_dcst_tars$mu_par[1], result_bcst_dexp_tars$mu_par[1], result_bexp_dexp_tars$mu_par[1]) mu_tars_tree<-as.data.frame(mu_tars_tree) mu_tars_tree colnames(mu_tars_tree)<-"mu" mu_aicws_tars_tree<-cbind(mu_tars_tree, aicws_tars) mu_aicws_tars_tree mod.av.mu<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.mu)<-"weighted.mu" for (i in 1:length(mu_aicws_tars_tree$w)){ mod.av.mu[i,1]<-mu_aicws_tars_tree[i,1]*mu_aicws_tars_tree[i,4] } mod.av.mu[,1] sum(mod.av.mu[,1]) #time varying parameter tim.var.par_tars_tree <- c(result_cst_tars$lamb_par[2], result_exp_tars$lamb_par[2], result_lin_tars$lamb_par[2],result_bcst_dcst_tars$lamb_par[2], result_bexp_dcst_tars$lamb_par[2], result_bcst_dexp_tars$lamb_par[2], result_bexp_dexp_tars$lamb_par[2]) tim.var.par_tars_tree<-as.data.frame(tim.var.par_tars_tree) tim.var.par_tars_tree tim.var.par_aicws_tars_tree<-cbind(tim.var.par_tars_tree, aicws_tars) tim.var.par_aicws_tars_tree mod.av.tim.var.par<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.tim.var.par)<-"weighted.tim.var.par" for (i in 1:length(tim.var.par_aicws_tars_tree$w)){ mod.av.tim.var.par[i,1]<-tim.var.par_aicws_tars_tree[i,1]*tim.var.par_aicws_tars_tree[i,4] } mod.av.tim.var.par[,1] sum(mod.av.tim.var.par[,1],na.rm=T) lam_mu_timvar_aicws_tars_tree<-cbind(lambdas_tars_tree,mu_tars_tree, tim.var.par_tars_tree, aicws_tars) str(lam_mu_timvar_aicws_tars_tree) lam_mu_timvar_aicws_tars_tree$taxon<-rep("tars",length(lam_mu_timvar_aicws_tars_tree$lambdas_tars_tree)) rownames(lam_mu_timvar_aicws_tars_tree)<-mod.names colnames(lam_mu_timvar_aicws_tars_tree)<-c("lambda","mu","time.var.par","fit","delta","w","taxon") mod_av_res_tars<-cbind(sum(mod.av.lambda[,1]),sum(mod.av.mu[,1]),sum(mod.av.tim.var.par[,1],na.rm=T)) #################################################### ##RPANDA models for NWM #################################################### tot_time<-max(node.age(NWM)$ages)+NWM$root.edge ######################################################################### # Fit the pure birth model (no extinction) with a constant speciation rate ######################################################################### f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){0} lamb_par<-c(0.15) mu_par<-c() result_cst_NWM <- fit_bd(NWM,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(NWM$tip.label)/142,cst.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_cst_NWM$model <- "pure birth with constant speciation rate" #plot_fit_bd(result_cst_NWM, tot_time) ################################################################# # Fit the pure birth model (no extinction) with exponential variation (decrease) # of the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] * exp(y[2] * t)} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.05) mu_par<-c() result_exp_NWM <- fit_bd(NWM,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(NWM$tip.label)/142,expo.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_exp_NWM$model <- "pure birth with exponential variation in speciation rate" #par(mfrow=c(1,3)) #plot_fit_bd(result_exp_NWM, tot_time) ################################################################# # Fit the a pure birth model (no extinction) with linear variation (decrease) of # the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] + y[2] * t} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.001) mu_par<-c() result_lin_NWM <- fit_bd(NWM,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(NWM$tip.label)/142,fix.mu=TRUE,dt=1e-3,cond="stem") result_lin_NWM$model <- "pure birth with linear variation in speciation rate" #plot_fit_bd(result_lin_NWM) ################################################################# # Fit a birth-death model with constant speciation and constant extinction ################################################################# f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){y[1]} lamb_par <- c(0.15) mu_par <-c(0.05) result_bcst_dcst_NWM <- fit_bd(NWM,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(NWM$tip.label)/142,cst.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bcst_dcst_NWM$model <- "birth-death with Constant speciation rate and constant extinction" #par(mfrow=c(1,4)) #plot_fit_bd(result_bcst_dcst_NWM) ################################################################# # Fit a birth-death model with exponential variation of the speciation # rate with time and constant extinction ################################################################# f.lamb<-function(t,y){y[1] * exp(y[2] * t)} f.mu <-function(t,y){y[1]} lamb_par <- c(0.15, 0.05) mu_par <-c(0.05) result_bexp_dcst_NWM <- fit_bd(NWM,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(NWM$tip.label)/142,expo.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bexp_dcst_NWM$model <- "birth-death with exponential variation in speciation rate and constant extinction" #par(mfrow=c(1,3)) #plot_fit_bd(result_bexp_dcst_NWM) ################################################################# # Fit a birth-death model with exponential variation (decrease) of the extinction # rate with time and constant speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1]} mu_par <- c(0.03, 0.01) lambda_par <-c(0.07) result_bcst_dexp_NWM <- fit_bd(NWM,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(NWM$tip.label)/142,cst.lamb=TRUE,expo.mu=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bcst_dexp_NWM$model <- "birth-death with exponential decrease in extinction rate and constant speciation" #par(mfrow=c(3,3)) #plot_fit_bd(result_bcst_dexp_NWM) #plot_dtt(result_bcst_dexp_NWM, tot_time, 142) ################################################################# # Fit a birth-death model with exponential variation of the extinction # rate with time and exponential speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1] * exp(y[2] * t)} mu_par <- c(0.05, 0.01) lambda_par <-c(0.15, 0.01) result_bexp_dexp_NWM <- fit_bd(NWM,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(NWM$tip.label)/142,expo.lamb=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bexp_dexp_NWM$model <- "birth-death with exponential extinction rate and exponential speciation" #par(mfrow=c(3,1)) #plot_fit_bd(result_bexp_dexp_NWM) ################################################################# # Find the best model ################################################################# aiccs_NWM <- c(result_cst_NWM$aicc, result_exp_NWM$aicc, result_lin_NWM$aicc,result_bcst_dcst_NWM$aicc,result_bexp_dcst_NWM$aicc, result_bcst_dexp_NWM$aicc, result_bexp_dexp_NWM$aicc) aicws_NWM<-aicw(aiccs_NWM) #write.table(aicws_3, file="NWM model comp rpanda.csv", sep=",") ###model averaging #lambda lambdas_NWM_tree <- c(result_cst_NWM$lamb_par, result_exp_NWM$lamb_par[1], result_lin_NWM$lamb_par[1],result_bcst_dcst_NWM$lamb_par[1], result_bexp_dcst_NWM$lamb_par[1], result_bcst_dexp_NWM$lamb_par[1], result_bexp_dexp_NWM$lamb_par[1]) lambdas_NWM_tree<-as.data.frame(lambdas_NWM_tree) lambdas_NWM_tree lambdas_aicws_NWM_tree<-cbind(lambdas_NWM_tree, aicws_NWM) lambdas_aicws_NWM_tree mod.names <- c(result_cst_NWM$model, result_exp_NWM$model, result_lin_NWM$model,result_bcst_dcst_NWM$model, result_bexp_dcst_NWM$model, result_bcst_dexp_NWM$model, result_bexp_dexp_NWM$model) rownames(lambdas_aicws_NWM_tree)<-mod.names d.aics<-cbind(lambdas_aicws_NWM_tree[1,2]-lambdas_aicws_NWM_tree[2,2],lambdas_aicws_NWM_tree[4,2]-lambdas_aicws_NWM_tree[5,2]) colnames(d.aics)<-c("d.aic.pb","d.aic.bd") mod.av.lambda<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.lambda)<-"weighted.lambda" for (i in 1:length(lambdas_aicws_NWM_tree$w)){ mod.av.lambda[i,1]<-lambdas_aicws_NWM_tree[i,1]*lambdas_aicws_NWM_tree[i,4] } mod.av.lambda[,1] sum(mod.av.lambda[,1]) #mu mu_NWM_tree <- c(0, 0, 0,result_bcst_dcst_NWM$mu_par[1], result_bexp_dcst_NWM$mu_par[1], result_bcst_dexp_NWM$mu_par[1], result_bexp_dexp_NWM$mu_par[1]) mu_NWM_tree<-as.data.frame(mu_NWM_tree) mu_NWM_tree colnames(mu_NWM_tree)<-"mu_NWM_tree" mu_aicws_NWM_tree<-cbind(mu_NWM_tree, aicws_NWM) mu_aicws_NWM_tree mod.av.mu<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.mu)<-"weighted.mu" for (i in 1:length(mu_aicws_NWM_tree$w)){ mod.av.mu[i,1]<-mu_aicws_NWM_tree[i,1]*mu_aicws_NWM_tree[i,4] } mod.av.mu[,1] sum(mod.av.mu[,1]) #time varying parameter tim.var.par_NWM_tree <- c(result_cst_NWM$lamb_par[2], result_exp_NWM$lamb_par[2], result_lin_NWM$lamb_par[2],result_bcst_dcst_NWM$lamb_par[2], result_bexp_dcst_NWM$lamb_par[2], result_bcst_dexp_NWM$lamb_par[2], result_bexp_dexp_NWM$lamb_par[2]) tim.var.par_NWM_tree<-as.data.frame(tim.var.par_NWM_tree) tim.var.par_NWM_tree tim.var.par_aicws_NWM_tree<-cbind(tim.var.par_NWM_tree, aicws_NWM) tim.var.par_aicws_NWM_tree mod.av.tim.var.par<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.tim.var.par)<-"weighted.tim.var.par" for (i in 1:length(tim.var.par_aicws_NWM_tree$w)){ mod.av.tim.var.par[i,1]<-tim.var.par_aicws_NWM_tree[i,1]*tim.var.par_aicws_NWM_tree[i,4] } mod.av.tim.var.par[,1] sum(mod.av.tim.var.par[,1],na.rm=T) lam_mu_timvar_aicws_NWM_tree<-cbind(lambdas_NWM_tree,mu_NWM_tree, tim.var.par_NWM_tree, aicws_NWM) str(lam_mu_timvar_aicws_NWM_tree) lam_mu_timvar_aicws_NWM_tree$taxon<-rep("NWM",length(lam_mu_timvar_aicws_NWM_tree$lambdas_NWM_tree)) lam_mu_timvar_aicws_NWM_tree rownames(lam_mu_timvar_aicws_NWM_tree)<-mod.names colnames(lam_mu_timvar_aicws_NWM_tree)<-c("lambda","mu","time.var.par","fit","delta","w","taxon") mod_av_res_NWM<-cbind(sum(mod.av.lambda[,1]),sum(mod.av.mu[,1]),sum(mod.av.tim.var.par[,1],na.rm=T)) #################################################### ##RPANDA models for OWM #################################################### tot_time<-max(node.age(anthropoids)$ages)+anthropoids$root.edge ######################################################################### # Fit the pure birth model (no extinction) with a constant speciation rate ######################################################################### f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){0} lamb_par<-c(0.15) mu_par<-c() result_cst_anthropoids <- fit_bd(anthropoids,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(anthropoids$tip.label)/123,cst.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_cst_anthropoids$model <- "pure birth with constant speciation rate" #plot_fit_bd(result_cst_anthropoids, tot_time) ################################################################# # Fit the pure birth model (no extinction) with exponential variation (decrease) # of the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] * exp(y[2] * t)} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.05) mu_par<-c() result_exp_anthropoids <- fit_bd(anthropoids,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(anthropoids$tip.label)/123,expo.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_exp_anthropoids$model <- "pure birth with exponential variation in speciation rate" #par(mfrow=c(1,3)) #plot_fit_bd(result_exp_anthropoids, tot_time) ################################################################# # Fit the a pure birth model (no extinction) with linear variation (decrease) of # the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] + y[2] * t} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.001) mu_par<-c() result_lin_anthropoids <- fit_bd(anthropoids,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(anthropoids$tip.label)/123,fix.mu=TRUE,dt=1e-3,cond="stem") result_lin_anthropoids$model <- "pure birth with linear variation in speciation rate" #plot_fit_bd(result_lin_anthropoids) ################################################################# # Fit a birth-death model with constant speciation and constant extinction ################################################################# f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){y[1]} lamb_par <- c(0.15) mu_par <-c(0.05) result_bcst_dcst_anthropoids <- fit_bd(anthropoids,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(anthropoids$tip.label)/123,cst.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bcst_dcst_anthropoids$model <- "birth-death with Constant speciation rate and constant extinction" #par(mfrow=c(1,4)) #plot_fit_bd(result_bcst_dcst_anthropoids) ################################################################# # Fit a birth-death model with exponential variation of the speciation # rate with time and constant extinction ################################################################# f.lamb<-function(t,y){y[1] * exp(y[2] * t)} f.mu <-function(t,y){y[1]} lamb_par <- c(0.15, 0.05) mu_par <-c(0.05) result_bexp_dcst_anthropoids <- fit_bd(anthropoids,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(anthropoids$tip.label)/123,expo.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bexp_dcst_anthropoids$model <- "birth-death with exponential variation in speciation rate and constant extinction" #par(mfrow=c(1,3)) #plot_fit_bd(result_bexp_dcst_anthropoids) ################################################################# # Fit a birth-death model with exponential variation (decrease) of the extinction # rate with time and constant speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1]} mu_par <- c(0.03, 0.01) lambda_par <-c(0.07) result_bcst_dexp_anthropoids <- fit_bd(anthropoids,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(anthropoids$tip.label)/123,cst.lamb=TRUE,expo.mu=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bcst_dexp_anthropoids$model <- "birth-death with exponential decrease in extinction rate and constant speciation" #par(mfrow=c(3,3)) #plot_fit_bd(result_bcst_dexp_anthropoids) #plot_dtt(result_bcst_dexp_anthropoids, tot_time, 123) ################################################################# # Fit a birth-death model with exponential variation of the extinction # rate with time and exponential speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1] * exp(y[2] * t)} mu_par <- c(0.05, 0.01) lambda_par <-c(0.15, 0.01) result_bexp_dexp_anthropoids <- fit_bd(anthropoids,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(anthropoids$tip.label)/123,expo.lamb=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bexp_dexp_anthropoids$model <- "birth-death with exponential extinction rate and exponential speciation" #par(mfrow=c(3,1)) #plot_fit_bd(result_bexp_dexp_anthropoids) ################################################################# # Find the best model ################################################################# aiccs_anthropoids <- c(result_cst_anthropoids$aicc, result_exp_anthropoids$aicc, result_lin_anthropoids$aicc,result_bcst_dcst_anthropoids$aicc,result_bexp_dcst_anthropoids$aicc, result_bcst_dexp_anthropoids$aicc, result_bexp_dexp_anthropoids$aicc) aicws_anthropoids<-aicw(aiccs_anthropoids) #write.table(aicws_3, file="anthropoids model comp rpanda.csv", sep=",") ###model averaging #lambda lambdas_anthropoids_tree <- c(result_cst_anthropoids$lamb_par, result_exp_anthropoids$lamb_par[1], result_lin_anthropoids$lamb_par[1],result_bcst_dcst_anthropoids$lamb_par[1], result_bexp_dcst_anthropoids$lamb_par[1], result_bcst_dexp_anthropoids$lamb_par[1], result_bexp_dexp_anthropoids$lamb_par[1]) lambdas_anthropoids_tree<-as.data.frame(lambdas_anthropoids_tree) lambdas_anthropoids_tree lambdas_aicws_anthropoids_tree<-cbind(lambdas_anthropoids_tree, aicws_anthropoids) lambdas_aicws_anthropoids_tree mod.names <- c(result_cst_anthropoids$model, result_exp_anthropoids$model, result_lin_anthropoids$model,result_bcst_dcst_anthropoids$model, result_bexp_dcst_anthropoids$model, result_bcst_dexp_anthropoids$model, result_bexp_dexp_anthropoids$model) rownames(lambdas_aicws_anthropoids_tree)<-mod.names d.aics<-cbind(lambdas_aicws_anthropoids_tree[1,2]-lambdas_aicws_anthropoids_tree[2,2],lambdas_aicws_anthropoids_tree[4,2]-lambdas_aicws_anthropoids_tree[5,2]) colnames(d.aics)<-c("d.aic.pb","d.aic.bd") mod.av.lambda<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.lambda)<-"weighted.lambda" for (i in 1:length(lambdas_aicws_anthropoids_tree$w)){ mod.av.lambda[i,1]<-lambdas_aicws_anthropoids_tree[i,1]*lambdas_aicws_anthropoids_tree[i,4] } mod.av.lambda[,1] sum(mod.av.lambda[,1]) #mu mu_anthropoids_tree <- c(0, 0, 0,result_bcst_dcst_anthropoids$mu_par[1], result_bexp_dcst_anthropoids$mu_par[1], result_bcst_dexp_anthropoids$mu_par[1], result_bexp_dexp_anthropoids$mu_par[1]) mu_anthropoids_tree<-as.data.frame(mu_anthropoids_tree) mu_anthropoids_tree colnames(mu_anthropoids_tree)<-"mu_anthropoids_tree" mu_aicws_anthropoids_tree<-cbind(mu_anthropoids_tree, aicws_anthropoids) mu_aicws_anthropoids_tree mod.av.mu<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.mu)<-"weighted.mu" for (i in 1:length(mu_aicws_anthropoids_tree$w)){ mod.av.mu[i,1]<-mu_aicws_anthropoids_tree[i,1]*mu_aicws_anthropoids_tree[i,4] } mod.av.mu[,1] sum(mod.av.mu[,1]) #time varying parameter tim.var.par_anthropoids_tree <- c(result_cst_anthropoids$lamb_par[2], result_exp_anthropoids$lamb_par[2], result_lin_anthropoids$lamb_par[2],result_bcst_dcst_anthropoids$lamb_par[2], result_bexp_dcst_anthropoids$lamb_par[2], result_bcst_dexp_anthropoids$lamb_par[2], result_bexp_dexp_anthropoids$lamb_par[2]) tim.var.par_anthropoids_tree<-as.data.frame(tim.var.par_anthropoids_tree) tim.var.par_anthropoids_tree tim.var.par_aicws_anthropoids_tree<-cbind(tim.var.par_anthropoids_tree, aicws_anthropoids) tim.var.par_aicws_anthropoids_tree mod.av.tim.var.par<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.tim.var.par)<-"weighted.tim.var.par" for (i in 1:length(tim.var.par_aicws_anthropoids_tree$w)){ mod.av.tim.var.par[i,1]<-tim.var.par_aicws_anthropoids_tree[i,1]*tim.var.par_aicws_anthropoids_tree[i,4] } mod.av.tim.var.par[,1] sum(mod.av.tim.var.par[,1],na.rm=T) lam_mu_timvar_aicws_anthropoids_tree<-cbind(lambdas_anthropoids_tree,mu_anthropoids_tree, tim.var.par_anthropoids_tree, aicws_anthropoids) str(lam_mu_timvar_aicws_anthropoids_tree) lam_mu_timvar_aicws_anthropoids_tree$taxon<-rep("anthropoids",length(lam_mu_timvar_aicws_anthropoids_tree$lambdas_anthropoids_tree)) rownames(lam_mu_timvar_aicws_anthropoids_tree)<-mod.names colnames(lam_mu_timvar_aicws_anthropoids_tree)<-c("lambda","mu","time.var.par","fit","delta","w","taxon") mod_av_res_anthropoids<-cbind(sum(mod.av.lambda[,1]),sum(mod.av.mu[,1]),sum(mod.av.tim.var.par[,1],na.rm=T)) #################################################### ##RPANDA models for OWM #################################################### tot_time<-max(node.age(OWM)$ages)+OWM$root.edge ######################################################################### # Fit the pure birth model (no extinction) with a constant speciation rate ######################################################################### f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){0} lamb_par<-c(0.15) mu_par<-c() result_cst_OWM <- fit_bd(OWM,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(OWM$tip.label)/123,cst.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_cst_OWM$model <- "pure birth with constant speciation rate" #plot_fit_bd(result_cst_OWM, tot_time) ################################################################# # Fit the pure birth model (no extinction) with exponential variation (decrease) # of the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] * exp(y[2] * t)} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.05) mu_par<-c() result_exp_OWM <- fit_bd(OWM,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(OWM$tip.label)/123,expo.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_exp_OWM$model <- "pure birth with exponential variation in speciation rate" #par(mfrow=c(1,3)) #plot_fit_bd(result_exp_OWM, tot_time) ################################################################# # Fit the a pure birth model (no extinction) with linear variation (decrease) of # the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] + y[2] * t} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.001) mu_par<-c() result_lin_OWM <- fit_bd(OWM,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(OWM$tip.label)/123,fix.mu=TRUE,dt=1e-3,cond="stem") result_lin_OWM$model <- "pure birth with linear variation in speciation rate" #plot_fit_bd(result_lin_OWM) ################################################################# # Fit a birth-death model with constant speciation and constant extinction ################################################################# f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){y[1]} lamb_par <- c(0.15) mu_par <-c(0.05) result_bcst_dcst_OWM <- fit_bd(OWM,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(OWM$tip.label)/123,cst.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bcst_dcst_OWM$model <- "birth-death with Constant speciation rate and constant extinction" #par(mfrow=c(1,4)) #plot_fit_bd(result_bcst_dcst_OWM) ################################################################# # Fit a birth-death model with exponential variation of the speciation # rate with time and constant extinction ################################################################# f.lamb<-function(t,y){y[1] * exp(y[2] * t)} f.mu <-function(t,y){y[1]} lamb_par <- c(0.15, 0.05) mu_par <-c(0.05) result_bexp_dcst_OWM <- fit_bd(OWM,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(OWM$tip.label)/123,expo.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bexp_dcst_OWM$model <- "birth-death with exponential variation in speciation rate and constant extinction" #par(mfrow=c(1,3)) #plot_fit_bd(result_bexp_dcst_OWM) ################################################################# # Fit a birth-death model with exponential variation (decrease) of the extinction # rate with time and constant speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1]} mu_par <- c(0.03, 0.01) lambda_par <-c(0.07) result_bcst_dexp_OWM <- fit_bd(OWM,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(OWM$tip.label)/123,cst.lamb=TRUE,expo.mu=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bcst_dexp_OWM$model <- "birth-death with exponential decrease in extinction rate and constant speciation" #par(mfrow=c(3,3)) #plot_fit_bd(result_bcst_dexp_OWM) #plot_dtt(result_bcst_dexp_OWM, tot_time, 123) ################################################################# # Fit a birth-death model with exponential variation of the extinction # rate with time and exponential speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1] * exp(y[2] * t)} mu_par <- c(0.05, 0.01) lambda_par <-c(0.15, 0.01) result_bexp_dexp_OWM <- fit_bd(OWM,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(OWM$tip.label)/123,expo.lamb=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bexp_dexp_OWM$model <- "birth-death with exponential extinction rate and exponential speciation" #par(mfrow=c(3,1)) #plot_fit_bd(result_bexp_dexp_OWM) ################################################################# # Find the best model ################################################################# aiccs_OWM <- c(result_cst_OWM$aicc, result_exp_OWM$aicc, result_lin_OWM$aicc,result_bcst_dcst_OWM$aicc,result_bexp_dcst_OWM$aicc, result_bcst_dexp_OWM$aicc, result_bexp_dexp_OWM$aicc) aicws_OWM<-aicw(aiccs_OWM) #write.table(aicws_3, file="OWM model comp rpanda.csv", sep=",") ###model averaging #lambda lambdas_OWM_tree <- c(result_cst_OWM$lamb_par, result_exp_OWM$lamb_par[1], result_lin_OWM$lamb_par[1],result_bcst_dcst_OWM$lamb_par[1], result_bexp_dcst_OWM$lamb_par[1], result_bcst_dexp_OWM$lamb_par[1], result_bexp_dexp_OWM$lamb_par[1]) lambdas_OWM_tree<-as.data.frame(lambdas_OWM_tree) lambdas_OWM_tree lambdas_aicws_OWM_tree<-cbind(lambdas_OWM_tree, aicws_OWM) lambdas_aicws_OWM_tree mod.names <- c(result_cst_OWM$model, result_exp_OWM$model, result_lin_OWM$model,result_bcst_dcst_OWM$model, result_bexp_dcst_OWM$model, result_bcst_dexp_OWM$model, result_bexp_dexp_OWM$model) rownames(lambdas_aicws_OWM_tree)<-mod.names d.aics<-cbind(lambdas_aicws_OWM_tree[1,2]-lambdas_aicws_OWM_tree[2,2],lambdas_aicws_OWM_tree[4,2]-lambdas_aicws_OWM_tree[5,2]) colnames(d.aics)<-c("d.aic.pb","d.aic.bd") mod.av.lambda<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.lambda)<-"weighted.lambda" for (i in 1:length(lambdas_aicws_OWM_tree$w)){ mod.av.lambda[i,1]<-lambdas_aicws_OWM_tree[i,1]*lambdas_aicws_OWM_tree[i,4] } mod.av.lambda[,1] sum(mod.av.lambda[,1]) #mu mu_OWM_tree <- c(0, 0, 0,result_bcst_dcst_OWM$mu_par[1], result_bexp_dcst_OWM$mu_par[1], result_bcst_dexp_OWM$mu_par[1], result_bexp_dexp_OWM$mu_par[1]) mu_OWM_tree<-as.data.frame(mu_OWM_tree) mu_OWM_tree colnames(mu_OWM_tree)<-"mu_OWM_tree" mu_aicws_OWM_tree<-cbind(mu_OWM_tree, aicws_OWM) mu_aicws_OWM_tree mod.av.mu<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.mu)<-"weighted.mu" for (i in 1:length(mu_aicws_OWM_tree$w)){ mod.av.mu[i,1]<-mu_aicws_OWM_tree[i,1]*mu_aicws_OWM_tree[i,4] } mod.av.mu[,1] sum(mod.av.mu[,1]) #time varying parameter tim.var.par_OWM_tree <- c(result_cst_OWM$lamb_par[2], result_exp_OWM$lamb_par[2], result_lin_OWM$lamb_par[2],result_bcst_dcst_OWM$lamb_par[2], result_bexp_dcst_OWM$lamb_par[2], result_bcst_dexp_OWM$lamb_par[2], result_bexp_dexp_OWM$lamb_par[2]) tim.var.par_OWM_tree<-as.data.frame(tim.var.par_OWM_tree) tim.var.par_OWM_tree tim.var.par_aicws_OWM_tree<-cbind(tim.var.par_OWM_tree, aicws_OWM) tim.var.par_aicws_OWM_tree mod.av.tim.var.par<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.tim.var.par)<-"weighted.tim.var.par" for (i in 1:length(tim.var.par_aicws_OWM_tree$w)){ mod.av.tim.var.par[i,1]<-tim.var.par_aicws_OWM_tree[i,1]*tim.var.par_aicws_OWM_tree[i,4] } mod.av.tim.var.par[,1] sum(mod.av.tim.var.par[,1],na.rm=T) lam_mu_timvar_aicws_OWM_tree<-cbind(lambdas_OWM_tree,mu_OWM_tree, tim.var.par_OWM_tree, aicws_OWM) str(lam_mu_timvar_aicws_OWM_tree) lam_mu_timvar_aicws_OWM_tree$taxon<-rep("OWM",length(lam_mu_timvar_aicws_OWM_tree$lambdas_OWM_tree)) rownames(lam_mu_timvar_aicws_OWM_tree)<-mod.names colnames(lam_mu_timvar_aicws_OWM_tree)<-c("lambda","mu","time.var.par","fit","delta","w","taxon") mod_av_res_OWM<-cbind(sum(mod.av.lambda[,1]),sum(mod.av.mu[,1]),sum(mod.av.tim.var.par[,1],na.rm=T)) #################################################### ##RPANDA models for apes #################################################### tot_time<-max(node.age(apes)$ages)+apes$root.edge ######################################################################### # Fit the pure birth model (no extinction) with a constant speciation rate ######################################################################### f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){0} lamb_par<-c(0.15) mu_par<-c() result_cst_apes <- fit_bd(apes,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(apes$tip.label)/23,cst.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_cst_apes$model <- "pure birth with constant speciation rate" #plot_fit_bd(result_cst_apes, tot_time) ################################################################# # Fit the pure birth model (no extinction) with exponential variation (decrease) # of the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] * exp(y[2] * t)} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.05) mu_par<-c() result_exp_apes <- fit_bd(apes,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(apes$tip.label)/23,expo.lamb=TRUE,fix.mu=TRUE,dt=1e-3,cond="stem") result_exp_apes$model <- "pure birth with exponential variation in speciation rate" #par(mfrow=c(1,3)) #plot_fit_bd(result_exp_apes, tot_time) ################################################################# # Fit the a pure birth model (no extinction) with linear variation (decrease) of # the speciation rate with time ################################################################# f.lamb <-function(t,y){y[1] + y[2] * t} f.mu<-function(t,y){0} lamb_par<-c(0.15, 0.001) mu_par<-c() result_lin_apes <- fit_bd(apes,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(apes$tip.label)/23,fix.mu=TRUE,dt=1e-3,cond="stem") result_lin_apes$model <- "pure birth with linear variation in speciation rate" #plot_fit_bd(result_lin_apes) ################################################################# # Fit a birth-death model with constant speciation and constant extinction ################################################################# f.lamb <-function(t,y){y[1]} f.mu<-function(t,y){y[1]} lamb_par <- c(0.15) mu_par <-c(0.05) result_bcst_dcst_apes <- fit_bd(apes,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(apes$tip.label)/23,cst.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bcst_dcst_apes$model <- "birth-death with Constant speciation rate and constant extinction" #par(mfrow=c(1,4)) #plot_fit_bd(result_bcst_dcst_apes) ################################################################# # Fit a birth-death model with exponential variation of the speciation # rate with time and constant extinction ################################################################# f.lamb<-function(t,y){y[1] * exp(y[2] * t)} f.mu <-function(t,y){y[1]} lamb_par <- c(0.15, 0.05) mu_par <-c(0.05) result_bexp_dcst_apes <- fit_bd(apes,tot_time,f.lamb,f.mu,lamb_par,mu_par, f=length(apes$tip.label)/23,expo.lamb=TRUE,cst.mu=TRUE,dt=1e-3,cond="stem") result_bexp_dcst_apes$model <- "birth-death with exponential variation in speciation rate and constant extinction" #par(mfrow=c(1,3)) #plot_fit_bd(result_bexp_dcst_apes) ################################################################# # Fit a birth-death model with exponential variation (decrease) of the extinction # rate with time and constant speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1]} mu_par <- c(0.03, 0.01) lambda_par <-c(0.07) result_bcst_dexp_apes <- fit_bd(apes,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(apes$tip.label)/23,cst.lamb=TRUE,expo.mu=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bcst_dexp_apes$model <- "birth-death with exponential decrease in extinction rate and constant speciation" #par(mfrow=c(3,3)) #plot_fit_bd(result_bcst_dexp_apes) #plot_dtt(result_bcst_dexp_apes, tot_time, 142) ################################################################# # Fit a birth-death model with exponential variation of the extinction # rate with time and exponential speciation ################################################################# f.mu<-function(t,y){y[1] * exp(y[2] * t)} f.lambda <-function(t,y){y[1] * exp(y[2] * t)} mu_par <- c(0.05, 0.01) lambda_par <-c(0.15, 0.01) result_bexp_dexp_apes <- fit_bd(apes,tot_time,f.lambda,f.mu,lambda_par,mu_par, f=length(apes$tip.label)/23,expo.lamb=TRUE,cst.mu=FALSE,dt=1e-3,cond="stem") result_bexp_dexp_apes$model <- "birth-death with exponential extinction rate and exponential speciation" #par(mfrow=c(3,1)) #plot_fit_bd(result_bexp_dexp_apes) ################################################################# # Find the best model ################################################################# aiccs_apes <- c(result_cst_apes$aicc, result_exp_apes$aicc, result_lin_apes$aicc,result_bcst_dcst_apes$aicc,result_bexp_dcst_apes$aicc, result_bcst_dexp_apes$aicc, result_bexp_dexp_apes$aicc) aicws_apes<-aicw(aiccs_apes) #write.table(aicws_3, file="apes model comp rpanda.csv", sep=",") ###model averaging #lambda lambdas_apes_tree <- c(result_cst_apes$lamb_par, result_exp_apes$lamb_par[1], result_lin_apes$lamb_par[1],result_bcst_dcst_apes$lamb_par[1], result_bexp_dcst_apes$lamb_par[1], result_bcst_dexp_apes$lamb_par[1], result_bexp_dexp_apes$lamb_par[1]) lambdas_apes_tree<-as.data.frame(lambdas_apes_tree) lambdas_apes_tree lambdas_aicws_apes_tree<-cbind(lambdas_apes_tree, aicws_apes) lambdas_aicws_apes_tree mod.names <- c(result_cst_apes$model, result_exp_apes$model, result_lin_apes$model,result_bcst_dcst_apes$model, result_bexp_dcst_apes$model, result_bcst_dexp_apes$model, result_bexp_dexp_apes$model) rownames(lambdas_aicws_apes_tree)<-mod.names d.aics<-cbind(lambdas_aicws_apes_tree[1,2]-lambdas_aicws_apes_tree[2,2],lambdas_aicws_apes_tree[4,2]-lambdas_aicws_apes_tree[5,2]) colnames(d.aics)<-c("d.aic.pb","d.aic.bd") mod.av.lambda<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.lambda)<-"weighted.lambda" for (i in 1:length(lambdas_aicws_apes_tree$w)){ mod.av.lambda[i,1]<-lambdas_aicws_apes_tree[i,1]*lambdas_aicws_apes_tree[i,4] } mod.av.lambda[,1] sum(mod.av.lambda[,1]) #mu mu_apes_tree <- c(0,0,0,result_bcst_dcst_apes$mu_par[1], result_bexp_dcst_apes$mu_par[1], result_bcst_dexp_apes$mu_par[1], result_bexp_dexp_apes$mu_par[1]) mu_apes_tree<-as.data.frame(mu_apes_tree) mu_apes_tree colnames(mu_apes_tree)<-"mu_apes_tree" mu_aicws_apes_tree<-cbind(mu_apes_tree, aicws_apes) mu_aicws_apes_tree mod.av.mu<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.mu)<-"weighted.mu" for (i in 1:length(mu_aicws_apes_tree$w)){ mod.av.mu[i,1]<-mu_aicws_apes_tree[i,1]*mu_aicws_apes_tree[i,4] } mod.av.mu[,1] sum(mod.av.mu[,1]) #time varying parameter tim.var.par_apes_tree <- c(result_cst_apes$lamb_par[2], result_exp_apes$lamb_par[2], result_lin_apes$lamb_par[2],result_bcst_dcst_apes$lamb_par[2], result_bexp_dcst_apes$lamb_par[2], result_bcst_dexp_apes$lamb_par[2], result_bexp_dexp_apes$lamb_par[2]) tim.var.par_apes_tree<-as.data.frame(tim.var.par_apes_tree) tim.var.par_apes_tree tim.var.par_aicws_apes_tree<-cbind(tim.var.par_apes_tree, aicws_apes) tim.var.par_aicws_apes_tree mod.av.tim.var.par<-data.frame(matrix(0,ncol=1,nrow=7)) colnames(mod.av.tim.var.par)<-"weighted.tim.var.par" for (i in 1:length(tim.var.par_aicws_apes_tree$w)){ mod.av.tim.var.par[i,1]<-tim.var.par_aicws_apes_tree[i,1]*tim.var.par_aicws_apes_tree[i,4] } mod.av.tim.var.par[,1] sum(mod.av.tim.var.par[,1],na.rm=T) lam_mu_timvar_aicws_apes_tree<-cbind(lambdas_apes_tree,mu_apes_tree, tim.var.par_apes_tree, aicws_apes) str(lam_mu_timvar_aicws_apes_tree) lam_mu_timvar_aicws_apes_tree$taxon<-rep("apes",length(lam_mu_timvar_aicws_apes_tree$lambdas_apes_tree)) rownames(lam_mu_timvar_aicws_apes_tree)<-mod.names colnames(lam_mu_timvar_aicws_apes_tree)<-c("lambda","mu","time.var.par","fit","delta","w","taxon") mod_av_res_apes<-cbind(sum(mod.av.lambda[,1]),sum(mod.av.mu[,1]),sum(mod.av.tim.var.par[,1],na.rm=T)) mod.avg.params<-rbind.data.frame(mod_av_res_tree,mod_av_res_streps,mod_av_res_tars,mod_av_res_NWM,mod_av_res_anthropoids, mod_av_res_OWM,mod_av_res_apes) mod.avg.params$clade<-c("primates","streps","tars","NWM","anthropoids","OWM","apes") colnames(mod.avg.params)<-c("lambda","mu","time-var.par","clade") mod.avg.params$treenum<-i str(lam_mu_timvar_aicws_tars_tree) mod.fits<-rbind.data.frame(lam_mu_timvar_aicws_tree_tree,lam_mu_timvar_aicws_streps_tree,lam_mu_timvar_aicws_tars_tree,lam_mu_timvar_aicws_NWM_tree,lam_mu_timvar_aicws_anthropoids_tree,lam_mu_timvar_aicws_OWM_tree,lam_mu_timvar_aicws_apes_tree) mod.avg.params.results<-rbind.data.frame(mod.avg.params.results,mod.avg.params) mod.fits.results<-rbind.data.frame(mod.fits.results, mod.fits) } #### look at results mod.avg.params.results str(mod.avg.params.results) std <- function(x) sd(x)/sqrt(length(x)) mod.avg.lam.means<-as.data.frame(t(tapply(mod.avg.params.results$lambda,mod.avg.params.results$clade,FUN=mean))) mod.avg.lam.se<-as.data.frame(t(tapply(mod.avg.params.results$lambda,mod.avg.params.results$clade,FUN=std))) mod.avg.mu.means<-as.data.frame(t(tapply(mod.avg.params.results$mu,mod.avg.params.results$clade,FUN=mean))) mod.avg.mu.se<-as.data.frame(t(tapply(mod.avg.params.results$mu,mod.avg.params.results$clade,FUN=std))) mod.avg.time.mean<-as.data.frame(t(tapply(mod.avg.params.results$`time-var.par`,mod.avg.params.results$clade,FUN=mean))) mod.avg.time.se<-as.data.frame(t(tapply(-mod.avg.params.results$`time-var.par`,mod.avg.params.results$clade,FUN=std))) str(mod.avg.time.mean) mod.fits.results[which(mod.fits.results$taxon=="tree"),] str(mod.fits.results) apes.mod.fits<-mod.fits.results[which(mod.fits.results$taxon=="streps"),] b_prim_mean <- function(t) (mod.avg.lam.means$primates-mod.avg.mu.means$primates)/exp(mod.avg.time.mean$primates*t) # primates b_prim_hci <- function(t) ((mod.avg.lam.means$primates+1.96*mod.avg.lam.se$primates)-(mod.avg.mu.means$primates+1.96*mod.avg.mu.se$primates))/exp((-mod.avg.time.mean$primate+1.96*mod.avg.time.se$primate)*t) # primates b_prim_lci <- function(t) ((mod.avg.lam.means$primates-1.96*mod.avg.lam.se$primates)-(mod.avg.mu.means$primates-1.96*mod.avg.mu.se$primates))/exp((-mod.avg.time.mean$primate-1.96*mod.avg.time.se$primate)*t) # primates b_streps_mean <- function(t) (mod.avg.lam.means$streps-mod.avg.mu.means$streps)/exp(mod.avg.time.mean$streps*t) # streps b_streps_hci <- function(t) ((mod.avg.lam.means$streps+1.96*mod.avg.lam.se$streps)-(mod.avg.mu.means$streps+1.96*mod.avg.mu.se$streps))/exp((-mod.avg.time.mean$primate+1.96*mod.avg.time.se$primate)*t) # streps b_streps_lci <- function(t) ((mod.avg.lam.means$streps-1.96*mod.avg.lam.se$streps)-(mod.avg.mu.means$streps-1.96*mod.avg.mu.se$streps))/exp((-mod.avg.time.mean$primate-1.96*mod.avg.time.se$primate)*t) # streps b_tars_mean <- function(t) (mod.avg.lam.means$tars-mod.avg.mu.means$tars)/exp(mod.avg.time.mean$tars*t) # tars b_tars_hci <- function(t) ((mod.avg.lam.means$tars+1.96*mod.avg.lam.se$tars)-(mod.avg.mu.means$tars+1.96*mod.avg.mu.se$tars))/exp((-mod.avg.time.mean$primate+1.96*mod.avg.time.se$primate)*t) # tars b_tars_lci <- function(t) ((mod.avg.lam.means$tars-1.96*mod.avg.lam.se$tars)-(mod.avg.mu.means$tars-1.96*mod.avg.mu.se$tars))/exp((-mod.avg.time.mean$primate-1.96*mod.avg.time.se$primate)*t) # tars b_OWM_mean <- function(t) (mod.avg.lam.means$OWM)/exp(mod.avg.time.mean$OWM*t) # OWM b_OWM_hci <- function(t) ((mod.avg.lam.means$OWM+1.96*mod.avg.lam.se$OWM))/exp((-mod.avg.time.mean$primate+1.96*mod.avg.time.se$primate)*t) # OWM b_OWM_lci <- function(t) ((mod.avg.lam.means$OWM-1.96*mod.avg.lam.se$OWM))/exp((-mod.avg.time.mean$primate-1.96*mod.avg.time.se$primate)*t) # OWM b_NWM_mean <- function(t) (mod.avg.lam.means$NWM)/exp(mod.avg.time.mean$NWM*t) # NWM b_NWM_hci <- function(t) ((mod.avg.lam.means$NWM+1.96*mod.avg.lam.se$NWM))/exp((-mod.avg.time.mean$primate+1.96*mod.avg.time.se$primate)*t) # NWM b_NWM_lci <- function(t) ((mod.avg.lam.means$NWM-1.96*mod.avg.lam.se$NWM))/exp((-mod.avg.time.mean$primate-1.96*mod.avg.time.se$primate)*t) # NWM b_apes_mean <- function(t) (mod.avg.lam.means$apes)/exp(mod.avg.time.mean$apes*t) # apes b_apes_hci <- function(t) (mod.avg.lam.means$apes+1.96*mod.avg.lam.se$apes)/exp((-mod.avg.time.mean$primate+1.96*mod.avg.time.se$primate)*t) # apes b_apes_lci <- function(t) (mod.avg.lam.means$apes-1.96*mod.avg.lam.se$apes)/exp((-mod.avg.time.mean$primate-1.96*mod.avg.time.se$primate)*t) # apes ###PLot it curve(b_prim_mean,0, 79.48951, xlab = "Time", ylab = "",lwd=2, col="black", las=1,ylim=c(0,1)) curve(b_prim_hci,0, 79.48951,lty=2,col="black",add=T) curve(b_prim_lci,0, 79.48951,lty=2,col="black",add=T) curve(b_streps_mean,0, 71.35462, lwd=2, col="red", add=T) curve(b_streps_hci,0, 71.35462,lty=2,col="red",add=T) curve(b_streps_lci,0, 71.35462,lty=2,col="red",add=T) curve(b_tars_mean,0, 64.17932, lwd=2, col="purple", add=T) curve(b_tars_hci,0, 64.17932,lty=2,col="purple",add=T) curve(b_tars_lci,0, 64.17932,lty=2,col="purple",add=T) curve(b_NWM_mean,0, 37.41346, lwd=2, col="orange", add=T) curve(b_NWM_hci,0, 37.41346,lty=2,col="orange",add=T) curve(b_NWM_lci,0, 37.41346,lty=2,col="orange",add=T) curve(b_OWM_mean,0, 21.28007, lwd=2, col="cyan", add=T) curve(b_OWM_hci,0, 21.28007,lty=2,col="cyan",add=T) curve(b_OWM_lci,0, 21.28007,lty=2,col="cyan",add=T) curve(b_apes_mean,0, 21.27997, lwd=2, col="magenta", add=T) curve(b_apes_hci,0, 21.27997,lty=2,col="magenta",add=T) curve(b_apes_lci,0, 21.27997,lty=2,col="magenta",add=T) legend("topright",c("primates","strepsirrhines","tarsiiforms","platyrhines","OW non-hominoid anthropoids","hominoids"),lty=c(1,1,1,1,1,1),lwd=c(2,2,2,2,2,2),col=c("black","red","purple","orange","cyan","magenta")) write.table(mod.fits.results,file="mod_fit_results.txt",sep="\t") write.table(mod.avg.params.results,file="mod_avg_results.txt",sep="\t") ####################################### ###means plots for figure #######################################3 pdf(file="fig1.pdf",height=6,width=11) par(mfrow=c(1,2)) curve(b_prim_mean,-79.48951,0, xlab = "Time before present (Ma)", ylab = "Diversification rate (n lineages / my)",lwd=2, col="black", las=1,ylim=c(0,1)) curve(b_streps_mean,-71.35462,0, lwd=2, col="red", add=T) curve(b_tars_mean,-64.17932,0, lwd=2, col="purple", add=T) curve(b_NWM_mean,-37.41346,0, lwd=2, col="orange", add=T) curve(b_OWM_mean,-21.28007,0,lwd=2, col="cyan", add=T) curve(b_apes_mean,-21.27997,0, lwd=2, col="magenta", add=T) abline(h=0,lty=2,col="darkred") legend("topleft",c("primates","strepsirrhines","tarsiiforms","platyrrhines","OW non-hominoid anthropoids","hominoids"),lty=c(1,1,1,1,1,1),lwd=c(2,2,2,2,2,2),col=c("black","red","purple","orange","cyan","magenta")) ######################plot the net diversification rates for fossils from PyRate ##all primates allprimates_mean=c(NA, -0.64,-0.637,-0.636,-0.637,-0.619,-0.599,-0.595,-0.593,-0.595,-0.597,-0.595,-0.596,-0.598,-0.596,-0.594,-0.591,-0.486,-0.147,-0.209,-0.092,-0.054,-0.15,-0.231,-0.217,-0.229,-0.201,-0.216,-0.223,-0.21,-0.181,-0.176,-0.181,-0.198,-0.192,-0.123,-0.143,-0.175,-0.2,-0.086,-0.01,0.054,0.065,0.058,0.054,0.053,0.054,0.061,0.068,0.072,0.078,0.117,0.202,0.241,0.252) streps_mean=c(NA, -9.814,-4.63,-3.751,-3.378,-2.854,-2.709,-2.564,-2.289,-2.044,-1.815,-1.632,-1.615,-1.601,-1.564,-1.544,-1.517,-1.492,-1.467,-1.424,-1.404,-1.355,-1.147,-0.443,-0.417,-0.299,-0.284,-0.281,-0.281,-0.281,-0.277,-0.242,-0.243,-0.238,-0.233,-0.188,-0.214,-0.204,-0.179,-0.1,-0.029,0.013,0.043,0.047,0.05,0.051,0.056,0.061,0.066,0.069,0.07,0.085,0.121,0.214,0.216) tars_mean=c(NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,-0.24,-0.239,-0.239,-0.238,-0.239,-0.24,-0.241,-0.242,-0.243,-0.243,-0.244,-0.244,-0.247,-0.244,-0.244,-0.244,-0.232,-0.229,-0.228,-0.218,-0.194,-0.184,-0.183,-0.181,-0.175,-0.152,-0.061,-0.034,-0.022,-0.021,-0.022,-0.022,-0.021,-0.021,-0.02,-0.014,0.016,0.128,0.174,0.184,0.187) platy_mean=c(NA, -7.346,-7.301,-7.184,-7.049,-6.903,-6.608,-6.412,-6.317,-6.119,-6.026,-5.882,-5.774,-5.304,-2.045,-2.83,-0.976,0.124,0.165,0.181,0.187,0.197,0.205,0.18,0.066,-0.064,-0.217,-0.384,-0.528,-0.593,-0.623,-0.67,-0.68,-0.685,-0.687,-0.681,-0.686,-0.688,-0.688,-0.687,-0.687,-0.687) anthro_mean=c(NA, -1.074,-1.072,-1.069,-1.797,-0.485,-0.15,-0.165,-0.093,-0.212,-0.326,-0.334,-0.334,-0.334,-0.333,-0.331,-0.332,-0.335,-0.34,-0.359,-0.406,-0.447,-0.463,-0.49,-0.516,-0.577,-0.613,-0.689,-0.719,-0.352,-0.043,-0.02,-0.015,-0.014,-0.014,-0.009,-0.009,-0.012,-0.01,-0.008,-0.01,-0.026,-0.041,-0.055,-0.069,-0.081,-0.085,-0.085,-0.085,-0.083,-0.082) hominoids_mean=c(NA, -0.812,1.34,-10.249,-11.519,-11.479,-11.436,-11.417,-11.341,-11.347,-11.453,-9.308,-1.732,-1.831,-1.89,-1.91,-1.986,-2.448,-2.264,-2.285,-2.277,-2.659,-2.756,-2.866,-2.96,-3.085,-3.217) trans=0.25 age=(0:(55-1))* -1 plot(age,age,type = 'n', ylim = c(-11.5, 2), xlim = c(-60,0), ylab = 'Diversification rate (n lineages / Ma)', xlab = 'Time (Ma)',las=1 ) abline(h=0,lty=2,col="darkred") lines(rev(age), rev(allprimates_mean), col = "black", lwd=3) lines(rev(age), rev(streps_mean), col = "red", lwd=2) lines(rev(age), rev(tars_mean), col = "purple", lwd=2) age=(0:(42-1))* -1 lines(rev(age), rev(platy_mean), col = "orange", lwd=2) age=(0:(51-1))* -1 lines(rev(age), rev(anthro_mean), col = "cyan", lwd=2) age=(0:(27-1))* -1 lines(rev(age), rev(hominoids_mean), col = "magenta", lwd=2) dev.off() ##################################### #####TESS with all primates ##################################### library(ape); library(geiger);library(phytools) setwd("C:/R2/projects/TESS") primates<-read.tree(file="C:/R2/projects/TESS/springer_trees.tre") par(mfrow=c(1,1)) plot(primates[[1]]) branching.times(primates[[1]]) trees_timescaled<-vector("list") class(trees_timescaled)<-"multiPhylo" for(i in 1:length(primates)){ trees_timescaled[[i]]<-rescale(primates[[i]],model="depth",branching.times(primates[[i]])[1]*100) } trees_timescaled branching.times(trees_timescaled[[1]]) library(TESS) tess1_primates<-tess.analysis(trees_timescaled[[1]],0.15 , 0.05, empiricalHyperPriors=TRUE, empiricalHyperPriorInflation = 10.0, empiricalHyperPriorForm = "lognormal", speciationRatePriorMean = 1.0, speciationRatePriorStDev = 5.0, extinctionRatePriorMean = 1.0, extinctionRatePriorStDev = 5.0, estimateNumberRateChanges = TRUE, numExpectedRateChanges = log(2), samplingProbability = 0.877, tInitialMassExtinction = 34, pInitialMassExtinction = 0.5, pMassExtinctionPriorShape1 = 2, pMassExtinctionPriorShape2 = 10, estimateMassExtinctionTimes = TRUE, numExpectedMassExtinctions = log(2), estimateNumberMassExtinctions = TRUE, MRCA = TRUE, CONDITION = "survival", BURNIN = 10000, MAX_ITERATIONS = 200000, THINNING = 100, OPTIMIZATION_FREQUENCY = 500, CONVERGENCE_FREQUENCY = 1000, MAX_TIME = Inf, MIN_ESS = 200, ADAPTIVE = TRUE, dir = "C:/R2/projects/TESS/allprimates1" , priorOnly = FALSE, verbose = TRUE) tess2_primates<-tess.analysis(trees_timescaled[[1]], 0.15, 0.05, empiricalHyperPriors=TRUE, empiricalHyperPriorInflation = 10.0, empiricalHyperPriorForm = "lognormal", speciationRatePriorMean = 1.0, speciationRatePriorStDev = 5.0, extinctionRatePriorMean = 1.0, extinctionRatePriorStDev = 5.0, estimateNumberRateChanges = TRUE, numExpectedRateChanges = log(2), samplingProbability = 0.877, tInitialMassExtinction = 34, pInitialMassExtinction = 0.5, pMassExtinctionPriorShape1 = 2, pMassExtinctionPriorShape2 = 10, estimateMassExtinctionTimes = TRUE, numExpectedMassExtinctions = log(2), estimateNumberMassExtinctions = TRUE, MRCA = TRUE, CONDITION = "survival", BURNIN = 10000, MAX_ITERATIONS = 200000, THINNING = 100, OPTIMIZATION_FREQUENCY = 500, CONVERGENCE_FREQUENCY = 1000, MAX_TIME = Inf, MIN_ESS = 200, ADAPTIVE = TRUE, dir = "C:/R2/projects/TESS/allprimates1" , priorOnly = FALSE, verbose = TRUE) # Process the output primatesOutput_1 <- tess.process.output(dir="C:/R2/projects/TESS/allprimates1", tree=trees_timescaled[[1]], numExpectedRateChanges=log(2), numExpectedMassExtinctions=log(2)) primatesOutput_2 <- tess.process.output(dir="C:/R2/projects/TESS/allprimates2", tree=trees_timescaled[[1]], numExpectedRateChanges=log(2), numExpectedMassExtinctions=log(2)) # Plot the output outputs_primates <- list(primatesOutput_1,primatesOutput_2) ###chain diagnostics par(mfrow=c(2,4)) tess.plot.multichain.diagnostics(outputs_primates, parameters=c("speciation rates", "speciation shift times", "extinction rates", "extinction shift times", "net-diversification rates", "relative-extinction rates", "mass extinction times"), diagnostics="Gelman-Rubin", gelman.crit=1.05, xlab="million years ago", col=NULL, xaxt="n", yaxt="s", pch=19) ##visualize outputs par(mfrow=c(4,3)) tess.plot.output(outputs_primates[[1]], fig.types=c("speciation rates", "speciation shift times", "speciation Bayes factors", "extinction rates", "extinction shift times", "extinction Bayes factors", "net-diversification rates", "relative-extinction rates", "mass extinction times", "mass extinction Bayes factors"), xlab="million years ago", col=NULL, col.alpha=50, xaxt="n", yaxt="s", pch=19, plot.tree=FALSE) par(mfrow=c(1,1)) tess.plot.output(outputs_primates[[1]], fig.types="speciation rates", xlab="million years ago", col=NULL, col.alpha=50, xaxt="n", yaxt="s", pch=19, plot.tree=FALSE) ############################ ##### compare lineage accumulation through time from phylogeny ##### to fossil species richness through time based on fossil record ########################## ###PALEOTREE library(paleotree) prim.fossils<-read.table(file="all primates.txt",sep="\t",head=T) str(prim.fossils) ##rarefy sp richness in fossils #first create 5 million year time interval bins to calculate the # occurrences of each species in each bin prim.fossils$bin1 <- cut(prim.fossils$base, breaks = c(0, seq(1, 56, by = 5)), labels = 1:12) str(prim.fossils) #calculate the # occurrence of each sp in each bin abund<-t(table(prim.fossils$taxon,prim.fossils$bin1)) head(abund) abund[1:6,1:10] ##raw sp # specnumber(abund) diversity(abund) mean(rowSums(abund)) #rarefy rare<-rarefy(abund,15) rarefied<-as.data.frame(rare) rarefied$int<-as.numeric(rownames(rarefied)) rarefied$intrev<--5*rarefied$int plot(log(rare)~intrev,data=rarefied,pch="",ylab="Rarefied species richness (log)",xlab="Time (Ma)") lines(log(rare)~intrev,data=rarefied) ## the raw primate datafile includes multiple records for species based on the number of occurrence records, use the following to omit duplicated species dup<-duplicated(prim.fossils$taxon) prim.fossils2<-prim.fossils[!dup,] prim.foss<-prim.fossils2 str(prim.foss) summary(prim.foss) rownames(prim.foss)<-prim.foss[,1] prim.foss.mat2<-as.matrix(prim.foss[,c(4,3)]) str(prim.foss.mat2) #plot the diversity curve based on the sampled ranges layout(1:2) logrich<-taxicDivCont(prim.foss.mat2,plotLogRich=T,int.length=1) rawrich<-taxicDivCont(prim.foss.mat2,plotLogRich=F,int.length=1) ###compare to phylogeny phylo_dtt<-phyloDiv(tree, int.length = 1, int.times = NULL, plot = TRUE, plotLogRich = FALSE, drop.ZLB = TRUE, timelims = NULL) str(phylo_dtt) ##simulate missing data for low preservation tax.div.timevar<-c() class(tax.div.timevar)<-"list" tax.div.timeinvar<-c() class(tax.div.timeinvar)<-"list" ###simulate 100 incomplete fossil records with time varying preservation for(i in 1:100){ prim.foss.incompl_timevar<-sampleRanges(prim.foss.mat2, r=0.5, alpha = 10, beta = 10, rTimeRatio = 2, modern.samp.prob = .5, min.taxa = 300, ranges.only = TRUE, minInt = 0.1, merge.cryptic = TRUE, randLiveHat = TRUE, alt.method = FALSE, plot = F) div.cont.timevar<-taxicDivCont(prim.foss.incompl_timevar,plotLogRich=F,int.length=1) tax.div.timevar[[i]]<-div.cont.timevar[,] } tax.div.timevar[[1]] rawdiv1<-tax.div.timevar[[1]] rawdiv2<-rawdiv1[,3] #plot them plot(log(int.div+1)~int.start,data=tax.div.timevar[[1]],type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(logrich[,3])))) for(i in 1:length(tax.div.timevar)){ lines(log(int.div+1)~int.start,data=tax.div.timevar[[i]],col="grey") } lines(log(int.div+1)~int.start,data=logrich,lwd=2) #plot them plot(log(int.div+1)~int.start,data=tax.div.timevar[[1]],type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(logrich[,3])))) for(i in 1:length(tax.div.timevar)){ lines(log(int.div+1)~int.start,data=tax.div.timevar[[i]],col="grey") } lines(log(int.div+1)~int.start,data=logrich,lwd=2) plot(mvalue) ######time invariant for(i in 1:100){ prim.foss.incompl_timeinvar<-sampleRanges(prim.foss.mat2, r=0.5, alpha = 1, beta = 1, rTimeRatio = 1, modern.samp.prob = .5, min.taxa = 300, ranges.only = TRUE, minInt = 0.1, merge.cryptic = TRUE, randLiveHat = TRUE, alt.method = FALSE, plot = F) div.cont.timeinvar<-taxicDivCont(prim.foss.incompl_timeinvar,plotLogRich=F,int.length=1) tax.div.timeinvar[[i]]<-div.cont.timeinvar[,] } #plot them plot(log(int.div+1)~int.start,data=tax.div.timeinvar[[1]],type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(logrich[,3])))) for(i in 1:length(tax.div.timeinvar)){ lines(log(int.div+1)~int.start,data=tax.div.timeinvar[[i]],col="grey80") } for(i in 1:length(tax.div.timevar)){ lines(log(int.div+1)~int.start,data=tax.div.timevar[[i]],col="grey20") } lines(log(int.div+1)~int.start,data=logrich,lwd=2) summary(prim.foss.mat2[,1]-prim.foss.mat2[,2]) layout(1:3) taxicDivCont(prim.foss.mat2,plotLogRich=T,int.length=5) incomplete_timevar<-taxicDivCont(prim.foss.incompl_timevar,plotLogRich = T) taxicDivCont(prim.foss.incompl_timeinvar,plotLogRich = T) ########### ### bootstrap the fossil data to subsample diversity prim.fossils<-as.data.frame(prim.foss.mat2) rownames(prim.fossils) samples<-c() class(samples)<-"list" for(i in 1:1000){ boot.data_rows<-sample(rownames(prim.fossils),size=583,replace=T) boot.data<-prim.fossils[boot.data_rows,] dups<-duplicated(rownames(boot.data)) boot.div<-taxicDivCont(boot.data,plotLogRich = F) samples[[i]]<-boot.div } length(samples) #plot them plot(log(int.div+1)~int.start,data=tax.div.timeinvar[[1]],type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(logrich[,3])))) for(i in 1:length(tax.div.timeinvar)){ lines(log(int.div+1)~int.start,data=tax.div.timeinvar[[i]],col="grey90") } for(i in 1:length(tax.div.timevar)){ lines(log(int.div+1)~int.start,data=tax.div.timevar[[i]],col="grey80") } for(i in 1:100){ lines(log(int.div+1)~int.start,data=samples[[i]],col="grey50",lty=2) } lines(log(int.div+1)~int.start,data=logrich,lwd=3) legend("topright",c("primates","strepsirrhines","tarsiiforms","platyrhines","OW non-hominoid anthropoids","hominoids"),lty=c(1,1,1,1,1,1),lwd=c(2,2,2,2,2,2),col=c("black","red","purple","orange","cyan","magenta")) ###compare to phylogeny tax.div.tree<-c() class(tax.div.tree)<-"list" for(i in 1:4){ phylo_dtt<-phyloDiv(trees[[i]], int.length = 1, int.times = NULL, plot = TRUE, plotLogRich = FALSE, drop.ZLB = TRUE, timelims = NULL) tax.div.tree[[i]]<-phylo_dtt } tax.div.tree1<-tax.div.tree[[1]] ###plot all together pdf(file="diversity through time.pdf",height=8,width=10) rev.int.div<-rev(tax.div.tree1[,3]) rev.int.start<-rev(tax.div.tree1[,1]) rev.int.start<-rev.int.start*-1 plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(rev.int.div))),las=1) ##sim time invar for(i in 1:length(tax.div.timeinvar)){ tax.div.timeinvar1<-tax.div.timeinvar[[i]] rev.int.div<-rev(tax.div.timeinvar1[,3]) rev.int.start<-rev(tax.div.timeinvar1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey80") } ##sim time var for(i in 1:length(tax.div.timevar)){ tax.div.timevar1<-tax.div.timevar[[i]] rev.int.div<-rev(tax.div.timevar1[,3]) rev.int.start<-rev(tax.div.timevar1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60") } ##boot for(i in 1:100){ sample1<-samples[[i]] rev.int.div<-rev(sample1[,3]) rev.int.start<-rev(sample1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey20",lty=2) } ##observed fossil rev.int.div<-rev(logrich[,3]) rev.int.start<-rev(logrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3) ##trees for(i in 1:4){ tax.div.tree1<-tax.div.tree[[i]] rev.int.div <-rev(tax.div.tree1[,3]) rev.int.start<-rev(tax.div.tree1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lty=4,lwd=3) } legend("topleft",c("Observed fossil","Bootstrapped fossil","Observed phylogeny","Simulated fossil, \ntime invariant preservation","Simulated fossil, \ntime varying preservation"), lty=c(1,2,4,1,1),lwd=c(4,2,4,2,2,2),col=c("black","grey20","black","grey80","grey60")) dev.off() ###calculate mean and confidence intervals on ltt from extant using ltt95 from phytools ltt_CI<-ltt95(trees,alpha=0.05, log=T, method="lineages",mode="mean") ltt_CI<-multiDiv(trees,int.length = 1, plot = TRUE, split.int = TRUE, drop.ZLB = TRUE, drop.cryptic = FALSE, extant.adjust = 0.01, plotLogRich = T, timelims = NULL, int.times = NULL, plotMultCurves = FALSE, multRainbow = TRUE, divPalette = NULL, divLineType = 1, main = NULL) ltt_CI$median.curve ltt.divs<-as.data.frame(ltt_CI$median.curve) ltt.times<-as.data.frame(ltt_CI$int.times) str(ltt.times) multlttres<-cbind(ltt.times$int.start,ltt.divs) str(multlttres) colnames(multlttres)<-c("time","median","low","high") plot(multlttres$median~-1*(multlttres$time)) ###plot without subsampled preservation bias pdf(file="diversity through time2.pdf",height=5,width=8) rev.int.div<-rev(tax.div.tree1[,3]) rev.int.start<-rev(tax.div.tree1[,1]) rev.int.start<-rev.int.start*-1 plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(tax.div.tree1[,3]))),las=1) ##boot for(i in 1:1000){ samples1<-samples[[i]] rev.int.div<-rev(samples1[,3]) rev.int.start<-rev(samples1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60",lty=2) } ##observed fossil rev.int.div<-rev(logrich[,3]) rev.int.start<-rev(logrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3) ##trees for(i in 1:4){ tax.div.tree1<-tax.div.tree[[i]] rev.int.div<-rev(tax.div.tree1[,3]) rev.int.start<-rev(tax.div.tree1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lty=4,lwd=3) } legend("topleft",c("Observed fossil","Bootstrapped fossil","Observed phylogeny"), lty=c(1,2,4),lwd=c(4,2,4),col=c("black","grey60","black")) dev.off() ###plot without subsampled preservation bias and with CIs for phylo instead of raw plots pdf(file="diversity through time2.pdf",height=5,width=9) tax.div.tree1<-multlttres[,2] rev.int.div<-rev(tax.div.tree1) rev.int.start<-rev(multlttres[,1]) rev.int.start<-rev.int.start*-1 plot(log(rev.int.div)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(tax.div.tree1))),las=1) ##boot for(i in 1:1000){ samples1<-samples[[i]] rev.int.div<-rev(samples1[,3]) rev.int.start<-rev(samples1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60",lty=2) } ##observed fossil rev.int.div<-rev(logrich[,3]) rev.int.start<-rev(logrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3) ##trees tax.div.tree1<-multlttres[,2] rev.int.div<-rev(tax.div.tree1) rev.int.start<-rev(multlttres[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lty=4,lwd=3) tax.div.tree1<-multlttres[,3] rev.int.div<-rev(tax.div.tree1) rev.int.start<-rev(multlttres[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey80",lty=4,lwd=3) tax.div.tree1<-multlttres[,4] rev.int.div<-rev(tax.div.tree1) rev.int.start<-rev(multlttres[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey80",lty=4,lwd=3) legend("topleft",c("Observed fossil","Bootstrapped fossil","Extant phylogeny","Phylogeny 95% quantiles"), lty=c(1,2,4,4),lwd=c(4,2,4,4),col=c("black","grey60","black","grey80")) dev.off() ###sensitivity of patterns to interval length ########### ### bootstrap the fossil data to subsample diversity prim.fossils<-as.data.frame(prim.foss.mat2) rownames(prim.fossils) summary(prim.foss.mat2[,1]-prim.foss.mat2[,2]) samples.int1<-c() class(samples.int1)<-"list" for(i in 1:1000){ boot.data_rows<-sample(rownames(prim.fossils),size=583,replace=T) boot.data<-prim.fossils[boot.data_rows,] boot.div.int<-taxicDivCont(boot.data,int.length=1,plotLogRich = F) samples.int1[[i]]<-boot.div.int } samples.int2<-c() class(samples.int2)<-"list" for(i in 1:1000){ boot.data_rows<-sample(rownames(prim.fossils),size=583,replace=T) boot.data<-prim.fossils[boot.data_rows,] boot.div.int<-taxicDivCont(boot.data,int.length=2,plotLogRich = F) samples.int2[[i]]<-boot.div.int } samples.int3<-c() class(samples.int3)<-"list" for(i in 1:1000){ boot.data_rows<-sample(rownames(prim.fossils),size=583,replace=T) boot.data<-prim.fossils[boot.data_rows,] boot.div.int<-taxicDivCont(boot.data,int.length=3,plotLogRich = F) samples.int3[[i]]<-boot.div.int } samples.int4<-c() class(samples.int4)<-"list" for(i in 1:1000){ boot.data_rows<-sample(rownames(prim.fossils),size=583,replace=T) boot.data<-prim.fossils[boot.data_rows,] boot.div.int<-taxicDivCont(boot.data,int.length=4,plotLogRich = F) samples.int4[[i]]<-boot.div.int } samples.int5<-c() class(samples.int5)<-"list" for(i in 1:1000){ boot.data_rows<-sample(rownames(prim.fossils),size=583,replace=T) boot.data<-prim.fossils[boot.data_rows,] boot.div.int<-taxicDivCont(boot.data,int.length=5,plotLogRich = F) samples.int5[[i]]<-boot.div.int } ###observed obs.int1<-taxicDivCont(prim.foss.mat2,int.length=1,plotLogRich = F) obs.int2<-taxicDivCont(prim.foss.mat2,int.length=2,plotLogRich = F) obs.int3<-taxicDivCont(prim.foss.mat2,int.length=3,plotLogRich = F) obs.int4<-taxicDivCont(prim.foss.mat2,int.length=4,plotLogRich = F) obs.int5<-taxicDivCont(prim.foss.mat2,int.length=5,plotLogRich = F) ####plot pdf(file="sensitivity of diversity patterns to interval length.pdf",height=11,width=8) layout(1:5) plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",xlim=c(-70,0),ylim=c(0,max(log(tax.div.tree1[,3]))),las=1,main="Diversity with interval length = 1 My") ##boot 1 for(i in 1:1000){ samples1<-samples.int1[[i]] rev.int.div<-rev(samples1[,3]) rev.int.start<-rev(samples1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60",lty=2) } ##observed fossil rev.int.div<-rev(obs.int1[,3]) rev.int.start<-rev(obs.int1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lwd=1) plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",xlim=c(-70,0),ylim=c(0,max(log(tax.div.tree1[,3]))),las=1,main="Diversity with interval length = 2 My") ##boot 2 for(i in 1:1000){ samples1<-samples.int2[[i]] rev.int.div<-rev(samples1[,3]) rev.int.start<-rev(samples1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60",lty=2) } ##observed fossil rev.int.div<-rev(obs.int2[,3]) rev.int.start<-rev(obs.int2[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lwd=1) plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",xlim=c(-70,0),ylim=c(0,max(log(tax.div.tree1[,3]))),las=1,main="Diversity with interval length = 3 My") ##boot 3 for(i in 1:1000){ samples1<-samples.int3[[i]] rev.int.div<-rev(samples1[,3]) rev.int.start<-rev(samples1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60",lty=2) } ##observed fossil rev.int.div<-rev(obs.int3[,3]) rev.int.start<-rev(obs.int3[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lwd=1) plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",xlim=c(-70,0),ylim=c(0,max(log(tax.div.tree1[,3]))),las=1,main="Diversity with interval length = 4 My") ##boot 4 for(i in 1:1000){ samples1<-samples.int4[[i]] rev.int.div<-rev(samples1[,3]) rev.int.start<-rev(samples1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60",lty=2) } ##observed fossil rev.int.div<-rev(obs.int4[,3]) rev.int.start<-rev(obs.int4[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lwd=1) plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",xlim=c(-70,0),ylim=c(0,max(log(tax.div.tree1[,3]))),las=1,main="Diversity with interval length = 5 My") ##boot 5 for(i in 1:1000){ samples1<-samples.int5[[i]] rev.int.div<-rev(samples1[,3]) rev.int.start<-rev(samples1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60",lty=2) } ##observed fossil rev.int.div<-rev(obs.int5[,3]) rev.int.start<-rev(obs.int5[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lwd=1) dev.off() ###################### ## LTTs/DTTs for each clade separately ########################### ###hominoids hominoid.fossils<-read.table(file="hominoids.txt",sep="\t",head=T) str(hominoid.fossils) dup<-duplicated(hominoid.fossils$taxon) hominoid.fossils2<-hominoid.fossils[!dup,] hominoid.foss<-hominoid.fossils2 str(hominoid.foss) summary(hominoid.foss) rownames(hominoid.foss)<-hominoid.foss[,1] summary(hominoid.foss[,4]-hominoid.foss[,3]) #prim.fossils<-read.table(file="data for paleotree.txt",sep="\t",head=T) hominoid.foss.mat2<-as.matrix(hominoid.foss[,c(4,3)]) #rownames(prim.foss.mat2)<-prim.foss[,1] str(hominoid.foss.mat2) #plot the diversity curve based on the sampled ranges hominoidlogrich<-taxicDivCont(hominoid.foss.mat2,plotLogRich=T,int.length=1) ########### ### bootstrap the fossil data to subsample diversity hominoid.fossils<-as.data.frame(hominoid.foss.mat2) rownames(hominoid.fossils) hominoid.samples<-c() class(hominoid.samples)<-"list" for(i in 1:1000){ boot.data_rows<-sample(rownames(hominoid.fossils),size=95,replace=T) boot.data<-hominoid.fossils[boot.data_rows,] boot.div<-taxicDivCont(boot.data,plotLogRich = F) hominoid.samples[[i]]<-boot.div } length(hominoid.samples) #plot them plot(log(int.div+1)~int.start,data=hominoidlogrich,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(hominoidlogrich[,3])+0.5)),las=1) for(i in 1:1000){ lines(log(int.div+1)~int.start,data=hominoid.samples[[i]],col="grey50",lty=2) } lines(log(int.div+1)~int.start,data=hominoidlogrich,lwd=3) ###compare to phylogeny tax.div.ape.tree<-c() class(tax.div.ape.tree)<-"list" for(i in 1:4){ apes_node<-findMRCA(trees[[i]],tips=c("Pan_paniscus","Hylobates_lar")) apes<-extract.clade(trees[[i]], apes_node,root.edge = 1) apes phylo_dtt<-phyloDiv(apes, int.length = 1, int.times = NULL, plot = F, plotLogRich = FALSE, drop.ZLB = TRUE, timelims = NULL) tax.div.ape.tree[[i]]<-phylo_dtt } tax.div.ape.tree1<-tax.div.ape.tree[[1]] ###plot without subsampled preservation bias pdf(file="ape diversity through time2.pdf",height=5,width=8) rev.int.div<-rev(tax.div.ape.tree1[,3]) rev.int.start<-rev(tax.div.ape.tree1[,1]) rev.int.start<-rev.int.start*-1 plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(tax.div.ape.tree1[,3]+1)+0.5)),las=1) ##boot for(i in 1:1000){ samples1<-hominoid.samples[[i]] rev.int.div<-rev(samples1[,3]) rev.int.start<-rev(samples1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60",lty=2) } ##observed fossil rev.int.div<-rev(hominoidlogrich[,3]) rev.int.start<-rev(hominoidlogrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3) ##trees for(i in 1:4){ tax.div.tree1<-tax.div.ape.tree[[i]] rev.int.div<-rev(tax.div.tree1[,3]) rev.int.start<-rev(tax.div.tree1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lty=4,lwd=3) } legend("topleft",c("Observed fossil","Bootstrapped fossil","Observed phylogeny"), lty=c(1,2,4),lwd=c(4,2,4),col=c("black","grey60","black")) dev.off() ###streps streps.fossils<-read.table(file="strepsirrhines.txt",sep="\t",head=T) str(streps.fossils) dup<-duplicated(streps.fossils$taxon) streps.fossils2<-streps.fossils[!dup,] streps.foss<-streps.fossils2 str(streps.foss) summary(streps.foss) rownames(streps.foss)<-streps.foss[,1] summary(streps.foss[,4]-streps.foss[,3]) #prim.fossils<-read.table(file="data for paleotree.txt",sep="\t",head=T) streps.foss.mat2<-as.matrix(streps.foss[,c(4,3)]) #rownames(prim.foss.mat2)<-prim.foss[,1] str(streps.foss.mat2) #plot the diversity curve based on the sampled ranges strepslogrich<-taxicDivCont(streps.foss.mat2,plotLogRich=T,int.length=1) ########### ### bootstrap the fossil data to subsample diversity streps.fossils<-as.data.frame(streps.foss.mat2) rownames(streps.fossils) streps.samples<-c() class(streps.samples)<-"list" for(i in 1:1000){ boot.data_rows<-sample(rownames(streps.fossils),size=161,replace=T) boot.data<-streps.fossils[boot.data_rows,] boot.div<-taxicDivCont(boot.data,plotLogRich = F) streps.samples[[i]]<-boot.div } length(streps.samples) #plot them plot(log(int.div+1)~int.start,data=strepslogrich,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(strepslogrich[,3])+0.5)),las=1) for(i in 1:1000){ lines(log(int.div+1)~int.start,data=streps.samples[[i]],col="grey50",lty=2) } lines(log(int.div+1)~int.start,data=strepslogrich,lwd=3) ###compare to phylogeny tax.div.streps.tree<-c() class(tax.div.streps.tree)<-"list" for(i in 1:4){ strepss_node<-findMRCA(trees[[i]],tips=c("Loris_tardigradus","Avahi_laniger")) strepss<-extract.clade(trees[[i]], strepss_node,root.edge = 1) strepss phylo_dtt<-phyloDiv(strepss, int.length = 1, int.times = NULL, plot = F, plotLogRich = FALSE, drop.ZLB = TRUE, timelims = NULL) tax.div.streps.tree[[i]]<-phylo_dtt } tax.div.streps.tree1<-tax.div.streps.tree[[1]] ###plot without subsampled preservation bias pdf(file="streps diversity through time2.pdf",height=5,width=8) rev.int.div<-rev(tax.div.streps.tree1[,3]) rev.int.start<-rev(tax.div.streps.tree1[,1]) rev.int.start<-rev.int.start*-1 plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(tax.div.streps.tree1[,3]+1)+0.5)),las=1) ##boot for(i in 1:1000){ samples1<-streps.samples[[i]] rev.int.div<-rev(samples1[,3]) rev.int.start<-rev(samples1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60",lty=2) } ##observed fossil rev.int.div<-rev(strepslogrich[,3]) rev.int.start<-rev(strepslogrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3) ##trees for(i in 1:4){ tax.div.tree1<-tax.div.streps.tree[[i]] rev.int.div<-rev(tax.div.tree1[,3]) rev.int.start<-rev(tax.div.tree1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lty=4,lwd=3) } legend("topleft",c("Observed fossil","Bootstrapped fossil","Observed phylogeny"), lty=c(1,2,4),lwd=c(4,2,4),col=c("black","grey60","black")) dev.off() ###tars tars.fossils<-read.table(file="tarsiiformes.txt",sep="\t",head=T) str(tars.fossils) dup<-duplicated(tars.fossils$taxon) tars.fossils2<-tars.fossils[!dup,] tars.foss<-tars.fossils2 str(tars.foss) summary(tars.foss) rownames(tars.foss)<-tars.foss[,1] summary(tars.foss[,4]-tars.foss[,3]) #prim.fossils<-read.table(file="data for paleotree.txt",sep="\t",head=T) tars.foss.mat2<-as.matrix(tars.foss[,c(4,3)]) #rownames(prim.foss.mat2)<-prim.foss[,1] str(tars.foss.mat2) #plot the diversity curve based on the sampled ranges tarslogrich<-taxicDivCont(tars.foss.mat2,plotLogRich=T,int.length=1) ########### ### bootstrap the fossil data to subsample diversity tars.fossils<-as.data.frame(tars.foss.mat2) rownames(tars.fossils) tars.samples<-c() class(tars.samples)<-"list" for(i in 1:1000){ boot.data_rows<-sample(rownames(tars.fossils),size=121,replace=T) boot.data<-tars.fossils[boot.data_rows,] boot.div<-taxicDivCont(boot.data,plotLogRich = F) tars.samples[[i]]<-boot.div } length(tars.samples) #plot them plot(log(int.div+1)~int.start,data=tarslogrich,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(tarslogrich[,3])+0.5)),las=1) for(i in 1:1000){ lines(log(int.div+1)~int.start,data=tars.samples[[i]],col="grey50",lty=2) } lines(log(int.div+1)~int.start,data=tarslogrich,lwd=3) ###compare to phylogeny tax.div.tars.tree<-c() class(tax.div.tars.tree)<-"list" for(i in 1:4){ tarss_node<-findMRCA(trees[[i]],tips=c("Tarsius_syrichta","Tarsius_tarsier")) tarss<-extract.clade(trees[[i]], tarss_node,root.edge = 1) tarss phylo_dtt<-phyloDiv(tarss, int.length = 1, int.times = NULL, plot = F, plotLogRich = FALSE, drop.ZLB = TRUE, timelims = NULL) tax.div.tars.tree[[i]]<-phylo_dtt } tax.div.tars.tree1<-tax.div.tars.tree[[1]] ###plot without subsampled preservation bias pdf(file="tars diversity through time2.pdf",height=5,width=8) rev.int.div<-rev(tarslogrich[,3]) rev.int.start<-rev(tarslogrich[,1]) rev.int.start<-rev.int.start*-1 plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(tarslogrich[,3]+1)+0.5)),las=1) ##boot for(i in 1:1000){ samples1<-tars.samples[[i]] rev.int.div<-rev(samples1[,3]) rev.int.start<-rev(samples1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60",lty=2) } ##observed fossil rev.int.div<-rev(tarslogrich[,3]) rev.int.start<-rev(tarslogrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3) ##trees for(i in 1:4){ tax.div.tree1<-tax.div.tars.tree[[i]] rev.int.div<-rev(tax.div.tree1[,3]) rev.int.start<-rev(tax.div.tree1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lty=4,lwd=3) } legend("topright",c("Observed fossil","Bootstrapped fossil","Observed phylogeny"), lty=c(1,2,4),lwd=c(4,2,4),col=c("black","grey60","black")) dev.off() ###platyrhine platyrhine.fossils<-read.table(file="platyrhines_species.txt",sep="\t",head=T) str(platyrhine.fossils) dup<-duplicated(platyrhine.fossils$taxon) platyrhine.fossils2<-platyrhine.fossils[!dup,] platyrhine.foss<-platyrhine.fossils2 str(platyrhine.foss) summary(platyrhine.foss) rownames(platyrhine.foss)<-platyrhine.foss[,1] summary(platyrhine.foss[,4]-platyrhine.foss[,3]) #prim.fossils<-read.table(file="data for paleotree.txt",sep="\t",head=T) platyrhine.foss.mat2<-as.matrix(platyrhine.foss[,c(4,3)]) #rownames(prim.foss.mat2)<-prim.foss[,1] str(platyrhine.foss.mat2) #plot the diversity curve based on the sampled ranges platyrhinelogrich<-taxicDivCont(platyrhine.foss.mat2,plotLogRich=T,int.length=1) ########### ### bootstrap the fossil data to subsample diversity platyrhine.fossils<-as.data.frame(platyrhine.foss.mat2) rownames(platyrhine.fossils) summary(platyrhine.fossils) platyrhine.samples<-c() class(platyrhine.samples)<-"list" for(i in 1:1000){ boot.data_rows<-sample(rownames(platyrhine.fossils),size=39,replace=T) boot.data<-platyrhine.fossils[boot.data_rows,] boot.div<-taxicDivCont(boot.data,plotLogRich = F) platyrhine.samples[[i]]<-boot.div } length(platyrhine.samples) #plot them plot(log(int.div+1)~int.start,data=platyrhinelogrich,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(platyrhinelogrich[,3])+0.5)),las=1) for(i in 1:1000){ lines(log(int.div+1)~int.start,data=platyrhine.samples[[i]],col="grey50",lty=2) } lines(log(int.div+1)~int.start,data=platyrhinelogrich,lwd=3) ###compare to phylogeny tax.div.platyrhine.tree<-c() class(tax.div.platyrhine.tree)<-"list" for(i in 1:4){ platyrhines_node<-findMRCA(trees[[i]],tips=c("Pithecia_pithecia","Cebus_apella")) platyrhines<-extract.clade(trees[[i]], platyrhines_node,root.edge = 1) platyrhines phylo_dtt<-phyloDiv(platyrhines, int.length = 1, int.times = NULL, plot = F, plotLogRich = FALSE, drop.ZLB = TRUE, timelims = NULL) tax.div.platyrhine.tree[[i]]<-phylo_dtt } tax.div.platyrhine.tree1<-tax.div.platyrhine.tree[[1]] layout(1:2) ###plot without subsampled preservation bias pdf(file="platyrhine diversity through time2.pdf",height=5,width=8) rev.int.div<-rev(platyrhinelogrich[,3]) rev.int.start<-rev(platyrhinelogrich[,1]) rev.int.start<-rev.int.start*-1 plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(tax.div.platyrhine.tree1[,3]+1)+0.5)),las=1) ##boot for(i in 1:1000){ samples1<-platyrhine.samples[[i]] rev.int.div<-rev(samples1[,3]) rev.int.start<-rev(samples1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60",lty=2) } ##observed fossil rev.int.div<-rev(platyrhinelogrich[,3]) rev.int.start<-rev(platyrhinelogrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3) ##trees for(i in 1:4){ tax.div.tree1<-tax.div.platyrhine.tree[[i]] rev.int.div<-rev(tax.div.tree1[,3]) rev.int.start<-rev(tax.div.tree1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lty=4,lwd=3) } legend("topleft",c("Observed fossil","Bootstrapped fossil","Observed phylogeny"), lty=c(1,2,4),lwd=c(4,2,4),col=c("black","grey60","black")) dev.off() ###owanthro owanthro.fossils<-read.table(file="catarhines sans hominoids.txt",sep="\t",head=T) str(owanthro.fossils) dup<-duplicated(owanthro.fossils$taxon) owanthro.fossils2<-owanthro.fossils[!dup,] owanthro.foss<-owanthro.fossils2 str(owanthro.foss) summary(owanthro.foss) rownames(owanthro.foss)<-owanthro.foss[,1] summary(owanthro.foss[,4]-owanthro.foss[,3]) #prim.fossils<-read.table(file="data for paleotree.txt",sep="\t",head=T) owanthro.foss.mat2<-as.matrix(owanthro.foss[,c(4,3)]) #rownames(prim.foss.mat2)<-prim.foss[,1] str(owanthro.foss.mat2) #plot the diversity curve based on the sampled ranges owanthrologrich<-taxicDivCont(owanthro.foss.mat2,plotLogRich=T,int.length=1) ########### ### bootstrap the fossil data to subsample diversity owanthro.fossils<-as.data.frame(owanthro.foss.mat2) rownames(owanthro.fossils) summary(owanthro.fossils) owanthro.samples<-c() class(owanthro.samples)<-"list" for(i in 1:1000){ boot.data_rows<-sample(rownames(owanthro.fossils),size=160,replace=T) boot.data<-owanthro.fossils[boot.data_rows,] boot.div<-taxicDivCont(boot.data,plotLogRich = F) owanthro.samples[[i]]<-boot.div } length(owanthro.samples) #plot them plot(log(int.div+1)~int.start,data=owanthrologrich,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(owanthrologrich[,3])+0.5)),las=1) for(i in 1:1000){ lines(log(int.div+1)~int.start,data=owanthro.samples[[i]],col="grey50",lty=2) } lines(log(int.div+1)~int.start,data=owanthrologrich,lwd=3) ###compare to phylogeny tax.div.owanthro.tree<-c() class(tax.div.owanthro.tree)<-"list" for(i in 1:4){ owanthros_node<-findMRCA(trees[[i]],tips=c("Papio_anubis","Pan_troglodytes")) owanthros<-extract.clade(trees[[i]], owanthros_node,root.edge = 1) owanthros phylo_dtt<-phyloDiv(owanthros, int.length = 1, int.times = NULL, plot = F, plotLogRich = FALSE, drop.ZLB = TRUE, timelims = NULL) tax.div.owanthro.tree[[i]]<-phylo_dtt } tax.div.owanthro.tree1<-tax.div.owanthro.tree[[1]] layout(1:2) ###plot without subsampled preservation bias pdf(file="owanthro diversity through time2.pdf",height=5,width=8) rev.int.div<-rev(owanthrologrich[,3]) rev.int.start<-rev(owanthrologrich[,1]) rev.int.start<-rev.int.start*-1 plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,max(log(tax.div.owanthro.tree1[,3]+1)+0.5)),las=1) ##boot for(i in 1:1000){ samples1<-owanthro.samples[[i]] rev.int.div<-rev(samples1[,3]) rev.int.start<-rev(samples1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="grey60",lty=2) } ##observed fossil rev.int.div<-rev(owanthrologrich[,3]) rev.int.start<-rev(owanthrologrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3) ##trees for(i in 1:4){ tax.div.tree1<-tax.div.owanthro.tree[[i]] rev.int.div<-rev(tax.div.tree1[,3]) rev.int.start<-rev(tax.div.tree1[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="black",lty=4,lwd=3) } legend("topleft",c("Observed fossil","Bootstrapped fossil","Observed phylogeny"), lty=c(1,2,4),lwd=c(4,2,4),col=c("black","grey60","black")) dev.off() ##### diversity of all primate clades pdf(file="primate clades diversity compared.pdf",height=5,width=8) rev.int.div<-rev(logrich[,3]) rev.int.start<-rev(logrich[,1]) rev.int.start<-rev.int.start*-1 plot(log(rev.int.div+1)~rev.int.start,type="n",ylab="Species richness (log)",xlab="Time (Ma)",ylim=c(0,5),las=1) ##observed fossil rev.int.div<-rev(logrich[,3]) rev.int.start<-rev(logrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3) rev.int.div<-rev(strepslogrich[,3]) rev.int.start<-rev(strepslogrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3,col="red") rev.int.div<-rev(tarslogrich[,3]) rev.int.start<-rev(tarslogrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3,col="purple") rev.int.div<-rev(platyrhinelogrich[,3]) rev.int.start<-rev(platyrhinelogrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3,col="orange") rev.int.div<-rev(hominoidlogrich[,3]) rev.int.start<-rev(hominoidlogrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,lwd=3,col="magenta") rev.int.div<-rev(owanthrologrich[,3]) rev.int.start<-rev(owanthrologrich[,1]) rev.int.start<-rev.int.start*-1 lines(log(rev.int.div+1)~rev.int.start,col="cyan",lwd=3) dev.off() legend("topright",c("primates","strepsirrhines","tarsiiforms","platyrrhines","OW non-hominoid anthropoids","hominoids"),lty=c(1,1,1,1,1,1),lwd=c(2,2,2,2,2,2),col=c("black","red","purple","orange","cyan","magenta")) ################################ ###prepare the data for pyrate ################################ source(file = "pyrate_utilities.r") extract.ages(file = "all primates.txt", replicates=10, cutoff=NULL, random=TRUE) extract.ages(file = "tarsiiformes.txt", replicates=10, cutoff=NULL, random=TRUE) extract.ages(file = "strepsirrhines_sp_logmass.txt", replicates=10, cutoff=NULL, random=T) extract.ages(file = "hominoids.txt", replicates=10, cutoff=NULL, random=FALSE) extract.ages(file = "platyrhines_sp_logmass.txt", replicates=10, cutoff=NULL, random=T) extract.ages(file = "catarhines_sans_hominoids.txt", replicates=10, cutoff=NULL, random=T) extract.ages(file = "anthropoids.txt", replicates=10, cutoff=NULL, random=T) ### example of code for MCMC diagnostics library(coda);library(HDInterval) prim<-read.table(file="/Users/jamesherrera/Documents/PyRate-master/pyrate_mcmc_logs/primates_mcmc.log",sep="\t",head=T) str(prim[,1:13]) plot(prim[,2]~prim[,1]) effectiveSize(prim[,2]) median(prim[,2]) hdi(prim[,2]) plot(prim[,3]~prim[,1]) effectiveSize(prim[,3]) plot(prim[,4]~prim[,1]) effectiveSize(prim[,4]) plot(prim[,5]~prim[,1]) effectiveSize(prim[,5]) plot(prim[,6]~prim[,1]) effectiveSize(prim[,6]) median(prim[,6]) hdi(prim[,6]) plot(prim[,7]~prim[,1]) effectiveSize(prim[,7]) median(prim[,7]) hdi(prim[,7])