R/plot_dm.detrend.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

#' Plot detrended dendrometer series
#'
#' @description
#' S3 plotting method for objects returned by [dm.detrend.fit()].
#'
#' The default plot compares:
#' \itemize{
#'   \item original daily dendrometer series,
#'   \item fitted curve reconstructed on the original daily scale,
#'   \item residuals (`observed - fitted_original`),
#'   \item detrended standardized series.
#' }
#'
#' @param x An object of class \code{"dm_detrended"} returned by
#'   [dm.detrend.fit()].
#' @param y Unused.
#' @param type Plot type. One of:
#'   \describe{
#'     \item{`"compare"`}{Default three-panel comparison of original daily vs
#'       fitted, residuals, and detrended standardized series.}
#'     \item{`"fit"`}{Original daily and fitted original-scale values only.}
#'     \item{`"residual"`}{Residuals only.}
#'     \item{`"detrended"`}{Detrended standardized series only.}
#'     \item{`"boxplot"`}{Distribution of detrended standardized values by series
#'       or season.}
#'   }
#' @param series Optional character vector of dendrometer series to plot.
#'   Default is \code{NULL}, meaning all available series are used.
#' @param seasons Optional character vector of vegetation-season labels to plot.
#'   Default is \code{NULL}, meaning all seasons are used.
#' @param x_axis Character string controlling the x-axis. One of:
#'   \describe{
#'     \item{`"default"`}{Uses date for `"compare"` and `"fit"`, and season day
#'       for `"residual"` and `"detrended"`.}
#'     \item{`"date"`}{Use actual calendar date.}
#'     \item{`"season_day"`}{Use vegetation-season day.}
#'     \item{`"doy"`}{Use calendar day-of-year.}
#'   }
#' @param facet_by Character string controlling faceting. One of:
#'   \describe{
#'     \item{`"series"`}{Facet by dendrometer series.}
#'     \item{`"season"`}{Facet by vegetation season.}
#'     \item{`"none"`}{No faceting.}
#'   }
#' @param ncol Optional integer giving the number of columns in faceted plots
#'   where [ggplot2::facet_wrap()] is used.
#' @param box_group For \code{type = "boxplot"}, grouping variable on the x-axis.
#'   One of \code{"series"} or \code{"season"}.
#' @param show_observed Logical. If \code{TRUE}, original daily observations are
#'   shown in plot types where relevant. Default is \code{TRUE}.
#' @param show_fitted Logical. If \code{TRUE}, fitted original-scale values are
#'   shown in plot types where relevant. Default is \code{TRUE}.
#' @param point_alpha Numeric alpha level used for observed points. Default is
#'   \code{0.7}.
#' @param line_width Numeric line width used for fitted, residual, and detrended
#'   lines. Default is \code{0.8}.
#' @param legend_position Character string specifying legend position.
#'   Default is \code{"right"}.
#' @param ... Further arguments passed to or from other methods.
#'
#' @return
#' A \code{ggplot2} object.
#'
#' @examples
#' \donttest{
#' fit1 <- dm.growth.fit(
#'   df = gf_nepa17,
#'   TreeNum = 1:2,
#'   method = "gompertz",
#'   year_mode = "yearly",
#'   verbose = FALSE
#' )
#'
#' det1 <- dm.detrend.fit(fit1)
#'
#' plot(det1)
#' plot(det1, type = "fit")
#' plot(det1, type = "residual")
#' plot(det1, type = "detrended")
#' plot(det1, type = "boxplot")
#' plot(det1, type = "compare", facet_by = "series")
#' plot(det1, type = "compare", facet_by = "season")
#' }
#'
#' @method plot dm_detrended
#' @importFrom dplyr filter mutate select %>%
#' @importFrom tibble as_tibble
#' @importFrom ggplot2 ggplot aes geom_point geom_line geom_hline geom_boxplot
#'   geom_jitter facet_wrap facet_grid theme_bw theme element_text labs
#' @export
plot.dm_detrended <- function(x,
                              y = NULL,
                              type = c("compare", "fit", "residual", "detrended", "boxplot"),
                              series = NULL,
                              seasons = NULL,
                              x_axis = c("default", "date", "season_day", "doy"),
                              facet_by = c("series", "season", "none"),
                              ncol = NULL,
                              box_group = c("series", "season"),
                              show_observed = TRUE,
                              show_fitted = TRUE,
                              point_alpha = 0.7,
                              line_width = 0.8,
                              legend_position = "right",
                              ...) {

  TIME <- doy <- season_day <- season_label <- series_name <- observed <- fitted_original <- NULL
  residual <- detrended_std <- x_plot <- facet_var <- colour_group <- group_id <- panel <- group_var <- NULL

  type <- match.arg(type)
  x_axis <- match.arg(x_axis)
  facet_by <- match.arg(facet_by)
  box_group <- match.arg(box_group)

  if (!inherits(x, "dm_detrended")) {
    stop("'x' must be an object of class 'dm_detrended'.")
  }

  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("Package 'ggplot2' is required for plot.dm_detrended().")
  }

  if (is.null(x$detrended_long) || !is.data.frame(x$detrended_long)) {
    stop("The object does not contain a valid 'detrended_long' table.")
  }

  dat <- tibble::as_tibble(x$detrended_long)

  if (!is.null(series)) {
    dat <- dat %>% dplyr::filter(.data$series %in% series)
  }
  if (!is.null(seasons)) {
    dat <- dat %>% dplyr::filter(.data$season_label %in% seasons)
  }

  if (nrow(dat) == 0) {
    stop("No rows remain after filtering by series and/or seasons.")
  }

  dat$season_label <- as.character(dat$season_label)
  dat$TIME_date <- as.Date(dat$TIME)

  x_cfg <- if (x_axis == "default") {
    if (type %in% c("compare", "fit")) {
      list(var = "TIME_date", label = "Date")
    } else {
      list(var = "season_day", label = "Season day")
    }
  } else if (x_axis == "date") {
    list(var = "TIME_date", label = "Date")
  } else if (x_axis == "doy") {
    list(var = "doy", label = "DOY")
  } else {
    list(var = "season_day", label = "Season day")
  }

  x_var <- x_cfg$var
  x_lab <- x_cfg$label

  grp_cfg <- dmdp_prepare_groups(dat, facet_by = facet_by)
  dat <- grp_cfg$data

  if (type == "compare") {
    p <- ggplot2::ggplot() +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = if (isTRUE(grp_cfg$use_colour)) legend_position else "none",
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 11)
      )

    obsfit_dat <- dat
    obsfit_dat$panel <- "Original vs fitted"

    resid_dat <- dat
    resid_dat$panel <- "Residual"

    detr_dat <- dat
    detr_dat$panel <- "Detrended standardized"

    if (isTRUE(show_observed)) {
      p <- p + ggplot2::geom_point(
        data = obsfit_dat[is.finite(obsfit_dat$observed), , drop = FALSE],
        mapping = ggplot2::aes(
          x = .data[[x_var]],
          y = .data$observed,
          colour = .data$colour_group,
          group = .data$group_id
        ),
        alpha = point_alpha,
        size = 1.4
      )
    }

    if (isTRUE(show_fitted)) {
      p <- p + ggplot2::geom_line(
        data = obsfit_dat[is.finite(obsfit_dat$fitted_original), , drop = FALSE],
        mapping = ggplot2::aes(
          x = .data[[x_var]],
          y = .data$fitted_original,
          colour = .data$colour_group,
          group = .data$group_id
        ),
        linewidth = line_width
      )
    }

    p <- p +
      ggplot2::geom_hline(
        data = data.frame(panel = "Residual"),
        mapping = ggplot2::aes(yintercept = 0),
        inherit.aes = FALSE,
        linetype = 2
      ) +
      ggplot2::geom_line(
        data = resid_dat[is.finite(resid_dat$residual), , drop = FALSE],
        mapping = ggplot2::aes(
          x = .data[[x_var]],
          y = .data$residual,
          colour = .data$colour_group,
          group = .data$group_id
        ),
        linewidth = line_width
      ) +
      ggplot2::geom_hline(
        data = data.frame(panel = "Detrended standardized"),
        mapping = ggplot2::aes(yintercept = 1),
        inherit.aes = FALSE,
        linetype = 2
      ) +
      ggplot2::geom_line(
        data = detr_dat[is.finite(detr_dat$detrended_std), , drop = FALSE],
        mapping = ggplot2::aes(
          x = .data[[x_var]],
          y = .data$detrended_std,
          colour = .data$colour_group,
          group = .data$group_id
        ),
        linewidth = line_width
      ) +
      ggplot2::labs(
        x = x_lab,
        y = NULL,
        colour = grp_cfg$legend_title,
        title = "Original daily, fitted, residual, and detrended standardized series"
      )

    if (facet_by == "none") {
      p <- p + ggplot2::facet_grid(panel ~ ., scales = "free_y")
    } else {
      p <- p + ggplot2::facet_grid(panel ~ facet_var, scales = "free_y")
    }

    return(p)
  }

  if (type == "fit") {
    p <- ggplot2::ggplot() +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = if (isTRUE(grp_cfg$use_colour)) legend_position else "none",
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 11)
      )

    if (isTRUE(show_observed)) {
      p <- p + ggplot2::geom_point(
        data = dat[is.finite(dat$observed), , drop = FALSE],
        mapping = ggplot2::aes(
          x = .data[[x_var]],
          y = .data$observed,
          colour = .data$colour_group,
          group = .data$group_id
        ),
        alpha = point_alpha,
        size = 1.4
      )
    }

    if (isTRUE(show_fitted)) {
      p <- p + ggplot2::geom_line(
        data = dat[is.finite(dat$fitted_original), , drop = FALSE],
        mapping = ggplot2::aes(
          x = .data[[x_var]],
          y = .data$fitted_original,
          colour = .data$colour_group,
          group = .data$group_id
        ),
        linewidth = line_width
      )
    }

    p <- p + ggplot2::labs(
      x = x_lab,
      y = "Growth",
      colour = grp_cfg$legend_title,
      title = "Original daily and fitted original-scale growth"
    )

    if (isTRUE(grp_cfg$facet_enabled)) {
      p <- p + ggplot2::facet_wrap(stats::as.formula("~ facet_var"), scales = "free_y", ncol = ncol)
    }

    return(p)
  }

  if (type == "residual") {
    p <- ggplot2::ggplot(
      dat[is.finite(dat$residual), , drop = FALSE],
      ggplot2::aes(
        x = .data[[x_var]],
        y = .data$residual,
        colour = .data$colour_group,
        group = .data$group_id
      )
    ) +
      ggplot2::geom_hline(yintercept = 0, linetype = 2) +
      ggplot2::geom_line(linewidth = line_width) +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = if (isTRUE(grp_cfg$use_colour)) legend_position else "none",
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 11)
      ) +
      ggplot2::labs(
        x = x_lab,
        y = "Residual",
        colour = grp_cfg$legend_title,
        title = "Residuals from original daily and fitted original-scale growth"
      )

    if (isTRUE(grp_cfg$facet_enabled)) {
      p <- p + ggplot2::facet_wrap(stats::as.formula("~ facet_var"), scales = "free_y", ncol = ncol)
    }

    return(p)
  }

  if (type == "detrended") {
    p <- ggplot2::ggplot(
      dat[is.finite(dat$detrended_std), , drop = FALSE],
      ggplot2::aes(
        x = .data[[x_var]],
        y = .data$detrended_std,
        colour = .data$colour_group,
        group = .data$group_id
      )
    ) +
      ggplot2::geom_hline(yintercept = 1, linetype = 2) +
      ggplot2::geom_line(linewidth = line_width) +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = if (isTRUE(grp_cfg$use_colour)) legend_position else "none",
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 11)
      ) +
      ggplot2::labs(
        x = x_lab,
        y = "Detrended standardized",
        colour = grp_cfg$legend_title,
        title = "Detrended standardized dendrometer series"
      )

    if (isTRUE(grp_cfg$facet_enabled)) {
      p <- p + ggplot2::facet_wrap(stats::as.formula("~ facet_var"), scales = "free_y", ncol = ncol)
    }

    return(p)
  }

  if (type == "boxplot") {
    if (box_group == "series") {
      dat$group_var <- dat$series
      x_lab2 <- "Series"
    } else {
      dat$group_var <- dat$season_label
      x_lab2 <- "Season"
    }

    p <- ggplot2::ggplot(
      dat[is.finite(dat$detrended_std), , drop = FALSE],
      ggplot2::aes(x = .data$group_var, y = .data$detrended_std)
    ) +
      ggplot2::geom_boxplot(outlier.shape = NA) +
      ggplot2::geom_jitter(width = 0.15, alpha = point_alpha, size = 1.2) +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = "none",
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 11),
        axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)
      ) +
      ggplot2::labs(
        x = x_lab2,
        y = "Detrended standardized",
        title = "Distribution of detrended standardized values"
      )

    if (facet_by == "series" && box_group == "season") {
      p <- p + ggplot2::facet_wrap(stats::as.formula("~ series"), scales = "free_y", ncol = ncol)
    } else if (facet_by == "season" && box_group == "series") {
      p <- p + ggplot2::facet_wrap(stats::as.formula("~ season_label"), scales = "free_y", ncol = ncol)
    }

    return(p)
  }

  stop("Unknown plot type.")
}


