[R] Confindence interval for Levenberg-Marquardt fit

[Joerg van den Hoff]

On Wed, Feb 21, 2007 at 09:41:29AM -0600, Douglas Bates wrote:
> On 2/21/07, joerg van den hoff <j.van_den_hoff at fzd.de> wrote:
> >On Wed, Feb 21, 2007 at 11:09:52AM +0000, Prof Brian Ripley wrote:
> >> Well, the algorithm used does not affect the confidence interval 
> >(provided
> >> it works correctly), but what is nls.ml (presumably in some package you
> >> have not mentioned) and why would I want to use an old-fashioned
> >> algorithm?
> 
> >is'nt this a bit strong? in what respect do you consider 
> >levenberg-marquardt
> >(going back to the 19-forties, I think) as old-fashioned (especially
> >in comparsion to the `nls' standard gauss-newton approach (both gentlemen
> >seem to have done their major work a bit longer ago :-))?
> >AFAICS levenberg-marquardt is generally appreciated for it's rapid
> >convergence achieved by a smooth transition from an
> >inverse-hessian approach to steepest descent. my own experience
> >with non-linear least squares minimization using this algorithm
> >are positive as well, but
> >I have not tried out the levenberg-marquardt
> >implementation in package `minpack.lm' (which originates from netlib.org)
> >and don't know if it's good. but in any case there sure are implementations
> >around (e.g. in the CERN MINUIT library) which have proven to be
> >of high quality.
> 
> The nls function provides implementations of the Gauss-Newton
> algorithm and the Golub-Pereyra algorithm and the NL2SOL algorithm.
> Check the possible values of the optional argument "algorithm" to
> nls().
I'm using `nls' rather frequently and I'm aware of this. maybe the formulation
"standard gauss-newton approach" was misleading (please attribute this to the
fact that I'm not a native speaker: I actually meant the default algorithm used
and was not implying that `nls' can only use some odd "standard"
procedure)
> 
> If you will allow me to wax philosophical for a moment, I submit that
> there are three stages to any iterative optimization:
>  - formulating initial values of the parameters
>  - determining what step to take from iteration i to iteration i+1
>  - determining when to declare convergence
> 
> For the most part optimization research focuses on the second step.
> The first and third steps are also important.  One of the advantages
> of the implementation of the Gauss-Newton and Golub-Pereyra algorithms
> in nls is that they use a convergence criterion (based only on the
> current parameter values) not a termination criterion (based on the
> changes in the parameter values or the changes in the objective
> function or both).  See chapter 2 of Bates and Watts, "Nonlinear
> Regression Analysis and Its Applications" (Wiley, 1988) for a
> discussion of the particular convergence criterion used.
I did this when starting to use `nls' seriously. being a physicist,
not a mathematician, my attitude is nevertheless a bit more pragmatic
(or simple-minded, if you like): I fully appreciate that these
things need a thorough and rigorous mathematical foundation but any
convergence criterion suits me just fine as long as (1) convergence
is actually achieved (i.e. the algorithm reports some results...)
and (2) the final least squares sum is sufficiently near the
achievable minimum value for that model (sufficiently near meaning
that the parameter estimation errors are much larger than any further
changes in the parameters when wandering around any longer near the
current point in parameter space -- possibly finding a slightly
"better mininmum") and (3) the results are obtained in sufficiently
short time.  It's not obivous to me (I did'nt encounter an
example case up to now), for instance, that a good termination
criterion is either more unreliable (yielding wrong answers way off
the true minimum) or less efficient (e.g. driving
the iteration count up) than the convergence criterion used in `nls'
(against which I do not have any objections -- apart from the minor
handicap that it can't converge on 'ideal' data which by and then
comes in the way as the help list shows)

an impression which I cannot prove by hard data (so it may be wrong) is,
that the `nls' gauss-newton approach tends to be more sensitive to
suboptimal start values (sending it towards outer space in the worst
case...) than levenberg-marquardt driven approaches and definitely
has on average a higher iteration count.

> 
> Another point that is often lost in a discussion of various
> optimization algorithms is that one does not compare algorithms - one
> compares implementations of algorithms which, at the very least, must
> involve all of the above steps.
sure
> 
> For many people in the nonlinear least squares field the NL2SOL
> algorithm by Dennis, Gay and Welch was considered the "gold standard"
> and David Gay's implementation is now available in nls.  This allows
> those who are interested to compare implementations of these
> algorithms.  Also, with R being an Open Source system anyone can add
> their own implementation or modify an existing implementation to try
> things out.
for one, I already stressed that I appreciate `nls'. it
is, for me, one of the most useful things in R. your second point
is true in principle, of course, but question is how many people
will dare and try to fiddle with the nls code (which seems rather
hermetic to the non-initiate). I'd love to see levenberg-marquardt among
the available alorithms but I'm far from sure whether I could manage
to integrate it adequately into `nls' myself (the author of the
minpack package seems to have had the same problem choosing to
implement it separately).

and while I'm at it: another thing for the wish-list (I think you
put it somewhere at position 10^(some_small_integer_larger_than_1)
on your TODO list a couple years ago :-)) would be to get access
to the parameter search path in parameter space (essentially the
trace output) as a part of the returned `nls'-object (maybe filling
that component only if trace = TRUE?). As a quick work around I
only managed to copy getPars() to a global variable by modifying
`nlsModel' accordingly but found no easy clean solution.

> 
> 
> 
> 
> >`nls' sure is a _very_ valuable function, but not necessarily the
> >"last word" with regards to the chosen algorithm(s).
> >
> >
> >>
> >> You could start nls at the solution you got from nls.ml and use confint()
> >> on that.
> >maybe one should look at profile.nls and confint.nls and see what 
> >information
> >of the usual `nls' object is actually used for the confidence intervall
> >computation and mimick this for the `nls.lm' output? at a (admittedly)
> >quick glance it seems that only parameters, std.errs. and the 
> >fitted/residual
> >values are needed which should all be provided by nls.lm as well.
> >maybe one could even try to map the nls.lm results into a structure of 
> >class
> >`nls' (although this would not be a clean solution, probably) in order
> >to use `confint.nls'?
> 
> Not really.  The profile.nls function requires, not surprisingly, the
> ability to profile the objective function.  (See the discussion of the
> profile sum of squares and the profile-t transformations in chapter 6
> of Bates and Watts).  Such a capability would need to be added to a
> nonlinear least squares implementation before methods for those
> generics could be developed.
I see.
> 
> >>
> >> On Wed, 21 Feb 2007, Michael Dondrup wrote:
> >>
> >> > Dear all,
> >> > I would like to use the Levenberg-Marquardt algorithm for non-linear
> >> > least-squares regression using function nls.lm. Can anybody help  me to
> >> >   find a a way to compute confidence intervals on the  fitted
> >> > parameters as it is possible for nls (using confint.nls, which does not
> >> > work for nls.lm)?
> >> >
> >> > Thank you for your help
> >> > Michael
> >>
> >>
> >
> >______________________________________________
> >R-help at stat.math.ethz.ch 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.
> >