R/plotting_functions.R

In EriqLaplus/OZAB: Implemtation of Ordinal Zero-Augmented Beta Models in Stan

#' Plot Cumulative Empirical Proportions of Species by Covariate Bins
#'
#' First found in Irvine, Rodhouse, Keren 2016
#'
#' @param data Dataframe in long format
#' @param covariate Name of Categorical Covariate to Plot By
#' @param cover_class_col Name of Cover Class column if it differs from Default 'Cover Class'
#'
#' @return ggplot2 object
#' @export
#'
#' @examples
#'   cumulative_plot(sagebrush, Fire)
#' @importFrom rlang .data
cumulative_plot <- function(data, covariate, cover_class_col = .data$`Cover Class`){

  data %>%
    ggplot2::ggplot(ggplot2::aes(x = {{ covariate }}, fill = {{ cover_class_col }})) +
      ggplot2::geom_bar(position = ggplot2::position_fill()) +
      ggplot2::ylab("Empirical Cumulative Proportions")
}

#' Plot Prior Sensitivity
#'
#' @param OZABfit resulting stan object from an OZAB modle
#' @param pattern Regex pattern to match covariate names on
#'
#' @return ggplot object
#' @export
#'
#' @examples
plot_prior_sensitivity <- function(OZABfit, pattern = NULL){

  # Extract Selected Posteriors
  if(!is.null(pattern)){
    param_names <- names(OZABfit@prior)[grepl(pattern, names(OZABfit@prior))]
  } else {
    param_names <- names(OZABfit@prior)
  }

  if(length(param_names) == 0) {
    stop(paste('No parameters found matching pattern:', pattern))
  }

  posteriors <-
    rstan::extract(OZABfit, param_names) %>%
    tibble::as_tibble() %>%
    tidyr::pivot_longer(everything(), names_to = 'parameter', values_to = 'value')

  limits <-
    posteriors %>%
    dplyr::summarise(min = min(value), max = max(value))

  x <- seq(from = limits[['min']],
           to = limits[['max']],
           length.out = nrow(posteriors))

  priors <- NULL

  # Compute Prior Densities
  for(param_name in param_names) {
    priors <- rbind(priors, dplyr::tibble(parameter = param_name, x = x, y = dnorm(x, mean = OZABfit@prior[[param_name]][1], sd = OZABfit@prior[[param_name]][2])))
  }

  # Plot Posteriors
  posteriors %>%
    ggplot2::ggplot(ggplot2::aes(x = value)) +
      ggplot2::geom_histogram(ggplot2::aes(y = stat(density))) +
      ggplot2::facet_wrap(~parameter, scales = "free_x") +
      ggplot2::geom_line(ggplot2::aes(x = x, y = y), data=filter(priors, parameter == parameter))
}