Detecting Multiple Step Changes Using Adaptive Regression Splines with Application to Neural Recordings

02/10/2018
by   Hazem Toutounji, et al.
ETH Zurich
0

Time series produced by dynamical systems as frequently the case in neuroscience are rarely stationary but often exhibit quite abrupt changes due to bifurcations or other dynamical phenomena. A plethora of methods for detecting such changes in time series statistics, commonly called change point analysis, have been developed over the years, in addition to test criteria to evaluate change point significance. Issues to consider when developing such methods include computational demands, difficulties arising from either limited amount of data or a large number of covariates, and arriving at statistical tests with sufficient power to detect as many change points as contained in potentially high-dimensional data sets. Here, a general method called Paired Adaptive Regressors for Cumulative Sum (PARCS) is developed for detecting multiple change points in multivariate time series. The method's flexibility to incorporate useful features from other change point detection techniques is highlighted. The advantages of PARCS over existing approaches are demonstrated through a series of simulation experiments, followed by a real data application to neural recordings from rat medial prefrontal cortex during learning.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 3

page 4

05/24/2021

Change Point Detection in Nonstationary Sub-Hourly Wind Time Series

In this paper, we present a change point detection method for detecting ...
10/08/2021

Subspace Change-Point Detection via Low-Rank Matrix Factorisation

Multivariate time series can often have a large number of dimensions, wh...
06/01/2020

Tonal harmony, the topology of dynamical score networks and the Chinese postman problem

We introduce the concept of dynamical score networks for the representat...
01/12/2021

Change-point detection using spectral PCA for multivariate time series

We propose a two-stage approach Spec PC-CP to identify change points in ...
11/06/2020

A computationally efficient, high-dimensional multiple changepoint procedure with application to global terrorism incidence

Detecting changepoints in datasets with many variates is a data science ...
01/15/2020

Detecting Changes in the Second Moment Structure of High-Dimensional Sensor-Type Data in a K-Sample Setting

The K sample problem for high-dimensional vector time series is studied,...
05/23/2021

Multiple Change Point Detection in Structured VAR Models: the VARDetect R Package

Vector Auto-Regressive (VAR) models capture lead-lag temporal dynamics o...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

Abstract

Time series, as frequently the case in neuroscience, are rarely stationary, but often exhibit abrupt changes due to attractor transitions or bifurcations in the dynamical systems producing them. A plethora of methods for detecting such change points in time series statistics have been developed over the years, in addition to test criteria to evaluate their significance. Issues to consider when developing change point analysis methods include computational demands, difficulties arising from either limited amount of data or a large number of covariates, and arriving at statistical tests with sufficient power to detect as many changes as contained in potentially high-dimensional time series. Here, a general method called Paired Adaptive Regressors for Cumulative Sum is developed for detecting multiple change points in the mean of multivariate time series. The method’s advantages over alternative approaches are demonstrated through a series of simulation experiments. This is followed by a real data application to neural recordings from rat medial prefrontal cortex during learning. Finally, the method’s flexibility to incorporate useful features from state-of-the-art change point detection techniques is discussed, along with potential drawbacks and suggestions to remedy them.

Keywords: change point, cumulative sum, adaptive regression splines, nonstationary, bootstrap test, block-permutation, behaviour, spike counts

1 Introduction

Stationary data are the exception rather than the rule in many areas of science (Aston  Kirch, 2012; Elsner ., 2004; Z. Fan ., 2015; Gärtner ., 2017; Latimer ., 2015; Paillard, 1998; Shah ., 2007; Stock  Watson, 2014). Time series statistics often change, sometimes abruptly, due to transitions in the underlying system dynamics, adaptive processes or external factors. In neuroscience, both behavioural time series (Durstewitz ., 2010; Powell  Redish, 2016; AC. Smith ., 2004) and their neural correlates (Durstewitz ., 2010; Gärtner ., 2017; Latimer ., 2015; Powell  Redish, 2016; Roitman  Shadlen, 2002) exhibit strongly nonstationary features which relate to important cognitive processes such as learning (Durstewitz ., 2010; Powell  Redish, 2016; AC. Smith ., 2004) and perceptual decision making (Hanks  Summerfield, 2017; Latimer ., 2015; Roitman  Shadlen, 2002). As such, identifying nonstationary features in behavioural and neural time series becomes necessary, both for interpreting the data in relation to the potential influences generating those features, and for removing those features from the data in order to perform statistical analyses that assume stationary observations (Hamilton, 1994; Shumway  Stoffer, 2010). Abrupt jumps in time series statistics form one important class of nonstationary events. These are often caused by bifurcations, which, in turn, may occur with gradual changes in parameters of the underlying system (Strogatz, 2001). Consequently, they are of wide interest to both statistical data analysis and the study of dynamical systems, and are commonly referred to as change points (CP; Chen  Gupta, 2012).

Detecting CPs has a long and varied history in statistics, and we will not attempt to exhaustively survey the different approaches, including regression models (Brown ., 1975; Quandt, 1958), Bayesian techniques (Chernoff  Zacks, 1964) and cumulative sum (CUSUM) statistics (Basseville, 1988; Page, 1954), to name but a few, within the limited scope of this article. Instead, we refer the reader to the excellent reviews on the topic (Aminikhanghahi  Cook, 2017; Bhattacharya, 1994; Chen  Gupta, 2012) and focus on the offline CUSUM class of methods (Hinkley, 19711) to which PARCS belongs (as opposed to sequential CUSUM methods, Page, 1954, that locate a CP online, while the time series is evolving), specifically methods that aim at detecting CPs in the mean of the time series. CUSUM-based methods are powerful, easy to implement, and are backed up by an extensive literature, theoretical results and various extensions to multiple CPs and multivariate scenarios, making them an ideal starting point. These methods assume that the time series is piecewise stationary in the statistic under consideration (e.g., piecewise constant mean) and rely on a cumulative sum transformation of the time series. Commonly, at-most-one-change (AMOC) is identified by maximum-type statistics (Kirch, 2007) at the extremum of the curve resulting from that transformation (Antoch ., 1995; Basseville, 1988).

Extending the CUSUM method to multiple CPs usually involves repetitive partitioning of the time series upon each detection (binary segmentation methods; Bai, 1997; Cho  Fryzlewicz, 2015; Fryzlewicz, 2014; Olshen ., 2004; Scott  Knott, 1974). This segmentation procedure, however, may hamper detection in later iterations as the reduction in number of observations depletes statistical power exponentially fast as more CPs are to be retrieved. In this article, we develop the PARCS (Paired Adaptive Regressors for Cumulative Sum) method which offers a straightforward extension that leverages the full time series in order to detect multiple CPs, thus providing a new solution to this issue. PARCS rests on the fact that a CUSUM transformation of the data relates to computing an integral transformation of the piecewise constant mean time series model, resulting in a piecewise linear mean function that bends at potential CPs and could be approximated by adaptive regression spline methods (Friedman, 1991; Friedman  Silverman, 1989; Stone ., 1997). Namely, rather than attempting to approximate the discontinuous time series mean directly (Efron ., 2004; Vert  Bleakley, 2010)

, the PARCS model is an approximation to the continuous CUSUM-transformed time series by a piecewise linear function. The bending points of the PARCS model are each defined by a pair of non-overlapping piecewise linear regression splines that are first selected by a two-stage iterative procedure.

The PARCS model is further refined by a nonparametric CP significance test based on bootstraps (Antoch  Hušková, 2001; Dumbgen, 1991; Hušková, 2004; Kirch, 2007; Matteson  James, 2014). While analytically derived parametric tests may usually be preferable over bootstrap-based tests due to better convergence and coverage of the tails, in the current CP setting closed form expressions for parametric tests are hard to come by and are usually replaced by approximations (Gombay  Horváth, 1996; Horváth, 1997)

. In this case, tests based on bootstraps are preferable since they are known to converge faster to the limit distribution of the test statistic

(often they are also not as conservative as parametric approximations for datasets of a relatively small size; Antoch  Hušková, 2001; Csörgö  Horváth, 1997; Kirch, 2007). In order to accommodate the possibility of temporally dependent noise in the data (Antoch ., 1997; Horváth, 1997; Picard, 1985), model selection is carried out by a nonparametric block-permutation bootstrap procedure (Davison  Hinkley, 1997; Hušková  Slabỳ, 2001; Kirch, 2007)

developed specifically for PARCS, which relies on a test statistic that quantifies the amount of bending at each candidate CP. Since model estimation is based on linear regression, PARCS is also effortlessly extended to spatially independent,

multivariate time series.

