R/plot.R

In bvhar: Bayesian Vector Heterogeneous Autoregressive Modeling

#' Density Plot for Minnesota Prior VAR Model
#' 
#' This function draws density plot for coefficient matrices of Minnesota prior VAR model.
#' 
#' @param object A `summary.normaliw` object
#' @param type The type of the plot. Trace plot (`trace`), kernel density plot (`dens`), and interval estimates plot (`area`).
#' @param pars Parameter names to draw.
#' @param regex_pars Regular expression parameter names to draw.
#' @param ... Other options for each [bayesplot::mcmc_trace()], [bayesplot::mcmc_dens()], and [bayesplot::mcmc_areas()].
#' @return A ggplot object
#' @importFrom bayesplot mcmc_trace mcmc_dens mcmc_areas
#' @export
autoplot.summary.normaliw <- function(object, type = c("trace", "dens", "area"), pars = character(), regex_pars = character(), ...) {
  type <- match.arg(type)
  bayes_plt <- switch(
    type,
    "trace" = mcmc_trace(x = object$param, pars = pars, regex_pars = regex_pars, ...),
    "dens" = mcmc_dens(x = object$param, pars = pars, regex_pars = regex_pars, ...),
    "area" = mcmc_areas(x = object$param, pars = pars, regex_pars = regex_pars, ...)
  )
  bayes_plt
}

#' Residual Plot for Minnesota Prior VAR Model
#' 
#' This function draws residual plot for covariance matrix of Minnesota prior VAR model.
#' 
#' @param object A `normaliw` object
#' @param hcol color of horizontal line = 0 (By default, grey)
#' @param hsize size of horizontal line = 0 (By default, 1.5)
#' @param ... additional options for geom_point
#' @return A ggplot object
#' @importFrom ggplot2 ggplot aes geom_point geom_hline facet_grid labs element_text element_blank
#' @importFrom tidyr pivot_longer
#' @export
autoplot.normaliw <- function(object, hcol = "grey", hsize = 1.5, ...) {
  X <- object$residuals |> as.data.frame()
  X[["id"]] <- 1:object$obs
  X <- 
    X |> 
    pivot_longer(-id, names_to = "name", values_to = "value")
  X |> 
    ggplot(aes(x = id, y = value)) +
    geom_hline(yintercept = 0, col = hcol, size = hsize) +
    geom_point(...) +
    facet_grid(
      name ~ .,
      scales = "free_y"
    ) +
    labs(
      x = element_blank(),
      y = element_blank()
    )
}

#' Make Data Form to Make Forecasting Plot
#' 
#' @param object A `predbvhar` object usually generated by [stats::predict()] function in this package
#' @return `gather_predbvhar` produces a `tibble` with columns named
#' * `forecast` - If the row is observed one or forecasted one
#' * `id` - Index of the row
#' * `variable` - Name of the variable
#' * `value_forecast` - Point forecasts
#' * `value_lower_joint` - Lower CI
#' * `value_upper_joint` - Upper CI
#' * `Model` - Fitting Model
#' @importFrom dplyr bind_rows mutate n left_join
#' @importFrom purrr reduce
#' @importFrom tidyr pivot_longer
#' @noRd
gather_predbvhar <- function(object) {
  Y <- 
    object$y |> 
    as.data.frame() |> 
    mutate(forecast = FALSE)
  PRED <- 
    object$forecast |> 
    as.data.frame() |> 
    mutate(forecast = TRUE)
  lapply(
    c("forecast", "lower_joint", "upper_joint"),
    function(comp) {
      PRED <- 
        object[[comp]] |> 
        as.data.frame() |> 
        mutate(forecast = TRUE)
      Y |> 
        bind_rows(PRED) |> 
        mutate(id = 1:n()) |> 
        pivot_longer(-c(id, forecast), names_to = "variable", values_to = paste("value", comp, sep = "_"))
    }
  ) |> 
    reduce(left_join, by = c("id", "variable", "forecast")) |> 
    mutate(Model = object$process)
}

