R help - Dec 2011 - nlrob problem

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.

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.

Pascal Niklaus

Please ignore my previous posting, in the meanwhile I have realised that I
did not look carefully enough how nlrob works (irls), although this is
pretty obvious from the code.

The only issue that remains is the case of zero weight for redescending M
estimators, in which case NaN's are produced. I think it would be better if
one could avoid the NaN's produced by dividing by zero when psi becomes
zero. 

Something along the line of replacing 

           resid <- -residuals(out)/sqrt(w)

by 

           resid <- ifelse(w>0,-residuals(out)/sqrt(w),0);

There probably is a more elegant way to do this. 

This would avoid that one has to use na.action=na.exclude for cases where
there are no NAs in the data set originally passed to nlrob (and the NAs
only occur 'internally' in the code of nlrob).

Hope I haven't missed something else this time...

Pascal Niklaus



--
View this message in context:
http://r.789695.n4.nabble.com/nlrob-problem-tp4215473p4229016.html
Sent from the R help mailing list archive at Nabble.com.

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