The article is structured as follows. Section 2.1 introduces the CUSUM method for AMOC detection. We then develop the PARCS method, presenting in Section 2.2 the procedure for inferring a nested model that allows for significance testing of multiple CPs, followed in Section 2.3 by an outline of the nonparametric permutation test procedure for refining the PARCS model further. Results in Section 3 illustrate that PARCS improves on several issues inherent in classical methods for change point analysis. In Section 3.1, we compare the PARCS approach to the CUSUM method in detecting a single CP, followed in Section 3.2 by a comparison with standard binary segmentation in detecting multiple CPs. We also demonstrate in Section3.3 that PARCS is successful in detecting CPs in spatially independent, multivariate time series. We then present in Section 3.4 an example from the neurosciences, in which neural and behavioural CPs are compared during operant rule-switching learning (Durstewitz ., 2010). Finally, we discuss in Section 4 the PARCS approach in relation to other state-of-the-art CP detection methods, along with drawbacks and potential extensions.

2 Methods

This section outlines the CUSUM method and the PARCS extension to multiple CPs, in addition to a nonparametric permutation technique to test for the statistical significance of CPs as identified by PARCS. For generality, the formulation assumes temporally dependent observations in the time series, independent observations being a special case.

2.1 CUSUM: Cumulative Sum of Differences to the Mean

A class of methods for identifying a single CP in the mean relies on computing a CUSUM transformation of the time series . A useful formulation that allows for dependent observations in the time series is given by the moving average (MA) step model (Antoch ., 1997; Horváth, 1997; Lombard  Hart, 1994; Kirch, 2007),

(1)

where a jump in the time series mean from baseline to occurs after time step , the change point. The step parameter or weight is positive (negative) when the time series mean increases (decreases) following . The largest integer such that noise coefficient defines a finite order of the MA process, which is 0 for temporally independent observations. We will assume that the MA process is stationary, which will always be the case if it is finite, with

independent and identically distributed (i.i.d.) random variables

(for an infinite process, points for may be considered unobserved, and coefficients have to fulfil certain conditions to make the process stationary, as given, for instance, in Shumway  Stoffer, 2010)

. The Gaussian noise assumption in the MA process can be relaxed, as long as the noise process has zero mean and finite, constant variance

(see Antoch ., 1997; Horváth, 1997; Lombard  Hart, 1994; Kirch, 2007, for theoretical results on the more general form of dependent noise). The discrete Heaviside step function, , is defined by,

Identifying the presence of a CP requires testing the null hypothesis,

, against the alternative, (Antoch ., 1995; Lombard  Hart, 1994). This begins by inferring the time of the step according to a CP locator statistic. A typical offline CP locator statistic is the maximum point of the weighted absolute cumulative sum of differences to the mean (Antoch  Hušková, 2001; Horváth, 1997),

(2)

where is the arithmetic mean of the time series (see Figure 1A). The first term on the right-hand side corrects for bias toward the centre, where more centrally-located points are down-weighed by an amount controlled by parameter . Other CUSUM-based locator statistics exist with different bias-correcting terms and cumulative sum transformations (Antoch ., 1997; Bhattacharya, 1994; Jirak, 2012; Kirch, 2007). As outlined in the Discussion, PARCS may be modified to include such bias-correcting terms as well. However, as we will demonstrate, PARCS can significantly reduce centre bias even without recourse to such a term. To show this, we will mostly deal with the generic case, , when comparing PARCS to the CUSUM transformation as defined in Eq. 2. This has the added advantage of avoiding having to select an optimal power or an optimal weight factor, a choice that usually depends on prior assumptions on the CP’s potential location (Bhattacharya, 1994). As such, and unless stated otherwise, the term CUSUM transformation will refer, thereof, to the cumulative sum of differences to the mean,

(3)

where the maximum value,

(4)

defines a test statistic by which it is decided whether to reject the null hypothesis.

Figure 1: Paired Adaptive Regressors for Cumulative Sum; (A,B) time series with (A) one or (B) two step changes and their corresponding CUSUM transformation ; (C) fitting by a piecewise linear model using two pairs of regressors and ; (D) the PARCS model fit to the CUSUM transformation of a time series , returning estimates of multiple CPs, and .

Given potentially dependent observations, , as defined by the model in Eq. 1, nonparametric bootstrap testing proceeds by block-permutation (Davison  Hinkley, 1997; Hušková, 2004; Kirch, 2007), such that temporal dependence in the data is preserved (see Section 3.2). The candidate CP is identified according to Eq. 2 and its associated test statistic is computed by Eq. 4. Estimates and are retrieved from the arithmetic means of before and after using the model in Eq. 1. By subtracting from the time series we arrive at a time series that provides an estimate of the null distribution. The stationary time series is split into blocks of size , chosen such that temporal dependencies are mostly preserved in the permuted time series (Davison  Hinkley, 1997). One way to do so is to select the block size to be larger than the order of the underlying MA process, (since the autocorrelation function of an MA process cuts off at order ; Davison  Hinkley, 1997). This requires identifying the order which can be determined from the -conform time series by inspecting its autocorrelation function (J. Fan  Yao, 2003) for different time lags . The autocorrelation function’s asymptotic distribution (Kendall  Stuart, 1983),

(5)

provides a test statistic for deciding the largest time lag at which to reject the null hypothesis , given some preset significance level . The resulting blocks are randomly permuted and each permutation is CUSUM-transformed according to Eq. 3 to compute an -conform sample of the test statistic in Eq. 4 (note that we do not know the true step parameter or the true CP , of course, such that this procedure will only yield an estimate of the distribution). A sufficiently large number of permutations results in samples of an -conform empirical distribution function (EDF) that weighs every sample equally. The candidate CP is detected when the test statistic as computed from the original time series satisfies , where is a preset significance level and the inverse of the EDF, defined as the largest value out of permutations (Davison  Hinkley, 1997; Durstewitz, 2017).

2.2 PARCS: Paired Adaptive Regressors for Cumulative Sum

The PARCS method for estimating multiple CPs rests on the fact that the integral of a piecewise constant function is piecewise linear. The AMOC model as defined in Eq. 1 assumes a piecewise stationary MA process, consisting of two segments with constant mean. A process consisting of segments generalises Eq. 1 to data containing at-most--change,

(6)

The CUSUM transformation of this process as given by Eq. 3 corresponds to the numerical integration of a piecewise stationary process . That is, is approximately (due to the noise) piecewise linear (exactly piecewise linear in the mean; see Figure 1B). If points at which bends were known, the latter can be fitted by a weighted sum of local piecewise linear basis functions or splines, centred at the knots ,

This fit corresponds to modelling the expected value of , conditioned on spline pair set , resulting in model inference,

which is a simple regression problem that can be solved by estimating the intercept and coefficients that minimise the mean-square-error,

(7)

However, in the multiple CP detection setting (assuming is known), optimal knot placement is not known a priori, but can be inferred by adaptively adjusting knot locations (Friedman, 1991; Friedman  Silverman, 1989; Stone ., 1997) to maximally satisfy the goodness-of-fit criterion in Eq. 7. In other words, and as shown in Figures 1C,D, the problem of identifying multiple CPs is replaced by the equivalent problem of inferring the order- PARCS model (or PARCS model),

(8)

with associated -tuple that best fits the CUSUM transformation of the time series. Regression coefficients in model 8 are real numbers, while knots in the present time series context are positive integers, excluding the first and last time steps, .

Fitting the PARCS model is based on a forward/backward spline selection strategy (PL. Smith, 1982) with added CP ranking stage and proceeds as outlined in Algorithm 1. Starting with an empty PARCS model, containing only the intercept , a forward sweep increases model complexity to a forward upper bound order by adding at each iteration the spline pair

, not yet contained in the model, that decreases residual mean-square-error the most. A reasonable heuristic for setting

is 2 to 3 times (assuming is known or given some liberal guess). This is followed by a backward pruning iteration, in which the spline pair whose removal increases residual mean-square-error the least is dropped from the model. Pruning removes those knots that were added at the beginning of the forward phase which became redundant as the model was refined by later additions (Friedman  Silverman, 1989). This stage continues until the number of knots reaches the preset final upper bound of model complexity , i.e., knots are pruned. Knots are then sorted in descending order according to the amount of explained variance. The ranking iteration returns a nested model by pruning the PARCS model further, down to the PARCS model. The first knot to be pruned, reducing the number of knots to , explains the least variance and is placed last as in the -tuple . The last knot to be pruned explains the most variance and is placed first as . Note that regression coefficients are re-estimated every time a knot is added to or removed from the PARCS model.

Input:
Output:
for  to  do // forward stage
     
      for  to  do // pruning stage
          
           for  to  do // ranking stage
               
Algorithm 1 Procedure for inferring the PARCS model with forward/backward spline selection (first/second loop) and CP ranking (third loop). Regression coefficients are computed by least squares estimation, conditioned on the set of knot locations of predefined size that minimises mean-square-error. Final knot locations are specified by eliminating spurious knots through block-permutation bootstrapping as described in Section 2.3.

The model can be effortlessly extended to the multiple response setting in the case of spatially independent time series (extension to a nondiagonal MA covariance matrix, Stone ., 1997, will be considered elsewhere). Given independent, piecewise stationary MA processes with common CPs ,

(9)

