Nothing
################################################################################
## ##
## RPANDA : Utils ##
## ##
## Julien Clavel - 01-02-2018 ##
## S3 methods, simulations, miscellaneous ##
## ##
################################################################################
# S3 generic method "ancestral" for reconstructing or retrieving ancestral states (see phyl.pca_pl.R) # should I use "predict" instead?
#ancestral <- function(object) UseMethod("ancestral")
# S3 for the fit_t_env class
ancestral.fit_t.env <- function(object, ...){
# extract objects
if(!inherits(object,"fit_t.env")) stop("only works with \"fit_t.env\" class objects. See ?fit_t_env")
# Ancestral state at the root
a <- object$root
names(a) = "root"
res <- a
return(res)
warning("only the root state is currently estimated for models of the class \"fit_t.env\"") # To remove later
}
# S3 for the fit_t_comp class?
# TODO
ancestral.fit_t.comp <- function(object, ...){
# extract objects
if(!inherits(object,"fit_t.comp")) stop("only works with \"fit_t.comp\" class objects. See ?fit_t_comp")
anc <- object$z0
names(anc) ="root"
return(anc)
warning("only the root state is currently estimated for models of the class \"fit_t.comp\"")
}
# Build a matrix with tip and internal covariances
.vcvPhyloInternal <- function(tree){
nbtip <- Ntip(tree)
dis <- dist.nodes(tree)
MRCA <- mrca(tree, full = TRUE)
M <- dis[as.character(nbtip + 1), MRCA]
dim(M) <- rep(sqrt(length(M)), 2)
return(M)
}
# Build the matrix square root and inverse matrix square root of the phylogenetic covariance matrix
.transformsqrt <- function(tree){
vcv_tr <- vcv.phylo(tree)
eig <- eigen(vcv_tr)
sqrtmValues <- sqrt(eig$values)
sqrtM1 <- tcrossprod(eig$vectors%*%diag(1/sqrtmValues), eig$vectors)
sqrtM <- tcrossprod(eig$vectors%*%diag(sqrtmValues), eig$vectors)
matrices <- list(sqrtM1=sqrtM1, sqrtM=sqrtM)
return(matrices)
}
# ---- Function to simulate random covariance matrices with a specified eigenstructure
# From Uyeda et al. 2015 - Systematic Biology 64(4):677-689.
Posdef <- function (p, ev = rexp(p, 1/100)) {
Z <- matrix(ncol=p, rnorm(p^2))
decomp <- qr(Z)
Q <- qr.Q(decomp)
R <- qr.R(decomp)
d <- diag(R)
ph <- d / abs(d)
O <- Q %*% diag(ph)
Z <- t(O) %*% diag(ev) %*% O
return(Z)
}
# --- Function to simulate multivariate normal distribution
# From the mvtnorm package
rmvnorm_util<-function (n, mean = rep(0, nrow(sigma)), sigma = diag(length(mean)),
method = c("eigen", "svd", "chol"), pre0.9_9994 = FALSE,
checkSymmetry = TRUE)
{
if (checkSymmetry && !isSymmetric(sigma, tol = sqrt(.Machine$double.eps),
check.attributes = FALSE)) {
stop("sigma must be a symmetric matrix")
}
if (length(mean) != nrow(sigma))
stop("mean and sigma have non-conforming size")
method <- match.arg(method)
R <- if (method == "eigen") {
ev <- eigen(sigma, symmetric = TRUE)
if (!all(ev$values >= -sqrt(.Machine$double.eps) * abs(ev$values[1]))) {
warning("sigma is numerically not positive semidefinite")
}
t(ev$vectors %*% (t(ev$vectors) * sqrt(pmax(ev$values,
0))))
}
else if (method == "svd") {
s. <- svd(sigma)
if (!all(s.$d >= -sqrt(.Machine$double.eps) * abs(s.$d[1]))) {
warning("sigma is numerically not positive semidefinite")
}
t(s.$v %*% (t(s.$u) * sqrt(pmax(s.$d, 0))))
}
else if (method == "chol") {
R <- chol(sigma, pivot = TRUE)
R[, order(attr(R, "pivot"))]
}
retval <- matrix(rnorm(n * ncol(sigma)), nrow = n, byrow = !pre0.9_9994) %*%
R
retval <- sweep(retval, 2, mean, "+")
colnames(retval) <- names(mean)
retval
}
## Wrapper to compute the inverse fourier transform as in pracma using fft from "stats"
## J. Clavel
ifft_panda <- function(x) return(fft(x, inverse=TRUE)*(1/length(x)))
## Compute the pseudoinverse from svd
## J. Clavel modified from corpcor
## most of the pseudoinverse are applied to square matrices in RPANDA. No need for fat/thin matrices implementations
pseudoinverse = function (m, tol)
{
# compute the svd
s = svd(m)
if( missing(tol) )
tol = max(dim(m))*max(s$d)*.Machine$double.eps
Positive = s$d > tol
if (length(s$d) == 0)
{
return(
array(0, dim(m)[2:1])
)
}
else
{
return(
s$v[, Positive, drop=FALSE] %*% (1/s$d[Positive] * t(s$u[, Positive, drop=FALSE]))
)
}
}
## detect extinct and remove them - from geiger 2.0.7
## J. Clavel
drop.extinct <- function (phy, tol=NULL) {
if (!"phylo" %in% class(phy)) {
stop("\"phy\" is not of class \"phylo\".");
}
if (is.null(phy$edge.length)) {
stop("\"phy\" does not have branch lengths.");
}
if (is.null(tol)) {
tol <- min(phy$edge.length)/100;
}
aa <- is.extinct(phy=phy, tol=tol);
if (length(aa) > 0) {
phy <- drop.tip(phy, aa); # use drop.tip from "ape" => Imports
}
return(phy);
}
# return tip.labels, so that tree ordering is not an issue
is.extinct <- function (phy, tol=NULL) {
if (!"phylo" %in% class(phy)) {
stop("\"phy\" is not of class \"phylo\".");
}
if (is.null(phy$edge.length)) {
stop("\"phy\" does not have branch lengths.");
}
if (is.null(tol)) {
tol <- min(phy$edge.length)/100;
}
phy <- reorder(phy);
xx <- numeric(Ntip(phy) + phy$Nnode);
for (i in 1:length(phy$edge[,1])) {
xx[phy$edge[i,2]] <- xx[phy$edge[i,1]] + phy$edge.length[i];
}
aa <- max(xx[1:Ntip(phy)]) - xx[1:Ntip(phy)] > tol;
if (any(aa)) {
return(phy$tip.label[which(aa)]);
} else {
return(NULL);
}
}
## EB transform
## part of the code used in fit_t_pl
## J. Clavel
transform_EB <- function(phy, beta, sigmasq){
# parents & descent
parent <- phy$edge[,1]
descendent <- phy$edge[,2]
if (beta!=0){
distFromRoot <- node.depth.edgelength(phy)
phy$edge.length = (exp(beta*distFromRoot[descendent])-exp(beta*distFromRoot[parent]))/beta
}
if(sigmasq!=0) phy$edge.length <- phy$edge.length*sigmasq
return(phy)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.