R/mean.detrended.dm.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

#' @title Calculate the mean detrended dendrometer series
#'
#' @description
#' Computes a mean detrended dendrometer series across multiple trees from the
#' output of \code{dm.detrend.fit()}.
#'
#' This is useful for creating one representative detrended series for a species,
#' site, or treatment group after detrending individual dendrometer series.
#'
#' Optionally, the function can:
#' \itemize{
#'   \item calculate a \strong{robust mean} using a trimmed mean across trees,
#'   \item remove temporal autocorrelation from the mean detrended series using
#'   \code{forecast::auto.arima()},
#'   \item rescale the autocorrelation-removed series so it stays non-negative
#'   and has mean = 1 within each vegetation season.
#' }
#'
#' @param detrended_dm An object of class \code{"dm_detrended"} returned by
#'   \code{dm.detrend.fit()}, or a data frame in the same wide format as
#'   \code{$detrended_data}.
#' @param series Optional character vector of tree/series names to include.
#'   Default is \code{NULL}, meaning all available detrended series are used.
#' @param ac1.remove Logical. If \code{TRUE}, removes temporal autocorrelation
#'   from the mean detrended series using \code{forecast::auto.arima()} applied
#'   separately within each season. Default is \code{TRUE}.
#' @param robust.mean Logical. If \code{TRUE}, calculates a trimmed mean across
#'   trees at each time step. Default is \code{TRUE}.
#' @param trim Proportion to trim from each tail when \code{robust.mean = TRUE}.
#'   Default is \code{0.15}.
#' @param seasonal_rescale Logical. If \code{TRUE}, the autocorrelation-removed
#'   series is shifted to non-negative values and rescaled to mean = 1 within
#'   each season. Default is \code{TRUE}.
#'
#' @return
#' A tibble of class \code{"mean_dm_detrended"} containing:
#' \itemize{
#'   \item metadata columns copied from the detrended input,
#'   \item \code{STD_DDM}: the mean detrended series,
#'   \item \code{RES_DDM}: the autocorrelation-removed mean detrended series
#'         (returned only when \code{ac1.remove = TRUE}).
#' }
#'
#' @examples
#' \donttest{
#' fit1 <- dm.growth.fit(
#'   df = gf_nepa17,
#'   TreeNum = 1:2,
#'   method = "gompertz",
#'   year_mode = "yearly",
#'   verbose = FALSE
#' )
#'
#' det1 <- dm.detrend.fit(fit1)
#'
#' m_det <- mean_detrended.dm(det1)
#' head(m_det, 10)
#' }
#'
#' @importFrom dplyr mutate group_by ungroup arrange select all_of %>%
#' @importFrom tibble as_tibble
#' @importFrom stats residuals
#' @export
mean_detrended.dm <- function(detrended_dm,
                              series = NULL,
                              ac1.remove = TRUE,
                              robust.mean = TRUE,
                              trim = 0.15,
                              seasonal_rescale = TRUE) {

  # -----------------------------
  # helpers
  # -----------------------------
  row_trimmed_mean <- function(x, trim = 0.15) {
    if (all(is.na(x))) return(NA_real_)
    mean(x, trim = trim, na.rm = TRUE)
  }

  row_simple_mean <- function(x) {
    if (all(is.na(x))) return(NA_real_)
    mean(x, na.rm = TRUE)
  }

  ac1_remove_one <- function(x) {
    if (all(is.na(x))) return(rep(NA_real_, length(x)))

    ok <- which(!is.na(x))
    if (length(ok) < 3) return(rep(NA_real_, length(x)))

    x_ok <- x[ok]

    fit <- tryCatch(
      forecast::auto.arima(x_ok),
      error = function(e) NULL
    )

    if (is.null(fit)) {
      return(rep(NA_real_, length(x)))
    }

    res <- rep(NA_real_, length(x))
    res[ok] <- as.numeric(stats::residuals(fit))
    res
  }

  rescale_positive_mean1 <- function(x, eps = 1e-6) {
    out <- rep(NA_real_, length(x))
    ok <- !is.na(x)

    if (!any(ok)) return(out)

    x_ok <- x[ok]
    xmin <- min(x_ok, na.rm = TRUE)
    shifted <- x_ok - xmin + eps
    m <- mean(shifted, na.rm = TRUE)

    if (!is.finite(m) || m <= 0) {
      out[ok] <- NA_real_
      return(out)
    }

    out[ok] <- shifted / m
    out
  }

  # -----------------------------
  # input
  # -----------------------------
  if (inherits(detrended_dm, "dm_detrended")) {
    dat <- tibble::as_tibble(detrended_dm$detrended_data)
  } else if (is.data.frame(detrended_dm)) {
    dat <- tibble::as_tibble(detrended_dm)
  } else {
    stop(
      "'detrended_dm' must be either:\n",
      "  1. an object of class 'dm_detrended', or\n",
      "  2. a data frame in the same format as $detrended_data."
    )
  }

  meta_cols <- intersect(
    c("TIME", "doy", "season_label", "season_start", "season_end", "season_day"),
    names(dat)
  )

  value_cols <- setdiff(names(dat), meta_cols)

  if (length(value_cols) == 0) {
    stop("No detrended series columns were found.")
  }

  if (!is.null(series)) {
    miss <- setdiff(series, value_cols)
    if (length(miss) > 0) {
      warning("These requested series were not found and were ignored: ",
              paste(miss, collapse = ", "))
    }
    value_cols <- intersect(value_cols, series)
    if (length(value_cols) == 0) {
      stop("None of the requested series were found.")
    }
  }

  if (!requireNamespace("forecast", quietly = TRUE) && isTRUE(ac1.remove)) {
    stop("Package 'forecast' is required when ac1.remove = TRUE.")
  }

  # -----------------------------
  # mean detrended series
  # -----------------------------
  mat <- as.matrix(dat[, value_cols, drop = FALSE])

  if (isTRUE(robust.mean)) {
    mean_series <- apply(mat, 1, row_trimmed_mean, trim = trim)
  } else {
    mean_series <- apply(mat, 1, row_simple_mean)
  }

  out <- dat[, meta_cols, drop = FALSE]
  out$STD_DDM <- as.numeric(mean_series)

  # -----------------------------
  # remove autocorrelation
  # -----------------------------
  if (isTRUE(ac1.remove)) {

    if ("season_label" %in% names(out)) {
      out <- out %>%
        dplyr::group_by(.data$season_label) %>%
        dplyr::arrange(.data$TIME, .by_group = TRUE) %>%
        dplyr::mutate(
          RES_DDM = ac1_remove_one(.data$STD_DDM)
        ) %>%
        dplyr::ungroup()
    } else {
      out$RES_DDM <- ac1_remove_one(out$STD_DDM)
    }

    if (isTRUE(seasonal_rescale)) {
      if ("season_label" %in% names(out)) {
        out <- out %>%
          dplyr::group_by(.data$season_label) %>%
          dplyr::mutate(
            RES_DDM = rescale_positive_mean1(.data$RES_DDM)
          ) %>%
          dplyr::ungroup()
      } else {
        out$RES_DDM <- rescale_positive_mean1(out$RES_DDM)
      }
    }
  }

  class(out) <- c("mean_dm_detrended", class(out))
  out
}
###----------------------------------------------------------------------------------------------