where , the corresponding multivariate CUSUM transformation is fitted by the multiple response, PARCS model, conditioned on common spline pairs,

using Algorithm 1. Returning CPs that are common to all variables is done by using the goodness-of-fit criterion in Eq. 7, averaged over all responses .

2.3 PARCS Model Selection by Block-Permutation Bootstrap

The piecewise linear PARCS formulation, Eq. 8, of the CUSUM transformation in Eq. 3 bends at the CPs. Due to the presence of noise in the original time series , some noise realisations may appear as slight bends in the CUSUM-transformed time series, leading PARCS to return false CPs. As such, the amount of bending at knot can be used as a test statistic for bootstrap significance testing that can refine the PARCS model further. No bending indicates either a constant fit, , or a smooth linear fit, (see also Figure 1C). Thus, a suitable test statistic that quantifies the amount of bending at is given by,

(10)

where for multivariate time series, the test statistic is the average over all time series.

Before describing the block-permutation bootstrap method for PARCS, we outline a procedure for identifying the order of the MA noise process, provided as pseudocode in Algorithm 2. First, an -conform time series is computed by regressing out the PARCS model of Eq. 8 from the CUSUM-transformed time series and then inverting the CUSUM transformation. This is followed by inspecting the autocorrelation function of for different time lags . The largest time lag at which the null hypothesis is rejected, given some preset significance level is then returned as the order , given some predefined upper bound of MA order, .

Input:
Output:
for  to  do
      if  then
           break
          
Algorithm 2 Identifying the order of the MA process, given some upper bound . The -conform time series is estimated before entering the loop. The loop increases the autocorrelation time lag and exits when the autocorrelation of is not significantly different from 0 anymore.

Given the -tuple CP set returned by Algorithm 1 and an estimate of the dependent normal noise order by Algorithm 2, a block-permutation bootstrap test returns the subset of significant CPs, as outlined in Algorithm 3. First, an -conform time series is computed. For each CP , starting with the one ranked highest, all CP-splines already deemed significant by the bootstrap test are regressed out of . A PARCS model with the remaining knots, including , is estimated and the test statistic , evaluated at according to Eq. 10, is computed. Knot is tested for significance against an -conform EDF, estimated through block-permutation bootstrapping: A total of bootstrap samples is generated from the -conform series by randomly permuting blocks of size . For each of these bootstrap samples test statistic is evaluated at knot location , yielding an EDF

which assigns equal probability

to each bootstrapped . A significant is then added to , or rejected as false discovery otherwise. The procedure repeats for the CP next in the rank order. Similar to Algorithm 1, regression coefficients are re-estimated every time a knot is added to or removed from the PARCS model.

Input:
Output:
for  to  do
      if  then
          
          
Algorithm 3 Block-permutation bootstrap procedure for PARCS, given block size . The -conform time series is estimated before entering the loop. The loop iterates over the rank-ordered CPs to test for each CP’s significance.

3 Results

We first evaluate the PARCS method on synthetic data in single and multiple CP detection settings, followed by a real data example on detecting behavioural and neural change points during rule learning.

3.1 Alleviating CUSUM Bias in AMOC Detection

We first compare the CUSUM method for detecting a single CP to the PARCS approach in order to evaluate the effect of each method on the centre bias in CP detection. Both white and MA Gaussian noise are considered. We also compare PARCS to the CUSUM locator statistic of Eq. 2 with (the maximum likelihood estimator of CP location under the assumption of i.i.d. Gaussian noise) and identify conditions under which one method is preferable over the other.

Univariate time series of length are simulated according to the step model in Eq. 1 with different levels of white Gaussian noise, , and different ground truth CP locations, . Baseline is set to and step parameter to . Note that in the step model with white Gaussian noise, increasing is equivalent to reducing . A single CP was identified by using the CUSUM method and estimating the PARCS model, both followed by bootstrap significance testing with permutations, nominal significance level , and blocks of size (since noise is independent in this example). Each parameter configuration was repeated for 1000 noise realisations.

We compare bias in CP detection toward the centre of the time series in both the CUSUM and PARCS methods. We measure this centre bias by , which is positive when estimate falls onto the side located toward the centre from , and is negative otherwise. As expected given the choice in Eq. 2, the CUSUM method shows a strong centre bias which increases for lower signal-to-noise ratio and more peripheral CPs (see Figure 2A). The CUSUM method’s power decreases for harder parameter settings (higher and more peripheral ) in that true CPs are missed in more of the realisations. The PARCS method results in a significant reduction in centre bias but does not eliminate it completely, and yields more misses relative to CUSUM if both are run at the same nominal level (see Figure 2B). Summary comparison between the two methods is shown in Figure 2C for two exemplary CP locations, .

Figure 2: Centre bias in PARCS compared to CUSUM for temporally independent noise; (A,B) bias,

, colour-coded as indicated by the colour bar; numbers indicate rounded type II error rates; (

C) bias s.e.m. for (solid) and (dashed); (D) centre bias distributions for and

; inset shows centre bias distributions as boxplots that mark the median and first and third quartiles; whiskers include points within 1.5 times the interquartile range; outliers are excluded; (

E-H) P-P plots comparing nominal (-axis) versus factual (-axis) true rejection rates in time series of length (E) , (F) , (G) , and (H) ; dotted vertical line, nominal ; dotted horizontal line, factual ; (I-L

) ROC curves depicting false discovery rate (type I error rate;

-axis) versus power (-axis) for different series lengths as in E-H; dotted vertical line, nominal ; In E-L, larger filled circles indicate the empirical rejection rates at a nominal , and empty circles indicate where the factual .

To fully appreciate the source of CUSUM centre bias and its reduction by PARCS, time series realisations with the two hardest parameter settings ( and ) are considered in Figure 2D, which compares the distribution of in the

of realisations in which both CUSUM and PARCS returned a CP. The histograms show a strongly skewed, heavy-tailed distribution for CUSUM, compared to a more symmetric distribution around 0 for PARCS, indicating only little bias. Most of the centre bias in PARCS is accounted for by outliers. This is illustrated by excluding outliers in the boxplots, which show a median of 1 time step centre bias in PARCS against median centre bias of 4 time steps in the case of CUSUM. Note that measuring centre bias as defined above does not differentiate between biased detections and false discoveries where, in extreme cases, a CP may be detected beyond the middle point

of the time series, corresponding to centre bias greater than . However, this scenario rarely occurred in the simulation results reported here.

While PARCS reduces centre bias, the simulation results above indicate that it behaves more conservatively than CUSUM at the same nominal level. In principle, false rejection rates (type II errors) may be reduced by adjusting the level, at the same time producing more false discoveries. In order to assess how well the nominal significance level agrees with the empirical type I error rate (false discoveries), 1000 white Gaussian noise realisations of length are simulated with . Conclusions drawn from this analysis are largely the same for larger signal-to-noise ratio (results not shown). A single CP was extracted using the CUSUM method and estimating the PARCS model on these time series conforming to the null hypothesis . Type I error rates at different nominal levels are shown in Figure 2E as probability-probability (P-P) plots, depicting the nominal, , against the empirical, , probabilities of accepting the null hypothesis when the null hypothesis is true. While the empirical type I error rate of CUSUM perfectly agrees with the nominal significance level, for PARCS, in contrast, the empirical rate of false discoveries tends toward 0% for and remains smaller than 1% for as large as . This entails that PARCS behaves highly conservatively, and that the nominal level may be adjusted considerably upward without strongly influencing the false discovery rate. On the other hand, despite being more conservative, Figure 2I shows that the receiver operating characteristic (ROC) curve for PARCS, depicting the method’s false discovery rate against its power for different nominal levels, consistently lies above that of CUSUM. For estimating the statistical power of each method, 1000 white Gaussian noise realisations with one CP at a random location in the range (and ) are simulated and type II error rates at different nominal levels are computed. This ROC analysis indicates that for every nominal level for CUSUM, there exists at least one nominal for PARCS such that PARCS has both higher power (fewer type II errors) and lower false detection rate (fewer type I errors), making it the preferable method. This point is explored in more detail later in the context of multiple CP detection in Section 3.2. For shorter time series, PARCS behaves similarly as for longer time series in our simulations, while for CUSUM type I error rates now start to fall below the nominal level as well (Figures 2F-H). The area under the ROC curves become smaller for shorter time series for both methods, but the ROC curve of PARCS consistently lies above that of CUSUM in those cases as well (Figures 2J-L).

Next, we examine the behaviour of the tests with dependent noise. In the case of temporally dependent noise, an appropriate block size for the bootstrap procedure could be determined by inspecting the autocorrelation function of (see Eq. 5 and Algorithm 2). 1000 noise realisations of length are drawn from an order-2 MA process with coefficients and . Since increasing noise variance in the temporally dependent case is not equivalent to decreasing the step parameter, we repeat the analysis with the same parameters from the temporally independent case above but varying the step parameter, and considering two levels of Gaussian noise, .

