Segmented regression

Question

I am working with data consisting of observations over time (abundance and harvest of different game species over several years). I am interested in finding if and when the direction is changed (i.e. an increase turns into a decrease). For this I find that segmented regression (package segmented) should be useful.

The first point is that this function creates as many breaks as is specified in the psi-part of the function. Since I do not know how many breaks there should be I adopted a strategy where I test with one, two, three and four breaks and then looked at the AIC values to pick the best function (lowest AIC). This works well for most datasets, but I get some peculiar-looking results also (looks like over-fitting). So my first question is if this is a proper way to use segmented regression.

Furthermore, the functions seems sensitive to the breakpoints suggested in the psi-part of the function. If I change these suggestions, the break-points decided by the function may change sligthly for the same data. This seems a bit subjective to me. How to handle this? I have tried to use the psi=NA but only get an error message.

And finally, by accident I ran the same regression twice, using exactly the same data and the same specifications. And the results changed slightly, both the estimates and the AIC value. Why is this and how to handle it?

Answer 1

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

The steps involved are:

First, define a vector breaks for a range of potential break points (V_expl, and further below, V_resp stand for explanatory variable and response variable respectively):

breaks <- V_expl[V_expl >= ... & V_expl <= ...]

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:

mse <- numeric(length(breaks))
for(i in 1:length(breaks)){
  piecewise <- lm(V_resp ~ V_expl*(V_expl < breaks[i]) + V_expl*(V_expl >=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.