A Short Introduction to splines2

B-spline Basis Functions

The bSpline() function generates the basis matrix for B-splines and extends the function bs() of the package splines by providing 1) the piece-wise constant basis functions when degree = 0, 2) the derivatives of basis functions for a positive derivs, 3) the integrals of basis functions if integral = TRUE, 4) periodic basis functions based on B-splines if periodic = TRUE.

One example of linear B-splines with three internal knots is as follows:

knots <- c(0.3, 0.5, 0.6)
x <- seq(0, 1, 0.01)
bsMat <- bSpline(x, knots = knots, degree = 1, intercept = TRUE)
plot(bsMat, mark_knots = "all")

B-splines of degree one with three internal knots placed at 0.3, 0.5, and 0.6.

Integrals and Derivatives of B-splines

For convenience, the package also provides functions ibs() and dbs() for constructing the B-spline integrals and derivatives, respectively. Two toy examples are as follows:

ibsMat <- ibs(x, knots = knots, degree = 1, intercept = TRUE)
op <- par(mfrow = c(1, 2))
plot(bsMat, mark_knots = "internal")
plot(ibsMat, mark_knots = "internal")
abline(h = c(0.15, 0.2, 0.25), lty = 2, col = "gray")

Piecewise linear B-splines (left) and their integrals (right).

bsMat <- bSpline(x, knots = knots, intercept = TRUE)
dbsMat <- dbs(x, knots = knots, intercept = TRUE)
plot(bsMat, mark_knots = "internal")
plot(dbsMat, mark_knots = "internal")

Cubic B-spline (left) and their first derivative (right).

We may also obtain the derivatives easily by the deriv() method as follows:

is_equivalent <- function(a, b) {
    all.equal(a, b, check.attributes = FALSE)
}
stopifnot(is_equivalent(dbsMat, deriv(bsMat)))

Periodic B-splines

The function bSpline() produces periodic spline basis functions following Piegl and Tiller (1997, chap. 12) when periodic = TRUE is specified. Different from the regular basis functions, the x is allowed to be placed outside the boundary and the Boundary.knots defines the cyclic interval. For instance, one may obtain the periodic cubic B-spline basis functions with cyclic interval (0, 1) as follows:

px <- seq(0, 3, 0.01)
pbsMat <- bSpline(px, knots = knots, Boundary.knots = c(0, 1),
                  intercept = TRUE, periodic = TRUE)
ipMat <- ibs(px, knots = knots, Boundary.knots = c(0, 1),
             intercept = TRUE, periodic = TRUE)
dp1Mat <- deriv(pbsMat)
dp2Mat <- deriv(pbsMat, derivs = 2)
par(mfrow = c(1, 2))
plot(pbsMat, ylab = "Periodic B-splines", mark_knots = "boundary")
plot(ipMat, ylab = "The integrals", mark_knots = "boundary")

plot(dp1Mat, ylab = "The 1st derivatives", mark_knots = "boundary")
plot(dp2Mat, ylab = "The 2nd derivatives", mark_knots = "boundary")

For reference, the corresponding integrals and derivatives are also visualized.

M-spline Basis Functions

M-splines (Ramsay 1988) can be considered the normalized version of B-splines with unit integral within boundary knots. An example given by Ramsay (1988) was a quadratic M-splines with three internal knots placed at 0.3, 0.5, and 0.6. The default boundary knots are the range of x, and thus 0 and 1 in this example.

msMat <- mSpline(x, knots = knots, degree = 2, intercept = TRUE)
par(op)
plot(msMat, mark_knots = "all")

Quadratic M-spline with three internal knots placed at 0.3, 0.5, and 0.6.

The derivative of the given order of M-splines can be obtained by specifying a positive integer to argument dervis of mSpline(). Similarly, for an existing mSpline object generated by mSpline(), one can use the deriv() method for derivaitives. For example, the first derivative of the M-splines given in the previous example can be obtained equivalently as follows:

dmsMat1 <- mSpline(x, knots = knots, degree = 2, intercept = TRUE, derivs = 1)
dmsMat2 <- deriv(msMat)
stopifnot(is_equivalent(dmsMat1, dmsMat2))

Periodic M-Splines

The mSpline() function produces periodic splines based on M-spline basis functions when periodic = TRUE is specified. The Boundary.knots defines the cyclic interval, which is the same with the periodic B-splines.

pmsMat <- mSpline(px, knots = knots, intercept = TRUE,
                  periodic = TRUE, Boundary.knots = c(0, 1))
plot(pmsMat, ylab = "Periodic Basis", mark_knots = "boundary")

Cubic periodic M-splines.

We may still specify the argument derivs in mSpline() or use the corresponding deriv() method to obtain the derivatives when periodic = TRUE.

dpmsMat <- deriv(pmsMat)
plot(dpmsMat, ylab = "The 1st derivatives", mark_knots = "boundary")

