Functions `naturalSpline()`

and `nsk()`

generate the natural cubic
spline basis functions, the corresponding derivatives or integrals (from the
left boundary knot). Both of them are different from `splines::ns()`

.
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 `nsk()`

are more interpretable compared to `naturalSpline()`

or
`splines::ns()`

.

## Usage

```
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,
...
)
```

## Arguments

- x
The predictor variable. Missing values are allowed and will be returned as they are.

- df
Degree of freedom that equals to the column number of returned matrix. One can specify

`df`

rather than`knots`

, then the function chooses`df - 1 - as.integer(intercept)`

internal knots at suitable quantiles of`x`

ignoring missing values and those`x`

outside of the boundary. Thus,`df`

must be greater than or equal to`2`

. If internal knots are specified via`knots`

, the specified`df`

will be ignored.- knots
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`

.- intercept
If

`TRUE`

, the complete basis matrix will be returned. Otherwise, the first basis will be excluded from the output.- Boundary.knots
Boundary points at which to anchor the splines. By default, they are the range of

`x`

excluding`NA`

. If both`knots`

and`Boundary.knots`

are 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.- trim
The fraction (0 to 0.5) of observations to be trimmed from each end of

`x`

before placing the default internal and boundary knots. This argument will be ignored if`Boundary.knots`

is specified. The default value is`0`

for 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.- derivs
A nonnegative integer specifying the order of derivatives of natural splines. The default value is

`0L`

for the spline basis functions.- integral
A logical value. The default value is

`FALSE`

. If`TRUE`

, this function will return the integrated natural splines from the left boundary knot.- ...
Optional arguments that are not used.

## Value

A numeric matrix of `length(x)`

rows and `df`

columns if `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.

## Details

The constructed spline basis functions from `naturalSpline()`

are
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 `naturalSpline()`

to encourage the use in a model formula.

The function `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 ```
intercept =
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
introduction.

## Examples

```
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)
#> [1] TRUE
all.equal(nsMat2, deriv(nsMat0))
#> [1] TRUE
all.equal(nsMat3, deriv(nsMat2))
#> [1] TRUE
all.equal(nsMat3, deriv(nsMat0, 2))
#> [1] 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[[2]])[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
#> [1] 114.5595 127.9333 143.3219 163.6409
## different interpretation when "intercept = FALSE"
unname(coef(fit1)[-1] - coef(fit1)[1]) # differences: yhat[2:4] - yhat[1]
#> [1] 13.37381 28.76248 49.08146
unname(coef(fit2))[-1] # ditto
#> [1] 13.37381 28.76248 49.08146
```