R/plot-funs.R

In FRBCesab/popbayes: Bayesian Model to Estimate Population Trends from Counts Series

#' Plot original vs converted counts
#' 
#' @description
#' This function plots a panel of two graphics for one count series 
#' (previously generated by [format_data()]):
#' * on the left side, scatter plot overlapping original (black points) and 
#'   converted counts (grey points);
#' * on the right side, scatter plot of converted counts with boundaries of the
#'   95% confident interval.
#' 
#' @param series a `character` string. The count series names (can be 
#'   retrieved by running [list_series()]).
#'   
#' @param title a `logical`. If `TRUE` (default) a title (series name) is 
#'   added.
#'   
#' @param path a `character` string. The directory in which count series have 
#'   been saved by the function [format_data()].
#'   
#' @param path_fig a `character` string. The directory where to save the plot 
#'   (if `save = TRUE`). This directory must exist and can be an absolute or a 
#'   relative path.
#'   
#' @param save a `logical`. If `TRUE` (default is `FALSE`) the plot is saved in 
#'   `path_fig`.
#'
#' @return No return value.
#' 
#' @export
#'
#' @examples
#' ## Load Garamba raw dataset ----
#' file_path <- system.file("extdata", "garamba_survey.csv", 
#'                          package = "popbayes")
#'                          
#' garamba <- read.csv(file = file_path)
#' 
#' ## Create temporary folder ----
#' temp_path <- tempdir()
#' 
#' ## Format dataset ----
#' garamba_formatted <- popbayes::format_data(
#'   data              = garamba, 
#'   path              = temp_path,
#'   field_method      = "field_method",
#'   pref_field_method = "pref_field_method",
#'   conversion_A2G    = "conversion_A2G",
#'   rmax              = "rmax")
#' 
#' ## Get series names ----
#' popbayes::list_series(path = temp_path)
#' 
#' ## Plot for Alcelaphus buselaphus at Garamba ----
#' popbayes::plot_series("garamba__alcelaphus_buselaphus", path = temp_path)

plot_series <- function(series, title = TRUE, path = ".", path_fig = ".",
                        save = FALSE) {
  
  
  if (!dir.exists(path)) {
    stop("The directory '", path, "' does not exist.")
  }
  
  if (missing(series)) {
    stop("Argument 'series' is required.")
  }
  
  if (!is.character(series) || length(series) != 1) {
    stop("Argument 'series' must be character string (one series name).")
  }
  
  if (!dir.exists(file.path(path, series))) {
    stop("The directory '", file.path(path, series), "' does not exist.")
  }
  
  if (!file.exists(file.path(path, series, paste0(series, "_data.RData")))) {
    stop("The file '", file.path(path, series, 
                                 paste0(series, "_data.RData")), 
         "' does not exist. Did you forget to run format_data()?")
  }
  
  
  ## Load formatted data ----
  
  data <- get(load(file.path(path, series, paste0(series, "_data.RData"))))
  
  data_copy <- data
  data <- data[[series]]
  
  
  ## Get plot extent ----
  
  counts_range <- c(0, max(c(data$"data_converted"$"upper_ci_conv", 
                             data$"data_original"$"count_orig")))
  dates_range <- range(data$"date")
  
  
  ## Export as PNG (if required) ----
  
  if (save) {
    
    if (!dir.exists(path_fig)) {
      stop("The directory '", path_fig, "' does not exist.")
    }
    
    filename <- paste0(names(data_copy[series]), "_data.png")
    grDevices::png(filename = file.path(path_fig, filename), width = 12, 
                   height = 6, pointsize = 14, res = 600, units = "in")
  }
  
  
  ## Graphical parameters ----
  
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  
  par(family = "serif", mgp = c(1.5, 0.1, 0), mar = c(1.5, 3.1, 1.9, 1.1),
      cex.axis = 0.75)
  par(mfrow = c(1, 2))
  
  
  ###
  ### Left plot - Original vs. Converted
  ###
  
  ## Empty plot ----
  
  plot(0, type = "n", xlim = dates_range, ylim = counts_range, axes = FALSE,
       xlab = "", ylab = "Original (black) and converted (grey) counts",
       font.lab = 1)
  
  
  ## Background ----
  
  grid(lwd = 0.25, lty = 1)
  axis(1, at = axTicks(1), lwd = 0)
  axis(2, at = axTicks(2), lwd = 0)
  box(lwd = 1.0, col = "white")
  box(lwd = 0.5)
  
  
  ## Deviation between original and converted data ----
  
  for (i in seq_len(nrow(data$"data_original"))) {
    
    lines(x = rep(data$"data_original"[i, "date"], 2),
          y = c(data$"data_original"[i, "count_orig"], 
                data$"data_converted"[i, "count_conv"]),
          lty = 4, lwd = 0.75)
  }
  
  
  ## Add data ----
  
  points(x = data$"data_original"$"date", 
         y = data$"data_original"$"count_orig", pch = 19, cex = 1.2)
  
  points(x = data$"data_converted"$"date", 
         y = data$"data_converted"$"count_conv", pch = 19, cex = 1.2)
  
  points(x = data$"data_converted"$"date", 
         y = data$"data_converted"$"count_conv", pch = 19, col = "#aaaaaa")
  
  
  ###
  ### Right plot - Converted with 95% CI
  ###
  
  
  ## Empty plot ----
  
  plot(0, type = "n", xlim = dates_range, ylim = counts_range, axes = FALSE,
       xlab = "", ylab = "Converted counts with 95% CI",
       font.lab = 1)
  
  
  ## Background ----
  
  grid(lwd = 0.25, lty = 1)
  axis(1, at = axTicks(1), lwd = 0)
  axis(2, at = axTicks(2), lwd = 0)
  box(lwd = 1.0, col = "white")
  box(lwd = 0.5)
  
  
  ## Add 95% CI ----
  
  for (i in seq_len(nrow(data$"data_converted"))) {
    
    lines(x = rep(data$"data_converted"[i, "date"], 2),
          y = c(data$"data_converted"[i, "lower_ci_conv"], 
                data$"data_converted"[i, "upper_ci_conv"]))
    
    
    lines(x = c(data$"data_converted"[i, "date"] - 
                  data$"data_converted"[i, "date"] * 0.0002,
                data$"data_converted"[i, "date"] + 
                  data$"data_converted"[i, "date"] * 0.0002),
          y = rep(data$"data_converted"[i, "lower_ci_conv"], 2))
    
    lines(x = c(data$"data_converted"[i, "date"] - 
                  data$"data_converted"[i, "date"] * 0.0002,
                data$"data_converted"[i, "date"] + 
                  data$"data_converted"[i, "date"] * 0.0002),
          y = rep(data$"data_converted"[i, "upper_ci_conv"], 2))
  }
  
  
  ## Add data ----
  
  points(x = data$"data_converted"$"date", 
         y = data$"data_converted"$"count_conv", pch = 19, cex = 1.2)
  
  points(x = data$"data_converted"$"date", 
         y = data$"data_converted"$"count_conv", pch = 19, col = "#aaaaaa")
  
  
  ## Add title ----
  
  if (title)
    mtext(paste(data$"location", "-", data$"species"), side = 3, line = 0.5, 
          adj = 1, at = par()$usr[2], cex = 1.2, font = 4)
  
  
  ## Restore user graphical parameters ----
  
  if (save)
    dev.off()
  
  invisible(NULL) 
}



