[R] [RsR] function in nls argument -- robust estimation

Stromberg, Arnold astro11 at email.uky.edu
Sun May 11 02:11:50 CEST 2008 

Martin, Kate, Fernando, et al.

Be careful bootstrapping robust estimators. The trouble is that when resampling is done with replacement, the outliers can be selected too many times which would ruin the standard error estimates. Martin is right that bootstrapping would be fine if there are not too many outliers. Otherwise, Jackknifing will likely work better, especially if you use a delete more than one version. For a zero weight for outliers M estimator, a high breakdown starting value like least trimmed squares would be a good idea.

arny

Arnold J. Stromberg
Professor and Chair
Department of Statistics
University of Kentucky
817 Patterson Office Tower
Lexington, KY 40506-0027
Phone: 859-257-6115
Fax: 859-323-1973
________________________________________
From: r-sig-robust-bounces at r-project.org [r-sig-robust-bounces at r-project.org] On Behalf Of Martin Maechler [maechler at stat.math.ethz.ch]
Sent: Saturday, May 10, 2008 5:07 AM
To: Katharine Mullen; elnano
Cc: r-help at r-project.org; R-SIG-robust at r-project.org
Subject: Re: [RsR] [R] function in nls argument -- robust estimation

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.

_______________________________________________
R-SIG-Robust at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-robust

More information about the R-help mailing list