R/mov.cor.dm.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

#' @title Running correlation between dendrometer data and climate
#'
#' @description
#' Calculates running correlations between a selected daily dendrometer summary
#' and one or more climate variables. The user can select the daily dendrometer
#' statistic, correlation method, optional bootstrap confidence intervals, and
#' lagged / antecedent climate transformations.
#'
#' @details
#' The dendrometer series is first aggregated to daily resolution. The daily
#' dendrometer statistic used for correlation is controlled by \code{dm_stat}:
#' \itemize{
#'   \item \code{"mean"}: daily mean dendrometer value
#'   \item \code{"min"}: daily minimum
#'   \item \code{"max"}: daily maximum
#'   \item \code{"median"}: daily median
#'   \item \code{"amplitude"}: daily amplitude (\code{max - min})
#'   \item \code{"change"}: day-to-day change in the daily mean
#' }
#'
#' Users can choose the climate variables to analyze via \code{clim_vars}.
#' Climate transformation settings can be given as:
#' \itemize{
#'   \item a single value applied to all selected climate variables
#'   \item an unnamed vector with one value per selected climate variable
#'   \item a named vector mapping each selected climate variable to its own setting
#' }
#'
#' This applies to \code{clim_fun}, \code{lag_days}, and \code{accum_days}.
#'
#' @param df A data frame with the first column containing date-time in the format
#'   \code{yyyy-mm-dd HH:MM:SS} (or convertible to \code{POSIXct}) and dendrometer
#'   series in the following columns.
#' @param Clim A data frame with the first column containing daily date
#'   (\code{yyyy-mm-dd}, \code{Date}, or convertible) and one or more climate
#'   variables in the following columns.
#' @param TreeNum Integer indicating the dendrometer series to analyze.
#' @param win_size Integer giving the running window size in days. Minimum is 18.
#' @param cor_method Correlation method: one of \code{"pearson"},
#'   \code{"kendall"}, or \code{"spearman"}.
#' @param boot Logical. If \code{TRUE}, bootstrap confidence intervals are
#'   computed for each running correlation.
#' @param R Integer number of bootstrap iterations.
#' @param boot.ci Numeric confidence level selector: one of \code{0.01},
#'   \code{0.05} (default), or \code{0.1}.
#' @param set_seed Integer seed for reproducibility of bootstrap results.
#' @param dm_stat Daily dendrometer statistic used for correlation. One of
#'   \code{"mean"}, \code{"min"}, \code{"max"}, \code{"median"},
#'   \code{"amplitude"}, or \code{"change"}.
#' @param clim_vars Optional character vector of climate variables to analyze.
#'   If \code{NULL}, all numeric climate variables are used.
#' @param lag_days Climate lag in days. Can be:
#'   \itemize{
#'     \item one value for all selected climate variables
#'     \item one value per selected climate variable
#'     \item a named numeric vector keyed by climate variable name
#'   }
#' @param accum_days Antecedent window length in days for climate transformation.
#'   Can be scalar, per-variable, or a named numeric vector.
#' @param clim_fun Climate transformation over the antecedent window. One of
#'   \code{"raw"}, \code{"mean"}, \code{"sum"}, \code{"max"},
#'   \code{"min"}, or \code{"median"}. Can be scalar, per-variable, or a named
#'   character vector.
#' @param min_complete Minimum number of complete paired observations required in
#'   a running window to calculate correlation. If \code{NULL}, defaults to
#'   \code{max(3, ceiling(0.8 * win_size))}.
#' @param p_adjust_method Method for p-value adjustment in the non-bootstrap
#'   output. Passed to \code{p.adjust()}. Use \code{"none"} to skip adjustment.
#'
#' @return
#' A list with class \code{"mov_cor_dm"} (and \code{"mov_cor_dm_boot"} if
#' bootstrapped) containing:
#' \itemize{
#'   \item \code{results}: named list of tibbles, one per climate variable
#'   \item \code{metadata}: analysis metadata
#'   \item \code{call}: the matched function call
#' }
#'
#' @examples
#' \donttest{
#' library(dendRoAnalyst)
#' data(gf_nepa17)
#' data(ktm_rain17)
#'
#' # one common climate transformation for all selected variables
#' out_corr <- mov.cor.dm(
#'   df = gf_nepa17,
#'   Clim = ktm_rain17,
#'   TreeNum = 1,
#'   win_size = 21,
#'   clim_fun = "raw"
#' )
#' print(out_corr)
#' summary(out_corr)
#'
#' # variable-specific climate transformations
#' out_varfun <- mov.cor.dm(
#'   df = gf_nepa17,
#'   Clim = ktm_rain17,
#'   TreeNum = 1,
#'   win_size = 21,
#'   clim_vars = c("rainfall"),
#'   clim_fun = c(rainfall = "sum"),
#'   lag_days = c(rainfall = 1),
#'   accum_days = c(rainfall = 7)
#' )
#'
#' # bootstrap confidence intervals
#' out_boot <- mov.cor.dm(
#'   df = gf_nepa17,
#'   Clim = ktm_rain17,
#'   TreeNum = 1,
#'   win_size = 21,
#'   boot = TRUE,
#'   R = 250
#' )
#' summary(out_boot)
#' }
#'
#' @importFrom lubridate ymd_hms ymd
#' @importFrom dplyr mutate group_by summarise rename select filter %>%
#' @importFrom tibble as_tibble tibble
#' @importFrom stats cor.test cor median sd p.adjust quantile
#' @importFrom boot boot
#' @export
mov.cor.dm <- function(df,
                       Clim,
                       TreeNum,
                       win_size,
                       cor_method = c("pearson", "kendall", "spearman"),
                       boot = FALSE,
                       R = 1000,
                       boot.ci = 0.05,
                       set_seed = 1,
                       dm_stat = c("mean", "min", "max", "median", "amplitude", "change"),
                       clim_vars = NULL,
                       lag_days = 0,
                       accum_days = 1,
                       clim_fun = "raw",
                       min_complete = NULL,
                       p_adjust_method = "BH") {
  TIME <- DATE <- dm <- TIME <- NULL

  cor_method <- match.arg(cor_method)
  dm_stat <- match.arg(dm_stat)

  allowed_clim_fun <- c("raw", "mean", "sum", "max", "min", "median")

  if (!is.numeric(TreeNum) || length(TreeNum) != 1 || is.na(TreeNum)) {
    stop("'TreeNum' must be a single numeric value.")
  }
  if (TreeNum > ncol(df) - 1) {
    stop("The tree number is not valid for the provided dendrometer dataset.")
  }

  if (!is.numeric(win_size) || length(win_size) != 1 || is.na(win_size) || win_size < 18) {
    stop("'win_size' must be a single numeric value >= 18.")
  }

  if (!is.logical(boot) || length(boot) != 1) {
    stop("'boot' must be TRUE or FALSE.")
  }

  if (!boot.ci %in% c(0.01, 0.05, 0.1)) {
    stop("'boot.ci' must be one of 0.01, 0.05, or 0.1.")
  }

  if (!is.numeric(R) || length(R) != 1 || is.na(R) || R <= 0) {
    stop("'R' must be a positive number.")
  }

  if (is.null(min_complete)) {
    min_complete <- max(3L, ceiling(0.8 * win_size))
  } else {
    if (!is.numeric(min_complete) || length(min_complete) != 1 || is.na(min_complete) || min_complete < 3) {
      stop("'min_complete' must be a single numeric value >= 3.")
    }
    min_complete <- as.integer(min_complete)
  }

  if (!p_adjust_method %in% c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")) {
    stop("'p_adjust_method' must be one of the supported p.adjust methods or 'none'.")
  }

  parse_datetime <- function(x) {
    if (inherits(x, "POSIXct")) return(as.POSIXct(x))
    if (inherits(x, "Date")) return(as.POSIXct(x))

    out <- suppressWarnings(lubridate::ymd_hms(x, quiet = TRUE))
    if (all(is.na(out))) {
      out_date <- suppressWarnings(lubridate::ymd(x, quiet = TRUE))
      if (!all(is.na(out_date))) out <- as.POSIXct(out_date)
    }
    if (all(is.na(out))) out <- suppressWarnings(as.POSIXct(x))
    out
  }

  parse_date_only <- function(x) {
    if (inherits(x, "Date")) return(as.Date(x))
    if (inherits(x, "POSIXct")) return(as.Date(x))

    out <- suppressWarnings(lubridate::ymd(x, quiet = TRUE))
    if (all(is.na(out))) {
      out_dt <- suppressWarnings(lubridate::ymd_hms(x, quiet = TRUE))
      if (!all(is.na(out_dt))) out <- as.Date(out_dt)
    }
    if (all(is.na(out))) out <- suppressWarnings(as.Date(x))
    out
  }

  resolve_per_var_arg <- function(arg, var_names, arg_name, allowed = NULL) {
    if (length(arg) == 1 && is.null(names(arg))) {
      out <- rep(arg, length(var_names))
      names(out) <- var_names
    } else if (!is.null(names(arg))) {
      miss <- setdiff(var_names, names(arg))
      if (length(miss) > 0) {
        stop(sprintf(
          "The following selected climate variables are missing in '%s': %s",
          arg_name, paste(miss, collapse = ", ")
        ))
      }
      out <- arg[var_names]
    } else if (length(arg) == length(var_names)) {
      out <- arg
      names(out) <- var_names
    } else {
      stop(sprintf(
        "'%s' must be either length 1, length equal to 'clim_vars', or a named vector.",
        arg_name
      ))
    }

    if (!is.null(allowed)) {
      bad <- setdiff(unique(as.character(out)), allowed)
      if (length(bad) > 0) {
        stop(sprintf(
          "Unsupported value(s) in '%s': %s",
          arg_name, paste(bad, collapse = ", ")
        ))
      }
    }

    out
  }

  transform_clim_series <- function(x, lag_days, accum_days, clim_fun) {
    n <- length(x)
    out <- rep(NA_real_, n)

    for (i in seq_len(n)) {
      end_idx <- i - lag_days

      if (clim_fun == "raw") {
        if (end_idx >= 1 && end_idx <= n) out[i] <- x[end_idx]
      } else {
        start_idx <- end_idx - accum_days + 1
        if (start_idx >= 1 && end_idx <= n) {
          ww <- x[start_idx:end_idx]
          if (!all(is.na(ww))) {
            out[i] <- switch(
              clim_fun,
              mean = mean(ww, na.rm = TRUE),
              sum = sum(ww, na.rm = TRUE),
              max = max(ww, na.rm = TRUE),
              min = min(ww, na.rm = TRUE),
              median = stats::median(ww, na.rm = TRUE)
            )
          }
        }
      }
    }

    out
  }

  summarize_daily_dm <- function(dm_df, stat) {
    daily <- dm_df %>%
      dplyr::mutate(DATE = as.Date(TIME)) %>%
      dplyr::group_by(DATE) %>%
      dplyr::summarise(
        dm_mean = mean(dm, na.rm = TRUE),
        dm_min = min(dm, na.rm = TRUE),
        dm_max = max(dm, na.rm = TRUE),
        dm_median = stats::median(dm, na.rm = TRUE),
        .groups = "drop"
      )

    daily$dm_amplitude <- daily$dm_max - daily$dm_min
    daily$dm_change <- c(NA_real_, diff(daily$dm_mean))

    daily_value <- switch(
      stat,
      mean = daily$dm_mean,
      min = daily$dm_min,
      max = daily$dm_max,
      median = daily$dm_median,
      amplitude = daily$dm_amplitude,
      change = daily$dm_change
    )

    tibble::tibble(TIME = as.Date(daily$DATE), dm_value = daily_value)
  }

  # dendrometer
  df[[1]] <- parse_datetime(df[[1]])
  if (all(is.na(df[[1]]))) {
    stop("Could not parse the first column of 'df' as date-time.")
  }

  dm_name <- names(df)[TreeNum + 1]

  dm_df <- tibble::as_tibble(df) %>%
    dplyr::select(1, TreeNum + 1) %>%
    dplyr::rename(
      TIME = 1,
      dm = 2
    )

  dm_daily <- summarize_daily_dm(dm_df, dm_stat)

  # climate
  Clim[[1]] <- parse_date_only(Clim[[1]])
  if (all(is.na(Clim[[1]]))) {
    stop("Could not parse the first column of 'Clim' as daily date.")
  }

  clim_tbl <- tibble::as_tibble(Clim) %>%
    dplyr::rename(TIME = 1) %>%
    dplyr::mutate(TIME = as.Date(TIME)) %>%
    dplyr::group_by(TIME) %>%
    dplyr::summarise(
      dplyr::across(
        .cols = where(is.numeric),
        .fns = ~ mean(.x, na.rm = TRUE)
      ),
      .groups = "drop"
    )

  all_clim_names <- names(clim_tbl)[-1]

  if (length(all_clim_names) == 0) {
    stop("'Clim' must contain at least one numeric climate variable after the date column.")
  }

  if (is.null(clim_vars)) {
    clim_vars <- all_clim_names
  } else {
    if (!is.character(clim_vars)) {
      stop("'clim_vars' must be a character vector.")
    }
    miss <- setdiff(clim_vars, all_clim_names)
    if (length(miss) > 0) {
      stop(sprintf(
        "The following climate variables were not found in 'Clim': %s",
        paste(miss, collapse = ", ")
      ))
    }
  }

  clim_fun_map <- resolve_per_var_arg(clim_fun, clim_vars, "clim_fun", allowed = allowed_clim_fun)
  lag_map <- resolve_per_var_arg(lag_days, clim_vars, "lag_days")
  accum_map <- resolve_per_var_arg(accum_days, clim_vars, "accum_days")

  lag_map <- as.integer(unname(lag_map))
  names(lag_map) <- clim_vars

  accum_map <- as.integer(unname(accum_map))
  names(accum_map) <- clim_vars

  if (any(is.na(lag_map)) || any(lag_map < 0)) {
    stop("All values in 'lag_days' must be non-negative integers.")
  }
  if (any(is.na(accum_map)) || any(accum_map < 1)) {
    stop("All values in 'accum_days' must be positive integers.")
  }

  raw_bad <- clim_vars[clim_fun_map == "raw" & accum_map != 1]
  if (length(raw_bad) > 0) {
    warning(
      "For variables using clim_fun = 'raw', accum_days is ignored. Affected variables: ",
      paste(raw_bad, collapse = ", ")
    )
  }

  clim_trans <- tibble::tibble(TIME = clim_tbl$TIME)
  for (nm in clim_vars) {
    clim_trans[[nm]] <- transform_clim_series(
      x = clim_tbl[[nm]],
      lag_days = unname(lag_map[nm]),
      accum_days = unname(accum_map[nm]),
      clim_fun = unname(clim_fun_map[nm])
    )
  }

  common_dates <- intersect(dm_daily$TIME, clim_trans$TIME)

  if (length(common_dates) == 0) {
    stop("'df' and 'Clim' must have time overlap after daily aggregation.")
  }

  if ((length(common_dates) != nrow(dm_daily)) || (length(common_dates) != nrow(clim_trans))) {
    warning("The temporal coverage of 'df' and 'Clim' is not the same. Output will be shortened to the common time period.")
  }

  dm_daily <- dm_daily %>% dplyr::filter(TIME %in% common_dates)
  clim_trans <- clim_trans %>% dplyr::filter(TIME %in% common_dates)

  dm_daily <- dm_daily[order(dm_daily$TIME), , drop = FALSE]
  clim_trans <- clim_trans[order(clim_trans$TIME), , drop = FALSE]

  n_row <- nrow(clim_trans)

  wl <- ifelse(win_size %% 2 == 0, win_size + 1, win_size)
  wl <- as.integer(wl)
  half_w <- (wl - 1) / 2

  if (n_row < wl) {
    stop("Not enough overlapping daily observations for the requested 'win_size'.")
  }

  if (boot) set.seed(set_seed)

  result_list <- lapply(clim_vars, function(cln) {
    cls <- clim_trans[[cln]]

    if (boot) {
      estimate <- data.frame(
        TIME = as.Date(dm_daily$TIME),
        corr = NA_real_,
        corr_lw_lim = NA_real_,
        corr_up_lim = NA_real_,
        n_complete = NA_real_
      )

      for (end_idx in seq.int(wl, n_row)) {
        start_idx <- end_idx - wl + 1
        mid_idx <- start_idx + half_w

        val_1 <- cls[start_idx:end_idx]
        val_2 <- dm_daily$dm_value[start_idx:end_idx]

        ok <- is.finite(val_1) & is.finite(val_2)
        estimate$n_complete[mid_idx] <- sum(ok)

        if (sum(ok) >= min_complete) {
          temp_data <- data.frame(x = val_1[ok], y = val_2[ok])

          b.corr <- suppressWarnings(
            boot::boot(
              temp_data,
              statistic = function(temp_data, i) {
                stats::cor(temp_data[i, "x"], temp_data[i, "y"], method = cor_method)
              },
              R = R
            )
          )

          estimate$corr[mid_idx] <- b.corr$t0

          cis <- c(0.01, 0.05, 0.1)
          loc.cis <- which(cis == boot.ci)
          cis_lw <- c(0.005, 0.025, 0.05)
          cis_up <- c(0.995, 0.975, 0.95)
          cis_app <- c(cis_lw[loc.cis], cis_up[loc.cis])

          qv <- stats::quantile(as.vector(b.corr$t), cis_app, na.rm = TRUE)
          estimate$corr_lw_lim[mid_idx] <- qv[1]
          estimate$corr_up_lim[mid_idx] <- qv[2]
        }
      }

      estimate <- tibble::as_tibble(estimate)
      estimate$significance <- with(estimate, ifelse(corr_lw_lim * corr_up_lim > 0, TRUE, FALSE))
      estimate$clim_var <- cln
      estimate$tree <- dm_name
      estimate

    } else {
      estimate <- data.frame(
        TIME = as.Date(dm_daily$TIME),
        corr = NA_real_,
        p_val = NA_real_,
        n_complete = NA_real_
      )

      for (end_idx in seq.int(wl, n_row)) {
        start_idx <- end_idx - wl + 1
        mid_idx <- start_idx + half_w

        val_1 <- cls[start_idx:end_idx]
        val_2 <- dm_daily$dm_value[start_idx:end_idx]

        ok <- is.finite(val_1) & is.finite(val_2)
        estimate$n_complete[mid_idx] <- sum(ok)

        if (sum(ok) >= min_complete) {
          a <- suppressWarnings(
            stats::cor.test(val_1[ok], val_2[ok], method = cor_method)
          )
          estimate$corr[mid_idx] <- unname(a$estimate)
          estimate$p_val[mid_idx] <- a$p.value
        }
      }

      estimate <- tibble::as_tibble(estimate)
      estimate$p_adj <- NA_real_
      okp <- !is.na(estimate$p_val)

      if (any(okp)) {
        if (p_adjust_method == "none") {
          estimate$p_adj[okp] <- estimate$p_val[okp]
        } else {
          estimate$p_adj[okp] <- stats::p.adjust(estimate$p_val[okp], method = p_adjust_method)
        }
      }

      estimate$significance <- ifelse(!is.na(estimate$p_adj) & estimate$p_adj < 0.05, TRUE, FALSE)
      estimate$clim_var <- cln
      estimate$tree <- dm_name
      estimate
    }
  })

  names(result_list) <- clim_vars

  climate_settings <- tibble::tibble(
    clim_var = clim_vars,
    clim_fun = unname(clim_fun_map[clim_vars]),
    lag_days = unname(lag_map[clim_vars]),
    accum_days = unname(accum_map[clim_vars])
  )

  out <- list(
    results = result_list,
    metadata = list(
      tree_num = TreeNum,
      tree_name = dm_name,
      clim_names = clim_vars,
      climate_settings = climate_settings,
      win_size = wl,
      cor_method = cor_method,
      boot = boot,
      R = if (boot) R else NA_integer_,
      boot_ci = if (boot) boot.ci else NA_real_,
      set_seed = if (boot) set_seed else NA_integer_,
      dm_stat = dm_stat,
      min_complete = min_complete,
      p_adjust_method = if (!boot) p_adjust_method else NA_character_,
      n_obs = n_row
    ),
    call = match.call()
  )

  class(out) <- unique(c(
    if (boot) "mov_cor_dm_boot",
    "mov_cor_dm",
    if (boot) "mov.cor.boot" else "mov.cor"
  ))

  out
}


