Piecewise linear regression in R (segmented.lm)

Question

I appreciate any help to make segmented.lm (or any other function) find the obvious breakpoints in this example:

data = list(x=c(50,60,70,80,90) , y= c(703.786,705.857,708.153,711.056,709.257))
plot(data, type='b')
require(segmented)
model.lm = segmented(lm(y~x,data = data),seg.Z = ~x, psi = NA)

It returns with the following error:

Error in solve.default(crossprod(x1), crossprod(x1, y1)) : system is computationally singular: reciprocal condition number = 1.51417e-20

If I change K:

model.lm = segmented(lm(y~x,data = data),seg.Z = ~x, psi = NA, control = seg.control(K=1))

I get another error:

Error in segmented.lm(lm(y ~ x, data = data), seg.Z = ~x, psi = NA, control = seg.control(K = 1)) : only 1 datum in an interval: breakpoint(s) at the boundary or too close each other

I hope that's just a sample data set because that's really far too few observations to try to estimate something like that. The problem when you don't specify K is that the default value is 10. Since you have fewer than observations, this default doesn't make sense and you can't differentiate between models. The problem with K=1 is that it's trying the break at just after 80 which leaves only one point in that second group making it impossible to calculate a slope. But you really need a much larger input data set. — MrFlick
Like @MrFlick said, you have too few observations here. Try: data = data.frame(x=c(50,60,70,80,90,100,110) , y= c(703.786,705.857,708.153,711.056,709.257, 707.4, 705.6)). Then this will work: model.lm = segmented(lm(y~x,data = data),seg.Z = ~x, psi = NA, control = seg.control(K=1)) — Jota
Thank you MrFlick and Frank. I will use a larger dataset eventually. But for this sample set, I could manage to get a model with: library("SiZer") model.pwl = piecewise.linear(x = data$x, y = data$y, middle = 1, CI = FALSE, bootstrap.samples = 1000, sig.level = 0.05) plot(model.pwl) — user1507435

Chris Ruehlemann Chris Ruehlemann · Accepted Answer · 2020-06-26T09:21:38

A nice objective method to determine the break point is described in Crawley (2007: 427).

First, define a vector breaks for a range of potential break points:

breaks <- data$x[data$x >= 70 & data$x <= 90]

Then run a for loop for piecewise regressions for all potential break points and yank out the minimal residual standard error (mse) for each model from the summary output:

mse <- numeric(length(breaks))
for(i in 1:length(breaks)){
  piecewise <- lm(data$y ~ data$y*(data$x < breaks[i]) + data$y*(data$x >= breaks[i]))
  mse[i] <- summary(piecewise)[6]
}
mse <- as.numeric(mse)

Finally, identify the break point with the least mse:

breaks[which(mse==min(mse))]

Hope this helps.

Piecewise linear regression in R (segmented.lm)

1 Answers