B-splines

A curve defined by f(x) is a spline of order k (degree k − 1), with knots ξ₁, ξ₂, …, ξ_m where ξ_j ≤ x_j + 1 and ξ_j < ξ_j + k∀j if f(x) is k − r − 1 times differentiable at any r-fold knot, and f(x) is a polynomial of order  ≤ k over each interval x ∈ [ξ_j, ξ_j + 1] for j = 0, …, m − 1.

In particular, B-splines, are defined as an affine combination:

where B_j, k, ξ(x) is the j^th basis spline function, ξ is a sequence of l interior knots and a total of 2k boundary knots, i.e., the cardinality of the knot sequence is |ξ| = 2k + l. The value of the k boundary knots is arbitrary, but a common choice is to use k-fold knots on the boundary:

ξ₁ = ξ₂ = ⋯ = ξ_k < ξ_k + 1 ≤ ⋯ ≤ ξ_k + l = ξ_{k + l + 1} = ⋯ = ξ_2k + l.

Alternative boundary knots can be used so long as the sequence ξ is non-decreasing. More on the implications of k-fold boundary knots follow in the next section.

The j^th B-spline is defined as:

where and

For a set of observations, x = x₁, …, x_n, the basis functions defined above generalize to a matrix:

Within the cpr package we can generate a basis matrix thusly:

x <- seq(0, 5.9999, length.out = 5000)
bmat <- bsplines(x, iknots = c(1, 1.5, 2.3, 4, 4.5), bknots = c(0, 6))
bmat
## Basis matrix dims: [5,000 x 9]
## Order: 4
## Number of internal knots: 5
## 
## First 6 rows:
## 
##           [,1]        [,2]         [,3]         [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 1.0000000 0.000000000 0.000000e+00 0.000000e+00    0    0    0    0    0
## [2,] 0.9964037 0.003593461 2.878634e-06 5.011451e-10    0    0    0    0    0
## [3,] 0.9928160 0.007172539 1.150485e-05 4.009161e-09    0    0    0    0    0
## [4,] 0.9892369 0.010737255 2.586411e-05 1.353092e-08    0    0    0    0    0
## [5,] 0.9856664 0.014287632 4.594188e-05 3.207329e-08    0    0    0    0    0
## [6,] 0.9821045 0.017823691 7.172363e-05 6.264314e-08    0    0    0    0    0

Note: the default order for bsplines is 4, and the default for the boundary knots is the range of x. However, relying on the default boundary knots can lead to unexpected behavior as, by definition the splines on the k-fold upper boundary is 0.

We can quickly view the plot of each of these spline functions as well.

plot(bmat, show_xi = TRUE, show_x = TRUE)

A few of many important properties of the basis functions:

B_j, k, xi(x) ∈ [0, 1], ∀x ∈ ℝ;

all(bmat >= 0)
## [1] TRUE
all(bmat <= 1)
## [1] TRUE

B_j, k, ξ(x) > 0  for  x ∈ (ξ_j, ξ_j + k);

and

$$\sum_j B_{j,k,\boldsymbol{xi}} \left(x\right) = \begin{cases} 1 & x \in \left[\xi_1, \xi_{2k+l}\right) \\ 0 & otherwise. \end{cases}$$

all.equal(rowSums(bmat), rep(1, nrow(bmat)))
## [1] TRUE

args(bsplines)
## function (x, iknots = NULL, df = NULL, bknots = range(x), order = 4L) 
## NULL
args(splines::bs)
## function (x, df = NULL, knots = NULL, degree = 3, intercept = FALSE, 
##     Boundary.knots = range(x), warn.outside = TRUE) 
## NULL

