[RsR] function in nls argument -- robust estimation

[Katharine Mullen]

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)