R/pdf_ranges.R

In mchevalier2/crestr: A Probabilistic Approach to Reconstruct Past Climates Using Palaeoecological Datasets

#' Calculate the climate tolerance of the taxa from their pdfs.
#'
#' Calculate the climate tolerance of the taxa from their pdfs.
#'
#' @inheritParams plot.crestObj
#' @param x A \code{\link{crestObj}} generated by either the \code{\link{crest.calibrate}},
#'        \code{\link{crest.reconstruct}} or \code{\link{crest}} functions.
#' @param climate Climate variables to be used to generate the plot. By default
#'        all the variables are included.
#' @param taxanames A list of taxa to use for the plot (default is all the
#'        recorded taxa).
#' @param orderby A string ('name', or one of the climate variables) to sort
#'        the output table
#' @return The set of coordinates ext projected in equal earth.
#' @export
#' @examples
#'
#' data(reconstr)
#' pdf_ranges(reconstr, climate='bio1')
#' pdf_ranges(reconstr, climate='bio12', orderby='bio1', uncertainties=c(0.2, 0.6, 0.95))
#'
pdf_ranges <- function( x,
                climate = x$parameters$climate,
                taxanames = x$input$taxa.name,
                uncertainties = x$parameters$uncertainties,
                orderby = NA
              ) {

    if(base::missing(x)) x

    if (is.crestObj(x)) {
        test <- is.na(x$modelling$pdfs)
        if( test[1] & length(test) == 1 ) {
            stop('The crestObj requires the climate space to be calibrated. Run crest.calibrate() on your data.\n')
        }

        for(clim in climate){
            if(! clim %in% x$parameters$climate) {
                stop(paste0("The climate variable should be one of the following options: '", paste(x$parameters$climate, collapse="', '"), "'.\n"))
                return(invisible(NA))
            }
        }

        if(!is.na(orderby)) {
            if(! orderby %in% c('name', climate)) {
                warning(paste0("The value orderby = '", orderby, "' is not accepted. It should be 'name' or one of the climate variables. The data were not sorted.\n"))
                orderby <- NA
            }
        }

        err <- c()
        for(tax in taxanames) {
            if(! tax %in% x$input$taxa.name) err <- c(err, tax)
        }
        if(length(err) > 0) {
            stop(paste0("The following taxa names are not available in your dataset: '", paste(err, collapse="', '"), "'\n\n"))
            return(invisible(NA))
        }

        rs <- list()
        for(uncer in uncertainties) {
            res <- as.data.frame(matrix(NA, ncol=3*length(climate), nrow=length(taxanames)))
            rownames(res) <- taxanames  ;
            for(clim in 1:length(climate)) {
                colnames(res)[(3*(clim-1)+1):(3*clim)] <- paste(climate[clim], c('tol_inf', 'tol_sup', 'range'), sep='_')
                for(tax in taxanames) {
                    if(!is.null(x$modelling$pdfs[[tax]])) {
                        pdf <- x$modellin$pdfs[[tax]][[climate[clim]]]$pdfpol
                        oo <- rev(order(pdf))
                        tmp1 <- x$modelling$xrange[[climate[clim]]][oo]
                        tmp2 <- pdf[oo]
                        oo <- order(tmp1)
                        pdfter <- cumsum(tmp2 / sum(tmp2))[oo]
                        w <- which(pdfter <= uncer)
                        res[tax, (3*(clim-1)+1):(3*(clim-1)+2)] <- x$modelling$xrange[[climate[clim]]][c(w[1], w[length(w)])]
                    }
                    res[, 3*clim] <- res[, (3*(clim-1)+2)] - res[, (3*(clim-1)+1)]
                }
            }
            if(!is.na(orderby)) {
                if(orderby == 'name') res <- res[order(taxanames), ]
                if(orderby %in% climate) res <- res[order(res[, 3*which(climate == orderby)]), ]
            }
            rs[[paste0('Range = ', 100*uncer, '%')]] <- res
        }
    } else {
        cat('This function only works with a crestObj.\n\n')
    }
    rs
}