#' Forecasting Lines and Region
#' 
#' This function draws h-step ahead forecasts their confidence (credible) region.
#' 
#' @param mapping Set of aesthetic mappings created by [ggplot2::aes()] or [ggplot2::aes_()].
#' If specified and `inherit.aes = TRUE` (the default), it is combined with the default mapping at the top level of the plot.
#' You must supply mapping if there is no plot mapping.
#' @param data The data to be displayed in this layer.
#' This should be generated by internal function, `gather_predbvhar()` that has specific columns.
#' @param stat The statistical transformation to use on the data for this layer, as a string.
#' @param position Position adjustment, either as a string, or the result of a call to a position adjustment function.
#' @param alpha_scale Scale of transparency parameter (`alpha`) between the two layers. `alpha` of CI ribbon = `alpha_scale` * `alpha` of path (By default, .3)
#' @param ci_param Parameter lists for [ggplot2::geom_ribbon()]. See details.
#' @param line_param Parameter lists for [ggplot2::geom_path()].
#' @param inherit.aes If `FALSE`, overrides the default aesthetics, rather than combining with them.
#' This is most useful for helper functions that define both data and aesthetics and shouldn't inherit behaviour from the default plot specification,
#' e.g. [ggplot2::borders()].
#' @details 
#' * `alpha_scale` is recommended to be smaller than 1.
#' * In case of `ci_param`, since it is controlling CI, try to specify `alpha` and `colour`.
#' @importFrom ggplot2 aes layer
#' @noRd
geom_predbvhar <- function(mapping = NULL, 
                           data = NULL, 
                           stat = "identity", 
                           position = "identity", 
                           alpha_scale = .3,
                           ci_param = list(),
                           line_param = list(),
                           inherit.aes = TRUE) {
  line_param$alpha <- ifelse(is.null(line_param$alpha), 1, line_param$alpha)
  ci_param$alpha <- line_param$alpha * alpha_scale
  ci_layer <- layer(
    geom = "ribbon",
    stat = "identity",
    data = data,
    mapping = aes(
      ymin = value_lower_joint, 
      ymax = value_upper_joint, 
      fill = Model
    ),
    position = position,
    params = ci_param,
    inherit.aes = inherit.aes
  )
  line_layer <- layer(
    geom = "path",
    stat = "identity",
    data = data,
    mapping = aes(colour = Model),
    position = position,
    params = line_param,
    inherit.aes = inherit.aes,
    show.legend = FALSE
  )
  list(ci_layer, line_layer)
}

#' Plot Forecast Result
#' 
#' Plots the forecasting result with forecast regions.
#' 
#' @param object A `predbvhar` object
#' @param type Divide variables using [ggplot2::facet_grid()] ("grid": default) or [ggplot2::facet_wrap()] ("wrap")
#' @param ci_fill color of CI
#' @param ci_alpha Transparency of CI
#' @param alpha_scale Scale of transparency parameter (`alpha`) between the two layers. `alpha` of CI ribbon = `alpha_scale` * `alpha` of path (By default, .5)
#' @param x_cut plot x axes from `x_cut` for visibility
#' @param viridis If `TRUE`, scale CI and forecast line using [ggplot2::scale_fill_viridis_d()] and [ggplot2::scale_colour_viridis_d], respectively.
#' @param viridis_option Option for viridis string. See `option` of [ggplot2::scale_colour_viridis_d]. Choose one of `c("A", "B", "C", "D", "E")`. By default, `D`.
#' @param NROW `nrow` of [ggplot2::facet_wrap()]
#' @param NCOL `ncol` of [ggplot2::facet_wrap()]
#' @param ... additional option for [ggplot2::geom_path()]
#' @return A ggplot object
#' @importFrom ggplot2 ggplot aes facet_grid geom_path labs element_blank scale_fill_viridis_d scale_colour_viridis_d
#' @importFrom dplyr filter
#' @export
autoplot.predbvhar <- function(object, 
                               type = c("grid", "wrap"), 
                               ci_alpha = .7,
                               alpha_scale = .3,
                               x_cut = 1,
                               viridis = FALSE,
                               viridis_option = "D",
                               NROW = NULL, 
                               NCOL = NULL, ...) {
  type <- match.arg(type)
  forecast_list <- 
    gather_predbvhar(object) |> 
    filter(id >= x_cut)
  p <- 
    forecast_list |> 
    ggplot(aes(x = id, y = value_forecast)) +
    geom_path(data = filter(forecast_list, forecast == FALSE))
  p <- 
    switch(
      type,
      "grid" = {
        p +
          geom_predbvhar(
            data = filter(forecast_list, forecast == TRUE),
            alpha_scale = alpha_scale,
            ci_param = list(alpha = ci_alpha, colour = NA),
            line_param = list(...)
          ) +
          facet_grid(variable ~ ., scales = "free_y") +
          labs(
            x = element_blank(),
            y = element_blank()
          )
      },
      "wrap" = {
        p +
          geom_predbvhar(
            data = filter(forecast_list, forecast == TRUE),
            alpha_scale = alpha_scale,
            ci_param = list(alpha = ci_alpha, colour = NA),
            line_param = list(...)
          ) +
          facet_wrap(variable ~ ., nrow = NROW, ncol = NCOL, scales = "free_y") +
          labs(
            x = element_blank(),
            y = element_blank()
          )
      }
    )
  # viridis----------------------------------
  if (viridis) {
    p +
      scale_fill_viridis_d(option = viridis_option) +
      scale_colour_viridis_d(option = viridis_option)
  } else {
    p
  }
}

