I'm struggling with specifying the right R syntax for natural and (cubic) B-splines, using
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
"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.
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?