I'm trying to implement a "change point" analysis, or a multiphase regression using nls() in R.

Here's some fake data I've made. The formula I want to use to fit the data is:

$y = \beta_0 + \beta_1x + \beta_2\max(0,x-\delta)$

What this is supposed to do is fit the data up to a certain point with a certain intercept and slope ($\beta_0$ and $\beta_1$), then, after a certain x value ($\delta$), augment the slope by $\beta_2$. That's what the whole max thing is about. Before the $\delta$ point, it'll equal 0, and $\beta_2$ will be zeroed out.

So, here's my function to do this:

changePoint <- function(x, b0, slope1, slope2, delta){ 
   b0 + (x*slope1) + (max(0, x-delta) * slope2)

And I try to fit the model this way

nls(y ~ changePoint(x, b0, slope1, slope2, delta), 
    data = data, 
    start = c(b0 = 50, slope1 = 0, slope2 = 2, delta = 48))

I chose those starting parameters, because I know those are the starting parameters, because I made the data up.

However, I get this error:

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

Have I just made unfortunate data? I tried fitting this on real data first, and was getting the same error, and I just figured that my initial starting parameters weren't good enough.


(At first I thought it could be a problem resulting from the fact that max is not vectorized, but that's not true. It does make it a pain to work with changePoint, wherefore the following modification:

changePoint <- function(x, b0, slope1, slope2, delta) { 
   b0 + (x*slope1) + (sapply(x-delta, function (t) max(0, t)) * slope2)

This R-help mailing list post describes one way in which this error may result: the rhs of the formula is overparameterized, such that changing two parameters in tandem gives the same fit to the data. I can't see how that is true of your model, but maybe it is.

In any case, you can write your own objective function and minimize it. The following function gives the squared error for data points (x,y) and a certain value of the parameters (the weird argument structure of the function is to account for how optim works):

sqerror <- function (par, x, y) {
  sum((y - changePoint(x, par[1], par[2], par[3], par[4]))^2)

Then we say:

optim(par = c(50, 0, 2, 48), fn = sqerror, x = x, y = data)

And see:

[1] 54.53436800 -0.09283594  2.07356459 48.00000006

Note that for my fake data (x <- 40:60; data <- changePoint(x, 50, 0, 2, 48) + rnorm(21, 0, 0.5)) there are lots of local maxima depending on the initial parameter values you give. I suppose if you wanted to take this seriously you'd call the optimizer many times with random initial parameters and examine the distribution of results.

  • $\begingroup$ This post by Bill Venables explains well the issues involved in this kind of analysis. $\endgroup$
    – Aaron
    Feb 23 '11 at 7:07
  • 6
    $\begingroup$ Instead of that (cumbersome) sapply call in your first code snippet, you can always just use pmax. $\endgroup$
    – cardinal
    Apr 14 '11 at 19:00

Just wanted to add that you can do this with many other packages. If you want to get an estimate of uncertainty around the change point (something nls cannot do), try the mcp package.

# Simulate the data
df = data.frame(x = 1:100)
df$y = c(rnorm(20, 50, 5), rnorm(80, 50 + 1.5*(df$x[21:100] - 20), 5))

# Fit the model
model = list(
  y ~ 1,  # Intercept
  ~ 0 + x  # Joined slope
fit = mcp(model, df)

Let's plot it with a prediction interval (green line). The blue density is the posterior distribution for the change point location:

# Plot it
plot(fit, q_predict = T)

You can inspect individual parameters in more detail using plot_pars(fit) and summary(fit).

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.