Asks whether a series shifted to a new level at some point (for example after a policy change) and, if so, where. It scans every admissible split of the series, reports the one with the largest difference between the before and after means, and gives it a p-value from the permutation distribution of the MAXIMUM over splits – not from the best split's own test, which is the standard way to find a change point in noise.
Arguments
- y
Numeric vector, the series in period order.
- x
Numeric or date periods, default
NULL(the position), strictly increasing and not missing whereyhas a value. Used only to label the break.- minseg
Whole number, default 2: the fewest periods allowed on either side of the break.
- nperm
Whole number, default 9999: the number of random permutations for the null distribution. The exact enumeration is used instead when \(n!\) is at most
nperm.- seed
Integer, default 1: the seed for the permutations, so the p-value is reproducible.
Value
A list of class rmbl_step_change with breakafter
(the period the series changes
after, from x), index (its position), before and
after (the segment means), difference (after minus
before), statistic (\(T_{\max}\)), p_value,
nperm (orderings or permutations used) and method.
Details
For a split after period \(k\), with \(n_a = k\) and
\(n_b = n - k\) periods and means \(\bar y_a\), \(\bar y_b\), the
statistic is the between-segment sum of squares
$$T_k = \frac{n_a n_b}{n} (\bar y_a - \bar y_b)^2,$$
over \(k = m, \dots, n - m\) with \(m\) = minseg. The
observed \(\max_k T_k\) is compared with the same maximum over
reorderings of the series, and \(p = (1 + \mathrm{count}(T^*_{\max} \ge
T_{\max}))/(1 + B)\) over \(B\) random reorderings. When \(n! \le\)
nperm (with the default, \(n \le 7\)) every ordering is
enumerated, the observed one among them, and the p-value is exact:
\(\mathrm{count}(T^*_{\max} \ge T_{\max}) / n!\);
otherwise nperm random permutations are drawn under seed,
and the caller's random number stream is restored afterwards.
The test assumes the periods are exchangeable under no change, so it ignores
autocorrelation and smooth trends; a steady trend can produce a
"step". With very few periods even a clean step cannot reach small
p-values (six periods allow only 720 orderings). Missing values are dropped.
A series with values beyond 1e150 is scanned in units of its largest absolute
value, so the statistic does not overflow; before, after and
difference are given back in the units of y, and
statistic stays in the rescaled units. Errors: infinite values; fewer
than 2 * minseg periods; numeric or date x that is not strictly
increasing. A progress bar is shown in an interactive session.
References
The permutation distribution of the maximum over splits, rather than the chosen split's own test, is what makes this a test of whether there is a break rather than a way of locating the largest wobble. See any treatment of the change-point problem, e.g. Coles, S. An Introduction to Statistical Modeling of Extreme Values (Springer), which discusses change-point detection alongside the threshold choices that raise the same multiple-comparison issue.
Examples
# A step down after the third period. With six periods there are only
# 720 orderings, so even a clean step has p near 0.1.
stepchg(c(100, 104, 98, 60, 63, 58), x = 2018:2023)
#> ╭── Step change ───────────────────────────────────────────────────────╮
#> │ Method exact over all 720 orderings │
#> │ Permutations 720 │
#> │ Break after 2020 │
#> │ Break at index 3 │
#> │ │
#> │ before after difference statistic p │
#> │ ------ ----- ---------- --------- ----- │
#> │ 100.7 60.33 -40.33 2,440 0.100 │
#> │ │
#> │ The best single split is after 2020: the mean moves from 100.7 to │
#> │ 60.33. The p-value compares that best split with the best split of │
#> │ shuffled series, so it allows for having searched every split (p = │
#> │ 0.100). │
#> ╰──────────────────────────────────────────────────────────────────────╯
# A longer series with the same step is clearly significant.
stepchg(c(100, 104, 98, 101, 97, 60, 63, 58, 61, 59), nperm = 999)$p_value
#> [1] 0.014
# Pure noise: the best split is still found, and is not significant.
set.seed(2)
stepchg(stats::rnorm(12), nperm = 999)$p_value
#> [1] 0.603
try(stepchg(c(1, 2, 3)))
#> Error : need at least 4 periods for a break with 2 either side