12
$\begingroup$

I have been reading the very helpful introduction on splines on http://freakonometrics.hypotheses.org/9184 and on http://www.stats.uwo.ca/faculty/braun/ss3859/chapters/splines/splines.pdf, as well as the examples so helpfully given by whuber on this forum. I understand that the function h(x) can be approximated by $$\sum\limits_{k = 0}^d \alpha_k(x-\alpha)^k+ \sum\limits_{i = 1}^j \beta_i(x-x_i)^d_+$$

I have trouble wrapping my head around on how this formula links to the actual basis functions that are generated in e.g. the splines packages in R. For example - following the R code below, the code creates 4 columns that i believe are the basis functions (each composed of 100 points), that look as per the graph below and/or the (head of) the matrix below. This might be simple - but how do these values link back to the formula per above?

To further clarify - in spline regression, in my understanding $x$ is still the vector of observations of our (1 dimensional) covariate, and that we add truncated power terms that are only positive for those observations that are greater than these knots. I would hence have expected a set of hinge functions that are 0 until the knot, and then continue with the slope as implied by the derivative $\beta$. however - I fail to understand why the basis functions are "triangles", and rise up with slope $\beta$ to the knot and then drop downwards with negative slope $-\beta$ after the knot. I presume this is fairly elementary - but I am just failing to wrap my head around this.

set.seed(1)
n=10
xr = seq(0,n,by=.1)
yr = sin(xr/2)+rnorm(length(xr))/2
db = data.frame(x=xr,y=yr)
plot(db)

attach(db)
library(splines)
B=bs(xr,knots=c(2,5,8),Boundary.knots=c(0,10),degre=1)
B

matplot(xr,B,type="l")

reg=lm(yr~B)
lines(xr,predict(reg),col="red")

> B
#             1          2          3    4
# [1,] 0.00000000 0.00000000 0.00000000 0.00
# [2,] 0.05000000 0.00000000 0.00000000 0.00
# [3,] 0.10000000 0.00000000 0.00000000 0.00
# [4,] 0.15000000 0.00000000 0.00000000 0.00
# [5,] 0.20000000 0.00000000 0.00000000 0.00
# [6,] 0.25000000 0.00000000 0.00000000 0.00
# [7,] 0.30000000 0.00000000 0.00000000 0.00
# [8,] 0.35000000 0.00000000 0.00000000 0.00
# [9,] 0.40000000 0.00000000 0.00000000 0.00

enter image description here

$\endgroup$
3
  • 2
    $\begingroup$ They don't. You have not reproduced the formula correctly. The final "+" is meaningless as you have typeset it: it belongs as a subscript and (presumably, because your source does not define it) indicates the positive part of a function. For some worked examples and illustrations, check out stats.stackexchange.com/questions/23137 and possibly (if you really like reading code) stats.stackexchange.com/questions/97845. $\endgroup$
    – whuber
    Jul 27, 2015 at 15:43
  • $\begingroup$ Sorry. i will adjust the notation and question . i am still trying to wrap my head around the concept of splines, and even after checking yourexamples and also the useful introduction on stats.uwo.ca/faculty/braun/ss3859/chapters/splines/splines.pdf i am still somewhat confused on probably some very basic concepts. i updated the question to show where i am failing to connect the dots. $\endgroup$
    – Wouter
    Jul 27, 2015 at 18:20
  • $\begingroup$ Thank you for improving your question (+1). I agree that although this may be basic, it is not necessarily easy! In particular, I am always challenged by such plots as the one you made because they are not symmetric (when the axis is reversed), but intuitively one would expect that kind of symmetry. $\endgroup$
    – whuber
    Jul 27, 2015 at 18:32

2 Answers 2

7
$\begingroup$

You can't get there from here. The basis splines in your graph do not emerge as a straightforward algebraic manipulation of the equation you have supplied -- at least, not straightforward to me. But that's not where they come from.

The basis functions come from various theoretical results about B-splines. The spline function is the smoothest function that passes close to (or that interpolates) the sampled function values (the knot points). It can be shown that the solution to this optimization lies in a finite dimentional function space composed of piecewise polynomials -- the degree of which depends on how much smoothness you want. The kink in the polynomials happens at the knot points.

So now we go in search of sensible basis functions. In your case, you have piecewise linear functions .. so the simplest piecewise linear function is a tent function ... and the tent pole has to occur at a knot point, because that's where we get the break in differentiability. The basis functions supplied by R are not the only choice, but they produce nice band matrices, cheap and easy to invert and otherwise manipulate.

Note that your basis functions must also respect the boundary conditions of the problem you have set yourself. The basis functions above will only give functions with $s(0)=0$. If I add the constant function to my basis, I can interpolate functions with $s(0)=c$.

