Hi, I'm trying to use Robust non-linear regression to fit dose response curves. Maybe I didnt look good enough, but I dind't find robust methods for NON linear regression implemented in R. A method that looked good to me but is unfortunately not (yet) implemented in R is described in http://www.graphpad.com/articles/RobustNonlinearRegression_files/frame.htm <http://www.graphpad.com/articles/RobustNonlinearRegression_files/frame.htm> in short: instead of using the premise that the residuals are gaussian they propose a Lorentzian distribution, in stead of minimizing the squared residus SUM (Y-Yhat)^2, the objective function is now SUM log(1+(Y-Yhat)^2/ RobustSD) where RobustSD is the 68th percentile of the absolute value of the residues my question is: is there a smart and elegant way to change to objective function from squared Distance to log(1+D^2/Rsd^2) ? or alternatively to write this as a weighted non-linear regression where the weights are recalculated during the iterations in nlme it is possible to specify weights, possibly that is the way to do it, but I didn't manage to get it working the weights should then be something like: SUM (log(1+(resid(.)/quantile(all_residuals,0.68))^2)) / SUM (resid(.)) the test data I use : x<-seq(-5,-2,length=50) x<-rep(x,4) y<-SSfpl(x,0,100,-3.5,1) y<-y+rnorm(length(y),sd=5) y[sample(1:length(y),floor(length(y)/50))]<-200 # add 2% outliers at 200 thanks a lot Hans Vermeiren [[alternative HTML version deleted]]
you might consider nlrq() in the quantreg package, which does median regression for nonlinear response functions.... url: www.econ.uiuc.edu/~roger Roger Koenker email rkoenker at uiuc.edu Department of Economics vox: 217-333-4558 University of Illinois fax: 217-244-6678 Champaign, IL 61820 On Nov 13, 2005, at 3:47 PM, Vermeiren, Hans [VRCBE] wrote:> Hi, > > I'm trying to use Robust non-linear regression to fit dose response > curves. > Maybe I didnt look good enough, but I dind't find robust methods > for NON > linear regression implemented in R. A method that looked good to me > but is > unfortunately not (yet) implemented in R is described in > http://www.graphpad.com/articles/RobustNonlinearRegression_files/ > frame.htm > <http://www.graphpad.com/articles/RobustNonlinearRegression_files/ > frame.htm> > > > in short: instead of using the premise that the residuals are > gaussian they > propose a Lorentzian distribution, > in stead of minimizing the squared residus SUM (Y-Yhat)^2, the > objective > function is now > SUM log(1+(Y-Yhat)^2/ RobustSD) > > where RobustSD is the 68th percentile of the absolute value of the > residues > > my question is: is there a smart and elegant way to change to > objective > function from squared Distance to log(1+D^2/Rsd^2) ? > > or alternatively to write this as a weighted non-linear regression > where the > weights are recalculated during the iterations > in nlme it is possible to specify weights, possibly that is the way > to do > it, but I didn't manage to get it working > the weights should then be something like: > > SUM (log(1+(resid(.)/quantile(all_residuals,0.68))^2)) / SUM (resid > (.)) > > the test data I use : > x<-seq(-5,-2,length=50) > x<-rep(x,4) > y<-SSfpl(x,0,100,-3.5,1) > y<-y+rnorm(length(y),sd=5) > y[sample(1:length(y),floor(length(y)/50))]<-200 # add 2% outliers > at 200 > > thanks a lot > > Hans Vermeiren > > > [[alternative HTML version deleted]] > > ______________________________________________ > 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 >
> -----Original Message----- > From: r-help-bounces at stat.math.ethz.ch [SMTP:r-help-bounces at stat.math.ethz.ch] On Behalf Of Vermeiren, Hans [VRCBE] > Sent: Sunday, November 13, 2005 7:48 PM > To: 'r-help at stat.math.ethz.ch' > Subject: [R] Robust Non-linear Regression > > Hi, > > I'm trying to use Robust non-linear regression to fit dose response curves. > Maybe I didnt look good enough, but I dind't find robust methods for NON > linear regression implemented in R. A method that looked good to me but is > unfortunately not (yet) implemented in R is described in > http://www.graphpad.com/articles/RobustNonlinearRegression_files/frame.htm > <http://www.graphpad.com/articles/RobustNonlinearRegression_files/frame.htm> > > > in short: instead of using the premise that the residuals are gaussian they > propose a Lorentzian distribution, > in stead of minimizing the squared residus SUM (Y-Yhat)^2, the objective > function is now > SUM log(1+(Y-Yhat)^2/ RobustSD) > > where RobustSD is the 68th percentile of the absolute value of the residues > > my question is: is there a smart and elegant way to change to objective > function from squared Distance to log(1+D^2/Rsd^2) ? >----------- I do not know about in-built robustness options in R but I have found that Dave Fournier's robust likelihood for nonlinear regression in ADMB does a pretty good job in detecting and counter-acting the influence of outliers (in my applications this has been used to counter-act the effect of reading errors in determination of the age of fish based on rings in bones). It relies on a likelihood function based on a mixture of a normal and another distribution with fatter tails. You can find the documentation in the ADMB manual at the ADMB website: http://otter-rsch.com/admodel.htm Ruben
Package 'sfsmisc' has had a function 'rnls()' for a while
which does robust non-linear regression via M-estimation.
[The name of the function is probably *really* a misnomer,
since the 'ls' part stands for "least squares"!]
Two weeks ago, there's been a small workshop
"Robustness and R" in Treviso (It),
http://www.dst.unive.it/rsr/
where we've talked about available and missing robustness
functionality ``in R''. One consequence of the workshop is the
new mailing list "R-SIG-Robust" {to which I CC this message}
and another planned and hopefully even more consequential
consequence will be collaboration on producing more widely
available robustness functionality for R. Do subscribe to the
list if you are interested.
>>>>> "Vermeiren" == Vermeiren, Hans [VRCBE]
<Vermeiren>
>>>>> on Sun, 13 Nov 2005 22:47:31 +0100 writes:
Vermeiren> Hi, I'm trying to use Robust non-linear
Vermeiren> regression to fit dose response curves. Maybe I
Vermeiren> didnt look good enough, but I dind't find robust
Vermeiren> methods for NON linear regression implemented in
Vermeiren> R. A method that looked good to me but is
Vermeiren> unfortunately not (yet) implemented in R is
Vermeiren> described in
Vermeiren>
http://www.graphpad.com/articles/RobustNonlinearRegression_files/frame.htm
Vermeiren> in short: instead of using the premise that the
Vermeiren> residuals are gaussian they propose a Lorentzian
Vermeiren> distribution, in stead of minimizing the squared
Vermeiren> residus SUM (Y-Yhat)^2, the objective function is
Vermeiren> now SUM log(1+(Y-Yhat)^2/ RobustSD)
Vermeiren> where RobustSD is the 68th percentile of the
Vermeiren> absolute value of the residues
Vermeiren> my question is: is there a smart and elegant way
Vermeiren> to change to objective function from squared
Vermeiren> Distance to log(1+D^2/Rsd^2) ?
no; not easily.
Vermeiren> or alternatively to write this as a weighted
Vermeiren> non-linear regression where the weights are
Vermeiren> recalculated during the iterations in nlme it is
Vermeiren> possible to specify weights, possibly that is the
Vermeiren> way to do it, but I didn't manage to get it
Vermeiren> working the weights should then be something
Vermeiren> like:
Vermeiren> SUM (log(1+(resid(.)/quantile(all_residuals,0.68))^2))
Vermeiren> / SUM (resid(.))
rnls() mentioned does use robust weights and IRLS (iteratively
reweighted LS) making use of nls() and rlm(),
similarly to your suggestion.
Vermeiren> the test data I use :
Vermeiren> x<-seq(-5,-2,length=50)
Vermeiren> x<-rep(x,4)
Vermeiren> y<-SSfpl(x,0,100,-3.5,1)
Vermeiren> y<-y+rnorm(length(y),sd=5)
Vermeiren> y[sample(1:length(y),floor(length(y)/50))]<-200 # add 2%
outliers at 200
Since you have only outliers in 'y' and none in 'x',
you could use the 'nlrq' (nonlinear regression quantiles)
package that Roger Koenker mentioned.
To really robustify such self starting models as the 4-parameter
logistic 'SSfpl' above, you would also need to provide a robust
initial estimator;
maybe that could be done pretty easily 'rlm()' instead of 'lm()'
and
using 'rnls()' instead of 'nls()' also for the
"initial" part in
something like
SSfpl.rob <-
selfStart(~ A + (B - A)/(1 + exp((xmid - input)/scal)),
initial = function( ...) { ...... },
parameters=
c("A","B","xmid","scal"))
{look at 'SSfpl() for the initial estimator}.
However, BTW, currently the "plinear" version fails for our robust
nonlinear procedure 'rnls()'.
Martin Maechler, ETH Zurich
thank you all for the valuable suggestions rnls() is indeed what I was looking for I've to apologize to Roger Koenker for not mentioning that I did try quantile regression (saw his answer in a previous post with a similar question, yes i did my homework) however, least medians regression gave not always satisfying results, I now understand that this is in fact due to variability in the concentrations (x-axis) (thanks to Martin Maechlers remark), my example dataset was in that sense a bit unfortunate regards Hans Vermeiren -----Original Message----- From: Martin Maechler [mailto:maechler at stat.math.ethz.ch] Sent: Monday, November 14, 2005 12:41 PM To: Vermeiren, Hans [VRCBE] Cc: 'r-help at stat.math.ethz.ch'; R-SIG-robust at stat.math.ethz.ch Subject: Re: [R] Robust Non-linear Regression Package 'sfsmisc' has had a function 'rnls()' for a while which does robust non-linear regression via M-estimation. Since you have only outliers in 'y' and none in 'x', you could use the 'nlrq' (nonlinear regression quantiles) package that Roger Koenker mentioned.