simulate.white.noise: simulate.white.noise

In sebastianduchene/NELSI: NELSI: Nucleotide EvoLution Simulator

simulate.white.noise

Description

simulate.white.noise simulates the rate along a phylogeneti tree according to a white noise process. In particular, it samples from a lognormal distribution with a given mean and standard deviation. For each branch a rate is drawn from a distribution with the given mean. The standard deviation is divided by the branch length divided by the mean of the branch lengths. As a result, the rate of a long branch will be drawn from a distribution with a lower standard deviation than that for a short branch (i.e. the rate is more uncertain for shorter branches).

Usage

simulate.white.noise(tree, params = NULL)

Arguments

tree	A phylogenetic tree of class 'phylo'. The branch lengths should be in units of time (a chronogram).
params	A list with the parameters for the simulation function, corresponding to the following two items:
mean.log	The mean of the lognormal distribution. It should be in log scale
sd.log	The standard deviation of the lognormal distribution. The actual value used for each branch will vary inversely with the branch length. See the description of the function

Details

None

Value

An object of class 'ratesim'. This is a list with two items:

phylogram	The phylogenetic tree with branch lengths in units of substitutions (phylogram)
tree.data.matrix	This is a matrix with the number of substitutions, rates, and times along every branch in the tree. See get.tree.data.matrix for more details

Note

none

Author(s)

Sebastian Duchene and David Duchene

References

A similar model is described in: Lepage, Thomas, et al. "A general comparison of relaxed molecular clock models." Molecular biology and evolution 24.12 (2007): 2669-2680.

See Also

simulate.uncor.lnorm

Examples

set.seed(1234525)

myTree <- rcoal(50)

rateTree <- simulate.white.noise(tree = myTree, params = NULL)
plot(rateTree, col.lineages = rainbow(50))

#See the histogram of the branch-wise rates
hist(rateTree$tree.data.matrix[, 5])

## The function is currently defined as
function (tree, params = NULL) 
{
    data.matrix <- get.tree.data.matrix(tree)
    means.rates <- data.matrix[, 7] / sum(data.matrix[, 7])
    branch.rates <- sapply(means.rates, function(x) rlnorm(1, log(x), x))
    data.matrix[, 5] <- branch.rates
    data.matrix[, 6] <- data.matrix[, 5] * data.matrix[, 7]
    tree$edge.length <- data.matrix[, 6]
    res <- list(tree, tree.data.matrix = data.matrix)
    class(res) <- "ratesim"
    return(res)
  }

sebastianduchene/NELSI documentation built on Sept. 5, 2024, 3:08 a.m.