R/broom-methods.R

In ggchangepoint: Combines Changepoint Analysis with 'ggplot2'

#' @importFrom generics tidy
#' @export
generics::tidy

#' @importFrom generics glance
#' @export
generics::glance

#' @importFrom generics augment
#' @export
generics::augment

#' @importFrom ggplot2 autoplot
#' @export
ggplot2::autoplot

#' @importFrom tibble as_tibble
#' @export
tibble::as_tibble

#' Tidy a ggcpt object
#'
#' Returns the changepoints tibble (one row per changepoint), or the
#' coefficients of every segment.
#'
#' @param x A \code{ggcpt} object.
#' @param what \code{"changepoints"} (the default) or
#'   \code{"coefficients"}: one row per segment and model term, with
#'   \code{estimate}, \code{std_error} and a Wald interval. For a formula
#'   fit the terms are the regression's; for a series they are the segment
#'   level (and its slope in time, for a change in slope), so a mean-shift
#'   result reports each segment's mean with its standard error. Engines
#'   that estimate the coefficients themselves (\pkg{segmented}'s slopes,
#'   \pkg{fastcpd}'s per-segment parameters) report their own estimates.
#' @param conf_level Confidence level for the coefficient intervals.
#' @param ... Additional arguments (ignored).
#' @return A tibble. For \code{"changepoints"}: columns \code{cp},
#'   \code{cp_value}, and any method-specific columns. For
#'   \code{"coefficients"}: \code{segment}, \code{start}, \code{end},
#'   \code{term}, \code{estimate}, \code{std_error}, \code{conf_low},
#'   \code{conf_high} (plus \code{coordinate} for a multivariate result).
#' @export
#' @examples
#' set.seed(1)
#' fit <- cpt_detect(c(rnorm(80), rnorm(80, 3)), method = "pelt")
#' tidy(fit)
#' tidy(fit, "coefficients")
tidy.ggcpt <- function(x, what = c("changepoints", "coefficients"),
                       conf_level = 0.95, ...) {
  what <- cpt_match_arg(what)
  if (what == "coefficients") {
    validate_scalar(conf_level, "conf_level", min = 0, max = 1,
                    min_open = TRUE, max_open = TRUE)
    return(tidy_coefficients(x, conf_level))
  }
  x$changepoints
}

#' Tidy the significance regions of a ggcpt object
#'
#' The interval-valued methods (currently \code{\link{nsp_wrapper}()})
#' return regions rather than points: each is guaranteed to contain at least
#' one changepoint at a prescribed \emph{global} significance level. A
#' region is a different object from a confidence interval around an
#' estimate, so it lives in its own slot and has its own accessor rather
#' than being folded into \code{tidy()}.
#'
#' @param x A \code{ggcpt} object.
#' @return A tibble with columns \code{start}, \code{end} (positions),
#'   \code{length}, and, when the result carries a time index,
#'   \code{start_index}/\code{end_index} on the original scale. Any further
#'   columns the engine supplied are carried through after those:
#'   \code{\link{nsp_wrapper}()} adds \code{value}, the region's statistic.
#'   A zero-row tibble when the result carries no regions.
#' @seealso \code{\link{nsp_wrapper}()}, \code{\link{geom_cpt_region}()},
#'   \code{\link{cpt_confint}()}.
#' @export
#' @examples
#' set.seed(2026)
#' fit <- as_ggcpt(50, c(rnorm(50), rnorm(50, 4)),
#'                 regions = data.frame(start = 45, end = 56))
#' cpt_regions(fit)
cpt_regions <- function(x) {
  if (!is_ggcpt(x)) {
    cpt_abort("`x` must be a ggcpt object.", class = "bad_argument")
  }
  reg <- x$regions
  if (is.null(reg) || nrow(reg) == 0) {
    # The same columns a non-empty return would have, so
    # rbind(cpt_regions(a), cpt_regions(b)) works when one of the two is
    # empty and the other is indexed. empty_confint() already does this,
    # and ?alarms advertises the property.
    out <- tibble::tibble(start = integer(), end = integer(),
                          length = integer())
    # `[[..., exact = TRUE]]`, not `$`: a tibble warns "Unknown or
    # uninitialised column" when asked for a column it does not have.
    if ("index_value" %in% names(x$data)) {
      iv <- x$data[["index_value", exact = TRUE]]
      out$start_index <- iv[0]
      out$end_index <- iv[0]
    }
    return(out)
  }
  out <- tibble::as_tibble(reg)
  out$length <- out$end - out$start + 1L
  front <- intersect(c("start", "end", "length", "start_index", "end_index"),
                     names(out))
  out[, c(front, setdiff(names(out), front)), drop = FALSE]
}

