trait-simulator: Simulate the Evolution of a Quantitative Trait

In MPSEM: Modelling Phylogenetic Signals using Eigenvector Maps

Simulate the Evolution of a Quantitative Trait

Description

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

Usage

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, ...)

Arguments

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 (row) to another (column) 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 in function (OUvar) or steepness in (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 other 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.

Details

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 argument 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, used to generate covariates drawn from a multi-normal distribution.

Value

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. It 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.

Functions

• EvolveOptimMarkovTree(): Trait Optima Simulator

  Simulates the evolution of trait optima along a phylogeny as a Markov process.

• TraitOUsimTree(): Trait Value Simulator

  Simulates the evolution of trait values along a phylogeny as a Ornstein–Uhlenbeck process.

• OUvar(): Ornstein–Uhlenbeck Variance Calculator

  Calculates the expected covariance matrix for a trait evolving following an Ornstein–Uhlenbeck process. This function is meant to be used with function TraitVarGraphSim.

• PEMvar(): Phylogenetic Eigenvector Maps Variance Calculator

  Calculates the covariance on the basis of the covariance model (power function) associated used in calculating Phylogenetic Eigenvector Maps. This function is meant to be used with function TraitVarGraphSim.

• TraitVarGraphSim(): Covariance-based Trait Evolution Simulator.

  Simulates trait evolution as covariates drawn from a multi-normal distribution whose covariance is estimated using an external function (functions OUvar, PEMvar provided with the package or any user-provided function).

Author(s)

Guillaume Guénard [aut, cre] (ORCID: <https://orcid.org/0000-0003-0761-3072>), Pierre Legendre [ctb] (ORCID: <https://orcid.org/0000-0002-3838-3305>) Maintainer: Guillaume Guénard <[e-mail]>

References

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

Guénard, G., Legendre, P., and Peres-Neto, P. 2013. Phylogenetic eigenvector maps: a framework to model and predict species traits. Methods in Ecology and Evolution 4: 1120-1131

Examples

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)
}

## Iteratively displays the simulated traits.
## Left-click on the display area to go to the next plot.
## To terminate: right-click (WIndows, X11), esc key (Mac), or hit the
## "finish" button (RStudio).

for(i in 1:10) {
  DisplayTreeEvol(tree2,y4[i,])
  if(is.null(locator(1)))
    break  ## Terminate: 
}

MPSEM documentation built on June 8, 2025, 1:56 p.m.