[R] predict(backSpline(x)): losing my marbles?

[Greg Snow]

Ben, 

I am pretty sure that you have found a bug, I am just not sure where exactly the bug would be.

The problem comes from the fact that your bsp object has the knots in decreasing order, if we reverse the order like this:

bsp3 <- bsp
bsp3$knots <- rev(bsp$knots)
bsp3$coefficients <- bsp$coefficients[20:1,]

lines(predict(bsp3,seq(0,6,length=101)), col='red')

Then the prediction looks like it should.

If we trace through the predict.polySpline function with x=5.6 as one example, then inside the prediction function the value of i (the interval between knots that x belongs to) is 18, well if the knots are in increasing order then that is the correct interval, this comes from the line: i <- as.numeric(cut(x, knots)) and the cut function treats the knots as being in increasing order, not decreasing, so that in a later step when the function is supposed to be finding how far x is from the closest knot below it (which should be 5.2366151 for your example) it grabs the 18 element of knots, which in its decending order is instead: 0.2873583, so the difference ends up being much too big, when the function squares and cubes this much too big value, we get values that are way off (as you noticed), the wrong set of coefficients is being used as well (also row 18, which correspond to 0.28 not 5.23).  I would say that the most surprising thing with your example is the range at which the prediction is actually decent, I would only expect this between the middle 2 knots.

I see 2 possible ways to fix the problem:

1. The predict function checks the order of the knots and either reorders the knots and coefficients when the knots are descending, or does the correct computations of the interval.

2. The backSpline function checks and makes sure that it always produces the results with the knots in ascending order (and the coefficients to match).

Sorry, but since it looks like you are not the author of backSpline, I don't think that you can claim the boneheadedness this time, your approach was reasonable, and you were just a victim of the reversed knots.

Unless I have missed something here (please double check), you should submit a bug report (unless someone fixes this before you get around to that). 

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> project.org] On Behalf Of Ben Bolker
> Sent: Tuesday, September 14, 2010 10:15 AM
> To: r-help at r-project.org
> Subject: [R] predict(backSpline(x)): losing my marbles?
> 
> 
>   I'm sure I'm doing something completely boneheaded here, but I've
> used this idiom
> (constructing an interpolation spline and using prediction from a
> backSpline to find
> an approximation profile confidence interval) many times before and
> haven't hit this
> particular problem:
> 
>   r2 <- c(1.04409027570601, 1.09953936543359, 1.15498845516117,
> 1.21043754488875,
>           1.26588663461633, 1.32133572434391, 1.37678481407149,
> 1.43223390379907,
>           1.48768299352664, 1.54313208325422, 1.5985811729818,
> 1.65403026270938,
>           1.70947935243696, 1.76492844216454, 1.82037753189212,
> 1.87582662161970,
>           1.93127571134727, 1.98672480107485, 2.04217389080243,
> 2.09762298053001
>           )
>   d2 <- c(6.1610616585333, 5.70079375491741, 5.2366151167289,
> 4.77263065800071,
>           4.31310259797181, 3.86232922249189, 3.42452047126494,
> 3.00367670365119,
>           2.6034766331926, 2.22717964214416, 1.87754657382891,
> 1.55678176465949,
>           1.26649764837839, 1.00770187223770, 0.780805622450771,
> 0.585650849661306,
>           0.421553364080296, 0.287358347766713, 0.18150469048976,
> 0.102094654969619
>           )
>   plot(d2,r2,type="b")
>   require(splines)
>   sp <- interpSpline(r2,d2)
>   psp <- predict(sp)
>   points(psp$y,psp$x,col=5)
>   bsp <- backSpline(sp)
>   lines(predict(bsp,seq(0,6,length=101)),col=2)
> 
>    The prediction from the spline (cyan dots) looks perfectly
> reasonable.
> The prediction from the inverted spline matches the curve over part of
> the range but
> goes crazy elsewhere.  I would have expected it to be reasonably close
> to this well-behaved
> curve over the whole range.
> 
>   I have looked at the docs for interpSpline, backSpline,
> predict.bSpline, ... and nothing jumps out at me.
> 
>    Please be gentle, if possible, in explaining to me what I'm missing
> ...
> 
>   thanks,
>     Ben Bolker
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
> and provide commented, minimal, self-contained, reproducible code.