# Selecting the number of breakpoints in segmented regression This function selects (and estimates) the number of breakpoints of the segmented relationship according to the BIC/AIC criterion or sequential hypothesis testing. ```r selgmented(olm, seg.Z, Kmax=2, type = c("score", "bic", "davies", "aic"), alpha = 0.05, control = seg.control(), refit = FALSE, stop.if = 5, return.fit = TRUE, bonferroni = FALSE, msg = TRUE, plot.ic = FALSE, th = NULL, G = 1, check.dslope = TRUE) ``` ## 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 it might be omitted, and will be taken as 1,2..., if `olm` includes a single numeric variable. - `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. ## 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. ## 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. ## Returns 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(s) Vito M. R. Muggeo ## See Also `segmented`, `pscore.test`, `davies.test` ## Examples ```r 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) #selection (Kmax=2) via the Score test (default) os <-selgmented(out.lm, type="bic", Kmax=3) #BIC-based selection ## Not run: ######################################## #Example 1: selecting a large number of breakpoints 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) ## End(Not run) ```

selgmented function