dwbapst / paleotree

R library for analyzing, time-scaling and simulating phylogenies of extinct/fossil lineages. Also plots diversity curves for stratigraphic range data and phylogenies, including combinations of these two data types.
Creative Commons Zero v1.0 Universal
21 stars 9 forks source link

cal3TimePaleoPhy not properly adjusting dates according to dateTreatment = minMax #19

Closed dwbapst closed 5 years ago

dwbapst commented 5 years ago

The following issue was reported via email in October 2018 by Armin Elsler.

The functions timePaleoPhy and cal3TimePaleoPhy allows for a method of defining dates as two-date minMax. This means the ages given are interpreted as minMax limits on a single occurrence. The limits themselves should not affect the tree as calculated. However, cal3TimePaleoPhy does not seem to be doing that, when minMax is used.

First, note that argument add.term isn't used, if option minMax is used with timePaleoPhy. The two options are incompatible - there's no terminal range with minMax, because each taxon is a single occurrence. However, add.term is effectively always true with cal3TimePaleoPh - terminal ranges are always added. But what is the 'terminal range' that cal3 is adding?

This leads to several questions, I have followed these questions with my expectations of how those situations should be resolved in paleotree, without checking which are legal, or if they work that way at present.

1) What is the behavior of cal3TimePaleoPhy when each taxon is supposed to be a single occurrence with a precise age? (This would be the default usage with firstLast, with precise first and last appearance times being identical.)

As Armin reports:

So, one option would be to use timePaleoPhy with the 'minMax' option. However, we cannot use the cal3 implementation of timePaleoPhy since it will automatically add the terminal range of the taxa (default option of cal3). If we use FAD.only = TRUE it will restrict the taxon to the lower bound of its oldest occurrence. So cal3timePaleoPhy with the 'minMax' option should (?) not be the correct choice (timePaleoPhy, on the other hand, would be fine).

In the above bug report, my understanding is that the 'terminal branch' is basically the duration of time between the randomly selected occurrence time, and the maximum age bound - these two are being interpreted as the taxon's first and last appearance times, and are equivalent to treating minMax as being identical to the randObs option.

Except, maybe it isn't?

More from Armin:

I have attached some piece of code (based on the examples you give in the help file of the cal3 method) that shows, that the "minMax" argument together with cal3timePaleoPhy will always give you the LAD for the tips. timePaleoPhy, on the other hand, will not (and therefore works as expected).

The issue that remains, however, is that "minMax" together with cal3timePaleoPhy will always (!) add the terminal branches. This is also mentioned in your help file:

"Also, cal3 will always add the terminal ranges of taxa."

So, the LAD of the "minMax" + cal3 option on the phylogeny will always coincide with the LAD as given in the "timeData" object. This is not the case for the standard timePaleoPhy function.

So, in this case, it sounds like minMax is doing something very different from it being an accidental clone of the randObs algorithm. Instead, the minMax option is pulling an age from between the minimum and maximum, this is being treated as a first occurrence time (which is sensible, for the sake of properly dating nodes on the tree, unlike randObs, where it is sensible to use the oldest age as the bound for determining branch ages). The problem is that the upper bound is then being treated as a precisely known last appearance time, and the remaining duration between the randomly selected first-appearance time and the latest possible age is treated as a 'terminal branch`.

Obviously, this is highly undesirable.

Finally, from Armin:

cal3timePaleoPhy definitely behaves differently compared to TimePaleoPhy when using the argument "minMax": Unlike TimePaleoPhy it will always (!) add the terminal ranges to the tips. Is this behavior intended?

No, it is not, and if this is true, it is extremely worrying.

Analysis of this issue, beginning with the file mentioned by Armin, and pinning down which scenarios violate the expectations listed above will continue below.

dwbapst commented 5 years ago

https://gist.github.com/dwbapst/909cc3c2e9ae8086d6030f047f1b3ed7

dwbapst commented 5 years ago

The following is output from R of running a script testing the 5 expectations given above.

