Message-ID: <Pine.GSO.4.56.0805101532110.12828@laurel.few.vu.nl>
Date: 2008-05-10T14:02:08Z
From: Katharine Mullen
Subject: [RsR] function in nls argument -- robust estimation
There has been recent discussion on R-help
(https://stat.ethz.ch/pipermail/r-help/2008-May/161611.html) regarding a
nonlinear least squares problem in which residuals above a certain
quantile are set to zero; the intent "is to avoid the strong influence of
points which push my modeled vs. measured values away from the 1:1 line".
Martin Maechler pointed out the connection of this approach to robust
nonlinear regression
(https://stat.ethz.ch/pipermail/r-help/2008-May/161786.html). Below is a
script for the problem under plain nls, nlrob, and nls.lm with the
residuals above a certain quantile zeroed out; maybe someone on this list
or the OP Fernando wants to comment on how the methods compare?
ss <-
read.table(url("ftp://ftp.bgc-jena.mpg.de/pub/outgoing/fmoyano/data_nls.lm_moyano.txt"),
skip=1)
ST <- ss[!is.na(ss[,3]),1]
SM <- ss[!is.na(ss[,3]),2]
SR <- ss[!is.na(ss[,3]),3]
## fit with nls
nlsRes <- nls(SR ~ (a*SM^2+b*SM+c) * exp(E*((1/(283.15-227.13)) -
(1/(ST+273.15-227.13)))),
start=list(a=-0.003, b=0.13, c=0.50, E=400),
control = list(printEval=TRUE))
## fit with nlrob
library(robustbase)
nlrobRes <- nlrob(SR ~ (a*SM^2+b*SM+c) * exp(E*((1/(283.15-227.13)) -
(1/(ST+273.15-227.13)))),
data=data.frame(SR),
start=list(a=-0.003, b=0.13, c=0.50, E=400),
trace=TRUE)
## note: adding control= list(printEval=TRUE) w/ and w/out trace=TRUE
##prints strange results
## fit with nls.lm and residual zeroing based on quantile
optim.f <- function(p, ST, SM, SR, q) {
res <- SR - m.f(p)
abserr <- abs(res)
qnum <- quantile(abserr, probs=q, na.rm=TRUE)
res[abserr > qnum] <- 0
res
}
m.f <- function(p) (p$a*SM^I(2)+p$b*SM+p$c)*exp(p$E*((1/(283.15-227.13))-
(1/(ST+273.15-227.13))))
library(minpack.lm)
q <- 0.90
p <- list(a=-0.003, b=0.13, c=0.50, E=400)
nls.lmRes <- nls.lm(par = p, fn = optim.f, ST = ST, SM = SM, SR = SR,
q = q, control = nls.lm.control(nprint=1))
plot(SR, type="p")
lines(fitted(nlsRes), col=3)
## note: the below line doesn't work; why is nlrobRes$fitted.values
## empty?
##lines(fitted(nlrobRes), col=2)
lines(SR-resid(nlrobRes), col=2)
lines(m.f(coef(nls.lmRes)), col=4)
legend(1000,5,legend=c("nls","nlrob","nls.lm w/0'ing"),lwd=2,col=c(3,2,4))
## observed vs predicted and 1:1 line
par(mfrow=c(3,1))
plot(SR~fitted(nlsRes), main="nls")
lines(seq(0,6,1),seq(0,6,1),col=4)
robustResult<-(SR-resid(nlrobRes))
plot(SR~robustResult,main="nlrob")
lines(seq(0,6,1),seq(0,6,1),col=4)
plot(SR~m.f(coef(nls.lmRes)),main="nls.lm w/0'ing")
lines(seq(0,6,1),seq(0,6,1),col=4)