Nothing
#' Calculate mean values for various matrix statistics
#'
#' Calculates: Mean Squared Correlation, ICV, Autonomy, ConditionalEvolvability, Constraints, Evolvability, Flexibility, Pc1Percent, Respondability.
#' @aliases Autonomy ConditionalEvolvability Constraints Evolvability Flexibility Pc1Percent Respondability
#' @param cov.matrix A covariance matrix
#' @param iterations Number of random vectors to be used in calculating the stochastic statistics
#' @param full.results If TRUE, full distribution of statistics will be returned.
#' @param parallel if TRUE computations are done in parallel. Some foreach backend must be registered, like doParallel or doMC.
#' @return dist Full distribution of stochastic statistics, only if full.resuts == TRUE
#' @return mean Mean value for all statistics
#' @export
#' @references Hansen, T. F., and Houle, D. (2008). Measuring and comparing evolvability
#' and constraint in multivariate characters. Journal of evolutionary
#' biology, 21(5), 1201-19. doi:10.1111/j.1420-9101.2008.01573.x
#' @author Diogo Melo Guilherme Garcia
#' @importFrom Matrix Matrix solve chol
#' @examples
#' cov.matrix <- cov(iris[,1:4])
#' MeanMatrixStatistics(cov.matrix)
#'
#' \dontrun{
#' #Multiple threads can be used with some foreach backend library, like doMC or doParallel
#' library(doMC)
#' registerDoMC(cores = 2)
#' MeanMatrixStatistics(cov.matrix, parallel = TRUE)
#' }
#' @keywords Autonomy
#' @keywords ConditionalEvolvability
#' @keywords Constraints
#' @keywords Evolvability
#' @keywords Flexibility
#' @keywords Pc1Percent
#' @keywords Respondability
MeanMatrixStatistics <- function (cov.matrix, iterations = 1000, full.results = FALSE, parallel = FALSE) {
matrix.stat.functions = list ('respondability' = Respondability,
'evolvability' = Evolvability,
'conditional.evolvability' = ConditionalEvolvability,
'autonomy' = Autonomy,
'flexibility' = Flexibility,
'constraints' = Constraints)
num.traits <- dim (cov.matrix) [1]
beta.mat <- array (rnorm (num.traits * iterations), c(num.traits, iterations))
beta.mat <- apply (beta.mat, 2, Normalize)
iso.vec <- Normalize (rep(1, num.traits))
null.dist <- abs (t (iso.vec) %*% beta.mat)
null.dist <- sort (null.dist)
crit.value <- null.dist [round (0.95 * iterations)]
if(full.results) cat ('critical value: ', crit.value, '\n')
MatrixStatisticsMap <- function (CurrentFunc) CurrentFunc(cov.matrix = cov.matrix, beta.mat = beta.mat)
stat.dist <- t(laply (matrix.stat.functions, MatrixStatisticsMap, .parallel = parallel))
stat.dist <- cbind (stat.dist, null.dist)
colnames (stat.dist) <- c('respondability',
'evolvability',
'conditional.evolvability',
'autonomy',
'flexibility',
'constraints',
'null.dist')
stat.mean <- colMeans (stat.dist[,-7])
integration <- c (CalcR2 (cov.matrix),
Pc1Percent (cov.matrix),
CalcICV(cov.matrix),
CalcEigenVar(cov.matrix, sd=TRUE, rel=FALSE))
names (integration) <- c ('MeanSquaredCorrelation', 'pc1.percent', 'ICV', 'EigenSd')
stat.mean <- c (integration, stat.mean)
if(full.results)
return (list ('dist' = stat.dist, 'mean' = stat.mean))
else
return (stat.mean)
}
#' @export
#' @importFrom Matrix Matrix solve chol nearPD diag
Autonomy <- function (cov.matrix, beta.mat = NULL, iterations = 1000){
cov.matrix <- Matrix(cov.matrix)
num.traits <- dim (cov.matrix) [1]
if(is.null(beta.mat)){
beta.mat <- array (rnorm (num.traits * iterations), c(num.traits, iterations))
beta.mat <- apply (beta.mat, 2, Normalize)
}
cov.matrix = tryCatch({cv = chol(cov.matrix); cov.matrix}, error = function(cond){
warning("matrix is singular, can't compute autonomy directly. Using nearPD, results could be wrong")
cov.matrix <- nearPD(cov.matrix)[[1]]
chol(cov.matrix)
})
(1/diag(t (beta.mat) %*% solve (cov.matrix, beta.mat))) / diag(t (beta.mat) %*% cov.matrix %*% beta.mat)
}
#' @export
#' @importFrom Matrix Matrix solve chol nearPD diag
ConditionalEvolvability <- function (cov.matrix, beta.mat = NULL, iterations = 1000){
cov.matrix <- Matrix(cov.matrix)
num.traits <- dim (cov.matrix) [1]
if(is.null(beta.mat)){
beta.mat <- array (rnorm (num.traits * iterations), c(num.traits, iterations))
beta.mat <- apply (beta.mat, 2, Normalize)
}
cov.matrix = tryCatch({chol(cov.matrix); cov.matrix}, error = function(cond){
warning("matrix is singular, can't compute conditional evolvability directly. Using nearPD, results could be wrong")
cov.matrix <- nearPD(cov.matrix)[[1]]
chol(cov.matrix)
})
return (1/diag(t (beta.mat) %*% solve (cov.matrix, beta.mat)))
}
#' @export
Constraints <- function (cov.matrix, beta.mat = NULL, iterations = 1000){
num.traits <- dim (cov.matrix) [1]
if(is.null(beta.mat)){
beta.mat <- array (rnorm (num.traits * iterations), c(num.traits, iterations))
beta.mat <- apply (beta.mat, 2, Normalize)
}
PC1 <- eigen(cov.matrix)$vectors[,1]
Cb <- apply(cov.matrix %*% beta.mat, 2, Normalize)
abs(apply(Cb, 2, `%*%`, PC1))
}
#' @export
Evolvability <- function (cov.matrix, beta.mat = NULL, iterations = 1000){
num.traits <- dim (cov.matrix) [1]
if(is.null(beta.mat)){
beta.mat <- array (rnorm (num.traits * iterations), c(num.traits, iterations))
beta.mat <- apply (beta.mat, 2, Normalize)
}
diag(t(beta.mat) %*% cov.matrix %*% beta.mat)
}
#' @export
Flexibility <- function (cov.matrix, beta.mat = NULL, iterations = 1000){
num.traits <- dim (cov.matrix) [1]
if(is.null(beta.mat)){
beta.mat <- array (rnorm (num.traits * iterations), c(num.traits, iterations))
beta.mat <- apply (beta.mat, 2, Normalize)
}
Cb <- apply(cov.matrix %*% beta.mat, 2, Normalize)
diag(t (beta.mat) %*% Cb)
}
#' @export
Pc1Percent <- function (cov.matrix) return (eigen (cov.matrix)$values [1] / sum (eigen (cov.matrix)$values))
#' @export
Respondability <- function (cov.matrix, beta.mat = NULL, iterations = 1000) {
num.traits <- dim (cov.matrix) [1]
if(is.null(beta.mat)){
beta.mat <- array (rnorm (num.traits * iterations), c(num.traits, iterations))
beta.mat <- apply (beta.mat, 2, Normalize)
}
apply(cov.matrix %*% beta.mat, 2, Norm)
}
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.