Veronique, Can you write down the likelihood function in R and send it to me (this is very easy to do, but I don't want to do your work)? Also send me the code for simulating the data. I will show you how to fit such models using optimization tools.
Best, Ravi ________________________________________ From: Vito Muggeo [vito.mug...@unipa.it] Sent: Tuesday, October 16, 2012 9:55 AM To: Bert Gunter Cc: Véronique Boucher Lalonde; r-help@r-project.org; Ravi Varadhan Subject: Re: [R] fit a "threshold" function with nls 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.lalo...@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@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 ==================================== ______________________________________________ R-help@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.