#' Plot estimated population trend
#' 
#' @description
#' This function plots a panel of two graphics for one BUGS model
#' (previously generated by [fit_trend()]):
#' * on the left side, the population trend estimated by the Bayesian model 
#'   (blue line) with the 95% CI (gray envelop). Dots (with intervals) 
#'   represent converted counts passed to the model (with the 95% CI);
#' * on the right side, a bar plot of estimated relative growth rates (r) by 
#'   date. Dark bars are real estimated r.
#' 
#' @param series a `character` string. The count series name (can be 
#'   retrieved by running [list_series()]).
#'   
#' @param title a `logical`. If `TRUE` (default) a title (series name) is 
#'   added.
#'   
#' @param path a `character` string. The directory in which count series (and 
#'   BUGS outputs) have been saved by the function [format_data()] (and by
#'   [fit_trend()]).
#'   
#' @param path_fig a `character` string. The directory where to save the plot 
#'   (if `save = TRUE`). This directory must exist and can be an absolute or a 
#'   relative path.
#'   
#' @param save a `logical`. If `TRUE` (default is `FALSE`) the plot is saved in 
#'   `path_fig`.
#'
#' @return No return value.
#' 
#' @export
#'
#' @examples
#' ## Load Garamba raw dataset ----
#' file_path <- system.file("extdata", "garamba_survey.csv", 
#'                          package = "popbayes")
#'                          
#' garamba <- read.csv(file = file_path)
#' 
#' ## Create temporary folder ----
#' temp_path <- tempdir()
#' 
#' ## Format dataset ----
#' garamba_formatted <- popbayes::format_data(
#'   data              = garamba, 
#'   path              = temp_path,
#'   field_method      = "field_method",
#'   pref_field_method = "pref_field_method",
#'   conversion_A2G    = "conversion_A2G",
#'   rmax              = "rmax")
#' 
#' ## Select one serie ----
#' a_buselaphus <- popbayes::filter_series(garamba_formatted, 
#'                                         location = "Garamba",
#'                                         species  = "Alcelaphus buselaphus")
#' \donttest{
#' ## Fit population trends (requires JAGS) ----
#' a_buselaphus_mod <- popbayes::fit_trend(a_buselaphus, path = temp_path)
#' 
#' ## Plot estimated population trend ----
#' popbayes::plot_trend(series = "garamba__alcelaphus_buselaphus", 
#'                      path   = temp_path)
#' }

