[R] nlrob problem

[Martin Maechler]

>>>>> 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