R/getSampleDCA.R

In paleoAM: Simulating Assemblage Models of Abundance for the Fossil Record

#' Get the Detrended correspondence Analysis (DCA) Score Value for A Single Sample
#' 
#' For a single simulated assemblage sample, projects its location in DCA space defined by the original abundance data.

#' @details
#' Detrended correspondence analysis (DCA) is a common method 
#' for producing an ordination of ecological data. 
#' It isn't the only such method.

#' @param simSample The assemblage data for a single sample (presumably from a simulation).

#' @param origAbundData The original matrix of abundance data, to be used to
#'  project the simulated data into the same detrended correspondence analysis (DCA) space.

#' @param useTransformedRelAbundance Should the DCA be analyzed 

#' @param projectIntoOrigDCA Should the new simulated data be projected in the DCA generated by analyzing the original data? This is \code{TRUE} by default, which is what most users will likely use, as it is the preferable way to consider how the new simulated data relates to the original data.

#' @param returnDCAforOrigAndSim Should the DCA score values for both the new simulated data and the original abundance data be returned? Default is \code{FALSE} as projecting the new data into the original DCA space means the original data should never be meaningfully different the original score values.

#' @param whichAxes Which dimensional score from the DCA should be used? By default this is 1. Unclear under what circumstances one would ever use a value other than 1, though, as detrending the correspondence analysis causes distortion along all scores other than the first axis of the ordination. Only the first score can be returned when \code{projectIntoOrigDCA = TRUE}.

#' @param powerRootTransform The power-root transform to be used on the abundance data before applying the DCA. By default this is 1, which means the data is not transformed at all. Note that the power-root transform is only performed if \code{useTransformedRelAbundance = TRUE}.

#' @return
#' A vector, containing either a single value, the DCA score value of the simulated sample, when \code{returnDCAforOrigAndSim = FALSE}, or a vector of the DCA scores for the original data and 

# @aliases

#' @seealso
#' This function is ultimately just a wrapper for using \code{decorana} in package \code{vegan}.

# @references

# @examples


#' @name getSampleDCA
#' @rdname getSampleDCA
#' @export
getSampleDCA <- function(
            simSample,
            origAbundData,
            useTransformedRelAbundance = TRUE,
            projectIntoOrigDCA = TRUE,
            returnDCAforOrigAndSim = FALSE,
            whichAxes = 1,
            powerRootTransform = 1
            ){

        if(projectIntoOrigDCA){
                
            if(length(whichAxes) != 1){
                stop("Cannot return any axis other than DCA-1 when 'projectIntoOrigDCA = TRUE' ")
                }
            if(whichAxes != 1){
                stop("Cannot return any axis other than DCA-1 when 'projectIntoOrigDCA = TRUE' ")
                }
            
            abundanceTable <- origAbundData
            ## Doing a DCA just on the original data
            if(useTransformedRelAbundance){
                # root transform
                abundanceTable <- powerRootTransformFun(abundanceTable, 
                    powerRootTransform = powerRootTransform)   
            }else{
                if(powerRootTransform != 1){stop(
                    "Non default value specified for powerRootTransform but transformed relative abundance is not enabled??")}
                }

            # do the DCA!
            dcaOut <- vegan::decorana(abundanceTable)
            
            newSample <- rbind(origAbundData,simSample)
            if(useTransformedRelAbundance){        
                newSample <- powerRootTransformFun(newSample, 
                    powerRootTransform = powerRootTransform)   
                }

            # project the new data in the DCA
            newSampleDCA <- predictDCA_sites_dca1(
                oldDCA = dcaOut,
                newSample = newSample
                )
            
            if(returnDCAforOrigAndSim){
                # get DCA value for sim assemblages AND empirical data
                output <- newSampleDCA[ , 1]
            }else{
                # get DCA value for the simulated mixed assemblage
                output <- newSampleDCA[nrow(newSample), 1]
                }
            
        }else{
            # take a simulated sample, combine with original data,
                # apply DCA and get the DCA-1 value

            # combined picked sample with original abundance data
            abundanceTable <- rbind(origAbundData, simSample)
            #    
            if(useTransformedRelAbundance){        
                # Transform the abundance table to relative abundances, 
                    # (and get the pair-wise Bray-Curtis distances... wait no)
                    # recalculate relative abundance table
                #relAbundanceTable <- t(apply(abundanceTable, 1, function(x) x/sum(x)))
                # root transform?
                abundanceTable <- powerRootTransformFun(abundanceTable, 
                    powerRootTransform = powerRootTransform)  
                }
            ## Doing a DCA on the Simulated Data
            dcaOut <- vegan::decorana(abundanceTable)
            dcaOut <- vegan::scores(dcaOut)
            
            if(returnDCAforOrigAndSim){
                # get DCA value for sim assemblages AND empirical data
                output <- dcaOut[ , whichAxes]
            }else{
                # get DCA value for the simulated mixed assemblage
                output <- dcaOut[nrow(abundanceTable), whichAxes]
                }
            }
        
        # get bray-curtis distances
        #bcdist <- vegan::vegdist(relAbundanceTable, method = "bray")

        return(output)
        }

powerRootTransformFun <- function(abundTable, powerRootTransform){
    # convert to relative abundance
    out <- apply(abundTable, 1, function(x) x/sum(x))
    out <- t(out) ^ (1/powerRootTransform)  
    return(out)
    }