natural splines in R with ns( )

Question

I'm struggling with specifying the right R syntax for natural and (cubic) B-splines, using ns() and bs() of package "splines".

To keep things simple, suppose I have a linear regression model, dependent variable Y and independent X. The X runs from 1 to 100, each value appearing once, the number of cases N being 100. I would like to use a natural cubic spline to model X's influence on Y. For the lowest region 1-20 of X, and also for the highest region 80-100, I want the regression effect of X to be linear; for regions 20-40, 40-60, 60-80, I want the effect of X to be cubic and different for each of the three regions. In R, I used this syntax to estimate the model:

modelns1 <- lm(y ~ ns(x, knots=c(40, 60), 
                   Boundary.knots=c(20, 80)))

Looking at the predicted values of this model, it seems to do what I want: linear in the two extreme regions and curved in the middle parts. But I'm not sure if I'm doing this right. Specifying

modelns2 <- lm(y ~ ns(x, knots=c(20, 40, 60, 80)))

also leads to predictions which are very close to linear in the two outer regions.

I think my confusion is about the meaning of "boundary knots". The description in the help-page about ns() says:

"boundary points at which to impose the natural boundary conditions and anchor the B-spline basis (default the range of the data)"

I reasoned that the natural boundary conditions refer to setting the constant and the first derivative equal for the first two regions, 1-20 and 20-40, and for the last two regions 60-80 and 80-100. Hence, I thought I should specify 20 and 80 as "Boundary.knots" in the ns() function. I don't know if this is right, though.

By default, ns() chooses the min and max of x, here 1 and 100, as boundary knots, the help-page says. So with the second modelns2, the natural "lower" boundary condition would refer to the regions x < 1 and 1-20. But this doesn't seem to make sense, since how could the conditions be applied to non-existing data x < 1?

Conclusion: I don't get it! Also, help for specifying a B-spline would be appreciated greatly. For such spline, the effect of X in the two outer regions would be cubic as well. I think I have to use:

modelbs <- lm(y ~ bs(x, knots=c(20, 40, 60, 80)))

The graph of predictions against X indeed shows a nonlinear trend for the two outer regions. But is this indeed the syntax to use?

Answer 1

After some studying and trial and error, I found answers to the questions. Also I found how one can obtain a "common" cubic spline at one boundary, and a natural spline at the other boundary, which can be useful if one has many observations at one end the data and only a few at the other end.

The following lines generate some data with a nonlinear relation between x and y.

# Generate x and y.
x <- 1:100
set.seed <- 123456
y <- sin(20*x*pi/180) + 0.01*x + 0.001*(x-50)^2 + 0.00002*(x-50)^3
# Show the graph.
library(ggplot2)
ggplot()+aes(x,y)+geom_line()

Next, estimate a natural spline with linearity in the x intervals 1-20 and 80-100, but cubic in intervals 20-40, 40-60 and 60-80:

library(splines)
mymodel1 <- lm(y ~ ns(x, knots=c(40,60), Boundary.knots=c(20,80)))
summary(mymodel1)
predm <- predict(mymodel1)
ggplot()+aes(x, predm)+geom_line()+geom_vline(xintercept=c(20,40,60,80))

The above model has 4 df, which equals the total number of knots, which is 4, including the two boundary knots specified in ns().

We can instead specify another natural spline model with the boundary knots lying at the minimum and maximum values of x, being 1 and 100. In that case, the linearity assumption of the natural spline pertains to the x values below 1 and above 100 (which are not observed, or may even not be possible):

# Natural spline with linearity below 1 and above 100.
mymodel2 <- lm(y ~ ns(x, knots=c(20,40,60,80), Boundary.knots=c(1,100)))
summary(mymodel2)
predm <- predict(mymodel2)
ggplot()+aes(x, predm)+geom_line()+geom_vline(xintercept=c(20,40,60,80))

In contrast to mymodel1, the second mymodel2 is not linear in the outer intervals 1-20 and 80-100, but cubic. However, the cubic polynomials in the intervals 1-20 and 80-100 have a second derivative at x=1 and x=100, respectively, which is fixed to be zero. This is, because the second derivative at x=1 has to be equal to the second derivatives of the linear function assumed to exist for x lower than 1, and this derivative is zero, obviously; also, the second derivative at x=100 is zero, because for the linear function above x=100, the second derivative is zero as well. For mymodel2 the number of df equals 6, again equal to the total number of knots specified in ns(), the boundary knots included.

Next, two models m1 and m2 are shown, m1 estimated with ns() and m2 estimated with bs(), both producing identical r-squares, df's and predictions.

A trick is applied in model m1 to have ns() fit a "common" cubic spline function between the lowest two knots 1 and 20, just like bs() would do. Note that the left boundary knot is specified as 0, and not as 1! To the left of this fake boundary there are no data, but ns() assumes a linear effect of x in that region. Between x=0 and the first knot x=1, where again no data exist, a cubic is assumed by ns() with zero second derivative at knot x=1. Further, in the interval 1-20, a "common" cubic is fitted, like bs() would do. A similar trick is applied at the right end of the data, with the fake boundary knot specified at 101. As a result, ns() in model m1 gives the same results in terms of df, r-square and predicted values as bs() does in model m2, which comes next.

The linearity constraints for the two extremest regions and/or the second derivatives constraints for the two next-to-extremest regions of x, which are typical for natural splines are circumvented here by specifying fake boundary knots in ns(). What remains is a "common" cubic spline function, with df equal to the number of internal knots (20,40,60,80) plus 4 = 4 + 4 = 8. As before, one could also say: the df equals the total number of knots given in ns(), the specified boundary knots included.

Btw, the df is the total number parameters in both models, including the intercept.

 m1 <- lm(y ~ ns(x, knots=c(1,20,40,60,80,100), Boundary.knots=c(0,101)))
 summary(m1)
 predm <- predict(m1)
 ggplot()+aes(x, predm)+geom_line()+geom_vline(xintercept=c(20,40,60,80))
    
 m2 <- lm(y ~ bs(x, knots=c(20,40,60,80)))
 summary(m2)
 predm <- predict(m2)
 ggplot()+aes(x, predm)+geom_line()+geom_vline(xintercept=c(20,40,60,80))

With the same trick as above, in the next model a common cubic spline is fitted at the left end, and a natural spline is fitted at the right of the x data, with a right boundary knot at x=100. There are 7 df, equal to the number of knots mentioned in ns(), again including the boundary knots specified.

m3 <- lm(y ~ ns(x, knots=c(1,20,40,60,80), Boundary.knots=c(0,100)))
summary(m3)
predm <- predict(m3)
ggplot()+aes(x, predm)+geom_line()+geom_vline(xintercept=c(20,40,60,80))