[R] function in nls argument -- robust estimation

Katharine Mullen kate at few.vu.nl
Sat May 10 15:29:00 CEST 2008


Dear Martin,

Thanks for the ideas regarding the relation of what Fernando is doing with
robust regression.  Indeed, it's an important point that he can't consider
the standard error estimates on his parameters correct.

I know from discussion off-list that he's happy with the results he has
now; nevertheless the robust regression route may be an interesting
alternative.  I'm posting a scipt to R-SIG-robust now that compares the 3
ways (nls, nlrob and nls.lm w/residuals above a certain quantile set to
zero).

best,
Kate

On Sat, 10 May 2008, Martin Maechler wrote:

> Hi Kate and Fernando,
>
> I'm late into this thread,
> but from reading it I get the impression that Fernando really
> wants to do *robust* (as opposed to least-squares) non-linear
> model fitting.  His proposal to set residuals to zero when they
> are outside a given bound is a very special case of an
> M-estimator, namely (if I'm not mistaken) the so-called "Huber
> skipped-mean", an M-estimator with psi-function
>    psi <- function(x, k) ifelse(abs(x) <= k, x, 0)
> It is known that this can be far from optimal, and either using
> Huber-psi or "a redescender" such as Tukey's biweight can be
> considerably better.
> Also note that the standard inference (std.errors, P-values, ...)
> that you'd get from summary(nlsfit) or anova(nls1, nl2) is
> *invalid* here, since you are effectively using *random* weighting.
>
> The nlrob() function in package 'robustbase'
> implements M-estimation of nonlinear models directly.
> Unfortunately, how to do correct inference in this situation
> is a hard problem, probably even an open research question in
> parts. I would expect that "the" bootstrap should work if you only
> have a few outliers.
>
> I don't have time at the moment to look at the example data and
> the model, and show you how to use it for nlrob();
> if you find a way to you it for nls() , then the same should
> work for nlrob().
>
> I'm CCing this to the specialists for "Robust Stats with R"
> mailing list, R-SIG-robust.
>
> Best regards,
> Martin Maechler
> ETH Zurich
>
> >>>>> "KateM" == Katharine Mullen <kate at few.vu.nl>
> >>>>>     on Fri, 9 May 2008 15:50:08 +0200 (CEST) writes:
>
>     KateM> You can take minpack.lm_1.1-0 (source code and MS Windows build,
>     KateM> respectively) from here:
>
>     KateM> http://www.nat.vu.nl/~kate/minpack.lm_1.1-0.tar.gz
>     KateM> http://www.nat.vu.nl/~kate/minpack.lm_1.1-0.zip
>
>     KateM> The bug that occurs when nprint = 0 is fixed.  Also fixed is another
>     KateM> problem suggested your example: when the argument par is a list, calling
>     KateM> summary on the output of nls.lm was not working.
>
>     KateM> I'll submit the new version to CRAN soon.
>
>     KateM> This disscusion has been fruitful - thanks for it.
>
>     KateM> On Fri, 9 May 2008, Katharine Mullen wrote:
>
>     >> You indeed found a bug.  I can reproduce it (which I should have tried to
>     >> do on other examples in the first place!).  Thanks for finding it.
>     >>
>     >> It will be fixed in version 1.1-0 which I will submit to CRAN soon.
>     >>
>     >> On Fri, 9 May 2008, elnano wrote:
>     >>
>     >> >
>     >> > Find the data (data_nls.lm_moyano.txt) here:
>     >> > ftp://ftp.bgc-jena.mpg.de/pub/outgoing/fmoyano
>     >> >
>     >> >
>     >> >
>     >> > Katharine Mullen wrote:
>     >> > >
>     >> > > Thanks for the details - it sounds like a bug.  You can either send me the
>     >> > > data in an email off-list or make it available on-line somewhere, so that
>     >> > > I and other people can download it.
>     >> > >
>     >> > >
>     >> > > ______________________________________________
>     >> > > 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.
>     >> > >
>     >> > >
>     >> >
>     >> > --
>     >> > View this message in context: http://www.nabble.com/function-in-nls-argument-tp17108100p17146812.html
>     >> > Sent from the R help mailing list archive at Nabble.com.
>     >> >
>     >> > ______________________________________________
>     >> > 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.
>     >> >
>     >>
>     >> ______________________________________________
>     >> 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.
>     >>
>
>     KateM> ______________________________________________
>     KateM> R-help at r-project.org mailing list
>     KateM> https://stat.ethz.ch/mailman/listinfo/r-help
>     KateM> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>     KateM> and provide commented, minimal, self-contained, reproducible code.
>



More information about the R-help mailing list