Selecting the number of breakpoints/changepoints in segmented/stepmented regression

This function selects (and estimates) the number of breakpoints/changepoints of the segmented/stepmented relationship according to the BIC/AIC criterion or sequential hypothesis testing (via the pseudo-Score test).

Usage

selgmented(olm, seg.Z, Kmax=10, type=c("bic", "aic", "score", "davies"), 
  alpha = 0.05, control = seg.control(), refit = FALSE, stop.if = 6, 
  return.fit = TRUE, bonferroni = FALSE, msg = TRUE, plot.ic = FALSE, th = NULL, 
  G = 1, check.dslope = TRUE)

stelpmented(olm, seg.Z, Kmax=10, type=c("bic", "aic", "score", "davies"), 
  alpha = 0.05, control = seg.control(), refit = FALSE, stop.if = 6, 
  return.fit = TRUE, bonferroni = FALSE, msg = TRUE, plot.ic = FALSE, th = NULL, 
  G = 1, check.dslope = TRUE, a=1, Cn=1)

Arguments

olm: A starting lm or glm object, or a simple numerical vector meaning the response variable.
seg.Z: A one-side formula for the segmented variable. Only one term can be included, and it can be omitted if olm is a (g)lm fit including just one numeric covariate. Also seg.Z might be omitted, and will be taken as 1,2..., if olm is a single numeric vector.
Kmax: The maximum number of breakpoints being tested. If type='bic' or type='aic', any integer value can be specified; otherwise at most Kmax=2 breakpoints can be tested via the Score or Davies statistics.
type: Which criterion should be used? Options "score" and "davies" allow to carry out sequential hypothesis testing with no more than 2 breakpoints (Kmax=2). Alternatively, the number of breakpoints can be selected via the BIC (or AIC) with virtually no upper bound for Kmax.
alpha: The fixed type I error probability when sequential hypothesis testing is carried out (i.e. type='score' or 'davies'). It is also used when type='bic' (or type='aic') and check.dslope=TRUE to remove the breakpoints based on the slope diffence t-value.
control: See seg.control.
refit: If TRUE, the final selected model is re-fitted using arguments in control, typically with bootstrap restarting. Set refit=FALSE to speed up computation (and possibly accepting near-optimal estimates). It is always TRUE if type='score' or type='davies'.
stop.if: An integer. If, when trying models with an increasing (when G=1) or decreasing (when G>1) number of breakpoints, stop.if consecutive fits exhibit higher AIC/BIC values, the search is interrupted. Set a large number, larger then Kmax say, if you want to assess the fits for all breakpoints 0, 1, 2, ..., Kmax. Ignored if type='score' or type='davies'.
return.fit: If TRUE, the fitted model (with the number of breakpoints selected according to type) is returned.
bonferroni: If TRUE, the Bonferroni correction is employed, i.e. alpha/Kmax (rather than alpha) is always taken as threshold value to reject or not. If FALSE, alpha is used in the second level of hypothesis testing. It is also effective when type="bic" (or 'aic') and check.dslope=TRUE, see Details.
msg: If FALSE the final fit is returned silently with the selected number of breakpoints, otherwise the messages about the selection procedure (i.e. the BIC values), and possible warnings are also printed.
plot.ic: If TRUE the information criterion values with respect to the number of breakpoints are plotted. Ignored if type='score' or type='davies' or G>1. Note that if check.dslope=TRUE, the final number of breakpoints could differ from the one selected by the BIC/AIC leading to an inconsistent plot of the information criterion, see Note below.
th: When a large number of breakpoints is being tested, it could happen that 2 estimated breakpoints are too close each other, and only one can be retained. Thus if the difference between two breakpoints is less or equal to th, one (the first) breakpoint is removed. Of course, th depends on the x scale: Integers, like 5 or 10, are appropriate if the covariate is the observation index. Default (NULL) means th=diff(range(x))/100. Set th=0 if you are willing to consider even breakpoints very clode each other. Ignored if type='score' or type='davies'.
G: Number of sub-intervals to consider to search for the breakpoints when type='bic' or 'aic'. See Details.
check.dslope: Logical. Effective only if type='bic' or 'aic'. After the optimal number of breakpoints has been selected (via AIC/BIC), should the \(t\) values of the slope differences be checked? If TRUE, the breakpoints corresponding to slope differences with a 'low' \(t\) values will be removed. Note the model is re-fitted at each removal and a new check is performed. Simulation evidence suggests that such strategy leads to better results. See Details.
a: An additional tuning parameter for the BIC. a=1 provides the classical BIC (see Details).
Cn: An additional tuning parameter for the BIC. Cn=1 provides the classical BIC (see Details).

Details