bs_mat <- splines::bs(x, knots = attr(bmat, "iknots"), Boundary.knots = attr(bmat, "bknots"))
str(attributes(bmat))
## List of 10
##  $ dim        : int [1:2] 5000 9
##  $ order      : num 4
##  $ df         : num 9
##  $ iknots     : num [1:5] 1 1.5 2.3 4 4.5
##  $ bknots     : num [1:2] 0 6
##  $ xi         : num [1:13] 0 0 0 0 1 1.5 2.3 4 4.5 6 ...
##  $ xi_star    : num [1:9] 0 0.333 0.833 1.6 2.6 ...
##  $ class      : chr [1:2] "cpr_bs" "matrix"
##  $ call       : language bsplines(x = x, iknots = c(1, 1.5, 2.3, 4, 4.5), bknots = c(0, 6))
##  $ environment:<environment: R_GlobalEnv>
str(attributes(bs_mat))
## List of 7
##  $ dim           : int [1:2] 5000 8
##  $ dimnames      :List of 2
##   ..$ : NULL
##   ..$ : chr [1:8] "1" "2" "3" "4" ...
##  $ degree        : int 3
##  $ knots         : num [1:5] 1 1.5 2.3 4 4.5
##  $ Boundary.knots: num [1:2] 0 6
##  $ intercept     : logi FALSE
##  $ class         : chr [1:3] "bs" "basis" "matrix"

bspline_eg <- bsplines(c(0, 1, 2, 5, 6), bknots = c(1, 5))
## Warning in bsplines(c(0, 1, 2, 5, 6), bknots = c(1, 5)): At least one x value <
## min(bknots)
## Warning in bsplines(c(0, 1, 2, 5, 6), bknots = c(1, 5)): At least one x value
## >= max(bknots)
bs_eg      <- splines::bs(c(0, 1, 2, 5, 6), Boundary.knots = c(1, 5), intercept = TRUE )
## Warning in splines::bs(c(0, 1, 2, 5, 6), Boundary.knots = c(1, 5), intercept =
## TRUE): some 'x' values beyond boundary knots may cause ill-conditioned bases

head(bspline_eg)
##          [,1]     [,2]     [,3]     [,4]
## [1,] 0.000000 0.000000 0.000000 0.000000
## [2,] 1.000000 0.000000 0.000000 0.000000
## [3,] 0.421875 0.421875 0.140625 0.015625
## [4,] 0.000000 0.000000 0.000000 0.000000
## [5,] 0.000000 0.000000 0.000000 0.000000
rowSums(bspline_eg)
## [1] 0 1 1 0 0

head(bs_eg)
##              1         2         3         4
## [1,]  1.953125 -1.171875  0.234375 -0.015625
## [2,]  1.000000  0.000000  0.000000  0.000000
## [3,]  0.421875  0.421875  0.140625  0.015625
## [4,]  0.000000  0.000000  0.000000  1.000000
## [5,] -0.015625  0.234375 -1.171875  1.953125
rowSums(bs_eg)
## [1] 1 1 1 1 1

Spline Spaces and Inserting a Knot

Consider two knot sequences ξ and ξ ∪ ξ′. Then, for a given polynomial order k, the spline space 𝕊_k, ξ ⊂ 𝕊_{k, ξ ∪ ξ′} (de Boor 2001, pg135). Given this relationship between spline spaces, and the convex sums generating spline functions, Boehm (1980) presented a method for inserting a knots into the knot sequence such that the spline function is unchanged. Specifically, for ξ′ = {ξ_j′|min (ξ) < ξ_j < max (ξ)  ∀j}, then there exist a θ_{ξ ∪ ξ′} such that B_k, ξ(x)θ_ξ = B_{k, ξ ∪ ξ′}(x)θ_{ξ ∪ ξ′}  ∀x ∈ [min (ξ), max (ξ)].

When inserting a singleton ξ′ into ξ, then θ_{ξ ∪ ξ′} = W_k, ξ(ξ′)θ_ξ where W_k, ξ(ξ′) is a (|θ| + 1) × |θ| lower bi-diagonal matrix $$ \boldsymbol{W}_{k,\boldsymbol{\xi}}\left(\xi'\right) = \begin{pmatrix} 1 & 0 & \cdots & 0 \\ 1 - \omega_{1, k, \boldsymbol{\xi}}\left(\xi'\right) & \omega_{1, k, \boldsymbol{\xi}} \left(\xi'\right) & \cdots & 0 \\ 0 & 1 - \omega_{2, k, \boldsymbol{\xi}}\left(\xi'\right) & \cdots 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 1 - \omega_{\left\lvert \boldsymbol{\theta} \right\rvert - 1, k, \boldsymbol{\xi}} \left(\xi'\right) & \omega_{\left\lvert \boldsymbol{\theta} \right\rvert - 1, k, \boldsymbol{\xi}} \left(\xi'\right) \\ 0 & 0 & 0 & 1 \end{pmatrix}, $$ with ω_j, k, ξ(x) as defined above in the de~Boor recursive algorithm. Through recursion we can insert as many knots into ξ without changing the value of the spline function.

