View source: R/OUphylregression.R

OU_RSS | R Documentation |

The `OU_RSS`

function calculates the residual sum of squares (RSS)
for given data under a multivariate OU model evolving on the phylogeny.
The user is recommended to install the suggested package
PCMBaseCpp which significantly speeds up the calculations (see Details).

```
OU_RSS(mY, phyltree, modelParams, M.error = NULL, do_centre = NA,
regimes = NULL, regimes.times = NULL, root.regime = NULL)
```

`mY` |
A matrix with the rows corresponding to the tip species while the columns correspond to the traits.
The rows should be named by species |

`phyltree` |
The phylogeny in |

`modelParams` |
List of model parameters of the BM/OUOU/OUBM model as |

`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 |

`do_centre` |
Should the data, |

`regimes` |
A vector or list of regimes. If vector then each entry corresponds to each of the branches of |

`regimes.times` |
A list of vectors for each tree node, it starts with |

`root.regime` |
The regime at the root of the tree. If not given, then it is taken as the regime that is present
on the daughter lineages stemming from the root and is the most frequent one in the |

The matrix algebra calculations are done using the likelihood function offerred 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.

For a given input data matrix, `mY`

, the function considers the stacking of it by rows
(i.e. stacking species by species). Let Y = vec(`mY`

), i.e. `Y<-c(t(mY))`

,
V be the between-species-between-traits variance-covariance matrix (under the parameters passed
in `modelParams`

) and `v`

a centring vector (if `do_centre`

is `NA`

, then
`v=0`

). The function calculates the value of the quadratic form

`(Y-v)^{T}V^{-1}(Y-v).`

A special centring is when `do_centre`

equals `"phylaverage"`

. In this situation
the centring vector is a phylogenetically weighted average, i.e.

`v=(D^{T}V^{-1}D)^{-1}D^{T}V^{-1}Y ,`

where denoting `1_{n}`

as a column vector of n ones and `Id_{k}`

as
the identity matrix with rows and columns equalling the number or columns of `mY`

,

`D=1_{n}\otimes Id_{k}.`

The value of the residual sum of squares, quadratic form with respect to the between-species-between-traits precision matrix. Also the used phylogeny is returned.

Krzysztof Bartoszek

```
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 a vector of regimes.
regimes<-c("small","small","large","small","small","large","large","large")
### Define SDE parameters to be able to simulate data under the OUOU model.
## 3D model
## OUOUparameters<-list(vY0=matrix(c(1,-1,0.5),nrow=3,ncol=1),
## A=rbind(c(9,0,0),c(0,5,0),c(0,0,1)),mPsi=cbind("small"=c(1,-1,0.5),"large"=c(-1,1,0.5)),
## Syy=rbind(c(1,0.25,0.3),c(0,1,0.2),c(0,0,1)))
## 2D model used to reduce running time on CRAN
OUOUparameters<-list(vY0=matrix(c(1,-1),nrow=2,ncol=1),
A=rbind(c(9,0),c(0,5)),mPsi=cbind("small"=c(1,-1),"large"=c(-1,1)),
Syy=rbind(c(1,0.25),c(0,1)))
### Now simulate the data.
OUOUdata<-simulOUCHProcPhylTree(phyltree,OUOUparameters,regimes,NULL)
OUOUdata<-OUOUdata[phyltree$tip.label,,drop=FALSE]
## below will return the RSS under the assumed OUOU model of evolution
OU_RSS(OUOUdata, phyltree, OUOUparameters, M.error=NULL,
do_centre="evolutionary_model", regimes = regimes)
```

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.