[RsR] [R] nlrob problem
Pascal A Niklaus <pascal.niklaus at ieu.uzh.ch>
on Mon, 19 Dec 2011 20:46:57 +0100 writes:
> Dear all, I am not sure if this mail is for R-help or
> should be sent to R-devel instead, and therefore post to
> both.
> While using nlrob from package 'robustbase', I ran into
> the following problem:
> For psi-functions that can become zero
> (e.g. psi.bisquare), weights in the internal call to nls
> can become zero. Example:
d <- data.frame(x=1:5,y=c(2,3,5,10,9)) d.nlrob <- nlrob(y ~ k*x,start=list(k=1),data=d,psi=psi.bisquare)
I think the problem is as follows: After the call to nls, a weighted residual vector is calculated by dividing by sqrt(w). This generates NaN's for weights of zero, which leads to problems in the subsequent call to nls during the next iteration:
for (iiter in 1:maxit) {
...
w <- psi(resid/Scale, ...)
...
data$w <- sqrt(w)
...
out <- nls( ....., data=data ....... )
...
resid <- -residuals(out)/sqrt(w) # NaN for w=0
...
}
I wonder whether this problem shouldn't be dealt with by setting 'w' to 0 for the NaN cases.
I can get a fit by calling nlrob with na.action=na.exclude, but I'd have intuitively assumed that na.action applies to the NAs in the data set passed to nlrob, not to the one internally generated by nlrob and passed to nls.
You are right. The next version of robustbase (0.8-1) will have this fixed. Note however, that for me, in your example and in other more interesting ones, as soon as I use a redescending psi() function, the IRLS iterations do *not* converge... but rather strangely ``cycle'' around the true solution. As a redescending psi() has a non-convex rho(), and non-linear problems depend quite a bit on "feasible"/"good" starting values, I currently think I would discourage from using such psi(). As others, some possibly more expert in robust non-linear fitting, may be interested, and have interesting feedback, I'm crossposting this reply to the R-SIG-robust mailing list. (Further replies should typically only go there!)
The second 'issue' is that the weights are passed as 'w', whereas the documentation of 'nls' says weights should be given as 'weights' :
data: an optional data frame in which to evaluate the variables in
?formula? and ?weights?. Can also be a list or an
environment, but not a matrix.
I think it would be good to mention in the documentation of 'nls' that weights can be passed as 'w' as well.
You have overlooked that nlrob uses a different *formula* for nls() --- nowadays unnecessarily --- and that formula has a 'w'. So this part of your report is not correct. Martin Maechler, ETH Zurich
Pascal Niklaus