# \donttest{
test = FALSE
if(test){
data(InfTemp)
set.seed(9999) # for reproducibility
# Let's start by simulating a trait under a climatic OU
beta = 0.6 # relationship to the climate curve
sim_theta = 4 # value of the optimum if the relationship to the climate
# curve is 0 (this corresponds to an 'intercept' in the linear relationship used below)
sim_sigma2 = 0.025 # variance of the scatter = sigma^2
sim_alpha = 0.36 # alpha value = strength of the OU; quite high here...
delta = 0.001 # time step used for the forward simulations => here its 1000y steps
tree <- pbtree(n=200, d=0.3) # simulate a bd tree with some extinct lineages
root_age = 60 # height of the root (almost all the Cenozoic here)
tree$edge.length <- root_age*tree$edge.length/max(nodeHeights(tree))
# here - for this contrived example - I scale the tree so that the root is at 60 Ma
# define a model - here we replicate the default model used in fit_t_env_ou
my_model <- function(t, env, param, theta0) theta0 + param[1]*env(t)
# simulate the traits
trait <- sim_t_env_ou(tree, sigma=sqrt(sim_sigma2), alpha=sim_alpha,
theta0=sim_theta, param=beta, model=my_model,
env_data=InfTemp, step=0.01, scale=TRUE, plot=TRUE)
## Fit the Environmental model (default)
result_fit <- fit_t_env_ou(phylo = tree, data = trait, env_data =InfTemp,
method = "Nelder-Mead", df=50, scale=TRUE)
plot(result_fit)
# We can also use the results from the previous fit to simulate a new dataset
trait2 <- sim_t_env_ou(tree, param=result_fit, step=0.001, plot=TRUE)
result_fit2 <- fit_t_env_ou(phylo = tree, data = trait2, env_data =InfTemp,
method = "Nelder-Mead", df=50, scale=TRUE)
result_fit2
}
# }
Run the code above in your browser using DataLab