R/gsp.R
In gmonette/spida2: Collection of tools developed for the Summer Programme in Data Analysis 2000-2012
Defines functions
sc
print.gsp
Documented in
sc
##
##
## General polynomial splines: August 8, 2008
##
## Modified GM 2019_02_22: added smoothness as list
##
## Revised June 16, 2010 to incorporate degree 0 polynomials
## and constraining evaluation of spline to 0 with 'intercept'
##
## Modified: June 15, 2013 - GM
## Added specified linear contraints to gsp and
## to sc expressed
## in terms of full polynomial parametrization
##
## Added 'PolyShift' function to specify constraints for
## a periodic spline. See example in the manual page.
##
## Added 'periodic' argument to gsp
##
##
#' General regression splines with variable degrees and ness, smoothing
#' splines
#'
#' These functions implement a general polynomial spline 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 \code{\link[spida2]{sc}} helps in the construction of linear hypothesis matrices
#' to estimate and test levels and derivatives of splines at arbitrary points
#' and the saltus of derivatives that have discontinuities at knots.
#'
#' 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 and smooth with a first derivative at each knot) with knots at 1 and 3
#' can be fitted 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 approach does not work for some more complex splines with different
#' degrees or different orders of continuity at the knots. An example is the
#' commonly used natural quadratic spline. A natural quadratic spline with
#' knots at -1,0 and 1 (where -1 and 1 are termed 'boundary knots') is linear
#' in the intervals (-Inf,-1) and (1,+Inf), and quadratic in the intervals
#' (-1,0) and (0,1). The spline is smooth of order 1 at each knot.
#'
#' Many
#' techniques for fitting splines generate a basis for the spline (columns of
#' the design matrix) that has good numerical properties but whose coefficients
#' do not have a simple interpretation.
#'
#' The \code{gsp} function makes it
#' easy to specify a spline with arbitrary degree in each interval and
#' arbitrary smoothness at each 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 each free saltus at each knot.
#'
#' The
#' \code{sc} function generates a matrix to estimate features of a fitted
#' spline that can be expressed as linear combinations of the spline
#' coefficients. Examples are various derivatives of the spline at any point,
#' left or right derivatives of different orders and the saltus in derivatives
#' at a knot.
#'
#' A disadvantage of \code{gsp} 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.
#'
#' \code{gsp}
#' generates a matrix of regressors for a spline with knots, degree of
#' polynomials in each interval and the degree of smoothness at each knot.
#' Typically, \code{gsp} is used to define a function that is then used in a
#' model equation. See the examples below.
#'
#' A function to fit a cubic spline with knots at 5 and 10 is generated with:
#'
#' Overview of utility functions:
#'
#' \code{Xmat = function( x, degree, D = 0, signif = 3)} design/estimation matrix for f[D](x) where f(x) is
#' polynomial of degree degree.
#'
#' \code{Xf = function( x, knots, degree = 3, D = 0, right = TRUE , signif = 3)}
#' uses
#' Xmat to form 'full' matrix with blocks determined by knots intervals
#'
#' \code{Cmat = function( knots, degree, smoothness, intercept = 0, signif = 3)}
#' linear constraints
#'
#' \code{Emat = function( knots, degree, smoothness , intercept = FALSE, signif = 3)}
#' estimates - not necessarily a basis
#'
#' \code{basis = function( X , coef = FALSE )} selects linear independent subset of
#' columns of X
#'
#' \code{spline.T = function( knots, degree, smoothness, intercept = 0, signif = 3 )}
#' full transformation of Xf to spline basis and constraints
#'
#' \code{spline.E = function( knots, degree, smoothness, intercept = 0, signif = 3 )}
#' transformation for spline basis (first r columns of spline.T)
#'
#' @aliases gsp smsp Xf Xmat qs lsp Cmat PolyShift
#' @param x value(s) where spline is evaluated
#' @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. If the spline should evaluate to
#' 0 in the interval, use the \code{intercept} argument to specify some value
#' in the interval at which the spline must evaluate to 0.
#' @param smoothness vector with the degree of smoothness at each knot (0 =
#' continuity, 1 = smoothness with continuous first derivative, 2 = continuous
#' second derivative, etc. The value -1 allows a discontinuity at the knot.
#' A scalar is recycled so its length equals the
#' number of knots. Alternatively,
#' a list of length equal to the number of knots. Each element of the list is a
#' vector of the orders of derivatives which are required to be smooth. THis allows
#' non-sequential constraints, e.g. to have the same first and second derivative
#' on either side of a knot but a possible discontinuity and change in
#' higher-order derivatives, the vector would be c(1,2). Note that if a list is used,
#' all elements must provide all desired constraints. That is the list argument corresponding
#' to `smoothness = c(1,2,-1)` is `smoothness=list(0:1, 0:2, -1)`. The default is
#' 1 less than the minimum of adjoining degrees at
#'
#'
#'
#'
#'
#' @param intercept value(s) of x at which the spline has value 0, i.e. the
#' value(s) of x for which yhat is estimated by the intercept term in the
#' model. The default is 0. If NULL, the spline is not constrained to evaluate
#' to 0 for any x.
#' @param periodic if TRUE generates a period spline on the base interval
#' [0,max(knots)]. A constraint is generated so that the coefficients generate
#' the same values to the right of max(knots) as they do to the right of 0.
#' Note that all knots should be strictly positive.
#' @param lin provides a matrix specifying additional linear contraints on the
#' 'full' parametrization consisting of blocks of polynomials of degree equal
#' to max(degree) in each of the length(knots)+1 intervals of the spline. See
#' below for examples of a spline that is 0 outside of its boundary knots.
#' @param signif number of significant digits used to label coefficients
#' @param exclude number of leading columns to drop from spline matrix: 0:
#' excludes the intercept column, 1: excludes the linear term as well. Terms
#' that are excluded from the spline matrix can be modeled explicitly.
#' @param a a 'periodic' knot at which a spline repeats itself
#' @param n the maximal 'order', i.e. maximal degree + 1, of a periodic spline
#' @return \code{gsp} returns a matrix generating a spline. \code{cs},
#' \code{qs} and \code{lsp} return matrices generating cubic, quadratic and
#' linear splines respectively.
#'
#' \code{smsp}, whose code is adapted from function in the package
#' \code{lmeSplines}, generates a smoothing spline to be used in the random
#' effects portion of a call to \code{lme}.
#'
#' @section Warning: The variables generated by \code{gsp} are designed so the
#' coefficients are interpretable as changes in derivatives at knots. The
#' resulting matrix is not designed to have optimal numerical properties.
#'
#' The intermediate matrices generated by \code{gsp} will contain \code{x}
#' raised to the power equal to the highest degree in \code{degree}. The values
#' of \code{x} should be scaled to avoid excessively large values in the spline
#' matrix and ensuing numerical problems.
#' @author Monette, G. \email{georges@@yorku.ca}
#' @seealso \code{\link[spida2]{wald}} \code{\link[spida2]{sc}}
#' @references %% ~put references to the literature/web site here ~
#' @keywords ~kwd1 ~kwd2
#' @examples
#' ## Fitting a quadratic spline
#' 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
#' sp <- function(x) gsp( x, knots = c(10,25,40), degree = c(1,2,2,1),
#' smoothness = c(1,1,1))
#'
#' fit <- lm( y ~ sp(age)*G, simd)
#'
#' require(lattice)
#' xyplot( predict(fit) ~ age , simd, groups = G,type = 'l')
#' summary(fit)
#'
#' ## Linear hypotheses
#' L <- list( "Overall test of slopes at 20" = rbind(
#' "Female slope at age 20" = c( F20 <- cbind( 0 , sc( sp, 20, D = 1), 0 , 0 * sc( sp, 20, D = 1))),
#' "Male slope at age 20" = c( M20 <- cbind( 0 , sc( sp, 20, D = 1), 0 , 1 * sc( sp, 20, D = 1))),
#' "Difference" = c(M20 - F20))
#' )
#' wald( fit, L)
#'
#' ## Right and left second derivatives at knots and saltus
#' L <- list( "Second derivatives and saltus for female curve at knot at 25" =
#' cbind( 0, sc( sp, c(25,25,25), D = 2, type =c(0,1,2)), 0,0,0,0))
#' L
#' wald( fit, L )
#'
#' # Quadratic natural spline
#' 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)))
#'
#' sp <- function(x) gsp( x, knots = c(10,25,40),
#' degree = c(1,2,2,1),
#' smoothness = c(1,1,1))
#'
#' fit <- lm( y ~ sp(age)*G, simd)
#' xyplot( predict(fit) ~ age , simd, groups = G, type = "l")
#' summary(fit) # convenient display
#'
#' # output:
#' #
#' # Call:
#' # lm(formula = y ~ sp(age) * G, data = simd)
#' #
#' # Residuals:
#' # Min 1Q Median 3Q Max
#' # -3.1111 -0.9429 0.0774 0.6633 3.6171
#' #
#' # Coefficients:
#' # Estimate Std. Error t value Pr(>|t|)
#' # (Intercept) 0.522681 0.623071 0.839 0.404
#' # sp(age)D1(0) -0.060361 0.056802 -1.063 0.291
#' # sp(age)C(10).2 0.009013 0.007116 1.267 0.208
#' # sp(age)C(25).2 -0.015459 0.013263 -1.166 0.247
#' # Gmale -0.092681 0.881156 -0.105 0.916
#' # sp(age)D1(0):Gmale 0.042878 0.080331 0.534 0.595
#' # sp(age)C(10).2:Gmale -0.006703 0.010063 -0.666 0.507
#' # sp(age)C(25).2:Gmale 0.007650 0.018757 0.408 0.684
#' #
#' # Residual standard error: 1.261 on 92 degrees of freedom
#' # Multiple R-squared: 0.05212, Adjusted R-squared: -0.02
#' # F-statistic: 0.7227 on 7 and 92 DF, p-value: 0.653
#' # ---
#'
#' L0 <- list(
#' "hat" = rbind(
#' "females at age=20" = c( 1, sc(sp,20), 0, 0 * sc(sp,20)),
#' "males at age=20" = c( 1, sc(sp,20), 1, 1* sc(sp,20))),
#' "male-female" = rbind(
#' "at 20" = c( 0 , 0*sc(sp,20), 1, 1*sc(sp,20))))
#' wald(fit, L0 )
#'
#' L1 <- list(
#' "D(yhat)/D(age)" =
#' rbind( "female at age = 25" = c(0, sc(sp,25,1), 0, 0*sc(sp,25,1)),
#' "male at x = 25" = c(0, sc(sp,25,1), 0, 1*sc(sp,25,1))))
#' wald( fit, L1)
#'
#' # A periodic spline function can be generated by forcing the coefficients
#' # beyond a periodic knot to repeat the pattern in a previous inteval. For
#' # example a periodic spline of period 1 can be created as follows:
#'
#' x <- seq(-3,3,.01)
#' y <- rnorm(x)
#'
#' per.sp <- function(x) gsp( x %% 1, knots = 1,
#' degree = c(3, 3),
#' smoothness = 1,
#' lin = cbind( diag(4), - PolyShift(1,4)))
#' fit <- lm( y ~ per.sp(x))
#' summary(fit)
#' xyplot( y + predict(fit) ~ x ,
#' panel = panel.superpose.2,
#' type = c("p","l"))
#' xyplot( predict(fit) ~ x , type = 'l')
#' #' # A periodic spline with additional knots can be created as shown below. Note
#' # that constraint matrix 'lin' expresses constraints for the 'full' polynomial
#' # parametrization, i.e. polynomials of degree 3 (thus of order 4 when
#' # including the constant term) in each of the 4 intervals. The constraint
#' # given in 'lin' forces the coefficients beyond the periodic knot at x = 1, to
#' # repeat the polynomial in the interval just to the right of x = 0.
#'
#' x <- seq(-3,3,.01)
#' y <- rnorm(x)
#'
#' per.sp <- function(x) gsp( x %% 1, knots = c(.333,.666,1),
#' degree = 3,
#' smoothness = 1,
#' lin = cbind( diag(4),0*diag(4),0*diag(4), - PolyShift(1,4)))
#' fit <- lm( y ~ per.sp(x))
#' summary(fit)
#' xyplot( y + predict(fit) ~ x ,
#' panel = panel.superpose.2,
#' type = c("p","l"))
#' xyplot( predict(fit) ~ x , type = 'l')
#'
#' # Cubic spline:
#'
#' sp <- function( x ) gsp( x, c(5,10), c(3,3,3), 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.
#'
#' # Natural cubic 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 <- function( x ) gsp( x, c(0,5,10,15), c(1,3,3,3,1), c(1,2,2,1))
#'
#' # Quadratic and linear splines, respectively:
#'
#' sp.quad <- function( x ) gsp( x, c(5,10), c(2,2,2), c(1,1))
#' sp.lin <- function( x ) gsp( x, c(5,10), c(1,1,1), c(0,0))
#'
#' # When the same degree is used for all
#' # intervals and knots, it suffices to give it once:
#'
#' sp.quad <- function( x ) gsp( x, c(5,10), 2, 1)
#' sp.lin <- function( x ) gsp( x, c(5,10), 1, 0)
#'
#' # 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 <- function( x ) gsp( x, c(5,10), c(3,3,3), 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 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.
#'
#' # The naming of coefficients should allow them to be easily interpreted. For
#' # example:
#'
#' zapsmall( gsp ( 0:10, c(3, 7) , c(2,3,2), c(1,1)))
#'
#' ## Output:
#' ##
#' ## D1(0) D2(0) C(3).2 C(3).3 C(7).2
#' ## f(0) 0 0.0 0.0 0.00000 0.0
#' ## f(1) 1 0.5 0.0 0.00000 0.0
#' ## f(2) 2 2.0 0.0 0.00000 0.0
#' ## f(3) 3 4.5 0.0 0.00000 0.0
#' ## f(4) 4 8.0 0.5 0.16667 0.0
#' ## f(5) 5 12.5 2.0 1.33333 0.0
#' ## f(6) 6 18.0 4.5 4.50000 0.0
#' ## f(7) 7 24.5 8.0 10.66667 0.0
#' ## f(8) 8 32.0 12.5 20.66667 0.5
#' ## f(9) 9 40.5 18.0 34.66667 2.0
#' ## f(10) 10 50.0 24.5 52.66667 4.5
#'
#' # 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 <- function(x) gsp ( x, c(3, 7) , c(2,3,2), c(1,1))
#' zd <- data.frame( x = seq(0,10, .5),
#' y = seq(0,10,.5)^2 + rnorm(21))
#' fit <- lm( y ~ sp( x ), zd)
#' summary(fit)
#' Ls <-cbind( 0, sc( sp, c(1,2,3,3,3,5,7,7,7,8),
#' D = 3, type = c(0,0,0,1,2,0,0,1,2,0)))
#' zapsmall(Ls )
#' ##
#' ## Output:
#' ## D1(0) D2(0) C(3).2 C(3).3 C(7).2
#' ## D3(1) 0 0 0 0 0 0
#' ## D3(2) 0 0 0 0 0 0
#' ## D3(3-) 0 0 0 0 0 0
#' ## D3(3+) 0 0 0 0 1 0
#' ## D3(3+)-D3(3-) 0 0 0 0 1 0
#' ## D3(5) 0 0 0 0 1 0
#' ## D3(7-) 0 0 0 0 1 0
#' ## D3(7+) 0 0 0 0 0 0
#' ## D3(7+)-D3(7-) 0 0 0 0 -1 0
#' ## D3(8) 0 0 0 0 0 0
#'
#' wald( fit, list( 'third derivatives' = Ls))
#'
#' ## Output:
#' ##
#' ## numDF denDF F.value p.value
#' ## third derivatives 1 15 2.184522 0.16009
#' ##
#' ## ? Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
#' ## D3(1) 0.000000 0.000000 15 2.408694 0.02932 0.000000 0.000000
#' ## D3(2) 0.000000 0.000000 15 2.408694 0.02932 0.000000 0.000000
#' ## D3(3-) 0.000000 0.000000 15 2.408694 0.02932 0.000000 0.000000
#' ## D3(3+) -0.535959 0.362621 15 -1.478013 0.16009 -1.308867 0.236950
#' ## D3(3+)-D3(3-) -0.535959 0.362621 15 -1.478013 0.16009 -1.308867 0.236950
#' ## D3(5) -0.535959 0.362621 15 -1.478013 0.16009 -1.308867 0.236950
#' ## D3(7-) -0.535959 0.362621 15 -1.478013 0.16009 -1.308867 0.236950
#' ## D3(7+) 0.000000 0.000000 Inf NaN NaN 0.000000 0.000000
#' ## D3(7+)-D3(7-) 0.535959 0.362621 15 1.478013 0.16009 -0.236950 1.308867
#' ## D3(8) 0.000000 0.000000 Inf NaN NaN 0.000000 0.000000
#'
#' # Note that some coefficients that are 0 by design may lead to invalid degrees
#' # of freedom and t-values.
#'
#' # sc generates a portion of a hypothesis matrix for the coefficients of
#' # a general spline constructed with gsp With:
#' #
#' # sc( sp, x, D, type ):
#' # sp is the spline function for which coefficients
#' # are required,
#' # x is the vector of values at which the spline is evaluated,
#' # D is the order of derivatives: 0 = value of spline,
#' # 1 = first derivative, etc.,
#' # type at knots: 0 limit from the left, 1 from the
#' # right, 2 is saltus (i.e. jump from left to right)
#'
#' # Warning: sc will not work correctly if the function defining the
#' # spline transforms the variable, e.g. sp <- function(x) gsp( x/100, knot=2 )
#'
#' # Example:
#'
#' 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)))
#'
#' sp <- function(x) gsp( x, knots = c(10,25,40),
#' degree = c(1,2,2,1),
#' smoothness = c(1,1,1))
#'
#' fit <- lm( y ~ sp(age)*G, simd)
#' xyplot( predict(fit) ~ age , simd, groups = G,type = "l")
#' summary(fit)
#'
#' L0 <- list(
#' "hat" = rbind(
#' "females at age=20" = c( 1, sc(sp,20), 0, 0 * sc(sp,20)),
#' "males at age=20" = c( 1, sc(sp,20), 1, 1* sc(sp,20))),
#' "male-female" = rbind(
#' "at 20" = c( 0 , 0*sc(sp,20), 1, 1*sc(sp,20))))
#' wald(fit, L0 )
#'
#' L1 <- list(
#' "D(yhat)/D(age)" =
#' rbind( "female at age = 25" = c(0, sc(sp,25,1), 0, 0*sc(sp,25,1)),
#' "male at x = 25" = c(0, sc(sp,25,1), 0, 1*sc(sp,25,1))))
#' wald( fit, L1)
#'
#' # A periodic spline function can be generated by forcing the coefficients
#' # beyond a periodic knot to repeat the pattern in a previous inteval. For
#' # example a periodic spline of period 1 can be created as follows:
#'
#' x <- seq(-3,3,.01)
#' y <- rnorm(x)
#'
#' per.sp <- function(x) gsp( x %% 1,
#' knots = 1,
#' degree = c(3, 3),
#' smoothness = 1,
#' lin = cbind( diag(4), - PolyShift(1,4)))
#' fit <- lm( y ~ per.sp(x))
#' summary(fit)
#' xyplot( y + predict(fit) ~ x ,
#' panel = panel.superpose.2,
#' type = c("p","l"))
#' xyplot( predict(fit) ~ x , type = 'l')
#'
#' # A periodic spline with additional knots can be created as shown below. Note
#' # that constraint matrix 'lin' expresses constraints for the 'full' polynomial
#' # parametrization, i.e. polynomials of degree 3 (thus of order 4 when
#' # including the constant term) in each of the 4 intervals. The constraint
#' # given in 'lin' forces the coefficients beyond the periodic knot at x = 1, to
#' # repeat the polynomial in the interval just to the right of x = 0.
#'
#' x <- seq(-3,3,.01)
#' y <- rnorm(x)
#'
#' per.sp <- function(x) gsp( x %% 1, knots = c(.333,.666,1),
#' degree = 3,
#' smoothness = 1,
#' lin = cbind( diag(4),0*diag(4),0*diag(4), - PolyShift(1,4)))
#' fit <- lm( y ~ per.sp(x))
#' summary(fit)
#' xyplot( y + predict(fit) ~ x ,
#' panel = panel.superpose.2,type = c("p","l"))
#' xyplot( predict(fit) ~ x , type = 'l')
#'
#'
#' ## Smoothing splines
#'
#' library(nlme)
#' library(MASS)
#' data(Spruce)
#' Spruce$all <- 1
#' range( Spruce$days)
#' sp <- function(x) smsp ( x, seq( 150, 675, 5))
#' spruce.fit1 <- lme(logSize ~ days, data=Spruce,
#' random=list(all= pdIdent(~sp(days) -1),
#' plot=~1, Tree=~1))
#' summary(spruce.fit1)
#' pred <- expand.grid( days = seq( 160, 670, 10), all = 1)
#' pred$logSize <- predict( spruce.fit1, newdata = pred, level = 1)
#' library( lattice )
#' xyplot( logSize ~ days, pred, type = 'l')
#'
#'
#' @export
gsp <- function (x, knots, degree = 3,
smoothness = pmax(pmin(degree[-1], degree[-length(degree)]) - 1, -1),
lin = NULL, periodic = FALSE, intercept = 0, signif = 3)
{
if (periodic) {
maxd <- max(degree)
maxk <- max(knots)
mid <- matrix(0, maxd + 1, (maxd + 1) * (length(knots) - 1))
const.per <- do.call(cbind, list(diag(maxd + 1), mid,
-PolyShift(maxk, maxd + 1)))
lin <- rbind(lin, const.per)
}
degree = rep(degree, length = length(knots) + 1)
smoothness = rep(smoothness, length = length(knots))
spline.attr <- list(knots = knots, degree = degree, smoothness = smoothness,
lin = lin, intercept = intercept, signif = signif)
if (is.null(x))
return(spline.attr)
if (periodic)
x <- x%%maxk
ret = Xf(x, knots, max(degree), signif = signif) %*% spline.E(knots,
degree, smoothness, lin = lin, intercept = intercept, signif = signif)
attr(ret, "spline.attr") <- spline.attr
class(ret) <- "gsp"
ret
}
print.gsp <- function(x, strip.attributes=TRUE, ...){
nms <- colnames(x)
ncol <- ncol(x)
xx <- x
if (strip.attributes) {
attributes(xx) <- NULL
xx <- matrix(xx, ncol=ncol)
colnames(xx) <- nms
}
else (xx <- unclass(xx))
print(zapsmall(xx), ...)
invisible(x)
}
#' Hypothesis matrix for general regression splines
#'
#' This function helps to build hypothesis matrices for general splines. See examples
#' in \code{\link[spida2]{gsp}}.
#'
#' @param sp a spline function generated by \code{\link[spida2]{gsp}}
#' @param x points where spline is tested
#' @param D (default 0) order of derivative to test, 0 is the value of the spline, 1 is the first derivative, etc.
#' @param type (default 1) if there is a discontinuity at a knot, type specifies whether to estimate the limit from the left (0), the limit from the right (1) or the saltus -- limit from the right minus the limit from left (2)
#' @seealso \code{\link[spida2]{wald}} \code{\link[spida2]{gsp}}
#' @export
sc <-
function( sp, x , D = 0, type = 1) {
a = sp(NULL)
D = rep(D, length.out = length(x))
type = rep ( type, length.out = length(x))
left = Xf( x, knots = a$knots, degree = max( a$degree), D = D, right = T)
right = Xf( x, knots = a$knots, degree = max( a$degree), D = D, right = F)
cleft = c(1,0,-1) [ match( type, c(0,1,2)) ]
cright = c(0,1,1) [ match( type, c(0,1,2)) ]
raw = left * cleft + right * cright
nam = rownames( raw )
nam = sub("^f","g", nam)
nam0 = sub( "\\)","-)", nam)
nam1 = sub( "\\)","+)", nam)
nam2 = paste( nam1 ,"-",nam0,sep = '')
rownames( raw ) =
ifelse( match( x , a$knots, 0) > 0,
cbind( nam0,nam1,nam2) [ cbind( seq_along(type), type+1)],
ifelse( type != 2, nam, '0'))
mod = raw %*% spline.E(
a$knots, a$degree, a$smoothness,
lin = a$lin,
intercept = a$intercept,
signif = a$signif)
mod
}
gmonette/spida2 documentation built on Aug. 20, 2023, 7:21 p.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 sc print.gsp
Documented in sc
##
##
## General polynomial splines: August 8, 2008
##
## Modified GM 2019_02_22: added smoothness as list
##
## Revised June 16, 2010 to incorporate degree 0 polynomials
## and constraining evaluation of spline to 0 with 'intercept'
##
## Modified: June 15, 2013 - GM
## Added specified linear contraints to gsp and
## to sc expressed
## in terms of full polynomial parametrization
##
## Added 'PolyShift' function to specify constraints for
## a periodic spline. See example in the manual page.
##
## Added 'periodic' argument to gsp
##
##
#' General regression splines with variable degrees and ness, smoothing
#' splines
#'
#' These functions implement a general polynomial spline 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 \code{\link[spida2]{sc}} helps in the construction of linear hypothesis matrices
#' to estimate and test levels and derivatives of splines at arbitrary points
#' and the saltus of derivatives that have discontinuities at knots.
#'
#' 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 and smooth with a first derivative at each knot) with knots at 1 and 3
#' can be fitted 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 approach does not work for some more complex splines with different
#' degrees or different orders of continuity at the knots. An example is the
#' commonly used natural quadratic spline. A natural quadratic spline with
#' knots at -1,0 and 1 (where -1 and 1 are termed 'boundary knots') is linear
#' in the intervals (-Inf,-1) and (1,+Inf), and quadratic in the intervals
#' (-1,0) and (0,1). The spline is smooth of order 1 at each knot.
#'
#' Many
#' techniques for fitting splines generate a basis for the spline (columns of
#' the design matrix) that has good numerical properties but whose coefficients
#' do not have a simple interpretation.
#'
#' The \code{gsp} function makes it
#' easy to specify a spline with arbitrary degree in each interval and
#' arbitrary smoothness at each 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 each free saltus at each knot.
#'
#' The
#' \code{sc} function generates a matrix to estimate features of a fitted
#' spline that can be expressed as linear combinations of the spline
#' coefficients. Examples are various derivatives of the spline at any point,
#' left or right derivatives of different orders and the saltus in derivatives
#' at a knot.
#'
#' A disadvantage of \code{gsp} 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.
#'
#' \code{gsp}
#' generates a matrix of regressors for a spline with knots, degree of
#' polynomials in each interval and the degree of smoothness at each knot.
#' Typically, \code{gsp} is used to define a function that is then used in a
#' model equation. See the examples below.
#'
#' A function to fit a cubic spline with knots at 5 and 10 is generated with:
#'
#' Overview of utility functions:
#'
#' \code{Xmat = function( x, degree, D = 0, signif = 3)} design/estimation matrix for f[D](x) where f(x) is
#' polynomial of degree degree.
#'
#' \code{Xf = function( x, knots, degree = 3, D = 0, right = TRUE , signif = 3)}
#' uses
#' Xmat to form 'full' matrix with blocks determined by knots intervals
#'
#' \code{Cmat = function( knots, degree, smoothness, intercept = 0, signif = 3)}
#' linear constraints
#'
#' \code{Emat = function( knots, degree, smoothness , intercept = FALSE, signif = 3)}
#' estimates - not necessarily a basis
#'
#' \code{basis = function( X , coef = FALSE )} selects linear independent subset of
#' columns of X
#'
#' \code{spline.T = function( knots, degree, smoothness, intercept = 0, signif = 3 )}
#' full transformation of Xf to spline basis and constraints
#'
#' \code{spline.E = function( knots, degree, smoothness, intercept = 0, signif = 3 )}
#' transformation for spline basis (first r columns of spline.T)
#'
#' @aliases gsp smsp Xf Xmat qs lsp Cmat PolyShift
#' @param x value(s) where spline is evaluated
#' @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. If the spline should evaluate to
#' 0 in the interval, use the \code{intercept} argument to specify some value
#' in the interval at which the spline must evaluate to 0.
#' @param smoothness vector with the degree of smoothness at each knot (0 =
#' continuity, 1 = smoothness with continuous first derivative, 2 = continuous
#' second derivative, etc. The value -1 allows a discontinuity at the knot.
#' A scalar is recycled so its length equals the
#' number of knots. Alternatively,
#' a list of length equal to the number of knots. Each element of the list is a
#' vector of the orders of derivatives which are required to be smooth. THis allows
#' non-sequential constraints, e.g. to have the same first and second derivative
#' on either side of a knot but a possible discontinuity and change in
#' higher-order derivatives, the vector would be c(1,2). Note that if a list is used,
#' all elements must provide all desired constraints. That is the list argument corresponding
#' to `smoothness = c(1,2,-1)` is `smoothness=list(0:1, 0:2, -1)`. The default is
#' 1 less than the minimum of adjoining degrees at
#'
#'
#'
#'
#'
#' @param intercept value(s) of x at which the spline has value 0, i.e. the
#' value(s) of x for which yhat is estimated by the intercept term in the
#' model. The default is 0. If NULL, the spline is not constrained to evaluate
#' to 0 for any x.
#' @param periodic if TRUE generates a period spline on the base interval
#' [0,max(knots)]. A constraint is generated so that the coefficients generate
#' the same values to the right of max(knots) as they do to the right of 0.
#' Note that all knots should be strictly positive.
#' @param lin provides a matrix specifying additional linear contraints on the
#' 'full' parametrization consisting of blocks of polynomials of degree equal
#' to max(degree) in each of the length(knots)+1 intervals of the spline. See
#' below for examples of a spline that is 0 outside of its boundary knots.
#' @param signif number of significant digits used to label coefficients
#' @param exclude number of leading columns to drop from spline matrix: 0:
#' excludes the intercept column, 1: excludes the linear term as well. Terms
#' that are excluded from the spline matrix can be modeled explicitly.
#' @param a a 'periodic' knot at which a spline repeats itself
#' @param n the maximal 'order', i.e. maximal degree + 1, of a periodic spline
#' @return \code{gsp} returns a matrix generating a spline. \code{cs},
#' \code{qs} and \code{lsp} return matrices generating cubic, quadratic and
#' linear splines respectively.
#'
#' \code{smsp}, whose code is adapted from function in the package
#' \code{lmeSplines}, generates a smoothing spline to be used in the random
#' effects portion of a call to \code{lme}.
#'
#' @section Warning: The variables generated by \code{gsp} are designed so the
#' coefficients are interpretable as changes in derivatives at knots. The
#' resulting matrix is not designed to have optimal numerical properties.
#'
#' The intermediate matrices generated by \code{gsp} will contain \code{x}
#' raised to the power equal to the highest degree in \code{degree}. The values
#' of \code{x} should be scaled to avoid excessively large values in the spline
#' matrix and ensuing numerical problems.
#' @author Monette, G. \email{georges@@yorku.ca}
#' @seealso \code{\link[spida2]{wald}} \code{\link[spida2]{sc}}
#' @references %% ~put references to the literature/web site here ~
#' @keywords ~kwd1 ~kwd2
#' @examples
#' ## Fitting a quadratic spline
#' 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
#' sp <- function(x) gsp( x, knots = c(10,25,40), degree = c(1,2,2,1),
#' smoothness = c(1,1,1))
#'
#' fit <- lm( y ~ sp(age)*G, simd)
#'
#' require(lattice)
#' xyplot( predict(fit) ~ age , simd, groups = G,type = 'l')
#' summary(fit)
#'
#' ## Linear hypotheses
#' L <- list( "Overall test of slopes at 20" = rbind(
#' "Female slope at age 20" = c( F20 <- cbind( 0 , sc( sp, 20, D = 1), 0 , 0 * sc( sp, 20, D = 1))),
#' "Male slope at age 20" = c( M20 <- cbind( 0 , sc( sp, 20, D = 1), 0 , 1 * sc( sp, 20, D = 1))),
#' "Difference" = c(M20 - F20))
#' )
#' wald( fit, L)
#'
#' ## Right and left second derivatives at knots and saltus
#' L <- list( "Second derivatives and saltus for female curve at knot at 25" =
#' cbind( 0, sc( sp, c(25,25,25), D = 2, type =c(0,1,2)), 0,0,0,0))
#' L
#' wald( fit, L )
#'
#' # Quadratic natural spline
#' 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)))
#'
#' sp <- function(x) gsp( x, knots = c(10,25,40),
#' degree = c(1,2,2,1),
#' smoothness = c(1,1,1))
#'
#' fit <- lm( y ~ sp(age)*G, simd)
#' xyplot( predict(fit) ~ age , simd, groups = G, type = "l")
#' summary(fit) # convenient display
#'
#' # output:
#' #
#' # Call:
#' # lm(formula = y ~ sp(age) * G, data = simd)
#' #
#' # Residuals:
#' # Min 1Q Median 3Q Max
#' # -3.1111 -0.9429 0.0774 0.6633 3.6171
#' #
#' # Coefficients:
#' # Estimate Std. Error t value Pr(>|t|)
#' # (Intercept) 0.522681 0.623071 0.839 0.404
#' # sp(age)D1(0) -0.060361 0.056802 -1.063 0.291
#' # sp(age)C(10).2 0.009013 0.007116 1.267 0.208
#' # sp(age)C(25).2 -0.015459 0.013263 -1.166 0.247
#' # Gmale -0.092681 0.881156 -0.105 0.916
#' # sp(age)D1(0):Gmale 0.042878 0.080331 0.534 0.595
#' # sp(age)C(10).2:Gmale -0.006703 0.010063 -0.666 0.507
#' # sp(age)C(25).2:Gmale 0.007650 0.018757 0.408 0.684
#' #
#' # Residual standard error: 1.261 on 92 degrees of freedom
#' # Multiple R-squared: 0.05212, Adjusted R-squared: -0.02
#' # F-statistic: 0.7227 on 7 and 92 DF, p-value: 0.653
#' # ---
#'
#' L0 <- list(
#' "hat" = rbind(
#' "females at age=20" = c( 1, sc(sp,20), 0, 0 * sc(sp,20)),
#' "males at age=20" = c( 1, sc(sp,20), 1, 1* sc(sp,20))),
#' "male-female" = rbind(
#' "at 20" = c( 0 , 0*sc(sp,20), 1, 1*sc(sp,20))))
#' wald(fit, L0 )
#'
#' L1 <- list(
#' "D(yhat)/D(age)" =
#' rbind( "female at age = 25" = c(0, sc(sp,25,1), 0, 0*sc(sp,25,1)),
#' "male at x = 25" = c(0, sc(sp,25,1), 0, 1*sc(sp,25,1))))
#' wald( fit, L1)
#'
#' # A periodic spline function can be generated by forcing the coefficients
#' # beyond a periodic knot to repeat the pattern in a previous inteval. For
#' # example a periodic spline of period 1 can be created as follows:
#'
#' x <- seq(-3,3,.01)
#' y <- rnorm(x)
#'
#' per.sp <- function(x) gsp( x %% 1, knots = 1,
#' degree = c(3, 3),
#' smoothness = 1,
#' lin = cbind( diag(4), - PolyShift(1,4)))
#' fit <- lm( y ~ per.sp(x))
#' summary(fit)
#' xyplot( y + predict(fit) ~ x ,
#' panel = panel.superpose.2,
#' type = c("p","l"))
#' xyplot( predict(fit) ~ x , type = 'l')
#' #' # A periodic spline with additional knots can be created as shown below. Note
#' # that constraint matrix 'lin' expresses constraints for the 'full' polynomial
#' # parametrization, i.e. polynomials of degree 3 (thus of order 4 when
#' # including the constant term) in each of the 4 intervals. The constraint
#' # given in 'lin' forces the coefficients beyond the periodic knot at x = 1, to
#' # repeat the polynomial in the interval just to the right of x = 0.
#'
#' x <- seq(-3,3,.01)
#' y <- rnorm(x)
#'
#' per.sp <- function(x) gsp( x %% 1, knots = c(.333,.666,1),
#' degree = 3,
#' smoothness = 1,
#' lin = cbind( diag(4),0*diag(4),0*diag(4), - PolyShift(1,4)))
#' fit <- lm( y ~ per.sp(x))
#' summary(fit)
#' xyplot( y + predict(fit) ~ x ,
#' panel = panel.superpose.2,
#' type = c("p","l"))
#' xyplot( predict(fit) ~ x , type = 'l')
#'
#' # Cubic spline:
#'
#' sp <- function( x ) gsp( x, c(5,10), c(3,3,3), 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.
#'
#' # Natural cubic 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 <- function( x ) gsp( x, c(0,5,10,15), c(1,3,3,3,1), c(1,2,2,1))
#'
#' # Quadratic and linear splines, respectively:
#'
#' sp.quad <- function( x ) gsp( x, c(5,10), c(2,2,2), c(1,1))
#' sp.lin <- function( x ) gsp( x, c(5,10), c(1,1,1), c(0,0))
#'
#' # When the same degree is used for all
#' # intervals and knots, it suffices to give it once:
#'
#' sp.quad <- function( x ) gsp( x, c(5,10), 2, 1)
#' sp.lin <- function( x ) gsp( x, c(5,10), 1, 0)
#'
#' # 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 <- function( x ) gsp( x, c(5,10), c(3,3,3), 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 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.
#'
#' # The naming of coefficients should allow them to be easily interpreted. For
#' # example:
#'
#' zapsmall( gsp ( 0:10, c(3, 7) , c(2,3,2), c(1,1)))
#'
#' ## Output:
#' ##
#' ## D1(0) D2(0) C(3).2 C(3).3 C(7).2
#' ## f(0) 0 0.0 0.0 0.00000 0.0
#' ## f(1) 1 0.5 0.0 0.00000 0.0
#' ## f(2) 2 2.0 0.0 0.00000 0.0
#' ## f(3) 3 4.5 0.0 0.00000 0.0
#' ## f(4) 4 8.0 0.5 0.16667 0.0
#' ## f(5) 5 12.5 2.0 1.33333 0.0
#' ## f(6) 6 18.0 4.5 4.50000 0.0
#' ## f(7) 7 24.5 8.0 10.66667 0.0
#' ## f(8) 8 32.0 12.5 20.66667 0.5
#' ## f(9) 9 40.5 18.0 34.66667 2.0
#' ## f(10) 10 50.0 24.5 52.66667 4.5
#'
#' # 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 <- function(x) gsp ( x, c(3, 7) , c(2,3,2), c(1,1))
#' zd <- data.frame( x = seq(0,10, .5),
#' y = seq(0,10,.5)^2 + rnorm(21))
#' fit <- lm( y ~ sp( x ), zd)
#' summary(fit)
#' Ls <-cbind( 0, sc( sp, c(1,2,3,3,3,5,7,7,7,8),
#' D = 3, type = c(0,0,0,1,2,0,0,1,2,0)))
#' zapsmall(Ls )
#' ##
#' ## Output:
#' ## D1(0) D2(0) C(3).2 C(3).3 C(7).2
#' ## D3(1) 0 0 0 0 0 0
#' ## D3(2) 0 0 0 0 0 0
#' ## D3(3-) 0 0 0 0 0 0
#' ## D3(3+) 0 0 0 0 1 0
#' ## D3(3+)-D3(3-) 0 0 0 0 1 0
#' ## D3(5) 0 0 0 0 1 0
#' ## D3(7-) 0 0 0 0 1 0
#' ## D3(7+) 0 0 0 0 0 0
#' ## D3(7+)-D3(7-) 0 0 0 0 -1 0
#' ## D3(8) 0 0 0 0 0 0
#'
#' wald( fit, list( 'third derivatives' = Ls))
#'
#' ## Output:
#' ##
#' ## numDF denDF F.value p.value
#' ## third derivatives 1 15 2.184522 0.16009
#' ##
#' ## ? Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
#' ## D3(1) 0.000000 0.000000 15 2.408694 0.02932 0.000000 0.000000
#' ## D3(2) 0.000000 0.000000 15 2.408694 0.02932 0.000000 0.000000
#' ## D3(3-) 0.000000 0.000000 15 2.408694 0.02932 0.000000 0.000000
#' ## D3(3+) -0.535959 0.362621 15 -1.478013 0.16009 -1.308867 0.236950
#' ## D3(3+)-D3(3-) -0.535959 0.362621 15 -1.478013 0.16009 -1.308867 0.236950
#' ## D3(5) -0.535959 0.362621 15 -1.478013 0.16009 -1.308867 0.236950
#' ## D3(7-) -0.535959 0.362621 15 -1.478013 0.16009 -1.308867 0.236950
#' ## D3(7+) 0.000000 0.000000 Inf NaN NaN 0.000000 0.000000
#' ## D3(7+)-D3(7-) 0.535959 0.362621 15 1.478013 0.16009 -0.236950 1.308867
#' ## D3(8) 0.000000 0.000000 Inf NaN NaN 0.000000 0.000000
#'
#' # Note that some coefficients that are 0 by design may lead to invalid degrees
#' # of freedom and t-values.
#'
#' # sc generates a portion of a hypothesis matrix for the coefficients of
#' # a general spline constructed with gsp With:
#' #
#' # sc( sp, x, D, type ):
#' # sp is the spline function for which coefficients
#' # are required,
#' # x is the vector of values at which the spline is evaluated,
#' # D is the order of derivatives: 0 = value of spline,
#' # 1 = first derivative, etc.,
#' # type at knots: 0 limit from the left, 1 from the
#' # right, 2 is saltus (i.e. jump from left to right)
#'
#' # Warning: sc will not work correctly if the function defining the
#' # spline transforms the variable, e.g. sp <- function(x) gsp( x/100, knot=2 )
#'
#' # Example:
#'
#' 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)))
#'
#' sp <- function(x) gsp( x, knots = c(10,25,40),
#' degree = c(1,2,2,1),
#' smoothness = c(1,1,1))
#'
#' fit <- lm( y ~ sp(age)*G, simd)
#' xyplot( predict(fit) ~ age , simd, groups = G,type = "l")
#' summary(fit)
#'
#' L0 <- list(
#' "hat" = rbind(
#' "females at age=20" = c( 1, sc(sp,20), 0, 0 * sc(sp,20)),
#' "males at age=20" = c( 1, sc(sp,20), 1, 1* sc(sp,20))),
#' "male-female" = rbind(
#' "at 20" = c( 0 , 0*sc(sp,20), 1, 1*sc(sp,20))))
#' wald(fit, L0 )
#'
#' L1 <- list(
#' "D(yhat)/D(age)" =
#' rbind( "female at age = 25" = c(0, sc(sp,25,1), 0, 0*sc(sp,25,1)),
#' "male at x = 25" = c(0, sc(sp,25,1), 0, 1*sc(sp,25,1))))
#' wald( fit, L1)
#'
#' # A periodic spline function can be generated by forcing the coefficients
#' # beyond a periodic knot to repeat the pattern in a previous inteval. For
#' # example a periodic spline of period 1 can be created as follows:
#'
#' x <- seq(-3,3,.01)
#' y <- rnorm(x)
#'
#' per.sp <- function(x) gsp( x %% 1,
#' knots = 1,
#' degree = c(3, 3),
#' smoothness = 1,
#' lin = cbind( diag(4), - PolyShift(1,4)))
#' fit <- lm( y ~ per.sp(x))
#' summary(fit)
#' xyplot( y + predict(fit) ~ x ,
#' panel = panel.superpose.2,
#' type = c("p","l"))
#' xyplot( predict(fit) ~ x , type = 'l')
#'
#' # A periodic spline with additional knots can be created as shown below. Note
#' # that constraint matrix 'lin' expresses constraints for the 'full' polynomial
#' # parametrization, i.e. polynomials of degree 3 (thus of order 4 when
#' # including the constant term) in each of the 4 intervals. The constraint
#' # given in 'lin' forces the coefficients beyond the periodic knot at x = 1, to
#' # repeat the polynomial in the interval just to the right of x = 0.
#'
#' x <- seq(-3,3,.01)
#' y <- rnorm(x)
#'
#' per.sp <- function(x) gsp( x %% 1, knots = c(.333,.666,1),
#' degree = 3,
#' smoothness = 1,
#' lin = cbind( diag(4),0*diag(4),0*diag(4), - PolyShift(1,4)))
#' fit <- lm( y ~ per.sp(x))
#' summary(fit)
#' xyplot( y + predict(fit) ~ x ,
#' panel = panel.superpose.2,type = c("p","l"))
#' xyplot( predict(fit) ~ x , type = 'l')
#'
#'
#' ## Smoothing splines
#'
#' library(nlme)
#' library(MASS)
#' data(Spruce)
#' Spruce$all <- 1
#' range( Spruce$days)
#' sp <- function(x) smsp ( x, seq( 150, 675, 5))
#' spruce.fit1 <- lme(logSize ~ days, data=Spruce,
#' random=list(all= pdIdent(~sp(days) -1),
#' plot=~1, Tree=~1))
#' summary(spruce.fit1)
#' pred <- expand.grid( days = seq( 160, 670, 10), all = 1)
#' pred$logSize <- predict( spruce.fit1, newdata = pred, level = 1)
#' library( lattice )
#' xyplot( logSize ~ days, pred, type = 'l')
#'
#'
#' @export
gsp <- function (x, knots, degree = 3,
smoothness = pmax(pmin(degree[-1], degree[-length(degree)]) - 1, -1),
lin = NULL, periodic = FALSE, intercept = 0, signif = 3)
{
if (periodic) {
maxd <- max(degree)
maxk <- max(knots)
mid <- matrix(0, maxd + 1, (maxd + 1) * (length(knots) - 1))
const.per <- do.call(cbind, list(diag(maxd + 1), mid,
-PolyShift(maxk, maxd + 1)))
lin <- rbind(lin, const.per)
}
degree = rep(degree, length = length(knots) + 1)
smoothness = rep(smoothness, length = length(knots))
spline.attr <- list(knots = knots, degree = degree, smoothness = smoothness,
lin = lin, intercept = intercept, signif = signif)
if (is.null(x))
return(spline.attr)
if (periodic)
x <- x%%maxk
ret = Xf(x, knots, max(degree), signif = signif) %*% spline.E(knots,
degree, smoothness, lin = lin, intercept = intercept, signif = signif)
attr(ret, "spline.attr") <- spline.attr
class(ret) <- "gsp"
ret
}
print.gsp <- function(x, strip.attributes=TRUE, ...){
nms <- colnames(x)
ncol <- ncol(x)
xx <- x
if (strip.attributes) {
attributes(xx) <- NULL
xx <- matrix(xx, ncol=ncol)
colnames(xx) <- nms
}
else (xx <- unclass(xx))
print(zapsmall(xx), ...)
invisible(x)
}
#' Hypothesis matrix for general regression splines
#'
#' This function helps to build hypothesis matrices for general splines. See examples
#' in \code{\link[spida2]{gsp}}.
#'
#' @param sp a spline function generated by \code{\link[spida2]{gsp}}
#' @param x points where spline is tested
#' @param D (default 0) order of derivative to test, 0 is the value of the spline, 1 is the first derivative, etc.
#' @param type (default 1) if there is a discontinuity at a knot, type specifies whether to estimate the limit from the left (0), the limit from the right (1) or the saltus -- limit from the right minus the limit from left (2)
#' @seealso \code{\link[spida2]{wald}} \code{\link[spida2]{gsp}}
#' @export
sc <-
function( sp, x , D = 0, type = 1) {
a = sp(NULL)
D = rep(D, length.out = length(x))
type = rep ( type, length.out = length(x))
left = Xf( x, knots = a$knots, degree = max( a$degree), D = D, right = T)
right = Xf( x, knots = a$knots, degree = max( a$degree), D = D, right = F)
cleft = c(1,0,-1) [ match( type, c(0,1,2)) ]
cright = c(0,1,1) [ match( type, c(0,1,2)) ]
raw = left * cleft + right * cright
nam = rownames( raw )
nam = sub("^f","g", nam)
nam0 = sub( "\\)","-)", nam)
nam1 = sub( "\\)","+)", nam)
nam2 = paste( nam1 ,"-",nam0,sep = '')
rownames( raw ) =
ifelse( match( x , a$knots, 0) > 0,
cbind( nam0,nam1,nam2) [ cbind( seq_along(type), type+1)],
ifelse( type != 2, nam, '0'))
mod = raw %*% spline.E(
a$knots, a$degree, a$smoothness,
lin = a$lin,
intercept = a$intercept,
signif = a$signif)
mod
}
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.