Errors in segmented package: breakpoints confusion

Question

Using the segmented package to create a piecewise linear regression I am seeing an error when I try to set my own breakpoints; it seems only when I try to set more than two.

(EDIT) Here is the code I am using:

# data
bullard <- structure(list(Rt = c(0, 4.0054, 25.1858, 27.9998, 35.7259, 39.0769, 
45.1805, 45.6717, 48.3419, 51.5661, 64.1578, 66.828, 111.1613, 
114.2518, 121.8681, 146.0591, 148.8134, 164.6219, 176.522, 177.9578, 
180.8773, 187.1846, 210.5131, 211.483, 230.2598, 262.3549, 266.2318, 
303.3181, 329.4067, 335.0262, 337.8323, 343.1142, 352.2322, 367.8386, 
380.09, 388.5412, 390.4162, 395.6409), Tem = c(15.248, 15.4523, 
16.0761, 16.2013, 16.5914, 16.8777, 17.3545, 17.3877, 17.5307, 
17.7079, 18.4177, 18.575, 19.8261, 19.9731, 20.4074, 21.2622, 
21.4117, 22.1776, 23.4835, 23.6738, 23.9973, 24.4976, 25.7585, 
26.0231, 28.5495, 30.8602, 31.3067, 37.3183, 39.2858, 39.4731, 
39.6756, 39.9271, 40.6634, 42.3641, 43.9158, 44.1891, 44.3563, 
44.5837)), .Names = c("Rt", "Tem"), class = "data.frame", row.names = c(NA, 
-38L))

library(segmented)

# create a linear model
out.lm <- lm(Tem ~ Rt, data=bullard)

o<-segmented(out.lm, seg.Z=~Rt, psi=list(Rt=c(200,300)), control=seg.control(display=FALSE))

Using the psi option, I have tried the following:

psi = list(x = c(150, 300)) -- OK
psi = list(x = c(100, 200)) -- OK
psi = list(x = c(200, 300)) -- OK
psi = list(x = c(100, 300)) -- OK
psi = list(x = c(120, 150, 300)) -- error 1 below
psi = list(x = c(120, 300)) -- OK
psi = list(x = c(120, 150)) -- OK
psi = list(x = c(150, 300)) -- OK
psi = list(x = c(100, 200, 300)) -- error 2 below

(1) Error in segmented.lm(out.lm, seg.Z = ~Rt, psi = list(Rt = c(120, 150, : only 1 datum in an interval: breakpoint(s) at the boundary or too close

(2) Error in diag(Cov[id, id]) : subscript out of bounds

I have already listed my data at this question, but as a guide the limits on the x data are about 0--400.

A second question that pertains to this one is: how do I actually fix the breakpoints using this segmented package?

Surmise

The data may be in the other question, but can you turn that into a reproducible example, showing exactly the function call you have been trying? – Sin

If you are using seg.lm.fit the docs for PSI say "appropriate matrix including the starting values of the breakpoints to be estimated" – Sin

Hi sean, thanks for response. i added my code, though I'm not sure how to easily convert my data into script (in-line so to speak), so i left in the read.delim call. – Surmise

Maybe you could add the output of dput(head(bullard)) to give a sense of what the data looks like. – Sin

@Sin ah that's lovely. thanks, didn't know how to do it. The dataset is small enough to have it all in there. – Surmise

Answer 1

The issue here seems to be poor error trapping in the segmented package. Having a look at the code for segmented.lm allows a bit of debugging. For example, in the case of psi = list(x = c(100, 200, 300)), an augmented linear model is fitted as shown below:

lm(formula = Tem ~ Rt + U1.Rt + U2.Rt + U3.Rt + psi1.Rt + psi2.Rt + 
    psi3.Rt, data = mf)

Call:
lm(formula = Tem ~ Rt + U1.Rt + U2.Rt + U3.Rt + psi1.Rt + psi2.Rt + 
    psi3.Rt, data = mf)

Coefficients:
(Intercept)           Rt        U1.Rt        U2.Rt        U3.Rt      psi1.Rt        
   15.34303      0.04149      0.04591    742.74186   -742.74499      1.02252       
   psi2.Rt      psi3.Rt  
        NA           NA  

As you can see, the fit has NA values which then result in a degenerate variance-covariance matrix (called Cov in the code). The function doesn't check for this and tries to pull out diagonal entries from Cov and fails with the error message shown. At least the first error, although perhaps not overly helpful, is caught by the function itself and suggests that the break-points are too close.

In the absence of better error trapping in the function, I think that all you can do is adopt a trial and error approach (and avoid break points which are too close). For example, psi = list(x = c(50, 200, 300)) seems to work ok.

Sin

This error really is frustrating when you know your breakpoints precisely. – Unloose

The five year old error seems still to be existing. Funny thing: when I repeat my segmented() call, the error somewhen dissapears! – Albritton

Answer 2

If you use while and tryCatch you can make the command repeat itself until it decides there is no error in the model @jaySf. I'm guessing this is down to the randomiser settings in the function, which can be seen in seg.control.

lm.model <- lm(xdat ~ ydat, data = x)
if.false <- F
while(if.false == F){
  tryCatch({
    s <- segmented(lm.model, seg.Z =~ydata, psi = NA)
    if.false <- T
  }, error = function(e){
  }, finally = {})
}

Bettor

There is a risk this runs forever, any way to try only a hundred times and then give up? – Snuffer

You may add a loop index, such as: lm.model <- lm(xdat ~ ydat, data = x) if.false <- F idx <- 1L while(if.false == F & idx <= 100L){ tryCatch({ s <- segmented(lm.model, seg.Z =~ydata, psi = NA) if.false <- T }, error = function(e){ }, finally = {}) idx <- idx + 1L } – Bimolecular