On 08/02/2012 05:00 AM, r-devel-request at r-project.org wrote:> Now I just have to grovel over the R code in ns() and bs() to figure > out how exactly they pick knots and handle boundary conditions, plus > there is some code that I don't understand in ns() that uses qr() to > postprocess the output from spline.des. I assume this is involved > somehow in imposing the boundary conditions... > > Thanks again everyone for your help, > -- NathanielThe ns and bs function post-process the spline bases to get an orthagonal basis matrix, this is the use of qr. I think this causes much more grief than it is worth, for the sake of a small increase in numeric stability. For instance when you plot the spline bases, they don't look anything like the basis functions one would expect. (Perhaps my background in numerical analysis was a hindrance here, since I know something about splines and thus have an expectation). Terry Therneau
On Thu, Aug 2, 2012 at 2:09 PM, Terry Therneau <therneau at mayo.edu> wrote:> > > On 08/02/2012 05:00 AM, r-devel-request at r-project.org wrote: >> >> Now I just have to grovel over the R code in ns() and bs() to figure >> out how exactly they pick knots and handle boundary conditions, plus >> there is some code that I don't understand in ns() that uses qr() to >> postprocess the output from spline.des. I assume this is involved >> somehow in imposing the boundary conditions... >> >> Thanks again everyone for your help, >> -- Nathaniel > > The ns and bs function post-process the spline bases to get an orthagonal > basis matrix, this is the use of qr. I think this causes much more grief > than it is worth, for the sake of a small increase in numeric stability. > For instance when you plot the spline bases, they don't look anything like > the basis functions one would expect. (Perhaps my background in numerical > analysis was a hindrance here, since I know something about splines and thus > have an expectation).You know, the white book does claim that S orthogonalizes the spline basis functions, but R doesn't seem to actually do that:> crossprod(bs(1:10, 4, intercept=TRUE))1 2 3 4 1 1.84104162 0.6154211 0.2570069 0.06430817 2 0.61542109 0.7710207 0.5787736 0.25700689 3 0.25700689 0.5787736 0.7710207 0.61542109 4 0.06430817 0.2570069 0.6154211 1.84104162> crossprod(ns(1:10, 4, intercept=TRUE))1 2 3 4 1 1.0659929 0.5314089 0.37125505 -0.17838443 2 0.5314089 1.0681360 0.23122018 0.14422236 3 0.3712550 0.2312202 1.30226657 -0.03809082 4 -0.1783844 0.1442224 -0.03809082 1.10144173 and they look quite reasonable when plotted to me. ns is doing something more complicated; I think it's computing the second derivatives of each basis function at the position of the boundary knots, and then using that to somehow transform the original basis. Since the key feature of ns is that the returned basis has a zero second derivative at the boundary knots, I'm sure this makes sense if one stares at it for long enough. -n
If R's bs() and ns() are like S+'s (they do give very similar results* and S+'s were written by Doug Bates), then bs() does not do any linear algebra (like qr()) on splineDesign's output. bs() needs to come up with a default set of knots (if the user didn't supply them), combine the Boundary.knots (repeated ord times) and knots into one vector to pass to splineDesign, and, if there are any x's outside of range(Boundary.knots), treat them specially. ns() needs to project splineDesign's basis functions onto the subspace of functions whose 2nd derivatives are zero at the endpoints. The projection can be done with qr() and friends (S+ uses qr()). (If you want basis functions for a periodic spline you could use a similar procedure to project to the subspace of functions whose first and second derivatives match at the upper and lower Boundary.knots.) * The only difference between the R and S+ versions of bs() that I've noticed is in how they deal with the x's that are outside of range(Boundary.knots). S+ extrapolates with cubics both above and below that range while R extrapolates with a cubic below the range and with a quadratic above the range. I don't know what the rationale for this is. Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com> -----Original Message----- > From: r-devel-bounces at r-project.org [mailto:r-devel-bounces at r-project.org] On Behalf > Of Terry Therneau > Sent: Thursday, August 02, 2012 6:10 AM > To: r-devel at r-project.org; Nathaniel Smith > Subject: Re: [Rd] Rd] Numerics behind splineDesign > > > > On 08/02/2012 05:00 AM, r-devel-request at r-project.org wrote: > > Now I just have to grovel over the R code in ns() and bs() to figure > > out how exactly they pick knots and handle boundary conditions, plus > > there is some code that I don't understand in ns() that uses qr() to > > postprocess the output from spline.des. I assume this is involved > > somehow in imposing the boundary conditions... > > > > Thanks again everyone for your help, > > -- Nathaniel > The ns and bs function post-process the spline bases to get an > orthagonal basis matrix, this is the use of qr. I think this causes > much more grief than it is worth, for the sake of a small increase in > numeric stability. For instance when you plot the spline bases, they > don't look anything like the basis functions one would expect. (Perhaps > my background in numerical analysis was a hindrance here, since I know > something about splines and thus have an expectation). > > Terry Therneau > > ______________________________________________ > R-devel at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel