The R package **splines2** is intended to be a user-friendly *supplementary* package to the base package **splines**.

## Features

The package **splines2** provides functions to construct basis matrices of

- B-splines
- M-splines
- I-splines
- convex splines (C-splines)
- periodic splines
- natural cubic splines
- generalized Bernstein polynomials
- their integrals (except C-splines) and derivatives of given order by closed-form recursive formulas

In addition to the R interface, **splines2** provides a C++ header-only library integrated with **Rcpp**, which allows the construction of spline basis functions directly in C++ with the help of **Rcpp** and **RcppArmadillo**. Thus, it can also be treated as one of the **Rcpp*** packages. A toy example package that uses the C++ interface is available here.

## Installation of CRAN Version

You can install the released version from CRAN.

`install.packages("splines2")`

## Development

The latest version of the package is under development at GitHub. If it is able to pass the automated package checks, one may install it by

```
if (! require(remotes)) install.packages("remotes")
remotes::install_github("wenjie2wang/splines2", upgrade = "never")
```

## Getting Started

The Online document provides a reference for all functions and contains the following vignettes:

## Performance

Since v0.3.0, the implementation of the main functions has been rewritten in C++ with the help of the **Rcpp** and **RcppArmadillo** packages. The computational performance has thus been boosted and comparable with the function `splines::splineDesign()`

.

Some quick micro-benchmarks are provided for reference as follows:

```
library(microbenchmark)
options(microbenchmark.unit="relative")
library(splines)
library(splines2)
set.seed(123)
x <- runif(1e3)
degree <- 3
ord <- degree + 1
internal_knots <- seq.int(0.1, 0.9, 0.1)
boundary_knots <- c(0, 1)
all_knots <- sort(c(internal_knots, rep(boundary_knots, ord)))
## check equivalency of outputs
my_check <- function(values) {
all(sapply(values[- 1], function(x) {
all.equal(unclass(values[[1]]), x, check.attributes = FALSE)
}))
}
```

For B-splines, function `splines2::bSpline()`

provides equivalent results with `splines::bs()`

and `splines::splineDesign()`

, and is about 3x faster than `bs()`

and 2x faster than `splineDesign()`

for this example.

```
## B-splines
microbenchmark(
"splines::bs" = bs(x, knots = internal_knots, degree = degree,
intercept = TRUE, Boundary.knots = boundary_knots),
"splines::splineDesign" = splineDesign(x, knots = all_knots, ord = ord),
"splines2::bSpline" = bSpline(
x, knots = internal_knots, degree = degree,
intercept = TRUE, Boundary.knots = boundary_knots
),
check = my_check,
times = 1e3
)
```

```
: relative
Unit
expr min lq mean median uq max neval::bs 3.9034 3.6314 3.5353 3.5735 3.5633 1.1702 1000
splines::splineDesign 2.3132 2.1168 2.1689 2.0677 2.0819 1.3073 1000
splines::bSpline 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1000 splines2
```

Similarly, for derivatives of B-splines, `splines2::dbs()`

provides equivalent results with `splines::splineDesign()`

, and is about 2x faster.

```
## Derivatives of B-splines
derivs <- 2
microbenchmark(
"splines::splineDesign" = splineDesign(x, knots = all_knots,
ord = ord, derivs = derivs),
"splines2::dbs" = dbs(x, derivs = derivs, knots = internal_knots,
degree = degree, intercept = TRUE,
Boundary.knots = boundary_knots),
check = my_check,
times = 1e3
)
```

```
: relative
Unit
expr min lq mean median uq max neval::splineDesign 2.8409 2.6553 2.6413 2.5336 2.5257 1.0834 1000
splines::dbs 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1000 splines2
```

The **splines** package does not contain an implementation for integrals of B-splines. Thus, we performed a comparison with package **ibs** (version `r packageVersion("ibs")`

), where the function `ibs::ibs()`

was also implemented in **Rcpp**.