#' @title Plot mean detrended dendrometer series
#'
#' @description
#' S3 plotting method for objects returned by \code{mean_detrended.dm()}.
#'
#' It can plot:
#' \itemize{
#'   \item the mean detrended series (\code{STD_DDM}),
#'   \item the autocorrelation-removed mean detrended series (\code{RES_DDM}),
#'   \item or both together.
#' }
#'
#' @param x An object of class \code{"mean_dm_detrended"} returned by
#'   \code{mean_detrended.dm()}.
#' @param y Unused.
#' @param type Plot type. One of:
#'   \describe{
#'     \item{`"series"`}{Plot the time series over calendar time, DOY, or season day.}
#'     \item{`"seasonal"`}{Overlay seasons on a season-day scale.}
#'     \item{`"boxplot"`}{Show boxplots of the detrended values by season or metric.}
#'   }
#' @param value Which variable to plot. One of:
#'   \describe{
#'     \item{`"both"`}{Plot both \code{STD_DDM} and \code{RES_DDM} when available.}
#'     \item{`"STD_DDM"`}{Plot only the mean detrended series.}
#'     \item{`"RES_DDM"`}{Plot only the autocorrelation-removed mean detrended series.}
#'   }
#' @param seasons Optional character vector of \code{season_label} values to retain.
#' @param x_axis Character string controlling the x-axis. One of:
#'   \describe{
#'     \item{`"default"`}{Uses calendar date for \code{type = "series"} and
#'       season day for \code{type = "seasonal"}.}
#'     \item{`"date"`}{Use actual calendar date.}
#'     \item{`"doy"`}{Use calendar day-of-year.}
#'     \item{`"season_day"`}{Use vegetation season day.}
#'   }
#' @param facet_by Character string controlling faceting. One of:
#'   \describe{
#'     \item{`"none"`}{No faceting.}
#'     \item{`"season"`}{Facet by \code{season_label}.}
#'     \item{`"metric"`}{Facet by metric (\code{STD_DDM}, \code{RES_DDM}).}
#'   }
#' @param ncol Optional integer giving the number of facet columns.
#' @param box_group For \code{type = "boxplot"}, grouping variable on the x-axis.
#'   One of \code{"metric"} or \code{"season"}.
#' @param alpha Numeric alpha transparency for lines/points.
#' @param line_width Numeric line width.
#' @param point_size Numeric point size.
#' @param legend_position Character legend position passed to ggplot2.
#' @param ... Further arguments passed to or from other methods.
#'
#' @return
#' A \code{ggplot2} object.
#'
#' @method plot mean_dm_detrended
#' @importFrom dplyr filter mutate %>%
#' @importFrom tidyr pivot_longer
#' @importFrom tibble as_tibble
#' @importFrom ggplot2 ggplot aes geom_line geom_point geom_hline geom_boxplot
#'   geom_jitter facet_wrap theme_bw theme element_text labs
#' @export
plot.mean_dm_detrended <- function(x,
                                   y = NULL,
                                   type = c("series", "seasonal", "boxplot"),
                                   value = c("both", "STD_DDM", "RES_DDM"),
                                   seasons = NULL,
                                   x_axis = c("default", "date", "doy", "season_day"),
                                   facet_by = c("none", "season", "metric"),
                                   ncol = NULL,
                                   box_group = c("metric", "season"),
                                   alpha = 0.8,
                                   line_width = 0.8,
                                   point_size = 1.4,
                                   legend_position = "right",
                                   ...) {

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

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

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

  dat <- tibble::as_tibble(x)

  if (!"TIME" %in% names(dat)) {
    stop("The input object must contain a 'TIME' column.")
  }

  if (!is.null(seasons)) {
    if (!"season_label" %in% names(dat)) {
      stop("The object has no 'season_label' column for seasonal filtering.")
    }
    dat <- dat %>% dplyr::filter(.data$season_label %in% seasons)
  }

  metric_cols <- intersect(c("STD_DDM", "RES_DDM"), names(dat))

  if (length(metric_cols) == 0) {
    stop("No plottable columns found. Expected 'STD_DDM' and/or 'RES_DDM'.")
  }

  if (value == "STD_DDM") {
    if (!"STD_DDM" %in% metric_cols) {
      stop("Column 'STD_DDM' is not available.")
    }
    metric_cols <- "STD_DDM"
  }

  if (value == "RES_DDM") {
    if (!"RES_DDM" %in% metric_cols) {
      stop("Column 'RES_DDM' is not available.")
    }
    metric_cols <- "RES_DDM"
  }

  dat$TIME_DATE <- as.Date(dat$TIME)

  if (!"doy" %in% names(dat)) {
    dat$doy <- lubridate::yday(dat$TIME_DATE)
  }
  if (!"season_day" %in% names(dat)) {
    dat$season_day <- seq_len(nrow(dat))
  }

  long_dat <- dat %>%
    tidyr::pivot_longer(
      cols = dplyr::all_of(metric_cols),
      names_to = "metric",
      values_to = "value"
    )

  if (nrow(long_dat) == 0) {
    stop("No data available to plot.")
  }

  if (x_axis == "default") {
    if (type == "seasonal") {
      x_var <- "season_day"
      x_lab <- "Season day"
    } else {
      x_var <- "TIME_DATE"
      x_lab <- "Date"
    }
  } else if (x_axis == "date") {
    x_var <- "TIME_DATE"
    x_lab <- "Date"
  } else if (x_axis == "doy") {
    x_var <- "doy"
    x_lab <- "DOY"
  } else {
    x_var <- "season_day"
    x_lab <- "Season day"
  }

  if (type == "series") {
    if (facet_by == "none") {
      p <- ggplot2::ggplot(
        long_dat,
        ggplot2::aes(
          x = .data[[x_var]],
          y = .data$value,
          colour = .data$metric,
          group = .data$metric
        )
      ) +
        ggplot2::geom_hline(yintercept = 1, linetype = 2) +
        ggplot2::geom_line(alpha = alpha, linewidth = line_width) +
        ggplot2::geom_point(alpha = alpha, size = point_size)
    } else if (facet_by == "season") {
      if (!"season_label" %in% names(long_dat)) {
        stop("The object has no 'season_label' column for faceting by season.")
      }
      p <- ggplot2::ggplot(
        long_dat,
        ggplot2::aes(
          x = .data[[x_var]],
          y = .data$value,
          colour = .data$metric,
          group = .data$metric
        )
      ) +
        ggplot2::geom_hline(yintercept = 1, linetype = 2) +
        ggplot2::geom_line(alpha = alpha, linewidth = line_width) +
        ggplot2::geom_point(alpha = alpha, size = point_size) +
        ggplot2::facet_wrap(stats::as.formula("~ season_label"), scales = "free_x", ncol = ncol)
    } else {
      p <- ggplot2::ggplot(
        long_dat,
        ggplot2::aes(
          x = .data[[x_var]],
          y = .data$value,
          colour = .data$metric,
          group = .data$metric
        )
      ) +
        ggplot2::geom_hline(yintercept = 1, linetype = 2) +
        ggplot2::geom_line(alpha = alpha, linewidth = line_width) +
        ggplot2::geom_point(alpha = alpha, size = point_size) +
        ggplot2::facet_wrap(stats::as.formula("~ metric"), scales = "free_y", ncol = ncol)
    }

    p <- p +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = legend_position,
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 11)
      ) +
      ggplot2::labs(
        x = x_lab,
        y = "Mean detrended value",
        colour = "Metric",
        title = "Mean detrended dendrometer series"
      )

    return(p)
  }

  if (type == "seasonal") {
    if (!"season_label" %in% names(long_dat)) {
      stop("The object has no 'season_label' column required for type = 'seasonal'.")
    }

    if (facet_by == "metric") {
      p <- ggplot2::ggplot(
        long_dat,
        ggplot2::aes(
          x = .data$season_day,
          y = .data$value,
          colour = .data$season_label,
          group = .data$season_label
        )
      ) +
        ggplot2::geom_hline(yintercept = 1, linetype = 2) +
        ggplot2::geom_line(alpha = alpha, linewidth = line_width) +
        ggplot2::facet_wrap(stats::as.formula("~ metric"), scales = "free_y", ncol = ncol) +
        ggplot2::labs(colour = "Season")
    } else if (facet_by == "season") {
      p <- ggplot2::ggplot(
        long_dat,
        ggplot2::aes(
          x = .data$season_day,
          y = .data$value,
          colour = .data$metric,
          group = .data$metric
        )
      ) +
        ggplot2::geom_hline(yintercept = 1, linetype = 2) +
        ggplot2::geom_line(alpha = alpha, linewidth = line_width) +
        ggplot2::facet_wrap(stats::as.formula("~ season_label"), scales = "free_y", ncol = ncol) +
        ggplot2::labs(colour = "Metric")
    } else {
      long_dat$colour_group <- interaction(long_dat$metric, long_dat$season_label, drop = TRUE)

      p <- ggplot2::ggplot(
        long_dat,
        ggplot2::aes(
          x = .data$season_day,
          y = .data$value,
          colour = .data$colour_group,
          group = .data$colour_group
        )
      ) +
        ggplot2::geom_hline(yintercept = 1, linetype = 2) +
        ggplot2::geom_line(alpha = alpha, linewidth = line_width) +
        ggplot2::labs(colour = "Metric | Season")
    }

    p <- p +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = legend_position,
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 11)
      ) +
      ggplot2::labs(
        x = "Season day",
        y = "Mean detrended value",
        title = "Seasonal overlay of mean detrended series"
      )

    return(p)
  }

  if (type == "boxplot") {
    if (box_group == "metric") {
      long_dat$group_var <- long_dat$metric
      x_lab2 <- "Metric"
    } else {
      if (!"season_label" %in% names(long_dat)) {
        stop("The object has no 'season_label' column required for boxplots by season.")
      }
      long_dat$group_var <- long_dat$season_label
      x_lab2 <- "Season"
    }

    p <- ggplot2::ggplot(
      long_dat[is.finite(long_dat$value), , drop = FALSE],
      ggplot2::aes(x = .data$group_var, y = .data$value, fill = .data$metric)
    ) +
      ggplot2::geom_boxplot(outlier.shape = NA) +
      ggplot2::geom_jitter(width = 0.15, alpha = alpha, size = point_size) +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = legend_position,
        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 = "Mean detrended value",
        fill = "Metric",
        title = "Distribution of mean detrended series"
      )

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

    return(p)
  }

  stop("Unknown plot type.")
}