Network-structure changepoints via non-negative matrix factorisation

Wraps fabisearch::detect.cps() (Ondrus, Olds and Cribben, 2024): factorised binary search for changes in the network structure of a high-dimensional series. Each candidate split is scored by how much better a rank-$r$ non-negative matrix factorisation fits the two halves separately than together, and significance is assessed by permutation.

Usage

fabisearch_wrapper(
  x,
  min_dist = 35,
  n_runs = 50,
  n_reps = 100,
  alpha = NULL,
  rank = NULL,
  n_core = 1,
  seed = NULL,
  ...
)

Arguments

x: A non-negative numeric matrix or data frame with rows as time points and columns as nodes. Non-negativity is a hard requirement of NMF, not a preference; see the section below.
min_dist: Minimum distance between changepoints. Defaults to 35 (the engine's default), lowered automatically when the series is too short for it.
n_runs: NMF runs per candidate split. Defaults to 50.
n_reps: Replicates on each side of the significance test: the split is refitted n_reps times and the rows permuted n_reps times. Defaults to 100; at least 2, because the test compares two samples of this size.
alpha: Significance level applied to the p-value each candidate split receives. Defaults to 0.05. That p-value is a two-sample test of the refitted losses against the permuted ones (a t-test unless testtype is passed through ...), adjusted by Benjamini-Hochberg across the candidate splits. It is not a permutation p-value, so it has no 1 / n_reps floor; the exact rank tests (testtype = "wilcox" or "ks") do have one, and the wrapper warns when it is above alpha.
rank: NMF rank, a positive whole number. NULL estimates it with fabisearch::opt.rank(), which is expensive; supplying a rank is much faster.
n_core: Cores for the permutation stage. Defaults to 1.
seed: Optional seed. The seed is scoped to this call: .Random.seed is saved and restored, so a seeded call inside a simulation loop does not pin the loop's own stream.
...: Additional arguments passed to fabisearch::detect.cps().

Value

A ggcpt object with change_in = "network". Multivariate input is reduced to one series per observation by taking the cross-sectional mean of the columns, and that is the series stored on the result: autoplot() draws it, tidy()'s cp_value reads it, and $segments$param_estimate and augment()'s .fitted/.resid are computed from it. It is not any one column of the input. The full input is kept in $data_wide for autoplot(type = "coordinates").

Non-negativity, cost, and the attached namespace

Three practical notes. (1) NMF is undefined for negative entries, so this wrapper refuses them rather than letting the engine fail deep inside a factorisation; shift or rescale the series first if it has negatives. (2) It is by far the most expensive engine here (n_runs times n_reps factorisations per candidate split), so the defaults are lowered in the examples. The engine's own progress output is suppressed, so a long call is silent until it returns. (3) fabisearch calls NMF's multi-run machinery, which resolves helpers through the search path and fails with "none of the packages are loaded" when NMF is merely loaded; this wrapper therefore attaches NMF for the duration of the call and detaches it again afterwards.

References

Ondrus M, Cribben I (2024). “fabisearch: a package for change point detection in and visualization of the network structure of multivariate high-dimensional time series in R.” Neurocomputing, 578, 127321. doi:10.1016/j.neucom.2024.127321 .

Examples

# \donttest{
# A change in *structure*, not in scale: two latent factors drive
# different halves of the node set before and after the change.
# Deliberately tiny -- this is by far the most expensive engine in the
# package (n_runs x n_reps factorisations per candidate split). Measured
# at 5-6 s across fresh sessions, against 27 s for the 2 x 25 /
# n_reps = 4
# version this replaced -- and this is the floor: `n_reps = 1` fails
# inside fabisearch with "not enough 'x' observations" (the permutation
# test needs two), and smaller matrices are not reliably cheaper because
# the search then evaluates more splits relative to `min_dist` (2 x 10 at
# min_dist = 8 measured 6.6 s). So this one example stays near CRAN's 5 s
# budget by necessity; `cran-comments.md` says so. Use the defaults on
# real data -- the settings here are for the budget, not for detection.
set.seed(2026)
block <- function(n, cols) {
  f <- abs(stats::rnorm(n)) + 0.5
  Y <- matrix(abs(stats::rnorm(n * 5)) * 0.2 + 0.1, n, 5)
  Y[, cols] <- Y[, cols] + f
  Y
}
Y <- rbind(block(12, 1:2), block(12, 3:5))
fabisearch_wrapper(Y, min_dist = 10, n_runs = 1, n_reps = 2,
                   alpha = 0.25, rank = 2)
#> ggcpt (changepoint detection result)
#>   Method:             fabisearch
#>   Change in:          network
#>   Changepoints found: 1
#>   CP convention:      left
#>   Penalty:            engine alpha = 0.25
#>   Series length:      24
#> 
#> Changepoints:
#> # A tibble: 1 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1    12    0.754
# }