R/plot.combinedPDFS.R

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

#' Plot representing how the \code{pdfs} combine to produce the reconstruction.
#'
#' Plot representing how the \code{pdfs} combine to produce the reconstruction.
#'
#' @inheritParams plot.crestObj
#' @param x A \code{\link{crestObj}} generated by the
#'        \code{\link{crest.reconstruct}} or \code{\link{crest}} functions.
#' @param ages An age range to subset the samples to plot. By default, all the
#'        samples will be selected.
#' @param samples The list of sample indexes for which the plot should be
#'        plotted. All samples will be plotted by default.
#' @param climate The climate variable to use to plot the variable. Default is
#'        first variable (\code{x$parameters$climate\[1\]}).
#' @param optima A boolean to indicate whether to plot the optimum (\code{TRUE})
#'        or the mean (\code{FALSE}) estimates.
#' @param only.present A boolean to only add the names of the taxa recorded in
#'        the sample (default \code{FALSE}).
#' @param only.selected A boolean to only add the names of the selected taxa
#'        (default \code{FALSE}).
#' @param xlim The climate range to plot the pdfs on. Default is the full range
#'        used to fit the \code{pdfs} (x$modelling$xrange).
#' @param col A range of colour values to colour the \code{pdfs}. Colours will
#'        be recycled to match the number of taxa.
#' @param filename An absolute or relative path that indicates where the diagram
#'        should be saved. Also used to specify the name of the file. Default:
#'        the file is saved in the working directory under the name
#'        \code{'samplePDFs.pdf'}.
#' @param width The width of the output file in inches (default 7.48in ~ 19cm).
#' @param height The height of the output file in inches (default 5in ~ 12.7cm).
#' @return No return value, this function is used to plot.
#' @export
#' @examples
#' \dontrun{
#'   data(crest_ex)
#'   data(crest_ex_pse)
#'   data(crest_ex_selection)
#'   reconstr <- crest(
#'     df = crest_ex, pse = crest_ex_pse, taxaType = 0,
#'     climate = c("bio1", "bio12"), bin_width = c(2, 20),
#'     shape = c("normal", "lognormal"),
#'     selectedTaxa = crest_ex_selection, dbname = "crest_example",
#'     leave_one_out = FALSE
#'    )
#' }
#' ## example using pre-saved reconstruction obtained with the previous command.
#' data(reconstr)
#' plot_combinedPDFs(reconstr, ages=c(6,9), climate='bio12')
#'
plot_combinedPDFs <- function( x, ages=range(x$inputs$x),
                               samples=which(x$inputs$x >= min(ages) & x$inputs$x <= max(ages)),
                               climate=x$parameters$climate[1],
                               optima=TRUE, xlim=NA, only.present=FALSE, only.selected=FALSE,
                               col=crestr::colour_theme(1),
                               save=FALSE, filename='samplePDFs.pdf',
                               as.png = FALSE, png.res=300,
                               width=7.48, height = 5
                              ) {

    if(base::missing(x)) x

    if (is.crestObj(x)) {
        if (length(x$reconstructions) == 0 ) {
            stop('The crestObj requires the fossil data. Run crest.reconstruct() on your data.\n\n')
            return(invisible())
        }

        if (length(ages) != 2 ) {
            stop("The parameter 'ages' should have a length of 2.\n\n")
            return(invisible())
        }

        if (length(climate) > 1 ) {
            warning('The function only works for one variable at a time. Only the first value was used.\n')
            climate <- climate[1]
        }

        if (max(samples) > length(x$inputs$x) ) {
            warning(paste0('The record contains ', length(x$inputs$x), " samples. The parameter 'samples' was adapted accordingly.\n"))
            samples <- samples[samples <= length(x$inputs$x)]
        }
        if(length(samples) == 0) {
            stop(paste0("You did not select samples to plot. Assign a vector of values comprised between 1 and ", length(x$inputs$x),"\n\n"))
            return(invisible(NA))
        }

        err <- c()
        for(clim in climate) {
            if(! clim %in% x$parameters$climate) err <- c(err, clim)
        }
        if(length(err) > 0) {
            stop(paste0("The following variable is not available in your crestObj: '", err, "'\n\n"))
            return(invisible(NA))
        }

        climVar <- accClimateVariables()
        new_clim <- climate
        if (!(climate %in% climVar[, 1] | climate %in% climVar[, 2])) {
            stop(paste0("The variable '", climate, "' is not an accepted value. Check the list of accepted values using 'accClimateVariables()'.\n\n"))
        } else {
            if (suppressWarnings(!is.na(as.numeric(climate)))) {
                new_clim <- as.character(climVar[which(climVar[, 1] == as.numeric(climate)), 2])
            }
        }

        if (is.na(xlim)[1]) {
            xlim <- range(x$modelling$xrange[[climate]])
        }

        if(!save) {
            par_usr <- graphics::par(no.readonly = TRUE)
            on.exit(graphics::par(par_usr))
            if(length(samples) > 1) graphics::par(ask=TRUE)
        }

        ymx <- max(x$reconstructions[[climate]]$likelihood[-1, ], na.rm=TRUE)
        for(tax in x$inputs$taxa.name) {
            if (x$inputs$selectedTaxa[tax, climate] > 0) {
                ymx <- max(ymx, ifelse(is.na(x$modelling$pdfs[[tax]][[climate]]), 0, max(x$modelling$pdfs[[tax]][[climate]]$pdfpol, na.rm=TRUE)))
            }
        }
        ymx <- ymx * 1.05

        COLS <- rep(col, length.out=length(x$inputs$taxa.name))
        names(COLS) <- x$inputs$taxa.name
        LTYS <- 1:length(x$inputs$taxa.name)%/%26+1
        names(LTYS) <- x$inputs$taxa.name

        var_to_plot <- ifelse(optima, 2, 3)
        climate_name <- accClimateVariables(climate)[3]

        continue <- TRUE
        for(s in samples) {
            if(save) {
                if(as.png) {
                    if(length(samples) > 1) {
                        grDevices::png(paste0(strsplit(filename, '.png')[[1]], '_sample_', s, '.png'), width = width, height = height, units='in', res=png.res)
                    } else {
                        grDevices::png(paste0(strsplit(filename, '.png')[[1]], '.png'), width = width, height = height, units='in', res=png.res)
                    }
                } else if(continue) {
                    grDevices::pdf(paste0(strsplit(filename, '.pdf')[[1]], '.pdf'), width=width, height=height)
                    continue <- FALSE
                }
            }


            graphics::layout(matrix(c(1,2), ncol=2), width=c(5,2), height=1)
            graphics::par(mar=c(2,2,0.5,0))
            graphics::par(ps=8, cex=1)

            graphics::plot(0,0, type='n', xaxs='i', yaxs='i', frame=FALSE, axes=FALSE, xlab='', ylab='', main='',
                 xlim=xlim, ylim=c(0, ymx))

            graphics::polygon(x$modelling$xrange[[climate]][c(1, 1:x$parameters$npoints, x$parameters$npoints)],
                    c(0, x$reconstructions[[climate]]$likelihood[s+1, ], 0),
                    col='black', border=NA)

            for(tax in x$inputs$taxa.name) {
                if (x$inputs$selectedTaxa[tax, climate] == 1 & x$modelling$weights[s, tax] > 0) {
                    graphics::polygon(x$modelling$xrange[[climate]][c(1, 1:x$parameters$npoints, x$parameters$npoints)],
                      c(0, x$modelling$pdfs[[tax]][[climate]]$pdfpol, 0),
                      col=crestr::makeTransparent(COLS[tax], alpha=0.3), border=NA)
                }
            }
            for(tax in x$inputs$taxa.name) {
                if (x$inputs$selectedTaxa[tax, climate] == 1 & x$modelling$weights[s, tax] > 0) {
                    graphics::polygon(x$modelling$xrange[[climate]][c(1, 1:x$parameters$npoints, x$parameters$npoints)],
                      c(0, x$modelling$pdfs[[tax]][[climate]]$pdfpol, 0),
                      #lwd=max(0.2, log10(1+10*x$modelling$weights[s, tax])),
                      lwd=max(0.2, 1.5*sqrt(x$modelling$weights[s, tax])),
                      lty=LTYS[tax], col=NA, border=COLS[tax])
                }
            }

            graphics::polygon(x$modelling$xrange[[climate]][c(1, 1:x$parameters$npoints, x$parameters$npoints)],
                    c(0, x$reconstructions[[climate]]$likelihood[s+1, ], 0),
                    col='black', border=NA)
            graphics::segments(x$reconstructions[[climate]]$optima[s, var_to_plot], 0, x$reconstructions[[climate]]$optima[s, var_to_plot], ymx, col='grey70', lwd=1.5, lty=3)

            graphics::rect(xlim[2]*0.99, ymx*0.98, xlim[2]*0.98 - graphics::strwidth(x$inputs$x[s], cex=10/8, lwd=2), ymx*0.99-graphics::strheight(x$inputs$x[s], cex=1.5, lwd=2)*2, col='white', border=NA)
            graphics::text(xlim[2]*0.99, ymx*0.98, x$inputs$x[s], cex=10/8, lwd=2, adj=c(1, 1))

            graphics::box(lwd=0.5)
            graphics::par(mgp=c(3,0.31,0))
            graphics::axis(2,at=c(0,ymx),labels=c("",""),tck=0,lwd=0.5)
            graphics::axis(2,lwd.ticks=0.5,lwd=0,tck=-0.01,cex.axis=6/8)
            graphics::mtext('Density of probability',side=2,line=1.1, cex=7/8)

            graphics::par(mgp=-c(3,0.1,0))
            graphics::mtext(climate_name,side=1,line=0.8, cex=7/8)
            graphics::axis(1,at=range(x$modelling$xrange[[climate]]),labels=c("",""),tck=0,lwd=0.5,pos=0)
            graphics::axis(1, lwd.ticks=0.5,lwd=0,pos=0,tck=-0.01,cex.axis=6/8)

            # taxa selected (ranked by weight), then taxa unselected, then taxa unavailable
            taxa <- rownames(x$inputs$selectedTaxa)[x$inputs$selectedTaxa[, climate] == 1 ]
            weights <- as.numeric(x$modelling$weights[s, taxa])
            if(only.present) {
                taxa <- taxa[weights > 0]
                weights <- weights[weights > 0]
            }
            ordered_tax <- taxa[order(weights, decreasing=TRUE)]

            if(!only.selected) {
                taxa <- rownames(x$inputs$selectedTaxa[x$inputs$selectedTaxa[, climate] == 0, ])
                weights <- as.numeric(x$modelling$weights[s, taxa])
                if(only.present) {
                    taxa <- taxa[weights > 0]
                    weights <- weights[weights > 0]
                }
                ordered_tax <- c(ordered_tax, taxa[order(weights, decreasing=FALSE)])
            }

            graphics::par(mar=c(2,0.1,0.5,0.1))
            graphics::par(ps=8, cex=1)
            graphics::plot( 0,0, type='n', xaxs='i', yaxs='i', frame=FALSE, axes=FALSE, xlab='', ylab='', main='',
                            xlim=c(0,1), ylim=c(1+length(ordered_tax), 1))

            if(length(ordered_tax) > 0) {
                for(tax in 1:length(ordered_tax)) {
                    graphics::segments(0, tax+0.5, 0.20, tax+0.5,
                    #lwd=max(0.2, log10(1+10*x$modelling$weights[s, ordered_tax[tax]])),
                    lwd=max(0.2, 1.5*sqrt(x$modelling$weights[s, ordered_tax[tax]])),
                    lty=LTYS[ordered_tax[tax]], col=ifelse(x$inputs$selectedTaxa[ordered_tax[tax], climate] > 0, COLS[ordered_tax[tax]], 'grey70'))
                    graphics::text(0.23, tax+0.5, paste(ordered_tax[tax], ' (',round(x$modelling$weights[s, ordered_tax[tax]],2),')',sep=''), adj=c(0, 0.45), cex = 7/8, col=ifelse(x$inputs$selectedTaxa[ordered_tax[tax], climate] > 0, 'black', 'grey70'))
                }
            }
            if(save & as.png) grDevices::dev.off()

        }

        if(save) {
            if(!as.png) grDevices::dev.off()
        }
    } else {
        cat('This function only works with a crestObj.\n\n')
    }
    invisible()
}