trait-simulator: Simulates the Evolution of a Quantitative Trait.
In MPSEM: Modeling Phylogenetic Signals using Eigenvector Maps

Description Usage Arguments Details Value Functions Author(s) References Examples

Functions to simulate the evolution of a quantitative trait along a phylogenetic tree inputted as an object of class ‘phylo’ (package ape) or graph-class object.

EvolveOptimMarkovTree(tp, tw, anc, p = 1, root = tp$edge[1, 1])

TraitOUsimTree(tp, a, sigma, opt, p = 1, root = tp$edge[1, 1])

OUvar(d, a = 0, theta = 1, sigma = 1)

PEMvar(d, a = 0, psi = 1)

TraitVarGraphSim(x, variance, distance = "distance", p = 1, ...)

`tp`	A rooted phylogenetic tree of class ‘phylo’ (see package ape).
`tw`	Transition matrix giving the probability that the optimum trait value changes from one state to another at vertices. All rows must sum to 1.
`anc`	Ancestral state of a trait (at the root).
`p`	Number of variates to generate.
`root`	Root node of the tree.
`a`	Selection rate (`OUvar`) or steepness (`PEMvar`).
`sigma`	Neutral evolution rate, i.e. mean trait shift by drift.
`opt`	An index vector of optima at the nodes.
`d`	Phylogenetic distances (edge lengths).
`theta`	Adaptive evolution rate, i.e. mean trait shift by natural selection.
`psi`	Mean evolution rate.
`x`	A `graph-class` object.
`variance`	Variance function (`OUvar`, `PEMvar`, or any suitable user-defined function).
`distance`	The name of the member of ‘x$edge’ where edge lengths can be found.
`...`	Additional parameters for the specified variance function.

Function EvolveOptimMarkovTree allows one to simulate the changes of optimum trait values as a Markov process. The index whereby the process starts, at the tree root, is set by parameter anc; this is the ancestral character state. From the root onwards to the tips, the optimum is given the opportunity to change following a multinomial random draw with transition probabilities given by the rows of matrix tw. The integers thus obtained can be used as indices of a vector featuring the actual optimum trait values corresponding to the simulated selection regimes. The resulting optimum trait values at the nodes are used by TraitOUsimTree as its parameters opt to simulate trait values at nodes and tips. Function TraitVarGraphSim uses a graph variance function (either OUvar or PEMvar) to reconstruct a covariance matrix that is used to generate covariates drawn from a multi-normal distribution.

Functions EvolveOptimMarkovTree and TraitOUsimTree return a matrix whose rows represent the vertices (nodes and tips) of the phylogenetic tree and whose columns stand for the n different trials the function was asked to perform. For EvolveQTraitTree, the elements of the matrix are integers, representing the selection regimes prevailing at the nodes and tips, whereas for TraitOUsimTree, the elements are simulated quantitative trait values at the nodes and tips. These functions are implemented in C language and therefore run swiftly even for large (10000+ species) trees.

Function TraitVarGraphSim returns p phylogenetic signals and is implemented using a rotation of a matrix of standard normal random (mean=0, variance=1) deviates. The rotation matrix is itself obtained by Choleski factorization of the trait covariance matrix expected for a given set of trees, variance function, and variance function parameters.

EvolveOptimMarkovTree: Simulates the evolution of trait optima along a phylogeny.
TraitOUsimTree: Simulates the evolution of trait values along a phylogeny.
OUvar: Describe here...
PEMvar: Describe here...
TraitVarGraphSim: Describe here...

Guillaume Guenard, with contribution from Pierre Legendre Maintainer: Guillaume Guenard <guillaume.guenard@gmail.com>

Butler, M. A. & King, A. A. 2004. Phylogenetic comparative analysis: a modeling approach for adaptive evolution. Am. Nat. 164: 683-695.

Guénard, G., Legendre, P., and Peres-Neto, P. 2013. Phylogenetic eigenvector maps (PEM): a framework to model and predict species traits. Meth. Ecol. Evol. 4: 1120–1131

opt <- c(-2,0,2) # Three trait optima: -2, 0, and 2
## Transition probabilities:
transit <- matrix(c(0.7,0.2,0.2,0.2,0.7,0.1,0.1,0.1,0.7),
                  length(opt),length(opt),dimnames=list(from=opt,to=opt))

## In this example, the trait has a probability of 0.7 to stay at a given
## optimum, a probability of 0.2 for the optimum to change from -2 to 0,
## from 0 to -2, and from 2 to -2, and a probability of 0.1 for the
## optimum to change from -2 to 2, from 0 to 2, and from 2 to 0.
nsp <- 25  # A random tree for 25 species.
tree2 <- rtree(nsp,tip.label=paste("Species",1:nsp,sep=""))
tree2$node.label=paste("N",1:tree2$Nnode,sep="")  # Node labels.

## Simulate 10 trials of optimum change.
reg <- EvolveOptimMarkovTree(tp=tree2,tw=transit,p=10,anc=2)
y1 <- TraitOUsimTree(tp=tree2,a=0,sigma=1,
                     opt=opt[reg[,1]],p=10)    ## Neutral
y2 <- TraitOUsimTree(tp=tree2,a=1,sigma=1,
                     opt=opt[reg[,1]],p=10)    ## Few selection.
y3 <- TraitOUsimTree(tp=tree2,a=10,sigma=1,
                     opt=opt[reg[,1]],p=10)    ## Strong selection.

## Display optimum change with colours.
displayOUprocess <- function(tp,trait,regime,mvalue) {
  layout(matrix(1:2,1,2))
  n <- length(tp$tip.label)
  ape::plot.phylo(tp,show.tip.label=TRUE,show.node.label=TRUE,root.edge=FALSE,
                  direction="rightwards",adj=0,
                  edge.color=rainbow(length(trait))[regime[tp$edge[,2]]])
  plot(y=1:n,x=mvalue[1:n],type="b",xlim=c(-5,5),ylab="",xlab="Trait value",yaxt="n",
       bg=rainbow(length(trait))[regime[1:n]],pch=21) 
  text(trait[regime[1:n]],y=1:n,x=5,col=rainbow(length(trait))[regime[1:n]])
  abline(v=0)
}

displayOUprocess(tree2,opt,reg[,1],y1[,1])  # Trait evolve neutrally,
displayOUprocess(tree2,opt,reg[,1],y2[,1])  # under weak selection,
displayOUprocess(tree2,opt,reg[,1],y3[,1])  # under strong selection.

x <- Phylo2DirectedGraph(tree2)
y4 <- TraitVarGraphSim(x, variance = OUvar, p=10, a=5)

DisplayTreeEvol <- function(tp,mvalue) {
  layout(matrix(1:2,1,2))
  n <- length(tp$tip.label)
  ape::plot.phylo(tp,show.tip.label = TRUE, show.node.label = TRUE,
                  root.edge = FALSE, direction = "rightwards", adj = 0)
  plot(y=1:n, x=mvalue[1:n], type="b", xlim=c(-5,5), ylab="",
       xlab="Trait value", yaxt="n", pch=21)
  abline(v=0)
}

## Recursively displays the simulated traits.
for(i in 1:10) {
  DisplayTreeEvol(tree2,y4[i,])
  if(is.null(locator(1)))
    break                  ## Stops recursive display on a mouse right-click.
}