R/gspline.R
In john-d-fox/gspline: Generalized Regression Splines
Defines functions
knots.gspline
print.gspline_matrix
print.gspline
gspline
Documented in
gspline
# Stable splines:
# - centering and scaling X
# - orthogonalize spline basis
# - orthogonalize spline basis wrt to intercept
#
# 2019_12_10 - GM
# - fixed ortho_basis so it returns the right size of matrix
# if the input is not of full rank ... thus not returning
# a basis, but a basis with appended columns of 0's to fill
# the matrix so it has the same shape as the input
#
#' General Parametric Regression Splines With Variable Degrees and Smoothness
#'
#' This function creates functions (closures) that implement general polynomial
#' splines with possibly different degrees in each interval and different
#' orders of smoothness at each knot, including the possibility of allowing a
#' discontinuity at a knot.
#'
#' The function returned by \code{gspline} can be used to generate columns of a
#' model matrix representing the spline, or to generate portions of a linear
#' hypothesis matrix for estimates and Wald tests of features of the spline,
#' such as derivatives of various orders, discontinuities, etc., using the
#' \code{\link{wald}} function.
#'
#' Many polynomial regression splines can be generated by 'plus' functions
#' although the resulting basis for the spline may be ill conditioned. For
#' example, a 'quadratic spline' (a spline that is quadratic in each interval
#' with matching first derivatives at each knot) with knots at 1 and 3 can be
#' fit with:
#'
#' \code{plus <- function(x, y) ifelse(x > 0, y, 0)}\cr \code{ lm(y ~ x +
#' I(x^2) + plus(x - 1, (x - 1)^2) + plus(x - 3, (x - 3)^2))}
#'
#' All 'standard' polynomial splines with the same degree in each interval and
#' continuity of order one less than the degree at each knot can be constructed
#' in this fashion. A convenient aspect of this parametrization of the spline
#' is that the estimated coefficients have a simple interpretation. The
#' coefficients for 'plus' terms represent the 'saltus' (jump) in the value of
#' a coefficient at the knot. Testing whether the true value of the saltus in a
#' coefficient is 0 is equivalent to a test for the need to retain the
#' corresponding knot.
#'
#' This simple approach does not work for some more complex splines with
#' different degrees or different orders of continuity at the knots.
#'
#' Many techniques for fitting splines generate a basis for the spline (columns
#' of the model matrix) that has good numerical properties but whose
#' coefficients do not have a simple interpretation.
#'
#' The \code{gspline} function generates functions representing splines with
#' arbitrary degrees in each interval and arbitrary smoothness at each knot.
#' Any set of orders of derivatives can be contrained to be continuous at any
#' particular knot. The parametrization produces coefficients that have a
#' simple interpretation. For a spline of degree p at x = 0, coefficients
#' correspond to the 1st, 2nd, ..., pth derivative at 0. Additional
#' coefficients correspond to free discontinuities at each knot.
#'
#' A disadvantage of spline functions generated by \code{gspline} is that the
#' spline basis may be poorly conditioned. The impact of this problem can be
#' mitigated by rescaling the $x$ variable so that it has an appropriate range.
#' For example, with a spline whose highest degree is cubic, ensuring that $x$
#' has a range between -10 and 10 should avert numerical problems. The
#' \code{gspline} function makes provision for stabilizing the spline basis
#' that it produces both by rescaling $x$ and by linearly transforming the basis.
#'
#' Let \code{sp} denote a spline function created by \code{gspline}. For
#' example \code{sp <- gspline(knots = c(-1,0,2), degree = c(2,3,3,2),
#' smoothness = c(1,-1,1))} creates a spline function \code{sp} that can be
#' used in linear formulas to fit a polynomial spline with knots at -1, 0 and
#' 2, with degrees 2, 3, 3, and 2 in the four intervals bounded by the knots,
#' and with "C 1" continuity (continuous first derivative) at the knots -1 and
#' 2, and a possible discontinuity in value at the knot 0.
#'
#' The linear formula to fit this spline has the form \code{y ~ 1 + sp(x)}
#' which can be augmented with other regressors in the usual way.
#'
#' Called with a single argument (e.g., \code{sp(x)}), \code{sp} returns a
#' matrix of regressors and can be used in linear model formulas, e.g .,
#' \code{y ~ sp(x)}, \code{y ~ A * sp(x)}, or \code{y ~ A / sp(x) -1}, if
#' \code{A} is a factor, to fit separate splines within each level of \code{A}.
#'
#' If the optional first \code{x} argument to \code{gspline} is specified, then
#' by default the closure resturned by \code{gspline} will produce a spline basis
#' for \code{x} that's more numerically stable. In the preceding example, suppose
#' that \code{x} resides in the data frame \code{D}. Then
#' \code{sp <- gspline(D\$x, knots = c(-1,0,2), degree = c(2,3,3,2), smoothness = c(1,-1,1))}
#' produces a spline-generating function \code{sp} that builds a more numerically
#' stable spline basis for \code{D\$x}.
#'
#' The spline-generating function \code{sp} can be called with additional
#' arguments beyond \code{x} to geneate portions of linear hypothesis matrices
#' that can be used to test or estimate various values or derivatives of the spline:
#'
#' \describe{
#'
#' \item{\code{D}}{\emph{argument of the closure created by} \code{gspline}:
#' value(s) or the order of derivatives to be estimated when creating a linear
#' hypothesis matrix. The default is \code{D=0} corresponding to the value of
#' the spline.}
#'
#' \item{\code{limit}}{\emph{argument of the closure created by}
#' \code{gspline}: specifies whether a derivative or value should be estimated
#' as as limit from the left \code{limit = -1}, from the right, \code{limit =
#' +1}, or the difference of the two limits, \code{limit == 0}. The default is
#' \code{limit == +1}.}
#' }
#'
#' For example, \code{sp(c(-1, -1, -1, 0, 2), D = c(2, 2, 2, 0, 1), limit =
#' c(-1, 1, 0, 0, 1))} will return a matrix of coefficients to estimate,
#' respectively, the limits from the left and from the right of the second
#' derivative at -1, the difference (jump or saltus) in this derivative at -1,
#' the size of the discontinuity in the value of the function at 0, and the
#' first derivative at 1, which being continuous there has the same value
#' whether a limit is taken from the right or from the left.
#'
#' A function to fit a cubic spline with knots at 5 and 10 is generated with
#' \code{sp <- gspline(knots = c(5, 10), degree = c(3, 3, 3), smoothness = c(2,
#' 2))}, or more briefly: \code{sp <- gspline(knots=c(5, 10), degree=3, smoothness=2)} since the
#' \code{degree} and \code{smoothness} arguments are recycled as required by
#' the number of knots.
#'
#' A natural cubic spline with knots at 5 and 10 and boundary knots at 0 and 20
#' can be created with \code{sp <- gspline(knots=c(0, 5, 10, 20), degree=c(1, 3, 3, 3, 1),
#' smoothness=c(2, 2, 2, 2))}.
#'
#' A step function with steps at 0, 1 and 2 would have the form: \code{sp <-
#' gspline(knots=c(0, 1, 2), degree=0, smoothness=-1)}.
#'
#' @param x an optional vector of values that provide an estimate of the range of
#' values over which the closure returned by \code{gspline} will be evaluated.
#' It is only used if \code{rescale} and/or \code{stable} are \code{TRUE}.
#' @param knots vector of knots.
#' @param degree vector giving the degree of the spline in each interval. Note
#' the number of intervals is equal to the number of knots + 1. A value of 0
#' corresponds to a constant in the interval.
#' @param smoothness vector with the degree of smoothness at each knot (0 =
#' continuity, 1 = smooth with continuous first derivative, 2 = continuous
#' second derivative, etc. The value -1 allows a discontinuity at the knot.
#' Alternatively, a list each of whose elements specifies the derivatives that
#' must be continuous at each knot. This allow the set of continuous
#' derivatives to skip a value. For example, to specify that a function can
#' have a discontinuity but must have the same limiting slope and second
#' derivative (curvature) on either side of the discontinuity, the smoothness
#' vector for the corresponding knot would have the value \code{c(1, 2)}, thus
#' leaving the value (denoted by 0) unconstrained.
#' @param constraints provides a vector or matrix specifying additional linear
#' contraints on the 'full' parametrization consisting of blocks of polynomials
#' of degree equal to \code{max(degree)} in each of the \code{length(knots) +
#' 1} intervals of the spline. See below for an example of a spline that is 0
#' outside of its boundary knots.
#' @param estimates provides a vector or matrix specifying additional linear
#' combination(s) of the parameters in the 'full' parametrization that should
#' be estimated by estimated coefficients of the spline.
#' @param periodic if \code{TRUE} a periodic spline is generated on the base
#' interval [0, \code{max(knots)}]. A constraint is imposed so that the
#' coefficients generate, in effect, the same fitted values and derivatives to
#' the right of \code{max(knots)} as they do to the right of 0. Note that all
#' knots must be strictly positive. Since the knot at 0 corresponds to the knot
#' at \code{max(knots)}, the knot at 0 must not be explicitly specified. Note
#' also that the degree in the interval to the right of \code{max(knots)} must
#' be equal to that in the interval to the left of \code{min(knots)}. If these
#' are specified to be equal, \code{gspline} uses the maximum for both
#' intervals. The continuity constraints associated with the last knot are, in
#' effect, continuity constraints at 0.
#' @param intercept intercept value(s) of x at which the spline has value 0,
#' i.e., the value(s) of x for which y-hat is estimated by the intercept term
#' in the model. The default is 0. If \code{NULL}, the spline is not
#' constrained to evaluate to 0 for any x. This will usually result in a model
#' matrix that that is collinear with an intercept term.
#' @param tol.continuity,tol.basis zero tolerances: largest postive number indistinguishable from 0.
#' @param stable if \code{TRUE} (the default if the \code{x} argument is specified and
#' inapplicable if it is not), transform the spline basis into a equivalent
#' orthonormal basis, after rescaling if \code{rescale = TRUE}.
#'
#' @param rescale if \code{TRUE} (the default is the value of \code{stable}),
#' rescale the \code{x} argument so it has a range between
#' -10 and 10 in order to improve the numerical properties of the spline
#' basis.
#'
#' @param ortho2intercept if \code{TRUE} (the default is the value of \code{stable}),
#' the spline basis is replaced with a basis for the
#' orthogonal complement to the 1-vector in the space spanned by the
#' 1-vector and the spline basis. This should improve
#' the numerical properties of the fitted model but will produce incorrect
#' results for some models that don't satisfy the principle of marginality,
#' for example if a model does not include an intercept term.
#'
#' The estimated coefficients are not directly interpretable if \code{stable} or \code{rescale}
#' is \code{TRUE}, but the \code{\link{wald}} function is designed to compute and display
#' interpretable coefficients in this situation.
#' To see interpretable estimated coefficients use \code{wald(model)}. Also, hypothesis
#' matrices for the \code{\link{wald}} function are expressed as linear combinations of
#' the \emph{interpretable} coefficients.
#'
#' @param debug if \code{TRUE} (default is \code{FALSE}) print additional information.
#'
#' @return \code{gspline} returns a closure (function) that creates portions of
#' model matrices or of linear hypothesis matrices for Wald tests for a general
#' polynomial spline.
#'
#' @author Georges Monette and John Fox
#'
#' @examples
#'
#' ## Fitting a quadratic spline
#' set.seed(143634) # for reproducibility
#' Simd <- data.frame(age = rep(1:50, 2), y = sin(2*pi*(1:100)/5) + rnorm(100),
#' G = rep(c('male','female'), c(50,50)))
#' # define a function generating the spline quadratic spline with linear extrapolation
#' sp <- gspline(Simd$age, knots = c(10, 25, 40), degree = c(1, 2, 2, 1),
#' smooth = c(1, 1,1))
#'
#' fit <- lm(y ~ sp(age)*G, data = Simd)
#'
#' if (require(lattice)){
#' xyplot(predict(fit) ~ age, data = Simd, groups = G, type = 'l')
#' summary(fit)
#' }
#'
#' fitg <- update(fit, y ~ G/sp(age) - 1)
#' summary(fitg)
#'
#' ## Linear hypotheses
#' L <- list("Overall test of slopes at 20" = rbind(
#' "Female slope at age 20" = c( F20 <- cbind( 0 , sp(20, D = 1), 0 , 0 * sp(20, D = 1))),
#' "Male slope at age 20" = c( M20 <- cbind( 0 , sp(20, D = 1), 0 , 1 * sp(20, D = 1))),
#' "Difference" = c(M20 - F20))
#' )
#' wald(fit, L)
#'
#' ## Right and left second derivatives at knots and saltus (jump)
#'
#' L <- list("Second derivatives and saltus for female curve at knot at 25" =
#' cbind( 0, sp(c(25,25,25), D = 2, limit =c(-1,1,0)), 0,0,0,0))
#' L
#' wald(fit, L)
#'
#' L0 <- list(
#' "hat" = rbind(
#' "females at age=20" = c( 1, sp(20), 0, 0 * sp(20)),
#' "males at age=20" = c( 1, sp(20), 1, 1* sp(20))),
#' "male-female" = rbind(
#' "at 20" = c( 0 , 0*sp(20), 1, 1*sp(20))))
#' # wald(fit, L0) # !FIXME
#'
#' L1 <- list(
#' "D(yhat)/D(age)" =
#' rbind( "female at age = 25" = c(0, sp(25,1), 0, 0*sp(25,1)),
#' "male at x = 25" = c(0, sp(25,1), 0, 1*sp(25,1))))
#' wald( fit, L1)
#'
#' # Cubic spline:
#'
#' sp <- gspline(knots=c(5, 10), degree=c(3, 3, 3), smoothness=c(2, 2))
#'
#' # The parameters indicate that a cubic polynomial is used in each of the three intervals
#' # and that the second derivative is continuous at each knot.
#'
#' # Cubic natural spline:
#' # is a cubic spline in bounded intervals with linear components
#' # in each unbounded interval and continuous first derivative at the
#' # two knots for unbounded intervals.
#'
#' sp.natural <- gspline(knots=c(0, 5, 10, 15), degree=c(1, 3, 3, 3, 1),
#' smoothness=c(1, 2, 2, 1))
#'
#' # Quadratic and linear splines and step functions, respectively:
#'
#' sp.quad <- gspline(knots=c(5, 10), degree=c(2, 2, 2), smoothness=c(1, 1))
#' sp.lin <- gspline(knots=c(5, 10), degree=c(1, 1, 1), smoothness=c(0, 0))
#' sp.step <- gspline(knots=c(5, 10), degree=c(0, 0, 0), smoothness=c(-1, -1))
#'
#' # When the same degree is used for all
#' # intervals and knots, it suffices to give it once:
#'
#' sp.quad <- gspline(knots=c(5, 10), degree=2, smoothness=1)
#' sp.lin <- gspline(knots=c(5, 10), degree=1, smoothness=0)
#' sp.step <- gspline(knots=c(5, 10), degree=0, smoothness=-1)
#'
#' # An easy way to specify a model in which a knot is dropped
#' # is to force a degree of continuity equal to the degree of adjoining
#' # polynomials, e.g. to drop the knot at 10, use:
#'
#' sp.1 <- gspline(knots=c(5,10), degree=c(3, 3, 3), smoothness=c(2, 3))
#'
#' # This is sometimes easier than laboriously rewriting the
#' # spline function for each null hypothesis.
#'
#' # Depending on the maximal degree of the spline, the range of \code{x} should not be
#' # excessive to avoid numerical problems. The spline matrix generated is 'raw'
#' # and values of max(abs(x))^max(degree) may appear in the matrix. For
#' # example, for a cubic spline, it might be desirable to rescale x and/or
#' # recenter x so abs(x) < 100 if that is not already the case. Note that the
#' # knots need to be correspondingly transformed. This is done automatically if
#' # the argument \code{rescale=TRUE}.
#'
#' # The naming of coefficients should allow them to be easily interpreted. For
#' # example:
#'
#' sp <- gspline(knots=c(3, 7), degree=c(2, 3, 2), smoothness=c(1, 1))
#' sp(0:10, c(1, 1))
#'
#'
#' # The coefficient for the first regressor is the first derivative at x = 0;
#' # for the second regressor, the second derivative at 0; the third, the saltus
#' # (change) in the second derivative at x = 3, the fourth, the saltus in the
#' # third derivative at x = 3 and, finally, the saltus in the second derivative
#' # at x = 7.
#'
#' sp <- gspline(knots=c(3, 7) , degree=c(2, 3, 2), smoothness=c(1,2))
#' set.seed(27349) # for reproducibility
#' Zd <- data.frame(x = seq(0, 10, 0.5),
#' y = seq(0, 10, 0.5)^2 + rnorm(21))
#' fit <- lm( y ~ sp(x), data=Zd)
#' summary(fit)
#' Ls <- cbind(0, sp(c(1, 2, 3, 3, 3, 5, 7, 7, 7, 8), D = 2,
#' limit = c(-1, -1, -1, 1, 0, -1, -1, 1, 0, -1)))
#' zapsmall(Ls)
#'
#' # Note that estimates of features that are continuous at a point
#' # are indicated without a direction for the limit, e.g. "D2(1)"
#' # is the second derivative at 1.
#'
#' wald(fit, list('second derivatives' = Ls))
#'
#' # Note that some coefficients that are 0 by design may lead to invalid degrees
#' # of freedom and t-values.
#' @export
gspline <- function(x = NULL,
knots = NULL,
degree = 3,
smoothness = pmax(pmin(degree[-1], degree[-length(degree)]) - 1, 0),
intercept = 0,
constraints = NULL,
estimates = NULL,
periodic = FALSE,
stable = !missing(x) && !periodic,
rescale = stable,
ortho2intercept = stable,
tol.basis = sqrt(.Machine$double.eps),
tol.continuity = 1e-14,
debug = FALSE) {
if (!missing(x) && missing(knots))
knots <- quantile(x, 1:3 / 4)
#orthonormalize <- stable
basis <- function(X) {
# basis is a misnomer, this should be called 'span'
# returns linearly independent columns with possible pivoting
q <- qr(X, tol.basis = tol.basis)
sel <- q$pivot[seq_len(q$rank)]
ret <- X[, sel, drop = FALSE]
colnames(ret) <- colnames(X)[sel]
if(ncol(ret) < nrow(ret)) {
zeros <- matrix(0, nrow(ret), nrow(ret) - ncol(ret))
colnames(zeros) <- colnames(X)[-sel]
ret <- cbind(ret, zeros)
}
ret
}
orth_basis <- function(X) {
# returns orthonormal basis using the SVD
ud <- svd(X, nv = 0)
ret <- ud$u[, ud$d > tol.basis, drop = FALSE]
ret <- c(ret, rep(0,length(X) - length(ret)))
X[] <- ret
X
}
Cmat <- function(knots,
degree,
smoothness,
lin = NULL,
intercept = 0,
signif = 3) {
dm <- max(degree)
# intercept
cmat <- NULL
if (!is.null(intercept)) {
cmat <- rbind(cmat, "Intercept" = Xf(intercept, knots, dm, D = 0))
}
# continuity constraints
for (i in seq_along(knots)) {
k <- knots[i]
sm <- smoothness[[i]]
if (max(sm) > -1) {
# sm = -1 corresponds to discontinuity
if (!is.list(smoothness))
sm <- 0:sm
dmat <-
Xf(k, knots, dm, D = sm, FALSE) - Xf(k, knots, dm, D = sm, TRUE)
rownames(dmat) <- paste0("C", sm, '|', signif(k, signif))
cmat <- rbind(cmat, dmat)
}
}
# reduced degree constraints
degree <- rep(degree, length.out = length(knots) + 1)
for (i in seq_along(degree)) {
di <- degree[i]
if (dm > di) {
dmat <- diag((length(knots) + 1) *
(dm + 1)) [(i - 1) * (dm + 1) + 1 + seq(di + 1, dm), , drop = FALSE]
rownames(dmat) <-
paste("I.", i, ".", seq(di + 1, dm), sep = '')
cmat <- rbind(cmat, dmat)
}
}
# add additional linear constraint
if (!is.null(lin))
cmat <- rbind(cmat, lin)
rk <- qr(cmat)$rank
spline.rank <- ncol(cmat) - rk
attr(cmat, "ranks") = c(
npar.full = ncol(cmat),
C.n = nrow(cmat),
C.rank = rk,
spline.rank = spline.rank
)
attr(cmat, "d") <- svd(cmat)$d
cmat
}
Emat <- function(knots,
degree,
smoothness ,
intercept = FALSE,
signif = 3) {
if (length(degree) < length(knots) + 1) {
degree <- rep(degree, length.out = length(knots) + 1)
}
dmin <- min(degree)
dmax <- max(degree)
imax <- length(degree)
# derivatives for changes to estimate
# Quick temporary fix if smoothness is a list
# - generate extra terms that will be dropped due to collinearity with constraints
if (is.list(smoothness))
smoothness <- rep(-1, length(knots))
zeroi <-
as.numeric(cut(0, c(-Inf, sort(knots), Inf))) # index of interval containing 0
dzero <- degree[zeroi] # degree of interval containing 0
cmat <-
Xf(0, knots, degree = dmax, D = seq(1, dzero)) # derivatives at 0 - may be discarded
if (imax > zeroi) {
# intervals to the right of the interval containing 0
for (i in (zeroi + 1):imax) {
d.right <- degree[i]
d.left <- degree[i - 1]
k <- knots[i - 1]
sm <- smoothness[i - 1]
if (d.right > sm) {
dmat <- Xf(k , knots, dmax, D = seq(sm + 1, d.right) , FALSE) -
Xf(k , knots, dmax, D = seq(sm + 1, d.right) , TRUE)
rownames(dmat) <-
paste0("C", seq(sm + 1, d.right), '|', signif(k, signif))
cmat <- rbind(cmat, dmat)
}
}
}
if (zeroi > 1) {
for (i in zeroi:2) {
d.right <- degree[i]
d.left <- degree[i - 1]
k <- knots[i - 1]
sm <- smoothness[i - 1]
if (d.left > sm) {
dmat <- Xf(k , knots, dmax, D = seq(sm + 1, d.left) , FALSE) -
Xf(k , knots, dmax, D = seq(sm + 1, d.left) , TRUE)
rownames(dmat) <-
paste0("C", seq(sm + 1, d.left), '|', signif(k, signif))
cmat <- rbind(cmat, dmat)
}
}
}
cmat
}
Xmat <- function(x,
degree,
D = 0,
signif = 3) {
# Returns rows of design matrix if D = 0 or linear hypothesis matrix for D-th derivative
if (length(D) < length(x))
D <- rep(D, length.out = length(x))
if (length(x) < length(D))
x <- rep(x, length.out = length(D))
xmat <- matrix(x, nrow = length(x), ncol = degree + 1)
expvec <- 0:degree
coeffvec <- rep(1, degree + 1)
expmat <- NULL
coeffmat <- NULL
for (i in 0:max(D)) {
expmat <- rbind(expmat, expvec)
coeffmat <- rbind(coeffmat, coeffvec)
coeffvec <- coeffvec * expvec
expvec <- ifelse(expvec > 0, expvec - 1, 0)
}
X <-
coeffmat[D + 1, , drop = FALSE] * xmat ^ expmat[D + 1, , drop = FALSE]
xlab <- signif(x, signif)
rownames(X) <- ifelse(D == 0,
paste('f(', xlab , ')', sep = ''),
paste0("D", D, "|", xlab))
colnames(X) <- paste("X", 0:(ncol(X) - 1), sep = "")
class(X) <- c('gspline_matrix', class(X))
X
}
Xf <- function(x,
knots,
degree = 3,
D = 0,
right = TRUE,
periodic = FALSE,
signif = 3) {
# Returns block diagonal matrix (if x ordered) with design matrix
# or linear hypothesis matrix for the 'full' model with contrained parameters.
# With the default, right == TRUE, if x is at a knot then it is included in
# in the lower interval. With right = FALSE, it is included in the higher
# interval. This is needed when building derivative constraints at the knot
if (periodic) {
period <- max(knots)
xx <- x %% period
if (right)
xx[xx == 0] <- period
x <- xx
knots <- knots[-length(knots)]
}
xmat <- Xmat (x, degree, D , signif)
k <- sort(knots)
g <- cut(x, c(-Inf, k, Inf), right = right)
ret <- do.call('cbind', lapply(seq_along(levels(g)), function(i) (g == levels(g)[i]) * xmat))
if (periodic) {
rownames(ret) <- paste0(rownames(ret), '/', signif(period, signif))
}
class(ret) <- c('gspline_matrix', class(ret))
ret
}
# Value and derivatives at 0 and all discontinuities at knots
Dmat <- function(knots,
degree,
periodic = FALSE,
signif = 3) {
dm <- max(degree)
cmat <- Xf(0, knots, dm, D = 0:dm, periodic = periodic)
n_knots <- length(knots)
for (i in seq_len(n_knots - 1)) {
k <- knots[i]
dmat <- Xf(k,
knots,
dm,
D = seq(0, dm),
FALSE,
periodic = periodic) -
Xf(k,
knots,
dm,
D = seq(0, dm),
TRUE,
periodic = periodic)
rownames(dmat) <-
paste0("C", seq(0, dm), '|', signif(k, signif))
if (periodic) {
rownames(dmat) <-
paste0(rownames(dmat), '/', signif(max(k), signif))
}
cmat <- rbind(cmat, dmat)
}
k <- knots[length(knots)]
if (periodic) {
dmat <-
Xf(0,
knots,
dm,
D = seq(0, dm),
FALSE,
periodic = periodic) -
Xf(k,
knots,
dm,
D = seq(0, dm),
TRUE,
periodic = periodic)
rownames(dmat) <-
paste("C",
seq(0, dm),
'|',
signif(k, signif),
'/',
signif(max(k), signif))
cmat <- rbind(cmat, dmat)
} else {
dmat <-
Xf(k,
knots,
dm,
D = seq(0, dm),
FALSE,
periodic = periodic) -
Xf(k,
knots,
dm,
D = seq(0, dm),
TRUE,
periodic = periodic)
rownames(dmat) <-
paste0("C", seq(0, dm), '|', signif(k, signif))
cmat <- rbind(cmat, dmat)
}
class(cmat) <- c('gspline_matrix', class(cmat))
cmat
}
# Parameters to constrain
#
# TODO: Change so degree can be a list
Pcon <- function(knots, degree, periodic) {
degree <- rep(degree, length.out = length(knots) + 1)
if (periodic) {
if (degree[length(degree)] != degree[1]) {
warning("For periodic splines, the degree of the last and first intervals should match")
}
knots <- knots[-length(knots)]
degree <- degree[-length(degree)]
}
dm <- max(degree)
cmat <- NULL
for (i in seq_along(degree)) {
di <- degree[i]
if (dm > di) {
dmat <- diag((length(knots) + 1) *
(dm + 1)) [(i - 1) * (dm + 1) + 1 + seq(di + 1, dm), , drop = FALSE]
rownames(dmat) <-
paste("I.", i, ".", seq(di + 1, dm), sep = '')
cmat <- rbind(cmat, dmat)
}
}
if (!is.null(cmat)) {
class(cmat) <- c('gspline_matrix', class(cmat))
}
cmat
}
#
# Parse knots, degree and smoothness to create constraint and estimation matrices:
# - Identify rows of Dmat that are used for constraints and rows used for estimates
if (periodic) {
if (min(knots) <= 0)
stop('For a periodic spline all knots must be positive')
}
if (any(knots != sort(knots)))
stop('Knots must be in ascending order')
degree <- rep(degree, length.out = length(knots) + 1)
max_degree <- max(degree)
smoothness <- rep(smoothness, length.out = length(knots))
if (!is.list(smoothness)) {
smoothness <-
lapply(smoothness, function(n)
if (n < 0)
- 1
else
0:n)
}
Dmat_smoothness_indices <- unlist(lapply(seq_along(smoothness),
function(i)
smoothness[[i]] + (max_degree + 1) * i + 1))
Dmat_smoothness_indices <-
Dmat_smoothness_indices[unlist(smoothness) > -1]
# Build constraints
constraint_mat <-
Dmat(knots, degree, periodic = periodic)[Dmat_smoothness_indices, , drop = FALSE]
constraint_mat <-
rbind(constraint_mat, Pcon(knots, degree, periodic = periodic))
if (!is.null(constraints))
constraint_mat <- rbind(constraint_mat, constraints)
if (!is.null(intercept)) {
constraint_mat <-
rbind(constraint_mat,
Xf(intercept, knots, max_degree, periodic = periodic))
Dmat_smoothness_indices <- c(1, Dmat_smoothness_indices)
}
cmat <- t(basis(t(constraint_mat)))
estimate_mat <-
rbind(estimates, Dmat(knots, degree, periodic = periodic)
[-Dmat_smoothness_indices, , drop = FALSE])
# the following is necessary in case some of the derivatives
# estimated at 0 have been coerced to 0 by other constraints
estimate_mat <- t(basis(t(rbind(
constraint_mat, estimate_mat
))))
estimate_mat <-
estimate_mat[-seq_len(nrow(constraint_mat)), , drop = FALSE]
# emat <- estimate_mat #FIXME not used JF
A <- rbind(constraint_mat, estimate_mat)
G <- solve(rbind(constraint_mat, estimate_mat),
diag(ncol(constraint_mat))[, -seq_len(nrow(constraint_mat)), drop = FALSE])
colnames(G) <- rownames(estimate_mat)
# create closure
fun <- function(x = NA, D = NULL, limit = 1) {
if (is.null(D)) {
ret <-
Xf(x, knots, max_degree, periodic = periodic) %*% G # model matrix
} else {
# Hypothesis matrix
D <- rep(D, length.out = length(x))
limit <- rep(limit, length.out = length(x))
left <- Xf(
x,
knots = knots,
degree = max_degree,
D = D,
periodic = periodic,
right = TRUE
)
right <- Xf(
x,
knots = knots,
degree = max_degree,
D = D,
periodic = periodic,
right = FALSE
)
cleft <- c(1, 0, -1)[match(limit, c(-1, 1, 0))]
cright <- c(0, 1, 1)[match(limit, c(-1, 1, 0))]
continuous <-
apply((right - left) %*% G, 1, function(x)
all(abs(x) < 10 * tol.continuity))
raw <- left * cleft + right * cright
ret <- raw %*% G
nam <- rownames(ret)
nam <- sub("^f", "g", nam)
nam_left <- sub("(/|$)", "-\\1", nam)
nam_right <- sub("(/|$)", "+\\1", nam)
nam_jump <- paste0(nam_right , "-", nam_left)
# If we're testing for a jump and the derivative is continuous we should show
# the identity of the jump tested.
#
# If we're testing for a derivative at a knot and the derivative is continuous, we should
# indicate that by removing the direction indicator
#
# So:
# If we're testing a jump we use the full name of the linear combination even
# if it happens to be at a point of continuity. Otherwise, used the directional notation only
# if the tested derivative is discontinuous
rownames(ret) <- ifelse(limit == 0, nam_jump,
ifelse(continuous, nam, ifelse(limit == 1, nam_right, nam_left)))
}
class(ret) <- c('gspline_matrix', class(ret))
ret
}
# Modify colnames of G matrix to show direction of limit
# if there is a discontinuity of derivatives at origin
if (!is.null(intercept)) {
dest <- fun(rep(intercept, max(degree) + 1),
D = 0:max(degree),
limit = -1)
tos <- grep('-$|-/', rownames(dest), value = TRUE)
froms <- sub('([0-9])-', '\\1', tos)
froms <- sub('(.*)', '^\\1$', froms)
froms <- sub('\\|', '\\\\|', froms)
for (i in seq_along(tos)) {
colnames(G) <- sub(froms[i], tos[i], colnames(G))
}
}
ffun <- if (!missing(x) && (stable || rescale)) {
if (rescale) {
if(!periodic) {
ran <- range(x)
a <- mean(ran)
b <- 20 / (ran[2] - ran[1]) # rescaled x goes from -10 to 10
} else { # if periodic
a <- 0
b <- 10 / knots[length(knots)]
}
fun_scaled <- gspline(
knots = (knots - a) * b,
degree = degree,
smoothness = smoothness,
periodic = periodic,
intercept = (intercept - a) * b,
stable = FALSE
)
}
function(x = NA, D = NULL, limit = 1) {
if (inwald() || !is.null(D)) {
ret <- fun(x, D, limit)
} else {
ret <- if (rescale)
fun_scaled((x - a) * b)
else
fun(x)
colnames(ret) <- NULL
# make basis orthogonal to intercept term GM:2019_10_11
if(ortho2intercept && !periodic)
ret <- sweep(ret, 2, colMeans(ret))
if (stable)
ret <- orth_basis(ret)
}
if (rescale)
attr(ret, 'rescaled') <- c(a = a, b = b)
if (stable)
attr(ret, 'orthonormalize') <- TRUE
if (ortho2intercept)
attr(ret, 'ortho2intercept') <- TRUE
class(ret) <- c('gspline_matrix', class(ret))
ret
}
} else
fun
class(ffun) <- c('gspline', class(ffun))
ffun
}
#' @method print gspline
#' @export
print.gspline <- function(x,
show = c('knots',
'degree',
'smoothness',
'G',
'constraint_mat',
'estimate_mat'),
...) {
cat('Spline function created by gspline\n')
envr <- environment(x)
ret <-
Filter(Negate(is.function), sapply(ls(environment(x)), get, environment(x)))
ret <-
lapply(ret, function(x)
if (is.numeric(x))
zapsmall(x)
else
x)
if (!is.null(show))
ret <- ret[show]
print(ret)
invisible(x)
}
#' @method print gspline_matrix
#' @export
print.gspline_matrix <- function(x, ...) {
header <- paste(if (!is.null(attr(x, "rescaled"))) "rescaled = TRUE",
if (!is.null(attr(x, "orthonormalize"))) "stable = TRUE",
if (!is.null(attr(x, "ortho2intercept"))) "ortho2intercept = TRUE",
sep="; ")
if (!is.null(header)) cat("\n", header, "\n")
xx <- zapsmall(x)
class(xx) <- "matrix"
attr(xx, "rescaled") <- NULL
attr(xx, "orthonormalize") <- NULL
attr(xx,"ortho2intercept") <- NULL
print(xx)
invisible(x)
}
#' @method knots gspline
#' @export
knots.gspline <- function(Fn, ...) {
environment(Fn)$knots
}
john-d-fox/gspline documentation built on Sept. 17, 2020, 3:19 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions knots.gspline print.gspline_matrix print.gspline gspline
Documented in gspline
# Stable splines:
# - centering and scaling X
# - orthogonalize spline basis
# - orthogonalize spline basis wrt to intercept
#
# 2019_12_10 - GM
# - fixed ortho_basis so it returns the right size of matrix
# if the input is not of full rank ... thus not returning
# a basis, but a basis with appended columns of 0's to fill
# the matrix so it has the same shape as the input
#
#' General Parametric Regression Splines With Variable Degrees and Smoothness
#'
#' This function creates functions (closures) that implement general polynomial
#' splines with possibly different degrees in each interval and different
#' orders of smoothness at each knot, including the possibility of allowing a
#' discontinuity at a knot.
#'
#' The function returned by \code{gspline} can be used to generate columns of a
#' model matrix representing the spline, or to generate portions of a linear
#' hypothesis matrix for estimates and Wald tests of features of the spline,
#' such as derivatives of various orders, discontinuities, etc., using the
#' \code{\link{wald}} function.
#'
#' Many polynomial regression splines can be generated by 'plus' functions
#' although the resulting basis for the spline may be ill conditioned. For
#' example, a 'quadratic spline' (a spline that is quadratic in each interval
#' with matching first derivatives at each knot) with knots at 1 and 3 can be
#' fit with:
#'
#' \code{plus <- function(x, y) ifelse(x > 0, y, 0)}\cr \code{ lm(y ~ x +
#' I(x^2) + plus(x - 1, (x - 1)^2) + plus(x - 3, (x - 3)^2))}
#'
#' All 'standard' polynomial splines with the same degree in each interval and
#' continuity of order one less than the degree at each knot can be constructed
#' in this fashion. A convenient aspect of this parametrization of the spline
#' is that the estimated coefficients have a simple interpretation. The
#' coefficients for 'plus' terms represent the 'saltus' (jump) in the value of
#' a coefficient at the knot. Testing whether the true value of the saltus in a
#' coefficient is 0 is equivalent to a test for the need to retain the
#' corresponding knot.
#'
#' This simple approach does not work for some more complex splines with
#' different degrees or different orders of continuity at the knots.
#'
#' Many techniques for fitting splines generate a basis for the spline (columns
#' of the model matrix) that has good numerical properties but whose
#' coefficients do not have a simple interpretation.
#'
#' The \code{gspline} function generates functions representing splines with
#' arbitrary degrees in each interval and arbitrary smoothness at each knot.
#' Any set of orders of derivatives can be contrained to be continuous at any
#' particular knot. The parametrization produces coefficients that have a
#' simple interpretation. For a spline of degree p at x = 0, coefficients
#' correspond to the 1st, 2nd, ..., pth derivative at 0. Additional
#' coefficients correspond to free discontinuities at each knot.
#'
#' A disadvantage of spline functions generated by \code{gspline} is that the
#' spline basis may be poorly conditioned. The impact of this problem can be
#' mitigated by rescaling the $x$ variable so that it has an appropriate range.
#' For example, with a spline whose highest degree is cubic, ensuring that $x$
#' has a range between -10 and 10 should avert numerical problems. The
#' \code{gspline} function makes provision for stabilizing the spline basis
#' that it produces both by rescaling $x$ and by linearly transforming the basis.
#'
#' Let \code{sp} denote a spline function created by \code{gspline}. For
#' example \code{sp <- gspline(knots = c(-1,0,2), degree = c(2,3,3,2),
#' smoothness = c(1,-1,1))} creates a spline function \code{sp} that can be
#' used in linear formulas to fit a polynomial spline with knots at -1, 0 and
#' 2, with degrees 2, 3, 3, and 2 in the four intervals bounded by the knots,
#' and with "C 1" continuity (continuous first derivative) at the knots -1 and
#' 2, and a possible discontinuity in value at the knot 0.
#'
#' The linear formula to fit this spline has the form \code{y ~ 1 + sp(x)}
#' which can be augmented with other regressors in the usual way.
#'
#' Called with a single argument (e.g., \code{sp(x)}), \code{sp} returns a
#' matrix of regressors and can be used in linear model formulas, e.g .,
#' \code{y ~ sp(x)}, \code{y ~ A * sp(x)}, or \code{y ~ A / sp(x) -1}, if
#' \code{A} is a factor, to fit separate splines within each level of \code{A}.
#'
#' If the optional first \code{x} argument to \code{gspline} is specified, then
#' by default the closure resturned by \code{gspline} will produce a spline basis
#' for \code{x} that's more numerically stable. In the preceding example, suppose
#' that \code{x} resides in the data frame \code{D}. Then
#' \code{sp <- gspline(D\$x, knots = c(-1,0,2), degree = c(2,3,3,2), smoothness = c(1,-1,1))}
#' produces a spline-generating function \code{sp} that builds a more numerically
#' stable spline basis for \code{D\$x}.
#'
#' The spline-generating function \code{sp} can be called with additional
#' arguments beyond \code{x} to geneate portions of linear hypothesis matrices
#' that can be used to test or estimate various values or derivatives of the spline:
#'
#' \describe{
#'
#' \item{\code{D}}{\emph{argument of the closure created by} \code{gspline}:
#' value(s) or the order of derivatives to be estimated when creating a linear
#' hypothesis matrix. The default is \code{D=0} corresponding to the value of
#' the spline.}
#'
#' \item{\code{limit}}{\emph{argument of the closure created by}
#' \code{gspline}: specifies whether a derivative or value should be estimated
#' as as limit from the left \code{limit = -1}, from the right, \code{limit =
#' +1}, or the difference of the two limits, \code{limit == 0}. The default is
#' \code{limit == +1}.}
#' }
#'
#' For example, \code{sp(c(-1, -1, -1, 0, 2), D = c(2, 2, 2, 0, 1), limit =
#' c(-1, 1, 0, 0, 1))} will return a matrix of coefficients to estimate,
#' respectively, the limits from the left and from the right of the second
#' derivative at -1, the difference (jump or saltus) in this derivative at -1,
#' the size of the discontinuity in the value of the function at 0, and the
#' first derivative at 1, which being continuous there has the same value
#' whether a limit is taken from the right or from the left.
#'
#' A function to fit a cubic spline with knots at 5 and 10 is generated with
#' \code{sp <- gspline(knots = c(5, 10), degree = c(3, 3, 3), smoothness = c(2,
#' 2))}, or more briefly: \code{sp <- gspline(knots=c(5, 10), degree=3, smoothness=2)} since the
#' \code{degree} and \code{smoothness} arguments are recycled as required by
#' the number of knots.
#'
#' A natural cubic spline with knots at 5 and 10 and boundary knots at 0 and 20
#' can be created with \code{sp <- gspline(knots=c(0, 5, 10, 20), degree=c(1, 3, 3, 3, 1),
#' smoothness=c(2, 2, 2, 2))}.
#'
#' A step function with steps at 0, 1 and 2 would have the form: \code{sp <-
#' gspline(knots=c(0, 1, 2), degree=0, smoothness=-1)}.
#'
#' @param x an optional vector of values that provide an estimate of the range of
#' values over which the closure returned by \code{gspline} will be evaluated.
#' It is only used if \code{rescale} and/or \code{stable} are \code{TRUE}.
#' @param knots vector of knots.
#' @param degree vector giving the degree of the spline in each interval. Note
#' the number of intervals is equal to the number of knots + 1. A value of 0
#' corresponds to a constant in the interval.
#' @param smoothness vector with the degree of smoothness at each knot (0 =
#' continuity, 1 = smooth with continuous first derivative, 2 = continuous
#' second derivative, etc. The value -1 allows a discontinuity at the knot.
#' Alternatively, a list each of whose elements specifies the derivatives that
#' must be continuous at each knot. This allow the set of continuous
#' derivatives to skip a value. For example, to specify that a function can
#' have a discontinuity but must have the same limiting slope and second
#' derivative (curvature) on either side of the discontinuity, the smoothness
#' vector for the corresponding knot would have the value \code{c(1, 2)}, thus
#' leaving the value (denoted by 0) unconstrained.
#' @param constraints provides a vector or matrix specifying additional linear
#' contraints on the 'full' parametrization consisting of blocks of polynomials
#' of degree equal to \code{max(degree)} in each of the \code{length(knots) +
#' 1} intervals of the spline. See below for an example of a spline that is 0
#' outside of its boundary knots.
#' @param estimates provides a vector or matrix specifying additional linear
#' combination(s) of the parameters in the 'full' parametrization that should
#' be estimated by estimated coefficients of the spline.
#' @param periodic if \code{TRUE} a periodic spline is generated on the base
#' interval [0, \code{max(knots)}]. A constraint is imposed so that the
#' coefficients generate, in effect, the same fitted values and derivatives to
#' the right of \code{max(knots)} as they do to the right of 0. Note that all
#' knots must be strictly positive. Since the knot at 0 corresponds to the knot
#' at \code{max(knots)}, the knot at 0 must not be explicitly specified. Note
#' also that the degree in the interval to the right of \code{max(knots)} must
#' be equal to that in the interval to the left of \code{min(knots)}. If these
#' are specified to be equal, \code{gspline} uses the maximum for both
#' intervals. The continuity constraints associated with the last knot are, in
#' effect, continuity constraints at 0.
#' @param intercept intercept value(s) of x at which the spline has value 0,
#' i.e., the value(s) of x for which y-hat is estimated by the intercept term
#' in the model. The default is 0. If \code{NULL}, the spline is not
#' constrained to evaluate to 0 for any x. This will usually result in a model
#' matrix that that is collinear with an intercept term.
#' @param tol.continuity,tol.basis zero tolerances: largest postive number indistinguishable from 0.
#' @param stable if \code{TRUE} (the default if the \code{x} argument is specified and
#' inapplicable if it is not), transform the spline basis into a equivalent
#' orthonormal basis, after rescaling if \code{rescale = TRUE}.
#'
#' @param rescale if \code{TRUE} (the default is the value of \code{stable}),
#' rescale the \code{x} argument so it has a range between
#' -10 and 10 in order to improve the numerical properties of the spline
#' basis.
#'
#' @param ortho2intercept if \code{TRUE} (the default is the value of \code{stable}),
#' the spline basis is replaced with a basis for the
#' orthogonal complement to the 1-vector in the space spanned by the
#' 1-vector and the spline basis. This should improve
#' the numerical properties of the fitted model but will produce incorrect
#' results for some models that don't satisfy the principle of marginality,
#' for example if a model does not include an intercept term.
#'
#' The estimated coefficients are not directly interpretable if \code{stable} or \code{rescale}
#' is \code{TRUE}, but the \code{\link{wald}} function is designed to compute and display
#' interpretable coefficients in this situation.
#' To see interpretable estimated coefficients use \code{wald(model)}. Also, hypothesis
#' matrices for the \code{\link{wald}} function are expressed as linear combinations of
#' the \emph{interpretable} coefficients.
#'
#' @param debug if \code{TRUE} (default is \code{FALSE}) print additional information.
#'
#' @return \code{gspline} returns a closure (function) that creates portions of
#' model matrices or of linear hypothesis matrices for Wald tests for a general
#' polynomial spline.
#'
#' @author Georges Monette and John Fox
#'
#' @examples
#'
#' ## Fitting a quadratic spline
#' set.seed(143634) # for reproducibility
#' Simd <- data.frame(age = rep(1:50, 2), y = sin(2*pi*(1:100)/5) + rnorm(100),
#' G = rep(c('male','female'), c(50,50)))
#' # define a function generating the spline quadratic spline with linear extrapolation
#' sp <- gspline(Simd$age, knots = c(10, 25, 40), degree = c(1, 2, 2, 1),
#' smooth = c(1, 1,1))
#'
#' fit <- lm(y ~ sp(age)*G, data = Simd)
#'
#' if (require(lattice)){
#' xyplot(predict(fit) ~ age, data = Simd, groups = G, type = 'l')
#' summary(fit)
#' }
#'
#' fitg <- update(fit, y ~ G/sp(age) - 1)
#' summary(fitg)
#'
#' ## Linear hypotheses
#' L <- list("Overall test of slopes at 20" = rbind(
#' "Female slope at age 20" = c( F20 <- cbind( 0 , sp(20, D = 1), 0 , 0 * sp(20, D = 1))),
#' "Male slope at age 20" = c( M20 <- cbind( 0 , sp(20, D = 1), 0 , 1 * sp(20, D = 1))),
#' "Difference" = c(M20 - F20))
#' )
#' wald(fit, L)
#'
#' ## Right and left second derivatives at knots and saltus (jump)
#'
#' L <- list("Second derivatives and saltus for female curve at knot at 25" =
#' cbind( 0, sp(c(25,25,25), D = 2, limit =c(-1,1,0)), 0,0,0,0))
#' L
#' wald(fit, L)
#'
#' L0 <- list(
#' "hat" = rbind(
#' "females at age=20" = c( 1, sp(20), 0, 0 * sp(20)),
#' "males at age=20" = c( 1, sp(20), 1, 1* sp(20))),
#' "male-female" = rbind(
#' "at 20" = c( 0 , 0*sp(20), 1, 1*sp(20))))
#' # wald(fit, L0) # !FIXME
#'
#' L1 <- list(
#' "D(yhat)/D(age)" =
#' rbind( "female at age = 25" = c(0, sp(25,1), 0, 0*sp(25,1)),
#' "male at x = 25" = c(0, sp(25,1), 0, 1*sp(25,1))))
#' wald( fit, L1)
#'
#' # Cubic spline:
#'
#' sp <- gspline(knots=c(5, 10), degree=c(3, 3, 3), smoothness=c(2, 2))
#'
#' # The parameters indicate that a cubic polynomial is used in each of the three intervals
#' # and that the second derivative is continuous at each knot.
#'
#' # Cubic natural spline:
#' # is a cubic spline in bounded intervals with linear components
#' # in each unbounded interval and continuous first derivative at the
#' # two knots for unbounded intervals.
#'
#' sp.natural <- gspline(knots=c(0, 5, 10, 15), degree=c(1, 3, 3, 3, 1),
#' smoothness=c(1, 2, 2, 1))
#'
#' # Quadratic and linear splines and step functions, respectively:
#'
#' sp.quad <- gspline(knots=c(5, 10), degree=c(2, 2, 2), smoothness=c(1, 1))
#' sp.lin <- gspline(knots=c(5, 10), degree=c(1, 1, 1), smoothness=c(0, 0))
#' sp.step <- gspline(knots=c(5, 10), degree=c(0, 0, 0), smoothness=c(-1, -1))
#'
#' # When the same degree is used for all
#' # intervals and knots, it suffices to give it once:
#'
#' sp.quad <- gspline(knots=c(5, 10), degree=2, smoothness=1)
#' sp.lin <- gspline(knots=c(5, 10), degree=1, smoothness=0)
#' sp.step <- gspline(knots=c(5, 10), degree=0, smoothness=-1)
#'
#' # An easy way to specify a model in which a knot is dropped
#' # is to force a degree of continuity equal to the degree of adjoining
#' # polynomials, e.g. to drop the knot at 10, use:
#'
#' sp.1 <- gspline(knots=c(5,10), degree=c(3, 3, 3), smoothness=c(2, 3))
#'
#' # This is sometimes easier than laboriously rewriting the
#' # spline function for each null hypothesis.
#'
#' # Depending on the maximal degree of the spline, the range of \code{x} should not be
#' # excessive to avoid numerical problems. The spline matrix generated is 'raw'
#' # and values of max(abs(x))^max(degree) may appear in the matrix. For
#' # example, for a cubic spline, it might be desirable to rescale x and/or
#' # recenter x so abs(x) < 100 if that is not already the case. Note that the
#' # knots need to be correspondingly transformed. This is done automatically if
#' # the argument \code{rescale=TRUE}.
#'
#' # The naming of coefficients should allow them to be easily interpreted. For
#' # example:
#'
#' sp <- gspline(knots=c(3, 7), degree=c(2, 3, 2), smoothness=c(1, 1))
#' sp(0:10, c(1, 1))
#'
#'
#' # The coefficient for the first regressor is the first derivative at x = 0;
#' # for the second regressor, the second derivative at 0; the third, the saltus
#' # (change) in the second derivative at x = 3, the fourth, the saltus in the
#' # third derivative at x = 3 and, finally, the saltus in the second derivative
#' # at x = 7.
#'
#' sp <- gspline(knots=c(3, 7) , degree=c(2, 3, 2), smoothness=c(1,2))
#' set.seed(27349) # for reproducibility
#' Zd <- data.frame(x = seq(0, 10, 0.5),
#' y = seq(0, 10, 0.5)^2 + rnorm(21))
#' fit <- lm( y ~ sp(x), data=Zd)
#' summary(fit)
#' Ls <- cbind(0, sp(c(1, 2, 3, 3, 3, 5, 7, 7, 7, 8), D = 2,
#' limit = c(-1, -1, -1, 1, 0, -1, -1, 1, 0, -1)))
#' zapsmall(Ls)
#'
#' # Note that estimates of features that are continuous at a point
#' # are indicated without a direction for the limit, e.g. "D2(1)"
#' # is the second derivative at 1.
#'
#' wald(fit, list('second derivatives' = Ls))
#'
#' # Note that some coefficients that are 0 by design may lead to invalid degrees
#' # of freedom and t-values.
#' @export
gspline <- function(x = NULL,
knots = NULL,
degree = 3,
smoothness = pmax(pmin(degree[-1], degree[-length(degree)]) - 1, 0),
intercept = 0,
constraints = NULL,
estimates = NULL,
periodic = FALSE,
stable = !missing(x) && !periodic,
rescale = stable,
ortho2intercept = stable,
tol.basis = sqrt(.Machine$double.eps),
tol.continuity = 1e-14,
debug = FALSE) {
if (!missing(x) && missing(knots))
knots <- quantile(x, 1:3 / 4)
#orthonormalize <- stable
basis <- function(X) {
# basis is a misnomer, this should be called 'span'
# returns linearly independent columns with possible pivoting
q <- qr(X, tol.basis = tol.basis)
sel <- q$pivot[seq_len(q$rank)]
ret <- X[, sel, drop = FALSE]
colnames(ret) <- colnames(X)[sel]
if(ncol(ret) < nrow(ret)) {
zeros <- matrix(0, nrow(ret), nrow(ret) - ncol(ret))
colnames(zeros) <- colnames(X)[-sel]
ret <- cbind(ret, zeros)
}
ret
}
orth_basis <- function(X) {
# returns orthonormal basis using the SVD
ud <- svd(X, nv = 0)
ret <- ud$u[, ud$d > tol.basis, drop = FALSE]
ret <- c(ret, rep(0,length(X) - length(ret)))
X[] <- ret
X
}
Cmat <- function(knots,
degree,
smoothness,
lin = NULL,
intercept = 0,
signif = 3) {
dm <- max(degree)
# intercept
cmat <- NULL
if (!is.null(intercept)) {
cmat <- rbind(cmat, "Intercept" = Xf(intercept, knots, dm, D = 0))
}
# continuity constraints
for (i in seq_along(knots)) {
k <- knots[i]
sm <- smoothness[[i]]
if (max(sm) > -1) {
# sm = -1 corresponds to discontinuity
if (!is.list(smoothness))
sm <- 0:sm
dmat <-
Xf(k, knots, dm, D = sm, FALSE) - Xf(k, knots, dm, D = sm, TRUE)
rownames(dmat) <- paste0("C", sm, '|', signif(k, signif))
cmat <- rbind(cmat, dmat)
}
}
# reduced degree constraints
degree <- rep(degree, length.out = length(knots) + 1)
for (i in seq_along(degree)) {
di <- degree[i]
if (dm > di) {
dmat <- diag((length(knots) + 1) *
(dm + 1)) [(i - 1) * (dm + 1) + 1 + seq(di + 1, dm), , drop = FALSE]
rownames(dmat) <-
paste("I.", i, ".", seq(di + 1, dm), sep = '')
cmat <- rbind(cmat, dmat)
}
}
# add additional linear constraint
if (!is.null(lin))
cmat <- rbind(cmat, lin)
rk <- qr(cmat)$rank
spline.rank <- ncol(cmat) - rk
attr(cmat, "ranks") = c(
npar.full = ncol(cmat),
C.n = nrow(cmat),
C.rank = rk,
spline.rank = spline.rank
)
attr(cmat, "d") <- svd(cmat)$d
cmat
}
Emat <- function(knots,
degree,
smoothness ,
intercept = FALSE,
signif = 3) {
if (length(degree) < length(knots) + 1) {
degree <- rep(degree, length.out = length(knots) + 1)
}
dmin <- min(degree)
dmax <- max(degree)
imax <- length(degree)
# derivatives for changes to estimate
# Quick temporary fix if smoothness is a list
# - generate extra terms that will be dropped due to collinearity with constraints
if (is.list(smoothness))
smoothness <- rep(-1, length(knots))
zeroi <-
as.numeric(cut(0, c(-Inf, sort(knots), Inf))) # index of interval containing 0
dzero <- degree[zeroi] # degree of interval containing 0
cmat <-
Xf(0, knots, degree = dmax, D = seq(1, dzero)) # derivatives at 0 - may be discarded
if (imax > zeroi) {
# intervals to the right of the interval containing 0
for (i in (zeroi + 1):imax) {
d.right <- degree[i]
d.left <- degree[i - 1]
k <- knots[i - 1]
sm <- smoothness[i - 1]
if (d.right > sm) {
dmat <- Xf(k , knots, dmax, D = seq(sm + 1, d.right) , FALSE) -
Xf(k , knots, dmax, D = seq(sm + 1, d.right) , TRUE)
rownames(dmat) <-
paste0("C", seq(sm + 1, d.right), '|', signif(k, signif))
cmat <- rbind(cmat, dmat)
}
}
}
if (zeroi > 1) {
for (i in zeroi:2) {
d.right <- degree[i]
d.left <- degree[i - 1]
k <- knots[i - 1]
sm <- smoothness[i - 1]
if (d.left > sm) {
dmat <- Xf(k , knots, dmax, D = seq(sm + 1, d.left) , FALSE) -
Xf(k , knots, dmax, D = seq(sm + 1, d.left) , TRUE)
rownames(dmat) <-
paste0("C", seq(sm + 1, d.left), '|', signif(k, signif))
cmat <- rbind(cmat, dmat)
}
}
}
cmat
}
Xmat <- function(x,
degree,
D = 0,
signif = 3) {
# Returns rows of design matrix if D = 0 or linear hypothesis matrix for D-th derivative
if (length(D) < length(x))
D <- rep(D, length.out = length(x))
if (length(x) < length(D))
x <- rep(x, length.out = length(D))
xmat <- matrix(x, nrow = length(x), ncol = degree + 1)
expvec <- 0:degree
coeffvec <- rep(1, degree + 1)
expmat <- NULL
coeffmat <- NULL
for (i in 0:max(D)) {
expmat <- rbind(expmat, expvec)
coeffmat <- rbind(coeffmat, coeffvec)
coeffvec <- coeffvec * expvec
expvec <- ifelse(expvec > 0, expvec - 1, 0)
}
X <-
coeffmat[D + 1, , drop = FALSE] * xmat ^ expmat[D + 1, , drop = FALSE]
xlab <- signif(x, signif)
rownames(X) <- ifelse(D == 0,
paste('f(', xlab , ')', sep = ''),
paste0("D", D, "|", xlab))
colnames(X) <- paste("X", 0:(ncol(X) - 1), sep = "")
class(X) <- c('gspline_matrix', class(X))
X
}
Xf <- function(x,
knots,
degree = 3,
D = 0,
right = TRUE,
periodic = FALSE,
signif = 3) {
# Returns block diagonal matrix (if x ordered) with design matrix
# or linear hypothesis matrix for the 'full' model with contrained parameters.
# With the default, right == TRUE, if x is at a knot then it is included in
# in the lower interval. With right = FALSE, it is included in the higher
# interval. This is needed when building derivative constraints at the knot
if (periodic) {
period <- max(knots)
xx <- x %% period
if (right)
xx[xx == 0] <- period
x <- xx
knots <- knots[-length(knots)]
}
xmat <- Xmat (x, degree, D , signif)
k <- sort(knots)
g <- cut(x, c(-Inf, k, Inf), right = right)
ret <- do.call('cbind', lapply(seq_along(levels(g)), function(i) (g == levels(g)[i]) * xmat))
if (periodic) {
rownames(ret) <- paste0(rownames(ret), '/', signif(period, signif))
}
class(ret) <- c('gspline_matrix', class(ret))
ret
}
# Value and derivatives at 0 and all discontinuities at knots
Dmat <- function(knots,
degree,
periodic = FALSE,
signif = 3) {
dm <- max(degree)
cmat <- Xf(0, knots, dm, D = 0:dm, periodic = periodic)
n_knots <- length(knots)
for (i in seq_len(n_knots - 1)) {
k <- knots[i]
dmat <- Xf(k,
knots,
dm,
D = seq(0, dm),
FALSE,
periodic = periodic) -
Xf(k,
knots,
dm,
D = seq(0, dm),
TRUE,
periodic = periodic)
rownames(dmat) <-
paste0("C", seq(0, dm), '|', signif(k, signif))
if (periodic) {
rownames(dmat) <-
paste0(rownames(dmat), '/', signif(max(k), signif))
}
cmat <- rbind(cmat, dmat)
}
k <- knots[length(knots)]
if (periodic) {
dmat <-
Xf(0,
knots,
dm,
D = seq(0, dm),
FALSE,
periodic = periodic) -
Xf(k,
knots,
dm,
D = seq(0, dm),
TRUE,
periodic = periodic)
rownames(dmat) <-
paste("C",
seq(0, dm),
'|',
signif(k, signif),
'/',
signif(max(k), signif))
cmat <- rbind(cmat, dmat)
} else {
dmat <-
Xf(k,
knots,
dm,
D = seq(0, dm),
FALSE,
periodic = periodic) -
Xf(k,
knots,
dm,
D = seq(0, dm),
TRUE,
periodic = periodic)
rownames(dmat) <-
paste0("C", seq(0, dm), '|', signif(k, signif))
cmat <- rbind(cmat, dmat)
}
class(cmat) <- c('gspline_matrix', class(cmat))
cmat
}
# Parameters to constrain
#
# TODO: Change so degree can be a list
Pcon <- function(knots, degree, periodic) {
degree <- rep(degree, length.out = length(knots) + 1)
if (periodic) {
if (degree[length(degree)] != degree[1]) {
warning("For periodic splines, the degree of the last and first intervals should match")
}
knots <- knots[-length(knots)]
degree <- degree[-length(degree)]
}
dm <- max(degree)
cmat <- NULL
for (i in seq_along(degree)) {
di <- degree[i]
if (dm > di) {
dmat <- diag((length(knots) + 1) *
(dm + 1)) [(i - 1) * (dm + 1) + 1 + seq(di + 1, dm), , drop = FALSE]
rownames(dmat) <-
paste("I.", i, ".", seq(di + 1, dm), sep = '')
cmat <- rbind(cmat, dmat)
}
}
if (!is.null(cmat)) {
class(cmat) <- c('gspline_matrix', class(cmat))
}
cmat
}
#
# Parse knots, degree and smoothness to create constraint and estimation matrices:
# - Identify rows of Dmat that are used for constraints and rows used for estimates
if (periodic) {
if (min(knots) <= 0)
stop('For a periodic spline all knots must be positive')
}
if (any(knots != sort(knots)))
stop('Knots must be in ascending order')
degree <- rep(degree, length.out = length(knots) + 1)
max_degree <- max(degree)
smoothness <- rep(smoothness, length.out = length(knots))
if (!is.list(smoothness)) {
smoothness <-
lapply(smoothness, function(n)
if (n < 0)
- 1
else
0:n)
}
Dmat_smoothness_indices <- unlist(lapply(seq_along(smoothness),
function(i)
smoothness[[i]] + (max_degree + 1) * i + 1))
Dmat_smoothness_indices <-
Dmat_smoothness_indices[unlist(smoothness) > -1]
# Build constraints
constraint_mat <-
Dmat(knots, degree, periodic = periodic)[Dmat_smoothness_indices, , drop = FALSE]
constraint_mat <-
rbind(constraint_mat, Pcon(knots, degree, periodic = periodic))
if (!is.null(constraints))
constraint_mat <- rbind(constraint_mat, constraints)
if (!is.null(intercept)) {
constraint_mat <-
rbind(constraint_mat,
Xf(intercept, knots, max_degree, periodic = periodic))
Dmat_smoothness_indices <- c(1, Dmat_smoothness_indices)
}
cmat <- t(basis(t(constraint_mat)))
estimate_mat <-
rbind(estimates, Dmat(knots, degree, periodic = periodic)
[-Dmat_smoothness_indices, , drop = FALSE])
# the following is necessary in case some of the derivatives
# estimated at 0 have been coerced to 0 by other constraints
estimate_mat <- t(basis(t(rbind(
constraint_mat, estimate_mat
))))
estimate_mat <-
estimate_mat[-seq_len(nrow(constraint_mat)), , drop = FALSE]
# emat <- estimate_mat #FIXME not used JF
A <- rbind(constraint_mat, estimate_mat)
G <- solve(rbind(constraint_mat, estimate_mat),
diag(ncol(constraint_mat))[, -seq_len(nrow(constraint_mat)), drop = FALSE])
colnames(G) <- rownames(estimate_mat)
# create closure
fun <- function(x = NA, D = NULL, limit = 1) {
if (is.null(D)) {
ret <-
Xf(x, knots, max_degree, periodic = periodic) %*% G # model matrix
} else {
# Hypothesis matrix
D <- rep(D, length.out = length(x))
limit <- rep(limit, length.out = length(x))
left <- Xf(
x,
knots = knots,
degree = max_degree,
D = D,
periodic = periodic,
right = TRUE
)
right <- Xf(
x,
knots = knots,
degree = max_degree,
D = D,
periodic = periodic,
right = FALSE
)
cleft <- c(1, 0, -1)[match(limit, c(-1, 1, 0))]
cright <- c(0, 1, 1)[match(limit, c(-1, 1, 0))]
continuous <-
apply((right - left) %*% G, 1, function(x)
all(abs(x) < 10 * tol.continuity))
raw <- left * cleft + right * cright
ret <- raw %*% G
nam <- rownames(ret)
nam <- sub("^f", "g", nam)
nam_left <- sub("(/|$)", "-\\1", nam)
nam_right <- sub("(/|$)", "+\\1", nam)
nam_jump <- paste0(nam_right , "-", nam_left)
# If we're testing for a jump and the derivative is continuous we should show
# the identity of the jump tested.
#
# If we're testing for a derivative at a knot and the derivative is continuous, we should
# indicate that by removing the direction indicator
#
# So:
# If we're testing a jump we use the full name of the linear combination even
# if it happens to be at a point of continuity. Otherwise, used the directional notation only
# if the tested derivative is discontinuous
rownames(ret) <- ifelse(limit == 0, nam_jump,
ifelse(continuous, nam, ifelse(limit == 1, nam_right, nam_left)))
}
class(ret) <- c('gspline_matrix', class(ret))
ret
}
# Modify colnames of G matrix to show direction of limit
# if there is a discontinuity of derivatives at origin
if (!is.null(intercept)) {
dest <- fun(rep(intercept, max(degree) + 1),
D = 0:max(degree),
limit = -1)
tos <- grep('-$|-/', rownames(dest), value = TRUE)
froms <- sub('([0-9])-', '\\1', tos)
froms <- sub('(.*)', '^\\1$', froms)
froms <- sub('\\|', '\\\\|', froms)
for (i in seq_along(tos)) {
colnames(G) <- sub(froms[i], tos[i], colnames(G))
}
}
ffun <- if (!missing(x) && (stable || rescale)) {
if (rescale) {
if(!periodic) {
ran <- range(x)
a <- mean(ran)
b <- 20 / (ran[2] - ran[1]) # rescaled x goes from -10 to 10
} else { # if periodic
a <- 0
b <- 10 / knots[length(knots)]
}
fun_scaled <- gspline(
knots = (knots - a) * b,
degree = degree,
smoothness = smoothness,
periodic = periodic,
intercept = (intercept - a) * b,
stable = FALSE
)
}
function(x = NA, D = NULL, limit = 1) {
if (inwald() || !is.null(D)) {
ret <- fun(x, D, limit)
} else {
ret <- if (rescale)
fun_scaled((x - a) * b)
else
fun(x)
colnames(ret) <- NULL
# make basis orthogonal to intercept term GM:2019_10_11
if(ortho2intercept && !periodic)
ret <- sweep(ret, 2, colMeans(ret))
if (stable)
ret <- orth_basis(ret)
}
if (rescale)
attr(ret, 'rescaled') <- c(a = a, b = b)
if (stable)
attr(ret, 'orthonormalize') <- TRUE
if (ortho2intercept)
attr(ret, 'ortho2intercept') <- TRUE
class(ret) <- c('gspline_matrix', class(ret))
ret
}
} else
fun
class(ffun) <- c('gspline', class(ffun))
ffun
}
#' @method print gspline
#' @export
print.gspline <- function(x,
show = c('knots',
'degree',
'smoothness',
'G',
'constraint_mat',
'estimate_mat'),
...) {
cat('Spline function created by gspline\n')
envr <- environment(x)
ret <-
Filter(Negate(is.function), sapply(ls(environment(x)), get, environment(x)))
ret <-
lapply(ret, function(x)
if (is.numeric(x))
zapsmall(x)
else
x)
if (!is.null(show))
ret <- ret[show]
print(ret)
invisible(x)
}
#' @method print gspline_matrix
#' @export
print.gspline_matrix <- function(x, ...) {
header <- paste(if (!is.null(attr(x, "rescaled"))) "rescaled = TRUE",
if (!is.null(attr(x, "orthonormalize"))) "stable = TRUE",
if (!is.null(attr(x, "ortho2intercept"))) "ortho2intercept = TRUE",
sep="; ")
if (!is.null(header)) cat("\n", header, "\n")
xx <- zapsmall(x)
class(xx) <- "matrix"
attr(xx, "rescaled") <- NULL
attr(xx, "orthonormalize") <- NULL
attr(xx,"ortho2intercept") <- NULL
print(xx)
invisible(x)
}
#' @method knots gspline
#' @export
knots.gspline <- function(Fn, ...) {
environment(Fn)$knots
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.