#' @rdname autoplot.predbvhar
#' @return A ggplot layer
#' @importFrom dplyr bind_rows mutate filter
#' @importFrom ggplot2 ggplot aes facet_grid labs element_blank last_plot ggplot_build geom_path
#' @export
autolayer.predbvhar <- function(object, 
                                ci_fill = "grey70", 
                                ci_alpha = .5,
                                alpha_scale = .3, ...) {
  aes_data <-
    last_plot() |>
    ggplot_build()
  aes_data <- aes_data$plot$data # same form as forecast_list above
  x_cut <- aes_data[["id"]][1]
  NEW_list <- 
    gather_predbvhar(object) |> # new forecast_list
    filter(id >= x_cut)
  geom_predbvhar(
    data = NEW_list |> filter(forecast == TRUE),
    alpha_scale = alpha_scale,
    ci_param = list(alpha = ci_alpha, colour = NA),
    line_param = list(...)
  )
}

#' Adding Test Data Layer
#' 
#' This function adds a layer of test dataset.
#' 
#' @param data Test data to draw, which has the same format with the train data.
#' @param colour Color of the line (By default, `red`).
#' @param ... Other arguments passed on the [ggplot2::geom_path()].
#' @return A ggplot layer
#' @importFrom ggplot2 aes geom_path
#' @importFrom dplyr filter mutate n
#' @importFrom tidyr pivot_longer
#' @export
geom_eval <- function(data, colour = "red", ...) {
  if (is.matrix(data)) {
    data <- as.data.frame(data)
  }
  if (!is.data.frame(data)) {
    stop("'data' should be a data frame or matrix.")
  }
  if (!all(apply(data, 2, is.numeric))) {
    stop("Every column must be numeric class.")
  }
  aes_id <-
    last_plot() |>
    ggplot_build()
  aes_id <-
    aes_id$plot$data |>
    filter(forecast == FALSE)
  aes_id <- aes_id$id # index
  new_data <- 
    data |> 
    mutate(id = 1:n() + max(aes_id)) |> 
    pivot_longer(-id, names_to = "variable", values_to = "value")
  geom_path(
    aes(x = id, y = value),
    data = new_data,
    colour = colour,
    ...
  )
}

#' Make Data Form to Make Loss Plot
#' 
#' @param object List of `predbvhar`.
#' @param y_test Test data to be compared. should be the same format with the train data and predict$forecast.
#' @param loss Loss function to be used (`mse`: MSE, `mae`: MAE, `mape`: MAPE, `mase`: MASE)
#' @importFrom dplyr mutate bind_rows
#' @importFrom tidyr pivot_longer
#' @noRd
gather_loss <- function(object, y_test, loss = c("mse", "mae", "mape", "mase")) {
  loss <- match.arg(loss)
  if (missing(y_test)) {
    stop("Provide 'y_test' data")
  }
  # Must be list(predbvhar)-------------
  if (is.predbvhar(object)) {
    object <- list(object)
  }
  if (!all(sapply(object, is.predbvhar)) && !all(sapply(object, is.bvharcv))) {
    stop("'object' should be the list of 'predbvhar' or 'bvharcv'")
  }
  # Model names-------------------------
  mod_name <- 
    object |> 
    lapply(function(PRED) PRED$process) |> 
    unlist()
  # error for each model----------------
  score_dt <- 
    switch(
      loss,
      "mse" = {
        object |> 
          lapply(mse, y_test)
      },
      "mae" = {
        object |> 
          lapply(mae, y_test)
      },
      "mape" = {
        object |> 
          lapply(mape, y_test)
      },
      "mase" = {
        object |> 
          lapply(mase, y_test)
      }
    )
  score_dt |> 
    bind_rows() |> 
    mutate(Model = mod_name) |> 
    pivot_longer(-Model, names_to = "name", values_to = "score") |> 
    mutate(Loss = toupper(loss))
}

