Script from Armin showing that "terminal branches" are being included in dating of trees

Cal3_minMax_problem.R

library(paleotree)


#Simulate some fossil ranges with simFossilRecord
set.seed(444)
record<-simFossilRecord(p=0.1, q=0.1, nruns=1,
                        nTotalTaxa=c(30,40), nExtant=0)
taxa<-fossilRecord2fossilTaxa(record)
#simulate a fossil record with imperfect sampling with sampleRanges
rangesCont <- sampleRanges(taxa,r=0.5)
#let's use taxa2cladogram to get the 'ideal' cladogram of the taxa
cladogram <- taxa2cladogram(taxa,plot=TRUE)
#this library allows one to use rate calibrated type time-scaling methods (Bapst, in prep.)
#to use these, we need an estimate of the sampling rate (we set it to 0.5 above)
likFun<-make_durationFreqCont(rangesCont)
srRes<-optim(parInit(likFun),likFun,lower=parLower(likFun),upper=parUpper(likFun),
             method="L-BFGS-B",control=list(maxit=1000000))
sRate <- srRes[[1]][2]
# we also need extinction rate and branching rate
# we can get extRate from getSampRateCont too
#we'll assume extRate=brRate (ala Foote et al., 1999); may not always be a good assumption
divRate<-srRes[[1]][1]


# ####Basic timescaling####
# #timePaleoPhy method using "minMax"
ttree.minmax <- timePaleoPhy(cladogram,rangesCont, 
                             dateTreatment="minMax",
                             ntrees=1,plot=TRUE)
# #Just by looking at the plot we already see that the tips do not necessarily extend up to their LAD
# #We can check this properly

# #Find node heights above the root for all tips
mynodeheights <- -(node.depth.edgelength(ttree.minmax) - ttree.minmax$root.time)[1:Ntip(ttree.minmax)]
names(mynodeheights) <- ttree.minmax$tip.label
mynodeheights

# #We exclude NAs from our rangeCont object and focus only on the LADs
rangesCont.LAD.only <- rangesCont[,2][complete.cases(rangesCont[,2])]

# #Let's match the order of mynodeheights with that of rangesCont.LAD.only
mynodeheights.ordered <- mynodeheights[match(names(rangesCont.LAD.only), names(mynodeheights))]

# #Compare our initial LADs and the LAD in our timescaled phylogeny
rangesCont.LAD.only
mynodeheights.ordered
all.equal(rangesCont.LAD.only, mynodeheights.ordered)
# #They are clearly different!



# ####Cal3 timescaling####
# #Now let's repeat this approach for cal3
# #cal3TimePaleoPhy method using "minMax" (note that we exclude ancestors)
ttree.cal3.minmax <- cal3TimePaleoPhy(cladogram,rangesCont,brRate=divRate,extRate=divRate,
                                 sampRate=sRate,dateTreatment="minMax",
                                 ntrees=1,anc.wt=0,plot=TRUE)




# #Find node heights above the root for all tips
mynodeheights.cal3 <- -(node.depth.edgelength(ttree.cal3.minmax) - ttree.cal3.minmax$root.time)[1:Ntip(ttree.cal3.minmax)]
names(mynodeheights.cal3) <- ttree.cal3.minmax$tip.label
mynodeheights.cal3


# #Let's match again the order of mynodeheights.cal3 with that of rangesCont.LAD.only
mynodeheights.cal3.ordered <- mynodeheights.cal3[match(names(rangesCont.LAD.only), names(mynodeheights.cal3))]

# #Compare our initial LADs and the LAD in our timescaled phylogeny
rangesCont.LAD.only
mynodeheights.cal3.ordered
all.equal(rangesCont.LAD.only, mynodeheights.cal3.ordered)
# #They are THE SAME!