[R] fit a "threshold" function with nls

[Vito Muggeo]

Véronique,
in addition to Bert's comments, I would like to bring to your attention 
that there are several packages that perform 
threshold/breakpoint/changepoint estimation in R, including

cumSeg, segmented, strucchange, and bcp for a Bayesian approach

Moreover some packages, such as cghFLasso, perfom smoothing with a L1 
penalty that can yield a step-function fit.

I hope this helps you,

vito


Il 15/10/2012 22.21, Bert Gunter ha scritto:
> Véronique:
>
> I've cc'ed this to a true expert (Ravi Varadhan) who is one of those
> who can give you a definitive response, but I **believe** the problem
> is that threshhold type function fits have  objective functions whose
> derivatives are discontinuous,and hence gradient -based methods can
> run into the sorts of problems that you see.
>
> **If** this is so, then you might do better using an explicit
> non-gradient optimizer = rss minimizer via one of the optim() suite of
> functions (maybe even the default Nelder-Mead); but this is definitely
> where the counsel of an expert would be valuable. Also check out the
> CRAN Optimization Task View for advice on optimization options.
>
> Cheers,
> Bert
>
> On Mon, Oct 15, 2012 at 9:43 AM, Véronique Boucher Lalonde
> <veronique.boucher.lalonde at gmail.com> wrote:
>> I am trying to model a dependent variable as a threshold function of
>> my independent variable.
>> What I mean is that I want to fit different intercepts to y following 2
>> breakpoints, but with fixed slopes.
>> I am trying to do this with using ifelse statements in the nls function.
>> Perhaps, this is not an appropriate approach.
>>
>> I have created a very simple example to illustrate what I am trying to do.
>>
>> #Creating my independent variable
>> x <- seq(0, 1000)
>> #Creating my dependent variable, where all values below threshold #1 are 0,
>> all values above threshold #2 are 0 and all values within the two
>> thresholds are 1
>> y <- ifelse(x < 300, 0, ifelse(x>700, 0, 1))
>> #Adding noise to my dependent variable
>> y <- y + rnorm(length(y), 0, 0.01)
>> #I am having trouble clearly explaining the model I want to fit but perhaps
>> the plot is self explanatory
>> plot(x, y)
>> #Now I am trying to adjust a nonlinear model to fit the two breakpoints,
>> with known slopes between the breakpoints (here, respectively 0, 1 and 0)
>> threshold.model <- nls(y~ifelse(x<b1, 0, ifelse(x>b2, 0, 1)),
>> start=list(b1=300, b2=700), trace=T)
>>
>> I have played around with this function quite a bit and systematically get
>> an error saying: singular gradient matrix at initial parameter estimates
>> I admit that I don't fully understand what a singular gradient matrix
>> implies.
>> But it seems to me that both my model and starting values are sensible
>> given the data, and that the break points are in fact estimable.
>> Could someone point to what I am doing wrong?
>>
>> More generally, I am interested in knowing
>> (1) whether I can use such ifelse statements in the function nls
>> (1) whether I can even use nls for this model
>> (3) whether I can model this with a function that would allow me to assume
>> that the errors are binomial, rather than normally distributed
>> (ultimately I want to use such a model on binomial data)
>>
>> I am using R version 2.15.1 on 64-bit Windows 7
>>
>> Any guidance would be greatly appreciated.
>> Veronique
>>
>>          [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
>
>

-- 
====================================
Vito M.R. Muggeo
Dip.to Sc Statist e Matem `Vianelli'
Università di Palermo
viale delle Scienze, edificio 13
90128 Palermo - ITALY
tel: 091 23895240
fax: 091 485726
http://dssm.unipa.it/vmuggeo