#' Compute Average Loss to Draw Horizontal Line
#' 
#' @param data Result made by `gather_loss`
#' @importFrom dplyr group_by mutate
#' @noRd
summarise_loss <- function(data) {
  data |> 
    group_by(Model, Loss) |> 
    mutate(average = mean(score))
}

#' Loss Lines
#' 
#' This function adds a layer of Loss given `predbvhar`.
#' 
#' @param mapping Set of aesthetic mappings created by [ggplot2::aes()] or [ggplot2::aes_()].
#' @param data The data to be displayed in this layer.
#' `predbvhar` object or list of `predbvhar`.
#' @param y_test Test data to be compared. should be the same format with the train data and predict$forecast.
#' @param loss Loss function to be used (`mse`: MSE, `mae`: MAE, `mape`: MAPE, `mase`: MASE)
#' @param mean_line Whether to draw average loss. By default, `FALSE`.
#' @param line_param Parameter lists for [ggplot2::geom_path()].
#' @param mean_param Parameter lists for average loss with [ggplot2::geom_hline()].
#' @param inherit.aes If `FALSE`, overrides the default aesthetics, rather than combining with them.
#' This is most useful for helper functions that define both data and aesthetics and shouldn't inherit behaviour from the default plot specification,
#' e.g. [ggplot2::borders()].
#' @param show.legend `logical`. Should this layer be included in the legend?
#' `NA`, the default, includes if any aesthetics are mapped.
#' `FALSE` never includes, and `TRUE` always includes.
#' It can also be a named logical vector to finely select the aesthetics to display.
#' @details 
#' Internal `gather_loss` function produces a `tibble` with columns named
#' 
#' \itemize{
#'   \item `Model` - The name of the model from `predbvhar$process`
#'   \item `name` - The variable name
#'   \item `score` - Values of losses
#' }
#'
#' Additionally, `summarise_loss` function adds grouped average of `score` by `Model` named `average`.
#' 
#' @seealso 
#' * [mse()] to compute MSE for given forecast result
#' * [mae()] to compute MAE for given forecast result
#' * [mape()] to compute MAPE for given forecast result
#' * [mase()] to compute MASE for given forecast result
#' @importFrom ggplot2 aes layer
#' @importFrom dplyr mutate bind_rows
#' @importFrom tidyr pivot_longer
#' @noRd
geom_loss <- function(mapping = NULL, 
                      data = NULL, 
                      y_test,
                      loss = c("mse", "mae", "mape", "mase"),
                      mean_line = FALSE,
                      line_param = list(),
                      mean_param = list(),
                      inherit.aes = TRUE,
                      show.legend = NA, ...) {
  loss_layer <- layer(
    geom = "line",
    stat = "identity",
    data = data,
    mapping = aes(
      x = name,
      y = score,
      colour = Model,
      group = Model
    ),
    position = "identity",
    params = line_param,
    inherit.aes = inherit.aes,
    show.legend = show.legend
  )
  if (mean_line) {
    mean_dt <- summarise_loss(data)
    mean_layer <- layer(
      geom = "hline",
      stat = "identity",
      data = mean_dt,
      mapping = aes(
        yintercept = average,
        colour = Model
      ),
      position = "identity",
      params = mean_param,
      inherit.aes = inherit.aes,
      show.legend = show.legend
    )
    list(mean_layer, loss_layer)
  } else {
    loss_layer
  }
}

