2

Possible Duplicate:
Change point analysis using R's nls()

I want to do a nonlinear regression with nls() but also include a specific type of segmented or piecewise regression. The Formula I want to implement is:

S ~ b0 + (A > T) * b1 * (A - T)

T should be the threshold value or breakpoint as identified by the nonlinear-segmented regression. I know that I can use "algorithm = plinear" but that does not work at all.

The data I have is:

A   S
0.000809371 1
0.003642171 3
0.009712455 4
0.010521827 2
0.004046856 4
0.015378054 5
0.000404686 0
0.000404686 0
0.000404686 0
0.000809371 0
0.000809371 3
0.037635765 3
0.008903084 2
0.016187426 5
0.043301364 1
0.000404686 1
0.002428114 1
0.003642171 1
0.013759312 4
0.051395077 9
0.394568501 9
0.005665599 1
0.013354626 1
0.028732681 3
0.026304567 2
0.004451542 1
0.050585705 2
0.00647497  1
0.010926512 0
0.013354626 1
1.695632841 4
0.013354626 2
Community
  • 1
Jens
  • 1,384
  • 1
  • 14
  • 31
  • @whuber : On second look, this isn't really a duplicate actually. The solution for this problem is not evident for the other. And the solution of the other question is far too complicated for this one, as shown in my edited answer. – Joris Meys Jul 04 '11 at 12:09
  • @Joris The functional form `S ~ b0 + (A > T) * b1 * (A - T)` is equivalent to `S ~ b0 + b1*max(0,A-T)`. Changing notation to $x$ = A, $\delta$ = T, $\beta_0$ = b0, $\beta_2$ = b2, and setting $\beta_1=0$ gives a special case of the problem solved in the duplicate post. If you have improvements to suggest for the solution there, then it would be best also to post your reply there so we can keep this common thread together. – whuber Jul 04 '11 at 16:31
  • @whuber : off course, my mistake. I'll add an answer there. – Joris Meys Jul 05 '11 at 08:10

1 Answers1

5

EDIT:

As requested by OP, here a minimal example of how to include breakpoints in a nls. Actually, after reviewing the question and answers that is considered a duplicate, I believe they go way around to come at a solution. It can be as simple as :

Data <- data.frame(
  A = rnorm(100),
  S = sample(0:9,100,TRUE)
)

X <- nls(S~ B0 + B1*as.numeric(A > Tval)*(A-Tval),
        start=list(B0=1,B1=1,Tval=1),data=Data)

So you make sure that the ID variable is numeric (i.e. 0 if FALSE and 1 if TRUE) and proceed as usual.

Consider constructing a dummy variable :

Data <- data.frame(
  A = rnorm(100),
  S = sample(0:9,100,TRUE)
)
Tval <- 1
Data$Id <- ifelse( Data$A > Tval , 1, 0 ) # the dummy variable

nls(S~ B0 + B1*Id*(A-Tval),data=Data)

Joris Meys
  • 5,475
  • 2
  • 32
  • 43
  • The O.P. states that the method needs to identify `Tval`, too. – whuber Jun 29 '11 at 12:44
  • @whuber : I see, overread it. it's anyway solved in the duplicate question. I answered it thinking it was still on SO. – Joris Meys Jun 30 '11 at 08:15
  • Hi, I did not understand how the other duplicate post relates to this one. I do not have any prior information on b0 and b1. So how should I use them in the formula? Joris, could you try to put that into the example you made? would really help me out! – Jens Jul 04 '11 at 11:57
  • @Jens you're right. I added an example of how to do that. – Joris Meys Jul 04 '11 at 12:06
  • I only get "singular gradient" or "singular gradient matrix at first parameter estimation (translated from german", no matter which data I use. is there a problem with the data? – Jens Jul 05 '11 at 10:09
  • @Jens : It's a classic. I've been looking at it, and the problem is that the equation as it stands can cause non-identifiable cases. This will lead the more sophisticated methods like nls() to fail. There are other options, but for that I suggest you read the other question. The solution there is definitely what you're looking for, as your case is just a special case of the more general problem presented there. I will add a few other options for piecewise later on, but I'm too busy right now to write it up. Cheers – Joris Meys Jul 05 '11 at 10:25
  • @Joris: Many thanks, I really appreciate your help and of course I understand that you cannot spent all day answering. I will try some more around with and the segmented package while hoping for an update. many thanks again – Jens Jul 05 '11 at 10:30