hmorlon
commented
6 years ago Hello Miao,
Because of the way the optimization is performed, and because speciation and extinction rates are always positive, we use absolute values of f.lamb and f.mu in the likelihood computation. So, the outputs of f.lamb and f.mu should also be interpreted as absolute values. This means that the estimated extinction rate through time is abs(f.mu). Here, for example, it is abs(-1.342-0.000335t)=1.342+0.000335t. So the estimated extinction rate at present is 1.342, and it increases from present to past (i.e. it decreases from past to present). Note that in other situations where mu_par[1]+mu_par[2]t is positive for some t values and negative for others, it should always be interpreted as abs(mu_par[1]+mu_par[2]t), which might not be a linear function (in can be linear increasing on one part of the time interval and linear decreasing in another part of the time interval for example). If this is confusing, maybe the easiest is to use the plot_fit_bd function to plot how the rates vary through time. If you do this, please use the github version of the codes >library(devtools)
install_github("hmorlon/PANDA") as there have been some modifications of the codes (in particular accounting for this absolute values issue) that have not yet been commited to CRAN. If you want to force a monotonous function for a linear dependency (i.e. either increasing or decreasing function of time), then instead of mu.l <- function(t,y){y[1] + y[2]*t} you should use a function that returns the max of this and 0. Also, if you want to avoid these issues, you might want to consider using exponential rather than linear dependencies, as exponential functions remain always >0. I hope this helps. Best. Hélène
Cactusolo
Hi Helene,
I run a few models with RPANDA, and I got a bunch negative values for mu (extinction rate). would you please shed some lights on how should I interpret negative mu values? is negative value correct?
here are the basic code I used:
lambda.cst <- function(t,y){y[1]} #t=time, y=lambda (λ) mu.l <- function(t,y){y[1] + y[2]*t} par <- c(0.05, 0.005, 0.001) bcst.dvar.l <- fit_bd(tree, tot_time, f.lamb=lambda.cst, f.mu=mu.l, lamb_par=par[1], mu_par=par[c(2,3)], cst.lamb=TRUE, cond="crown", f=fraction, dt=1e-3)and got results like this:$bcst.dvar.l` model : [1] "birth death"LH : [1] -58309.65
aicc : [1] 116625.3
lamb_par : [1] 1.406383
mu_par : [1] -1.3424997816 -0.0003351743 `
I have read your RPANDA tutorial on Morlon et al. (2014) in MEE journal, but the mu_par[1] and mu[2] are not explained, except res$lamb_par[1] and res$lamb_par[2]. So would you please explain what does it mean if I got a negative value for mu_par[1]? for mu_par[2], does it mean increased the speciation rate?
Thanks!
Miao