Thank you all for your help. I think I now understand the issue.
I tried to write a likelihood function for my binomial model.
Please excuse my ignorance if I am not doing this right; I do not have any
statistical background.
#Example data
x <- seq(0, 1000)
y <- ifelse(x < 300, 0, ifelse(x>700, 0, 1))
#Function assuming binomial errors
fn <- function(par) {
y.pred = ifelse(x<par[1], 0, ifelse(x>par[2], 0, 1))
-sum(dbinom(1-abs(y-y.pred), prob = mean(y.pred == y), size = 1, log = T))
}
optim(par=c(300,700), fn)
Would the "Nelder-Mead" optimization method be appropriate?
On Tue, Oct 16, 2012 at 11:28 AM, Ravi Varadhan <ravi.varad...@jhu.edu>wrote:
> 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
> ====================================
>
[[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.