Now convince yourself that the tent functions shown above are a basis for the requisite spline space. Consider what happens when you add linear combinations of the functions you illustrated above: they will all have linear portions, with kinks at the knot points. None of these can be obtained from the others. Finally, you need to show that you have the right number of them (I can't remember the formula, off hand, for the dimension of the spline space in terms of the number of knots and the degree of the polys).

Smoother results would be obtained by increasing the degree of the polynomials -- you could have piecewise quadratics, or cubics (the usual choice) $\ldots$ and then your basis functions will look like a sequence of bells centered about the knot points.

The truncated polynomials from your equation can also be used to build the spline smoother or interpolant, but they do not have the attractive numeric properties of the tent functions ... so that's why R does not supply them.

I would not attempt to learn about spline functions from the references cited above. Ramsay and Hooker's Functional Data Analysis with R and Matlab ties in the theory with implementations in R. You could also dig up the original papers by Kimeldorf and Wahba on smoothing and interpolating splines.

$\endgroup$
6
$\begingroup$

I know this is an old question but maybe my answer can be useful to someone in the future.

It exists a connection between truncated polynomial functions bases (TPF) and B-splines. This is discussed, for example, in the following papers:

  1. Paul H.C. Eilers and Brian D. Marx - Splines, Knots, and Penalties, WIREs Computational Statistics (2010).
  2. Harald Binder, Willi Sauerbrei - Increasing the usefulness of additive spline models by knot removal, CSDA (2008)

My answer will follow ref.1 (further details can be found there). In order to keep notation simple, I will use B-splines defined on equally spaced knots. I will indicate with $x$ a vector of abscissa points and with $t_{1} < t_{2} <...$ a series of equally spaced knots. Then, the TPF-bases of order $p$ are obtained by computing: $$ f^{p}_{ij}(x_{i}, t_{j}) = (x_{i} - t_{j})^{p} I (x_{i} > t_{j}) $$ TPFs of order 1 look like this

enter image description here

How can we "transform" these broken lines in "triangles"? The operator that makes this possible is the difference operator: B-splines can be computed as (appropriately scaled) differences of adjacent TPFs. A B-spline basis of degree 1 (one of the triangles in the original question) can be obtained as: $$ B^{1}_{j}(x) = \frac{f_{j-2}^{1}(x) - 2 f^{1}_{j-1}(x) + f^{1}_{j}(x)}{h} $$

where $h = t_{2} - t_{1} = \mbox{const}$ is the distance between 2 consecutive knots (which is constant here because I am assuming equally spaced knots). By applying this operation to the 3 TPFs above, we obtain:

enter image description here

The formula can be, of course, generalized to any degree $p$: $$ B^{p}_{j}(x) = (-1)^{p+1} \frac{\Delta^{p+1} f_{j}^{p}(x)}{h^{p} p!} $$ where $\Delta^{p+1}$ is the $p+1$th order difference operator. Once again, this holds only for equally spaced knots. If the distance between successive knots is not constat, more complex transformations are available.

The code to reproduce my examples is rather simple (see below). Please note that, differently from the original quesiton, I will use the R function splineDesign() to create the B-splines. The bs() funciton is actually based on splineDesign() (see help page of bs()).

# Libraries
library(colorout)
library(scales)

# Get some data
set.seed(1)
n = 100
x = seq(0, 1, len = n)

# B-splines stuffs
xr    = max(x)
xl    = min(x)
ndx   = 4
deg   = 1
dx    = (xr - xl) / ndx
knots = seq(xl - deg * dx, xr + deg * dx, by = dx)
B     = splineDesign(knots = knots, x = x, ord = deg + 1, 
                derivs = 0, outer.ok = TRUE)

# Compute tpf basis and B-slines from them
tpower = function(x, t, p)  (x - t) ^ p * (x > t)
P = outer(x, knots, tpower, deg)

# Plot TPFs 
matplot(x, P[, 3:5], pch = 16, col = alpha(1:3, 0.5), 
         ylab = "")
matlines(x, P[, 3:5], pch = 16, lty = 1, 
                col = alpha(1:3, 0.5))

# Compute B-splines as differences of TPFs
Bp = t(apply(P, 1, diff, diff = 2)) / dx

# Plot single B-spline and compare with splineDesign
plot(x, Bp[, 3], pch = 16, ylab=" ")
lines(x, B[, 3], pch = 16, lty = 1)

# Compare all with splineDesign B-splines
matplot(Bp,col = alpha(1:ncol(Bp), 0.5), pch = 16)
matlines(B,col = alpha(1:ncol(Bp), 0.5), lty = 1)

Hope this is somehow helpful and it replies (at least partially) the original question.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.