#' @title Print method for running dendrometer-climate correlation objects
#'
#' @param x Object of class \code{"mov_cor_dm"}.
#' @param ... Unused.
#'
#' @return The object, invisibly.
#' @export
print.mov_cor_dm <- function(x, ...) {
  cat("Running dendrometer-climate correlation\n")
  cat("Tree: ", x$metadata$tree_name, "\n", sep = "")
  cat("Window size: ", x$metadata$win_size, " days\n", sep = "")
  cat("Correlation method: ", x$metadata$cor_method, "\n", sep = "")
  cat("Dendrometer statistic: ", x$metadata$dm_stat, "\n", sep = "")
  cat("Bootstrap: ", x$metadata$boot, "\n", sep = "")
  cat("Climate variables: ", paste(x$metadata$clim_names, collapse = ", "), "\n\n", sep = "")

  cat("Climate settings\n")
  print(x$metadata$climate_settings)

  invisible(x)
}


#' @title Summary method for running dendrometer-climate correlation objects
#'
#' @description
#' Summarizes running dendrometer-climate correlations for both bootstrapped and
#' non-bootstrapped outputs, and appends climate-specific settings to the
#' summary tables.
#'
#' @param object Object of class \code{"mov_cor_dm"}.
#' @param absolute Logical. If \code{TRUE} (default), the strongest correlation is
#'   identified by absolute magnitude.
#' @param top_n Integer. Number of top windows to return per climate variable.
#' @param ... Unused.
#'
#' @return A list of class \code{"summary_mov_cor_dm"} containing:
#' \itemize{
#'   \item \code{table}: summary table for each climate variable
#'   \item \code{top_windows}: top running windows per climate variable
#'   \item \code{metadata}: metadata from the original object
#'   \item \code{call}: original function call
#' }
#'
#' @export
summary.mov_cor_dm <- function(object, absolute = TRUE, top_n = 5, ...) {

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

  if (!is.numeric(top_n) || length(top_n) != 1 || is.na(top_n) || top_n < 1) {
    stop("'top_n' must be a single positive number.")
  }
  top_n <- as.integer(top_n)

  boot_flag <- isTRUE(object$metadata$boot)
  climate_settings <- object$metadata$climate_settings

  summarize_one <- function(tab, boot = FALSE, absolute = TRUE, top_n = 5) {
    tab <- tibble::as_tibble(tab)

    ok <- !is.na(tab$corr)
    if (!any(ok)) {
      sm <- tibble::tibble(
        clim_var = unique(tab$clim_var)[1],
        tree = unique(tab$tree)[1],
        n_windows = 0L,
        n_significant = 0L,
        prop_significant = NA_real_,
        peak_time = as.Date(NA),
        peak_corr = NA_real_,
        mean_corr = NA_real_,
        median_corr = NA_real_,
        min_corr = NA_real_,
        max_corr = NA_real_
      )

      if (boot) {
        sm$peak_lw <- NA_real_
        sm$peak_up <- NA_real_
      } else {
        sm$peak_p <- NA_real_
        sm$peak_p_adj <- NA_real_
      }

      return(list(summary = sm, top = tab[0, , drop = FALSE]))
    }

    tab_ok <- tab[ok, , drop = FALSE]

    peak_idx <- if (absolute) {
      which.max(abs(tab_ok$corr))
    } else {
      which.max(tab_ok$corr)
    }

    sm <- tibble::tibble(
      clim_var = unique(tab_ok$clim_var)[1],
      tree = unique(tab_ok$tree)[1],
      n_windows = nrow(tab_ok),
      n_significant = sum(tab_ok$significance, na.rm = TRUE),
      prop_significant = mean(tab_ok$significance, na.rm = TRUE),
      peak_time = tab_ok$TIME[peak_idx],
      peak_corr = tab_ok$corr[peak_idx],
      mean_corr = mean(tab_ok$corr, na.rm = TRUE),
      median_corr = stats::median(tab_ok$corr, na.rm = TRUE),
      min_corr = min(tab_ok$corr, na.rm = TRUE),
      max_corr = max(tab_ok$corr, na.rm = TRUE)
    )

    if (boot) {
      sm$peak_lw <- tab_ok$corr_lw_lim[peak_idx]
      sm$peak_up <- tab_ok$corr_up_lim[peak_idx]

      ord_idx <- order(!tab_ok$significance, -abs(tab_ok$corr))
      top_tab <- tab_ok[ord_idx, , drop = FALSE]
    } else {
      sm$peak_p <- tab_ok$p_val[peak_idx]
      sm$peak_p_adj <- tab_ok$p_adj[peak_idx]

      p_adj_sort <- ifelse(is.na(tab_ok$p_adj), Inf, tab_ok$p_adj)
      ord_idx <- order(!tab_ok$significance, p_adj_sort, -abs(tab_ok$corr))
      top_tab <- tab_ok[ord_idx, , drop = FALSE]
    }

    top_tab <- utils::head(top_tab, top_n)

    list(summary = sm, top = top_tab)
  }

  res_list <- lapply(
    object$results,
    summarize_one,
    boot = boot_flag,
    absolute = absolute,
    top_n = top_n
  )

  summary_table <- dplyr::bind_rows(lapply(res_list, `[[`, "summary"))
  top_windows <- dplyr::bind_rows(lapply(res_list, `[[`, "top"))

  if (!is.null(climate_settings) && nrow(climate_settings) > 0) {
    idx_sum <- match(summary_table$clim_var, climate_settings$clim_var)
    summary_table$clim_fun <- climate_settings$clim_fun[idx_sum]
    summary_table$lag_days <- climate_settings$lag_days[idx_sum]
    summary_table$accum_days <- climate_settings$accum_days[idx_sum]

    idx_top <- match(top_windows$clim_var, climate_settings$clim_var)
    top_windows$clim_fun <- climate_settings$clim_fun[idx_top]
    top_windows$lag_days <- climate_settings$lag_days[idx_top]
    top_windows$accum_days <- climate_settings$accum_days[idx_top]
  }

  out <- list(
    table = summary_table,
    top_windows = top_windows,
    metadata = object$metadata,
    call = object$call
  )

  class(out) <- "summary_mov_cor_dm"
  out
}


#' @title Print method for summaries of running dendrometer-climate correlation objects
#'
#' @param x Object of class \code{"summary_mov_cor_dm"}.
#' @param ... Unused.
#'
#' @return The summary object, invisibly.
#' @export
print.summary_mov_cor_dm <- function(x, ...) {
  cat("Summary of running dendrometer-climate correlations\n")
  cat("Tree: ", x$metadata$tree_name, "\n", sep = "")
  cat("Window size: ", x$metadata$win_size, " days\n", sep = "")
  cat("Method: ", x$metadata$cor_method, "\n", sep = "")
  cat("Dendrometer statistic: ", x$metadata$dm_stat, "\n", sep = "")
  cat("Bootstrap: ", x$metadata$boot, "\n\n", sep = "")

  cat("Overall summary\n")
  print(x$table)

  cat("\nTop windows\n")
  print(x$top_windows)

  invisible(x)
}