--- title: "The ggchangepoint Feature Tour: A Complete Map of the Package Surface" author: "Youzhi Yu
University of Chicago" bibliography: vignette_reference.bib output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{The ggchangepoint Feature Tour: A Complete Map of the Package Surface} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 8, fig.height = 5, message = FALSE, warning = FALSE, fig.alt = "ggchangepoint feature tour plot" ) library(ggchangepoint) library(ggplot2) theme_set(theme_light()) # Optional engines live in Suggests; every chunk that needs one is guarded so # the vignette builds on any installation. has_fpop <- requireNamespace("fpop", quietly = TRUE) has_wbs <- requireNamespace("wbs", quietly = TRUE) has_breakfast <- requireNamespace("breakfast", quietly = TRUE) has_not <- requireNamespace("not", quietly = TRUE) has_mosum <- requireNamespace("mosum", quietly = TRUE) has_idetect <- requireNamespace("IDetect", quietly = TRUE) has_stepR <- requireNamespace("stepR", quietly = TRUE) has_cpop <- requireNamespace("cpop", quietly = TRUE) has_bcp <- requireNamespace("bcp", quietly = TRUE) has_ocp <- requireNamespace("ocp", quietly = TRUE) # Rbeast (<= 1.0.2) can crash the vignette-building subprocess on Windows # (its C internals have known memory-state issues); run it elsewhere only. has_rbeast <- requireNamespace("Rbeast", quietly = TRUE) && .Platform$OS.type != "windows" has_cpm <- requireNamespace("cpm", quietly = TRUE) has_kcprs <- requireNamespace("kcpRS", quietly = TRUE) has_cptnonpar <- requireNamespace("CptNonPar", quietly = TRUE) has_decafs <- requireNamespace("DeCAFS", quietly = TRUE) has_snseg <- requireNamespace("SNSeg", quietly = TRUE) has_inspect <- requireNamespace("InspectChangepoint", quietly = TRUE) has_geomcp <- requireNamespace("changepoint.geo", quietly = TRUE) has_struc <- requireNamespace("strucchange", quietly = TRUE) has_segmented <- requireNamespace("segmented", quietly = TRUE) has_envcpt <- requireNamespace("EnvCpt", quietly = TRUE) has_fastcpd <- requireNamespace("fastcpd", quietly = TRUE) has_plotly <- requireNamespace("plotly", quietly = TRUE) ``` # Abstract **ggchangepoint** provides a unified, tidy, `ggplot2`-native interface to changepoint detection in R: one dispatcher (`cpt_detect()`) covering 31 methods across six methodological families, one result class (`ggcpt`) with a stable tidy contract, and one visualisation entry point (`autoplot()`) that draws everything a method reports — including confidence intervals and posterior probabilities [@wickham2016ggplot2; @robinson2017broom]. This vignette is the *feature tour*: it visits **every exported function** in the package at the point where it belongs in the workflow, so a reader can map the full surface in one sitting. The companion vignettes develop the methodology in depth (`vignette("introduction")`) and treat method comparison and evaluation (`vignette("comparison")`). ```{r data} set.seed(2026) x <- c(rnorm(100, 0, 1), rnorm(100, 6, 1)) # one mean shift at t = 100 x3 <- c(rnorm(100), rnorm(100, 4), rnorm(100, 1)) # shifts at 100, 200 ``` # The result object and its methods Every detector returns a `ggcpt` object: a list carrying the tidy `changepoints` tibble (`cp` = last index of the left segment, `cp_value` = series value there, plus any method-specific columns), a `segments` table, the `data`, the raw engine `fit`, and metadata (method, `change_in`, penalty, convention, runtime). `new_ggcpt()` is the low-level constructor and `is_ggcpt()` the class test; most users never call either directly. ```{r class} res <- cpt_detect(x, method = "pelt", change_in = "mean") is_ggcpt(res) print(res) ``` ```{r constructor} manual <- new_ggcpt( changepoints = tibble::tibble(cp = 100L, cp_value = x[100]), data = tibble::tibble(index = seq_along(x), value = x), method = "manual" ) is_ggcpt(manual) ``` The `broom` verbs give one row per changepoint (`tidy()`), a one-row model summary (`glance()`), and the data augmented with segment ids, fitted levels, residuals, and a changepoint flag (`augment()`): ```{r broom} tidy(res) glance(res) head(augment(res)) ``` The remaining S3 surface: a human-readable digest, tibble/data-frame coercion, a one-line format, and a base-`plot()` fallback that delegates to `autoplot()`. ```{r s3} summary(res) as_tibble(res) head(as.data.frame(res)) format(res) ``` ```{r plot-fallback} plot(res) ``` # Unified detection `cpt_detect()` is the recommended entry point: pick a `method`, say what the change is in (`change_in` = `"mean"`, `"var"`, `"meanvar"`, `"slope"`, or `"distribution"`), and optionally set a `penalty`. Incompatible `method`/`change_in` combinations error — they are never silently substituted. ```{r detect} cpt_detect(x, method = "binseg", change_in = "mean") ``` Anything else in `...` reaches the underlying wrapper, and takes precedence over the value the dispatcher would otherwise derive from `change_in` — here the NOT contrast is set directly rather than inherited: ```{r detect-dots, eval = has_not} tidy(cpt_detect(x, method = "not", contrast = "pcwsLinMean")) ``` `cpt_methods()` is the live capability table: every method the package knows, its engine package, what it can detect, and whether the engine is installed. Four rows carry status `"planned"` rather than `"available"`: their engines (`gfpop`, `robseg`, `FOCuS`, `hdbinseg`) are not on CRAN, so they are documented as future work and are *not* wired to `cpt_detect()`. ```{r methods} print(cpt_methods(), n = Inf) ``` `cpt_penalty()` constructs standard penalty values for the engines that take numeric penalties; see its help for the per-engine penalty semantics. ```{r penalty} cpt_penalty("BIC", n = 200) cpt_penalty("MBIC", n = 200) cpt_penalty("Manual", value = 10) ``` # Engine wrappers Each engine also has a direct wrapper exposing its native arguments. Every wrapper from 0.2.0 onwards returns a `ggcpt` object; only the two original wrappers below predate the class and still return a bare tibble. ## The classical core (0.1.0) `cpt_wrapper()` and `ecp_wrapper()` are the original interface to the *changepoint*/*changepoint.np* [@killick2012pelt; @killick2014changepoint; @haynes2017computationally] and *ecp* [@matteson2014nonparametric; @james2014ecp] engines; they return bare tibbles for backward compatibility. ```{r classic} cpt_wrapper(x, change_in = "mean", cp_method = "PELT") ecp_wrapper(x, algorithm = "divisive", seed = 1) ``` ## The search and pruning wave (0.2.0) Seven multiscale, randomised, and functional-pruning engines [@fryzlewicz2014wild; @fryzlewicz2020detecting; @baranowski2019narrowest; @eichinger2018mosum; @anastasiou2022idetect; @fryzlewicz2018tail; @maidstone2017optimal], one call each: ```{r fpop, eval = has_fpop} fpop_wrapper(x, penalty = 2 * log(length(x))) ``` ```{r wbs, eval = has_wbs} wbs_wrapper(x, n_intervals = 2000, seed = 1) ``` ```{r wbs2, eval = has_breakfast} wbs2_wrapper(x) ``` ```{r not, eval = has_not} not_wrapper(x, contrast = "pcwsConstMean", seed = 1) ``` ```{r mosum, eval = has_mosum} mosum_wrapper(x) mosum_wrapper(x3, multiscale = TRUE) ``` ```{r idetect, eval = has_idetect} idetect_wrapper(x, seed = 1) ``` ```{r tguh, eval = has_breakfast} tguh_wrapper(x) ``` ## The engine wave (0.4.0) **Multiscale inference with confidence.** `smuce_wrapper()` implements SMUCE [@frick2014smuce], which bounds the probability of over-estimating the number of changes at the level `alpha` and returns a confidence interval for every changepoint location (`ci_lower`/`ci_upper`); `family = "hsmuce"` switches to the heteroskedastic extension HSMUCE [@pein2017hsmuce], which estimates a variance per segment: ```{r smuce, eval = has_stepR} res_smuce <- smuce_wrapper(x, alpha = 0.5) tidy(res_smuce) ``` **Change in slope.** `cpop_wrapper()` performs exact penalised estimation of a continuous piecewise-linear mean [@fearnhead2019cpop; @fearnhead2024cpop]: ```{r cpop, eval = has_cpop} y_slope <- cumsum(c(rep(0.4, 100), rep(-0.3, 100))) + rnorm(200) res_cpop <- cpop_wrapper(y_slope) tidy(res_cpop) ``` ```{r cpop-fallback, include = FALSE, eval = !has_cpop} # y_slope is reused by the segmented example below y_slope <- cumsum(c(rep(0.4, 100), rep(-0.3, 100))) + rnorm(200) ``` **Bayesian detection.** `bcp_wrapper()` implements the Barry–Hartigan product-partition model [@barry1993bayesian; @erdman2007bcp]; `bocpd_wrapper()` runs Bayesian online changepoint detection over the run-length posterior [@adams2007bocpd]; `beast_wrapper()` wraps the BEAST Bayesian model-averaging ensemble [@zhao2019beast]. The first and third keep the locations whose posterior probability clears `prob_threshold` and record it in a `posterior_prob` column; BOCPD returns the maximum a posteriori changepoint set. ```{r bcp, eval = has_bcp} res_bcp <- bcp_wrapper(x, seed = 1) tidy(res_bcp) ``` ```{r bocpd, eval = has_ocp} res_bocpd <- bocpd_wrapper(x) tidy(res_bocpd) ``` ```{r beast, eval = has_rbeast} res_beast <- beast_wrapper(x, seed = 1) tidy(res_beast) ``` **Sequential and kernel nonparametrics.** `cpm_wrapper()` runs distribution-free sequential tests and reports, alongside each estimated location, the `detection_time` at which a stream monitor would have flagged it [@ross2015cpm]; `kcp_wrapper()` applies kernel change point analysis to running statistics (mean, variance, autocorrelation, correlation) [@arlot2019kernel; @cabrieto2018kcprs]; `npmojo_wrapper()` detects distributional changes under serial dependence [@mcgonigle2023npmojo]: ```{r cpm, eval = has_cpm} tidy(cpm_wrapper(x, cpm_type = "Mann-Whitney")) ``` ```{r kcp, eval = has_kcprs} tidy(kcp_wrapper(x, running_stat = "mean", nperm = 100, seed = 1)) ``` ```{r npmojo, eval = has_cptnonpar} tidy(npmojo_wrapper(x)) ``` **Robustness to drift and dependence.** `decafs_wrapper()` detects abrupt changes when the signal also drifts and the noise is autocorrelated [@romano2022decafs]; `sn_wrapper()` segments a chosen `parameter` (mean, variance, autocorrelation, or bivariate correlation) by self-normalisation, avoiding long-run variance estimation entirely [@zhao2022snseg]; `envcpt_wrapper()` reports changepoints only when a changepoint model beats the trend and autoregressive alternatives [@beaulieu2018envcpt]: ```{r decafs, eval = has_decafs} tidy(decafs_wrapper(x)) ``` ```{r sn, eval = has_snseg} tidy(sn_wrapper(x3, parameter = "mean")) ``` ```{r envcpt, eval = has_envcpt} res_env <- envcpt_wrapper(x, models = c("mean", "meancpt", "trendcpt")) tidy(res_env) res_env$penalty$type # criterion and winning model ``` **High-dimensional and multivariate.** These wrappers take a matrix or data frame with one row per time point. `inspect_wrapper()` finds sparse mean changes by projection [@wang2018inspect]; `ocd_wrapper()` monitors a high-dimensional stream online [@chen2022ocd]; `geomcp_wrapper()` maps each observation to a distance and an angle and segments both [@grundy2020geomcp]; multivariate `ecp` input flows through `cpt_detect()` unchanged. The univariate wrappers, by contrast, reject a multi-column argument rather than silently flattening it. ```{r hd-data} set.seed(1) X <- cbind(a = c(rnorm(80), rnorm(80, 4)), b = c(rnorm(80), rnorm(80, -3)), c = rnorm(160)) ``` ```{r inspect, eval = has_inspect} res_hd <- inspect_wrapper(X) tidy(res_hd) ``` ```{r ocd, eval = FALSE} # Online detection: the reported locations are declaration times (the change # plus the detection delay), in a `declared_at` column. Monte Carlo threshold # calibration makes this the slowest wrapper, so it is shown but not run # here. It needs at least two coordinates and rejects a univariate series. res_ocd <- ocd_wrapper(X, mc_reps = 100) tidy(res_ocd) ``` ```{r geomcp, eval = has_geomcp} tidy(geomcp_wrapper(X)) ``` ```{r ecp-mv} tidy(cpt_detect(X, method = "ecp", seed = 1)) ``` **Regression structure.** `strucchange_wrapper()` dates Bai–Perron breaks, either in a bare series or in the coefficients of a formula [@bai1998estimating; @bai2003computation; @zeileis2002strucchange]; `segmented_wrapper()` fits a continuous broken-line regression, so the change it reports is a kink in the trend rather than a jump in the level [@muggeo2003segmented; @muggeo2008segmented]. Both carry a confidence interval for every break: ```{r strucchange, eval = has_struc} tidy(strucchange_wrapper(x)) ``` ```{r segmented, eval = has_segmented} tidy(segmented_wrapper(y_slope, npsi = 1, seed = 1)) ``` **The modern PELT family.** `fastcpd_wrapper()` exposes the *fastcpd* engine, whose `family` argument covers changes in the mean, the variance, or both, as well as changes in a fitted AR/ARMA/GARCH model [@li2024fastcpd]: ```{r fastcpd, eval = has_fastcpd} tidy(fastcpd_wrapper(x, family = "mean")) ``` # The penalty path Rather than guessing one penalty, `cpt_crops()` computes *every* optimal segmentation over a penalty interval (the CROPS algorithm, via *changepoint* [@killick2014changepoint]) and returns a `ggcpt_path` object with its own `print()`, `tidy()`, and three `autoplot()` types: ```{r crops} path <- cpt_crops(x3) path tidy(path) ``` ```{r crops-elbow} autoplot(path) # the cost elbow ``` ```{r crops-path} autoplot(path, type = "path") # changepoints vs penalty ``` ```{r crops-seg} autoplot(path, type = "segmentations") # the candidate models themselves ``` # The visualisation layer `autoplot()` renders any `ggcpt`. Options: `show_segments` (fitted segment means), `show_fit` (the engine's own fitted signal, where provided — SMUCE, DeCAFS, CPOP, segmented, bcp, BEAST), `show_ci` (changepoint-location confidence intervals, where provided — SMUCE/HSMUCE, strucchange, segmented), `show_points`/`show_line`, an `index` for a date axis, and the `cptline_*` styling arguments (`cptline_color`, `cptline_alpha`, `cptline_type`, `cptline_linewidth`). Asking for an overlay the result cannot supply warns rather than failing silently. ```{r autoplot} autoplot(res, show_segments = TRUE, cptline_color = "firebrick", cptline_type = "dashed", cptline_linewidth = 0.8) ``` ```{r autoplot-ci, eval = has_stepR} autoplot(res_smuce, show_ci = TRUE, show_fit = TRUE) ``` Multivariate results facet automatically: ```{r autoplot-mv, eval = has_inspect, fig.height = 6} autoplot(res_hd) ``` The composable layers work inside any ggplot pipeline: `geom_changepoint()` draws vertical rules, `geom_cpt_segment()` draws segment levels, `geom_cpt_ci()` draws horizontal interval whiskers, and `stat_changepoint()` runs detection *inside* the plot. `theme_ggcpt()` is a publication theme and `annotate_segments()` shades alternating segments. ```{r geoms} cp_tbl <- tidy(res) df <- data.frame(index = seq_along(x), value = x) ggplot(df, aes(index, value)) + annotate_segments(cp = cp_tbl$cp, n = length(x)) + geom_line() + geom_changepoint(data = cp_tbl, aes(xintercept = cp), color = "red") + geom_cpt_segment( data = res$segments, aes(x = start, xend = end, y = param_estimate, yend = param_estimate), inherit.aes = FALSE, color = "darkred", linewidth = 1 ) + theme_ggcpt() ``` ```{r stat} ggplot(df, aes(index, value)) + geom_line() + stat_changepoint(method = "pelt", color = "blue") ``` ```{r geom-ci, eval = has_stepR} ci_tbl <- tidy(res_smuce) ci_tbl$y_pos <- min(x) - 1 ggplot(df, aes(index, value)) + geom_line(color = "grey60") + geom_cpt_ci( data = ci_tbl, aes(y = y_pos, xmin = ci_lower, xmax = ci_upper), width = 0.6, color = "blue", inherit.aes = FALSE ) + geom_changepoint(data = ci_tbl, aes(xintercept = cp), color = "blue") ``` `ggcptplot()` and `ggecpplot()` are the original one-call plots for the two classical engines: ```{r original-plots} ggcptplot(x, change_in = "mean", cp_method = "PELT") ``` ```{r ggecpplot} ggecpplot(x, algorithm = "divisive", seed = 1) ``` The Bayesian engines get the field's signature displays: `ggcpt_posterior()` shows the posterior mean over the series and the per-location changepoint probability (for `bcp` and BEAST results); `ggcpt_runlength()` shows the BOCPD run-length posterior as a heatmap. ```{r posterior, eval = has_bcp, fig.height = 6} ggcpt_posterior(res_bcp) ``` ```{r runlength, eval = has_ocp} ggcpt_runlength(res_bocpd) ``` Finally, `ggcpt_interactive()` turns any result (or any ggplot built from one) into a `plotly` widget with values on hover. The widget itself is not embedded here, only its class: ```{r interactive, eval = has_plotly} class(ggcpt_interactive(res)) # requires plotly ``` # Comparing methods `ggcpt_compare()` runs several detectors on the same series and renders them faceted (default) or overlaid, honouring `future::plan()` for parallel execution; `ggcpt_compare_table()` returns the tidy union of the same runs. A method that finds nothing keeps its panel and contributes an `NA` row — "no changepoints" is a result, not a missing one. ```{r compare, fig.height = 6} ggcpt_compare(x, methods = c("pelt", "binseg", "amoc")) ``` ```{r compare-overlay} ggcpt_compare(x, methods = c("pelt", "binseg"), layout = "overlay") ``` ```{r compare-table} ggcpt_compare_table(x, methods = c("pelt", "binseg", "amoc")) ``` # Batch detection and stability `cpt_batch()` runs one detector over many series (matrix, data frame, or list of vectors), also under `future::plan()`, and returns a tidy `ggcpt_batch` tibble with list-columns for the per-series changepoints and the `ggcpt` objects themselves: ```{r batch} XB <- cbind(shifted = x, pure_noise = rnorm(200)) batch <- cpt_batch(XB, method = "pelt") batch tidy(batch) autoplot(batch) ``` `cpt_stability()` quantifies how fragile a segmentation is: it resamples residuals within the fitted segments, re-runs the detector, and reports the re-detection frequency at every location — a model-agnostic confidence signal for engines with no native intervals: ```{r stability} st <- cpt_stability(x, method = "pelt", B = 30, seed = 1) st autoplot(st) ``` # Evaluation against ground truth `cpt_metrics()` scores predictions against a known truth: precision, recall and F1 under one-to-one matching, the covering metric and adjusted Rand index in the conventions of @van2020evaluation, Hausdorff distance, annotation error, and matched MAE/RMSE. `cpt_metrics_annotated()` averages over multiple annotators, and `ggcpt_eval()` draws the agreement (true positives, false positives, and misses): ```{r metrics} truth <- 100 pred <- tidy(res)$cp cpt_metrics(pred, truth, n = length(x), margin = 5) cpt_metrics_annotated(pred, list(100, 101, 99), n = length(x)) ``` ```{r eval-plot} ggcpt_eval(pred, truth, x, margin = 5) ``` # Simulation and canonical signals `cpt_simulate()` (alias `rcpt()`) generates series with known changepoints in mean, variance, both, or slope, under Gaussian, Student-t, AR(1), or random-walk noise; the truth travels with the data in the `true_changepoints` and `true_segments` attributes. Five canonical test signals from the literature ship ready-made: `signal_blocks()` (the Donoho–Johnstone blocks signal [@donoho1994ideal]), `signal_fms()`, `signal_mix()`, `signal_teeth()`, and `signal_stairs()`. ```{r simulate} sim <- cpt_simulate(300, changepoints = c(100, 200), change_in = "mean", params = c(0, 5, 1), seed = 1) attr(sim, "true_changepoints") sim2 <- rcpt(300, changepoints = 150, params = c(0, 3), seed = 2) # the alias attr(sim2, "true_changepoints") signals <- list(blocks = signal_blocks(1024, seed = 1), fms = signal_fms(500, seed = 1), mix = signal_mix(500, seed = 1), teeth = signal_teeth(400, seed = 1), stairs = signal_stairs(500, seed = 1)) vapply(signals, function(s) length(attr(s, "true_changepoints")), integer(1)) ``` ```{r blocks-plot} blocks <- signals$blocks ggplot(blocks, aes(index, value)) + geom_line(color = "grey50") + geom_vline(xintercept = attr(blocks, "true_changepoints"), color = "blue", linewidth = 0.3) + labs(title = "The blocks test signal with its true changepoints") ``` # Citing the methodology `cpt_cite()` returns the verified reference(s) behind a result or a method name, so a write-up can cite the right paper without leaving R; called with no argument it returns the whole method-to-reference table. ```{r cite} cpt_cite("pelt") cpt_cite(res) ``` # Closing note This tour visited every exported function in the package. For the framework's design, the mathematics of the wrapped methods, and worked analyses, see `vignette("introduction", package = "ggchangepoint")`; for method comparison and accuracy evaluation in depth, see `vignette("comparison", package = "ggchangepoint")`. # References