Figure 3: Centre bias in PARCS compared to CUSUM for temporally dependent noise with (top) and (bottom); (A,B) bias, , colour-coded as indicated by the colour bar; numbers indicate rounded rounded type II error rates; (C) bias s.e.m. for (solid) and (dashed); (D) centre bias distributions for and ; inset shows centre bias distributions as boxplots that mark the median and first and third quartiles; whiskers include points within 1.5 times the interquartile range; outliers are excluded.

Figure 3 shows results of the comparison for (top row) and

(bottom row). Similar to the white noise case, the CUSUM method’s centre bias increases for smaller signal-to-noise ratio (smaller

), and PARCS, in comparison, consistently reduces centre bias. For the same nominal level, PARCS misses more of the true CPs than CUSUM for peripheral CPs when (top panels of Figures 3A,B), but the type II error rate of the two methods is more comparable in the high noise case, (bottom panels of Figures 3A,B), despite the PARCS method’s more conservative behaviour (far lower type I error rates) in this setting as well. A summary comparison between the two methods is shown in Figure 3C for two exemplary CP locations, . Despite a significant reduction when using PARCS, both centre bias distributions in the 62% of realisations with a CP identified by the two methods, and with the two hardest parameter settings (, and ) remain strongly skewed (bottom panel of Figure 3D).

So far, we compared PARCS to the CUSUM statistic with . It is intuitive when developing PARCS to choose for the CUSUM transformation in Eq. 2, since this directly corresponds to the numerical integral of the time series upon which the PARCS approach is based (but see Section 4). Besides, under certain conditions, the CUSUM method using the test statistic with is more sensitive than that with (Antoch ., 1995). However, the CUSUM statistic with returns the maximum likelihood estimator of CP location in an AMOC scenario when noise in the step model of Eq. 1

is i.i.d. and normally distributed, leading theoretically to the strongest centre bias reduction under those conditions

(Antoch  Hušková, 2001). We therefore also compare PARCS to this maximum likelihood CUSUM estimator here, henceforth referred to as CUSUM.

Figure 4: Bias s.e.m. in PARCS compared to CUSUM with time series of length (A) , (B) , and (C) ; noise is temporally independent with .

Univerariate time series of length are simulated according to the step model Eq. 1 with different ground truth CP locations, . We consider only the scenario with largest white Gaussian noise variance, , in this analysis, for which PARCS showed the largest centre bias. A single CP was identified by using the CUSUM method and estimating the PARCS model, both followed by bootstrap significance testing. Other parameters are as in the previous analyses above. CUSUM results in a significant reduction in centre bias compared to PARCS in three of the most peripheral ground truth CPs, , but does not eliminate it completely (see Figure 4A). For the same nominal level, CUSUM also shows lower type II error rates compared to PARCS for these CPs as reported in Table 1, but recall that PARCS has a far lower type I error rate than CUSUM for the same choice of nominal (cf. Figures 2E-H and Kirch, 2007). The two methods are comparable in the quality of their detections for all other ground truth CP locations, , with PARCS having a slight advantage.

In order to assess how well the two methods fare in the small sample size limit, and to characterise the convergence behaviour of the bootstrap procedure in each method, we repeat the same analysis for shorter series lengths, . Ground truth CPs are set to the same relative location within the time series as in the simulations. As summarised in Table 1, detection rates deteriorate as series length decreases, as does the bias relative to series length (where the relative location within the series with respect to the periphery is more relevant than the absolute CP location; see Figures 4B,C). Especially for , PARCS performs mostly better than CUSUM, giving higher detection rates (see Table 1) and smaller centre bias (Figure 4C) in the majority of ground truth CPs, although it is still more conservative with near 0% type I error rate (given the bootstrap resolution; see Figure 2G). As we show next, this is a particularly important advantage of PARCS over the CUSUM-based methods when detecting multiple CPs, since CUSUM-based techniques rely on dissecting the time series into smaller segments in this case, reducing sample size at each iteration.

method type II error rate
CUSUM 07 03 02 01 01 03 12
100
PARCS 17 03 01 01 01 03 16
CUSUM 30 20 16 13 16 26 37
50
PARCS 44 22 12 08 10 19 41
CUSUM 53 38 35 32 38 44 59
26
PARCS 68 41 32 24 29 37 58
Table 1: Type II error rates in PARCS compared to CUSUM for different lengths of the time series; underline, method with higher detection rate; nominal level, 0.05 for both methods.

3.2 Detecting Multiple CPs in Univariate Data

For the scenario with multiple CPs, we assess the performance of PARCS in comparison to the CUSUM method with standard binary segmentation (Bai, 1997; Scott  Knott, 1974) for univariate data with white Gaussian noise. Standard binary segmentation is known to mislocate CPs in some scenarios, but modifications to the segmentation procedure have been proposed for solving this problem (Fryzlewicz, 2014). We show that PARCS provides an alternative approach. We then discuss a fundamental practical problem in statistical testing when using segmentation methods in general that is avoided by PARCS. Through comparison with standard binary segmentation, we illustrate conditions under which using such methods becomes infeasible. We then consider temporally dependent noise in univariate time series.

The binary segmentation method (Bai, 1997; Scott  Knott, 1974) for detecting multiple CPs proceeds as follows (pseudocode can be found in Fryzlewicz, 2014): If, according to a CUSUM test criterion, a CP is detected over the full time series, the series is partitioned at . The procedure is repeated on the resulting left and right segments, potentially returning two additional CPs, and , respectively. The procedure is terminated when no more CPs are detected after subsequent partitioning. Similar to the single CP case, we use bootstrap testing in deciding the significance of a CP at each stage. In the present context, we will refer to this CUSUM-based binary segmentation method simply by ‘CUSUM’.

Three different processes with white Gaussian noise, , of the form given by Eq. 6 and of length are simulated for 1000 realisations each. A baseline and two CPs at time steps and are set in all three scenarios. Weights are set to , and (see insets in top row of Figure 5). For the present comparison, binary segmentation is terminated after at most one partitioning, as this completely suffices to compare the methods (note that this allows CUSUM to detect up to three potential CPs, with only two present in the series). Similarly, the PARCS model is estimated, and both methods use the corresponding permutation bootstrap test with and . In order to compare type I and type II error rates between the two methods, we set nominal levels to 0.05 and 0.30 for CUSUM and PARCS, respectively, which is expected to return about the same factual type I error rates of 5% for both methods in series of length , according to Figures 2E-G.

Figure 5: Comparing PARCS to CUSUM with binary segmentation for multiple CP detection; stacked histograms of correct detection rates for CUSUM’s (top) and PARCS’ (bottom) over 1000 realisations; transparent bars show candidate CPs excluded by the permutation test; dashed grey, ground truth CPs; top inset, deterministic component of time series for the respective column’s scenario; left, centre, and right panels refer to first, second, and third scenario, respectively.

A first look at Figure 5 suggests that both CUSUM and PARCS detect CPs that are close to the ground truth. A more detailed comparison with respect to type I and type II error rates and the quality of CP detections is provided in Table 2. The quality of detections using accuracy scores is defined as the correct detection rate within a range from the ground truth CP location, adjusted for type I errors by an additive term of , where is the factual (empirical) level. This way, the accuracy score is an overall performance measure that takes into account both type I and type II errors, and how far off the detected CP is from the true one.

The first scenario (left panels in Figure 5) has the hardest parameter setting, since is both more peripheral and smaller in magnitude than . CUSUM first detects the easier CP, followed by detecting at the left hand segment as , but with a lower accuracy score as confirmed in Table 2. Similarly, PARCS returns and as the first and second rank CPs, respectively, with accuracy scores higher than those of CUSUM. The relatively low accuracy scores at detecting in both methods are due to its peripheral location and small magnitude. Both type I and type II error rates are markedly lower in PARCS.

method error rate accuracy score
type I type II
CUSUM 10 / 14 / 41 08 / 02 / 01 72 / 92 / 77 91 / 78 / 45
100
PARCS 02 / 03 / 02 04 / 00 / 01 80 / 96 / 95 96 / 74 / 76
CUSUM 12 / 18 / 19 27 / 21 / 11 34 / 56 / 61 74 / 38 / 16
50
PARCS 04 / 04 / 03 13 / 02 / 05 51 / 82 / 82 85 / 52 / 52
CUSUM 12 / 22 / 18 41 / 44 / 16 28 / 44 / 71 77 / 26 / 25
26
PARCS 06 / 07 / 06 24 / 09 / 10 37 / 75 / 76 79 / 47 / 47
Table 2: Comparing PARCS to CUSUM with binary segmentation for multiple CP detection for different lengths of the time series; error rates and accuracy scores are rounded; triplet, scenarios 1 / 2 / 3; underline, method with lower error rate or higher accuracy score; nominal levels, 0.05 and 0.30 for CUSUM and PARCS, respectively.