#' Glance at a ggcpt object
#'
#' Returns a one-row summary of a changepoint detection result.
#'
#' @param x A \code{ggcpt} object.
#' @param ... Additional arguments (ignored).
#' @return A one-row tibble with columns: \code{n}, \code{n_changepoints},
#'   \code{method}, \code{change_in}, \code{penalty_type}, \code{penalty_value},
#'   \code{cp_convention}, \code{total_cost} (\code{NA} when the engine does
#'   not expose a cost), \code{runtime} (elapsed seconds when measured by
#'   \code{cpt_detect()}, otherwise \code{NA}) and \code{engine_version}
#'   (the version of the engine package that produced the result, recorded
#'   when it was made, or \code{NA} for a result with no engine package,
#'   such as one built by \code{\link{as_ggcpt}()}; see
#'   \code{\link{cpt_verify}()}) and \code{family} (the distribution family
#'   the fit was asked for with \code{cpt_detect(family = )}, or \code{NA}
#'   for the method's own default).
#'
#' @details
#' \code{total_cost} is reported on whatever scale the engine itself uses, so
#' it is meaningful when comparing penalties within one method and not when
#' comparing one method against another. For the \pkg{changepoint} engines it
#' is the unpenalised \eqn{-2\log L} of the chosen segmentation. Four cases
#' there are \code{NA} rather than filled with a number that would not mean
#' the same thing:
#' \itemize{
#'   \item \code{"binseg"} and \code{"segneigh"}, whose \code{cpt.range}
#'     fits report the raw within-segment cost instead: for one and the same
#'     segmentation that is 219.7 where a PELT fit reports 659.9;
#'   \item \code{"np"}, because \pkg{changepoint.np} defines no
#'     \code{logLik} method;
#'   \item a change in \emph{mean} under the default \code{"MBIC"} penalty.
#'     Loading \pkg{changepoint.np} (which this package imports, so it is
#'     always loaded) replaces \pkg{changepoint}'s \code{logLik} method for
#'     \code{cpt} objects with one that errors on exactly that combination.
#'     Any other penalty (\code{"BIC"}, \code{"AIC"}, a numeric value)
#'     reports normally, as do \code{change_in = "var"} and
#'     \code{"meanvar"}.
#' }
#' @export
#' @examples
#' set.seed(2026)
#' fit <- cpt_detect(c(rnorm(60), rnorm(60, 3)), method = "pelt")
#' glance(fit)
glance.ggcpt <- function(x, ...) {
  total_cost <- NA_real_
  if (!is.null(x$fit)) {
    if (inherits(x$fit, "cpt") && !inherits(x$fit, "cpt.range") &&
        identical(cpt_test_stat(x$fit), "Normal")) {
      # changepoint's logLik() returns c(`-2*logLik`, `-2*logLik + pen`);
      # report the unpenalised first element, which is already a cost (do not
      # negate it). Verified against changepoint's own output: under a BIC
      # penalty the pair is (544.4, 555.0), i.e. the second element is the
      # first plus the penalty, and `change_in = "var"` and `"meanvar"` both
      # report on the same scale.
      #
      # The two guards above keep values that are NOT on that scale out of
      # the column. A `cpt.range` fit (BinSeg, SegNeigh) returns the raw
      # within-segment cost instead -- for one and the same segmentation,
      # cpt at 100, BinSeg reported 178.1 where PELT reported 544.4 -- and
      # it is that path, not the test statistic, which warns "Not changed
      # to be -2*logLik". A non-Normal test statistic does not warn at all:
      # CUSUM and CSS *error* with "logLik is only valid for distributional
      # assumptions", which is why they are excluded by name rather than
      # relied on to warn. cpt.np() fits are skipped for a third reason:
      # they carry the `cpt` class but describe no distributional
      # changepoint type, so logLik() rejects them with "Unknown
      # changepoint type".
      #
      # `logLik` is deliberately unqualified: the method for `cpt` is an S4
      # method owned by changepoint, which this package @imports, so plain
      # dispatch finds it. `stats::logLik` reaches only the S3 generic and
      # never dispatches to it. The tryCatch is not belt-and-braces: loading
      # changepoint.np replaces that method with one that errors on a
      # change in mean under the MBIC penalty (see the details section).
      total_cost <- tryCatch(
        {
          utils::capture.output(ll <- logLik(x$fit))
          as.numeric(ll)[1]
        },
        error = function(e) NA_real_
      )
    } else if (is.list(x$fit) && !is.data.frame(x$fit)) {
      # Exact [[ ]] subsetting: $ would partial-match unrelated elements
      # (e.g. DeCAFS's costFunction). And NOT a data frame: one IS a list,
      # and fchange_run() stores a frame directly as $fit, so a frame that
      # ever acquired a `value` column would make `cand` the whole COLUMN
      # -- whose last element the reconciliation below would then report as
      # the model's total cost. (The length > 1 case is kept for fpop,
      # which genuinely exposes a per-position cost vector whose terminal
      # element IS the total; a data frame column is not that.)
      cand <- x$fit[["cost"]] %||% x$fit[["loss"]] %||% x$fit[["value"]]
      if (!is.null(cand) && is.numeric(cand)) total_cost <- cand
    }
    # Some engines (e.g. fpop) expose a per-position cost vector; glance is
    # one row per model, so keep the terminal (total) cost only.
    total_cost <- suppressWarnings(as.numeric(total_cost))
    if (length(total_cost) == 0) {
      total_cost <- NA_real_
    } else if (length(total_cost) > 1) {
      total_cost <- total_cost[length(total_cost)]
    }
  }

  penalty_value <- if (is.list(x$penalty)) x$penalty$value else NA_real_
  penalty_value <- suppressWarnings(as.numeric(penalty_value %||% NA_real_))
  if (length(penalty_value) != 1) penalty_value <- NA_real_

  runtime <- x$runtime %||% NA_real_

  # glance() promises exactly one row. A hand-built object may carry a
  # zero-length `method` or `change_in` (that used to be new_ggcpt()'s
  # default), and tibble() would then recycle every other column down to
  # zero rows -- an empty summary rather than a summary saying "unknown".
  tibble::tibble(
    n = nrow(x$data),
    n_changepoints = nrow(x$changepoints),
    method = scalar_chr(x$method),
    change_in = scalar_chr(x$change_in),
    penalty_type = if (is.list(x$penalty)) {
      scalar_chr(x$penalty$type)
    } else {
      NA_character_
    },
    penalty_value = penalty_value,
    cp_convention = x$cp_convention %||% "left",
    total_cost = total_cost,
    runtime = runtime,
    engine_version = scalar_chr(x$versions$engine_version %||% NA_character_),
    family = scalar_chr(x$family %||% NA_character_)
  )
}

