---
title: "Introduction to splines2"
author: Wenjie Wang
date: "`r Sys.Date()`"
bibliography:
- ../inst/bib/splines2.bib
vignette: >
  %\VignetteIndexEntry{Introduction to splines2}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
output: rmarkdown::html_vignette
---

```{r setup, echo = FALSE}
knitr::opts_knit$set(global.par = TRUE)
knitr::opts_chunk$set(fig.width = 7, fig.height = 4)
```

```{r set-par, echo = FALSE}
library(graphics)
par(mar = c(2.5, 2.5, 0.2, 0.1), mgp = c(1.5, 0.5, 0))
```

The **splines2** package (version `r packageVersion("splines2")`) provides
functions constructing a variety of regression spline basis functions 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.


## B-splines with their integrals and derivatives {#bSpline}

Function `bSpline()` provides B-spline basis matrix and allows `degree = 0` for
piece-wise constant basis function, which extends the `bs()` function in package
**splines** with a better computational performance.  One example of linear
B-splines with three internal knots is given as follows:

```{r bSpline, fig.cap="B-splines of degree one with three internal knots placed at 0.3, 0.5, and 0.6."}
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
@boor1978practical is implemented in function `ibs()` and `dbs()`, respectively.
Two toy examples are given as follows:

```{r ibs, fig.cap="Piecewise linear B-splines (left) and their integrals (right)."}
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")
```

```{r dbs, fig.cap="Cubic B-spline (left) and their first derivative (right)."}
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:

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

## M-splines using `mSpline()` {#mSpline}

M-splines [@ramsay1988monotone] can be considered as a normalized version of
B-splines with unit integral within boundary knots.  An example given by
@ramsay1988monotone 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.

```{r reset-par-mSpline, echo = FALSE}
par(mfrow = c(1, 1))
```

```{r mSpline, fig.cap = "Quadratic M-spline with three internal knots placed at 0.3, 0.5, and 0.6."}
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:

```{r mSpline-derivs}
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 splines based on M-spline basis
functions when `periodic = TRUE` is specified.  The `Boundary.knots` defines the
start and end points of the cyclic period.  The construction follows periodic
B-splines discussed in @piegl1997nurbs [chapter 12].

```{r pms-basis, fig.cap = "Cubic periodic M-splines."}
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 Basis")
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`.

```{r pms-deriv, fig.cap = "The first derivatives of the periodic M-splines."}
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-splines by specifying
`integral = TRUE`.  The integral is defined to be integrated from the left
boundary knot.

```{r pms-integral, fig.cap = "The integrals of the periodic M-splines."}
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()` {#iSpline}

I-splines [@ramsay1988monotone] 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 basis functions with non-negative (non-positive) coefficients *plus
a constant*, where the coefficient of the constant is unconstrained.

The example given by @ramsay1988monotone 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.

```{r iSpline, fig.cap = "I-splines of degree two with three internal knots placed at 0.3, 0.5, and 0.6."}
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.

```{r msMat}
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.

```{r dmsMat}
dmsMat3 <- deriv(isMat, 2)
stopifnot(is_equivalent(dmsMat1, dmsMat3))
```

## C-splines using `cSpline` {#cSpline}

Convex splines [@meyer2008inference] called C-splines are scaled integrals of
I-splines with unit maximum value at the right boundary knot.
@meyer2008inference 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 basis functions 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.


```{r cSpline-scaled, fig.cap = "C-splines of degree two with three internal knots placed at 0.3, 0.5, and 0.6."}
csMat1 <- cSpline(x, knots = knots, degree = 2, intercept = TRUE)
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.

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

## Generalized Bernstein Polynomials

The Bernstein polynomials have also been applied to shape-constrained regression
analysis [@wang2012csda].  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:

```{r bp-1, fig.cap = "Bernstein polynomials of degree 4 over [0, 1] (left) and the generalized version over [- 1, 1] (right)."}
x1 <- seq.int(0, 1, 0.01)
x2 <- seq.int(- 1, 1, 0.01)
bpMat1 <- bernsteinPoly(x1, degree = 4, intercept = TRUE)
bpMat2 <- bernsteinPoly(x2, degree = 4, intercept = TRUE)
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.

```{r bp-2, fig.height=6, fig.cap = "The integrals (upper panel) and the first derivatives (lower panel) of Bernstein polynomials of degree 4."}
ibpMat1 <- bernsteinPoly(x1, degree = 4, intercept = TRUE, integral = TRUE)
ibpMat2 <- bernsteinPoly(x2, degree = 4, intercept = TRUE, integral = TRUE)
dbpMat1 <- bernsteinPoly(x1, degree = 4, intercept = TRUE, derivs = 1)
dbpMat2 <- bernsteinPoly(x2, degree = 4, intercept = TRUE, derivs = 1)
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.

```{r bp-deriv}
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)))
```

## Natural Cubic Splines

The function `naturalSpline()` returns nonnegative basis functions (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 basis are
nonnegative within boundary.  When `integral = TRUE`, `naturalSpline()` returns
integral of each natural spline basis.

```{r ns-basis, fig.cap = "Nonnegative natural cubic splines (left) and corresponding integrals (right)."}
nsMat <- naturalSpline(x, knots = knots, intercept = TRUE)
insMat <- naturalSpline(x, knots = knots, intercept = TRUE, integral = TRUE)
par(mfrow = c(1, 2))
matplot(x, nsMat, type = "l", ylab = "Basis")
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 basis functions.

```{r ns-deriv, fig.cap = "The derivatives of natural cubic splines."}
d1nsMat <- naturalSpline(x, knots = knots, intercept = TRUE, derivs = 1)
d2nsMat <- deriv(nsMat, 2)
par(mfrow = c(1, 2))
matplot(x, d1nsMat, type = "l", ylab = "The 1st derivatives")
matplot(x, d2nsMat, type = "l", ylab = "The 2nd derivatives")
```


## Evaluation on New Values by `predict` {#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:

```{r predict}
new_x <- c(0.275, 0.525, 0.8)
names(new_x) <- paste0("x=", new_x)
predict(isMat, new_x)
```

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.



## Reference

