R/eems.plots.R

In reems: Estimating Effective Migration Surfaces from Single Nucleotide Polymorphism Data

#' A function to plot effective migration and diversity surfaces from EEMS output
#'
#' Given a vector of EEMS output directories, this function generates several figures
#' to visualize EEMS results. It is a good idea to examine all these figures, which is
#' why they are generated by default.
#' \itemize{
#'   \item \code{plotpath-mrates01}: effective migration surface. This contour plot visualizes
#'   the estimated effective migration rates \code{m}, on the log10 scale after mean centering.
#'   \item \code{plotpath-mrates02}: posterior probabilities \code{P(m > 0 | diffs)} and
#'   \code{P(m < 0 | diffs)} for each location in the habitat. Since migration rates are
#'   visualized on the log10 scale after mean centering, 0 corresponds to the overall mean
#'   migration rate. This contour plot emphasizes regions with effective migration that is
#'   significantly higher/lower than the overall average.
#'   \item \code{plotpath-qrates01}: effective diversity surface. This contour plot visualizes
#'   the estimated effective diversity rates \code{q}, on the log10 scale after mean centering.
#'   \item \code{plotpath-qrates02}: posterior probabilities \code{P(q > 0 | diffs)} and
#'   \code{P(q < 0 | diffs)}. Similar to \code{plotpath-mrates02} but applied to the effective
#'   diversity rates.
#'   \item \code{plotpath-rdist01}: scatter plot of the observed vs the fitted between-deme
#'   component of genetic dissimilarity, where one point represents a pair of sampled demes.
#'   \item \code{plotpath-rdist02}: scatter plot of the observed vs the fitted within-deme
#'   component of genetic dissimilarity, where one point represents a sampled deme.
#'   \item \code{plotpath-rdist03}: scatter plot of observed genetic dissimilarities between
#'   demes vs observed geographic distances between demes.
#'   \item \code{plotpath-pilogl01}: posterior probability trace
#' }
#' The \code{mrates} and \code{qrates} figures visualize (properties of) the effective migration
#' and diversity rates across the habitat. The other figures can help to check that the MCMC
#' sampler has converged (the trace plot \code{pilogl}) and that the EEMS model fits the data
#' well (the scatter plots of genetic dissimilarities \code{rdist}).
#'
#' The function \code{eems.plots} will work given the results from a single EEMS run (one
#' directory in \code{mcmcpath}) but it is better to run EEMS several times, randomly initializing
#' the MCMC chain each time. In other words, simulate several realizations of the Markov chain
#' and let each realization start from a different state in the parameter space (by using a
#' different random seed).
#'
#' Detail about the within-deme and between-deme components of genetic dissimilarity: Let
#' \code{D(a,b)} be the dissimilarity between one individual from deme \code{a} and another
#' individual from deme \code{b}. Then the within-deme component for \code{a} and \code{b}
#' is simply \code{D(a,a)} and \code{D(b, b)}, respectively. The between-deme component is
#' \code{D(a,b) - [D(a,a) + D(b,b)] / 2} and it represents dissimilarity that is due to the spatial
#' structure of the population and is not a consequence of the local diversity in the two demes.
#' @param mcmcpath A vector of EEMS output directories, for the same dataset. Warning: There is
#' minimal checking that the given  directories are for the same dataset.
#' @param plotpath The full path and the file name for the graphics to be generated.
#' @param longlat A logical value indicating whether the coordinates are given as pairs
#' (longitude, latitude) or (latitude, longitude).
#' @param plot.width The width of the graphics region for the two rate contour plots, in inches.
#' The default value is 10.
#' @param plot.height The height of the graphics region, in inches. The default value is 10.
#' @param out.png A logical value which, if set,  forces output graphics to be generated as PNGs
#' (if `TRUE`) or PDFs (if `FALSE`). If left unset, the format depends on the nature of the plot. 
#' @param res Resolution, in dots per inch; used only for PNG images.
#' The default is 600.
#' @param xpd A logical value indicating whether to clip plotting to the figure region
#' (\code{xpd} = TRUE, which is the default) or clip plotting to the plot region
#' (\code{xpd} = FALSE).
#' @param add.grid A logical value indicating whether to add the population grid or not.
#' @param col.grid The color of the population grid. Defaults to \code{gray80}.
#' @param lwd.grid The line width of the population grid. Defaults to 1.
#' @param add.outline A logical value indicating whether to add the habitat outline or not.
#' @param col.outline The color of the habitat outline. Defaults to \code{white}.
#' @param lwd.outline The line width of the habitat outline. Defaults to 2.
#' @param add.demes A logical value indicating whether to add the observed demes or not.
#' @param col.demes The color of the demes. Defaults to \code{black}.
#' @param pch.demes The symbol, specified as an integer, or the character to be used for
#' plotting the demes. Defaults to 19.
#' @param min.cex.demes The minimum size of the deme symbol/character.
#' @param max.cex.demes The maximum size of the deme symbol/character. Defaults to 1 and 3,
#' respectively. If \code{max.cex.demes} > \code{min.cex.demes}, then demes with more samples
#' also have bigger size: the deme with the fewest samples has size \code{min.cex.demes} and
#' the deme with the most samples has size \code{max.cex.demes}.
#' @param projection.in The input cartographic projection, specified as a PROJ.4 string.
#' @param projection.out The output cartographic projection, specified as a PROJ.4 string.
#' @param add.map A logical value indicating whether to add a high-resolution geographic map.
#' Requires the \code{rworldmap} and \code{rworldxtra} packages. It also requires that
#' \code{projection.in} is specified.
#' @param col.map The color of the geographic map. Default is \code{gray60}.
#' @param lwd.map The line width of the geographic map. Defaults to 2.
#' @param eems.colors The EEMS color scheme as a vector of colors, ordered from low to high.
#' Defaults to a DarkOrange to Blue divergent palette with six orange shades, white in the middle,
#' six blue shades. Acknowledgement: The default color scheme is adapted from the \code{dichromat}
#' package.
#' @param m.colscale A fixed range for log10-transformed migration rates. If the estimated rates
#' fall outside the specified range, then the color scale is ignored. By default, no range is
#' specified for either type of rates.
#' @param q.colscale A fixed range for log10-transformed diversity rates.
#' @param add.colbar A logical value indicating whether to add the color bar (the key that shows
#' how colors map to rates) to the right of the plot. Defaults to TRUE.
#' @param remove.singletons Remove demes with a single observation from the diagnostic scatter
#' plots. Defaults to TRUE.
#' @param add.abline Add the line \code{y = x} to the diagnostic scatter plots of observed vs
#' fitted genetic dissimilarities.
#' @param add.r.squared Add the R squared coefficient to the diagnostic scatter plots of observed
#' vs fitted genetic dissimilarities.
#' @param add.title A logical value indicating whether to add the main title in the contour plots.
#' Defaults to TRUE.
#' @param prob.levels A vector of probabilities for plotting the posterior probability contours of
#' \code{P(m > 0 | diffs)} and \code{P(m < 0 | diffs)}. Defaults to \code{c(0.9, 0.95)}.
#' @param m.plot.xy Statements which add graphical elements (e.g. points) on top of the migration
#' sufrace.
#' @param q.plot.xy Statements which add graphical elements (e.g. points) on top of the diversity
#' surface.
#' @param xy.coords Additional coordinates at which to estimate the migration and diversity rates.
#' @returns None 
#' @references Light A and Bartlein PJ (2004). The End of the Rainbow? Color Schemes for Improved
#' Data Graphics. EOS Transactions of the American Geophysical Union, 85(40), 385.
#' @examples
#' # Use the provided example or supply the path to your own EEMS run.
#' extdata_path <- system.file("extdata", package = "reems")
#' eems_results <- file.path(extdata_path, "EEMS-example")
#' # Create a temporary output directory for the sake of the example
#' outdir <- file.path(tempdir(), "plot_out")
#' dir.create(outdir, showWarnings = FALSE) 
#' name_figures <- file.path(outdir, "EEMS-example")
#'
#' # Produce the set of EEMS figures, with default values for all optional parameters.
#' eems.plots(
#'   mcmcpath = eems_results,
#'   plotpath = paste0(name_figures, "-default"),
#'   longlat = TRUE,
#'   out.png = FALSE
#' )
#' 
#' # Delete the temporary output directory to tidy up. 
#' unlink(outdir, recursive = TRUE, force = TRUE)
#' @seealso \code{\link{eems.voronoi.samples}, \link{eems.posterior.draws}, \link{eems.population.grid}}
#' @export

