Natural Cubic Spline Basis for Polynomial SplinesSource:
nsk() generate the natural cubic
spline basis functions, the corresponding derivatives or integrals (from the
left boundary knot). Both of them are different from
However, for a given model fitting procedure, using different variants of
spline basis functions should result in identical prediction values. The
coefficient estimates of the spline basis functions returned by
are more interpretable compared to
naturalSpline( x, df = NULL, knots = NULL, intercept = FALSE, Boundary.knots = NULL, trim = 0, derivs = 0L, integral = FALSE, ... ) nsp( x, df = NULL, knots = NULL, intercept = FALSE, Boundary.knots = NULL, trim = 0, derivs = 0L, integral = FALSE, ... ) nsk( x, df = NULL, knots = NULL, intercept = FALSE, Boundary.knots = NULL, trim = 0, derivs = 0L, integral = FALSE, ... )
The predictor variable. Missing values are allowed and will be returned as they are.
Degree of freedom that equals to the column number of returned matrix. One can specify
knots, then the function chooses
df - 1 - as.integer(intercept)internal knots at suitable quantiles of
xignoring missing values and those
xoutside of the boundary. Thus,
dfmust be greater than or equal to
2. If internal knots are specified via
knots, the specified
dfwill be ignored.
The internal breakpoints that define the splines. The default is
NULL, which results in a basis for ordinary polynomial regression. Typical values are the mean or median for one knot, quantiles for more knots. For periodic splines, the number of knots must be greater or equal to the specified
degree - 1.
TRUE, the complete basis matrix will be returned. Otherwise, the first basis will be excluded from the output.
Boundary points at which to anchor the splines. By default, they are the range of
NA. If both
Boundary.knotsare supplied, the basis parameters do not depend on
x. Data can extend beyond
Boundary.knots. For periodic splines, the specified bounary knots define the cyclic interval.
The fraction (0 to 0.5) of observations to be trimmed from each end of
xbefore placing the default internal and boundary knots. This argument will be ignored if
Boundary.knotsis specified. The default value is
0for backward compatibility, which sets the boundary knots as the range of
x. If a positive fraction is specified, the default boundary knots will be equivalent to
quantile(x, probs = c(trim, 1 - trim), na.rm = TRUE), which can be a more sensible choice in practice due to the existence of outliers. The default internal knots are placed within the boundary afterwards.
A nonnegative integer specifying the order of derivatives of natural splines. The default value is
0Lfor the spline basis functions.
A logical value. The default value is
TRUE, this function will return the integrated natural splines from the left boundary knot.
Optional arguments that are not used.
A numeric matrix of
length(x) rows and
df is specified or
length(knots) + 1 +
as.integer(intercept) columns if
knots are specified instead.
Attributes that correspond to the arguments specified are returned for
usage of other functions in this package.
The constructed spline basis functions from
nonnegative within boundary with the second derivatives being zeros at
boundary knots. The implementation utilizes the close-form null space that
can be derived from the recursive formula for the second derivatives of
B-splines. The function
nsp() is an alias of
to encourage the use in a model formula.
nsk() produces another variant of natural cubic spline
matrix such that only one of the basis functions is nonzero and takes a
value of one at every boundary and internal knot. As a result, the
coefficients of the resulting fit are the values of the spline function at
the knots, which makes it easy to interpret the coefficient estimates. In
other words, the coefficients of a linear model will be the heights of the
function at the knots if
intercept = TRUE. If
FALSE, the coefficients will be the change in function value between each
knot. This implementation closely follows the function
nsk() of the
survival package (version 3.2-8). The idea corresponds directly to
the physical implementation of a spline by passing a flexible strip of wood
or metal through a set of fixed points, which is a traditional way to create
smooth shapes for things like a ship hull.
The returned basis matrix can be obtained by transforming the corresponding
B-spline basis matrix with the matrix
H provided in the attribute of
the returned object. Each basis is assumed to follow a linear trend for
x outside of boundary. A similar implementation is provided by
splines::ns, which uses QR decomposition to find the null space of
the second derivatives of B-spline basis at boundary knots. See
Supplementray Materials of Wang and Yan (2021) for a more detailed
library(splines2) x <- seq.int(0, 1, 0.01) knots <- c(0.3, 0.5, 0.6) ## naturalSpline() nsMat0 <- naturalSpline(x, knots = knots, intercept = TRUE) nsMat1 <- naturalSpline(x, knots = knots, intercept = TRUE, integral = TRUE) nsMat2 <- naturalSpline(x, knots = knots, intercept = TRUE, derivs = 1) nsMat3 <- naturalSpline(x, knots = knots, intercept = TRUE, derivs = 2) op <- par(mfrow = c(2, 2), mar = c(2.5, 2.5, 0.2, 0.1), mgp = c(1.5, 0.5, 0)) plot(nsMat0, ylab = "basis") plot(nsMat1, ylab = "integral") plot(nsMat2, ylab = "1st derivative") plot(nsMat3, ylab = "2nd derivative") par(op) # reset to previous plotting settings ## nsk() nskMat <- nsk(x, knots = knots, intercept = TRUE) plot(nskMat, ylab = "nsk()", mark_knots = "all") abline(h = 1, col = "red", lty = 3) ## use the deriv method all.equal(nsMat0, deriv(nsMat1), check.attributes = FALSE) #>  TRUE all.equal(nsMat2, deriv(nsMat0)) #>  TRUE all.equal(nsMat3, deriv(nsMat2)) #>  TRUE all.equal(nsMat3, deriv(nsMat0, 2)) #>  TRUE ## a linear model example fit1 <- lm(weight ~ -1 + nsk(height, df = 4, intercept = TRUE), data = women) fit2 <- lm(weight ~ nsk(height, df = 3), data = women) ## the knots (same for both fits) knots <- unlist(attributes(fit1$model[])[c('Boundary.knots', 'knots')]) ## predictions at the knot points predict(fit1, data.frame(height = sort(unname(knots)))) #> 1 2 3 4 #> 114.5595 127.9333 143.3219 163.6409 unname(coef(fit1)) # equal to the coefficient estimates #>  114.5595 127.9333 143.3219 163.6409 ## different interpretation when "intercept = FALSE" unname(coef(fit1)[-1] - coef(fit1)) # differences: yhat[2:4] - yhat #>  13.37381 28.76248 49.08146 unname(coef(fit2))[-1] # ditto #>  13.37381 28.76248 49.08146