The **splines2** package (version 0.4.1) provides functions constructing a variety of regression spline bases that are not available from the **splines** package shipped with base R. Most functions have a very similar user interface with the function `splines::bs()`

. To be more specific, it provides functions to construct basis matrix of

- B-splines
- M-splines
- I-splines
- C-splines
- periodic M-splines
- natural cubic splines
- generalized Bernstein polynomials

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.

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)
<- c(0.3, 0.5, 0.6)
knots <- seq(0, 1, 0.01)
x <- bSpline(x, knots = knots, degree = 1, intercept = TRUE)
bsMat 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:

```
<- ibs(x, knots = knots, degree = 1, intercept = TRUE)
ibsMat 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")
```

```
<- bSpline(x, knots = knots, intercept = TRUE)
bsMat <- dbs(x, knots = knots, intercept = TRUE)
dbsMat 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:

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

`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.

```
<- mSpline(x, knots = knots, degree = 2, intercept = TRUE)
msMat 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:

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

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.

```
<- seq.int(0, 3, 0.01)
x1 <- mSpline(x1, knots = knots, degree = 3, intercept = TRUE,
pmsMat 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`

.

```
<- deriv(pmsMat)
dpmsMat 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.

```
<- mSpline(x1, knots = knots, degree = 3, intercept = TRUE,
ipmsMat 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")
```

`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.

```
<- iSpline(x, knots = knots, degree = 2, intercept = TRUE)
isMat 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.

```
<- deriv(isMat, 2)
dmsMat3 stopifnot(is_equivalent(dmsMat1, dmsMat3))
```

`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-splines will be returned. If `scale = TRUE`

(by default), each C-spline basis is scaled to have unit height at the right boundary knot.

```
<- cSpline(x, knots = knots, degree = 2, intercept = TRUE)
csMat1 matplot(x, csMat1, type = "l", ylab = "y")
abline(h = 1, v = knots, lty = 2, col = "gray")
```

Similarly, the `deriv()`

method can be used to obtain the derivatives. A nested call of `deriv()`

is supported for derivatives of order greater than one. However, the argument `derivs`

of the `deriv()`

method can be specified directly for a better computational performance. For example, the first and second derivatives can be obtained by the following equivalent approaches, respectively.

```
<- cSpline(x, knots = knots, degree = 2, intercept = TRUE, scale = FALSE)
csMat2 stopifnot(is_equivalent(isMat, deriv(csMat2)))
stopifnot(is_equivalent(msMat, deriv(csMat2, 2)))
stopifnot(is_equivalent(msMat, deriv(deriv(csMat2))))
```

The Bernstein polynomials have also been applied to shape-constrained regression analysis (Wang and Ghosh 2012). The \(i\)-th basis of the generalized Bernstein polynomials of degree \(n\) over \([a, b]\) is defined as follows: \[ B_i^n(x)=\frac{1}{(b-a)^n}{n\choose i}(x-a)^i (b-x)^{n-i},~i\in\{0,\ldots,n\}, \] where \(a\le x\le b\). Obviously, it reduces to regular Bernstein polynomials defined over \([0, 1]\) when \(a = 0\) and \(b = 1\).

We may obtain the basis matrix of the generalized using the function `bernsteinPoly()`

. For example, the Bernstein polynomials of degree 4 over \([0, 1]\) and is generated as follows:

```
<- seq.int(0, 1, 0.01)
x1 <- seq.int(- 1, 1, 0.01)
x2 <- bernsteinPoly(x1, degree = 4, intercept = TRUE)
bpMat1 <- bernsteinPoly(x2, degree = 4, intercept = TRUE)
bpMat2 par(mfrow = c(1, 2))
matplot(x1, bpMat1, type = "l", ylab = "y")
matplot(x2, bpMat2, type = "l", ylab = "y")
```

In addition, we may specify `integral = TRUE`

or `derivs = 1`

in `bernsteinPoly()`

for their integrals or first derivatives, respectively.

```
<- bernsteinPoly(x1, degree = 4, intercept = TRUE, integral = TRUE)
ibpMat1 <- bernsteinPoly(x2, degree = 4, intercept = TRUE, integral = TRUE)
ibpMat2 <- bernsteinPoly(x1, degree = 4, intercept = TRUE, derivs = 1)
dbpMat1 <- bernsteinPoly(x2, degree = 4, intercept = TRUE, derivs = 1)
dbpMat2 par(mfrow = c(2, 2))
matplot(x1, ibpMat1, type = "l", ylab = "Integrals")
matplot(x2, ibpMat2, type = "l", ylab = "y")
matplot(x1, dbpMat1, type = "l", ylab = "y")
matplot(x2, dbpMat2, type = "l", ylab = "y")
```

Similarly, we may also use the `deriv()`

method to get derivatives of an existing `bernsteinPoly`

object.

```
stopifnot(is_equivalent(dbpMat1, deriv(bpMat1)))
stopifnot(is_equivalent(dbpMat2, deriv(bpMat2)))
stopifnot(is_equivalent(dbpMat1, deriv(ibpMat1, 2)))
stopifnot(is_equivalent(dbpMat2, deriv(ibpMat2, 2)))
```

The function `naturalSpline()`

returns nonnegative bases (within boundary) for natural cubic splines by utilizing the close-form null space derived from the second derivatives of cubic B-splines. While `splines::ns()`

uses QR decomposition to find the null space of the second derivatives of B-spline basis at boundary knots with no guarantee that the resulting bases are nonnegative within boundary. When `integral = TRUE`

, `naturalSpline()`

returns integral of each natural spline basis.

```
<- naturalSpline(x, knots = knots, intercept = TRUE)
nsMat <- naturalSpline(x, knots = knots, intercept = TRUE, integral = TRUE)
insMat par(mfrow = c(1, 2))
matplot(x, nsMat, type = "l", ylab = "Bases")
matplot(x, insMat, type = "l", ylab = "Integrals")
```

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

Similar to `bernsteinPoly()`

, one may specify the argument `derivs`

in `naturalSpline()`

or use the corresponding `deriv()`

method to obtain the derivatives of spline bases.

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

`predict`

The methods for **splines2** objects dispatched by generic function `predict`

will be useful if we want to evaluate the spline object at possibly new \(x\) values. For instance, we may evaluate the value of I-splines object in previous example at 0.275, 0.525, and 0.8, respectively, as follows:

```
<- c(0.275, 0.525, 0.8)
new_x names(new_x) <- paste0("x=", new_x)
predict(isMat, 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
```

Technically speaking, the methods take all information needed, such as `knots`

, `degree`

, `intercept`

, etc., from attributes of the original objects and call the corresponding function automatically for those new \(x\) values. Therefore, the `predict`

methods will not be applicable if those attributes are somehow lost after some operations.

De Boor, Carl. 1978. *A Practical Guide to Splines*. Vol. 27. New York: springer-verlag. https://doi.org/10.1002/zamm.19800600129.

Meyer, Mary C. 2008. “Inference Using Shape-Restricted Regression Splines.” *The Annals of Applied Statistics* 2 (3): 1013–33. https://doi.org/10.1214/08-aoas167.

Ramsay, J. O. 1988. “Monotone Regression Splines in Action.” *Statistical Science*, 425–41. https://doi.org/10.1214/ss/1177012761.

Wang, Jiangdian, and Sujit K. Ghosh. 2012. “Shape Restricted Nonparametric Regression with Bernstein Polynomials.” *Computational Statistics & Data Analysis* 56 (9): 2729–41.