The first derivatives of the periodic M-splines.

Furthermore, we can obtain the integrals of the periodic M-splines by specifying integral = TRUE. The integral is integrated from the left boundary knot.

ipmsMat <- mSpline(px, knots = knots, intercept = TRUE,
                   periodic = TRUE, Boundary.knots = c(0, 1), integral = TRUE)
plot(ipmsMat, ylab = "Integrals", mark_knots = "boundary")
abline(h = seq.int(0, 3), lty = 2, col = "gray")

The integrals of the periodic M-splines.

Nonnegative Natural Cubic Basis Functions

The package provides two variants of the natural cubic splines that can be constructed by naturalSpline() and nsk(), respectively, both of which are different from splines::ns().

The naturalSpline() function returns nonnegative basis functions (within the boundary) for natural cubic splines by utilizing a closed-form null space derived from the second derivatives of cubic B-splines. When integral = TRUE, the function naturalSpline() returns the integral of each natural spline basis.

nsMat <- naturalSpline(x, knots = knots, intercept = TRUE)
insMat <- naturalSpline(x, knots = knots, intercept = TRUE, integral = TRUE)
par(mfrow = c(1, 2))
plot(nsMat, ylab = "Basis")
plot(insMat, ylab = "Integrals")

Nonnegative natural cubic splines (left) and corresponding integrals (right).

stopifnot(is_equivalent(nsMat, deriv(insMat)))

Similarly, one may directly specify the argument derivs in naturalSpline() or use the corresponding deriv() method to obtain the derivatives of spline basis functions.

d1nsMat <- naturalSpline(x, knots = knots, intercept = TRUE, derivs = 1)
d2nsMat <- deriv(nsMat, 2)
matplot(x, d1nsMat, type = "l", ylab = "The 1st derivatives")
matplot(x, d2nsMat, type = "l", ylab = "The 2nd derivatives")

The derivatives of natural cubic splines.

Natural Cubic Basis Functions with Unit Heights at Knots

The function nsk() produces another variant of natural cubic splines, where only one of the spline basis functions is nonzero with unit height at every boundary and internal knot. As a result, the coefficients of the basis functions are the values of the spline function at the knots, which makes it more straightforward to interpret the coefficient estimates. This idea originated from the function nsk() of the survival package (introduced in version 3.2-8). The implementation of the nsk() of splines2 essentially follows the survival::nsk() function. One noticeable argument for nsk() is trim (equivalent to the argument b of survival::nsk()). One may specify trim = 0.05 to exclude 5% of the data from both sides when setting the boundary knots, which can be a more sensible choice in practice due to possible outliers. The trim argument is also available for naturalSpline(), which however is zero by default for backward compatibility. An illustration of the basis functions generated by nsk() is as follows:

nskMat <- nsk(x, knots = knots, intercept = TRUE)
par(op)
plot(nskMat, ylab = "nsk()", mark_knots = "all")
abline(h = 1, col = "red", lty = 3)

We can visually verify that only one basis function takes a value of one at each knot.

Update Spline’s Specification by update()

The update() function is an S3 method to generate spline basis functions with new x, degree, knots, or df specifications. The first argument is an existing splines2 object and additional named arguments will be passed to the corresponding functions to update the spline basis functions. Suppose we want to add two more knots to nskMat for natural cubic spline basis functions and exclude 5% of the data from both sides in total when placing the boundary knots. We can utilize the update() method as follows:

nskMat2 <- update(nskMat, knots = c(knots, 0.2, 0.4), trim = 0.025)
knots(nskMat2)

## [1] 0.2 0.3 0.4 0.5 0.6

stopifnot(all.equal(quantile(x, c(0.025, 0.975), names = FALSE),
                    knots(nskMat2, "boundary")))

Evaluation by predict()

The predict() method for splines2 objects allows one to evaluate the spline function if a coefficient vector is specified via the coef argument. In addition, it internally calls the update() method to update the basis functions before computing the spline function, which can be useful to get the derivatives of the spline function. If the coef argument is not specified, the predict() method will be equivalent to the update() method. For instance, we can compute the first derivative of the I-spline function from the previous example with a coefficient vector seq(0.1, by = 0.1, length.out = ncol(isMat)) at \(x = (0.275, 0.525, 0.8)\) as follows:

new_x <- c(0.275, 0.525, 0.8)
names(new_x) <- paste0("x=", new_x)
(isMat2 <- predict(isMat, newx = new_x)) # the basis functions at new x

##                 1         2         3         4        5     6
## x=0.275 0.9994213 0.7730556 0.2310764 0.0000000 0.000000 0.000
## x=0.525 1.0000000 1.0000000 0.9765625 0.2696429 0.000625 0.000
## x=0.8   1.0000000 1.0000000 1.0000000 0.9428571 0.580000 0.125