For an example, insert a knot ξ′ = 3 into the control polygon defined above.

cp1 <- insert_a_knot(cp0, xi_prime = 3)
plot(cp0, cp1, color = TRUE, show_spline = TRUE)

Assessing influence of ξ_j

Here we derive a metric for assessing how much influence ξ_j ∈ ξ has on B_k, ξ(x)θ_ξ. Using the relationship defined by Boehm (1980), we can derive a metric for the influence of ξ_j on the spline function.

Start with an defined k, ξ, and θ_k, ξ. The relationship θ_ξ = W_{k, ξ ∖ ξj}(ξ_j)θ_{ξ ∖ ξj} holds if ξ_j has zero influence. However, in practice, we would expect that the relationship is θ_ξ = W_{k, ξ ∖ ξj}(ξ_j)θ_{ξ ∖ ξj} + d for a set of deviations d which would be equal to 0 if ξ_j has no influence on the spline.

We can estimate d via least squares. To simplify the notation in the following we drop some of the subscripts and parenthesis, that is, let W = W_{k, ξ ∖ ξj}(ξ_j).

$$ \begin{align} \boldsymbol{d} &= \boldsymbol{\theta}_{\boldsymbol{\xi}} - \boldsymbol{W} \boldsymbol{\theta}_{\boldsymbol{\xi}\backslash\xi_{j}} \\ &= \boldsymbol{\theta}_{\boldsymbol{\xi}} - \boldsymbol{W} \left(\boldsymbol{W}^{T}\boldsymbol{W}\right)^{-1}\boldsymbol{W}^{T} \boldsymbol{\theta}_{\boldsymbol{\xi}} \\ &= \left(\boldsymbol{I} - \boldsymbol{W} \left(\boldsymbol{W}^{T}\boldsymbol{W}\right)^{-1}\boldsymbol{W}^{T}\right) \boldsymbol{\theta}_{\boldsymbol{\xi}} \\ \end{align} $$

Finally, we define the influence of ξ_j on CP_{k, ξ, θξ} as w_j = ∥d∥₂².

The influence of knots on the spline used in the above section.

x <- influence_of_iknots(cp0)
summary(x)
##   j iknot  influence influence_rank chisq chisq_rank p_value os_p_value
## 1 5   1.0 1.64661178              5    NA         NA      NA         NA
## 2 6   1.5 0.29066719              2    NA         NA      NA         NA
## 3 7   2.3 0.31205029              3    NA         NA      NA         NA
## 4 8   4.0 0.07702981              1    NA         NA      NA         NA
## 5 9   4.5 0.41987740              4    NA         NA      NA         NA

Let’s look at the following plots to explore the influence of ξ₇ (the third interior knot in a k = 4 order spline) on the spline. In panel (a) we see the original control polygon and spline along with the coarsened control polygon and spline. Note that the there are fewer control points, and the spline deviates from the original. In panel (b) we see that the restored control polygon is the result of inserting ξ₇ into the coarsened control polygon of panel (a). These plots are also a good example of the local support and strong convexity of the control polygons as there are only k + 1 = 5 control points which are impacted by the removal and re-insertion of ξ₇. Lastly, in panel (c) we see all three control polygons plotted together.