The second and third scenarios (centre and right panels in Figure 5, respectively) are easier in terms of ground truth parameter settings, since the more peripheral CP has the larger magnitude. In the second scenario (centre panels in Figure 5), the two methods are comparable with regards to their overall accuracy scores as defined above (many detections lie outside the accuracy score range, especially for ), but PARCS has the lower type I and type II error rates. CUSUM has a higher rate of false discoveries than in the first scenario. Its first detection is the higher magnitude CP, which comes with a lower accuracy than PARCS, followed by a detection at the right hand segment with a higher accuracy than PARCS. The relatively low accuracy rates for detecting in both methods are due to its small magnitude. The third scenario (right panels in Figure 5) is an example of a setting in which standard binary segmentation may fail in correctly allocating CPs (see Fryzlewicz, 2014, for a binary segmentation approach that solves this problem). While the performance of PARCS remains about the same as in the second scenario, CUSUM’s first detection diverges from either of the two ground truth CPs in a large number of realisations (see top-right panel in Figure 5). The large type I error rate (more than three times that in the first and second scenarios) markedly reduces the accuracy scores, which are substantially lower than for PARCS for both ground truth CPs (see Table 2). Another factor behind the high type I error rates of any iterative procedure including binary segmentation is that the same CP could be detected again at later iterations when an earlier detection is slightly biased.

We now compare the impact of shorter series on the performance of binary segmentation methods and PARCS. We consider 1000 noise realisations from each of the three scenarios with . Two ground truth CPs are set to the same relative location within the time series as for . As seen in Table 2, both type I and type II error rates increase in both methods with the decrease in series length, with the exception of type I error rates for CUSUM in the third scenario. PARCS is consistently the superior method, having both higher statistical power and less false discoveries than CUSUM. While accuracy scores predictably decrease with shorter series length, comparison between the two methods remains qualitatively similar to the case in the first and third scenario, while there is a marked change in the second scenario: While CUSUM is slightly superior in accurately detecting for , PARCS progressively surpasses CUSUM in accuracy as series length decreases. This behaviour is a result of deterioration in the power of the bootstrap test statistic for shorter series. Not only does the overall sample size decrease, but CUSUM in the second scenario with is tasked after a potential first detection of with bootstrapping the CUSUM test statistic on segments as short as 4 or 5 time steps only, which is statistically infeasible, either when using bootstraps or approximate parametric tests (Cho  Fryzlewicz, 2015; Fryzlewicz, 2014; Olshen ., 2004). We stress that this is a fundamental drawback to any method that relies on partitioning, and is not specific to standard binary segmentation.

The limitations of binary segmentation methods become more obvious when noise is temporally dependent. For instance, given a time series of length with parameters as in the second scenario and an order-2 MA noise process, blocks of size are required for proper block-permutation. If CUSUM first detected accurately, , the left hand segment would be only 20 time steps long. This allows for only 7 blocks, yielding 5040 possible permutations. This number drops to 720 permutations had been detected only 2 time steps further to the left, which makes it hard to approximate the EDF of the CUSUM statistic reliably. In addition, specifying the block size first requires estimating the MA process order by approximating the -conform time series (see Algorithm 2), but potential CPs are not known a priori, due to the recursive nature of binary segmentation methods.

Given these considerations, we focus on PARCS only as we now move over to the case of detecting multiple CPs in series with temporally dependent noise. 1000 noise realisations of length are drawn from an order-2 MA process with and coefficients and . Other parameters are as in the previous analysis. In Figure 6, we illustrate distributions of correct detections for the second scenario only (with the other two scenarios qualitatively comparable to their counterparts in Figure 5), but error and accuracy rates on these scenarios are reported as well. After removing the three step changes following PARCS CP detection, the majority of residual time series (70% as shown in bottom panel in Figure 6A) had an autocorrelation that cuts off at the correct order of the ground truth MA noise process, i.e. with being the last coefficient that lies outside the 95% confidence bounds (dashed lines; top panel in Figure 6A; results for the first and third scenarios are comparable). Block-permutation bootstrap testing with nominal and is carried out on these series with blocks of size 3 and, for other time series, according to the estimated order in Figure 6A (with an upper bound of 10 on block size). Exactly two CPs are detected in more than 99.5% of realisations in all scenarios. Figure 6B shows that the distribution of correct detections for the second scenario is largely concentrated around the ground truth CPs. Accuracy scores in each scenario are, respectively, 96%, 99%, and 99% for and 99%, 89%, and 89% for . Note also the oscillation in the example realisation in Figure 6C, which results from dependent noise with a negative MA coefficient .

Figure 6: Multiple CP detection in temporally dependent data from the second scenario; (A) estimating MA order; (top) average autocorrelation over time series realisations for different time lags

; dashed grey, 95% confidence interval; (bottom) ratio of 1000 realisations with a given estimated order; (

B) stacked histograms of correct detection rates over 1000 realisations; transparent bars show candidate CPs excluded by the permutation test; dashed grey, ground truth CPs; inset, deterministic component of time series; (C) deterministic component of the time series (grey) superimposed on an exemplary time series (blue); bottom panel shows a close up over 16 data points around ; red line highlights oscillation due to the negative MA coefficient.

3.3 Detecting Multiple CPs in Multivariate Data

The PARCS method’s ability to detect multiple CPs in spatially independent, multivariate time series is demonstrated in Figure 7 on 1000 realisations of length with covariates and white Gaussian noise, . Parameters are set to , , , and . The scaling parameter controls signal-to-noise ratio in the time series and is initially set to 1.0. Given these parameter values, the two CPs are not represented in all covariates of the time series, as exemplified in Figure 7A, rendering CPs harder to detect from the averaged univariate time series (with steps differing in sign across the covariates partially cancelling each other, resulting in small weights, and , for the resulting univariate time series).

Figure 7: Multiple CP detection in spatially independent, multivariate data; (A) deterministic component of the different covariates in the time series (grey) superimposed on an exemplary time series (blue); (B) stacked histograms of correct detection rates over 1000 realisations; transparent bars show candidate CPs excluded by the permutation test; dashed grey, ground truth CPs; (C) stacked histograms of correct detection rates over 1000 realisations given different signal-to-noise ratios (controlled by ), and binned with 5 time step windows; rates are logarithmically scaled as indicated by the colour bar.

Following PARCS CP detection, augmented by a bootstrap test with nominal , and , exactly two CPs are detected in 99.9% of realisations. Accuracy scores are 99.8% and 98% for and , respectively. The lower variance in detections, as seen in Figure 7B, is due to the higher average absolute weight compared to .

We then test the method’s performance for smaller signal-to-noise ratios with . Figure 7C shows that correct detection rates within a range from the ground truth CPs for different values of remain above 50%, even for magnitudes as small as . These rates are a result of PARCS leveraging CP information from multiple covariates simultaneously, rather than depleting the signal through averaging.

3.4 Detecting Neural Events that Reflect Learning

A previous study by one of the authors and colleagues exemplifies the practical value of change point detection in neuroscience (Durstewitz ., 2010). These authors demonstrated that acquiring a new behavioural rule in rats is accompanied by sudden jumps in behavioural performance, which in turn is reflected in the activity of neural units recorded simultaneously in the medial prefrontal cortex (mPFC). In the current section, we revisit part of these data to showcase PARCS in a real data scenario.

Before moving to the demonstration, it is important to note that the data in question are not normally distributed and potentially include linear trends (Durstewitz ., 2010) not accounted for by the step models in Eqs. 1, 6 and 9. As such, some preprocessing may be necessary for a statistical analysis that is more consistent with the step model assumptions (this may include detrending and potentially some mild smoothing with Gaussian kernels; see Durstewitz ., 2010)

. However, to keep the present demonstration simple, PARCS was applied directly to the data with minimal preprocessing, which only involves square-root-transforming the neural count data for bringing them closer to a Gaussian distribution and stabilising the variance

(Kihlberg ., 1972).

In order to show that PARCS can still return reasonable CP estimates under these non-Gaussian conditions, we first test its performance on simulated spike count data, before applying it to the empirical data. We simulate 1000 realisations of length with covariates according to a Poisson process. Parameters are set to , , , and . This choice of parameters results in average firing rates that are comparable in their means to the white Gaussian noise case (cf. Figures 7A and 8

A) and to the low firing rates often observed in mPFC neurons. One obvious diversion from Gaussian assumptions in the case of Poisson noise is that the variance is not constant anymore, but is equal to the means within each of the segments separated by true CPs. Following square-root-transforming the data and PARCS

CP detection, augmented by a bootstrap test with nominal , and , exactly two CPs are detected in 92% of realisations. Accuracy scores are 98% and 70% for and , respectively. The lower accuracy scores compared to the white Gaussian noise scenario are due to the lower signal-to-noise ratios, resulting from the increase in noise variance with firing rates (cf. Figures 7B and 8B). Nevertheless, these results sufficiently justify the use of PARCS in the present context.

