# Tanja Stadler, Feb 3, 2015 # Some examples to simulate phylogenies and estimate speciation and extinction rates library(TreeSim) library(TreePar) set.seed(10) ################################################################### # Simulation example - time-dependent rates # # First we simulate a tree: # Number of species nspecies <- 100 # At time 1 in the past, we have a rate shift: time <- c(0,1) #=(t0,t1=t2) from slide 11 in Macroevolution lecture # Half of the present day species are sampled (rho[1]=0.5): rho <- c(0.5,1) # speciation rates (between t[i],t[i+1] we have speciation rate lambda[i]): lambda <- c(2,5) #=(lambda0,lambda1=lambda2) from slide 11 in Macroevolution lecture # extinction rates (between t[i],t[i+1] we have extinction rate mu[i]): mu <- c(1.5,0) #=(mu0,mu1=mu2) from slide 11 in Macroevolution lecture # Simulation of a tree: tree<-sim.rateshift.taxa(nspecies,1,lambda,mu,frac=rho,times=time,complete=FALSE) # Extracting the speciation times x: x<-sort(getx(tree[[1]]),decreasing=TRUE) # When estimating the rate shift times t based on branching times x, # we allow the shift times to be 0.6, 0.8, 1, 1.2, .. ,2.4: start <- 0.4 end <- 2 grid <- 0.2 # We fix rho and estimate time, lambda, mu: res <- bd.shifts.optim(x,c(rho,1),grid,start,end)[[2]] res # res[[2]] tells us about the maximum likelihood estimate given one rate shift: # - log lik = 60.6940763. # rate shift at time 1.0. # turnover (extinction/speciation) = 0.68 more recent than 1.0, # and = 0.19 more ancestral than 1.0. # net diversification (speciation-extinction) rate = 0.81 more recent than 1.0, # and = 3.66 more ancestral than 1.0. #values used for simulation: mu/lambda lambda-mu #test if 1 shift explain the tree significantly better than 0 shifts: #if test>0.95 then 1 shift is significantly better than 0 shifts at a 5% error i<-1 test<-pchisq(2*(res[[i]][1]-res[[i+1]][1]),3) test #test if 2 shifts explain the tree significantly better than 1 shift: i<-2 test<-pchisq(2*(res[[i]][1]-res[[i+1]][1]),3) test # Plot parameter estimates bd.shifts.plot(list(res),1,2.1,0,5) # Plot parameters used for simulation lines(c(-max(x),-time[2]),c((lambda-mu)[2],(lambda-mu)[2]),col="red") lines(c(-time[2],0),c((lambda-mu)[1],(lambda-mu)[1]),col="red") ################################################################### # Simulation example - Test for Diversity-Dependent speciation set.seed(1) tree<-sim.rateshift.taxa(10,1,c(2,0.1),c(0,0.05),frac=c(1,1),times=time,complete=FALSE) x<-sort(getx(tree[[1]]),decreasing=TRUE) # The tree size is too small to make reliable estimations! We just analyse this tree as an examplefor getting answers quickly! resDD<-bd.densdep.optim(x,discrete=FALSE,continuous=TRUE)[[2]] # Slow! Package expoTree on CRAN performs same calculations much faster!! Same method also within Package DDD. resShifts <- bd.shifts.optim(x,c(rho,1),0.1,0.1,2.1)[[2]][[2]] resDD resShifts # Best model where AIC smallest AICDD <- 2*3+2*resDD\$value AICShifts <- 2*5+2*resShifts[1] AICDD AICShifts ################################################################### # Data example - bird orders # # Data analysis in Stadler & Bokma, 2013: # Number of species in each order from Sibley and Monroe (1990) data(bird.orders) plot(bird.orders) # Many species are missing - ie only one per order S <- c(10, 47, 69, 214, 161, 17, 355, 51, 56, 10, 39, 152, 6, 143, 358, 103, 319, 23, 291, 313, 196, 1027, 5712) # number of species per order: cbind(bird.orders\$tip.label,S) groupscut<-get.groups(bird.orders,S,96.43) #all orders established at 96.43 Ma x<-branching.times(bird.orders) # transforming molecular timescale into calendar timescale x<-x/0.207407 # allowing one shift in rates: resbirds<-bd.shifts.optim(x,sampling=c(1,1),grid=1,start=96,end=135,survival=1,groups=groupscut)[[2]] i<-1 test<-pchisq(2*(res[[i]][1]-res[[i+1]][1]),3) test # Sampling model is important: resbirdsWrongSamp<-bd.shifts.optim(x,sampling=c(length(S)/sum(S),1),grid=1,start=96,end=135,survival=1,groups=0)[[2]] resbirdsCompleteSamp<-bd.shifts.optim(x,sampling=c(1,1),grid=1,start=96,end=135,survival=1,groups=0)[[2]] # Check the difference (especially look at turnover estimates!): resbirds resbirdsWrongSamp resbirdsCompleteSamp