> # testing issues with minMax date treatment
> 
> library(paleotree)
Loading required package: ape
> 
> test_identical_tip_ages <- function(trees){
+ # a function to test if tip ages are different
+ #
+ # checks
+ if(!inherits(trees,"multiPhylo")){
+ stop("trees is not class multiphylo")
+ }
+ if(length(trees)<2){
+ stop("fewer than two trees given")
+ }
+ ###########
+ # get tip ages
+ node_dates <- lapply(trees,dateNodes)
+ tip_dates <- lapply(1:length(node_dates),function(i){
+ x <- node_dates[[i]]
+ x <- x[1:Ntip(trees[[i]])]
+ names(x) <- trees[[i]]$tip.label
+ x <- x[order(names(x))]
+ return(x)
+ })
+ # test if tip ages are nearly the same
+ all_nearly_identical <- all(sapply(tip_dates,
+ function(x) isTRUE(all.equal(x,tip_dates[[1]]))
+ ))
+ return(all_nearly_identical)
+ }
> 
> ################################
> 
> #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)
> 
> # set up a set of ranges with identical FADs and LADs
> # just repeat FADs
> rangesCont_pointocc <- rangesCont
> rangesCont_pointocc[,2] <- rangesCont[,1]
> 
> #this library allows one to use
> # rate calibrated type time-scaling methods (Bapst, 2014.)
> #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]
> ##################################################################################
> # 1) What is the behavior of `cal3TimePaleoPhy` when each taxon is supposed
> # to be a single occurrence with a precise age? 
> # This would be the default usage with `firstLast`, with precise
> # first and last appearance times being identical (`rangesCont_pointocc`),
> # and `FAD.only = FALSE`.
> # - Expectation: Tip ages shouldn't change from one tree to another,
> # as there is no tip age uncertainty expressed.
> 
> # #cal3TimePaleoPhy method using "firstLast" 
> # ancestors excluded
> # FAD.only = FALSE
> # FADs and LADs are identical (point occurrences)
> ttrees_1 <- cal3TimePaleoPhy(
+ cladogram, rangesCont_pointocc,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="firstLast",
+ FAD.only=FALSE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Warning: Following taxa dropped from tree: t2, t6, t24, t30, t31
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_1)
[1] TRUE
> 
> # Result
> # Expectation met, treatment 'firstLast' has no impact on resulting tip ages.
> 
> ##################################################################################
> # 2) Is this behavior identical to `cal3TimePaleoPhy` with `minMax`,
> # using identical minimum and maximum ages?
> # This would be the default usage with `firstLast`, with precise
> # first and last appearance times being identical(`rangesCont_pointocc`),
> # and `FAD.only = FALSE`.
> # - Expectation: Tip ages shouldn't change from one tree to another,
> # as (again) there is no tip age uncertainty expressed.
> 
> # #cal3TimePaleoPhy method using "minMax" 
> # ancestors excluded
> # FAD.only = FALSE
> # FADs and LADs are identical (point occurrences)
> ttrees_2<- cal3TimePaleoPhy(
+ cladogram, rangesCont_pointocc,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="minMax",
+ FAD.only=FALSE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Warning: Following taxa dropped from tree: t2, t6, t24, t30, t31
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_2)
[1] TRUE
> 
> # RESULT
> # Expectation met, treatment 'minMax' has no impact on resulting tip ages.
> #################################################################################
> # 3a) What is the difference in behavior of 
> # `cal3TimePaleoPhy` with `firstLast` versus `minMax`?
> # when the oldest and youngest ages given for each taxon are non-identical?
> # - Expectation: Tip ages for `firstLast` should always be equal
> # to the given last appearance times, with no variation across analyses.
> 
> # #cal3TimePaleoPhy method using "firstLast" 
> # ancestors excluded
> # FAD.only = FALSE
> # FADs and LADs can differ
> ttrees_3a <- cal3TimePaleoPhy(
+ cladogram, rangesCont,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="firstLast",
+ FAD.only=FALSE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Warning: Following taxa dropped from tree: t2, t6, t24, t30, t31
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_3a)
[1] TRUE
> 
> # Result
> # Expectation met, treatment 'firstLast' has no impact on resulting tip ages.
> #################################################################################
> # 3b) What is the difference in behavior of 
> # `cal3TimePaleoPhy` with `firstLast` versus `minMax`?
> # when the oldest and youngest ages given for each taxon are non-identical?
> # - Expectation: Tip ages for `minMax` should be somewhere between the given
> # minimum and maximum bounds, with variation across runs.
> 
> # #cal3TimePaleoPhy method using "minMax" 
> # ancestors excluded
> # FAD.only = FALSE
> # FADs and LADs can differ
> ttrees_3b <- cal3TimePaleoPhy(
+ cladogram, rangesCont,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="minMax",
+ FAD.only=FALSE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Warning: Following taxa dropped from tree: t2, t6, t24, t30, t31
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_3b)
[1] TRUE
> 
> # RESULT
> # Tip ages do not differ. Problematic.
> ################################################################################
> # 4a) When `FAD.only = TRUE`, what is the difference in behavior
> #  of `cal3TimePaleoPhy` with `firstLast` versus `minMax`?
> # With FADs and LADs *not* being identical.
> # - Expectation: Tip ages for `firstLast` should always be equal
> # to the given first appearance times, with no variation across analyses.
> 
> # #cal3TimePaleoPhy method using "firstLast" 
> # ancestors excluded
> # FAD.only = TRUE
> # FADs and LADs can differ
> ttrees_4a <- cal3TimePaleoPhy(
+ cladogram, rangesCont,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="firstLast",
+ FAD.only=TRUE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Warning: Following taxa dropped from tree: t2, t6, t24, t30, t31
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_4a)
[1] TRUE
> 
> # Result
> # Expectation met, treatment 'firstLast' has no impact on resulting tip ages.
> 
> ################################################################################
> # 4b) When `FAD.only = TRUE`, what is the difference in behavior
> #  of `cal3TimePaleoPhy` with `firstLast` versus `minMax`?
> # With FADs and LADs *not* being identical.
> # - Expectation: Tip ages for `minMax` should be somewhere between the given
> # minimum and maximum bounds, with variation across runs.
> 
> # #cal3TimePaleoPhy method using "minMax" 
> # ancestors excluded
> # FAD.only = TRUE
> # FADs and LADs can differ
> ttrees_4b<- cal3TimePaleoPhy(
+ cladogram, rangesCont,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="minMax",
+ FAD.only=TRUE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Error in cal3TimePaleoPhy(cladogram, rangesCont, brRate = divRate, extRate = divRate,  : 
  FAD.only = TRUE and dateTreatment = 'minMax' are conflicting, as there are no FADs, as dates are simply point occurrences
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_4b)
Error in test_identical_tip_ages(trees = ttrees_4b) : 
  object 'ttrees_4b' not found
> 
> # RESULT
> # Returns an error, as it finds minMax incongruent with FAD.only = TRUE.
> 
> ###########################################################################
> # 5a) What is the difference in behavior with `randObs` between 
> # `FAD.only = FALSE` and `FAD.only = TRUE` ?
> # With FADs and LADs *not* being identical.
> # - Expectation: When `FAD.only` is `FALSE`, tip ages from `randObs`
> # should be somewhere between the given minimum and maximum bounds,
> # with variation across runs, much like `minMax`. 
> 
> # #cal3TimePaleoPhy method using "randObs" 
> # ancestors excluded
> # FAD.only = FALSE
> # FADs and LADs can differ
> ttrees_5a<- cal3TimePaleoPhy(
+ cladogram, rangesCont,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="randObs",
+ FAD.only=FALSE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Warning: Following taxa dropped from tree: t2, t6, t24, t30, t31
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_5a)
[1] FALSE
> 
> # RESULT
> # Expectation met, there is variation as expected.
> 
> ###########################################################################
> # 5b) What is the difference in behavior with `randObs` between 
> # `FAD.only = FALSE` and `FAD.only = TRUE` ? 
> # With FADs and LADs *not* being identical.
> # - Expectation: When `FAD.only` is `TRUE`, tip agesfrom `randObs`
> # should be equal to the given first appearance times, with no variation
> # across runs, as with `firstLast`.
> 
> # #cal3TimePaleoPhy method using "randObs" 
> # ancestors excluded
> # FAD.only = TRUE
> # FADs and LADs can differ
> ttrees_5b<- cal3TimePaleoPhy(
+ cladogram, rangesCont,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="randObs",
+ FAD.only=TRUE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Error in cal3TimePaleoPhy(cladogram, rangesCont, brRate = divRate, extRate = divRate,  : 
  FAD.only = TRUE and dateTreatment = 'randObs' are conflicting arguments
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_5b)
Error in test_identical_tip_ages(trees = ttrees_5b) : 
  object 'ttrees_5b' not found
> 
> # RESULT
> # Returns an error, as it finds randObs incongruent with FAD.only = TRUE.
dwbapst commented 5 years ago

The expectations fail with 3b, in that dateTreatment = 'minMax' is treating the age given in the second column as a precisely known LAD, not an upper bound, and treating the point occurrence pulled from the range as a FAD.

So, yes, Armin was exactly correct with his bug report.

dwbapst commented 5 years ago

Note that randObs is not similar to minMax - in fact its nearly the opposite scenario from the current. randObs will treat the FADs as precisely known, and treat the LAD as a free floating random variable within the bounds of the ranges given . minMax instead treats the FAD as freely floating, affecting when branching events can occur, while treating LADs as fixed.

dwbapst commented 5 years ago

As of v3.2.7:

> # testing consistency issues with date treatment under cal3 (05-31-19)
> # making sure date treatment issues as reported by Armin Essler in October 2018
> # are avoided in current code base 
> 
> library(paleotree)
Loading required package: ape
> 
> ###########################################################
> 
> test_identical_tip_ages <- function(trees){
+ # a function to test if tip ages are different
+ #
+ # checks
+ if(!inherits(trees,"multiPhylo")){
+ stop("trees is not class multiphylo")
+ }
+ if(length(trees)<2){
+ stop("fewer than two trees given")
+ }
+ ###########
+ # get tip ages
+ node_dates <- lapply(trees,dateNodes)
+ tip_dates <- lapply(1:length(node_dates),function(i){
+ x <- node_dates[[i]]
+ x <- x[1:Ntip(trees[[i]])]
+ names(x) <- trees[[i]]$tip.label
+ x <- x[order(names(x))]
+ return(x)
+ })
+ # test if tip ages are nearly the same
+ all_nearly_identical <- all(sapply(tip_dates,
+ function(x) isTRUE(all.equal(x,tip_dates[[1]]))
+ ))
+ return(all_nearly_identical)
+ }
> 
> ################################
> 
> #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)
> 
> # set up a set of ranges with identical FADs and LADs
> # just repeat FADs
> rangesCont_pointocc <- rangesCont
> rangesCont_pointocc[,2] <- rangesCont[,1]
> 
> 
> #this library allows one to use
> # rate calibrated type time-scaling methods (Bapst, 2014.)
> #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]
>  ##################################################################################
> # 1) What is the behavior of `cal3TimePaleoPhy` when each taxon is supposed
> # to be a single occurrence with a precise age? 
> # This would be the default usage with `firstLast`, with precise
> # first and last appearance times being identical (`rangesCont_pointocc`),
> # and `FAD.only = FALSE`.
> # - Expectation: Tip ages shouldn't change from one tree to another,
> # as there is no tip age uncertainty expressed.
> 
> # #cal3TimePaleoPhy method using "firstLast" 
> # ancestors excluded
> # FAD.only = FALSE
> # FADs and LADs are identical (point occurrences)
> ttrees_1 <- cal3TimePaleoPhy(
+ cladogram, rangesCont_pointocc,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="firstLast",
+ FAD.only=FALSE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Warning: Following taxa dropped from tree:
 t2, t6, t24, t30, t31
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_1)
[1] TRUE
> 
> # Result
> # Expectation met, treatment 'firstLast' has no impact on resulting tip ages.
> 
#############################################################################
> # 2) Is this behavior identical to `cal3TimePaleoPhy` with `minMax`,
> # using identical minimum and maximum ages?
> # This would be the default usage with `firstLast`, with precise
> # first and last appearance times being identical(`rangesCont_pointocc`),
> # and `FAD.only = FALSE`.
> # - Expectation: Tip ages shouldn't change from one tree to another,
> # as (again) there is no tip age uncertainty expressed.
> 
> # #cal3TimePaleoPhy method using "minMax" 
> # ancestors excluded
> # FAD.only = FALSE
> # FADs and LADs are identical (point occurrences)
> ttrees_2<- cal3TimePaleoPhy(
+ cladogram, rangesCont_pointocc,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="minMax",
+ FAD.only=FALSE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Warning: Following taxa dropped from tree:
 t2, t6, t24, t30, t31
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_2)
[1] TRUE
> 
> # RESULT
> # Expectation met, treatment 'minMax' has no impact on resulting tip ages.
> 
#############################################################################
> # 3a) What is the difference in behavior of 
> # `cal3TimePaleoPhy` with `firstLast` versus `minMax`?
> # when the oldest and youngest ages given for each taxon are non-identical?
> # - Expectation: Tip ages for `firstLast` should always be equal
> # to the given last appearance times, with no variation across analyses.
> 
> # #cal3TimePaleoPhy method using "firstLast" 
> # ancestors excluded
> # FAD.only = FALSE
> # FADs and LADs can differ
> ttrees_3a <- cal3TimePaleoPhy(
+ cladogram, rangesCont,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="firstLast",
+ FAD.only=FALSE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Warning: Following taxa dropped from tree:
 t2, t6, t24, t30, t31
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_3a)
[1] TRUE
> 
> # Result
> # Expectation met, treatment 'firstLast' has no impact on resulting tip ages.
> 
#############################################################################
> # 3b) What is the difference in behavior of 
> # `cal3TimePaleoPhy` with `firstLast` versus `minMax`?
> # when the oldest and youngest ages given for each taxon are non-identical?
> # - Expectation: Tip ages for `minMax` should be somewhere between the given
> # minimum and maximum bounds, with variation across runs.
> 
> # #cal3TimePaleoPhy method using "minMax" 
> # ancestors excluded
> # FAD.only = FALSE
> # FADs and LADs can differ
> ttrees_3b <- cal3TimePaleoPhy(
+ cladogram, rangesCont,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="minMax",
+ FAD.only=FALSE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Warning: Following taxa dropped from tree:
 t2, t6, t24, t30, t31
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_3b)
[1] FALSE
> 
> # RESULT
> # Tip ages DO differ among the trees, as expected. Huzzah!
> 
############################################################################
> # 4a) When `FAD.only = TRUE`, what is the difference in behavior
> #  of `cal3TimePaleoPhy` with `firstLast` versus `minMax`?
> # With FADs and LADs *not* being identical.
> # - Expectation: Tip ages for `firstLast` should always be equal
> # to the given first appearance times, with no variation across analyses.
> 
> # #cal3TimePaleoPhy method using "firstLast" 
> # ancestors excluded
> # FAD.only = TRUE
> # FADs and LADs can differ
> ttrees_4a <- cal3TimePaleoPhy(
+ cladogram, rangesCont,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="firstLast",
+ FAD.only=TRUE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Warning: Following taxa dropped from tree:
 t2, t6, t24, t30, t31
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_4a)
[1] TRUE
> 
> # Result
> # Expectation met, treatment 'firstLast' has no impact on resulting tip ages.
> 
############################################################################
> # 4b) When `FAD.only = TRUE`, what is the difference in behavior
> #  of `cal3TimePaleoPhy` with `firstLast` versus `minMax`?
> # With FADs and LADs *not* being identical.
> # - Expectation: Tip ages for `minMax` should be somewhere between the given
> # minimum and maximum bounds, with variation across runs.
> 
> # #cal3TimePaleoPhy method using "minMax" 
> # ancestors excluded
> # FAD.only = TRUE
> # FADs and LADs can differ
> ttrees_4b<- cal3TimePaleoPhy(
+ cladogram, rangesCont,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="minMax",
+ FAD.only=TRUE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Error in cal3TimePaleoPhy(cladogram, rangesCont, brRate = divRate, extRate = divRate,  : 
  FAD.only = TRUE is inapplicable with dateTreatment = 'minMax' or 'randObs'
'minMax' assumes dates of observation of tips are point occurrences
   And thus taxa are not considered to have ranges that have FADs
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_4b)
Error in test_identical_tip_ages(trees = ttrees_4b) : 
  object 'ttrees_4b' not found
> 
> # RESULT
> # Returns an error, as it finds minMax incongruent with FAD.only = TRUE.
> 
> ###########################################################################
> # 5a) What is the difference in behavior with `randObs` between 
> # `FAD.only = FALSE` and `FAD.only = TRUE` ?
> # With FADs and LADs *not* being identical.
> # - Expectation: When `FAD.only` is `FALSE`, tip ages from `randObs`
> # should be somewhere between the given minimum and maximum bounds,
> # with variation across runs, much like `minMax`. 
> 
> # #cal3TimePaleoPhy method using "randObs" 
> # ancestors excluded
> # FAD.only = FALSE
> # FADs and LADs can differ
> ttrees_5a<- cal3TimePaleoPhy(
+ cladogram, rangesCont,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="randObs",
+ FAD.only=FALSE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Warning: Following taxa dropped from tree:
 t2, t6, t24, t30, t31
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_5a)
[1] FALSE
> 
> # RESULT
> # Expectation met, there is variation as expected.
> 
> ###########################################################################
> # 5b) What is the difference in behavior with `randObs` between 
> # `FAD.only = FALSE` and `FAD.only = TRUE` ? 
> # With FADs and LADs *not* being identical.
> # - Expectation: When `FAD.only` is `TRUE`, tip agesfrom `randObs`
> # should be equal to the given first appearance times, with no variation
> # across runs, as with `firstLast`.
> 
> # #cal3TimePaleoPhy method using "randObs" 
> # ancestors excluded
> # FAD.only = TRUE
> # FADs and LADs can differ
> ttrees_5b<- cal3TimePaleoPhy(
+ cladogram, rangesCont,
+ brRate=divRate, extRate=divRate,
+ sampRate=sRate, dateTreatment="randObs",
+ FAD.only=TRUE, ntrees=2,
+ anc.wt=0,plot=FALSE)
Error in cal3TimePaleoPhy(cladogram, rangesCont, brRate = divRate, extRate = divRate,  : 
  FAD.only = TRUE is inapplicable with dateTreatment = 'minMax' or 'randObs'
'randObs' assumes dates of observation unknown point occurrences within a taxon's range
   And thus there is no sensible reason for restricting dating to FADs if you chose
> 
> # test if tip ages have changed
> test_identical_tip_ages(trees = ttrees_5b)
Error in test_identical_tip_ages(trees = ttrees_5b) : 
  object 'ttrees_5b' not found
> 
> # RESULT
> # Returns an error, as it finds randObs incongruent with FAD.only = TRUE.

The behavior for 3b is now as expected, and all other cases are unimpacted.

dwbapst commented 5 years ago

The above test is now implemented via testthat so it will be checked on an on-going basis:

https://github.com/dwbapst/paleotree/blob/developmentBranch/tests/testthat/test_cal3_dateTreatment_Armin.R