---
title: "ggchangepoint: A Unified Tidy Interface for Changepoint Analysis in R"
author: "Youzhi Yu<br><span style='font-size:85%;'>University of Chicago</span>"
bibliography: vignette_reference.bib
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{ggchangepoint: A Unified Tidy Interface for Changepoint Analysis in R}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  fig.width = 8,
  fig.height = 5,
  dpi = 72,
  message = FALSE,
  warning = FALSE,
  comment = "#>",
  fig.alt = "ggchangepoint plot of a time series with its detected changepoints"
)
library(ggchangepoint)
library(ggplot2)
theme_set(theme_light())

# Optional (Suggests) engines: gate the chunks that need them so the
# vignette builds on any installation.
has_stepR       <- requireNamespace("stepR", quietly = TRUE)
has_cpop        <- requireNamespace("cpop", quietly = TRUE)
has_bcp         <- requireNamespace("bcp", quietly = TRUE)
has_ocp         <- requireNamespace("ocp", quietly = TRUE)
has_fpop        <- requireNamespace("fpop", quietly = TRUE)
has_wbs         <- requireNamespace("wbs", quietly = TRUE)
has_not         <- requireNamespace("not", quietly = TRUE)
has_mosum       <- requireNamespace("mosum", quietly = TRUE)
has_idetect     <- requireNamespace("IDetect", quietly = TRUE)
has_breakfast   <- requireNamespace("breakfast", quietly = TRUE)
has_cpm         <- requireNamespace("cpm", quietly = TRUE)
has_decafs      <- requireNamespace("DeCAFS", quietly = TRUE)
has_inspect     <- requireNamespace("InspectChangepoint", quietly = TRUE)
has_strucchange <- requireNamespace("strucchange", quietly = TRUE)
has_segmented   <- requireNamespace("segmented", quietly = TRUE)
has_fastcpd     <- requireNamespace("fastcpd", quietly = TRUE)
has_envcpt      <- requireNamespace("EnvCpt", quietly = TRUE)
```

# Abstract

**ggchangepoint** provides a unified, tidy interface to changepoint detection
across the methodological spectrum. It introduces a single S3 result class,
`ggcpt`, with `broom`-style methods (`tidy()`, `glance()`, `augment()`)
[@robinson2017broom], a central dispatcher `cpt_detect()` covering 50
detection methods across six algorithmic families, and native `ggplot2`
[@wickham2016ggplot2] visualisation through `autoplot()` and a set of
composable geoms. Where a method quantifies its uncertainty (the
simultaneous confidence intervals of SMUCE [@frick2014smuce], the break-date
intervals of Bai-Perron [@bai1998estimating], the posterior distributions of
Bayesian detectors [@barry1993bayesian; @adams2007bocpd]), the result object
carries that uncertainty and the plotting layer can draw it. The package
further supplies a penalty-path diagnostic (CROPS), batch detection over
panels of series, bootstrap stability diagnostics, accuracy metrics aligned
with current benchmarking conventions [@van2020evaluation], ground-truth
simulation, and per-method citations. This article sets out the statistical
background, the design of the package, and each method family in turn, with
worked examples throughout.

# Introduction

Changepoint analysis (locating the instants at which the stochastic
behaviour of an ordered sequence changes) is one of the oldest problems in
statistics, dating back at least to the continuous-inspection schemes of
@page1954continuous, and one of its most active: recent surveys catalogue
dozens of methods [@truong2020selective; @aminikhanghahi2017survey]. Its
applications span virtually every domain that produces sequential data,
including genomics [@picard2005statistical], finance [@athey2022detecting],
climate science [@haslett1989space], and signal processing
[@lavielle2005using].

The R ecosystem mirrors this breadth. Penalised optimal partitioning lives
in **changepoint** [@killick2014changepoint] and **fpop**
[@maidstone2017optimal]; wild binary segmentation in **wbs** and
**breakfast** [@fryzlewicz2014wild; @fryzlewicz2020detecting]; multiscale
inference in **stepR** [@frick2014smuce]; Bayesian analysis in **bcp**
[@erdman2007bcp] and **ocp** [@adams2007bocpd]; structural breaks in
**strucchange** [@zeileis2002strucchange]; and so on. Each of these packages
is excellent at what it does, and each returns a different object, follows a
different indexing convention, and draws (or does not draw) its own plots.

An analyst who wants to *compare* a PELT segmentation with a Bayesian
posterior and a multiscale confidence set (a routine task in applied work)
must therefore learn several APIs, reconcile several conventions, and write
custom plotting code for each. **ggchangepoint** removes that friction. Its
design goals are:

1. **One vocabulary.** A single front door, `cpt_detect(x, method,
   change_in, penalty, ...)`, whose arguments mean the same thing for every
   engine.
2. **One result type.** Every detector returns a `ggcpt` object with a
   stable tidy contract, whatever the upstream engine returned.
3. **One rendering path.** Every result (point estimates, confidence
   intervals, fitted signals, posteriors, penalty paths) draws with
   `autoplot()` and extends with ordinary `ggplot2` layers.
4. **Wrap, don't reinvent.** All detection is delegated to the
   peer-reviewed upstream engines; optional engines live in `Suggests` and
   are loaded only when requested.

# The changepoint problem

Let $y_{1:n} = (y_1, \dots, y_n)$ be an ordered sequence. A segmentation with
$m$ changepoints is an ordered set $\tau_{1:m}$ satisfying
$0 = \tau_0 < \tau_1 < \dots < \tau_m < \tau_{m+1} = n$, which partitions the
data into the $m + 1$ segments $y_{(\tau_{i-1}+1):\tau_i}$,
$i = 1, \dots, m+1$. Both the number of changepoints and their locations are
unknown, and estimating them jointly is what makes the problem hard.
Throughout the package a changepoint $\tau$ is reported as the **last index
of the left segment** (the convention of the **changepoint** package), so the
admissible locations are $1, \dots, n-1$; results from engines using the
opposite convention are shifted on the way in, and the convention is recorded
on every result object.

## Penalised cost minimisation

The classical formulation chooses the segmentation minimising a penalised
cost,
$$
\min_{m,\ \tau_{1:m}} \; \sum_{i=1}^{m+1}
  \mathcal{C}\!\left(y_{(\tau_{i-1}+1):\tau_i}\right) \;+\; \beta m,
$$
where $\mathcal{C}$ is a segment cost (for a change in mean under Gaussian
noise the residual sum of squares, more generally twice the negative
maximised log-likelihood of the segment) and $\beta > 0$ is the price of
each additional changepoint. Writing $k$ for the number of parameters a
changepoint introduces, the familiar choices are $\beta = 2k$ (AIC) and
$\beta = k \log n$ (BIC, or SIC in the changepoint literature)
[@yao1988estimating], alongside the strengthened and modified variants
discussed below. Solved naively by dynamic programming, the minimisation
costs $O(n^2)$; **PELT** [@killick2012pelt] prunes candidate changepoints to
reach linear expected cost while remaining exact, and **FPOP**
[@maidstone2017optimal] reaches comparable speed by functional pruning.

A practical consequence of writing the objective this way is easy to miss:
$\beta$ and $\mathcal{C}$ must live on the same scale. For a change in mean
the **changepoint** engines evaluate the Normal cost with the noise standard
deviation fixed at 1, and **fpop** penalises the residual sum of squares
directly, so multiplying the data by a constant multiplies the cost while
leaving $\beta$ untouched. On a 200-point series with a single changepoint
whose jump is five standard deviations, `pelt` recovers exactly one
changepoint at $\sigma = 1$ but returns 39 at $\sigma = 3$ and 141 at
$\sigma = 10$ (means over 20 draws, since a single draw is not stable at
these settings). Standardise the series, pass a penalty on the data's own
scale (say `2 * log(n) * var(diff(x)) / 2`), or use `change_in = "meanvar"`,
which estimates a variance per segment. Methods that estimate the noise
level as part of their procedure (SMUCE, the WBS family, CPOP, bcp, BEAST
and the nonparametric engines) return the same segmentation whatever the
units. Three exceptions are worth knowing: `geomcp` runs PELT on its mapped
series and inherits its sensitivity, DeCAFS floors its noise estimate near
0.03, and BOCPD's default prior is on the data's own scale, so the last two
miss changes in a series measured in very small units.

## Search-based and multiscale methods

A complementary family locates changepoints by scanning test statistics.
**Binary segmentation** [@scott1974cluster; @vostrikova1981detecting]
recursively splits the series at the maximal CUSUM statistic; **wild binary
segmentation** [@fryzlewicz2014wild] and its successor WBS2
[@fryzlewicz2020detecting] draw random subintervals so that short segments
are not masked; **narrowest-over-threshold** (NOT)
[@baranowski2019narrowest] favours the narrowest interval on which the
contrast exceeds a threshold, which generalises cleanly to changes in slope;
**MOSUM** [@eichinger2018mosum] scans a moving-sum statistic at a fixed
bandwidth, or across a range of bandwidths; Isolate-Detect
[@anastasiou2022idetect] isolates each changepoint in an expanding interval;
and TGUH [@fryzlewicz2018tail] performs a tail-greedy bottom-up merge.
**SMUCE** [@frick2014smuce] occupies a special place: it estimates the step
function with the fewest jumps that still passes a *simultaneous multiscale
test* at level $\alpha$, and in doing so delivers confidence intervals for
every changepoint location, uncertainty statements most competitors cannot
make. HSMUCE [@pein2017hsmuce] extends this to heterogeneous noise.

## Beyond the mean

Changes need not be in the mean: the `change_in` argument accepts `"mean"`,
`"var"`, `"meanvar"`, `"slope"`, `"distribution"`, `"covariance"`,
`"network"`, `"regression"` and `"seasonality"`, the nine values
`cpt_methods()` lists in its `supports` column, each routable to the
methods that declare it.
Nonparametric engines (energy statistics [@matteson2014nonparametric],
nonparametric cost functions [@haynes2017computationally], kernel running
statistics [@arlot2019kernel; @cabrieto2018kcprs], joint characteristic
functions [@mcgonigle2023npmojo], self-normalisation [@zhao2022snseg])
detect distributional change without likelihood assumptions; Bayesian
engines [@barry1993bayesian; @adams2007bocpd; @zhao2019beast] return
posteriors instead of point sets; high-dimensional engines
[@wang2018inspect; @chen2022ocd; @grundy2020geomcp] aggregate evidence
across coordinates; and regression engines [@bai1998estimating;
@muggeo2003segmented] date breaks in model coefficients. The tour below
visits each family with runnable code.

# Design of the package

## The `ggcpt` result contract

Every detector returns an object of class `ggcpt` containing:

- `changepoints`: a tibble with one row per changepoint. Columns `cp`
  (location, "left" convention) and `cp_value` (the data value at `cp`) are
  always present; engines add `ci_lower`/`ci_upper` (SMUCE, HSMUCE,
  strucchange, segmented, bfast, taylor, mcp), `posterior_prob` (bcp,
  BEAST), `detection_time` (CPM), `strength` (inspect), `declared_at`
  (ocd), or `mapping` (geomcp) when they have more to say. `nsp` reports
  its uncertainty as a significance region instead: `region_start`,
  `region_end` and the `regions` slot.
- `segments`: a tibble of the induced segments (`seg_id`, `start`, `end`,
  `n`, `param_estimate`).
- `data`: the analysed series as a tibble (`index`, `value`), plus a
  `fitted` column when the engine estimates a signal (SMUCE, HSMUCE, CPOP,
  bcp, BEAST, DeCAFS, segmented, mcp, bfast; these are the engines
  `cpt_methods()` marks in its `fitted` column).
- `method`, `change_in`, `penalty` (a `list(type, value)` descriptor),
  `cp_convention` (always `"left"`), `runtime` (elapsed seconds, timed by
  `cpt_detect()` and `NA` when a wrapper is called directly), and `fit` (the
  untouched upstream object, for experts).

Multivariate results additionally carry a `data_wide` tibble with one column
per coordinate, which `autoplot()` renders as faceted small-multiples. Two
methods plot a derived series instead, because for them the coordinates are
not what a reader wants to see: `network` takes a sequence of adjacency
matrices and reports mean edge weight (one facet per matrix entry would be
unreadable), and `hdreg` plots the response it regressed on the covariates.
Both carry no `data_wide`.

## Tidy methods and the plotting layer

The class implements the full complement of generics R users expect:

```{r contract}
set.seed(2022)
x <- c(rnorm(100, 0, 1), rnorm(100, 10, 1))
res <- cpt_detect(x, method = "pelt", change_in = "mean")
res
tidy(res)
glance(res)
head(augment(res))
```

`autoplot()` draws the series, the changepoint rules and, on request,
the fitted segment means (`show_segments`), the engine's fitted signal
(`show_fit`), and changepoint-location confidence intervals (`show_ci`):

```{r contract-plot, fig.alt = "Series with its changepoint rules and the fitted segment levels drawn as horizontal steps"}
autoplot(res, show_segments = TRUE)
```

Composable layers (`geom_changepoint()`, `geom_cpt_segment()`,
`geom_cpt_ci()`, `stat_changepoint()`), a theme (`theme_ggcpt()`), and
segment shading (`annotate_segments()`) let the same results be built into
bespoke graphics; `summary()`, `as_tibble()`, `as.data.frame()`,
`format()`, and `plot()` complete the S3 surface.

## Design principles

The four goals above are recorded in the package as principles **P1. Wrap,
don't reinvent** (bind to peer-reviewed CRAN engines), **P2. Tidy in, tidy
out** (stable column names across all methods), **P3. ggplot2 all the way
down** (every result renders and extends), and **P4. One vocabulary** (`x`,
`method`, `change_in`, `penalty`, `...`). Three further principles govern how
the interface evolves:

- **P5. Progressive disclosure**: beginners call `cpt_detect()` +
  `autoplot()`; experts reach the upstream fit via `$fit`.
- **P6. No surprises**: the 0.1.0 functions still work unchanged.
- **P7. Document everything you ship**: every export is introduced in the
  README and a vignette.

Release 0.4.0 adds an eighth: **P8. Carry the uncertainty**. Where a method
quantifies uncertainty, the `ggcpt` object records it and `autoplot()` can
draw it.

Release 0.5.0 adds a ninth: **P9. Be extensible from the outside**. A
detector this package does not wrap, cannot wrap, or has never heard of can
join the same grammar through `as_ggcpt()` and `cpt_register_method()`, and
is labelled as user-supplied wherever it appears. See
`vignette("extending", package = "ggchangepoint")`.

# The unified dispatcher

`cpt_detect()` dispatches by method name; `cpt_methods()` reports every
method the package knows, its engine, what it can detect, and whether the
engine is installed:

```{r methods-table}
cpt_methods()
```

Requests are validated against this capability matrix: asking a mean-only
engine for a variance change is an error with the legal alternatives named,
never a silent substitution. Univariate methods likewise refuse multi-column
input rather than flattening it.

Penalty semantics differ across engines, and `cpt_penalty()` documents and
constructs the standard values:

```{r penalty}
cpt_penalty("BIC", n = 200)
cpt_penalty("MBIC", n = 200)
cpt_penalty("Hannan-Quinn", n = 200)
cpt_penalty("sSIC", n = 200)
```

Two of these warrant a word. `"sSIC"` is the strengthened Schwarz criterion
$k (\log n)^{\alpha}$, with $\alpha = 1.01$ by default
[@fryzlewicz2014wild], marginally heavier than BIC, and the criterion the
search-based engines apply internally. `"MBIC"` in `cpt_penalty()` returns
$0.5 (k+1) \log n + \log \binom{n}{k}$: a BIC-type term plus the
combinatorial cost of placing $k$ changepoints among $n$ observations. It is
deliberately stronger than `"BIC"`, but it is *not* the modified BIC of Zhang
and Siegmund (2007), whose penalty
$1.5 k \log n + 0.5 \sum_i \log(\ell_i / n)$ depends on the segment lengths
$\ell_i$ and therefore cannot be written as a function of $n$ and $k$ alone;
nor is it the quantity the **changepoint** package computes for its own
character penalty `"MBIC"`.

Character penalties (`"MBIC"`, `"BIC"`, ...) pass through to the
`changepoint`-family engines natively and are resolved to numeric values for
the functional-pruning engines (`fpop`, `cpop`, `decafs`). Search-based
engines (WBS, NOT, MOSUM, ...) select their own models and ignore the
argument, as do the engines tuned by a significance level, a
posterior-probability threshold, or an average run length (SMUCE, bcp, BEAST,
CPM, SNSeg).

# A tour of the method families

Throughout we use simulated series with known truth, so that results can be
checked by eye. `cpt_simulate()` draws series with prescribed changepoints,
and five canonical test signals ship as ready-made generators:
`signal_blocks()` (the Donoho-Johnstone blocks signal [@donoho1994ideal]),
`signal_fms()`, `signal_teeth()`, `signal_stairs()`, and `signal_mix()`. The
comparison vignette puts them to work. The sections that follow work through
six families (penalised and optimal partitioning, multiscale and search,
Bayesian, nonparametric and sequential, multivariate and high-dimensional, and
regression-based) plus two concerns that cut across all of them: change in
slope, and robustness to drift, autocorrelation and model ambiguity.

```{r tour-data}
set.seed(2026)
x_mean  <- c(rnorm(100), rnorm(100, 4))                # mean shift at 100
x_multi <- c(rnorm(100), rnorm(100, 3), rnorm(100, -1)) # shifts at 100, 200
x_slope <- cumsum(c(rep(0.4, 100), rep(-0.3, 100))) + rnorm(200) # kink at 100
```

## Penalised and optimal partitioning

PELT [@killick2012pelt], binary segmentation [@scott1974cluster], segment
neighbourhoods [@auger1989segment], and at-most-one-change (AMOC)
[@hinkley1970inference] come from the **changepoint** package
[@killick2014changepoint];
FPOP [@maidstone2017optimal] from **fpop**:

```{r penalised}
tidy(cpt_detect(x_multi, method = "pelt"))
tidy(cpt_detect(x_multi, method = "binseg"))
```

```{r penalised-fpop, eval = has_fpop}
tidy(cpt_detect(x_multi, method = "fpop"))
```

The Achilles heel of penalised methods is the choice of $\beta$. Rather
than committing to one value, `cpt_crops()` computes *every* optimal
segmentation as $\beta$ ranges over an interval (the CROPS algorithm of
Haynes, Eckley and Fearnhead (2017), as implemented by **changepoint**)
and turns penalty selection into a diagnostic:

```{r crops, fig.alt = "CROPS elbow plot: segmentation cost against the number of changepoints"}
path <- cpt_crops(x_multi)
path
autoplot(path)
```

The default plot puts segmentation cost against model size, and the usual
reading takes the model beyond which the cost stops falling appreciably. The
sweep over the default interval $[\log n,\, 10 \log n]$ admits only a handful
of distinct segmentations here, and the most parsimonious of them already
recovers the two true changepoints.
`autoplot(path, type = "segmentations")` shows the candidate models
themselves, and `autoplot(path, type = "path")` the map from penalty to model
size:

```{r crops-segmentations, fig.alt = "The series faceted by CROPS solution, each panel showing that solution's changepoints"}
autoplot(path, type = "segmentations")
```

The modern **fastcpd** engine [@li2024fastcpd] brings the same penalised
formulation to a wide family of models (mean, variance, mean-and-variance,
and AR/ARMA/GARCH model changes) with sequential-gradient-descent speed:

```{r fastcpd, eval = has_fastcpd}
tidy(fastcpd_wrapper(x_multi, family = "mean"))
```

## Multiscale and search methods

The randomised and multiscale searchers are one call each, whether through
`cpt_detect()` or through the wrapper directly:

```{r search-wbs, eval = has_wbs}
tidy(wbs_wrapper(x_multi, seed = 1))
```

```{r search-not, eval = has_not}
tidy(not_wrapper(x_multi, seed = 1))
```

```{r search-mosum, eval = has_mosum}
tidy(mosum_wrapper(x_multi))
```

```{r search-others, eval = has_idetect && has_breakfast}
tidy(idetect_wrapper(x_multi, seed = 1))
tidy(wbs2_wrapper(x_multi))
tidy(tguh_wrapper(x_multi))
```

SMUCE [@frick2014smuce] is the family's inferential flagship: its level
$\alpha$ bounds the probability of overestimating the number of changepoints
(the default is `alpha = 0.5`, **stepR**'s own recommendation for
estimation rather than testing), and every location comes with a confidence
interval, stored in `ci_lower`/`ci_upper` and drawn by `show_ci = TRUE` as
whiskers near the foot of the panel (the step fit is drawn by
`show_fit = TRUE`):

```{r smuce, eval = has_stepR, fig.alt = "SMUCE step fit with changepoint-location confidence intervals drawn as horizontal whiskers"}
res_smuce <- smuce_wrapper(x_multi)
tidy(res_smuce)
autoplot(res_smuce, show_ci = TRUE, show_fit = TRUE)
```

For heterogeneous noise, `smuce_wrapper(x, family = "hsmuce")` (or
`cpt_detect(x, method = "hsmuce")`) runs HSMUCE [@pein2017hsmuce].

## Changes in slope

A kink in the trend is not a jump in the level, and running a mean-change
detector on a trending series over-detects notoriously. CPOP
[@fearnhead2019cpop; @fearnhead2024cpop] solves the change-in-slope problem
*exactly* under an $L_0$ penalty, returning a continuous piecewise-linear
fit:

```{r cpop, eval = has_cpop, fig.alt = "Series with the CPOP piecewise-linear fit overlaid and its changepoints marked"}
res_cpop <- cpop_wrapper(x_slope)
tidy(res_cpop)
autoplot(res_cpop, show_fit = TRUE)
```

NOT with its linear contrast [@baranowski2019narrowest] offers a
search-based alternative; the dispatcher routes
`cpt_detect(x, method = "not", change_in = "slope")` to it automatically:

```{r not-slope, eval = has_not}
tidy(cpt_detect(x_slope, method = "not", change_in = "slope"))
```

## Bayesian detection

The Barry-Hartigan product partition model [@barry1993bayesian], via the
**bcp** package [@erdman2007bcp], returns a *posterior probability of a
changepoint at every location* along with posterior segment means. Locations
clearing `prob_threshold` populate the changepoints tibble (with their
probabilities), and `ggcpt_posterior()` draws the classic two-panel
display:

```{r bcp, eval = has_bcp, fig.alt = "Two-panel Bayesian display: the series with its posterior mean above, per-location posterior changepoint probability below"}
res_bcp <- bcp_wrapper(x_mean, seed = 2026)
tidy(res_bcp)
ggcpt_posterior(res_bcp)
```

Bayesian *online* changepoint detection [@adams2007bocpd] instead tracks the
posterior over the current **run length** (the time elapsed since the last
change), updating it recursively as each observation arrives. Its signature
graphic is the run-length heatmap, in which a change shows up as the
posterior mass falling back to a run length of zero:

```{r bocpd, eval = has_ocp, fig.alt = "Run-length heatmap: posterior probability of each run length over time"}
res_bocpd <- bocpd_wrapper(x_mean)
tidy(res_bocpd)
ggcpt_runlength(res_bocpd)
```

A third Bayesian engine, BEAST [@zhao2019beast] via **Rbeast**, averages
over models rather than conditioning on one, and is wired as
`cpt_detect(x, method = "beast")` (or `beast_wrapper()`); it too reports
`posterior_prob` and renders with `ggcpt_posterior()`.

## Nonparametric and sequential detection

When no parametric form is trustworthy, the nonparametric cost approach of
**changepoint.np** [@haynes2017computationally] and the energy-statistics
E-Divisive of **ecp** [@matteson2014nonparametric; @james2014ecp] detect
general distributional change:

```{r nonparam}
set.seed(2022)
tidy(cpt_detect(x_mean, method = "np"))
tidy(cpt_detect(x_mean, method = "ecp", seed = 1))
```

The **cpm** package [@ross2015cpm] recasts detection as a stream of
two-sample tests (Mann-Whitney for location, Mood for scale, Lepage,
Kolmogorov-Smirnov, Cramér-von Mises, and parametric variants), run here over
the whole series in one pass to mimic an online monitor. Its results
distinguish where a change *happened* (`cp`) from when it was *detected*
(`detection_time`), the lag inherent in sequential monitoring:

```{r cpm, eval = has_cpm}
tidy(cpm_wrapper(x_mean, cpm_type = "Mann-Whitney"))
```

Three further nonparametric engines are wired and worth knowing: kernel
change-point analysis on running statistics (`kcp_wrapper()`, engine
**kcpRS**), which detects changes in running means, variances,
autocorrelations, or correlations [@arlot2019kernel; @cabrieto2018kcprs];
NP-MOJO (`npmojo_wrapper()`, engine **CptNonPar**), which detects changes in
the marginal or lagged joint distribution while remaining valid under serial
dependence [@mcgonigle2023npmojo]; and self-normalised segmentation
(`sn_wrapper()`, engine **SNSeg**), which avoids long-run variance
estimation altogether and tests changes in means, variances,
autocorrelations, or bivariate correlations [@zhao2022snseg].

## Robustness to drift, autocorrelation, and model ambiguity

The most common failure of mean-change detection in practice is not a subtle
statistical one: it is running a Gaussian-mean detector on data whose
baseline drifts or whose noise is autocorrelated, and then reporting a
changepoint wherever the model is wrong. DeCAFS [@romano2022decafs] models
exactly this regime (abrupt changes superimposed on random-walk drift and
AR(1) noise) and separates the two:

```{r decafs, eval = has_decafs, fig.alt = "Series with the DeCAFS fit overlaid, which separates gradual drift from abrupt change"}
res_decafs <- decafs_wrapper(x_mean)
tidy(res_decafs)
autoplot(res_decafs, show_fit = TRUE)
```

EnvCpt [@beaulieu2018envcpt] attacks the same confusion by model selection:
it fits up to twelve competing descriptions (constant mean or linear trend,
each with or without changepoints, and with white-noise, AR(1) or AR(2)
errors) and reports changepoints only if a changepoint model wins on an
information criterion:

```{r envcpt, eval = has_envcpt}
res_env <- envcpt_wrapper(x_mean, models = c("mean", "meancpt", "trendcpt"))
glance(res_env)
```

The winning model's name is recorded in the penalty descriptor
(`penalty_type` above, here `AIC: meancpt`), so "no changepoints, it's just
autocorrelation" is a first-class answer.

## Multivariate and high-dimensional detection

Multivariate methods accept a matrix (rows are time points) directly. The
energy-statistics E-Divisive of **ecp** was built for this
[@matteson2014nonparametric]; for high-dimensional data whose change is
confined to a sparse subset of coordinates, `inspect` [@wang2018inspect]
finds an optimal sparse projection of the CUSUM matrix and reports the
projected evidence (`strength`). Multivariate results render as faceted
small-multiples with shared changepoint rules:

```{r inspect, eval = has_inspect, fig.alt = "Faceted small-multiples, one panel per coordinate, sharing the detected changepoint rules"}
set.seed(2026)
X <- cbind(a = c(rnorm(80), rnorm(80, 3)),
           b = c(rnorm(80), rnorm(80, -2)),
           c = rnorm(160))
