Recursion

The problem is to find break-points in 2-dimensional ordered data {(x_i, y_i), i = 1, ..., n}, eg., a time series. This can be formulated as a dynamic programming recursion:

$$ \begin{equation} S_j = \max_{i\le j} (S_{i-\jump} + \text{score}(i,j)) - P\;, (\#eq:recur) \end{equation} $$

where the scoring function score(i, j) quantifies how well a segment between points i and j is defined.

The break-point penalty term P sets bounds on segment lengths and should be chosen close to the optimal value of the scoring function (Section @ref(selectp)). Higher P will yield longer segments, and especially for non-monotonic data P lower than the optimal value for the scoring function can work better.

The binary jump parameter $\jump \in \{0,1\}$ determines whether the break-point i is both, the end of the previous and start of the current segment ($\jump=0$), or the previous segments ends at i − 1 ($\jump=1$). The latter choice allows discontinuous jumps between adjacent segments (Fig. @ref(fig:jumps)).

Scoring Functions

Segmentation into linear segments can be achieved by using a goodness-of-fit measure. Linear models are often evaluated by the coefficient of determination R² (“r-squared”), and we can use this directly as a scoring function (type="r2" in dpseg):

$$ \begin{equation} \text{score}(i,j) = R^2 - 1 = \frac{\sum (\hat y_i - \bar y)^2}{\sum (y_i - \bar y)^2}-1\;, (\#eq:scorer2) \end{equation} $$

where ŷ_i is the linear model and ȳ the mean of the original data, see Section @ref(appii) and equation @ref(eq:r2) for details. Its square root, the Pearson correlation can also be used as scoring function (type="cor"). Subtraction of 1 simply aligns different scoring functions in dpseg at score ≤ 0 and thereby allows the use of a consistent default penalty parameter of P = 0.

R² depends on the slope (eq. @ref(eq:r2)) and will score poorly for segments without slope (Fig. @ref(fig:lines) and @ref(fig:custom), Section @ref(lines)). The negative variance of the residuals r_i = y_i − ŷ_i between data and the regression line does not depend on the slope (eq. @ref(eq:var)) and is for many cases the better choice and the default in dpseg (type="var"):

$$ \begin{equation} \text{score}(i,j) = -\Var(r) = - \frac{1}{n-1} \sum r_i^2\;, (\#eq:score) \end{equation} $$

See Section @ref(scoring) for a short discussion on the choice of the built-in scoring functions. Section @ref(lines) discusses the handling of straight horizontal or vertical lines (with zero variance) in the data. And Sections @ref(matrix) and @ref(generic) introduce custom scoring options.

Incremental calculation of the variances (Section @ref(appii)) allows for a computationally efficient implementation (Section @ref(benchmark)) of the recursion in dpseg.

Backtracing

During calculation of S_j the position i_max < j that yielded the maximal score is stored as a vector. The segmentation can be easily retrieved by starting at the last position j = n, the end of the last segment, looking up its i_max, the start of the segment, and proceed from $j=i_\text{max} - \jump$ (the end of the next segment). If jumps were allowed ($\jump=1$) the previous segment ended one position ahead of the next segment at i − 1, and this must be accounted for in backtracing.

simple_backtrace <- function (imax, jumps = 0) {
    end <- length(imax) # end of the first segment = end of the data
    ends <- c(end)      # initiate vector
    while (end > 1) {
        end <- imax[end] - jumps # end of previous segment
        ends <- c(end, ends)     # note the order, new end is prepended
    }
    ends
}

Segment Length Restrictions

Depending on the scoring function, short segments may not be meaningful, eg. goodness-of-fit measures for linear regression such as the variance of residuals are 0 for segments of length ℓ = 2 and an “optimal” segmentation would split the data in pairs with perfect “fits”. Thus, we can restrict score function maximization for only $i\le(j-\lmin+1)$.

When handling data sets much larger then the expected segment length, one can also define a maximal segment length, and thereby save memory and computation time by only considering $i\ge(j-\lmax +1)$.

Combining both restrictions, and with $\Ell_x=\ell_x+1$ the recursion becomes:

$$ \begin{equation} S_j = \max_{\substack{i\le j-\Ell_{\min}\\i\ge j-\Ell_{\max}}} (S_{i-\jump} + \text{score}(i,j)) - P\;. (\#eq:recurl) \end{equation} $$

Basic

The main function dpseg takes x and y vectors of equal length as input and returns an object of class dpseg. This results object can be inspected by print and plot methods. A predict method returns y values and specified novel x values at the straight line segments.

The results object is a simple list (S3 object) with the list item segments as the main result: a table (data frame) that lists the start and end x-values of the segments, the start and end indices in the data vectors, the linear regression coefficients and goodness-of-fit measures for the segments (intercept, slope, r-squared, variance of residuals).

library(dpseg)

type <- "var"  # use the (default) scoring, -var(residuals(lm(y~x)))
jumps <- FALSE # allow discrete jumps between segments?
P <- 1e-4      # break-point penalty, use higher P for longer segments

# get example data `oddata` - bacterial growth measured as optical density OD
x <- oddata$Time
y <- log(oddata[,"A2"]) # note: exponential growth -> log(y) is linear

# NOTE: the scoring function results are stored as a matrix for re-use below
segs <- dpseg(x=x, y=y, jumps=jumps, P=P, type=type, store.matrix=TRUE)

## using 'recursion_linreg_c' with type 'var'
## calculating recursion for 229 datapoints

print(segs)

## 
## Dynamic programming-based segmentation of 229 xy data points:
## 
##         x1       x2 start end  intercept     slope        r2         var
## 1 1.022778 1.362778     1   3  0.7194953 -4.383977 0.9909440 0.005076021
## 2 1.362778 1.703194     3   5 -9.3479269  3.016111 0.9626144 0.010235471
## 3 1.703194 2.213611     5   8 -6.2212774  1.153347 0.9766701 0.001532996
## ...
##         x1       x2 start end  intercept       slope        r2          var
## 7 15.68903 17.91181    87 100 -3.4426960  0.18899698 0.9940106 1.101408e-04
## 8 17.91181 21.32958   100 120 -0.3825874  0.01895679 0.8638382 6.374447e-05
## 9 21.32958 39.94000   120 229  0.4767462 -0.01989881 0.9919736 9.501867e-05
## 
## Found: 9  segments.
## Parameters:
##       type: var; minl: 3; maxl: 229; P: 1e-04; jumps: 0

plot(segs)
lines(predict(segs),lty=2, lwd=3, col="yellow")

Segmentation of a growth curve (E. coli in M9 minimal medium) by dpseg. Vertical lines indicate segment starts (dashed red) and ends (solid black). The predict method returns the piecewise linear model (dashed yellow).

Exploratory & Educational Movies

For both, educational purposes or detailed evaluation of novel scoring functions (Section @ref(generic)) or parameters, the movie function allows to visualize the performance of the dynamic programming as it progresses through the data:

# use arguments frames and res to control file size
movie(segs, format="gif", file.name="movie", path="pkg/vignettes/figures",
      frames=seq(1,length(x),3), delay=.3,res=50)

This will plot each step of the algorithm and visualize the development of the scoring function. For this function to work, the argument store.matrix=TRUE must be set in the call to dpseg.

Animation of the progress of the algorithm through the data (gray circles and right y-axis). The black vertical line is the current position j, and the black circles are the values of the scoring function score(i, j) for all i < j. The blue line was the i_max where the maximal value of the recursion S_j was found. The thin red line indicates $j-\lmin$. Colors indicate the final segmentation after backtracing.

Fine-tuning Segment Length: The Penalty Parameter P

The minimal and maximal segment length parameters $\lmin$ and $\lmax$ (arguments minl and maxl of dpseg) are the easiest way to restrict the algorithm (and memory usage) to certain segment lengths. However, this overrules potentially better segmentations with lower variance.

A more meaningful way to restrict segment lengths for given data is to explicitly allow higher variance of segments. This can be achieved with the break-point penalty parameter P (argument P). This parameter can be directly used to tune segmentation. It can be chosen in the order of magnitude of the tolerated variances. A higher P will allow higher variance of the individual segments and will yield longer segments. P < 0 will reward break-points and yield more shorter segments.

To choose an optimal P for a given application, the package offers two functions: estimateP and scanP. estimateP makes use of the excellent performance of base R’s smooth.spline implementation and reports the variance of the residuals; argument plot=TRUE allows to evaluate results:

p <- estimateP(x=x, y=y, plot=TRUE)

plot(dpseg(x=x, y=y, jumps=jumps, P=round(p,3)))

## using 'recursion_linreg_c' with type 'var'
## calculating recursion for 229 datapoints

The reported value is a good starting point. For our example data, values an order of magnitude lower than this estimate, P=estimateP(x=x,y=y)/10 achieved a satisfying segmentation, but this will strongly depend on the type of data and noise around expected segments.

Alternative estimators can be easily defined, eg. estimateP simply calls smooth.spline:

simple_estimateP <- function (x, y, ...) {
    var(smooth.spline(x, y, ...)$y - y)
}

The convenience function scanP calculates segmentations for a vector of P values and returns (plots) the resulting numbers of segments and the median of segment variances of residuals. Based on these values one can select a value of P that splits the data into a reasonable number of segments with acceptable variance of residuals:

scanP is more efficient than scanning over dpseg calls since it calculates the scoring function matrix only once and uses this matrix for scans over P (Section @ref(matrix)).

## NOTE: dpseg is slower for many segments!
sp <- scanP(x=x, y=y, P=seq(-.01,.1,length.out=50), plot=TRUE)

## running dpseg 50 times: ..................................................

A higher P will yield fewer and longer segments. P should be chosen close to the optimal value of the scoring function.

Scoring Functions: −Var(r) vs. R²

The performance of the different built-in scoring functions depends on the type of data and the goal of the user.

The default scoring function, the negative variance of residuals (type="var", $-\Var(r)$), is independent of the slope of segments (eq. @ref(eq:var)) and is recommended in most cases.

Only when looking explicitly for segmentation of zig-zag like data with only increases and decreases, the coefficient of determination (type="r2", R²) or its square root, Pearson’s correlation (type="cor"), may be better choices. In dpseg, these are not used in the sense of a goodness-of-fit measure but in their mathematical meaning as a measure of (linear) correlation: how well can x explain variation in y? Lines with low slope score worse than steeper lines (see slope b₁ in eq. @ref(eq:r2)). This behaviour is exemplified in the Sections @ref(example), @ref(lines), and @ref(generic). Figure @ref(fig:rcppdynprogdata) provides an example where the solution is stable over a broader range of penalty P when scoring by the Pearson correlation.

Data Example: P < 0 and Correlation-based Scoring.

Using the non-monotonic, sine-based example data from the RcppDynProg R package (Mount and Zumel 2019), we can first use scanP to test different penalties for different scoring functions. Here, dpseg yields the intended segmentation with the default scoring function (type="var") and a negative penalty, while the Pearson correlation (type="cor") worked well with default penalty P = 0 and appears to have a broader range of penalties that achieve the intented split into 7 segments:

## example from https://winvector.github.io/RcppDynProg/articles/SegmentationL.html
set.seed(2018)
d <- data.frame(x = 0.05*(1:(300))) # ordered in x
d$y_real <- sin((0.3*d$x)^2) 
d$y <- d$y_real + 0.25*rnorm(length(d$y_real))

par(mfrow=c(2,2))
scanP(x=d$x, y=d$y, P=seq(-0.005,0.005,length.out=50), verb=0) # note: no messages by verb=0
scanP(x=d$x, y=d$y, P=seq(-0.05,0.05,length.out=50), type="cor", verb=0)
plot(dpseg(x=d$x, y=d$y, P=-.0025, verb=0)) 
plot(dpseg(x=d$x, y=d$y, type="cor", verb=0))

Segmentation of the test data from the RcppDynProg package with the default scoring function, negative variance of residuals (left), and with Pearson correlation (right).

Vertical and Horizontal Lines with Zero Variance.

Experimental data may contain short horizontal or vertical data stretches, e.g., as artefacts of low measurement sensitivity. For such segments the calculation of the scoring function involves division by 0 and undefined values.

For a horizontal straight line with $\Var(y)=0$ (top panel in Fig. @ref(fig:lines)), i.e., a short segment of x values with equal y values, the R² and Pearson correlation measures are undefined (eq. @ref(eq:r2)). dpseg will set them to 0, their minimal value, consistent with their general behaviour to favor sloped and ignore horizontal segments. The variance of residuals can be calculated (eq. @ref(eq:var)) as $\Var(r)=0$ and the segments can be detected with the default scoring function (top left panel in Fig. @ref(fig:lines)).

For a vertical line with $\Var(x)=0$ (lower panel in Fig. @ref(fig:lines)), i.e., different y values at equal x, all measures are undefined (eq. @ref(eq:r2) and @ref(eq:var)). dpseg will set the scoring function to −∞, effectively ignoring them in segmentation. The data points will be part of adjacent or spanning segments. If such cases are present in the data between intended segments, it may be advisable to use jumps=TRUE to allow for split segments and discrete jumps between adjacent segments (see Section @ref(jumps)).

For both cases, dpseg will issue a summarized warning if it was relevant for the used scoring function, and providing the counts of occurrences of such cases during the recursion:

y <- c(1:5,rep(5,10),6:10) # several equal y values
x <- 1:length(y) 
par(mfrow=c(2,2))
plot(dpseg(x=x, y=y, type="var", verb=0)) # Var(y)=0, horizontal line
plot(dpseg(x=x, y=y, type="r2", verb=0))

## Warning in recursion(x = x, y = y, maxl = maxl, jumps = jumps, P = P, minl =
## minl, : Syy==0: equal y-values or numeric cancellation in segment, score was
## set to -1 (r2=cor=0). This occurred 45 times.

plot(dpseg(x=y, y=x, type="var", verb=0)) # Var(x)=0, vertical line, flip x<->y

## Warning in recursion(x = x, y = y, maxl = maxl, jumps = jumps, P = P, minl =
## minl, : Sxx==0: equal x-values or numeric cancellation in segment, score was
## set to -infinity. This occurred 45 times.

plot(dpseg(x=y, y=x, type="r2", verb=0))

## Warning in recursion(x = x, y = y, maxl = maxl, jumps = jumps, P = P, minl =
## minl, : Sxx==0: equal x-values or numeric cancellation in segment, score was
## set to -infinity. This occurred 45 times.

Effect of straight vertical (top) or horizontal (bottom) lines.

Discrete Jumps between Segments

Per default, dpseg assumes that adjacent segments are linked, that is: the last data point of a segment is also the first data point of the next adjacent segment. Setting the argument jumps=TRUE allows to obtain segments that do not share their start and end points (see $\jump$ in eq. @ref(eq:recur)). This can be useful if the data contains discrete jumps between adjacent segments (Fig. @ref(fig:jumps)).

x <- c(1:10) 
y <- c(1:4,7:12) # discrete jump between two lines
par(mfcol=c(1,2))
plot(dpseg(x=x, y=y, type="var", jumps=FALSE, verb=0))
plot(dpseg(x=x, y=y, type="var", jumps=TRUE, verb=0))

Handling discrete jumps in data: the default behaviour adds a third segment across the discrete jump between segments (left panel). Allowing for such jumps by setting argument jumps=TRUE correctly identifies the jump and the two adjacent segments (right panel).

Pre-Calculated Scoring Matrix

When the scoring function can be pre-calculated, a simple look-up in the triangular matrix (banded by $\lmin$ and $\lmax$) suffices for the recursion. This option is available in dpseg by passing a scoring matrix S_ij as argument y (instead of a numeric vector).

For example, we can use the scoring matrix generated above (Fig. @ref(fig:dpsegdemo), stored due to argument store.matrix=TRUE) and test different parameters P, $\jump$, $\lmin'>\lmin$ and $\lmax'<\lmax$ more efficiently (Section @ref(benchmark)). This approach is also used in the scanP function (Section @ref(selectp)).

## use the scoring matrix from a previous run for generic recursion,
## with store.matrix=TRUE, and test different parameters.
segm <- dpseg(y=segs$SCR, jumps=jumps, P=2*P, minl=5, maxl=50)

## setup use with pre-calculated scoring matrix
## calculating recursion for 229 datapoints

print(segm)

## 
## Dynamic programming-based segmentation of 229 xy data points:
## 
##   x1 x2 start end
## 1  1 37     1  37
## 2 37 46    37  46
## 3 46 87    46  87
## ...
##    x1  x2 start end
## 5 100 131   100 131
## 6 131 180   131 180
## 7 180 229   180 229
## 
## Found: 7  segments.
## Parameters:
##       type: matrix; minl: 5; maxl: 50; P: 2e-04; jumps: 0

Note, that in this case the algorithm does not use or know of the original x, y data and the result therefore contains only the segment break-point indices in the original data vectors. The predict and plot functions will not work. We can use the convenience function addLm to add the original x, y data and linear regression coefficients (via lm) for the calculated segments:

## add lm-based regression coefficients and original x/y data
segm <- addLm(segm, x=oddata$Time, y=log(oddata[,"A2"]))
plot(segm)

Compared to the first run (Fig. @ref(fig:dpsegdemo)) both, the higher $\lmin=5$ and higher P = 0.0002, contributed to get fewer segments, while the lower $\lmax=50$ split the last segment in two.

Custom Scoring Functions

Alternative scoring functions can be tested easily by either providing a scoring matrix S_ij as described above, or by providing the score(i,j) function directly. The latter can be achieved by defining a score function with the signature as in the following example, testing to use the coefficient of determination R² (r-squared) for segmentation:

Since our scoring function doesn’t provide linear regression parameters, we can use the argument add.lm=TRUE to add intercept and slope data via base R’s lm, required for the predict and plot methods.

score_rsq <- function(i, j, x, y,...) summary(lm(y[i:j]~x[i:j]))$r.squared
segn <- dpseg(x=x, y=y, jumps=jumps, P=.99, scoref=score_rsq, add.lm=TRUE)
plot(segn) 

Note, that as above a meaningful penalty term P should be close to the optimal value of the scoring function, here R² = 1. Further note, that the R² measure depends on the slope and misses a short almost horizontal segment around x = 20 that is detected by the −Var(r) default scoring function.

Least-Squares Method

Consider n measured data pairs {(x_i, y_i), i = 1, ..., n}, with the dependent variable y_i (eg. a measured value) and independent variable x_i (eg. the time of measurement), for which we suspect a linear relation with intercept β₀ and slope β₁, and with added random measurement errors, in regression analysis denoted the residuals r_i:

The goal is to find a straight line, denoted the regression line, that best describes this linear relation:

The gold standard approach to find the regression line is to minimize the sum of squares of the residuals: $\min\limits_{\beta_0,\beta_1} \sum r_i^2$. We treat the data as constants and look for parameters b₀ and b₁, estimators of the real parameters β₀ and β₁, that minimize the deviation of the data from this line.

$$ \begin{equation} f(\beta_0, \beta_1) = \sum_{i=1}^n r_i^2 = \sum_{i=1}^n (y_i-\beta_0-\beta_1 x_i)^2 (\#eq:beta) \end{equation} $$

Note that for clarity we shorten the sum symbol from now on as: $\sum = \sum_{i=1}^n$. It is important not to mix such estimators with the real parameters, since we do not know how accurate our estimation will be. The real parameters β₀ and β₁ reflect the actual physical process that underlies the data we measured.

To minimize the residuals, we get the partial derivatives and set them to 0.

$$ \begin{align} \frac{\partial f}{\partial \beta_0} \vert_{b_0} &= \sum 2\cdot (y_i-b_0-\beta_1 x_i)\cdot 1 &=0\\ &= 2 \sum y_i - 2 n b_0 - 2 \beta_1 \sum x_i &=0\\ b_0 &= \frac{1}{n}\sum y_i - \beta_1 \frac{1}{n} \sum x_i\\\hline \frac{\partial f}{\partial \beta_1} \vert_{b_1} &= \sum 2\cdot (y_i-\beta_0-b_1 x_i)\cdot x_i &=0\\ &= \sum x_i y_i - \beta_0 \sum x_i - b_1 \sum x_i^2 &=0 \end{align} $$

Combining both minimization criteria, with b₀ = β₀ and b₁ = β₁:

$$ \begin{align} \sum x_i y_i - \left(\frac{1}{n}\sum y_i - b_1 \frac{1}{n} \sum x_i\right) \sum x_i - b_1 \sum x_i^2 =0\\ \sum x_i y_i - \frac{1}{n}\sum y_i \sum x_i + b_1 \frac{1}{n} \left(\sum x_i\right)^2 - b_1 \sum x_i^2 =0\\ \sum x_i y_i - \frac{1}{n}\sum y_i \sum x_i = b_1 \left(\sum x_i^2 - \frac{1}{n} \left(\sum x_i\right)^2\right)\\ b_1 = \frac{\sum x_i y_i - \frac{1}{n}\sum y_i \sum x_i}{\sum x_i^2 - \frac{1}{n} \left(\sum x_i\right)^2} \end{align} $$

Introducing the arithmetic mean, $\bar{x} = \frac{1}{n} \sum_{i=1}^n x_i$, and using the centered form of the variance, we obtain the well known equation for linear regression:

$$ \begin{align} b_0 &= \bar y - b_1 \bar x\\ b_1 &= \frac{\sum(x_i-\bar x )(y_i-\bar y)}{\sum(x_i-\bar x)^2}\;. (\#eq:emin) \end{align} $$

This can be understood intuitively, especially the centered form. It contains the sum of the squared differences of the x_i data points from their arithmetic mean x̄, as well as a similar expression with the products of x and y data. These “sums of squares (products)” (SQ) are a measure of the spread of the data around their mean, known as variance (covariance).

Required Differentiation Rules:

• (f(x) + g(x))′ = f′(x) + g′(x),
• (f(x)²)′ = 2 ⋅ f(x) ⋅ f′(x).

Required Summation Rules:

• $\sum_{i=1}^n (a_i+b_i) = \sum_{i=1}^n a_i + \sum_{i=1}^n b_i$,
• $\sum_{i=1}^n(c a_i) = c \sum_{i=1}^n a_i$,
• $\sum_{i=1}^n c = n c$,
• also good to know, no rule for $\sum_{i=1}^n (a_i b_i)$.

Incremental Calculation

The last step to obtain equation @ref(eq:emin) is known in German as Verschiebungssatz von Steiner (Steiner translation theorem), in French as théorème de König-Huygens. In English this transformation has no special name and the two forms are just known as the centered (left) and expanded (right) forms of the variance:

We can apply this to all sum of squares required to calculate various statistical measures of the data:

$$ \begin{align} S_{xx} \def & \sum(x_i - \bar x)^2 &=& \sum x_i^2 - \frac{1}{n} (\sum x_i)^2\\ S_{yy} \def & \sum(y_i - \bar y)^2 &=& \sum y_i^2 - \frac{1}{n} (\sum y_i)^2\\ S_{xy} \def & \sum(y_i - \bar y)(x_i - \bar x) &=& \sum y_i x_i - \frac{1}{n} \sum y_i\sum x_i (\#eq:sos) \end{align} $$

Their expanded forms consist of additions of simple sums of the data, and will thus allow to calculate the sum of squares incrementally by adding up the sum expressions ∑x_i², ∑x_i, ∑y_i², ∑y_i, and ∑x_iy_i.

An important caveat is that the expanded forms can become 0 for large x̄², known as “catastrophic cancellation”.

Useful Definitions

For completeness and context, we provide some standard definitions in statistics literature.

Incremental Calculation of Scoring Functions

Let’s keep these definitions in mind but return to linear regression. The aim was to obtain expressions that we can calculate incrementally. We have found a regression line that minimizes the sum of squares of the residuals, the distance of actual measured values y from the regression line ŷ_i.