[R] function in nls argument -- robust estimation

Fernando Moyano nanomoyano at yahoo.com
Tue May 20 12:39:25 CEST 2008


Hi Kate and others,
thanks for the info.
Btw, you sent the different
methods to analyze the data: nls, nls.lm and nlrob. Comparing the
results visually nlrob performed better then nls, but nls.lm (using the
0.9 quantile of residuals) was still better than nlrob. My data may
have a rather large amount of contamination, so that an M-estimator
with a higher breakdown point should be used (least trimmed squares?).
I haven't found this in R and wouldn't know how to implement it. But I
can live with my results. Then remains the question of obtaining the
parameter st. errors. Jackknife was suggested. Is there an R function I
could use for that?

cheers,
Fernando



Katharine Mullen wrote:
> 
> 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.
>>
> 
> ______________________________________________
> 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-tp17108100p17337520.html
Sent from the R help mailing list archive at Nabble.com.



More information about the R-help mailing list