#' Compare Lists of Models
#' 
#' Draw plot of test error for given models
#' 
#' @param mod_list Lists of forecast results (`predbvhar` objects)
#' @param y Test data to be compared. should be the same format with the train data and predict$forecast.
#' @param type Loss function to be used (`mse`: MSE, `mae`: MAE, `mape`: MAPE, `mase`: MASE)
#' @param mean_line Whether to draw average loss. By default, `FALSE`.
#' @param line_param Parameter lists for [ggplot2::geom_path()].
#' @param mean_param Parameter lists for average loss with [ggplot2::geom_hline()].
#' @param viridis If `TRUE`, scale CI and forecast line using [ggplot2::scale_fill_viridis_d()] and [ggplot2::scale_colour_viridis_d], respectively.
#' @param viridis_option Option for viridis string. See `option` of [ggplot2::scale_colour_viridis_d]. Choose one of `c("A", "B", "C", "D", "E")`. By default, `D`.
#' @param NROW `nrow` of [ggplot2::facet_wrap()]
#' @param NCOL `ncol` of [ggplot2::facet_wrap()]
#' @param ... Additional options for `geom_loss` (`inherit.aes` and `show.legend`)
#' @return A ggplot object
#' @seealso 
#' * [mse()] to compute MSE for given forecast result
#' * [mae()] to compute MAE for given forecast result
#' * [mape()] to compute MAPE for given forecast result
#' * [mase()] to compute MASE for given forecast result
#' @importFrom ggplot2 ggplot labs element_blank scale_colour_viridis_d facet_wrap
#' @importFrom dplyr bind_rows mutate
#' @export
gg_loss <- function(mod_list, 
                    y, 
                    type = c("mse", "mae", "mape", "mase"), 
                    mean_line = FALSE,
                    line_param = list(),
                    mean_param = list(),
                    viridis = FALSE, 
                    viridis_option = "D",
                    NROW = NULL,
                    NCOL = NULL, ...) {
  # Input data for geom_loss----------------
  fct_arrange <- toupper(type)
  data <- 
    type |> 
    lapply(
      function(loss) {
        gather_loss(mod_list, y, loss)
      }
    ) |> 
    bind_rows() |> 
    mutate(Loss = factor(Loss, levels = fct_arrange))
  # plot------------------------------------
  p <- 
    ggplot() +
    geom_loss(
      data = data,
      y_test = y,
      loss = type,
      mean_line = mean_line,
      line_param = line_param,
      mean_param = mean_param,
      ...
    ) +
    labs(
      x = element_blank(),
      y = element_blank()
    )
  # viridis--------------------------------
  if (viridis) {
    p <- 
      p +
      scale_colour_viridis_d(option = viridis_option)
  }
  # facet----------------------------------
  if (length(type) == 1) {
    p
  } else {
    p +
      facet_wrap(Loss ~ ., nrow = NROW, ncol = NCOL, scales = "free_y")
  }
}

#' Plot Impulse Responses
#' 
#' Draw impulse responses of response ~ impulse in the facet.
#' 
#' @param object A `bvharirf` object
#' @param ... Other arguments passed on the [ggplot2::geom_path()].
#' @seealso [irf()]
#' @return A ggplot object
#' @importFrom tidyr unite
#' @importFrom ggplot2 ggplot aes facet_wrap vars labs geom_path scale_x_continuous labs element_blank
#' @export
autoplot.bvharirf <- function(object, ...) {
  irf_df <- object$df_long
  irf_df |> 
    unite("term", impulse, response, c("impulse", "response"), sep = "->") |> 
    ggplot(aes(x = period, y = value)) +
    geom_path(...) +
    scale_x_continuous(breaks = 0:(object$lag_max)) +
    # facet_grid(response ~ impulse) + # y ~ x: impulse (x) -> response (y)
    facet_wrap(vars(term), scales = "free_y") + # better be transposed
    labs(
      x = element_blank(),
      y = element_blank()
    )
}