ggpubr::ggarrange(
  ggpubr::ggarrange(
      plot(x, j = 3, coarsened = TRUE, restored = FALSE, color = TRUE, show_spline = TRUE) +
        ggplot2::theme(legend.position = "none")
    , plot(x, j = 3, coarsened = FALSE, restored = TRUE, color = TRUE, show_spline = TRUE) +
      ggplot2::theme(legend.position = "none")
    , labels = c("(a)", "(b)")
    , nrow = 1
  )
  , plot(x, j = 3, coarsened = TRUE, restored = TRUE, color = TRUE, show_spline = TRUE) +
    ggplot2::theme(legend.position = "bottom")
  , labels = c("", "(c)")
  , nrow = 2
  , ncol = 1
  , heights = c(1, 2)
)

Next, consider the influence of ξ₈, the fourth interior knot. By the influence metric defined above, this is the least influential knot in the sequence. This can be seen easily as the spline between the original and the coarsened spline are very similar, this is despite the apparent large difference in the magnitude of the control point ordinates between the original and coarsened control polygons. However, when re-inserting ξ₈ the recovered control polygon is very similar to the original, hence the low influence of ξ₈.

ggpubr::ggarrange(
  ggpubr::ggarrange(
      plot(x, j = 4, coarsened = TRUE, restored = FALSE, color = TRUE, show_spline = TRUE) +
        ggplot2::theme(legend.position = "none")
    , plot(x, j = 4, coarsened = FALSE, restored = TRUE, color = TRUE, show_spline = TRUE) +
      ggplot2::theme(legend.position = "none")
    , labels = c("(a)", "(b)")
    , nrow = 1
  )
  , plot(x, j = 4, coarsened = TRUE, restored = TRUE, color = TRUE, show_spline = TRUE) +
    ggplot2::theme(legend.position = "bottom")
  , labels = c("", "(c)")
  , nrow = 2
  , ncol = 1
  , heights = c(1, 2)
)

If you were required to omit an internal knot, it would be preferable to omit ξ₈ over ξ₅, ξ₆, ξ₇, or ξ₉ as that will have the least impact on the spline approximation of the original functional form.

Example with known knots

Apply CPR to the initial_cp from the above example.

cpr0 <- cpr(initial_cp)

cpr0
## A list of control polygons
## List of 7
##  - attr(*, "ioik")=List of 7
##  - attr(*, "class")= chr [1:2] "cpr_cpr" "list"

There are 7 control polygons within the cpr_cpr object, length(initial_cp$iknots) + 1. The indexing is set such at that the i^th element has i − 1 internal knots.

Before exploring the results, let’s just verify that the results of the call to cpr are the same as the manual results found about. There are some differences in the metadata of the objects, but the important parts, like the control polygons, are the same.

all.equal( cpr0[["cps"]][[7]][["cp"]],  initial_cp[["cp"]])
## [1] "target is NULL, current is data.frame"

# some attributes are different with the last cp due to how the automation
# creates the call vs how the call was created manually.
call_idx <- which(names(cpr0[["cps"]][[6]]) == "call")
all.equal( cpr0[["cps"]][[6]][-call_idx], first_reduction_cp [-call_idx])
## [1] "target is NULL, current is list"
all.equal( cpr0[["cps"]][[5]][-call_idx], second_reduction_cp[-call_idx])
## [1] "target is NULL, current is list"
all.equal( cpr0[["cps"]][[4]][-call_idx], third_reduction_cp [-call_idx])
## [1] "target is NULL, current is list"
all.equal( cpr0[["cps"]][[3]][-call_idx], fourth_reduction_cp[-call_idx])
## [1] "target is NULL, current is list"
all.equal( cpr0[["cps"]][[2]][-call_idx], fifth_reduction_cp [-call_idx])
## [1] "target is NULL, current is list"
all.equal( cpr0[["cps"]][[1]][-call_idx], sixth_reduction_cp [-call_idx], check.attributes = FALSE)
## [1] "target is NULL, current is list"

In the manual process we identified third_reduction_cp as the preferable model. For the cpr0 object we can quickly see a similar result as we did for the manual process.

s0 <- summary(cpr0)
knitr::kable(s0, row.names = TRUE)

