2. Piecewise Polynomials and Regresson Splines
This section presents a succession of increasingly complex segmented polynomial regression models, culminating in traditional cubic regression splines. We show how the models can be fit by manually coding the regressors in the model (i.e., the model matrix). Later on, we’ll explain how to fit these models using the gspline() (“generalized spline”) function, and, where appropriate, the bs() (“B-spline”) and ns() (“natural spline”) functions in the splines package, which is part of the standard R distribution. For simplicity, we work in this section with a single numeric predictor, \(x\), but everything that we do generalizes straightforwardly to numeric predictors in more complex regression models. The exposition here borrows heavily from, but also extends, Fox (2016, Sec. 17.2). A complementary vignette explains the theory underlying our implementation of regression splines.
Let’s start by randomly generating some nonlinear data:
The line on the graph represents the “true” population regression curve.
Piecewise polynomials begin by dividing the range of \(x\) into non-overlapping intervals at \(k\) ordered \(x\)-values \(t_j\) called knots: \((-\inf, t_1], (t_1, t_2], \ldots, (t_k, \inf)\). Then a degree-\(p\) polynomial is fit by least-squares regression to the data in each interval. The simplest case is an degree-0 polynomial, which generates a piecewise-constant fit by computing the mean \(y\)-value in each interval. For the example, we set \(k = 2\) knots at \(t_1 = 10/3\) and \(t_2 = 2\times10/3\), which are the 1/3 and 2/3 quantiles of the uniform[0, 10] distribution from which the \(x\)-values were generated. For real data, we can pick knots using sample \(x\)-quantiles or (as we discuss in another vignette) theoretically interesting transitional points in the data.
plot(x, y, main="Piecewise Constant Fit")
curve(cos(1.25*(x + 1)) + x/5, add=TRUE, lty=2, lwd=2)
t <- c(10/3, 20/3) # knots
abline(v=t, lty=3)
mtext(c(expression(t[1]), expression(t[2])), side=1, line=0.5, at=t)
group <- cut(x, breaks=c(-Inf, t, Inf))
means <- tapply(y, group, mean) # within-interval means
x0 <- c(0, t, 10) # knots augmented with boundaries
for (i in 1:3) lines(x0[c(i, i+1)], rep(means[i],2), lwd=2)
## [1] 0.4853018
RMSE is the root-mean-squared error of estimating \(E(y)\) for the cases in the data. It provides a rough guide to how well the model recovers the population regression curve, with the caveats that (1) it is an estimate because we have only \(n = 200\) sampled cases, and (2) we haven’t penalized the RMSE for the number of parameters in the model.
A next step is to generalize the model to a piecewise-linear fit (that is, a piecewise degree-1 polynomial):
mods2 <- by(data.frame(x=x, y=y, group=group), group,
function(df) lm(y ~ x, data=df)) # within-group linear regressions
plot(x, y, main="Piecewise Linear Fit")
curve(cos(1.25*(x + 1)) + x/5, add=TRUE, lty=2, lwd=2)
abline(v=t, lty=3)
mtext(c(expression(t[1]), expression(t[2])), side=1, line=0.5, at=t)
for (j in 1:3){
x1 <- x0[c(j, j + 1)]
y1 <- predict(mods2[[j]], list(x=x1))
lines(x1, y1, lwd=2) # draw within-group regression lines
}
## [1] 0.1336994
We convert the piecewise-linear fit into a linear regression spline by constraining the regression lines on the two sides of each knot to be equal at the knot. We can do this by fitting the linear model \[
y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + \varepsilon_i
\] where \(x_{i1} = x_i\), \[
x_{i2} = \left\{ {\begin{array}{ll}
0 & \mathrm{for} \: x_i \le t_1 \\
x_i - t_1 & \mathrm{for} \: x_i > t_1 \\
\end{array} } \right.
\] and \[
x_{i3} = \left\{ {\begin{array}{ll}
0 & \mathrm{for} \: x_i \le t_2 \\
x_i - t_2 & \mathrm{for} \: x_i > t_2 \\
\end{array} } \right.
\] Notice that the regression-spline coefficients each have a straightforward interpretation: \(\beta_0\) is the expectation of the response \(y\) when \(x = 0\); \(\beta_1\) is the slope of the regression line in the first interval; \(\beta_2\) is the change in slope between the first and second interval; and \(\beta_3\) is the change in slope between the second and third interval.
Fitting the linear regression spline to our artificial data produces the following result:
x1 <- x # define spline regressors
x2 <- (x - t[1]) * (x > t[1])
x3 <- (x - t[2]) * (x > t[2])
mod3 <- lm(y ~ x1 + x2 + x3) # linear regression spline model
plot(x, y, main="Linear Regression Spline")
curve(cos(1.25*(x + 1)) + x/5, add=TRUE, lty=2, lwd=2)
abline(v=t, lty=3)
mtext(c(expression(t[1]), expression(t[2])), side=1, line=0.5, at=t)
x01 <- x0
x02 <- (x0 - t[1]) * (x0 > t[1])
x03 <- (x0 - t[2]) * (x0 > t[2])
lines(x0, predict(mod3, list(x1=x01, x2=x02, x3=x03)), lwd=2) # draw fitted lines
## [1] 0.2292751
Because of the 2 constraints that the regression lines join at the 2 knots, the linear regression spline uses 2 fewer parameters than the unconstrained within-interval linear regressions: 4 parameters for the linear regression spline, \(3 \times 2 = 6\) parameters for the unconstrained linear regressions.
This approach generalizes readily to higher-degree polyomials, the most common of which is the third-degree or cubic regression spline. Consider first unconstrained within-interval cubic regressions:
mods4 <- by(data.frame(x=x, y=y, group=group), group, # within-interval cubic regressions
function(df) lm(y ~ poly(x, degree=3, raw=TRUE), data=df))
plot(x, y, main="Piecewise Cubic Fit")
curve(cos(1.25*(x + 1)) + x/5, add=TRUE, lty=2, lwd=2)
abline(v=t, lty=3)
mtext(c(expression(t[1]), expression(t[2])), side=1, line=0.5, at=t)
xx <- seq(0, 10, length=500)
for (j in 1:3){
xj <- xx[xx >= x0[j] & xx < x0[j + 1]]
yj <- predict(mods4[[j]], list(x=xj))
lines(xj, yj, lwd=2) # draw within-interval cubics
}
## [1] 0.0110144
This model uses \(4 \times 3 = 12\) parameters.
As before, we can constrain the cubic polynomials to join at the knots, generating the model \[\begin{align*}
y_i = \beta_0 &+ \beta_{11} x_{i1} + \beta_{12} x_{i1}^2 + \beta_{13} x_{i1}^3 \\
&+ \beta_{21} x_{i2} + \beta_{22} x_{i2}^2 + \beta_{23} x_{i2}^3 \\
&+ \beta_{31} x_{i3} + \beta_{32} x_{i3}^2 + \beta_{33} x_{i2}^3 + \varepsilon_i
\end{align*}\]
mod5 <- lm(y ~ poly(x1, degree=3, raw=TRUE) +
poly(x2, degree=3, raw=TRUE) +
poly(x3, degree=3, raw=TRUE))
plot(x, y, main="Cubic Regression Spline with Order-0 Smoothness")
curve(cos(1.25*(x + 1)) + x/5, add=TRUE, lty=2, lwd=2)
abline(v=t, lty=3)
mtext(c(expression(t[1]), expression(t[2])), side=1, line=0.5, at=t)
x01 <- xx
x02 <- (xx - t[1]) * (xx > t[1])
x03 <- (xx - t[2]) * (xx > t[2])
lines(xx, predict(mod5, list(x1=x01, x2=x02, x3=x03)), lwd=2)
## [1] 0.009149368
Because this model uses only 1 intercept, it has 2 fewer parameters than the unconstrained within-interval cubic regressions. The cubic spline is continuous but, as close inspection reveals, it isn’t smooth at the knots, where the slope of the fitted regression is free to change. We say that it has order-0 smoothness.
If we eliminate the linear terms \(\beta_{21} x_{i2}\) and \(\beta_{31} x_{i3}\) from the model we force the slope—that is, the first derivative—of the regression function to be equal on both sides of each knot, producing order-1 smoothness: \[\begin{align*}
y_i = \beta_0 + \beta_{11} x_{i1} &+ \beta_{12} x_{i1}^2 + \beta_{13} x_{i1}^3 \\
&+ \beta_{22} x_{i2}^2 + \beta_{23} x_{i2}^3 \\
&+ \beta_{32} x_{i3}^2 + \beta_{33} x_{i2}^3 + \varepsilon_i
\end{align*}\]
mod6 <- lm(y ~ poly(x1, degree=3, raw=TRUE) + I(x2^2) + I(x3^2) + I(x2^3) + I(x3^3))
plot(x, y, main="Cubic Regression Spline with Order-1 Smoothness")
curve(cos(1.25*(x + 1)) + x/5, add=TRUE, lty=2, lwd=2)
abline(v=t, lty=3)
mtext(c(expression(t[1]), expression(t[2])), side=1, line=0.5, at=t)
lines(xx, predict(mod6, list(x1=x01, x2=x02, x3=x03)), lwd=2)
## [1] 0.00902555
Finally, removing both the linear and the quadratic terms \(\beta_{22} x_{i2}^2\) and \(\beta_{32} x_{i3}^2\) forces the slope (first derivative) and curvature (second derivative) to be equal on both sides of each knot, producing order-2 smoothess and, except for parametrization (discussed below), a traditional cubic regression spline: \[
y_i = \beta_0 + \beta_{11} x_{i1} + \beta_{12} x_{i1}^2 + \beta_{13} x_{i1}^3 + \beta_{23} x_{i2}^3 + \beta_{33} x_{i3}^3 + \varepsilon_i
\] This model uses \(4 + 2 = 7\), or more generally \(4 + k\), parameters, considerably fewer than the \(4 \times 3 = 12\), or more generally \(4 \times (k + 1)\) parameters used by the unconstained piecewise cubic model.
mod7 <- lm(y ~ poly(x1, degree=3, raw=TRUE) + I(x2^3) + I(x3^3))
plot(x, y, main="Cubic Regression Spline with Order-2 Smoothness")
curve(cos(1.25*(x + 1)) + x/5, add=TRUE, lty=2, lwd=2)
abline(v=t, lty=3)
mtext(c(expression(t[1]), expression(t[2])), side=1, line=0.5, at=t)
lines(xx, predict(mod7, list(x1=x01, x2=x02, x3=x03)), lwd=2)
## [1] 0.02948439
5. How gspline() and bs() Construct Bases for Fitting Regression Splines
We’ve demonstrated how (1) directly coding spline regressors, (2) using gspline(), and (3) using bs() construct equivalent regression-spline bases that produce, within rounding error, identical fits to the data for cases, such as a standard cubic regressions spline, to which all three apply. In this section, we look more closely at the spline bases constructed directly and by gspline() and bs().
We return to the cubic regression spline with two knots at \(t_1\) and \(t_2\) discussed in Section 2: \[
y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + \varepsilon_i
\] where \(x_{i1} = x_i\), \[
x_{i2} = \left\{ {\begin{array}{ll}
0 & \mathrm{for} \: x_i \le t_1 \\
x_i - t_1 & \mathrm{for} \: x_i > t_1 \\
\end{array} } \right.
\] and \[
x_{i3} = \left\{ {\begin{array}{ll}
0 & \mathrm{for} \: x_i \le t_2 \\
x_i - t_2 & \mathrm{for} \: x_i > t_2 \\
\end{array} } \right.
\]
Here are the bases constructed directly, by gspline(), and by bs(), respectively named X, X.gsp, and X.bs, for a simple \(x\) vector with the integers from 0 to 10:
x <- 0:10
t <- c(10/3, 20/3)
x1 <- x # define spline regressors
x2 <- (x - t[1]) * (x > t[1])
x3 <- (x - t[2]) * (x > t[2])
X <- cbind(x1, x1^2, x1^3, x2^3, x3^3)
colnames(X) <- c("x1", "x1^2", "x1^3", "x2^3", "x3^3")
round(X, 5)
## x1 x1^2 x1^3 x2^3 x3^3
## [1,] 0 0 0 0.00000 0.00000
## [2,] 1 1 1 0.00000 0.00000
## [3,] 2 4 8 0.00000 0.00000
## [4,] 3 9 27 0.00000 0.00000
## [5,] 4 16 64 0.29630 0.00000
## [6,] 5 25 125 4.62963 0.00000
## [7,] 6 36 216 18.96296 0.00000
## [8,] 7 49 343 49.29630 0.03704
## [9,] 8 64 512 101.62963 2.37037
## [10,] 9 81 729 181.96296 12.70370
## [11,] 10 100 1000 296.29630 37.03704
## Spline function created by gspline
## $knots
## [1] 3.333333 6.666667
##
## $degree
## [1] 3 3 3
##
## $smoothness
## $smoothness[[1]]
## [1] 0 1 2
##
## $smoothness[[2]]
## [1] 0 1 2
##
##
## $G
## D1|0 D2|0 D3|0 C3|3.33 C3|6.67
## X0 0 0.0 0.00000 0.00000 0.00000
## X1 1 0.0 0.00000 0.00000 0.00000
## X2 0 0.5 0.00000 0.00000 0.00000
## X3 0 0.0 0.16667 0.00000 0.00000
## X0 0 0.0 0.00000 -6.17284 0.00000
## X1 1 0.0 0.00000 5.55556 0.00000
## X2 0 0.5 0.00000 -1.66667 0.00000
## X3 0 0.0 0.16667 0.16667 0.00000
## X0 0 0.0 0.00000 -6.17284 -49.38272
## X1 1 0.0 0.00000 5.55556 22.22222
## X2 0 0.5 0.00000 -1.66667 -3.33333
## X3 0 0.0 0.16667 0.16667 0.16667
##
## $constraint_mat
## X0 X1 X2 X3 X0 X1 X2 X3
## C0|3.33 -1 -3.33333 -11.11111 -37.03704 1 3.33333 11.11111 37.03704
## C1|3.33 0 -1.00000 -6.66667 -33.33333 0 1.00000 6.66667 33.33333
## C2|3.33 0 0.00000 -2.00000 -20.00000 0 0.00000 2.00000 20.00000
## C0|6.67 0 0.00000 0.00000 0.00000 -1 -6.66667 -44.44444 -296.29630
## C1|6.67 0 0.00000 0.00000 0.00000 0 -1.00000 -13.33333 -133.33333
## C2|6.67 0 0.00000 0.00000 0.00000 0 0.00000 -2.00000 -40.00000
## f(0) 1 0.00000 0.00000 0.00000 0 0.00000 0.00000 0.00000
## X0 X1 X2 X3
## C0|3.33 0 0.00000 0.00000 0.0000
## C1|3.33 0 0.00000 0.00000 0.0000
## C2|3.33 0 0.00000 0.00000 0.0000
## C0|6.67 1 6.66667 44.44444 296.2963
## C1|6.67 0 1.00000 13.33333 133.3333
## C2|6.67 0 0.00000 2.00000 40.0000
## f(0) 0 0.00000 0.00000 0.0000
##
## $estimate_mat
## X0 X1 X2 X3 X0 X1 X2 X3 X0 X1 X2 X3
## D1|0 0 1 0 0 0 0 0 0 0 0 0 0
## D2|0 0 0 2 0 0 0 0 0 0 0 0 0
## D3|0 0 0 0 6 0 0 0 0 0 0 0 0
## C3|3.33 0 0 0 -6 0 0 0 6 0 0 0 0
## C3|6.67 0 0 0 0 0 0 0 -6 0 0 0 6
## D1|0 D2|0 D3|0 C3|3.33 C3|6.67
## f(0) 0 0.0 0.00000 0.00000 0.00000
## f(1) 1 0.5 0.16667 0.00000 0.00000
## f(2) 2 2.0 1.33333 0.00000 0.00000
## f(3) 3 4.5 4.50000 0.00000 0.00000
## f(4) 4 8.0 10.66667 0.04938 0.00000
## f(5) 5 12.5 20.83333 0.77160 0.00000
## f(6) 6 18.0 36.00000 3.16049 0.00000
## f(7) 7 24.5 57.16667 8.21605 0.00617
## f(8) 8 32.0 85.33333 16.93827 0.39506
## f(9) 9 40.5 121.50000 30.32716 2.11728
## f(10) 10 50.0 166.66667 49.38272 6.17284
## 1 2 3 4 5
## [1,] 0.00000 0.00000 0.00000 0.00000 0.000
## [2,] 0.54225 0.11025 0.00450 0.00000 0.000
## [3,] 0.55800 0.34200 0.03600 0.00000 0.000
## [4,] 0.33075 0.54675 0.12150 0.00000 0.000
## [5,] 0.12800 0.58800 0.28200 0.00200 0.000
## [6,] 0.03125 0.46875 0.46875 0.03125 0.000
## [7,] 0.00200 0.28200 0.58800 0.12800 0.000
## [8,] 0.00000 0.12150 0.54675 0.33075 0.001
## [9,] 0.00000 0.03600 0.34200 0.55800 0.064
## [10,] 0.00000 0.00450 0.11025 0.54225 0.343
## [11,] 0.00000 0.00000 0.00000 0.00000 1.000
## attr(,"degree")
## [1] 3
## attr(,"knots")
## [1] 3.333333 6.666667
## attr(,"Boundary.knots")
## [1] 0 10
## attr(,"intercept")
## [1] FALSE
## attr(,"class")
## [1] "bs" "basis" "matrix"
We’ll discuss boundary knots in the next section when we take up natural splines. The other attributes associated with the matrix returned by bs() are self-explanatory; for further information see ?bs.
The coding in X is entirely straightforward, directly reflecting the definitions of the variables \(x_1\), \(x_2\), and \(x_3\) above. Comparison of X and X.gsp suggests that the columns of the latter are multiples of those of the former as we can easily verify:
## x1 x1^2 x1^3 x2^3 x3^3
## [1,] 0 NaN NaN NaN NaN
## [2,] 1 2 6 0 0
## [3,] 1 2 6 0 0
## [4,] 1 2 6 0 0
## [5,] 1 2 6 6 0
## [6,] 1 2 6 6 0
## [7,] 1 2 6 6 0
## [8,] 1 2 6 6 6
## [9,] 1 2 6 6 6
## [10,] 1 2 6 6 6
## [11,] 1 2 6 6 6
The NaNs are produced by dividing 0 by 0 and so are innocuous here. We can produce X.gsp from X by post-multiplying the latter by the diagonal weight matrix \(\mathrm{diag}(1, 1/2, 1/6, 1/6, 1/6)\):
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0.0 0.00000 0.00000 0.00000
## [2,] 1 0.5 0.16667 0.00000 0.00000
## [3,] 2 2.0 1.33333 0.00000 0.00000
## [4,] 3 4.5 4.50000 0.00000 0.00000
## [5,] 4 8.0 10.66667 0.04938 0.00000
## [6,] 5 12.5 20.83333 0.77160 0.00000
## [7,] 6 18.0 36.00000 3.16049 0.00000
## [8,] 7 24.5 57.16667 8.21605 0.00617
## [9,] 8 32.0 85.33333 16.93827 0.39506
## [10,] 9 40.5 121.50000 30.32716 2.11728
## [11,] 10 50.0 166.66667 49.38272 6.17284
We explain below why gspline() uses this scaling of the columns.
Demonstrating the equivalence of the bases produces by bs() and gspline() is a bit harder because the coefficients for the individual columns of the B-spline basis aren’t straightforwardly interpretable (which, by the way, is the essential advantage of using gspline()), but an indirect way to do this is to project the columns of X.bs onto the column space of X.gsp; if the two span the same subspace, then the projection should be equal to X.bs. We can easily perform the projection by least-squares regression, extracting the projected columns as the “fitted values”:
## [1] TRUE
To understand the coding used by gspline() let’s rewrite the simple cubic-spline model with 2 knots in terms of \(x\), omitting the subscript \(i\) for cases and taking the expectation \(\mu = E(y)\) of the response: \[
\mu = \left\{ {\begin{array}{ll}
\beta_0 + \beta_{11}x + \beta_{12}x^2 + \beta_{13}x^3 & \mathrm{for} \: x \le t_1 \\
\beta_0 + \beta_{11}x + \beta_{12}x^2 + \beta_{13}x^3 + \beta_{23}(x - t_1)^3 & \mathrm{for} \: t_1 < x \le t_2 \\
\beta_0 + \beta_{11}x + \beta_{12}x^2 + \beta_{13}x^3 + \beta_{23}(x - t_1)^3 + \beta_{33}(x - t_2)^3 & \mathrm{for} \: x > t_2 \\
\end{array} } \right.
\] Now differentiate \(\mu\) repeatedly with respect to \(x\) in each of the 3 intervals generated by the 2 knots: \[
\begin{array}{ll}
\dfrac{d\mu}{dx} &= \left\{ {\begin{array}{ll}
\beta_{11} + 2\beta_{12}x + 3\beta_{13}x^2 & \mathrm{for} \: x \le t_1 \\
\beta_{11} + 2\beta_{12}x + 3\beta_{13}x^2 + 3\beta_{23}(x - t_1)^2 & \mathrm{for} \: t_1 < x \le t_2 \\
\beta_{11} + 2\beta_{12}x + 3\beta_{13}x^2 + 3\beta_{23}(x - t_1)^2 + 3\beta_{33}(x - t_2)^2 & \mathrm{for} \: x > t_2 \\
\end{array} } \right.\\
\\
\dfrac{d^2\mu}{dx^2} &= \left\{ {\begin{array}{ll}
2\beta_{12} + 6\beta_{13}x & \mathrm{for} \: x \le t_1 \\
2\beta_{12} + 6\beta_{13}x + 6\beta_{23}(x - t_1) & \mathrm{for} \: t_1 < x \le t_2 \\
2\beta_{12} + 6\beta_{13}x + 6\beta_{23}(x - t_1) + 6\beta_{33}(x - t_2) & \mathrm{for} \: x > t_2 \\
\end{array} } \right.\\
\\
\dfrac{d^3\mu}{dx^3} &= \left\{ {\begin{array}{ll}
6\beta_{13} & \mathrm{for} \: x \le t_1 \\
6\beta_{13} + 6\beta_{23} & \mathrm{for} \: t_1 < x \le t_2 \\
6\beta_{13} + 6\beta_{23} + 6\beta_{33} & \mathrm{for} \: x > t_2 \\
\end{array} } \right.\\
\end{array}
\] The change in the third derivative at each knot is consequently 6 times the coefficient of the corresponding cubic regressor: \[
\begin{array}{ll}
\dfrac{d^3\mu}{dx^3}\biggr\rvert_{x = t_{1+}} - \dfrac{d^3\mu}{dx^3}\biggr\rvert_{x = t_{1}} &= 6\beta_{23}\\
\dfrac{d^3\mu}{dx^3}\biggr\rvert_{x = t_{2+}} - \dfrac{d^3\mu}{dx^3}\biggr\rvert_{x = t_{1}} &= 6\beta_{33}\\
\end{array}
\] where \(\dfrac{d^3\mu}{dx^3}\biggr\rvert_{x = t_{j+}}\) is the third derivative evaluated just to the right of knot \(t_j\).
Moreover, had the spline been a quadratic rather than a cubic regression spline, then the differences in the second derivative at the knots would have been 2 times the corresponding coefficients of the quadratic spline regressors, \(2\beta_{22}\) and \(2\beta_{23}\); and, finally, had the spline been a linear regression spline, then the differences in the first derivative at the knots would simply have been equal to the corresponding coefficients of the linear spline regressors, \(\beta_{12}\) and \(\beta_{13}\).
These observations explain the coding used by gspline(), which therefore produces coefficients that directly give differences in the highest-order derivative at the knots (e.g., of the third derivative for a cubic regression spline). This result is reflected in the names of the columns of the model matrix generated by gspline(). In the example, D1|0, D2|0, and D3|0 are respectively the first, second, and third derivatives at the left boundary (i.e., \(x = 0\)); C3|3.33 is the change in the third derivative (symbolized by “3”) at the first interior knot (\(t_1 = 10/3 \approx 3.33\)), and C3|6.67 is the change in the third derivative at the second interior knot (\(t_1 = 20/3 \approx 6.67\)).
5.1 Numerical Conditioning
Although the spline regression bases constructed directly, by gspline(), and by bs() are equivalent in the sense that they span the same subspace, these bases have different numerical properties. One way to see this is by computing the correlations among the regressors in each case:
## x1 x1^2 x1^3 x2^3 x3^3
## x1 1.000 0.963 0.909 0.813 0.645
## x1^2 0.963 1.000 0.986 0.929 0.782
## x1^3 0.909 0.986 1.000 0.977 0.862
## x2^3 0.813 0.929 0.977 1.000 0.940
## x3^3 0.645 0.782 0.862 0.940 1.000
## D1|0 D2|0 D3|0 C3|3.33 C3|6.67
## D1|0 1.000 0.963 0.909 0.813 0.645
## D2|0 0.963 1.000 0.986 0.929 0.782
## D3|0 0.909 0.986 1.000 0.977 0.862
## C3|3.33 0.813 0.929 0.977 1.000 0.940
## C3|6.67 0.645 0.782 0.862 0.940 1.000
## 1 2 3 4 5
## 1 1.000 0.313 -0.484 -0.456 -0.297
## 2 0.313 1.000 0.272 -0.484 -0.451
## 3 -0.484 0.272 1.000 0.313 -0.371
## 4 -0.456 -0.484 0.313 1.000 0.027
## 5 -0.297 -0.451 -0.371 0.027 1.000
The correlations for the first two bases are the same because, as we just explained, the columns of X.gsp are multiples of those of X, and some are quite large, while the correlations among the columns of X.bs are much smaller.
Even more directly, checking the condition numbers of the three basis matrices, each augmented with a column of 1s for the regression intercept, shows that both X and X.gsp are much more ill-conditioned than X.bs:
## [1] 4736.711
## [1] 1554.19
## [1] 10.42067
Notice that the condition numbers of the augmented X and X.gsp matrices are not identical because of the different scalings of the columns.
These observations suggest that if we’re interested in regression splines simply for fitting smooth nonlinear regressions, and if the model is expressible as a B-spline, then there’s an advantage in using the B-spline basis in preference to the basis constructed by gspline(). As we mentioned, the advantages of gspline() are in supporting a wider variety of piecewise polynomial and regression spline models and, crucially, in producing interpretable coefficients and hypothesis tests. We illustrate uses of the latter in another vignette.
Finally, modern least-squares algorithms, based for example on the QR and singular-value decompositions, decrease the numerical advantages of having relatively uncorrelated regressors. These methods effectively orthogonalize the model matrix in the process of the computating least-squares fits.
Natural Splines
Natural splines are B-splines with additional boundary knots, typically placed at the extremes of \(x\), and the additional constraints that the fitted regression is linear to the left of the lower boundary knot and to the right of the upper boundary knots. Derivatives up to order 1 less than the degree of the interior spline polynomials also match at the boundaries. Natural splines in effect place 2 additional constraints on the model and therefore are representable using 2 fewer parameters than corresponding B-splines: A cubic natural regression spline with \(k\) knots uses \(1 + k\) parameters (and regressors), excluding the intecept, 2 fewer than the corresponding B-spline.
Because of the condition that the regression is linear beyond the boundaries, natural splines tend to produce less wild extrapolations than B-splines. Natural splines may also behave more reasonably in the interior near the boundaries, but are less flexible than B-splines. We can compensate for the lower flexibility of natural splines within the range of the data by introducing 1 or 2 additional interior knots, perhaps spending the degrees of freedom more wisely than for a B-spline with the same number of parameters. The ns() function in the standard R splines package only fits natural cubic regression splines.
To illustrate, we continue the example from the preceding section, with \(x = (0, 1, 2, \ldots, 10)\), \(k = 2\) interior knots at 10/3 and 20/3, and boundary knots at 0 and 10:
## 1 2 3
## [1,] 0.00000 0.00000 0.00000
## [2,] -0.06960 0.22231 -0.14821
## [3,] -0.10044 0.40933 -0.27289
## [4,] -0.05376 0.52578 -0.35052
## [5,] 0.10238 0.54085 -0.35923
## [6,] 0.32047 0.47608 -0.29655
## [7,] 0.50176 0.38672 -0.17248
## [8,] 0.54762 0.32813 0.00275
## [9,] 0.41887 0.32738 0.21774
## [10,] 0.16843 0.36771 0.45936
## [11,] -0.14286 0.42857 0.71429
## attr(,"degree")
## [1] 3
## attr(,"knots")
## [1] 3.333333 6.666667
## attr(,"Boundary.knots")
## [1] 0 10
## attr(,"intercept")
## [1] FALSE
## attr(,"class")
## [1] "ns" "basis" "matrix"
The Boundary.knots argument to ns() defaults to the range of x, and so we didn’t have to specify it explicitly but have done so for clarity. The knots argument gives the interior knots.
We can specify an equivalent basis using gspline() by specifying linear (i.e., degree-1) polynomials beyond the boundaries along with interior cubics, and by matching derivatives up to order-2:
sp_ns <- gspline(knots=c(0, 10/3, 20/3, 10),
degree=c(1, 3, 3, 3, 1),
smoothness=c(2, 2, 2, 2))
sp_ns
## Spline function created by gspline
## $knots
## [1] 0.000000 3.333333 6.666667 10.000000
##
## $degree
## [1] 1 3 3 3 1
##
## $smoothness
## $smoothness[[1]]
## [1] 0 1 2
##
## $smoothness[[2]]
## [1] 0 1 2
##
## $smoothness[[3]]
## [1] 0 1 2
##
## $smoothness[[4]]
## [1] 0 1 2
##
##
## $G
## D1|0 C3|0 C3|3.33
## X0 0 0.00000 0.00000
## X1 1 0.00000 0.00000
## X2 0 0.00000 0.00000
## X3 0 0.00000 0.00000
## X0 0 0.00000 0.00000
## X1 1 0.00000 0.00000
## X2 0 0.00000 0.00000
## X3 0 0.16667 0.00000
## X0 0 0.00000 -6.17284
## X1 1 0.00000 5.55556
## X2 0 0.00000 -1.66667
## X3 0 0.16667 0.16667
## X0 0 148.14815 92.59259
## X1 1 -66.66667 -38.88889
## X2 0 10.00000 5.00000
## X3 0 -0.33333 -0.16667
## X0 0 -185.18519 -74.07407
## X1 1 33.33333 11.11111
## X2 0 0.00000 0.00000
## X3 0 0.00000 0.00000
##
## $constraint_mat
## X0 X1 X2 X3 X0 X1 X2 X3 X0 X1 X2
## C0|0 -1 0 0 0 1 0.0000 0.0000 0.0000 0 0.0000 0.0000
## C1|0 0 -1 0 0 0 1.0000 0.0000 0.0000 0 0.0000 0.0000
## C2|0 0 0 -2 0 0 0.0000 2.0000 0.0000 0 0.0000 0.0000
## C0|3.33 0 0 0 0 -1 -3.3333 -11.1111 -37.0370 1 3.3333 11.1111
## C1|3.33 0 0 0 0 0 -1.0000 -6.6667 -33.3333 0 1.0000 6.6667
## C2|3.33 0 0 0 0 0 0.0000 -2.0000 -20.0000 0 0.0000 2.0000
## C0|6.67 0 0 0 0 0 0.0000 0.0000 0.0000 -1 -6.6667 -44.4444
## C1|6.67 0 0 0 0 0 0.0000 0.0000 0.0000 0 -1.0000 -13.3333
## C2|6.67 0 0 0 0 0 0.0000 0.0000 0.0000 0 0.0000 -2.0000
## C0|10 0 0 0 0 0 0.0000 0.0000 0.0000 0 0.0000 0.0000
## C1|10 0 0 0 0 0 0.0000 0.0000 0.0000 0 0.0000 0.0000
## C2|10 0 0 0 0 0 0.0000 0.0000 0.0000 0 0.0000 0.0000
## I.1.2 0 0 1 0 0 0.0000 0.0000 0.0000 0 0.0000 0.0000
## I.1.3 0 0 0 1 0 0.0000 0.0000 0.0000 0 0.0000 0.0000
## I.5.2 0 0 0 0 0 0.0000 0.0000 0.0000 0 0.0000 0.0000
## I.5.3 0 0 0 0 0 0.0000 0.0000 0.0000 0 0.0000 0.0000
## f(0) 1 0 0 0 0 0.0000 0.0000 0.0000 0 0.0000 0.0000
## X3 X0 X1 X2 X3 X0 X1 X2 X3
## C0|0 0.0000 0 0.0000 0.0000 0.0000 0 0 0 0
## C1|0 0.0000 0 0.0000 0.0000 0.0000 0 0 0 0
## C2|0 0.0000 0 0.0000 0.0000 0.0000 0 0 0 0
## C0|3.33 37.0370 0 0.0000 0.0000 0.0000 0 0 0 0
## C1|3.33 33.3333 0 0.0000 0.0000 0.0000 0 0 0 0
## C2|3.33 20.0000 0 0.0000 0.0000 0.0000 0 0 0 0
## C0|6.67 -296.2963 1 6.6667 44.4444 296.2963 0 0 0 0
## C1|6.67 -133.3333 0 1.0000 13.3333 133.3333 0 0 0 0
## C2|6.67 -40.0000 0 0.0000 2.0000 40.0000 0 0 0 0
## C0|10 0.0000 -1 -10.0000 -100.0000 -1000.0000 1 10 100 1000
## C1|10 0.0000 0 -1.0000 -20.0000 -300.0000 0 1 20 300
## C2|10 0.0000 0 0.0000 -2.0000 -60.0000 0 0 2 60
## I.1.2 0.0000 0 0.0000 0.0000 0.0000 0 0 0 0
## I.1.3 0.0000 0 0.0000 0.0000 0.0000 0 0 0 0
## I.5.2 0.0000 0 0.0000 0.0000 0.0000 0 0 1 0
## I.5.3 0.0000 0 0.0000 0.0000 0.0000 0 0 0 1
## f(0) 0.0000 0 0.0000 0.0000 0.0000 0 0 0 0
##
## $estimate_mat
## X0 X1 X2 X3 X0 X1 X2 X3 X0 X1 X2 X3 X0 X1 X2 X3 X0 X1 X2 X3
## D1|0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## C3|0 0 0 0 -6 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 0
## C3|3.33 0 0 0 0 0 0 0 -6 0 0 0 6 0 0 0 0 0 0 0 0
## D1|0 C3|0 C3|3.33
## f(0) 0 0.00000 0.00000
## f(1) 1 0.16667 0.00000
## f(2) 2 1.33333 0.00000
## f(3) 3 4.50000 0.00000
## f(4) 4 10.66667 0.04938
## f(5) 5 20.83333 0.77160
## f(6) 6 36.00000 3.16049
## f(7) 7 57.14815 8.20370
## f(8) 8 84.14815 16.14815
## f(9) 9 115.14815 26.09259
## f(10) 10 148.14815 37.03704
## [1] TRUE
Notice that gspline() doesn’t distinguish between interior and boundary knots, and that the distinction is imposed implicitly by the distribution of the \(x\)-values and the specification of the degree and smoothness arguments. Because all of the values of smoothness are equal to 2, we could have specified smoothness=2, but we entered the argument as a 4-element vector to emphasize the number of constraints placed on the derivatives.