#' Plot the Result of BVAR and BVHAR MCMC
#' 
#' Draw BVAR and BVHAR MCMC plots.
#' 
#' @param object A `bvharsp` object
#' @param type The type of the plot. Posterior coefficient (`coef`), Trace plot (`trace`), kernel density plot (`dens`), and interval estimates plot (`area`).
#' @param pars Parameter names to draw.
#' @param regex_pars Regular expression parameter names to draw.
#' @param ... Other options for each [bayesplot::mcmc_trace()], [bayesplot::mcmc_dens()], and [bayesplot::mcmc_areas()].
#' @return A ggplot object
#' @importFrom bayesplot mcmc_trace mcmc_dens mcmc_areas
#' @export
autoplot.bvharsp <- function(object,
                             type = c("coef", "trace", "dens", "area"),
                             pars = character(),
                             regex_pars = character(), ...) {
  type <- match.arg(type)
  bayes_plt <- switch(
    type,
    "coef" = autoplot.summary.bvharsp(object, point = TRUE, ...),
    "trace" = mcmc_trace(x = object$param, pars = pars, regex_pars = regex_pars, ...),
    "dens" = mcmc_dens(x = object$param, pars = pars, regex_pars = regex_pars, ...),
    "area" = mcmc_areas(x = object$param, pars = pars, regex_pars = regex_pars, ...)
  )
  # additional processing later (title, labs)--------------
  bayes_plt
}

#' Make Data Form to Heatmap
#' 
#' @param object List of `summary.bvharsp`.
#' @importFrom tibble rownames_to_column
#' @importFrom tidyr pivot_longer separate_wider_regex
#' @importFrom dplyr mutate case_when
#' @noRd
gather_heat <- function(object) {
  heat_coef <-
    object$coefficients |>
    as.data.frame() |>
    rownames_to_column("term") |>
    pivot_longer(-term, names_to = "x", values_to = "value") |>
    mutate(x = factor(x, levels = colnames(object$coefficients)))
  is_vhar <- gsub(pattern = "(?=\\_).*", replacement = "", object$process, perl = TRUE) == "VHAR"
  if (object$type == "const") {
    heat_coef <-
      heat_coef |> 
      mutate(term = ifelse(
        term == "const",
        paste("const", term, sep = "_"),
        term
      ))
  }
  heat_coef <-
    heat_coef |> 
    separate_wider_regex(
      term,
      patterns = c(y = ".*", "_", ord = ".*$")
    ) |> 
    mutate(y = factor(y, levels = c(rev(colnames(object$coefficients)), "const")))
  # VHAR model--------------------------------
  if (is_vhar) {
    heat_coef <- 
      heat_coef |> 
      mutate(
        ord = case_when(
          ord == "day" ~ "Daily",
          ord == "week" ~ "Weekly",
          ord == "month" ~ "Monthly",
          .default = ord
        ),
        ord = factor(ord, levels = c("Daily", "Weekly", "Monthly", "const"))
      )
  }
  heat_coef
}

#' Plot the Heatmap of SSVS Coefficients
#' 
#' Draw heatmap for SSVS prior coefficients.
#' 
#' @param object A `summary.bvharsp` object
#' @param point Use point for sparsity representation
#' @param ... Other arguments passed on the [ggplot2::geom_tile()].
#' @return A ggplot object
#' @importFrom ggplot2 ggplot aes geom_tile geom_point scale_x_discrete guides guide_colourbar labs element_blank facet_grid scale_colour_gradient2
#' @export
autoplot.summary.bvharsp <- function(object, point = FALSE, ...) {
  heat_coef <- gather_heat(object)
  p <- 
    heat_coef |> 
    ggplot(aes(x = x, y = y))
  if (point) {
    p <- 
      p +
      geom_tile(fill = NA, colour = "#403d3d") +
      geom_point(aes(colour = value, size = abs(value))) +
      scale_colour_gradient2(low = "#d73027", mid = "#f0f0f0", high = "#4575b4") +
      guides(
        colour = guide_colourbar(title = element_blank()),
        size = "none"
      )
  } else {
    p <- 
      p +
      geom_tile(aes(fill = value), ...)
  }
  p <-
    p +
    scale_x_discrete(position = "top") +
    labs(
      x = element_blank(),
      y = element_blank()
    )
  if (object$p == 1) {
    return(p)
  }
  # plot for VHAR or p > 1-----------
  p +
    facet_grid(ord ~ ., switch = "x", scales = "free_y")
}