eems.plots <- function(mcmcpath,
                       plotpath,
                       longlat,
                       plot.width = 10,
                       plot.height = 10,
                       out.png = NULL,
                       res = 600,
                       xpd = TRUE,
                       add.grid = FALSE,
                       col.grid = "gray80",
                       lwd.grid = 1,
                       add.demes = FALSE,
                       col.demes = "black",
                       pch.demes = 19,
                       min.cex.demes = 1,
                       max.cex.demes = 3,
                       add.outline = FALSE,
                       col.outline = "gray90",
                       lwd.outline = 2,
                       projection.in = NULL,
                       projection.out = NULL,
                       add.map = FALSE,
                       col.map = "gray60",
                       lwd.map = 2,
                       eems.colors = NULL,
                       prob.levels = c(0.9, 0.95),
                       add.colbar = TRUE,
                       m.colscale = NULL,
                       q.colscale = NULL,
                       remove.singletons = TRUE,
                       add.abline = FALSE,
                       add.r.squared = FALSE,
                       add.title = TRUE,
                       m.plot.xy = NULL,
                       q.plot.xy = NULL,
                       xy.coords = NULL) {
  plot.params <- list(
    eems.colors = eems.colors, m.colscale = m.colscale, q.colscale = q.colscale,
    add.map = add.map, add.grid = add.grid, add.outline = add.outline, add.demes = add.demes,
    col.map = col.map, col.grid = col.grid, col.outline = col.outline, col.demes = col.demes,
    lwd.map = lwd.map, lwd.grid = lwd.grid, lwd.outline = lwd.outline, pch.demes = pch.demes,
    min.cex.demes = min.cex.demes, proj.in = projection.in, add.colbar = add.colbar,
    max.cex.demes = max.cex.demes, proj.out = projection.out, add.title = add.title,
    prob.levels = prob.levels
  )
  plot.params <- check.plot.params(plot.params)

  # A vector of EEMS output directories, for the same dataset.
  # Assume that if eemsrun.txt exists, then all EEMS output files exist.
  mcmcpath <- mcmcpath[file.exists(file.path(mcmcpath, "eemsrun.txt"))]
  if (!length(mcmcpath)) {
    stop("Please provide at least one existing EEMS output directory, mcmcpath")
  }

  dimns <- read.dimns(mcmcpath[1], longlat, coord = xy.coords)

  save.params <- list(height = plot.height, width = plot.width, res = res, out.png = out.png)

  message("Processing the following EEMS output directories :")
  message(mcmcpath)

  
  # Plot filled contour of estimated effective migration rates
  mrates.raster <- plot_with_graphics_device(
      paste0(plotpath, "-mrates"),
      save.params,
      average.eems.contours(
          mcmcpath, dimns, longlat, plot.params,
          is.mrates = TRUE, plot.xy = m.plot.xy
      ),
      def.out = TRUE,
      par.args = list(las = 1, font.main = 1, xpd = xpd)
  )
  xym.values <- mrates.raster$xyz.values
  colnames(xym.values) <- c("x", "y", "m")

   # Plot filled contour of estimated effective diversity rates
  qrates.raster <- plot_with_graphics_device(
      paste0(plotpath, "-qrates"),
      save.params, 
      average.eems.contours(
          mcmcpath, dimns, longlat, plot.params,
          is.mrates = FALSE, plot.xy = q.plot.xy
      ),
      def.out = TRUE,
      par.args = list(las = 1, font.main = 1, xpd = xpd)
  )
  xyq.values <- qrates.raster$xyz.values
  colnames(xyq.values) <- c("x", "y", "q")

  if (!add.colbar) {
    if (out.png) {
      save.params$height <- 6
      save.params$width <- 1.5
    } else {
      save.params$height <- 12
      save.params$width <- 3
    }
    plot_with_graphics_device(
        paste0(plotpath, "-mkey"),
        save.params,
        myfilled.legend(
            levels = mrates.raster$eems.levels,
            col = mrates.raster$eems.colors,
            key.axes = axis(4, tick = FALSE, hadj = 1, line = 4, cex.axis = 2),
            key.title = mtext(expression(paste(log, "(", italic(m), ")", sep = "")),
                              side = 3, cex = 2.5, line = 1.5, font = 1
                              )
        ),
        def.out = TRUE,
        par.args = list(las = 1, font.main = 1, xpd = xpd, mar = c(0, 1, 5, 8))
    )
    plot_with_graphics_device(
        paste0(plotpath, "-qkey"),
        save.params,
        myfilled.legend(
            levels = qrates.raster$eems.levels,
            col = qrates.raster$eems.colors,
            key.axes = axis(4, tick = FALSE, hadj = 1, line = 6, cex.axis = 2),
            key.title = mtext(expression(paste(log, "(", italic(q), ")", sep = "")),
                              side = 3, cex = 2.5, line = 1.5, font = 1
                              )
        ),
        def.out = TRUE,
        par.args = list(las = 1, font.main = 1, xpd = xpd, mar = c(0, 1, 5, 8))
    )
  }

  save.params$height <- 6
  save.params$width <- 6.5
  # Plot trace plot of posterior probability to check convergence
  plot_with_graphics_device(
      paste0(plotpath, "-pilogl"),
      save.params,
      plot.logposterior(mcmcpath),
      def.out = FALSE,
      par.args = list(las = 0, font.main = 1, mar = c(5, 5, 4, 5) + 0.1, xpd = TRUE)
  )
  

  # Plot scatter plots of observed vs fitted genetic differences
  dist.points <- plot_with_graphics_device(
      paste0(plotpath, "-rdist"),
      save.params,
      dist.scatterplot(
          mcmcpath, longlat, plot.params,
          remove.singletons = remove.singletons,
          add.abline = add.abline, add.r.squared = add.r.squared
      ),
      def.out = FALSE,
      par.args = list(las = 1, font.main = 1, mar = c(5, 5, 4, 2) + 0.1)
  )

  B.component <- dist.points$B.component
  W.component <- dist.points$W.component
  G.component <- dist.points$G.component
  save(B.component, W.component, G.component, xym.values, xyq.values,
    file = paste0(plotpath, "-rdist.RData")
  )
}