More on bounded Brownian motion & important papers you haven't read yet

Several days ago now, I wrote a post about the most important phylogenetic comparative methods article you’ve probably never read.

To recap in a couple of sentences, the post refers an article by Florian Boucher & Vincent Démery that was published in the journal Systematic Biology in 2016 describing a new continuous trait evolution model of Brownian motion with reflexive bounds (BBM), as well as (arguably, more importantly) a method to calculate the likelihood numerically, which involves discretizing the state space of the continuous trait. The article thus identifies an important and (in retrospect, obvious) valuable equivalency (in the limit as the number of discrete state categories goes to \(\infty\)) between continuous & discrete-state trait evolution models. This framework was subsequently taken a step further in a 2018 article by Boucher et al. describing a generalization they refer to as the Fokker-Planck-Kolmogorov or FPK model. Together, both of these articles have accumulated a meager 72 citations (in total) over the 8 & 6 years since their publication, respectively, which I think really understates the potential significance of this approach to modeling & understanding trait evolution for continuous characters. More on this in future posts!

The other day, however, a colleague who was interested in BBM asked me the simple question “what about single-sided bounds?” In other words, what about a trait whose evolution is (for example) bounded on the lower, but not the upper end?

To think about this, though, I’m going to digress backwards a tiny bit and add an observation to my previous post on this subject in which I note that the BBM model will always explain our data better (will have a better model fit: i.e., higher log-likelihood) than its unbounded counterpart.

To understand why this might be, we should just imagine a scenario in which our trait evolved under Brownian motion with no bounds. In this case, the probability of observing the data that we have (in fact) observed will always be higher if there exist upper & lower bounds on trait evolution lying just beyond our uppermost and lowermost observations. This is because they make having observed these values of the trait more probable (because they add probability in the tails) than if the bounds don’t exist – while the probability of our data in the center of the distribution remain largely unaffected.

This also helps us to understand why bounds, even though they aren’t optimized, add parameters to our fitted model. (Another way of understanding this is that they are not given a priori – they must be observed from our data!) On the other hand, it is a completely reasonable endeavor to ask if the bounds, which add model complexity, add significant explanatory power to our fitted model.

Let’s see this by first simulating unbounded trait evolution and then fitting two models to these data: one without bounds & the other with bounds equal to the observed range of the trait. (I’m going to do visualize my simulation using a technique illustrated here, but this is just for fun. We could’ve likewise simply run fastBM on the original tree.)

library(phytools)
packageVersion("phytools")

## [1] '2.1.13'

tree<-pbtree(n=100,scale=100)
h<-max(nodeHeights(tree))
mt<-make.era.map(tree,seq(0,h,length.out=100))
st<-map.to.singleton(mt)
y<-fastBM(st,sig2=0.8,a=5,internal=TRUE)
par(mar=c(5.1,4.1,1.1,1.1))
phenogram(st,y,ftype="off",mar=c(5.1,4.1,1.1,1.1),
  lwd=5,las=1,cex.axis=0.8)
grid()
phenogram(st,y,ftype="off",add=TRUE,lwd=3,
  color="white")
phenogram(st,y,ftype="off",add=TRUE,lwd=1,
  color="blue")

x<-y[tree$tip.label]

Now, I’m going to fit a bounded & unbounded trait evolution models, and show that the former has a higher log-likelihood – even though the latter was our generating scenario. (The new function option I’m using here, lik.func="eigen", does faster calculation of the likelihood via eigen-decomposition of the Q matrix. I’m still trying to figure out if this makes it more or less precise!)

unbounded_fit<-bounded_bm(tree,x,lik.func="eigen",parallel=TRUE,
  root="mle",levs=400)
unbounded_fit

## Object of class "bounded_bm" based on
##  	a discretization with k = 400 levels.
## 
## Unwrapped (i.e., bounded) model
## 
## Set or estimated bounds: [ -70.506162 , 81.686123 ]
## 
## Fitted model parameters:
## 	 sigsq: 0.880901 
##      x0: 4.638779 
## 
## Log-likelihood: -286.695767 
## 
## R thinks it has found the ML solution.

(Here, since we don’t specify bounds, the function puts upper- and lower-bounds such that they evenly span 3 \(\times\) the observed range of values in our data. This is arbitrary, but as an approximation this range should be well outside the observed range of our data.)

We can see that this approximates geiger::fitContinuous for model="BM" (and this approximation would be even more accurate for increased levs).

geiger::fitContinuous(tree,x,model="BM")

## GEIGER-fitted comparative model of continuous data
##  fitted 'BM' model parameters:
## 	sigsq = 0.886079
## 	z0 = 4.767101
## 
##  model summary:
## 	log-likelihood = -287.039597
## 	AIC = 578.079193
## 	AICc = 578.202905
## 	free parameters = 2
## 
## Convergence diagnostics:
## 	optimization iterations = 100
## 	failed iterations = 0
## 	number of iterations with same best fit = 100
## 	frequency of best fit = 1.000
## 
##  object summary:
## 	'lik' -- likelihood function
## 	'bnd' -- bounds for likelihood search
## 	'res' -- optimization iteration summary
## 	'opt' -- maximum likelihood parameter estimates

