EvolQG: Evolutionary Quantitative Genetics

#' 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)
}