res_hd <- inspect_wrapper(X)
tidy(res_hd)
autoplot(res_hd)
```

Two further engines complete the family. `geomcp_wrapper()` (engine
**changepoint.geo**) maps each observation to its distance from and its angle
to a reference point, then segments the two mapped series, catching changes in
magnitude and in orientation respectively [@grundy2020geomcp]. And
`ocd_wrapper()` (engine **ocd**) monitors a high-dimensional stream *online*
with worst-case detection-delay guarantees [@chen2022ocd]. Because detection
there is sequential, the locations it reports are *declaration times* (the
change plus the detection delay), recorded in `declared_at`; the wrapper
estimates the pre-change baseline from an initial training window and resets
after each declaration so that several changes can be found. Like the method
itself, it needs at least two coordinates and refuses a single series.

## Structural breaks in regression

Econometric practice dates breaks in regression coefficients. The
Bai-Perron estimator [@bai1998estimating; @bai2003computation], via
**strucchange** [@zeileis2002strucchange], returns break dates *with
confidence intervals*; called on a bare series it dates mean shifts, and
called with a formula it dates breaks in arbitrary regressions:

```{r strucchange, eval = has_strucchange, fig.alt = "Series with the Bai-Perron breakpoints and their confidence intervals drawn as horizontal whiskers"}
res_bp <- strucchange_wrapper(x_mean)
tidy(res_bp)
autoplot(res_bp, show_ci = TRUE)
```

Where the regression function is continuous (a kink rather than a jump),
**segmented** [@muggeo2003segmented; @muggeo2008segmented] estimates
broken-line relationships with standard errors for the breakpoints:

```{r segmented, eval = has_segmented, fig.alt = "Series with the broken-line fit, its breakpoint, and the breakpoint's confidence interval"}
res_seg <- segmented_wrapper(x_slope, npsi = 1, seed = 1)
tidy(res_seg)
autoplot(res_seg, show_fit = TRUE, show_ci = TRUE)
```

# Beyond detection

## Batch detection over many series

Applied work rarely stops at one series. `cpt_batch()` runs one detector over
every column of a matrix or data frame (or every element of a list) and
returns a tibble with one row per series, carrying both the tidy changepoints
and the full `ggcpt` object in list-columns. With **future** and
**future.apply** installed it honours a non-sequential `future::plan()`, using
parallel-safe RNG:

```{r batch, fig.alt = "Small-multiples of a panel of series, each with its own detected changepoints"}
set.seed(2026)
panel <- cbind(shifted = x_mean, quiet = rnorm(200))
batch <- cpt_batch(panel, method = "pelt")
batch
tidy(batch)
autoplot(batch)
```

## Stability diagnostics

Most engines report a point set with no measure of its fragility.
`cpt_stability()` resamples residuals *within* the fitted segments (so the
estimated regime structure is preserved), re-runs the detector on each
replicate, and reports how often each location is re-detected: a cheap,
model-agnostic confidence signal available for *every* engine, including the
many that ship no intervals of their own:

```{r stability, fig.alt = "Bootstrap detection-frequency profile across the series, with the original changepoints marked"}
st <- cpt_stability(x_mean, method = "pelt", B = 50, seed = 1)
st
autoplot(st)
```

## Evaluation, interactivity, and citations

When ground truth is known, `cpt_metrics()` computes precision, recall and
F1 under one-to-one matching, the covering metric, Hausdorff distance, and
adjusted Rand index, following the conventions of the modern benchmarking
literature [@van2020evaluation]; `ggcpt_eval()` draws the agreement, and
`ggcpt_compare()` juxtaposes methods. These are the subject of the
companion vignette `vignette("comparison", package = "ggchangepoint")`.

Any result renders as an interactive HTML widget with
`ggcpt_interactive(res)` (engine **plotly**, in `Suggests`); the static
`autoplot()` path is untouched.

Finally, because every method here is someone's published work, `cpt_cite()`
returns the reference(s) behind a result, so analyses can cite the right
paper without leaving R:

```{r cite}
cpt_cite("pelt")
```

# Discussion

ggchangepoint does not contribute a new detection algorithm; it contributes a
*surface*. The value of a common contract compounds with the number of
methods behind it: the same `tidy()` pipeline, the same plot, and the same
evaluation code now span penalised, multiscale, nonparametric, Bayesian,
high-dimensional, and regression-based detection: 50 methods in this
release. Five more, whose engines are not currently on CRAN
(graph-constrained gfpop [@hocking2020gfpop], robust segmentation under
outliers [@fearnhead2019changepoint], FOCuS, sparsified binary segmentation,
and random-forest classification [@londschien2023changeforest]), are listed as
*planned* in `cpt_methods()` and will slot into the same wrapper pattern once
their engines return; until then they are not callable.

Two practical notes. First, wrapped engines run with sensible defaults, but
every wrapper forwards `...` to its engine and the raw fit is always in
`$fit`; the package is a front door, not a cage. Second, detection quality
belongs to the engines; the package's own additions (metrics, stability,
penalty paths) are deliberately engine-agnostic, so conclusions drawn with
them transfer between methods.

# Acknowledgements

This package stands on the shoulders of the authors of the wrapped engines
and of the R [@rcore], ggplot2 [@wickham2016ggplot2], and broom
[@robinson2017broom] projects.

# References