#' Comparisons of Effect Sizes from Modularity Analyses
#'
#' Function performs an analysis to compare the effect sizes of two or more CR effects
#'
#' The function statistically compares the effect sizes of two or more CR analyses.  Typically, this
#' function might be used to compare levels of modularity between two or more samples, each measuring the degree 
#' of morphological modularity in each.  Alternatively, the approach can compare the degree of modular signal as 
#' expressed by alternative modular hypotheses. 
#' 
#' The analysis calculates effect sizes as standard deviates, z, and performs two-sample z-tests, using the pooled 
#' standard error from the sampling distributions of the CR analyses. The method follows that of Adams and Collyer (2016) used 
#' to compare patterns of integration via PLS across datasets
#'  
#' To use this function, simply perform \code{\link{modularity.test}}, or \code{\link{phylo.modulatiry}} on as many samples or 
#' alternative modular hypotheses as desired.  Any number of objects of class CR can be input.
#' 
#' SOMEWHERE NEED TO TALK ABOUT CR.null = FALSE/TRUE.  IF TRUE, for including a 'no modularity' model in the comparisons
#' 
#' @param ... saved analyses of class CR
#' @param CR.null A logical value to indicate whether a Null CR model (no modularity) should also be included in analysis.
#' @param two.tailed A logical value to indicate whether a two-tailed test (typical and default) should be performed.
#' @keywords analysis
#' @export
#' @author Dean Adams and Michael Collyer
#' @return An object of class compare.CR, returns a list of the following
#' \item{sample.z}{A vector of effect sizes for each sample.}
#' \item{sample.r.sd}{A vector of standard deviations for each sampling distribution.}
#' \item{pairwise.z}{A matrix of pairwise, two-sample z scores between all pairs of effect sizes.}
#' \item{pairwise.p}{A matrix of corresponding P-values.}
#' @references Adams, D.C. and M.L. Collyer. 2016.  On the comparison of the strength of morphological integration across morphometric 
#' datasets. Evolution. 70:2623-2631.
#' @references Adams, D.C.  PAPER TO BE WRITTEN LATER
#' @examples
#' 
#' TO BE ADDED LATER
#' 
 compare.CR <- function(...,CR.null = TRUE, two.tailed = TRUE){
   dots <- list(...)
   tails <- if(two.tailed) 2 else 1
   if(length(dots) == 1) n <- length(dots[[1]]) else n <- length(dots)
   if(n == 1) stop("At least two objects of class CR are needed")
   if(length(dots) == 1) {
     list.names <- names(dots[[1]]) 
     dots <- lapply(1:n, function(j) dots[[1]][[j]])
     names(dots) <- list.names
     } else list.names <- names(dots)
   
   if(length(dots) < 2) stop("At least two objects of class CR are needed")
     
   if(is.null(list.names)) list.names <- paste("CR", 1:n, sep = ".")
   names(dots) <- list.names
        
   is.CR <- function(x) class(x) == "CR"
   sdn <- function(x) sqrt(sum((x-mean(x))^2)/length(x))
   list.check <- sapply(1:length(dots), function(j) any(is.CR(dots[[j]])))
   if(any(list.check == FALSE)) stop("Not all objects are class CR")
   k <- length(list.check)
   if(is.null(list.names)) list.names <- as.list(substitute(list(...)))[-1L]
   k.combn <- combn(k,2)
   list.sds <- sapply(1:k, function(j) sdn(scale(dots[[j]]$random.CR[-1])))
   list.zs <- sapply(1:k, function(j) RRPP:::effect.size(dots[[j]]$random.CR, center=TRUE))  #NOTE: pull out RRPP::: once in geomorph
   if (CR.null == TRUE){
      k <- k + 1
      k.combn <- combn(k,2)
      list.sds <- c(0,list.sds)
      list.zs <- c(0,list.zs)
   }
   #Changed this: no 'abs', and forced min as r1, then other as r2
   z12 <- sapply(1:ncol(k.combn), function(j){
   a <- k.combn[1,j]; b <- k.combn[2,j]
     r1 <- list.zs[a]; r2 <- list.zs[b]; sd1 <- list.sds[a]; sd2 <- list.sds[b]
     (r1-r2)/sqrt(sd1^2+sd2^2)
   })
   z12.p <- sapply(1:length(z12), function(j) pnorm(abs(z12[[j]]), lower.tail = FALSE)* tails)  
   d <- rep(0,k); names(d) <- list.names
   D <-dist(d)
   z12.pw <- p12.pw <- D
   for(i in 1:length(z12)) z12.pw[i] <-z12[i]
   for(i in 1:length(z12)) p12.pw[i] <-z12.p[i]
   names(list.zs) <- names(list.sds) <-list.names
   pairwise.z <- as.matrix(z12.pw)
   pairwise.P <- as.matrix(p12.pw)
   diag(pairwise.P) <- 1
   if (CR.null == TRUE){
      names(list.zs) <- names(list.sds) <-c("No_Modules",list.names)
      pairwise.z <- as.matrix(z12.pw); rownames(pairwise.z) <- colnames(pairwise.z)<-c("No_Modules",list.names)
      pairwise.P <- as.matrix(p12.pw); rownames(pairwise.P) <- colnames(pairwise.P)<-c("No_Modules",list.names)
      diag(pairwise.P) <- 1
   }
   comment <- c("NOTE: more negative effects represent stronger modular signal!")
   out <- list(comment=comment, sample.z = list.zs,
               sample.r.sd = list.sds,
               pairwise.z = pairwise.z,
               pairwise.P = pairwise.P)
   class(out) <- "compare.CR"
   out
 }
 