[R] Bug in lowess

[Frank E Harrell Jr]

Prof Brian Ripley wrote:
> Frank Harrell wrote:
> 
> [...]
> 
>> Thank you Brian.  It seems that no matter what is the right answer, the
>> answer currently returned on my system is clearly wrong.  lowess()$y
>> should be constrained to be within range(y).
> 
> Really?  That assertion is offered without proof and I am pretty sure is 
> incorrect.  Consider
> 
>> x <- c(1:10, 20)
>> y <- c(1:10, 5) + 0.01*rnorm(11)
>> lowess(x,y)
> $x
>  [1]  1  2  3  4  5  6  7  8  9 10 20
> 
> $y
>  [1]  0.9983192  1.9969599  2.9960805  3.9948224  4.9944158  5.9959855
>  [7]  6.9964400  7.9981434  8.9990607 10.0002567 19.9946117
> 
> Remember that lowess is a local *linear* fitting method, and may give 
> zero weight to some data points, so (as here) it can extrapolate.

Brian - thanks - that's a good example though not typical of the kind I 
see from patients.

> 
> After reading what src/appl/lowess.doc says should happen with zero 
> weights, I think the answer given on Frank's system probably is the 
> correct one.  Rounding error is determining which of the last two points 
> is given zero robustness weight: on a i686 system both of the last two 
> are, and on mine only the last is. As far as I can tell in 
> infinite-precision arithmetic both would be zero, and hence the value at 
> x=120 would be found by extrapolation from those (far) to the left of it.
> 
> I am inclined to think that the best course of action is to quit with a 
> warning when the MAD of the residuals is effectively zero.  However, we 
> need to be careful not to call things 'bugs' that we do not understand 
> well enough.  This might be a design error in lowess, but it is not 
> AFAICS a bug in the implementation.

Yes it appears to be a weakness in the underlying algorithm.

Thanks
Frank