#' Augment a ggcpt object
#'
#' Returns the original data with added columns: \code{seg_id}, \code{.fitted},
#' \code{.resid}, and \code{is_changepoint}.
#'
#' @param x A \code{ggcpt} object.
#' @param ... Additional arguments (ignored).
#' @return A tibble with one row per observation: the data as the result
#'   carries it, plus four added columns. The data half depends on the
#'   result: \code{index} and \code{value} for a univariate one,
#'   \code{index} plus \strong{one column per coordinate} (named as the
#'   input's columns were) for a multivariate one, and an extra
#'   \code{fitted} column for the engines that supply their own fitted
#'   signal. The added four are always the same:
#'   \describe{
#'     \item{\code{seg_id}}{which segment the observation falls in,
#'       counting from 1.}
#'     \item{\code{.fitted}}{the segment's \code{param_estimate}: the
#'       segment \strong{mean}, for every method in the package, or the
#'       engine's own fitted signal where there is one. See the details
#'       below for the multivariate case.}
#'     \item{\code{.resid}}{\code{value - .fitted}, against the
#'       univariate series the result carries.}
#'     \item{\code{is_changepoint}}{\code{TRUE} at each detected
#'       location, under the result's \code{cp_convention}.}
#'   }
#'   Measured on a \code{pelt} fit the columns are \code{index},
#'   \code{value}, \code{seg_id}, \code{.fitted}, \code{.resid},
#'   \code{is_changepoint}.
#'
#' @details
#' For a multivariate result every coordinate is returned, but the
#' changepoints are shared across them, so \code{seg_id} and
#' \code{is_changepoint} apply to the whole row. Each coordinate also gets
#' its own \code{.fitted_<name>} and \code{.resid_<name>} (its segment
#' means and residuals), while \code{.fitted} and
#' \code{.resid} describe the \strong{univariate series the result
#' carries}: \code{$data$value}, the same series
#' \code{$segments$param_estimate} summarises, so \code{.resid} is always
#' \code{value - .fitted}. For most multivariate engines that series is
#' the first coordinate; \code{fmean}, \code{fcov}, \code{kwc} and
#' \code{fabisearch} store the cross-sectional mean \code{rowMeans()}
#' instead, and for those \code{.fitted}/\code{.resid} describe that mean
#' rather than any one column. Either way the two columns agree with each
#' other, which is what makes \code{.resid} a residual. When an engine
#' supplies its own fitted signal that signal is used for \code{.fitted} in
#' place of the segment means, and rides along in a \code{fitted} column of
#' its own, so for those engines the two columns agree. The engines that
#' do this are exactly the ones \code{\link{cpt_methods}()} marks in its
#' \code{fitted} column: \code{smuce}, \code{hsmuce}, \code{cpop},
#' \code{bcp}, \code{bocpd}, \code{beast}, \code{decafs},
#' \code{segmented}, \code{envcpt}, \code{mcp} and \code{bfast}; a
#' formula fit carries its piecewise regression's fitted values the same
#' way.
#' @export
#' @examples
#' set.seed(2026)
#' fit <- cpt_detect(c(rnorm(60), rnorm(60, 3)), method = "pelt")
#' head(augment(fit))
augment.ggcpt <- function(x, ...) {
  # For a multivariate result use the wide frame (index + one column per
  # coordinate) so all coordinates are kept; otherwise use the univariate
  # single-column frame.
  use_wide <- !is.null(x$data_wide)
  data <- if (use_wide) tibble::as_tibble(x$data_wide) else x$data

  if (!use_wide) {
    # Rename columns by position in a way that works with >2 columns
    nms <- names(data)
    if (length(nms) >= 2) {
      names(data)[1:2] <- c("index", "value")
    }
  }

  # The index used to flag changepoints, and the value vector used for
  # `.resid`.
  index_col <- data[["index"]]
  # `.resid` has to be computed against the series `param_estimate` was
  # computed FROM, which is always `x$data$value`. Reading coordinate one
  # out of `data_wide` instead was right for twelve of the sixteen
  # multivariate wrappers and wrong for four: `fmean`, `fcov`, `kwc` and
  # `fabisearch` store `rowMeans(X)` as their univariate series (see
  # wrap-functional.R and wrap-highdim2.R), so build_segments() derived
  # `param_estimate` from the row averages while this subtracted it from
  # coordinate one -- two different quantities, and `.resid` was not a
  # residual at all.
  #
  # Measured on a 120x6 `fmean` fit: `.resid` equalled `X[, 1] - .fitted`
  # and did NOT equal `value - .fitted`. Using `data$value` is identical
  # for the twelve wrappers where the two agree, so this changes nothing
  # for them and makes the other four correct.
  #
  # Taken from `x$data`, not from the local `data`: on the multivariate path
  # that has been replaced by the wide frame, which has no `value` column at
  # all -- reading it there gave a zero-length vector and `$<-` refused to
  # recycle it to n rows.
  # Two conventions live in the next twenty lines and they agree only
  # because `$data$index` is `seq_len(n)`: `.fitted`/`seg_id` are written by
  # ROW POSITION (`idx <- seq(s$start, s$end)`, and segment bounds are
  # positions), while `is_changepoint` is written by VALUE MATCH against
  # `$changepoints$cp`. attach_index() keeps `index` positional and adds the
  # user's scale as `index_value` precisely so that stays true -- if either
  # ever became the user's scale, one of the two would silently be wrong.
  value_vec <- as.numeric(x$data$value)

  data$seg_id <- NA_integer_
  data$.fitted <- NA_real_
  data$.resid <- NA_real_
  data$is_changepoint <- FALSE

  if (nrow(x$segments) > 0) {
    for (i in seq_len(nrow(x$segments))) {
      s <- x$segments[i, ]
      idx <- seq(s$start, s$end)
      data$seg_id[idx] <- s$seg_id
      data$.fitted[idx] <- s$param_estimate
    }
    data$.resid <- value_vec - data$.fitted
  }

  # Prefer the engine's fitted signal when the object carries one.
  engine_fitted <- x$data[["fitted"]]
  if (!is.null(engine_fitted) && length(engine_fitted) == nrow(data)) {
    data$.fitted <- as.numeric(engine_fitted)
    data$.resid <- value_vec - data$.fitted
  }

  # And for a multivariate result, every coordinate's own segment means and
  # residuals: `.fitted`/`.resid` describe only the one univariate series
  # the result carries, which was the whole of the diagnostics a
  # multivariate fit got.
  if (use_wide && nrow(x$segments) > 0) {
    coords <- setdiff(names(x$data_wide), c("index", "index_value"))
    for (cl in coords) {
      v <- as.numeric(x$data_wide[[cl]])
      fit_c <- rep(NA_real_, length(v))
      for (i in seq_len(nrow(x$segments))) {
        idx <- seq(x$segments$start[i], x$segments$end[i])
        fit_c[idx] <- mean(v[idx], na.rm = TRUE)
      }
      data[[paste0(".fitted_", cl)]] <- fit_c
      data[[paste0(".resid_", cl)]] <- v - fit_c
    }
  }

  if (nrow(x$changepoints) > 0) {
    data$is_changepoint[index_col %in% x$changepoints$cp] <- TRUE
  }

  tibble::as_tibble(data)
}