stopifnot(all.equal(predict(isMat, newx = new_x), update(isMat, x = new_x)))
## compute the first derivative of the I-spline function in different ways
beta <- seq(0.1, by = 0.1, length.out = ncol(isMat))
deriv_ispline1 <- predict(isMat, newx = new_x, coef = beta, derivs = 1)
deriv_ispline2 <- predict(update(isMat, x = new_x, derivs = 1), coef = beta)
deriv_ispline3 <- c(predict(deriv(isMat), newx = new_x) %*% beta)
stopifnot(all.equal(deriv_ispline1, deriv_ispline2))
stopifnot(all.equal(deriv_ispline2, deriv_ispline3))

Visualization by plot()

As one may notice in the previous examples, we may visualize the spline basis functions easily with the plot() method. By default, the spline basis functions are visualized at 101 equidistant grid points within the range of x, which can be tweaked by arguments from, to, and n. In addition, we can visualize the placement of knots by vertical lines through the argument mark_knots. The available options are "none", "internal", "boundary", and "all".

Including Spline Basis Functions in Model Formulas

It is common to directly include spline basis functions in a model formula. To avoid a lengthy model formula, the package provides alias functions that are summarized in the following table:

One may create new alias functions. For example, we can create a new alias function simply named b() for B-splines and obtain equivalent models as follows:

b <- bSpline # create an alias for B-splines
mod1 <- lm(weight ~ b(height, degree = 1, df = 3), data = women)
iknots <- with(women, knots(bSpline(height, degree = 1, df = 3)))
mod2 <- lm(weight ~ bSpline(height, degree = 1, knots = iknots), data = women)
pred1 <- predict(mod1, head(women, 10))
pred2 <- predict(mod2, head(women, 10))
stopifnot(all.equal(pred1, pred2))

Nevertheless, there is a possible pitfall when using a customized wrapper function for spline basis functions along with a data-dependent placement of knots. When we make model predictions for a given new data, the placement of the internal/boundary knots can be different from the original placement that depends on the training set. As a result, the spline basis functions generated for prediction may not be the same as the counterparts used in the model fitting. A simple example is as follows:

## generates quadratic spline basis functions based on log(x)
qbs <- function(x, ...) {
    splines2::bSpline(log(x), ..., degree = 2)
}
mod3 <- lm(weight ~ qbs(height, df = 5), data = women)
mod4 <- lm(weight ~ bsp(log(height), degree = 2, df = 5), data = women)
stopifnot(all.equal(unname(coef(mod3)), unname(coef(mod4)))) # the same coef
pred3 <- predict(mod3, head(women, 10))
pred4 <- predict(mod4, head(women, 10))
all.equal(pred3, pred4)

## [1] "Mean relative difference: 0.07185939"

pred0 <- predict(qbs(women$height, df = 5),
                 newx = head(log(women$height), 10),
                 coef = coef(mod3)[- 1]) + coef(mod3)[1]
stopifnot(all.equal(pred0, pred4, check.names = FALSE))

Although the coefficient estimates are the same, the prediction results by using the predict.lm() differ. Using an alias function in the model formula produces correct results.

To resolve this issue, we can create an S3 method for stats::makepredictcall() as follows:

## generates quadratic spline basis functions based on log(x) with a new class
qbs <- function(x, ...) {
    res <- splines2::bSpline(log(x), ..., degree = 2)
    class(res) <- c("qbs", class(res))
    return(res)
}
## a utility to help model.frame() create the right matrices
makepredictcall.qbs <- function(var, call) {
    if (as.character(call)[1L] == "qbs" ||
    (is.call(call) && identical(eval(call[[1L]]), qbs))) {
        at <- attributes(var)[c("knots", "Boundary.knots", "intercept",
                                "periodic", "derivs", "integral")]
        call <- call[1L:2L]
        call[names(at)] <- at
    }
    call
}
## the same example
mod3 <- lm(weight ~ qbs(height, df = 5), data = women)
mod4 <- lm(weight ~ bsp(log(height), degree = 2, df = 5), data = women)
stopifnot(all.equal(unname(coef(mod3)), unname(coef(mod4)))) # the same coef
pred3 <- predict(mod3, head(women, 10))
pred4 <- predict(mod4, head(women, 10))
all.equal(pred3, pred4) # should be TRUE this time

## [1] TRUE

Extract Specifications by $

The basis specifications are saved as attributes of the returned splines2 objects, which means that we can extract one of the specifications by attr(). Alternatively, we can treat splines2 objects as lists and use the corresponding $ method. For example, it is straightforward to extract the specified trim of nskMat2 by attr(nskMat2, "trim") or simply nskMat2$trim.

c(nskMat2$trim, attr(nskMat2, "trim"))

## [1] 0.025 0.025