Figure 8: Multiple CP detection in spatially independent, multivariate, Poisson data; (A) deterministic component of the different covariates in the time series (grey) superimposed on an exemplary time series (blue); -axis, square-root-transformed spike counts; (B) stacked histograms of correct detection rates over 1000 realisations; transparent bars show candidate CPs excluded by the permutation test; dashed grey, ground truth CPs.
Figure 9: Comparing behavioural and mPFC neural CPs; (A) blue, square-root-transformed spike count data in the three seconds following cue onset from 6 representative mPFC units of one rat; grey, mean as estimated by inverting the neural multiple response PARCS model. Note potential CP in top-centre unit which was not detected by PARCS since it did not contribute strongly to population-wide CPs; dashed lines, behavioural CPs from the same animal; (B) blue, lever press at each trial; this animal is rewarded for pressing the right lever during the spatial rule; grey, probability of pressing right lever as estimated by inverting the behavioural PARCS model; dashed lines, neural CPs from the same animal (see A); (C) relating behavioural and neural CPs; blue, behavioural CPs with higher weight; , correlation coefficients as computed over all 12 data points (black) and over those where behavioural CPs have the higher weight (blue); -values, significance levels of corresponding ; black and blue lines, respective least-square linear regression fits to the two sets of data points; red and yellow circles, neural and behavioural CP pairs from the exemplary animal in A and B, respectively.

We now turn to the experimentally obtained dataset. Six animals were trained on a two-choice deterministic operant rule switching task which proceeds as follows: At the beginning of the session, the animal follows a previously acquired behavioural rule whereby it responds to a visual cue by a lever press for attaining a reward (visual rule). Unknown to the animal, reward contingencies are switched after 20 trials to a novel spatial rule, in which attaining the reward requires pressing a certain baited lever (right or left), regardless of the visual cue. The session is terminated when the animal reaches a preset criterion that indicates that the new rule behaviour has been learnt. In addition to the binary behavioural data of lever presses over trials, spike counts emitted by mPFC units during the 3 seconds following cue onset were collected through single unit recording techniques. Neural and behavioural data from one animal are shown in Figure 9A and 9B, respectively. Trials corresponding to the steady state visual and spatial rule (first and last 20 trials, respectively) are not considered in the analysis.

Time series with one or two significant CPs were described in the original study byDurstewitz . (2010), so PARCS models are estimated for each animal for both the multivariate neural data (multiple response PARCS model; Figure 9A) and the univariate behavioural data (Figure 9B), in addition to PARCS models for the behavioural data as summarised in Figure 9C. As shown in Figure 9B, one neural CP in that exemplary animal matches its behavioural counterpart. The second neural CP, while not as close to its behavioural counterpart, is only 7 trials apart, and the two are highly correlated across animals, as shown in Figure 9C, concurring with the original findings of Durstewitz . (2010). Besides the significant correlation, the corresponding black linear regression line lies very close to the diagonal, indicating that neural and behavioural CPs are not only correlated, but are almost equal. Moreover, those authors report that data from many animals contain at most a single CP (also note the low weight of one of the CPs as estimated from one of the animals using PARCS). Comparing the PARCS behavioural CP to its neural counterpart (blue circles in Figure 9C) shows that correlation remains high and significant. A sample size of is usually recommended for evaluating the significance of a correlation. However, the corresponding linear regression line is also close to the diagonal, in further support for the reliability of this result, and in agreement with the original results despite different procedures: For the neural data in the original study, CUSUM-based detection was performed on a multivariate discrimination statistic defined across the whole neural population, while here, the model was determined directly from the multiple spike count data.

4 Discussion

In the current article, we introduced PARCS, a method for detecting multiple step changes, or CPs, in potentially multivariate, temporally dependent data, supported by a bootstrap-based nonparametric test. We also showed that PARCS substantially reduces centre bias in estimating CPs compared to the most basic specification of the CUSUM method, and presented conditions under which it compares to or outperforms the maximum-likelihood CUSUM statistic. Furthermore, we demonstrated that PARCS may achieve higher sensitivity (statistical power) than CUSUM-based methods while at the same time having lower type I errors in multiple CP scenarios, mainly because PARCS can make use of the full time series while CUSUM-based methods rely on segmenting the time series for detecting multiple CPs. We finally confirmed previous results pertaining to the acquisition of a new behavioural rule and the role of the medial prefrontal cortex in this process.

As already apparent from some of our simulation studies, the basic PARCS method as introduced here leaves room for improvement. In the presence of a single CP, we showed that PARCS strongly reduces the amount of bias toward the centre that results from the direct application of the most basic form of the CUSUM locator statistic. Theoretically-grounded modifications to the CUSUM transformation that reduce this amount of bias rely on down-weighing more centrally-located points (Kirch, 2007). As shown with PARCS, this problem is not quite as severe. Nevertheless, since PARCS approximates the CUSUM transformation using a regression model, similar down-weighing could be incorporated into the PARCS procedure as well by using weighted least squares instead of regular least squares (Hastie ., 2009), which is a straightforward amendment. Furthermore, the PARCS method currently requires a liberal guess of the number of CPs in advance, followed by refinements through nonparametric bootstrap testing. It is desirable, however, especially when no prior information on is available, to have statistical tests as termination criteria for the forward and backward stages. In adaptive regression spline methods (Friedman, 1991; Friedman  Silverman, 1989; Stone ., 1997), there is strong empirical evidence (Hinkley, 1969, 19712) backed by theoretical results (Feder, 1975) that the difference in residual mean-square-error between two nested models that differ in one additional knot is well approximated, albeit conservatively, by a scaled

statistic on 4 degrees of freedom

(Friedman, 1991). This led to one nonparametric termination recipe that is based on generalised cross-validation (Craven  Wahba, 1979). Another approach is to infer the piecewise linear regression model with the aid of a parametric test for specifying the number and location of knots, without recourse to iterative procedures (Liu ., 1997). Unfortunately, neither approach is directly applicable to PARCS, since they both require assumptions that are not met in the CUSUM-transformed time series. The CUSUM transformation of the time series is a nonstationary ARMA process. Deriving reasonable generalised cross-validation (Craven  Wahba, 1979; Friedman, 1991; Friedman  Silverman, 1989), F-ratio (Durstewitz, 2017; Hastie ., 2009) or parametric (Liu ., 1997) test statistics require currently unknown corrections to those tests which account for nonstationarity and the particular form of the ARMA model underlying the CUSUM-transformed data.

When multiple CPs are present in the data, PARCS can outperform standard binary segmentation (Bai, 1997; Scott  Knott, 1974). Other segmentation methods also solve the problem of mislocating CPs inherent in the standard procedure (Cho  Fryzlewicz, 2015; Fryzlewicz, 2014; Olshen ., 2004). Wild binary segmentation (WBS; Fryzlewicz, 2014), for instance, relies on sampling local CUSUM transformations of randomly chosen segments of the time series. The candidate CP with the largest value among sampled CUSUM curves is returned to be tested against a criterion, followed by binary segmentation. WBS is preferable to PARCS in that its test statistic and termination criterion when noise is independent are backed up by rigorous theory, and may be the favourable method when segments are large enough for the test statistic to converge. If series of only limited length are available, however, WBS may run into similar problems as standard binary segmentation for CUSUM, since each detection is still followed by partitioning the data further. WBS also, to the best of our knowledge, currently lacks a thorough analysis on the behaviour of its test statistic for dependent data. It is tempting to speculate on the potential for a hybrid method that capitalises on the desirable features of both methods. Computational demands arise in WBS from the need to choose segment range parameters by sampling few thousand CUSUM curves to which PARCS may offer an easy and efficient workaround: Fryzlewicz (2014) demonstrated that the optimal WBS segment choice is the segment bounded by the two CPs closest to the target CP from each side. PARCS could thus provide an informed selection of boundaries by returning candidate CPs in the data and use these to demarcate segments, rather than random sampling as in WBS.

In dealing with multivariate data, recent methods tackled the computational demands of having a large number of covariates and sparse CP representations (Cho  Fryzlewicz, 2015; Wang  Samworth, 2018). These methods rely on low dimensional projections of the multivariate CUSUM curve that preserve the CPs and follow this projection by a binary segmentation method. Since PARCS for multivariate time series is also based on the CUSUM transformation, it is straightforward to leverage the computational savings provided by such projection methods in reducing the dimensionality of the PARCS input, while avoiding the drawbacks of binary segmentation methods. This may offer a route for extending PARCS to the important case of multivariate CP detection in mutually dependent time series with spatial dependence, a configuration which these projection methods also consider (Cho  Fryzlewicz, 2015; Wang  Samworth, 2018). Alternatively, nondiagonal covariance structure in multivariate series may be accounted for by extending the PARCS formulation to the multivariate regression spline realm (Friedman, 1991; Stone ., 1997).