`%||%` <- function(x, y) if (is.null(x)) y else x

# Internal: coerce a metadata field to exactly one string, so a
# zero-length or over-long value cannot change the shape of glance()'s row.
#' @noRd
scalar_chr <- function(x) {
  if (length(x) == 0L) return(NA_character_)
  as.character(x)[1L]
}

# Internal: the test statistic of a changepoint `cpt` fit, or NA when the
# object does not carry one. Only "Normal" fits have a -2*logLik on the scale
# glance() reports.
#' @noRd
cpt_test_stat <- function(fit) {
  tryCatch(as.character(changepoint::test.stat(fit))[1],
           error = function(e) NA_character_)
}

#' Summary of a ggcpt object
#'
#' Provides a human-readable digest of a changepoint detection result,
#' including the segment table with levels and lengths (the level is the
#' segment mean whatever `change_in` says; see \code{\link{new_ggcpt}}),
#' total cost,
#' penalty, and runtime.
#'
#' @param object A \code{ggcpt} object.
#' @param ... Additional arguments (ignored).
#' @return A list with class \code{summary.ggcpt} containing the summary.
#' @exportS3Method base::summary
#' @examples
#' set.seed(2026)
#' summary(cpt_detect(c(rnorm(60), rnorm(60, 3)), method = "pelt"))
summary.ggcpt <- function(object, ...) {
  structure(
    list(
      method = object$method,
      change_in = object$change_in,
      n_changepoints = nrow(object$changepoints),
      cp_convention = object$cp_convention %||% "left",
      n_obs = nrow(object$data),
      penalty = object$penalty,
      segments = object$segments,
      changepoints = object$changepoints,
      runtime = object$runtime %||% NA_real_
    ),
    class = "summary.ggcpt"
  )
}