plot_trend <- function(series, title = TRUE, path = ".", path_fig = ".", 
                       save = FALSE) {
  
  if (!dir.exists(path)) {
    stop("The directory '", path, "' does not exist.")
  }
  
  if (missing(series)) {
    stop("Argument 'series' is required.")
  }
  
  if (!is.character(series) || length(series) != 1) {
    stop("Argument 'series' must be a character string (one series name).")
  }
  
  if (!dir.exists(file.path(path, series))) {
    stop("The directory '", file.path(path, series), "' does not exist.")
  }
  
  if (!file.exists(file.path(path, series, paste0(series, "_data.RData")))) {
    stop("The file '", file.path(path, series, 
                                 paste0(series, "_data.RData")), 
         "' does not exist. Did you forget to run format_data()?")
  }
  
  if (!file.exists(file.path(path, series, paste0(series, "_bugs.RData")))) {
    stop("The file '", file.path(path, series, 
                                 paste0(series, "_bugs.RData")), 
         "' does not exist. Did you forgot to run fit_trend()?")
  }
  
  if (!is.logical(title)) {
    stop("Argument 'title' must be TRUE or FALSE.")
  }
  
  if (!is.logical(save)) {
    stop("Argument 'save' must be TRUE or FALSE.")
  }
  
  
  ## Load formatted data ----
  
  data <- get(load(file.path(path, series, paste0(series, "_data.RData"))))
  
  data_copy <- data
  
  data <- data[[1]]$"data_converted"
  
  
  ## Load BUGS outputs ----
  
  outputs <- get(load(file.path(path, series, paste0(series, "_bugs.RData"))))
  outputs <- outputs[[1]]$"BUGSoutput"$"summary"
  
  
  ## Select estimated N ----
  
  outputs_n <- outputs[grep("^N\\[", rownames(outputs)), c("mean", "2.5%", 
                                                           "97.5%")]
  rownames(outputs_n) <- NULL
  colnames(outputs_n) <- c("count_jags", "lower_ci_jags", "upper_ci_jags")
  
  
  ## Create table N for graphic ----
  
  data_n_for_graph <- data.frame("date"     = data$"date", 
                                 "count"    = data$"count_conv",
                                 "lower_ci" = data$"lower_ci_conv", 
                                 "upper_ci" = data$"upper_ci_conv")
  
  data_n_for_graph <- cbind(data_n_for_graph, outputs_n)
  
  
  ## Select estimated r ----
  
  outputs_r <- outputs[grep("^r\\[", rownames(outputs)), "mean"]
  outputs_r <- c(outputs_r, NA)
  outputs_r <- data.frame("date" = data$"date", "r_jags" = outputs_r)
  rownames(outputs_r) <- NULL
  
  
  ## Create table r for graphic ----
  
  dates_seq <- data.frame(
    "date" = seq(min(outputs_r$"date"), max(outputs_r$"date")))
  
  outputs_r <- merge(outputs_r, dates_seq, by = "date", all = TRUE)
  
  outputs_r$"color" <- ifelse(is.na(outputs_r$"r_jags"), "#888888", "black")
  outputs_r$"lwd"   <- ifelse(is.na(outputs_r$"r_jags"), 1, 1.2)
  
  for (i in seq_len(nrow(outputs_r))) {
    if (is.na(outputs_r[i, "r_jags"])) {
      outputs_r[i, "r_jags"] <- outputs_r[i - 1, "r_jags"]
    }
  }
  
  data_r_for_graph <- outputs_r
  
  
  ## Get plot extent ----
  
  dates_range  <- range(data_n_for_graph$"date")
  
  counts_range <- c(0, max(c(data_n_for_graph$"upper_ci", 
                             data_n_for_graph$"upper_ci_jags")))
  
  growth_range <- c(min(data_r_for_graph$"r_jags"),
                    max(data_r_for_graph$"r_jags"))
  
  growth_range[1] <- ifelse(growth_range[1] > 0, 0, growth_range[1])
  growth_range[2] <- ifelse(growth_range[2] < 0, 0, growth_range[2])
  growth_range[2] <- growth_range[2] + 
    (growth_range[2] - growth_range[1]) * 0.05
  
  
  ## Export as PNG (if required) ----
  
  if (save) {
    
    if (!dir.exists(path_fig)) {
      stop("The directory '", path_fig, "' does not exist.")
    }
    
    filename <- paste0(names(data_copy), "_trend.png")
    grDevices::png(filename = file.path(path_fig, filename), width = 12, 
                   height = 6, pointsize = 14, res = 600, units = "in")
  }
  
  
  ## Graphical parameters ----
  
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  
  par(family = "serif", mgp = c(1.5, 0.1, 0), mar = c(1.5, 3.1, 1.9, 1.1),
      cex.axis = 0.75)
  par(mfrow = c(1, 2))
  
  
  
  ###
  ### Left plot - Population trend
  ###
  
  
  
  ## Empty plot ----
  
  plot(0, type = "n", xlim = dates_range, ylim = counts_range, axes = FALSE,
       xlab = "", ylab = "Population size", font.lab = 1)
  
  
  ## Background ----
  
  grid(lwd = 0.25, lty = 1)
  axis(1, at = axTicks(1), lwd = 0)
  axis(2, at = axTicks(2), lwd = 0)
  box(lwd = 1.0, col = "white")
  box(lwd = 0.5)
  
  
  ## Add JAGS CI polygon ----
  
  polygon(x = c(data_n_for_graph$"date", 
                data_n_for_graph$"date"[rev(seq_len(nrow(data_n_for_graph)))]),
          y = c(data_n_for_graph[ , "lower_ci_jags"],
                data_n_for_graph[rev(seq_len(nrow(data_n_for_graph))), 
                                 "upper_ci_jags"]),
          col = "#bbbbbb", border = "#bbbbbb", lwd = 1)
  
  
  ## Add JAGS mean prediction ----
  
  lines(x = data_n_for_graph$"date", y = data_n_for_graph$"count_jags", 
        col = "steelblue", lwd = 2)
  
  
  ## Add Observed (converted) 95% CI ----
  
  for (i in seq_len(nrow(data_n_for_graph))) {
    
    lines(x = rep(data_n_for_graph[i, "date"], 2),
          y = c(data_n_for_graph[i, "lower_ci"], 
                data_n_for_graph[i, "upper_ci"]))
    
    
    lines(x = c(data_n_for_graph[i, "date"] - 
                  data_n_for_graph[i, "date"] * 0.0002,
                data_n_for_graph[i, "date"] + 
                  data_n_for_graph[i, "date"] * 0.0002),
          y = rep(data_n_for_graph[i, "lower_ci"], 2))
    
    lines(x = c(data_n_for_graph[i, "date"] - 
                  data_n_for_graph[i, "date"] * 0.0002,
                data_n_for_graph[i, "date"] + 
                  data_n_for_graph[i, "date"] * 0.0002),
          y = rep(data_n_for_graph[i, "upper_ci"], 2))
  }
  
  
  ## Add Observed (converted) counts ----
  
  points(x = data_n_for_graph$"date", y = data_n_for_graph$"count", pch = 19, 
         cex = 1.2)
  
  
  
  ###
  ### Right plot - Estimated r barplot
  ###
  
  
  
  ## Empty plot ----
  
  plot(0, type = "n", xlim = dates_range, ylim = growth_range, axes = FALSE,
       xlab = "", ylab = "Relative growth rate (r)",
       font.lab = 1)
  
  
  ## Background ----
  
  grid(lwd = 0.25, lty = 1)
  axis(1, at = axTicks(1), lwd = 0)
  axis(2, at = axTicks(2), lwd = 0)
  
  
  ## Add data ----
  
  points(x = data_r_for_graph$"date", y = data_r_for_graph$"r", type = "h", 
         col = data_r_for_graph$"color", lwd = data_r_for_graph$"lwd")
  
  
  ## Add average value (top-left text) ----
  
  r_mean <- outputs[grep("^meanr$", rownames(outputs)), c("mean", "2.5%", 
                                                          "97.5%")]
  r_mean <- round(r_mean, 3)
  
  text(x = par()$usr[1], y = growth_range[2],
       labels = bquote(bar(r) == .(r_mean["mean"]) ~ .(paste0("[", 
                                                              r_mean["2.5%"], 
                                                              ", ", 
                                                              r_mean["97.5%"], 
                                                              "]"))), pos = 4)
  
  
  ## Add y = 0 line ----
  
  abline(h = 0, lwd = 1.2, col = "white")
  abline(h = 0, lwd = 0.5)
  box(lwd = 1.0, col = "white")
  box(lwd = 0.5)
  
  
  ## Add title ----
  
  if (title)
    mtext(paste(data_copy[[1]]$"location", "-", data_copy[[1]]$"species"), 
          side = 3, line = 0.5, adj = 1, at = par()$usr[2], cex = 1.2, 
          font = 4)
  
  
  ## Restore user graphical parameters ----
  
  if (save)
    dev.off()
  
  invisible(NULL) 
}