Finally, when analysing the neural and behavioural data during the rule switching task, we mentioned that data may also contain trends that are not accounted for by step change time series models (Durstewitz ., 2010). Caution must be made when analysing real data using CP detection methods in that these methods, PARCS included, assume a step change model underlying the generation of the data and hence may attempt to approximate trends and other nonstationary features by a series of step changes, a point made more explicit by Fryzlewicz (2014) (Durstewitz ., 2010, therefore removed trends around candidate CPs first). Hence, to avoid wrong conclusions with respect to the source and type of nonstationarity in experimental time series, it may be necessary to either augment change point detection by adequate preprocessing (Durstewitz ., 2010) or to generalise time series models for CP detection to include other forms of nonstationarity.

Author Contributions

HT, DD conceived study; HT developed and implemented methods; HT carried out simulations; HT analysed data; HT, DD interpreted results; HT prepared figures; HT wrote manuscript; HT, DD revised and finalised manuscript.

Funding

This research was funded by grants to DD by the German Research Foundation (DFG) (SPP1665, DU 354/8-2) and through the German Ministry for Education and Research (BMBF) via the e:Med framework (01ZX1311A & 01ZX1314E).

Acknowledgements

The authors thank Dr. Georgia Koppe and Dr. Eleonora Russo for discussions and Dr. Jeremy Seamans for providing the neural and behavioural data.

Data Availability Statement

Method implementation will be freely available on an online repository upon publication.

