BrownianMotionModel {mvSLOUCH}

Estimate parameters under a Brownian motion model of evolution

Description

The BrownianMotionModel function uses maximum likelihood to fit parameters of a Brownian motion model evolving on the phylogeny. The user is recommended to install the suggested package PCMBaseCpp which significantly speeds up the calculations (see Details).

Usage

BrownianMotionModel(phyltree, mData, predictors = NULL, M.error = NULL, 
min_bl = 0.0003)

Arguments

phyltree

The phylogeny in phylo format. The tree can be obtained from e.g. a nexus file by the read.nexus() function from the ape package. The "standard" ape node indexing is assumed: for a tree with n tips, the tips should have indices 1:n and the root index n+1. The root.edge field is ignored.

mData

A matrix with the rows corresponding to the tip species while the columns correspond to the traits. The rows should be named by species
(field phyltree$tip.label), if not, then a warning is thrown and the order of the species is assumed to be the same as the order in which the species are in the phylogeny (i.e. correspond to the node indices 1:n, where n is the number of tips). The columns should be named by traits, otherwise a warning is thrown and generic names are generated.

predictors

A vector giving the numbers of the columns from data which are to be considered predictor ones, i.e. conditioned on in the program output. If not provided, then none will be considered to be predictors.

M.error

An optional measurement error covariance structure. The measurement errors between species are assumed independent. The program tries to recognize the structure of the passed matrix and accepts the following possibilities :

• a single number that is a common measurement error for all tips and species,

• a m element vector with each value corresponding to a variable, measurement errors are independent between variables and each species is assumed to have the same measurement errors,

• a m x m ((number of variables) x (number of variables)) matrix, all species will have the same measurement error,

• a list of length n (number of species), each list element is the covariance structure for the appropriate (numbering according to tree) species, either a single number (each variable has same variance), vector (of length m for each variable), or m x m matrix, the order of the list has to correspond to the order of the nodes in the phyltree object,

• NULL no measurement error.

From version 2.0.0 of mvSLOUCH it is impossible to pass a single joint measurement error matrix for all the species and traits.

min_bl

Value to which PCMBase's PCMBase.Threshold.Skip.Singular should be set. It indicates that branches of length shorter than min_bl should be skipped in likelihood calculations. Short branches can result in singular covariance matrices for the transition density along a branch. The user should adjust this value if a lot of warnings are raised by PCMBase about singularities during the likelihood calculations. However, this does not concern tip branches-these cannot be skipped and hence should be long enough so that numerical issues are not raised.

Details

The likelihood calculations are done by the PCMBase package. However, there is a C++ backend, PCMBaseCpp. If it is not available, then the likelihood is calculated slower using pure R. However, with the calculations in C++ up to a 100-fold increase in speed is possible (more realistically 10-20 times). The PCMBaseCpp package is available from https://github.com/venelin/PCMBaseCpp.

This function estimates the parameters of a multivariate Brownian motion model defined by the SDE,

dX(t) = \Sigma dB(t), X(0)=X_{0}

evolving on a phylogenetic tree.

Without measurement error the parameters are obtained analytically via a GLS procedure. If measurement error is present, then the parameters are optimized over using optim(). The initial conditions for the optimization are motivated by Bartoszek \& Sagitov (2015)'s univariate results.

From version 2.0.0 of mvSLOUCH the data has to be passed as a matrix. To underline this the data parameter's name has been changed to mData.

The phyltree_paths() function enhances the tree for usage by mvSLOUCH. Hence, to save time, it is advisable to first do phyltree<-mvSLOUCH::phyltree_paths(phyltree) and only then use it with BrownianMotionModel().

From version 2.0.0 of mvSLOUCH the parameter calcCI has been removed. The package now offers the possibility of bootstrap confidence intervals, see function parametric.bootstrap.

Value

ParamsInModel

A list with estimated model parameters. The elements are vX0 : the ancestral trait, and Sxx where t\Sigma_{xx}\Sigma_{xx}^{T} is the Brownian motion's covariance matrix at time t.

ParamSummary

A list with summary statistics with elements, StS the infinitesimal covariance matrix \Sigma_{xx}\Sigma_{xx}^{T}, LogLik the log–likelihood, dof the degrees of freedom, m2loglik is -2log–likelihood, aic is the Akaike information criterion, aic.c is the Akaike information criterion corrected for small sample size, sic is the Schwarz information criterion, bic is the Bayesian information criterion (which is the same as the Schwarz information criterion) and RSS is the residual sum of squares.

Author(s)

Krzysztof Bartoszek

References

Bartoszek, K. and Pienaar, J. and Mostad. P. and Andersson, S. and Hansen, T. F. (2012) A phylogenetic comparative method for studying multivariate adaptation. Journal of Theoretical Biology 314:204-215.

Bartoszek, K. and Sagitov S. (2015) A consistent estimator of the evolutionary rate. Journal of Theoretical Biology 371:69-78.

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

Felsenstein, J. (1985) Phylogenies and the comparative method. American Naturalist 125:1-15.

Hansen, T.F. and Bartoszek, K. (2012) Interpreting the evolutionary regression: the interplay between observational and biological errors in phylogenetic comparative studies. Systematic Biology 61(3):413-425.

Mitov, V. and Bartoszek, K. and Asimomitis, G. and Stadler, T. (2020) Fast likelihood calculation for multivariate Gaussian phylogenetic models with shifts Theoretical Population Biology 131:66-78.

Pienaar et al (in prep) An overview of comparative methods for testing adaptation to external environments.

See Also

brown,mvBM, PCMLik, SummarizeBM, simulBMProcPhylTree, parametric.bootstrap

Examples

RNGversion(min(as.character(getRversion()),"3.6.1"))
set.seed(12345, kind = "Mersenne-Twister", normal.kind = "Inversion")
### We will first simulate a small phylogenetic tree  using functions from ape. 
### For simulating the tree one could also use alternative functions, e.g. sim.bd.taxa 
### from the TreeSim package
phyltree<-ape::rtree(5)

## The line below is not necessary but advisable for speed
phyltree<-phyltree_paths(phyltree)

### Define Brownian motion parameters to be able to simulate data under 
### the Brownian motion model.
BMparameters<-list(vX0=matrix(0,nrow=3,ncol=1),
Sxx=rbind(c(1,0,0),c(0.2,1,0),c(0.3,0.25,1)))

### Now simulate the data.
BMdata<-simulBMProcPhylTree(phyltree,X0=BMparameters$vX0,Sigma=BMparameters$Sxx)
BMdata<-BMdata[phyltree$tip.label,,drop=FALSE]

### Recover the parameters of the Brownian motion.
BMestim<-BrownianMotionModel(phyltree,BMdata)

## Not run: 
### And finally obtain bootstrap confidence intervals for some parameters
BMbootstrap<-parametric.bootstrap(estimated.model=BMestim,phyltree=phyltree,
values.to.bootstrap=c("vX0","StS"),M.error=NULL,numboot=2)

## End(Not run)
RNGversion(as.character(getRversion()))