[RsR] robust fitting using Lorentzian merit function
Martin Maechler
m@ech|er @end|ng |rom @t@t@m@th@ethz@ch
Thu Oct 16 16:13:39 CEST 2014
(Joachim has agreed to make this a public reply)
>>>>> Joachim Audenaert <Joachim.Audenaert using pcsierteelt.be>
>>>>> on Wed, 15 Oct 2014 13:43:13 +0200 writes:
> Dear Dr. Maechler, Thank you very much for your quick
> response.
> I am quite new to R, so please tell me if I gave data or
> code in a wrong way, all code is in bold.
(I'm reading E-mail in a "plain text" interface -- yes, old-fashioned --
so I don't see any "bold", but luckily your e-mail program
sends both, "HTML like" and "plain text like").
> this is my example dataset using dput, it is a table with
> FR and N0 values.
[I append a more complete R script later which includes the data]
> I want to make a regression model of this with nls, with
> the following formula RogersII <-
> nls(FR~(N0*attackR)/(1+(Th*N0*attackR)),start=list(attackR=1.5,Th=0.5),control=list(maxiter=10000))
> Because this is a biological experiment, there are often
> outliers. To get rid of them according to a method
> described in:
> http://www.biomedcentral.com/1471-2105/7/123; I should
> first fit a robust model, which minimizes the residuals
> with the assumption that they follow a Lorentizan
> distribution in stead of a normal distribution. This
> because with a Lorentzian distribution outliers have less
> impact on the model.
Well, it is a good idea to use robust procedures which do work
also when the errors are not normal.
This indeed has been a main focus of robust statistics since the 1980s
(even earlier: first cornerstone publication Huber (1964) .. Annals of Stats).
However, robust statistics methodology has generalized the
problem and the solution by using M-estimators (and later MM-,
CM-, S-, and more), all of these being more reasonable, than
assuming a "Lorentzian" -- which we typically call Cauchy
distribution or t_1 (t-distrib. with 1 degree of freedom).
> Currently, to make a robust fit I use:
> RogersII_R <-
> nlrob(FR~(N0*attackR)/(1+(Th*N0*attackR)),data=N0,FR,start=list(attackR=1.5,Th=0.5))
(the 'data=..' cannot be correct, but see also my R script below).
Using nlrob() from package robustbase is a very good idea, indeed...
> But this fit does not meet the criterion of residuals
> following a lorentzian distribution. Is it possible to
> extend the possibilities of nlrob, so the residuals can
> minimlized by assuming they follow a lorentzian
> distribution in stead of a normal distribution?
Using nlrob() you will use a method which is more sophisticated
than assuming a t_1 (aka Cauchy aka lorentzian) distribution for
the errors. Also it is *not* a two stage method (robust fit ->
outliers -> least squares) as in the biomed paper above, but
rather a "all-in-one" method which provides asymptotically
correct standard errors and inference where as
a two-stage method as the above leads to clearly *invalid* inference, i.e. all
P-values from the final least squares fit will be wrong "by
definition"...
In short: do use nlrob() and keep its result.
But then, as you asked, and I want to be nice (:-):
It *is* possible to use nlrob() with the correct 'psi = ...'
argument to directly compute the M-estimator which corresponds
to the MLE (maximum likelihood) for "t_1" aka Cauchy aka
lorentzian error distribution.
I don't have time to give you the details at the moment, but
hope to get there before the weekend.
Here's the promised R script with your original example and some
of the experiments I did.
Two more remarks:
1) You will notice that your data example is exactly a situation
where an "MM" estimate (or "CM" or ..) is clearly considerably
better than "M" because you have so called "leverage points".
2) Actually the parametrization you use in *non*linear
regression is important.
Here I believe it would be better to use iTh <- 1/Th
as parameter. I have not yet demonstrated it though ...
I need to send this e-mail out, as I'm busy for a while with
other matters.
Best regards,
Martin Maechler
------------------------------------------------------------------------------
## From: Joachim Audenaert <Joachim.Audenaert using pcsierteelt.be>
## To: Martin Maechler <maechler using stat.math.ethz.ch>
## [........]
## This is my example dataset using dput, it is a table with FR and N0 values.
d.NF <- data.frame(
N0 = c(10, 4, 1, 3, 5, 4, 7, 2, 7, 19,
12, 7, 15, 17, 10, 16, 17, 11, 10, 18, 15, 26, 20,
25, 30, 31, 28, 25, 36, 24, 46, 28, 40, 58, 55, 60,
47, 31, 60, 71, 18, 17, 69, 64, 74, 53, 135, 113,
161, 102, 156, 78),
FR = c(10, 4, 1, 3, 5, 1, 6, 2,
6, 18, 9, 7, 14, 15, 3, 16, 8, 11, 10, 15, 15, 25,
20, 25, 24, 31, 28, 21, 31, 24, 34, 24, 29, 32, 36,
39, 38, 28, 43, 41, 13, 17, 37, 42, 43, 33, 36, 24,
36, 40, 43, 21))
str(d.NF)
## I want to make a regression model of this with nls, with the following
## formula
## {Joachim called his models "RogersII"}
fm <- nls(FR ~ (N0*attackR)/(1+(Th*N0*attackR)),
data = d.NF, start = list(attackR = 1.5, Th = 0.5))
summary(fm)
## Because this is a biological experiment, there are often outliers. To get
## rid of them according to a method described in:
## http://www.biomedcentral.com/1471-2105/7/123; I should first fit a robust
## model, which minimizes the residuals with the assumption that they follow
## a Lorentizan distribution in stead of a normal distribution. This because
## with a Lorentzian distribution outliers have less impact on the model.
## Currently, to make a robust fit I use:
require(robustbase)
fm.R <-
nlrob(FR ~ (N0*attackR)/(1+(Th*N0*attackR)),
data = d.NF,
start = list(attackR = 1.5,Th = 0.5))
summary(fm.R)
## MM: Compute predictions on grid :
d.new <- with(d.NF, data.frame(N0 = seq(min(N0), max(N0), length=512+1)))
d.new[, "ls.fit"] <- predict(fm, d.new)
d.new[,"rob.fit"] <- predict(fm.R, d.new)
## Visualize
plot(FR ~ N0, data=d.NF, ylim = c(0,90)); abline(h=0,v=0, lty=3, col="gray")
cols <- c(LS = adjustcolor("brown",.6),
M = adjustcolor("red", .6))
lines( ls.fit ~ N0, d.new, col=cols["LS"], lty=2, lwd=2)
lines(rob.fit ~ N0, d.new, col=cols["M"], lty=1, lwd=2)
legend("bottomright", c("L.squares", "robust ('M', def.)"), bty="n",
lwd=2, lty=2:1, col = cols)
## MM: also try the other (new) methods for nlrob():
(meths <- eval(formals(nlrob)$method))
if(FALSE)# fails: they need (lower,upper) differently for different meth
allRm <- lapply(meths, function(.M.) update(fm.R, method = .M.))
allRm <- as.list(setNames(nm=meths))
allRm[["M"]] <- fm.R
## and the other methods:
allRm[-1] <- lapply(meths[-1],
function(.M.)
update(fm.R, method = .M.,
lower = c(attackR = 1, Th = 0),
upper = c(attackR = 3, Th = 0.2)))
t(sapply(allRm, coef))
## attackR Th
## M 1.482408 0.017703029
## MM 1.218437 0.012441474
## tau 1.166282 0.011655022
## CM 1.128051 0.011370435
## mtl 1.077946 0.006229957
## and the plots --- these fits are even more robust:
fit.rob <- sapply(allRm, function(m) predict(m, d.new))
matlines(d.new[,"N0"], fit.rob[,c("MM","tau","CM")], col="blue2")
## the asymptotes for N0 --> oo (Inf) :
abline(h = 1/sapply(allRm[c("M","MM","tau","CM")], coef)["Th",],
col=c(cols["M"], rep("blue2",3)), lty=1:4)
## More experiments: Try 'M' with different psi function:
## Redescending psi: starting value matters
summary(fm.R2 <- update(fm.R, psi = .Mwgt.psi1("lqq"))
## -> complete non-sense (local minimum)
## Try better "MM" with good start (and show iterations):
fm.R3 <- update(fm.R, psi = .Mwgt.psi1("lqq"), method="MM",
lower = c(attackR = 1, Th = 0),
upper = c(attackR = 3, Th = 0.2),
trace=TRUE)
## much more reasonable
More information about the R-SIG-Robust
mailing list