Re: [R] Closed form for regression splines - solution

From: Stephen A Roberts <Steve.Roberts_at_manchester.ac.uk>
Date: Thu 08 Dec 2005 - 03:08:19 EST

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@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)
}



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Dec 08 03:14:48 2005

This archive was generated by hypermail 2.1.8 : Thu 08 Dec 2005 - 05:31:26 EST