evolvability (version 2.0.1)

simulate_rate: Simulating evolutionary rate model

Description

simulate_rate Simulates three different evolutionary rates models. In the two first, 'predictor_BM' and 'predictor_gBM', the evolution of y follows a Brownian motion with variance linear in x, while the evolution of x either follows a Brownian motion or a geometric Brownian motion, respectively. In the third model, 'recent_evol', the residuals of the macroevolutionary predictions of y have variance linear in x.

Usage

simulate_rate(
  tree,
  startv_x = NULL,
  sigma_x = NULL,
  a,
  b,
  sigma_y = NULL,
  x = NULL,
  model = "predictor_BM"
)

Value

An object of class

'simulate_rate', which is a list with the following components:

tips				A data frame of x and y values for the tips.	percent_negative_roots
			The percent of iterations with negative roots in the rates of y (not given for model = 'recent_evol').	compl_dynamics

Arguments

tree

A phylo object. Must be ultrametric and scaled to unit length.

startv_x

The starting value for x (usually 0 for 'predictor_BM' and 'recent_evol', and 1 for 'predictor_gBM').

sigma_x

The evolutionary rate of x.

a

A parameter of the evolutionary rate of y.

b

A parameter of the evolutionary rate of y.

sigma_y

The evolutionary rate for the macroevolution of y (Brownian motion process) used in the 'recent_evolution' model.

x

Optional fixed values of x (only for the 'recent_evol' model), must equal number of tips in the phylogeny, must correspond to the order of the tip labels.

model

Either a Brownian motion model of x 'predictor_BM', geometric Brownian motion model of x 'predictor_gBM', or 'recent_evol'.

Author

Geir H. Bolstad

Details

See the vignette 'Analyzing rates of evolution' for an explanation of the evolutionary models. The data of the tips can be analyzed with the function rate_gls. Note that a large part of parameter space will cause negative roots in the rates of y (i.e. negative a+bx). In these cases the rates are set to 0. A warning message is given if the number of such instances is larger than 0.1%. For model 1 and 2, it is possible to set `a='scaleme'`, if this chosen then `a` will be given the lowest possible value constrained by a+bx>0.

Examples

# Also see the vignette 'Analyzing rates of evolution'.
if (FALSE) {
# Generating a tree with 50 species
set.seed(102)
tree <- ape::rtree(n = 50)
tree <- ape::chronopl(tree, lambda = 1, age.min = 1)

### model = 'predictor_BM' ###
sim_data <- simulate_rate(tree, startv_x = 0, sigma_x = 0.25, a = 1, b = 1, 
model = "predictor_BM")
head(sim_data$tips)
par(mfrow = c(1, 3))
plot(sim_data)
plot(sim_data, response = "y")
plot(sim_data, response = "x")
par(mfrow = c(1, 1))

### model = 'predictor_gBM' ###
sim_data <- simulate_rate(tree, startv_x = 1, sigma_x = 1, a = 1, b = 0.1, 
model = "predictor_gBM")
head(sim_data$tips)
par(mfrow = c(1, 3))
plot(sim_data)
plot(sim_data, response = "y")
plot(sim_data, response = "x")
par(mfrow = c(1, 1))

### model = 'recent_evol' ###
sim_data <- simulate_rate(tree,
  startv_x = 0, sigma_x = 1, a = 1, b = 1, sigma_y = 1,
  model = "recent_evol"
)
head(sim_data$tips)
}