bounded_fit<-bounded_bm(tree,x,lik.func="eigen",parallel=TRUE,
  root="mle",levs=400,lims=range(x))
bounded_fit

## Object of class "bounded_bm" based on
##  	a discretization with k = 400 levels.
## 
## Unwrapped (i.e., bounded) model
## 
## Set or estimated bounds: [ -19.775401 , 30.955361 ]
## 
## Fitted model parameters:
## 	 sigsq: 0.90261 
##      x0: 4.765605 
## 
## Log-likelihood: -284.36066 
## 
## R thinks it has found the ML solution.

lmtest::lrtest(unbounded_fit,bounded_fit)

## Likelihood ratio test
## 
## Model 1: unbounded_fit
## Model 2: bounded_fit
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1   2 -286.70                       
## 2   4 -284.36  2 4.6702     0.0968 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This model comparison happens to be non-significant, but it would seem to be more frequently so than expected under \(\alpha = 0.05\), suggesting that it might be necessary to generate a null distribution via simulation (or consider alternative ways of assessing model support).

Now let’s try the same simulation – but with a trait lower boundary of 0 & no upper boundary at all. Evolution under this scenario might look as follows.

tree<-pbtree(n=100,scale=100)
h<-max(nodeHeights(tree))
mt<-make.era.map(tree,seq(0,h,length.out=100))
st<-map.to.singleton(mt)
y<-fastBM(st,sig2=0.8,a=5,internal=TRUE,bounds=c(0,Inf))
par(mar=c(5.1,4.1,1.1,1.1))
phenogram(st,y,ftype="off",fsize=0.4,
  mar=c(5.1,4.1,1.1,1.1),lwd=5,las=1,cex.axis=0.8,
  spread.cost=c(1,0),ylim=c(-5,max(y)))
grid()
phenogram(st,y,ftype="off",add=TRUE,lwd=3,
  color="white")
phenogram(st,y,ftype="off",add=TRUE,lwd=1,
  color="blue")
segments(x=0,y=0,x1=h,y0=0,
  lwd=6,col=palette()[2])
segments(x=0,y=0,x1=h,y0=0,
  lwd=2,lty="dashed")

x<-y[tree$tip.label]

Here, I’ll first fit a model with just a lower bound, then one with both a lower & upper bound, and finally one that is unbounded. I’m going to still use lik.func="eigen", but I’ll revert to root="nuisance" and levs=200, the defaults.

lower_bounded<-bounded_bm(tree,x,lims=c(0,Inf),lik.func="eigen",
  parallel=TRUE)
lower_bounded

## Object of class "bounded_bm" based on
##  	a discretization with k = 200 levels.
## 
## Unwrapped (i.e., bounded) model
## 
## Set or estimated bounds: [ 0 , 47.531448 ]
## 
## Fitted model parameters:
## 	 sigsq: 0.802814 
##      x0: 4.301189 
## 
## Log-likelihood: -253.098453 
## 
## R thinks it has found the ML solution.

both_bounded<-bounded_bm(tree,x,lims=c(0,max(x)),lik.func="eigen",
  parallel=TRUE)
both_bounded

## Object of class "bounded_bm" based on
##  	a discretization with k = 200 levels.
## 
## Unwrapped (i.e., bounded) model
## 
## Set or estimated bounds: [ 0 , 23.815825 ]
## 
## Fitted model parameters:
## 	 sigsq: 0.803155 
##      x0: 4.351448 
## 
## Log-likelihood: -251.514354 
## 
## R thinks it has found the ML solution.

unbounded<-bounded_bm(tree,x,lik.func="eigen",parallel=TRUE)
unbounded

## Object of class "bounded_bm" based on
##  	a discretization with k = 200 levels.
## 
## Unwrapped (i.e., bounded) model
## 
## Set or estimated bounds: [ -23.615421 , 47.531448 ]
## 
## Fitted model parameters:
## 	 sigsq: 0.624178 
##      x0: 7.692622 
## 
## Log-likelihood: -266.082377 
## 
## R thinks it has found the ML solution.

AIC(unbounded,lower_bounded,both_bounded)

##               df      AIC
## unbounded      2 536.1648
## lower_bounded  3 512.1969
## both_bounded   4 511.0287

This shows us how difficult it might be to distinguish between trait evolution with just a bound on one side of the observed distribution of the trait, compared to evolution with both upper & lower bounds. On the other hand, we can easily reject our null hypothesis of unbounded Brownian motion evolution for this trait!

More on this soon.