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.