The function uses properly the functions segmented, pscore.test or davies.test to select the 'optimal' number of breakpoints 0,1,...,Kmax. If type='bic' or 'aic', the procedure stops if the last stop.if fits have increasing values of the information criterion. Moreover, a breakpoint is removed if too close to other, actually if the difference between two consecutive estimates is less then th. Finally, if check.dslope=TRUE, breakpoints whose corresponding slope difference estimate is `small' (i.e. \(p\)-value larger then alpha or alpha/Kmax) are also removed.

When \(G>1\) the dataset is split into \(G\) groups, and the search is carried out separately within each group. This approach is fruitful when there are many breakpoints not evenly spaced in the covariate range and/or concentrated in some sub-intervals. G=3 or 4 is recommended based on simulation evidence.

Note Kmax is always tacitely reduced in order to have at least 1 residual df in the model with Kmax changepoints. Namely, if \(n=20\), the maximal segmented model has 2*(Kmax + 1) parameters, and therefore the largest Kmax allowed is 8.

When type='score' or 'davies', the function also returns the 'overall p-value' coming from combing the single p-values using the Fisher method. The pooled p-value, however, does not affect the final result which depends on the single p-values only.

Arguments a and Cn in stelpmented() modify the classical penalty term in the BIC, namely \(p*(log(n)^a)*Cn\) where \(p\) is the number of model parameters. I don't know which are the most appropriate values for them.

Note

If check.dslope=TRUE, there is no guarantee that the final model has the lowest AIC/BIC. Namely the model with the best A/BIC could have `non-significant' slope differences which will be removed (with the corresponding breakpoints) by the final model. Hence the possible plot (obtained via plot.ic=TRUE) could be misleading. See Example 1 below.

Value

The returned object depends on argument return.fit. If FALSE, the returned object is a list with some information on the compared models (i.e. the BIC values), otherwise a classical 'segmented' object (see segmented for details) with the component selection.psi including the A/BIC values and
- if refit=TRUE, psi.no.refit that represents the breakpoint values before the last fit (with boot restarting)
- if G>1, cutvalues including the cutoffs values used to split the data.

References

Muggeo V (2020) Selecting number of breakpoints in segmented regression: implementation in the R package segmented https://www.researchgate.net/publication/343737604

Author

Vito M. R. Muggeo

Examples


## breakpoint selection (segmented)

set.seed(12)
xx<-1:100
zz<-runif(100)
yy<-2+1.5*pmax(xx-35,0)-1.5*pmax(xx-70,0)+15*pmax(zz-.5,0)+rnorm(100,0,2)
dati<-data.frame(x=xx,y=yy,z=zz)
out.lm<-lm(y~x,data=dati)

os <-selgmented(out.lm, type="score", Kmax=2) #selection (Kmax=2) via the Score test 
#> Error in eval(mf): object 'dati' not found

os <-selgmented(out.lm, type="bic", Kmax=5) #BIC-based selection
#> No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 
#> 
#> BIC to detect no. of breakpoints:
#>        0        1        2        3        4        5 
#> 716.3031 696.9431 545.1816 547.3707 555.7056 563.0675 
#> 
#> No. of selected breakpoints: 2  


## changepoint (stepmented)
yy<-2+1.5*(xx>30)-1.5*(xx>65)+1.5*(xx>80)+rnorm(100,0,.5)


os <-stelpmented(yy, Kmax=5)  
#> No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 
#> 
#> BIC to detect no. of breakpoints:
#>          0          1          2          3          4          5 
#>  -20.42470  -61.54908  -65.58944 -113.82651 -108.05171 -105.89439 
#> 
#> No. of selected breakpoints: 3  

#selection accounting for an additional linear covariate: 
out.lm<-lm(yy~zz)
os <-stelpmented(out.lm, ~xx, Kmax=5) 
#> Error in eval(parse(text = nomeX)): object 'xx' not found




if (FALSE) { # \dontrun{
########################################
#Example 1: selecting a large number of breakpoints in segmented regression

b <- c(-1,rep(c(1.5,-1.5),l=15))
psi <- seq(.1,.9,l=15)
n <- 2000
x <- 1:n/n
X <- cbind(x, outer(x,psi,function(x,y)pmax(x-y,0)))
mu <- drop(tcrossprod(X,t(b)))
set.seed(113)
y<- mu + rnorm(n)*.022 
par(mfrow=c(1,2))

#select number of breakpoints via the BIC (and plot it)
o<-selgmented(y, Kmax=20, type='bic', plot.ic=TRUE, check.dslope = FALSE) 
plot(o, res=TRUE, col=2, lwd=3)
points(o)
# select via the BIC + check on the slope differences (default)
o1 <-selgmented(y, Kmax=20, type='bic', plot.ic=TRUE) #check.dslope = TRUE by default
#note the plot of BIC is misleading.. But the number of psi is correct 
plot(o1, add=TRUE, col=3)
points(o1, col=3, pch=3)

##################################################
#Example 2: a large number of breakpoints not evenly spaced.  

b <- c(-1,rep(c(2,-2),l=10))
psi <- seq(.5,.9,l=10)
n <- 2000
x <- 1:n/n
X <- cbind(x, outer(x,psi,function(x,y)pmax(x-y,0)))
mu <- drop(tcrossprod(X,t(b)))
y<- mu + rnorm(n)*.02 

#run selgmented with G>1. G=3 or 4 recommended. 
#note G=1 does not return the right number of breaks  
o1 <-selgmented(y, type="bic", Kmax=20, G=4)
} # }