The additional columns in this summary, loglik_elbow and rse_elbow, indicate a location in the plot for either the loglik or rse by model index (degrees of freedom) where the trade-off between additional degrees of freedom and improvement in the fix statistic is negligible. See plot below. This is determined by finding the breakpoint such that a continuous, but not differentiable at breakpoint, quadratic spline fits the plot with minimal residual standard error.

ggpubr::ggarrange(
  ggpubr::ggarrange(
      plot(cpr0, color = TRUE)
    , plot(cpr0, show_cp = FALSE, show_spline = TRUE, color = TRUE)
    , ncol = 1
    , common.legend = TRUE
  )
  ,
  ggpubr::ggarrange(
      plot(s0, type = "rse")
    , plot(s0, type = "rss")
    , plot(s0, type = "loglik")
    , plot(s0, type = "wiggle")
    , plot(s0, type = "fdsc")
    , plot(s0, type = "Pr(>w_(1))")
  , common.legend = TRUE
  )
  , widths = c(2, 3)
)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

Example when knots are unknown

In practice it is be extremely unlikely to know where knots should be placed. Analytic solutions are difficult, if not impossible to derive (Jupp 1978). However, an optimal solution may not be necessary.

From de Boor (2001) (page 106) “…a B-spline doesn’t change much if one changes its k + 1 knots a little bit. Therefore, if one has multiple knots, then it is very easy to find B-spline almost like with simple knots: Simply replace knot of multiplicity r > 1 by r simple knots nearby.”

That is,

B_k, ξ(x)θ_ξ ≈ B_k, ξ′(x)θ_ξ′ where |ξ_j − ξ_j′| < δ  for small  δ > 0.

So, in the case when we are looking for a good set of knots (a good set of knots should be parsimonious, and provide a good model fit) we start with an initial knot sequence with a high cardinality, this is not without precedent Eilers and Marx (2010). We then apply the CPR algorithm to find a good set of knots.

For example we will use 50 internal knots. Not surprisingly we have a fit that is more “connect-the-dots” than a smooth fit.

initial_cp <- cp(y ~ bsplines(x, df = 54, bknots = c(0, 6)), data = DF)

ggpubr::ggarrange(
  plot(initial_cp, show_cp = TRUE, show_spline = FALSE) + original_data_ggplot_layers
  ,
  plot(initial_cp, show_cp = FALSE, show_spline = TRUE) + original_data_ggplot_layers
)

Apply CPR to the initial_cp and look at the summary. Only the first 10 of 51 rows are provided here.

cpr1 <- cpr(initial_cp)

x <- summary(cpr1)
knitr::kable(head(x, 10))

From this, the preferable model is suggested to be index 5, the model with four internal knots. Inspection of the rse by index plot, I would argue from a manual selection that index 5 is preferable overall.

plot(x, type = "rse")

Let’s compare the models in indices 3, 4, and 5.

ggpubr::ggarrange(
  plot(cpr1[[3]], cpr1[[4]], cpr1[[5]], show_cp = TRUE, show_spline = FALSE, color = TRUE)
  ,
  plot(cpr1[[3]], cpr1[[4]], cpr1[[5]], show_cp = FALSE, show_spline = TRUE, color = TRUE)
  ,
  common.legend = TRUE
)

The noticeable differences are, for the most part, located on the left side of the support between between 0 and 1. For the fifth index, there is change in convexity of the spline. Knowing that the models at index 3 and 4 are good fits too, then it would be easy not select index 5 due to extra “wiggle” in the spline.

In practice I would likely pick the model in index 4 to have a smooth (small wiggle) and low rse. Compare the selected model to the original data.

ggpubr::ggarrange(
    plot(cpr1[[3]], show_cp = FALSE, show_spline = TRUE) +
    ggplot2::ggtitle("Model Index 3") +
    original_data_ggplot_layers
  , plot(cpr1[[4]], show_cp = FALSE, show_spline = TRUE) +
    ggplot2::ggtitle("Model Index 4") +
    original_data_ggplot_layers
  , plot(cpr1[[5]], show_cp = FALSE, show_spline = TRUE) +
    ggplot2::ggtitle("Model Index 5") +
    original_data_ggplot_layers
  , nrow = 1
  , legend = "none"
)