[R] Closed form for regression splines - solution

[Stephen A Roberts]

Greetings,

The question was about getting closed form equations from a bs() representation of a curve so enabling publication of a fitted curve for uise outside R. In case anyone else is interested, the solution lies in exploiting the fact that we can get accurate derivatives at the breakpoints and from these reconstruct the individual polynomial pieces.

Prototype function below - maybe the basis for a useful addition to the splines library?

Steve.

###########################################################################
require(splines)

# pieces
#
# Prototype function to compute the coefficients of the individual 
# polynomial pieces making up a cubic spline
# Works for bs() and ns() with intercept=FALSE
#
# Based on the observation that the derivatives at the breakpoints can be
# reliably computed using splineDesign and these define the functions.
# See De Boor SIAM J Num Anal 14 441-472 (1977) 
#
# Arguments:
#  call: spline function call used in modelling function (eg bs(x,4)) could 
#        be taken from attr(fit$terms,"predvars")
#  data: frame to evaluate call in - usually data argument of modelling 
#        function call
#  coef: spline function coefficients (no intercept assumed!)
#
# Notes:
#  Currently assumes intercept=FALSE used in bs and ns call - the usual for #  modelling use
#  - relatively trivial to modify to allow this, but messy.
#
# Future:
#  Needs bulletproofing and deeper testing and code could be more elegant
#  Would like to add  a text representation or maybe an R function.
#
# Value:
#  A list with components
#       npieces: no of pieces
#       degree: degree of fitted polynomial pieces
#       cutpoints: the ranges of the pieces: 
#                  piece i covers cutpoint[i] to cutpoint[i+1]
#       pars: A list with one vector per piece giving coefficients of each                         #             piece for 1,(x-c),(x-c)^2.... 
#               where c=cutpoint[i] for the ith piece.
#
# Steve Roberts Dec 2005 
# steve.Roberts at manchester.ac.uk
############################################################################

pieces <- function(call,coef,data=NULL)
{
    basis <- eval(call,list(data,sys.parent))
    degree <- attr(basis,"degree")
    
    #knot positions as in bs and ns
    if (class(basis)[1]=="bs")     
        Aknots <- sort(c(rep(attr(basis,"Boundary.knots"), degree+1), 
                  attr(basis,"knots"))) 
    else if (class(basis)[1]=="ns") 
        Aknots <- sort(c(rep(attr(basis,"Boundary.knots"), 4), 
                   attr(basis,"knots")))
    else 
        stop( paste("Unrecognised spline type", class(basis)[1]))
    
    cutpoints <- sort(unique(c(min(attr(basis,"Boundary.knots")),
                 attr(basis,"knots"))))
    npieces<-length(cutpoints)
    
    #evaluate derivatives and collect in array dd[order,piece]
    dd<-matrix(NA,degree+1,length(cutpoints))
    if (class(basis)[1]=="bs")
    {
    for (j in 0:degree)
    dd[j+1,] <- splineDesign(Aknots,cutpoints,ord=degree+1,outer=T, 
                deriv=rep(j,npieces))[,-1,drop=F] %*% coef
    }
    else if (class(basis)[1]=="ns") 
    {
    
    for (j in 0:degree)
    {#Code from ns()
        sdes <- splineDesign(Aknots, cutpoints, 4,outer=T, 
                deriv=rep(j,npieces))[, -1, drop = FALSE]
        const <- splineDesign(Aknots, attr(basis,"Boundary.knots"), 4, 
                 c(2, 2))[, -1, drop = FALSE]
        qr.const <- qr(t(const))
        sdes <- as.matrix((t(qr.qty(qr.const, t(sdes))))[, -(1:2)])
        dd[j+1,] <- sdes %*% coef
    }   
    }
    
    #add upper boundary knot to output list for convenience
    cutpoints <- c(cutpoints,max(attr(basis,"Boundary.knots") ))
    
    #create output coefficient lists
    pars <- list()
    for ( i in 1:npieces)
    {
    pars[[i]] <- dd[,i]/factorial(0:degree)
    names(pars[[i]]) <- c("1",paste("x^",1:degree,sep=""))
    names(pars)[i] <- paste(cutpoints[i],"to",cutpoints[i+1])
    }
    
    list(npieces=npieces,degree=degree,cutpoints=cutpoints,pars=pars)
}

########################
# example

x <- seq(10,30,len=10)
y <- sin(x/2)+1
fit.bs <- lm( y~bs(x,5) )
fit.ns <- lm( y~ns(x,5) )
coef.bs <- coef(fit.bs)[-1]
coef.ns <- coef(fit.ns)[-1]

pieces.bs <- pieces(attr(fit.bs$terms,"predvars")[[3]],coef.bs)
pieces.ns <- pieces(attr(fit.ns$terms,"predvars")[[3]],coef.ns)

#verify in plot
plot(y~x)
xx <- seq(min(x),max(x),len=3000)
lines(predict(fit.ns,new=data.frame(x=xx))~xx, lty=2)
lines(predict(fit.bs,new=data.frame(x=xx))~xx, lty=1)
abline(v=pieces.bs$cutpoints)
abline(v=pieces.ns$cutpoints,lty=2)

for (i in 1:pieces.bs$npieces)
{
xx <- seq(pieces.bs$cutpoints[i],pieces.bs$cutpoints[i+1],len=1000)
mm <- matrix(NA,length(xx),pieces.bs$degree+1)
for (j in 0:pieces.bs$degree) mm[,j+1] <- (xx-pieces.bs$cutpoints[i])^j 
yy <- coef(fit.bs)[1] + mm %*% pieces.bs$pars[[i]]
lines(yy~xx, col=i+1)
}
for (i in 1:pieces.ns$npieces)
{
xx <- seq(pieces.ns$cutpoints[i],pieces.ns$cutpoints[i+1],len=1000)
mm <- matrix(NA,length(xx),pieces.ns$degree+1)
for (j in 0:pieces.ns$degree) mm[,j+1] <- (xx-pieces.ns$cutpoints[i])^j 
yy <- coef(fit.ns)[1] + mm %*% pieces.ns$pars[[i]]
lines(yy~xx, col=i+1, lty=2)
}