R/extractCorrelation.R

In Caetanods/ratematrix: Bayesian Estimation of the Evolutionary Rate Matrix

##' Function extracts the posterior distribution of evolutionary correlation among traits.
##'
##' Returns a list with length equal to the number of regimes. Each list element is composed by a matrix with trait correlation types in the columns and the evolutionary correlations for each sample at the rows. \cr
##' \cr
##' One can plot the correlation values using boxplots and compare their distribution. Pairwise statistical tests across the samples is also possible.
##' @title Extract the posterior distribution of evolutionary correlation
##' @param post a posterior distribution object as returned by the function 'readMCMC' or a merged posterior generated by 'mergePosterior'.
##' @return a list with the posterior distribution of evolutionary correlations among traits. If the data is a 2x2 matrix then the object will be a matrix and each regime will be a column of this matrix.
##' @author Daniel Caetano and Luke Harmon
##' @export
##' @examples
##' \donttest{
##' data( centrarchidae )
##' handle <- ratematrixMCMC(data=centrarchidae$data, phy=centrarchidae$phy.map
##'                          , gen=10000, dir=tempdir())
##' posterior <- readMCMC(handle, burn = 0.2, thin = 1)
##' ## Get the correlations:
##' cor.list <- extractCorrelation(post = posterior)
##' ## Plot the results:
##' class( cor.list ) ## In this case we have a matrix.
##' boxplot(cor.list, names = c("Regime 1", "Regime 2"))
##' ## Same example with more traits.
##' data( anoles )
##' handle <- ratematrixMCMC(data=anoles$data[,1:3], phy=anoles$phy
##'                          , gen=5000, dir=tempdir())
##' anole_post <- readMCMC(handle, burn = 0.2, thin = 1)
##' ## Get the correlations:
##' cor.anole <- extractCorrelation(post = anole_post)
##' ## Plot the results:
##' class( cor.anole ) ## In this case we have a list!
##' names( cor.anole ) ## Each element is a regime.
##' ## We can plot the regimes in separate.
##' boxplot(cor.anole$island)
##' boxplot(cor.anole$mainland)
##' }
extractCorrelation <- function(post){
    ## Check the class of the input.
    if( !inherits(post, what = c("ratematrix_multi_chain", "ratematrix_single_chain") ) ) stop("Argument 'post' need to be a posterior distribution object. See 'readMCMC' help page.")
  
    if( inherits(post, what = "ratematrix_single_chain" ) ){
        ## Get some information from the data.
        n.regimes <- 1
        regime.names <- NULL
        cor.list <- list()
        cor.list[[1]] <- lapply(post$matrix, cov2cor )
    } else{
        ## Get some information from the data.
        n.regimes <- length( post$matrix )
        regime.names <- names( post$matrix )
        cor.list <- lapply(1:n.regimes, function(x) lapply(post$matrix[[x]], cov2cor ) )

    }

    ## These objects are the same for any number of regimes:
    traits <- colnames( post$root )
    n.samples <- length(cor.list[[1]])
    n.traits <- ncol( cor.list[[1]][[1]] )
    n.upper <- sum( upper.tri(cor.list[[1]][[1]]) )

    ## Create storage objects.
    out.list <- lapply(1:n.regimes, function(x) matrix(nrow = n.samples, ncol = n.upper) )
    upper.id <- matrix(nrow = n.upper, ncol = 2)
    upper.name <- vector("character", length = n.upper)

    ## Find the indexes for the upper triangular elements.
    ## Also generate informative names.
    count <- 1
    for( n.row in 1:n.traits ){
        for( n.col in 1:n.traits ){
            if( n.col > n.row ){
                upper.id[count,] <- c(n.row,n.col)
                upper.name[count] <- paste0(traits[n.col], "_x_", traits[n.row])
                count <- count + 1
            }
        }
    }

    ## Loop through the thing.
    for( regime in 1:n.regimes ){
        for( upper in 1:n.upper ){
            out.list[[regime]][,upper] <- sapply(1:n.samples, function(y) cor.list[[regime]][[y]][upper.id[upper,1],upper.id[upper,2]])
        }
        colnames( out.list[[regime]] ) <- upper.name
    }
    names( out.list ) <- regime.names

    ## Return a list-type object if number of traits is larger than 2.
    if( n.traits > 2 ){
        return( out.list )
    } else{
        out.mat <- matrix(nrow = n.samples, ncol = n.regimes)
        for( i in 1:n.regimes){
            out.mat[,i] <- as.numeric( out.list[[i]] )
        }
        colnames( out.mat ) <- regime.names
        return( out.mat )
    }
    
}