#' @keywords internal
dmdp_prepare_groups <- function(dat, facet_by = c("series", "season", "none")) {
  facet_by <- match.arg(facet_by)

  dat$group_id <- interaction(dat$series, dat$season_label, drop = TRUE)

  n_series <- length(unique(dat$series))
  n_seasons <- length(unique(dat$season_label))

  if (facet_by == "series") {
    dat$facet_var <- dat$series
    if (n_seasons > 1) {
      dat$colour_group <- factor(dat$season_label)
      use_colour <- TRUE
      legend_title <- "Season"
    } else {
      dat$colour_group <- factor("curve")
      use_colour <- FALSE
      legend_title <- NULL
    }
    return(list(
      data = dat,
      facet_enabled = TRUE,
      use_colour = use_colour,
      legend_title = legend_title
    ))
  }

  if (facet_by == "season") {
    dat$facet_var <- dat$season_label
    if (n_series > 1) {
      dat$colour_group <- factor(dat$series)
      use_colour <- TRUE
      legend_title <- "Series"
    } else {
      dat$colour_group <- factor("curve")
      use_colour <- FALSE
      legend_title <- NULL
    }
    return(list(
      data = dat,
      facet_enabled = TRUE,
      use_colour = use_colour,
      legend_title = legend_title
    ))
  }

  dat$facet_var <- "All"
  if (n_series > 1 && n_seasons > 1) {
    dat$colour_group <- factor(paste(dat$series, dat$season_label, sep = " | "))
    use_colour <- TRUE
    legend_title <- "Series | Season"
  } else if (n_series > 1) {
    dat$colour_group <- factor(dat$series)
    use_colour <- TRUE
    legend_title <- "Series"
  } else if (n_seasons > 1) {
    dat$colour_group <- factor(dat$season_label)
    use_colour <- TRUE
    legend_title <- "Season"
  } else {
    dat$colour_group <- factor("curve")
    use_colour <- FALSE
    legend_title <- NULL
  }

  list(
    data = dat,
    facet_enabled = FALSE,
    use_colour = use_colour,
    legend_title = legend_title
  )
}