#' @rdname summary.ggcpt
#' @param x A \code{summary.ggcpt} object (for \code{print()}).
#' @export
print.summary.ggcpt <- function(x, ...) {
  cat("ggcpt Summary\n")
  # Same column as print.ggcpt(): the two are two views of one object, and
  # they used to indent their values differently.
  cat_field("Method", x$method)
  cat_field("Change in", x$change_in)
  cat_field("Changepoints found", x$n_changepoints)
  cat_field("CP convention", x$cp_convention)
  cat_field("Series length", x$n_obs)
  cat_field("Penalty", format_penalty(x$penalty))
  cat_field("Runtime (seconds)", format(x$runtime, digits = 4))
  if (nrow(x$segments) > 0) {
    cat("\nSegments:\n")
    print(x$segments, n = min(nrow(x$segments), 10))
  }
  if (nrow(x$changepoints) > 0) {
    cat("\nChangepoints:\n")
    print(x$changepoints, n = min(nrow(x$changepoints), 10))
  }
  invisible(x)
}

#' Coerce, format, and plot ggcpt objects
#'
#' Convenience S3 methods for working with \code{ggcpt} objects: coerce the
#' changepoints to a tibble or data frame, render a one-line summary string, or
#' produce the default plot (a base-graphics fallback that delegates to
#' \code{\link{autoplot.ggcpt}}).
#'
#' @param x A \code{ggcpt} object.
#' @param row.names,optional Passed to \code{\link[base]{as.data.frame}}.
#' @param .name_repair Ignored (the changepoints tibble already has valid,
#'   unique names); present for signature compatibility with the generic.
#' @param ... Additional arguments passed to methods.
#' @return \code{as_tibble()} and \code{as.data.frame()} return the changepoints
#'   table; \code{format()} returns a length-one character string; \code{plot()}
#'   draws the default plot and returns the \code{ggplot} object invisibly.
#' @seealso \code{\link{ggcpt_plot_methods}} for \code{plot()} on the other
#'   result classes, and \code{\link{tidy.ggcpt}()} for the tidy accessors.
#' @name ggcpt_methods
#' @examples
#' set.seed(2022)
#' res <- cpt_detect(c(rnorm(50), rnorm(50, 5)), method = "pelt")
#' as_tibble(res)
#' as.data.frame(res)
#' format(res)
#' @export
as_tibble.ggcpt <- function(x, ..., .name_repair = NULL) {
  x$changepoints
}

#' @rdname ggcpt_methods
#' @export
as.data.frame.ggcpt <- function(x, row.names = NULL, optional = FALSE, ...) {
  as.data.frame(x$changepoints, row.names = row.names, optional = optional, ...)
}

#' @rdname ggcpt_methods
#' @export
format.ggcpt <- function(x, ...) {
  paste0(
    "ggcpt [", x$method, "] ",
    nrow(x$changepoints), " changepoint(s) on ",
    nrow(x$data), " observations"
  )
}

#' @rdname ggcpt_methods
#' @exportS3Method base::plot
plot.ggcpt <- function(x, ...) {
  plot_via_autoplot(x, ...)
}