##### DATING USING CAL3 ##### # Load in geological stages Ranges<-read.table("Ranges.txt", header=TRUE, row.names=1) # Load in First and Last Appearance dates taxa<-read.table("Taxa.txt", header=TRUE, row.names=1) # Get into the right format data<-list(Ranges, taxa) # ML optimisation for sampling probability and extinction rates SPres<-getSampProbDisc(data) # Establish means for interval length intLength<--apply(data[[1]],1,diff) meanInt<-mean(intLength) # Calculate sampling rate sRate<-sProb2sRate(SPres[[2]][2], int.length=meanInt) # Calculate diversification and extinction rate divRate<-SPres[[2]][1]/meanInt # Load trees for cal3 trees<-read.nexus("trees.nex") # Get consensus of trees contree<-consensus(trees) # Date trees with cal3 datedtrees<-bin_cal3TimePaleoPhy(CFO, data, brRate=divRate, extRate=divRate, sampRate=sRate, ntrees=50, plot=FALSE) ##### FITTING OF EVOLUTIONARY MODELS ##### # Read in raw masses data<-read.table("Masses.txt", header=F, row.names=1) # If necessary, read dated trees from other analysis # datedtrees<-read.nexus("datedtrees.nex") # Log the data data<-log(data) # Create dummy matrix weights<-matrix(ncol=50, nrow=7) # Fit models to each dated tree for (i in 1:50) { all.data<-treedata(datedtrees[[i]], data)$data bm.fit <- fitContinuous_paleo(datedtrees[[i]], all.data[,1], model = "BM",meserr = 0.15) ou.fit <- fitContinuous_paleo(datedtrees[[i]], all.data[,1], model = "OU",meserr = 0.15) white.fit <- fitContinuous_paleo(datedtrees[[i]], all.data[,1], model = "white",meserr = 0.15) trend.fit<- fitContinuous_paleo(datedtrees[[i]], all.data[,1], model = "trend",meserr = 0.15) shift.fit <- fitContinuous_paleo(datedtrees[[i]], all.data[,1], model = "timeshift",meserr = 0.15, shift.time = 66) release.fit <- fitContinuous_paleo(datedtrees[[i]], all.data[,1], model = "release",meserr = 0.15, shift.time = 66) releaseradiate.fit <- fitContinuous_paleo(datedtrees[[i]], all.data[,1], model = "releaseradiate",meserr = 0.15, shift.time = 66) aic.c <- c(bm.fit$Trait1$aicc, trend.fit$Trait1$aicc, ou.fit$Trait1$aicc, white.fit$Trait1$aicc, shift.fit$Trait1$aicc, release.fit$Trait1$aicc, releaseradiate.fit$Trait1$aicc); names(aic.c) <- c("Brownian\nMotion", "Directional \nTrend","Ornstein \nUhlenbeck", "White \nNoise", "K-Pg \nshift", "K-Pg \nrelease", "K-Pg release \n& radiate"); wts<-aicw(aic.c) weights[,i]<-wts[,3] }