References

  • Aminikhanghahi  Cook (2017) Aminikhanghahi2017Aminikhanghahi, S.  Cook, DJ.  2017. A survey of methods for time series change point detection A survey of methods for time series change point detection. Knowl. Inf. Syst.512339–367. 10.1007/s10115-016-0987-z
  • Antoch . (1997) Antoch1997Antoch, J., Hušková, M.  Prášková, Z.  1997. Effect of dependence on statistics for determination of change Effect of dependence on statistics for determination of change. J. Stat. Plan. Inference602291–310. 10.1016/S0378-3758(96)00138-3
  • Antoch . (1995) Antoch1995Antoch, J., Hušková, M.  Veraverbeke, N.  1995. Change-point problem and bootstrap Change-point problem and bootstrap. J. Nonparametr. Statist.52123–144. 10.1080/10485259508832639
  • Antoch  Hušková (2001) Antoch2001Antoch, J.  Hušková, M.  2001. Permutation tests in change point analysis Permutation tests in change point analysis. Stat. Probab. Lett.53137–46. 10.1016/S0167-7152(01)00009-8
  • Aston  Kirch (2012) Aston2012Aston, JAD.  Kirch, C.  2012. Evaluating stationarity via change-point alternatives with applications to fMRI data Evaluating stationarity via change-point alternatives with applications to fMRI data. Ann. Appl. Stat.641906–1948. 10.1214/12-AOAS565
  • Bai (1997) Bai1997Bai, J.  1997. Estimating multiple breaks one at a time Estimating multiple breaks one at a time. Econ. Theory133315–352. 10.1017/S0266466600005831
  • Basseville (1988) Basseville1988Basseville, M.  1988. Detecting changes in signals and systems–a survey Detecting changes in signals and systems–a survey. Automatica243309–326. 10.1016/0005-1098(88)90073-8
  • Bhattacharya (1994) Bhattacharya1994Bhattacharya, PK.  1994. Some aspects of change-point analysis Some aspects of change-point analysis. Lect. Notes Monogr. Ser.2328–56. 10.1214/lnms/1215463112
  • Brown . (1975) Brown1975Brown, RL., Durbin, J.  Evans, JM.  1975. Techniques for testing the constancy of regression relationships over time Techniques for testing the constancy of regression relationships over time. J. R. Stat. Soc. Series B Stat. Methodol.149–192.
  • Chen  Gupta (2012) Chen2012Chen, J.  Gupta, AK.  2012. Parametric statistical change point analysis: With applications to genetics, medicine, and finance Parametric statistical change point analysis: With applications to genetics, medicine, and finance. Basel, SwitzerlandSpringer Science+Business Media, LLC.
  • Chernoff  Zacks (1964) Chernoff1964Chernoff, H.  Zacks, S.  1964. Estimating the current mean of a normal distribution which is subjected to changes in time Estimating the current mean of a normal distribution which is subjected to changes in time. Ann. Math. Stat.353999–1018. 10.1214/aoms/1177700517
  • Cho  Fryzlewicz (2015) Cho2015Cho, H.  Fryzlewicz, P.  2015. Multiple-change-point detection for high dimensional time series via sparsified binary segmentation Multiple-change-point detection for high dimensional time series via sparsified binary segmentation. J. R. Stat. Soc. Series B Stat. Methodol.772475–507. 10.1111/rssb.12079
  • Craven  Wahba (1979) Craven1978Craven, P.  Wahba, G.  1979. Smoothing Noisy Data with Spline Functions: Estimating the Correct Degree of Smoothing by the Method of Generalized Cross-Validation Smoothing noisy data with spline functions: Estimating the correct degree of smoothing by the method of generalized cross-validation. Numer. Math.314317–403. 10.1007/BF01404567
  • Csörgö  Horváth (1997) Csorgo1997Csörgö, M.  Horváth, L.  1997. Limit theorems in change-point analysis Limit theorems in change-point analysis ( 18). New York, NYJohn Wiley & Sons Inc.
  • Davison  Hinkley (1997) Davison1997Davison, AC.  Hinkley, DV.  1997. Bootstrap methods and their application Bootstrap methods and their application. New York, NYCambridge University Press.
  • Dumbgen (1991) Dumbgen1991Dumbgen, L.  1991. The asymptotic behavior of some nonparametric change-point estimators The asymptotic behavior of some nonparametric change-point estimators. Ann. Stat.1471–1495. 10.1214/aos/1176348257
  • Durstewitz (2017) Durstewitz2017Durstewitz, D.  2017. Advanced Data Analysis in Neuroscience: Integrating Statistical and Computational Models Advanced data analysis in neuroscience: Integrating statistical and computational models. Springer International Publishing.
  • Durstewitz . (2010) Durstewitz2010Durstewitz, D., Vittoz, NM., Floresco, SB.  Seamans, JK.  2010. Abrupt Transitions between Prefrontal Neural Ensemble States Accompany Behavioral Transitions during Rule Learning Abrupt transitions between prefrontal neural ensemble states accompany behavioral transitions during rule learning. Neuron663438–448. 10.1016/j.neuron.2010.03.029
  • Efron . (2004) Efron2004Efron, B., Hastie, T., Johnstone, I.  Tibshirani, R.  2004. Least angle regression Least angle regression. Ann. Stat.322407–499. 10.1214/009053604000000067
  • Elsner . (2004) Elsner2004Elsner, JB., Niu, X.  Jagger, TH.  2004.

    Detecting shifts in hurricane rates using a Markov chain Monte Carlo approach Detecting shifts in hurricane rates using a Markov chain Monte Carlo approach.

    J. Clim.17132652–2666. 10.1175/1520-0442(2004)017¡2652:DSIHRU¿2.0.CO;2
  • J. Fan  Yao (2003) Fan2003Fan, J.  Yao, Q.  2003. Nonlinear time series: nonparametric and parametric methods Nonlinear time series: nonparametric and parametric methods. New York, NYSpringer.
  • Z. Fan . (2015) Fan2015Fan, Z., Dror, RO., Mildorf, TJ., Piana, S.  Shaw, DE.  2015. Identifying localized changes in large systems: Change-point detection for biomolecular simulations Identifying localized changes in large systems: Change-point detection for biomolecular simulations. Proc. Natl. Acad. Sci. U.S.A.112247454–9. 10.1073/pnas.1415846112
  • Feder (1975) Feder1975Feder, PI.  1975. The log likelihood ratio in segmented regression The log likelihood ratio in segmented regression. Ann. Stat.84–97.
  • Friedman (1991) Friedman1991Friedman, JH.  1991. Multivariate adaptive regression splines Multivariate adaptive regression splines. Ann. Stat.1911–67. 10.1214/aos/1176347963
  • Friedman  Silverman (1989) Friedman1989Friedman, JH.  Silverman, BW.  1989. Flexible Parsimonious Smoothing and Additive Modeling Flexible parsimonious smoothing and additive modeling. Technometrics3113–21. 10.1080/00401706.1989.10488470
  • Fryzlewicz (2014) Fryzlewicz2014Fryzlewicz, P.  2014. Wild binary segmentation for multiple change-point detection Wild binary segmentation for multiple change-point detection. Ann. Stat.4262243–2281. 10.1214/14-AOS1245
  • Gärtner . (2017) Gartner2017Gärtner, M., Duvarci, S., Roeper, J.  Schneider, G.  2017. Detecting joint pausiness in parallel spike trains Detecting joint pausiness in parallel spike trains. J. Neurosci. Methods28569–81. 10.1016/j.jneumeth.2017.05.008
  • Gombay  Horváth (1996) Gombay1996Gombay, E.  Horváth, L.  1996. On the rate of approximations for maximum likelihood tests in change-point models On the rate of approximations for maximum likelihood tests in change-point models. J. Multivar. Anal.561120–152. 10.1006/jmva.1996.0007
  • Hamilton (1994) Hamilton1994Hamilton, JD.  1994. Time series analysis Time series analysis. Princeton, NJPrinceton University Press.
  • Hanks  Summerfield (2017) Hanks2017Hanks, TD.  Summerfield, C.  2017. Perceptual Decision Making in Rodents, Monkeys, and Humans Perceptual decision making in rodents, monkeys, and humans. Neuron93115–31. 10.1016/j.neuron.2016.12.003
  • Hastie . (2009) Hastie2009Hastie, T., Tibshirani, R.  Friedman, J.  2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction The elements of statistical learning: Data mining, inference, and prediction (second ). New York, NYSpringer.
  • Hinkley (1969) Hinkley1969Hinkley, DV.  1969. Inference about the intersection in two-phase regression Inference about the intersection in two-phase regression. Biometrika563495–504. 10.1093/biomet/56.3.495
  • Hinkley (19711) Hinkley1971aHinkley, DV.  19711. Inference about the change-point from cumulative sum tests Inference about the change-point from cumulative sum tests. Biometrika583509–523. 10.1093/biomet/58.3.509
  • Hinkley (19712) Hinkley1971bHinkley, DV.  19712. Inference in two-phase regression Inference in two-phase regression. J. Am. Stat. Ass.66336736–743. 10.1080/01621459.1971.10482337
  • Horváth (1997) Horvath1997Horváth, L.  1997. Detection of changes in linear sequences Detection of changes in linear sequences. Ann. Inst. Stat. Math.492271–283. 10.1023/A:1003110912735
  • Hušková (2004) Huskova2004Hušková, M.  2004. Permutation principle and bootstrap in change point analysis Permutation principle and bootstrap in change point analysis. L. Horváth  B. Szyszkowicz (), Asymptotic Methods in Stochastics: Festschrift for Miklós Csörgő Asymptotic methods in stochastics: Festschrift for Miklós Csörgő ( 44,  273–291). American Mathematical Soc.
  • Hušková  Slabỳ (2001) Huvskova2001Hušková, M.  Slabỳ, A.  2001. Permutation tests for multiple changes Permutation tests for multiple changes. Kybernetika375605–622.
  • Jirak (2012) Jirak2012Jirak, M.  2012. Change-point analysis in increasing dimension Change-point analysis in increasing dimension. J. Multivar. Anal.111136–159. 10.1016/j.jmva.2012.05.007
  • Kendall  Stuart (1983) Kendall1983Kendall, MG.  Stuart, A.  1983. The Advanced Theory of Statistics The advanced theory of statistics ( 3). London, UKGriffin.
  • Kihlberg . (1972) Kihlberg1972Kihlberg, J., Herson, J.  Schotz, W.  1972. Square root transformation revisited Square root transformation revisited. Appl. Statist.2076–81. 10.1214/aos/1176343000
  • Kirch (2007) Kirch2007Kirch, C.  2007. Block permutation principles for the change analysis of dependent data Block permutation principles for the change analysis of dependent data. J. Stat. Plan. Inference13772453–2474. 10.1016/j.jspi.2006.09.026
  • Latimer . (2015) Latimer2015Latimer, KW., Yates, JL., Meister, MLR., Huk, AC.  Pillow, JW.  2015. Single-trial spike trains in parietal cortex reveal discrete steps during decision-making Single-trial spike trains in parietal cortex reveal discrete steps during decision-making. Science3496244184–187. 10.1126/science.aaa4056
  • Liu . (1997) Liu1997Liu, J., Wu, S.  Zidek, JV.  1997. On segmented multivariate regression On segmented multivariate regression. Stat. Sin.497–525.
  • Lombard  Hart (1994) Lombard1994Lombard, F.  Hart, J.  1994. The analysis of change-point data with dependent errors The analysis of change-point data with dependent errors. Lect. Notes Monogr. Ser.23194–209. doi:10.1214/lnms/1215463125
  • Matteson  James (2014) Matteson2014Matteson, DS.  James, NA.  2014. A nonparametric approach for multiple change point analysis of multivariate data A nonparametric approach for multiple change point analysis of multivariate data. J. Am. Stat. Assoc.109505334–345. 10.1080/01621459.2013.849605
  • Olshen . (2004) Olshen2004Olshen, AB., Venkatraman, E., Lucito, R.  Wigler, M.  2004. Circular binary segmentation for the analysis of array-based DNA copy number data Circular binary segmentation for the analysis of array-based DNA copy number data. Biostatistics54557–572. 10.1093/biostatistics/kxh008
  • Page (1954) Page1954Page, E.  1954. Continuous inspection schemes Continuous inspection schemes. Biometrika411/2100–115. 10.2307/2333009
  • Paillard (1998) Paillard1998Paillard, D.  1998. The timing of Pleistocene glaciations from a simple multiple-state climate model The timing of pleistocene glaciations from a simple multiple-state climate model. Nature3916665378. 10.1038/34891
  • Picard (1985) Picard1985Picard, D.  1985. Testing and estimating change-points in time series Testing and estimating change-points in time series. Adv. Appl. Probab.174841–867. 10.1017/S0001867800015433
  • Powell  Redish (2016) Powell2016Powell, NJ.  Redish, AD.  2016. Representational changes of latent strategies in rat medial prefrontal cortex precede changes in behaviour Representational changes of latent strategies in rat medial prefrontal cortex precede changes in behaviour. Nat. Commun.712830. 10.1038/ncomms12830
  • Quandt (1958) Quandt1958Quandt, RE.  1958. The estimation of the parameters of a linear regression system obeying two separate regimes The estimation of the parameters of a linear regression system obeying two separate regimes. J. Am. Stat. Assoc.53284873–880. 10.1080/01621459.1958.10501484
  • Roitman  Shadlen (2002) Roitman2002Roitman, JD.  Shadlen, MN.  2002. Response of neurons in the lateral intraparietal area during a combined visual discrimination reaction time task Response of neurons in the lateral intraparietal area during a combined visual discrimination reaction time task. J. Neurosci.22219475–9489.
  • Scott  Knott (1974) Scott1974Scott, AJ.  Knott, M.  1974.

    A cluster analysis method for grouping means in the analysis of variance A cluster analysis method for grouping means in the analysis of variance.

    Biometrics303507–512. 10.2307/2529204
  • Shah . (2007) Shah2007Shah, SP., Lam, WL., Ng, RT.  Murphy, KP.  2007. Modeling recurrent DNA copy number alterations in array CGH data Modeling recurrent DNA copy number alterations in array CGH data. Bioinformatics2313i450–i458. 10.1093/bioinformatics/btm221
  • Shumway  Stoffer (2010) Shumway2010Shumway, RH.  Stoffer, DS.  2010. Time series analysis and its applications: With R examples Time series analysis and its applications: With R examples (3rd ). New York, NYSpringer.
  • AC. Smith . (2004) Smith2004Smith, AC., Frank, LM., Wirth, S., Yanike, M., Hu, D., Kubota, Y.Brown, EN.  2004. Dynamic Analysis of Learning in Behavioral Experiments Dynamic analysis of learning in behavioral experiments. J. Neurosci.242447–461.
  • PL. Smith (1982) Smith1982Smith, PL.  1982. Curve fitting and modeling with splines using statistical variable selection techniques Curve fitting and modeling with splines using statistical variable selection techniques NASA Report 166034. Langley Research Center, Hampton, VA.
  • Stock  Watson (2014) Stock2014Stock, JH.  Watson, MW.  2014. Estimating turning points using large data sets Estimating turning points using large data sets. J. Econom.178368–381. 10.1016/j.jeconom.2013.08.034
  • Stone . (1997) Stone1997Stone, CJ., Hansen, MH., Kooperberg, C.  Truong, YK.  1997. Polynomial splines and their tensor products in extended linear modeling: 1994 Wald memorial lecture Polynomial splines and their tensor products in extended linear modeling: 1994 Wald memorial lecture. Ann. Stat.2541371–1470. 10.1214/aos/1031594728
  • Strogatz (2001) Strogatz2001Strogatz, SH.  2001. Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering. Perseus Books Publishing.
  • Vert  Bleakley (2010) Vert2010Vert, JP.  Bleakley, K.  2010. Fast detection of multiple change-points shared by many signals using group LARS Fast detection of multiple change-points shared by many signals using group LARS. JD. Lafferty, CKI. Williams, J. Shawe-Taylor, RS. Zemel  A. Culotta (), Advances in Neural Information Processing Systems 23 Advances in neural information processing systems 23 ( 2343–2351). Curran Associates, Inc.
  • Wang  Samworth (2018) Wang2018Wang, T.  Samworth, RJ.  2018. High dimensional change point estimation via sparse projection High dimensional change point estimation via sparse projection. J. R. Stat. Soc. Series B Stat. Methodol.80157–83. 10.1111/rssb.12243