Nothing
General Parametric Splines in carEx"
In carEx: Supplemental and Experimental Functions
knitr::opts_chunk$set(echo = TRUE, comment='', message = FALSE)
library(carEx)
library(latticeExtra)
sp <- gspline(1)
Dmat <- environment(sp)$Dmat
Emat <- environment(sp)$Emat
basis <- environment(sp)$basis
Xmat <- environment(sp)$Xmat
Pcon <- environment(sp)$Pcon
Xf <- environment(sp)$Xf
Cmat <-
function( knots, degree, smoothness, lin = NULL, intercept = 0, signif = 3) {
# GM 2013-06-13
# add lin: contraints,
# generates constraint matrix
# GM: 2019_02_22: added facility to input smoothness as a list
# for smoothness constraints that don't have the form 0:smoothness[[i]]
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, F ) -
Xf( k, knots, dm, D = sm, T )
rownames( dmat ) <- paste( "C(",signif(k, signif),").",
sm, sep = '')
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 = F]
rownames( dmat ) = paste( "I.", i,".",seq(di+1,dm),sep = '')
cmat = rbind( cmat, dmat)
}
}
# add additional linear constraints
if( ! is.null(lin)) cmat <- rbind(cmat,lin) # GM:2013-06-13
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
}
Introduction
The parametric polynomial splines implemented in the 'carEx' package are
piecewise polynomial functions on $k+1$ intervals formed by $k$ knots partitioning
the real line:
$$(-\infty, t_1], (t_1, t_2], ..., (t_{i-1}, t_i], ..., (t_k, \infty)$$
with degree $d_i$ on the $i$th interval $(t_{i-1}, t_i]$, $i = 1,...,k+1$,
and order of continuity $c_i$ at the $i$th knot, $i = 1,...,k$.
The order of continuity refers to the highest order
for which the derivatives of the polynomial on the interval to the left and to the right of
a knot, $t_i$,
have the same limits at $t_i$.
For all orders above $c_i$, derivatives, if any, are not constained to have the same limit.
Such a spline is parametrized by three vectors:
a vector of knots, $t_1 < t_2 < ... < t_k$, of length $k > 0$,
a vector of polynomial degrees, $d_1, d_2, ... , d_{k+1}$, of length $k+1$, and
a vector of orders of continuity or 'smoothness', $c_1, c_2, ... , c_{k}$, of length $k$.
Theory
We first describe the general principles that underly the implemention of splines in this package.
Let $X_f$ be a $n \times q$ matrix for a model whose coefficients are subject
to $c$ linearly independent constraints given by a $c \times q$ matrix $C$. That is,
the linear space for the model is:
$$\mathcal{M} = { \eta = X_f \phi : \phi \in \mathbb{R}^q, C \phi = 0 }$$
We wish to construct a $n \times p$ design matrix $X$ with $p=q-c$ so that
$$\mathcal{M} = { \eta = X \beta : \beta \in \mathbb{R}^p}$$
Supppose further that we want the parameters $\beta$ to provide $p$ specified linearly independent
function of $\phi$ represented by the rows of the $p \times q$ matrix $E$ whose rows are
linearly independent of the rows of $C$ to ensure that they are not equal to 0 on
$\mathcal{M}$.
Consider the $q \times q$ partitioned matrix $\left[\begin{array}{cc} C \ E \end{array} \right]$.
Since its rows are linearly independent, it is invertible and has a conformably
partitioned inverse:
$$\left[\begin{array}{cc} F & G \end{array} \right] = \left[\begin{array}{cc} C \ E \end{array} \right]^{-1}$$
Thus $FC + GE = I$, $CF = I$, etc.
Consider the model matrix $X = X_f G$.
We show that $\mathcal{M} = {X\beta : \beta \in \mathbb{R}^p}$ and that
for any $\phi \in \mathbb{R}^q$, such that $C\phi =0$, $\beta = E \phi$.
Suppose $C \phi = 0$. Then
$$\phi = \left[\begin{array}{cc} F & G \end{array} \right] \left[\begin{array}{cc} C \ E \end{array} \right] \phi = \left[\begin{array}{cc} F & G \end{array} \right] \left[\begin{array}{cc} 0 \ E \phi \end{array} \right] = GE \phi$$
Thus, with $\beta = E \phi$, we have
$$ X_f \phi = X_f GE \phi = X \beta$$
We therefore have a 1-1 correspondence between $\beta \in \mathbb{R}^p$ and
${\phi \in \mathbb{R}^q : C \phi = 0}$ given by $\beta = E \phi$ and $\phi = G\beta$.
If $X$ is of full rank, we can obtain the least-squares estimator $\hat{\beta} = (X'X)^{-1}X'Y$.
We can then estimate any linear function $\psi = L\phi$ of $\phi$ under the constraint $C \phi = 0$ with the estimator $\hat{\psi} = A\hat{\beta}$ with
$$A = LG$$
Thus, the matrix $G$ serves as a post-multiplier to transform $X_f$ into a model matrix $X=X_f G$
that can be used in a linear model. The matrix $G$ also serves as a post-multiplier to
transform any general linear hypothesis matrix expressed in terms of $\phi$ into a
general linear hypothesis matrix in terms of of $\beta$.
Application to Splines
Our goal is to generate model matrices for splines in a way that produces interpretable
coefficients and lends itself to easily estimating and testing properties of the spline
that are linear functions of parameters: slope, curvature, discontinuities, etc.
Given $k$ knots, $-\infty = t_0 < t_1 < ... < t_k < t_{k+1} = \infty$, the spline
in the $i$th interval, $(t_{i-1},t_i]$, is a polynomial of degree $d_i$, a non-negative integer
with the value 0 signifying a constant over the corresponding interval.
The order of smoothness $c_i$ at $t_i$ is either a non-negative integer or -1
to allow a discontinuity. (TODO: control direction of discontinuity)
Generating a model matrix for some piecewise polynomial functions is simple.
For example, if the degrees, $d_i$, are
non-decreasing and the order of continuity is a constant $c$ less than $\min(d_i)$, one can add terms using 'plus' functions at each knot. For example, a
quadratic spline (degree 2, continuity 1) with one knot at 1 can be generated with a model
matrix with three columns, in addition to the intercept term:
$$x, x^2, (x -1)+^2$$
where
$$(y)+ = \left{ \begin{array}{ll}
0 & \textrm{if } y < 0 \
y & \textrm{otherwise }
\end{array}\right.$$
A spline that is quadratic on the interval $(-\infty,1]$ and cubic on $(1,\infty)$ with
continuity of order 1, $c_1 = 1$, at $t_1 = 1$, can be generated by the columns:
$$x, x^2, (x -1)+^2, (x -1)+^3$$
However, if one allows the degree of the polynomial or the order of smoothness to
vary in diffenret parts of the spline, the approach above works only in special cases.
Generating model matrices in more general situations, for example with degrees that are
not monotone, nor monotone increasing as the index radiates from a central value, is more
challenging. The approach described here works for any pattern of degrees, $d_i$ and smoothness
constraints, $c_i$.
We start by constucting a matrix, $X_f$, for a spline in which the polynomial degree in
each interval is the maximal value, $\max(d_i)$. We then construct constraints for the coefficients
of this model to produce the desired spline.
As an example, consider a spline, $\mathcal{S}$, with knots at 3 and 7, polynomial degrees,
$(2, 3, 2)$,
and smoothness, $(1, 2)$, meaning that $\mathcal{S}$ is smooth of order 1 at
$x=3$, and smooth of order 2 at $x=7$.
Columns of the full matrix $X_f$ contain the intercept,
linear and quadratic and cubic terms in each
interval of the spline.
To create an instance of $X_f$ we need to specify the values over which
the matrix is evaluated. Evaluating $X_f$ at $x = 0, 1, ... 9$, we obtain
the following matrix, which happens here to be block diagonal because of the
ordering of the $x$ values:
Xf(0:9, knots = c(3,7), degree = 3)
The model for the unconstrained maximal polynomial is
$X_f \phi: \phi \in \mathbb{R}^{12}$.
We impose three types of constraints on $\phi$.
- $X_f \phi$ should evaluate to 0 at $x=0$ so an intercept
term in the model will have the correct interpretation,
- the limits of the value and of the first derivative of the
spline must be the same when approaching the first knot
from the right or from the left, and the limits of the
value, the first and second derivatives should be the same
when approaching the second knot from the right or
from the left, and
- the degree of the polynomial in the first and third
intervals must be reduced to 2.
The constraint marix, $C$ is created by the 'Cmat' function:
Cmat(knots = c(3, 7), degree = c(2, 3, 2), smooth = c(1, 2))
The row labels of the constraint matrix show the role of each row.
For example, "f(0)" is
the value of the spline when $x=0$ which is constrained to 0 so that
an intercept term in a linear model can have its usual interpretation,
"C(3).0" ensures continuity
at $x=3$, "C(7).2" forces continuity of the
second derivative at $x = 7$, "I.1.3" constrains the cubic term
to be 0 in the first interval, etc.
Attributes give the length of the $\phi$ vector as 'npar.full',
the number of constraints as 'C.n', the rank of the constraint
matrix as 'C.rank' and the rank of the spline, omitting
the intercept term, as 'spline.rank'.
The 'd' attribute contains the vector of singular values of the
constraint matrix.
The following is the matrix $E$ of estimable functions created by the 'Emat'
function:
Emat(knots = c(3, 7), degree = c(2, 3, 2), smooth = c(1, 2))
The row labels signify the first derivative at $x=0$, 'D1(0)',
the second derivative at $x=0$, 'D2(0)', the saltus in the second
derivative at $x = 3$ and the saltus in the third derivative
at $x=3$.
The full rank model for the spline is generated by a
matrix $X= X_f G$ as described in the previous section.
The spline modelling function
is a closure generated by the gspline function.
sp <- gspline(knots = c(3, 7), degree = c(2, 3, 2), smoothness = c(1, 2))
sp(0:9)
produce a matrix $X = X_fG$ that will generate the desired spline parametrized by linear estimable coefficients.
The closure created by the gspline function can be used in a linear model formulas.
We illustrate its use with a small example. Note that the spline function can be used in
any linear model formula. It can, for example, be modelled as interacting with other predictors.
df <- data.frame(x = 0:10)
set.seed(123)
df <- within(df, y <- -2* (x-5) + .1 * (x-5)^3 + rnorm(x))
df <- rbind(df, data.frame(x = seq(0,10,.1), y = NA))
df <- df[order(df$x),]
plot(y~x, df, pch = 16)
fit <- lm(y ~ sp(x), data = df)
summary(fit)
lines(df$x , predict(fit, df))
Linear hypotheses
Linear hypotheses about a spline may be easy to formulate in terms
of its 'full' parameter vector
$\phi$ but challenging in terms of the 'working' parameters, $\beta$. For example, the
derivative or curvature of the spline over a range of values is easily expressed in terms
of $\phi$. To do this We use the relationship between linear hypotheses in terms of $\phi$ with those
in terms of $\beta$ to generate linear hypotheses based on $\hat{\beta}$. Namely the
least-squares estimator of
$\psi = L \phi$ under the contraint $C \phi =0$ is
$\hat{\psi} = A \hat{\beta}$ where $A = LG$.
Given a spline function sp created by the gspline function:
sp <- gspline(knots = c(3,7), degree = c(2,3,2), smoothness = c(1,2))
sp(0:9)
The sp function will generate a hypothesis matrix to query values and derivatives
of the spline.
sp(c(2, 3, 7), D = 1)
Denoting the matrix above by $A$, $A \hat{\beta}$ will estimate the first derivative
of the spline at $x=2$ and its limit from the right at the knots $x = 3, 7$.
The limit parameter to the spline function is used to select whether the value
estimated is a limit from the right, from the left, or the saltus (jump) in value
if discontinuous. For example, at $x=3$ where the spline has a discontinuous
second derivatives:
sp(c(3, 3, 3), D = 2, limit = c(-1,0,1))
Using the 'wald' function it is possible to graph these estimates as a function of
of $x$.
# customized panel function
gpanel.fit <- function (x, y, fit, lower, upper, subscripts, ..., type, group.number,
alpha, col, col.line, col.symbol, border = F, font, fontface) {
if (!missing(fit)) {
if (missing(col))
col <- "blue"
if (!missing(group.number)) {
col <- col.line
}
if (!missing(subscripts)) {
fit <- fit[subscripts]
}
dd <- data.frame(x = x, fit = fit)
dd <- dd[order(dd$x), ]
panel.lines(dd$x, dd$fit, ..., col = col, type = "l")
}
if (!missing(lower)) {
if (missing(alpha) || alpha == 1)
alpha <- 0.3
if (missing(col))
col <- "blue"
if (!missing(group.number)) {
col <- col.symbol
}
if (!missing(subscripts)) {
upper <- upper[subscripts]
lower <- lower[subscripts]
}
dd <- data.frame(x = x, lower = lower, upper = upper)
dd <- dd[order(dd$x), ]
panel.polygon(c(dd$x, rev(dd$x)), c(dd$upper, rev(dd$lower)),
border = border, col = col, alpha = alpha, ...)
}
}
xpred <- seq(0,10, .05)
A.1 <- cbind(0, sp(xpred, D = 1))
ww.1 <- as.data.frame(wald(fit, A.1))
A.2 <- cbind(0, sp(xpred, D = 2))
ww.2 <- as.data.frame(wald(fit, A.2))
plot(xpred, ww.1$coef, type = 'l')
plot(xpred, ww.2$coef, type = 'l')
library(latticeExtra)
ww.1$x <- xpred
xyplot(coef ~ x, ww.1, type = 'l',
lower = ww.1$L2, upper = ww.1$U2,
ylab = 'first derivative',
subscripts = TRUE) +
layer(gpanel.fit(...))
head(ww.1)
Periodic splines
We show how periodic splines can be used to fit periodic patterns such as seasonal
patterns. We use the monthly U.S. unemployment rates from January 1995 to February 2019.
The 'crash' in November 2008 creates a discontinuity in the series which we can use
to illustrate how to a discontinuity at an a priori value of time. We model the series
as consisting of two components, a long-term secular component and a periodic
annual component. We also illustrate how to use a parsimonious secular spline to model
secular interactions with the seasonal spline.
unemp <- Unemployment
head(unemp)
library(lattice)
library(latticeExtra)
xyplot(unemployment ~ date, unemp) + layer(panel.abline(v = as.Date('2008-12-15')))
toyear <- function(x) {
# number of years from January 1, 2000
(as.numeric(x) - as.numeric(as.Date('2000-01-01')))/365.25
}
unemp <- within(
unemp,
{
year <- toyear(date)
month <- as.numeric(format(date, '%m'))
})
summary(unemp)
quintiles <- quantile(unemp$year, 1:4/5)
sp2 <- gspline(quintiles, 2, 1) # quadratic spline
sp3 <- gspline(quintiles, 3, 2) # cubic spline
quintiles_with_crash <- sort(c(quintiles, toyear(as.Date('2008-12-15'))))
sp2d <- gspline(quintiles_with_crash, 2, c(1,1,-1,1,1))
sp3d <- gspline(quintiles_with_crash, 3, c(2,2,-1,2,2))
fit2 <- lm(unemployment ~ sp2(year), unemp)
summary(fit2)
unemp$fit2 <- predict(fit2)
fit3 <- lm(unemployment ~ sp3(year), unemp)
summary(fit3)
unemp$fit3 <- predict(fit3)
fit2d <- lm(unemployment ~ sp2d(year), unemp)
summary(fit2d)
unemp$fit2d <- predict(fit2d)
fit3d <- lm(unemployment ~ sp3d(year), unemp)
summary(fit3d)
unemp$fit3d <- predict(fit3d)
pp <- xyplot(unemployment ~ date, unemp, type = 'b',
col = 'gray',
key = list(
corner = c(1,1),
text = list(lab = c('quadratic','cubic','quadratic disc.','cubic disc.')),
lines = list(col= c('red','blue','red','blue'),
lty = c(1,1,3,3))
)) +
layer(panel.lines(x, unemp$fit3, col = 'blue')) +
layer(panel.lines(x, unemp$fit2, col = 'red')) +
layer(panel.lines(x, unemp$fit3d, col = 'blue', lty = 2)) +
layer(panel.lines(x, unemp$fit2d, col = 'red', lty = 2))
pp
Periodic spline
We add a periodic spline component as a function of months
using a cubic spline with period 12 and four internal knot at
$12 \times (1/5 2/5 3/5 4/5)$.
per3 <- gspline(12 * 1:5/5, 3, 2, periodic = TRUE)
fitper3 <- lm(unemployment ~ sp3d(year) + per3(month), unemp)
summary(fitper3)
unemp$fitper3 <- predict(fitper3)
pp <- xyplot(unemployment ~ date, unemp, type = 'b',
col = 'gray') +
layer(panel.lines(x, unemp$fitper3, col = 'blue'))
pp
We can examine the monthly periodic spline fit
pred <- data.frame(month = seq(0,24,.01), year = 0)
pred$fitper3 <- predict(fitper3, newdata = pred)
xyplot(fitper3 ~ month, pred, type = 'l',
xlim = c(0,24),
scales = list( x = list( at =3 * 0:8)),
ylab = 'seasonal unemployment') +
layer(panel.abline(v = 12))
Try the carEx package in your browser
Any scripts or data that you put into this service are public.
carEx documentation built on June 28, 2019, 3:01 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.
knitr::opts_chunk$set(echo = TRUE, comment='', message = FALSE) library(carEx) library(latticeExtra) sp <- gspline(1) Dmat <- environment(sp)$Dmat Emat <- environment(sp)$Emat basis <- environment(sp)$basis Xmat <- environment(sp)$Xmat Pcon <- environment(sp)$Pcon Xf <- environment(sp)$Xf Cmat <- function( knots, degree, smoothness, lin = NULL, intercept = 0, signif = 3) { # GM 2013-06-13 # add lin: contraints, # generates constraint matrix # GM: 2019_02_22: added facility to input smoothness as a list # for smoothness constraints that don't have the form 0:smoothness[[i]] 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, F ) - Xf( k, knots, dm, D = sm, T ) rownames( dmat ) <- paste( "C(",signif(k, signif),").", sm, sep = '') 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 = F] rownames( dmat ) = paste( "I.", i,".",seq(di+1,dm),sep = '') cmat = rbind( cmat, dmat) } } # add additional linear constraints if( ! is.null(lin)) cmat <- rbind(cmat,lin) # GM:2013-06-13 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 }
Introduction
The parametric polynomial splines implemented in the 'carEx' package are piecewise polynomial functions on $k+1$ intervals formed by $k$ knots partitioning the real line: $$(-\infty, t_1], (t_1, t_2], ..., (t_{i-1}, t_i], ..., (t_k, \infty)$$ with degree $d_i$ on the $i$th interval $(t_{i-1}, t_i]$, $i = 1,...,k+1$, and order of continuity $c_i$ at the $i$th knot, $i = 1,...,k$.
The order of continuity refers to the highest order for which the derivatives of the polynomial on the interval to the left and to the right of a knot, $t_i$, have the same limits at $t_i$. For all orders above $c_i$, derivatives, if any, are not constained to have the same limit.
Such a spline is parametrized by three vectors: a vector of knots, $t_1 < t_2 < ... < t_k$, of length $k > 0$, a vector of polynomial degrees, $d_1, d_2, ... , d_{k+1}$, of length $k+1$, and a vector of orders of continuity or 'smoothness', $c_1, c_2, ... , c_{k}$, of length $k$.
Theory
We first describe the general principles that underly the implemention of splines in this package.
Let $X_f$ be a $n \times q$ matrix for a model whose coefficients are subject to $c$ linearly independent constraints given by a $c \times q$ matrix $C$. That is, the linear space for the model is: $$\mathcal{M} = { \eta = X_f \phi : \phi \in \mathbb{R}^q, C \phi = 0 }$$ We wish to construct a $n \times p$ design matrix $X$ with $p=q-c$ so that $$\mathcal{M} = { \eta = X \beta : \beta \in \mathbb{R}^p}$$ Supppose further that we want the parameters $\beta$ to provide $p$ specified linearly independent function of $\phi$ represented by the rows of the $p \times q$ matrix $E$ whose rows are linearly independent of the rows of $C$ to ensure that they are not equal to 0 on $\mathcal{M}$.
Consider the $q \times q$ partitioned matrix $\left[\begin{array}{cc} C \ E \end{array} \right]$. Since its rows are linearly independent, it is invertible and has a conformably partitioned inverse: $$\left[\begin{array}{cc} F & G \end{array} \right] = \left[\begin{array}{cc} C \ E \end{array} \right]^{-1}$$ Thus $FC + GE = I$, $CF = I$, etc.
Consider the model matrix $X = X_f G$. We show that $\mathcal{M} = {X\beta : \beta \in \mathbb{R}^p}$ and that for any $\phi \in \mathbb{R}^q$, such that $C\phi =0$, $\beta = E \phi$.
Suppose $C \phi = 0$. Then
$$\phi = \left[\begin{array}{cc} F & G \end{array} \right] \left[\begin{array}{cc} C \ E \end{array} \right] \phi = \left[\begin{array}{cc} F & G \end{array} \right] \left[\begin{array}{cc} 0 \ E \phi \end{array} \right] = GE \phi$$ Thus, with $\beta = E \phi$, we have $$ X_f \phi = X_f GE \phi = X \beta$$ We therefore have a 1-1 correspondence between $\beta \in \mathbb{R}^p$ and ${\phi \in \mathbb{R}^q : C \phi = 0}$ given by $\beta = E \phi$ and $\phi = G\beta$.
If $X$ is of full rank, we can obtain the least-squares estimator $\hat{\beta} = (X'X)^{-1}X'Y$. We can then estimate any linear function $\psi = L\phi$ of $\phi$ under the constraint $C \phi = 0$ with the estimator $\hat{\psi} = A\hat{\beta}$ with $$A = LG$$ Thus, the matrix $G$ serves as a post-multiplier to transform $X_f$ into a model matrix $X=X_f G$ that can be used in a linear model. The matrix $G$ also serves as a post-multiplier to transform any general linear hypothesis matrix expressed in terms of $\phi$ into a general linear hypothesis matrix in terms of of $\beta$.
Application to Splines
Our goal is to generate model matrices for splines in a way that produces interpretable coefficients and lends itself to easily estimating and testing properties of the spline that are linear functions of parameters: slope, curvature, discontinuities, etc.
Given $k$ knots, $-\infty = t_0 < t_1 < ... < t_k < t_{k+1} = \infty$, the spline in the $i$th interval, $(t_{i-1},t_i]$, is a polynomial of degree $d_i$, a non-negative integer with the value 0 signifying a constant over the corresponding interval.
The order of smoothness $c_i$ at $t_i$ is either a non-negative integer or -1 to allow a discontinuity. (TODO: control direction of discontinuity)
Generating a model matrix for some piecewise polynomial functions is simple. For example, if the degrees, $d_i$, are non-decreasing and the order of continuity is a constant $c$ less than $\min(d_i)$, one can add terms using 'plus' functions at each knot. For example, a quadratic spline (degree 2, continuity 1) with one knot at 1 can be generated with a model matrix with three columns, in addition to the intercept term: $$x, x^2, (x -1)+^2$$ where $$(y)+ = \left{ \begin{array}{ll} 0 & \textrm{if } y < 0 \ y & \textrm{otherwise } \end{array}\right.$$ A spline that is quadratic on the interval $(-\infty,1]$ and cubic on $(1,\infty)$ with continuity of order 1, $c_1 = 1$, at $t_1 = 1$, can be generated by the columns: $$x, x^2, (x -1)+^2, (x -1)+^3$$ However, if one allows the degree of the polynomial or the order of smoothness to vary in diffenret parts of the spline, the approach above works only in special cases.
Generating model matrices in more general situations, for example with degrees that are not monotone, nor monotone increasing as the index radiates from a central value, is more challenging. The approach described here works for any pattern of degrees, $d_i$ and smoothness constraints, $c_i$.
We start by constucting a matrix, $X_f$, for a spline in which the polynomial degree in each interval is the maximal value, $\max(d_i)$. We then construct constraints for the coefficients of this model to produce the desired spline.
As an example, consider a spline, $\mathcal{S}$, with knots at 3 and 7, polynomial degrees, $(2, 3, 2)$, and smoothness, $(1, 2)$, meaning that $\mathcal{S}$ is smooth of order 1 at $x=3$, and smooth of order 2 at $x=7$. Columns of the full matrix $X_f$ contain the intercept, linear and quadratic and cubic terms in each interval of the spline.
To create an instance of $X_f$ we need to specify the values over which the matrix is evaluated. Evaluating $X_f$ at $x = 0, 1, ... 9$, we obtain the following matrix, which happens here to be block diagonal because of the ordering of the $x$ values:
Xf(0:9, knots = c(3,7), degree = 3)
The model for the unconstrained maximal polynomial is $X_f \phi: \phi \in \mathbb{R}^{12}$.
We impose three types of constraints on $\phi$.
- $X_f \phi$ should evaluate to 0 at $x=0$ so an intercept term in the model will have the correct interpretation,
- the limits of the value and of the first derivative of the spline must be the same when approaching the first knot from the right or from the left, and the limits of the value, the first and second derivatives should be the same when approaching the second knot from the right or from the left, and
- the degree of the polynomial in the first and third intervals must be reduced to 2.
The constraint marix, $C$ is created by the 'Cmat' function:
Cmat(knots = c(3, 7), degree = c(2, 3, 2), smooth = c(1, 2))
The row labels of the constraint matrix show the role of each row. For example, "f(0)" is the value of the spline when $x=0$ which is constrained to 0 so that an intercept term in a linear model can have its usual interpretation, "C(3).0" ensures continuity at $x=3$, "C(7).2" forces continuity of the second derivative at $x = 7$, "I.1.3" constrains the cubic term to be 0 in the first interval, etc.
Attributes give the length of the $\phi$ vector as 'npar.full', the number of constraints as 'C.n', the rank of the constraint matrix as 'C.rank' and the rank of the spline, omitting the intercept term, as 'spline.rank'.
The 'd' attribute contains the vector of singular values of the constraint matrix.
The following is the matrix $E$ of estimable functions created by the 'Emat' function:
Emat(knots = c(3, 7), degree = c(2, 3, 2), smooth = c(1, 2))
The row labels signify the first derivative at $x=0$, 'D1(0)', the second derivative at $x=0$, 'D2(0)', the saltus in the second derivative at $x = 3$ and the saltus in the third derivative at $x=3$.
The full rank model for the spline is generated by a matrix $X= X_f G$ as described in the previous section.
The spline modelling function
is a closure generated by the gspline function.
sp <- gspline(knots = c(3, 7), degree = c(2, 3, 2), smoothness = c(1, 2)) sp(0:9)
produce a matrix $X = X_fG$ that will generate the desired spline parametrized by linear estimable coefficients.
The closure created by the gspline function can be used in a linear model formulas.
We illustrate its use with a small example. Note that the spline function can be used in
any linear model formula. It can, for example, be modelled as interacting with other predictors.
df <- data.frame(x = 0:10) set.seed(123) df <- within(df, y <- -2* (x-5) + .1 * (x-5)^3 + rnorm(x)) df <- rbind(df, data.frame(x = seq(0,10,.1), y = NA)) df <- df[order(df$x),] plot(y~x, df, pch = 16) fit <- lm(y ~ sp(x), data = df) summary(fit) lines(df$x , predict(fit, df))
Linear hypotheses
Linear hypotheses about a spline may be easy to formulate in terms of its 'full' parameter vector $\phi$ but challenging in terms of the 'working' parameters, $\beta$. For example, the derivative or curvature of the spline over a range of values is easily expressed in terms of $\phi$. To do this We use the relationship between linear hypotheses in terms of $\phi$ with those in terms of $\beta$ to generate linear hypotheses based on $\hat{\beta}$. Namely the least-squares estimator of $\psi = L \phi$ under the contraint $C \phi =0$ is $\hat{\psi} = A \hat{\beta}$ where $A = LG$.
Given a spline function sp created by the gspline function:
sp <- gspline(knots = c(3,7), degree = c(2,3,2), smoothness = c(1,2)) sp(0:9)
The sp function will generate a hypothesis matrix to query values and derivatives
of the spline.
sp(c(2, 3, 7), D = 1)
Denoting the matrix above by $A$, $A \hat{\beta}$ will estimate the first derivative
of the spline at $x=2$ and its limit from the right at the knots $x = 3, 7$.
The limit parameter to the spline function is used to select whether the value
estimated is a limit from the right, from the left, or the saltus (jump) in value
if discontinuous. For example, at $x=3$ where the spline has a discontinuous
second derivatives:
sp(c(3, 3, 3), D = 2, limit = c(-1,0,1))
Using the 'wald' function it is possible to graph these estimates as a function of of $x$.
# customized panel function gpanel.fit <- function (x, y, fit, lower, upper, subscripts, ..., type, group.number, alpha, col, col.line, col.symbol, border = F, font, fontface) { if (!missing(fit)) { if (missing(col)) col <- "blue" if (!missing(group.number)) { col <- col.line } if (!missing(subscripts)) { fit <- fit[subscripts] } dd <- data.frame(x = x, fit = fit) dd <- dd[order(dd$x), ] panel.lines(dd$x, dd$fit, ..., col = col, type = "l") } if (!missing(lower)) { if (missing(alpha) || alpha == 1) alpha <- 0.3 if (missing(col)) col <- "blue" if (!missing(group.number)) { col <- col.symbol } if (!missing(subscripts)) { upper <- upper[subscripts] lower <- lower[subscripts] } dd <- data.frame(x = x, lower = lower, upper = upper) dd <- dd[order(dd$x), ] panel.polygon(c(dd$x, rev(dd$x)), c(dd$upper, rev(dd$lower)), border = border, col = col, alpha = alpha, ...) } } xpred <- seq(0,10, .05) A.1 <- cbind(0, sp(xpred, D = 1)) ww.1 <- as.data.frame(wald(fit, A.1)) A.2 <- cbind(0, sp(xpred, D = 2)) ww.2 <- as.data.frame(wald(fit, A.2)) plot(xpred, ww.1$coef, type = 'l') plot(xpred, ww.2$coef, type = 'l') library(latticeExtra) ww.1$x <- xpred xyplot(coef ~ x, ww.1, type = 'l', lower = ww.1$L2, upper = ww.1$U2, ylab = 'first derivative', subscripts = TRUE) + layer(gpanel.fit(...)) head(ww.1)
Periodic splines
We show how periodic splines can be used to fit periodic patterns such as seasonal patterns. We use the monthly U.S. unemployment rates from January 1995 to February 2019.
The 'crash' in November 2008 creates a discontinuity in the series which we can use to illustrate how to a discontinuity at an a priori value of time. We model the series as consisting of two components, a long-term secular component and a periodic annual component. We also illustrate how to use a parsimonious secular spline to model secular interactions with the seasonal spline.
unemp <- Unemployment head(unemp) library(lattice) library(latticeExtra) xyplot(unemployment ~ date, unemp) + layer(panel.abline(v = as.Date('2008-12-15'))) toyear <- function(x) { # number of years from January 1, 2000 (as.numeric(x) - as.numeric(as.Date('2000-01-01')))/365.25 } unemp <- within( unemp, { year <- toyear(date) month <- as.numeric(format(date, '%m')) }) summary(unemp) quintiles <- quantile(unemp$year, 1:4/5) sp2 <- gspline(quintiles, 2, 1) # quadratic spline sp3 <- gspline(quintiles, 3, 2) # cubic spline quintiles_with_crash <- sort(c(quintiles, toyear(as.Date('2008-12-15')))) sp2d <- gspline(quintiles_with_crash, 2, c(1,1,-1,1,1)) sp3d <- gspline(quintiles_with_crash, 3, c(2,2,-1,2,2)) fit2 <- lm(unemployment ~ sp2(year), unemp) summary(fit2) unemp$fit2 <- predict(fit2) fit3 <- lm(unemployment ~ sp3(year), unemp) summary(fit3) unemp$fit3 <- predict(fit3) fit2d <- lm(unemployment ~ sp2d(year), unemp) summary(fit2d) unemp$fit2d <- predict(fit2d) fit3d <- lm(unemployment ~ sp3d(year), unemp) summary(fit3d) unemp$fit3d <- predict(fit3d) pp <- xyplot(unemployment ~ date, unemp, type = 'b', col = 'gray', key = list( corner = c(1,1), text = list(lab = c('quadratic','cubic','quadratic disc.','cubic disc.')), lines = list(col= c('red','blue','red','blue'), lty = c(1,1,3,3)) )) + layer(panel.lines(x, unemp$fit3, col = 'blue')) + layer(panel.lines(x, unemp$fit2, col = 'red')) + layer(panel.lines(x, unemp$fit3d, col = 'blue', lty = 2)) + layer(panel.lines(x, unemp$fit2d, col = 'red', lty = 2)) pp
Periodic spline
We add a periodic spline component as a function of months using a cubic spline with period 12 and four internal knot at $12 \times (1/5 2/5 3/5 4/5)$.
per3 <- gspline(12 * 1:5/5, 3, 2, periodic = TRUE) fitper3 <- lm(unemployment ~ sp3d(year) + per3(month), unemp) summary(fitper3) unemp$fitper3 <- predict(fitper3) pp <- xyplot(unemployment ~ date, unemp, type = 'b', col = 'gray') + layer(panel.lines(x, unemp$fitper3, col = 'blue')) pp
We can examine the monthly periodic spline fit
pred <- data.frame(month = seq(0,24,.01), year = 0) pred$fitper3 <- predict(fitper3, newdata = pred) xyplot(fitper3 ~ month, pred, type = 'l', xlim = c(0,24), scales = list( x = list( at =3 * 0:8)), ylab = 'seasonal unemployment') + layer(panel.abline(v = 12))
Try the carEx package in your browser
Any scripts or data that you put into this service are public.
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.