and their integrals (except C-spliness) and derivatives of given order by close-form recursive formulas. Compared to the splines package, splines2 allows piecewise constant basis for B-splines and provides a more user-friendly interface for their derivatives with consistent handling on NA’s. Most of the implementations had been (re)written in C++ with the help of Rcpp and RcppArmadillo since v0.3.0, which boosted the computational performance.
In the remaining of this vignette, we introduce the basic usage of most functions in the package through examples. The details of function syntax are available in the package manual and thus will not be discussed.
B-splines with their integrals and derivatives
Function bSpline() provides B-spline basis matrix and allows degree = 0 for piece-wise constant bases, which extends the bs() function in package splines with a better computational performance. One example of linear B-splines with two internal knots is given as follows:
library(splines2)
knots <- c(0.3, 0.5, 0.6)
x <- seq(0, 1, 0.01)
bsMat <- bSpline(x, knots = knots, degree = 1, intercept = TRUE)
matplot(x, bsMat, type = "l", ylab = "y")
abline(v = knots, lty = 2, col = "gray")
The close-form recursive formula of B-spline integrals and derivatives given by De Boor (1978) is implemented in function ibs() and dbs(), respectively. Two toy examples are given as follows:
ibsMat <- ibs(x, knots = knots, degree = 1, intercept = TRUE)
par(mfrow = c(1, 2))
matplot(x, bsMat, type = "l", ylab = "y")
abline(v = knots, h = 1, lty = 2, col = "gray")
matplot(x, ibsMat, type = "l", ylab = "y")
abline(v = knots, h = c(0.15, 0.2, 0.25), lty = 2, col = "gray")
bsMat <- bSpline(x, knots = knots, intercept = TRUE)
dbsMat <- dbs(x, knots = knots, intercept = TRUE)
par(mfrow = c(1, 2))
matplot(x, bsMat, type = "l", ylab = "y")
abline(v = knots, lty = 2, col = "gray")
matplot(x, dbsMat, type = "l", ylab = "y")
abline(v = knots, lty = 2, col = "gray")
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)))
M-splines using mSpline()
M-splines (Ramsay 1988) can be considered as a 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 boundary knots by default are the range of the data x, thus 0 and 1 in this example.
msMat <- mSpline(x, knots = knots, degree = 2, intercept = TRUE)
matplot(x, msMat, type = "l", ylab = "y")
abline(v = knots, lty = 2, col = "gray")
The derivative of given order of M-splines can be obtained by specifying a positive integer to argument dervis of mSpline(). Also, for an existing mSpline object generated by mSpline(), the deriv() method can be used conveniently. For example, the first derivative of the M-splines given in last 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 function mSpline() produces periodic spline basis (based on M-splines) when periodic = TRUE is specified. The Boundary.knots defines the start and end points of the cyclic period.
x1 <- seq.int(0, 3, 0.01)
pmsMat <- mSpline(x1, knots = knots, degree = 3, intercept = TRUE,
periodic = TRUE, Boundary.knots = c(0, 1))
matplot(x1, pmsMat, type = "l", xlab = "x", ylab = "Periodic Bases")
abline(v = seq.int(0, 3), lty = 2, col = "gray")
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)
matplot(x1, dpmsMat, type = "l", xlab = "x", ylab = "The 1st derivatives")
abline(v = seq.int(0, 3), lty = 2, col = "gray")
Furthermore, we can obtain the integrals of the periodic M-spline bases by specifying integral = TRUE. The integral is defined to be integrated from the left boundary knot.
ipmsMat <- mSpline(x1, knots = knots, degree = 3, intercept = TRUE,
periodic = TRUE, Boundary.knots = c(0, 1), integral = TRUE)
matplot(x1, ipmsMat, type = "l", xlab = "x", ylab = "Integrals")
abline(v = seq.int(0, 3), h = seq.int(0, 3), lty = 2, col = "gray")
I-splines using iSpline()
I-splines (Ramsay 1988) are simply the integral of M-splines and thus monotonically non-decreasing with unit maximum value. A monotonically non-decreasing (non-increasing) function can be fitted by a linear combination of I-spline bases with non-negative (non-positive) coefficients plus a constant, where the coefficient of the constant is unconstrained.
The example given by Ramsay (1988) was the I-splines corresponding to that quadratic M-splines with three internal knots placed at 0.3, 0.5, and 0.6. Notice that the degree of I-splines is defined from the associated M-splines instead of their own polynomial degree.
isMat <- iSpline(x, knots = knots, degree = 2, intercept = TRUE)
matplot(x, isMat, type = "l", ylab = "y")
abline(h = 1, v = knots, lty = 2, col = "gray")
The corresponding M-spline basis matrix can be obtained easily as the first derivatives of the I-splines by the deriv() method.
stopifnot(is_equivalent(msMat, deriv(isMat)))
We may specify the derivs = 2 in the deriv() method for the second derivatives of the I-splines, which are equivalent to the first derivatives of the corresponding M-splines.
dmsMat3 <- deriv(isMat, 2)
stopifnot(is_equivalent(dmsMat1, dmsMat3))
C-splines using cSpline
Convex splines (Meyer 2008) called C-splines are scaled integrals of I-splines with unit maximum value at the right boundary knot. Meyer (2008) applied C-splines to shape-restricted regression analysis. The monotone (non-decreasing) property of I-spines ensures the convexity of C-splines. A convex regression function can be estimated using linear combinations of the C-spline bases with non-negative coefficients, plus an unconstrained linear combination of a constant and an identity function \(g(x)=x\). If the underlying regression function is both increasing and convex, the coefficient on the identity function is restricted to be non-negative as well.
We may specify the argument scale = FALSE in function cSpline() to disable the scaling of the integrals of I-splines. Then the actual integrals of the corresponding I-spl