```
## integrals of B-splines
set.seed(123)
coef_sp <- rnorm(length(all_knots) - ord)
microbenchmark(
"ibs::ibs" = ibs::ibs(x, knots = all_knots, ord = ord, coef = coef_sp),
"splines2::ibs" = as.numeric(
splines2::ibs(x, knots = internal_knots, degree = degree,
intercept = TRUE, Boundary.knots = boundary_knots) %*%
coef_sp
),
check = my_check,
times = 1e3
)
```

```
: relative
Unit
expr min lq mean median uq max neval::ibs 22.49 20.654 21.372 21.103 21.016 25.442 1000
ibs::ibs 1.00 1.000 1.000 1.000 1.000 1.000 1000 splines2
```

The function `ibs::ibs()`

returns the integrated B-splines instead of the integrals of spline basis functions. Thus, we applied the same coefficients to the basis functions from `splines2::ibs()`

for equivalent results, which was still much faster than `ibs::ibs()`

.

For natural cubic splines (based on B-splines), `splines::ns()`

uses the QR decomposition to find the null space of the second derivatives of B-spline basis functions at boundary knots, while `splines2::naturalSpline()`

utilizes the closed-form null space derived from the second derivatives of cubic B-splines, which produces nonnegative basis functions (within boundary) and is more computationally efficient.

```
microbenchmark(
"splines::ns" = ns(x, knots = internal_knots, intercept = TRUE,
Boundary.knots = boundary_knots),
"splines2::naturalSpline" = naturalSpline(
x, knots = internal_knots, intercept = TRUE,
Boundary.knots = boundary_knots
),
times = 1e3
)
```

```
: relative
Unit
expr min lq mean median uq max neval::ns 5.3088 5.0697 4.9431 4.7628 4.7805 0.61586 1000
splines::naturalSpline 1.0000 1.0000 1.0000 1.0000 1.0000 1.00000 1000 splines2
```

The function `mSpline()`

produces periodic spline basis functions (based on M-splines) when `periodic = TRUE`

is specified. The `splines::periodicSpline()`

returns a periodic interpolation spline (based on B-splines) instead of basis matrix. Thus, we performed a comparison with package **pbs** (version `r packageVersion("pbs")`

), where the function `pbs::pbs()`

produces a basis matrix of periodic B-spline by using `splines::spline.des()`

(a wrapper function of `splines::splineDesign()`

).

```
microbenchmark(
"pbs::pbs" = pbs::pbs(x, knots = internal_knots, degree = degree,
intercept = TRUE, periodic = TRUE,
Boundary.knots = boundary_knots),
"splines2::mSpline" = mSpline(
x, knots = internal_knots, degree = degree, intercept = TRUE,
Boundary.knots = boundary_knots, periodic = TRUE
),
times = 1e3
)
```

```
: relative
Unit
expr min lq mean median uq max neval::pbs 3.6309 3.4413 3.5077 3.2826 3.2763 3.3523 1000
pbs::mSpline 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1000 splines2
```

## Session Information for Benchmarks

```
4.2.0 (2022-04-22)
R version : x86_64-pc-linux-gnu (64-bit)
Platform: Arch Linux
Running under
: default
Matrix products: /usr/lib/libopenblasp-r0.3.20.so
BLAS: /usr/lib/liblapack.so.3.10.1
LAPACK
:
locale1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
[4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
[
:
attached base packages1] splines stats graphics grDevices utils datasets methods base
[
:
other attached packages1] splines2_0.4.6 microbenchmark_1.4.9
[
namespace (and not attached):
loaded via a 1] Rcpp_1.0.9 codetools_0.2-18 ibs_1.4 digest_0.6.29 magrittr_2.0.3
[6] evaluate_0.16 rlang_1.0.4 stringi_1.7.8 cli_3.3.0.9000 rmarkdown_2.14
[11] tools_4.2.0 stringr_1.4.0 xfun_0.32 yaml_2.3.5 fastmap_1.1.0
[16] compiler_4.2.0 pbs_1.1 htmltools_0.5.3 knitr_1.39 [
```