Nothing
#' Compare matrices via Selection Response Decomposition
#'
#' Based on Random Skewers technique, selection response vectors are
#' expanded in direct and indirect components by trait and compared via
#' vector correlations.
#' @param cov.x Covariance matrix being compared. cov.x can be a matrix or a list.
#' @param cov.y Covariance matrix being compared. Ignored if cov.x is a list.
#' @param iterations Number of random vectors used in comparison
#' @param parallel if TRUE computations are done in parallel. Some foreach back-end must be registered, like doParallel or doMC.
#' @param ... additional parameters passed to other methods
#' @param x Output from SRD function, used in plotting
#' @param matrix.label Plot label
#' @details Output can be plotted using PlotSRD function
#' @return List of SRD scores means, confidence intervals, standard
#' deviations, centered means e centered standard deviations
#' @return pc1 scored along the pc1 of the mean/SD correlation matrix
#' @return model List of linear model results from mean/SD correlation. Quantiles, interval and divergent traits
#' @note If input is a list, output is a symmetric list array with pairwise comparisons.
#' @export
#' @rdname SRD
#' @references Marroig, G., Melo, D., Porto, A., Sebastiao, H., and Garcia, G.
#' (2011). Selection Response Decomposition (SRD): A New Tool for
#' Dissecting Differences and Similarities Between Matrices. Evolutionary
#' Biology, 38(2), 225-241. doi:10.1007/s11692-010-9107-2
#' @author Diogo Melo, Guilherme Garcia
#' @seealso \code{\link{RandomSkewers}}
#' @examples
#' cov.matrix.1 <- cov(matrix(rnorm(30*10), 30, 10))
#' cov.matrix.2 <- cov(matrix(rnorm(30*10), 30, 10))
#' colnames(cov.matrix.1) <- colnames(cov.matrix.2) <- sample(letters, 10)
#' rownames(cov.matrix.1) <- rownames(cov.matrix.2) <- colnames(cov.matrix.1)
#' srd.output <- SRD(cov.matrix.1, cov.matrix.2)
#'
#' #lists
#' m.list <- RandomMatrix(10, 4)
#' srd.array.result = SRD(m.list)
#'
#' #divergent traits
#' colnames(cov.matrix.1)[as.logical(srd.output$model$code)]
#'
#' #Plot
#' plot(srd.output)
#'
#' ## For the array generated by SRD(m.list) you must index the idividual positions for plotting:
#' plot(srd.array.result[1,2][[1]])
#' plot(srd.array.result[3,4][[1]])
#'
#' \dontrun{
#' #Multiple threads can be used with some foreach backend library, like doMC or doParallel
#' library(doMC)
#' registerDoMC(cores = 2)
#' SRD(m.list, parallel = TRUE)
#' }
#' @keywords SRD
#' @keywords RandomSkewers
#' @keywords selectionresponsedecomposition
SRD <- function (cov.x, cov.y, ...) UseMethod("SRD")
#' @rdname SRD
#' @method SRD default
#' @export
SRD.default <- function (cov.x, cov.y, iterations = 1000, ...) {
size <- dim (cov.x)[1]
if(iterations < 2*size) stop("Iterations must be larger than twice the number of traits.")
r2s <- array (0, c(size,iterations))
beta <- apply (array (rnorm (size*iterations, mean = 0, sd = 1), c(size,iterations)), 2, Normalize)
for (I in 1:iterations){
beta.matrix <- diag (beta[,I])
dz1 <- apply (cov.x %*% beta.matrix, 1, Normalize)
dz2 <- apply (cov.y %*% beta.matrix, 1, Normalize)
r2s[,I] <- colSums (dz1 * dz2)
}
# results
mean.r2 <- apply (r2s, 1, mean)
sample.conf <- function (x, lower = TRUE){
ox <- x[order(x)]
lox <- length (ox)
if (lower)
crit <- round (0.025 * lox)
else
crit <- round (0.975 * lox)
return (ox[crit])
}
low.r2 <- apply (r2s, 1, sample.conf, lower = TRUE)
up.r2 <- apply (r2s, 1, sample.conf, lower = FALSE)
sd.r2 <- apply (r2s,1,sd)
cmean.r2 <- scale (mean.r2, scale = FALSE)
csd.r2 <- scale (sd.r2, scale = FALSE)
cent <- cbind (cmean.r2,csd.r2)
pca.cent <- princomp (cent, cor = TRUE, scores = TRUE)
pc1 <- pca.cent$scores[,1]
if (pca.cent$loadings[1,1] < 0)
pc1 <- - pc1
pc1 <- pc1 / sd (pc1)
pc1.quant <- quantile (pc1,
probs = c (1,5,10,20,25,30,40,50,60,70,75,80,90,95,99)/100,
names = FALSE)
pc1.int <- - 1.96 / sqrt (length (pc1))
pc1.sig <- ifelse (pc1 < pc1.int, 1, 0)
model <- list ("quantiles" = pc1.quant,
"interval" = pc1.int,
"code" = pc1.sig)
output <- cbind (mean.r2, low.r2, up.r2, sd.r2, cmean.r2, csd.r2)
colnames (output) <- c("ARC","IC-","IC+","SD","CMEAN","CSD")
rownames (output) <- rownames (cov.x)
srd.out <-
list ("output" = output,
"pc1" = pc1,
"model" = model)
class(srd.out) <- 'SRD'
return (srd.out)
}
#' @rdname SRD
#' @method SRD list
#' @export
SRD.list <- function (cov.x, cov.y = NULL, iterations = 1000, parallel = FALSE, ...){
n.matrix <- length(cov.x)
if(is.null(names(cov.x))) {names(cov.x) <- 1:n.matrix}
matrix.names <- names (cov.x)
CompareToN <- function(n) llply(cov.x[(n+1):n.matrix],
function(x) {SRD(x, cov.x[[n]], iterations = iterations)},
.parallel = parallel)
comparisons <- alply(1:(n.matrix-1), 1, CompareToN, .parallel = parallel)
corrs <- array(list(), c(n.matrix, n.matrix))
for(i in 1:(n.matrix-1)){
corrs[i,(i+1):n.matrix] <- comparisons[[i]]
}
corrs[lower.tri(corrs)] <- t(corrs)[lower.tri(corrs)]
colnames(corrs) <- rownames(corrs) <- names(cov.x)
return (corrs)
}
#' @rdname SRD
#' @export
plot.SRD <- function (x, matrix.label = "", ...) {
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
output <- x
layout (array (c(1,1,2,2),c(2,2)))
par (mar = c(4.0, 4.0, 8, 0.4))
mean.r2 <- output$output[,1]
low.r2 <- output$output[,2]
up.r2 <- output$output[,3]
c.mean.r2 <- output$output[,5]
c.sd.r2 <- output$output[,6]
if (is.null (rownames (output$output)))
dists <- 1:length (mean.r2)
else
dists <- rownames (output$output)
### plot scores
### b l t r
pc.pch <- output$model$code + 20
plot (mean.r2, type = "p", lty = 2, pch = pc.pch,
ylab = "", xlab = "", xaxt = "n", ylim = c(-1,1), ...)
for (i in 1:length (mean.r2))
{
abline (v = i, lty = 1, col = rgb (0.8,0.8,0.8))
}
arrows (x0 = 1:length (mean.r2),
y0 = low.r2, y1 = up.r2,
angle = 90, length = 0.05, code = 3)
abline (h = mean (mean.r2), lty = 3)
axis (3, 1:length (mean.r2), dists, las = 2, cex.axis = 1.3)
mtext (side = 2, at = mean (mean.r2), text = round (mean (mean.r2), 2), las = 2, cex=1.8)
### plot av sd
### b l t r
par (mar = c(4.0, 0.0, 8, 4.6))
output$model$code[output$model$code == 1] <- 1.2
output$model$code[output$model$code == 0] <- 1
pc.cex <- output$model$code
plot (c.sd.r2 ~ c.mean.r2, pch = pc.pch, xlab = "", ylab = "",
yaxt = "n", main = matrix.label, cex.main= 3, cex.axis=1.3, ...)
abline (v = 0, lty = 2)
abline (h = 0, lty = 2)
text (c.mean.r2, c.sd.r2, labels = dists, pos = 4, cex = pc.cex)
axis (4, las = 2)
}
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.