smooth.construct.lmpi.smooth.spec: Locally shape-constrained P-spline based constructor (LSCOP-spline): locally increasing splines up to a change point.

Description

scam can use locally shape-constrained smooths that are monotone increasing up to a specified change point and become unconstrained beyond that point. This smooth is specified via model terms like s(x,bs="lmpi",xt=list(xc=xc)), where xc sets a change point. The construction of the smooths uses B-spline bases with mildly non-linear re-parametrization of the basis coefficients over the shape-constrained interval. The 'wiggliness' of the smooths is controlled by discrete penalties applied directly to the basis coefficients, in particular, by penalizing the squared differences between adjacent spline coefficients. lipl.smooth.spec works similarly to lpmi.smooth.spec, constructing a monotone increasing smooth up to a specified change point xc but plateaus beyond it.

The method is build by the mgcv constructor function for smooth terms, smooth.construct.

Usage

# S3 method for lmpi.smooth.spec
smooth.construct(object, data, knots)
# S3 method for lipl.smooth.spec
smooth.construct(object, data, knots)

Value

An object of class "lmpi.smooth".

Arguments

object: A smooth specification object, generated by an s term in a SCAM formula.
data: A data frame or list containing the data required by this term, with names given by object$term. The by variable is the last element.
knots: An optional list containing the knots supplied for basis setup. If it is NULL then the knot locations are generated automatically.

Author

Natalya Pya <nat.pya@gmail.com>

with contributions from Jens Lichter

Details

A smooth term of the form s(x,k,bs="lmpi",xt=list(xc=xc),m=2), where xc specifies a change point. k denotes the basis dimension and m+1 is the order of the B-spline basis. The default is a cubic spline m=2. The smooth uses a first order difference penalty on the coefficients starting from the second coefficient over the constrained part, and a second order difference penalty (a second order P-spline penalty) over the unconstrained part. The first order difference penalty over the constrained part shares with a second order P-spline penalty the feature of 'smoothing towards a straight line'. lipl.smooth.spec smooth sets zero values for the coefficients beyond the change point, which keeps the coefficient values constant (no change in the B-spline's coefficients beyond xc).

The default basis dimension is k=10. The basis dimensions for constrained and unconstrained parts are set proportionally to the share of each part from the total range of the covariate value, but with a minimum of 5 basis functions on each side.

The smooth currently does not support user-supplied knots. The knots are placed evenly throughout the constrained and unconstrained ranges of the covariate values (inner knots) and outside the covariate values range (outer knots). The change point xc is set as an inner knot m times. Having multiple knots at the change point guarantees the correct imposing of the shape constraints up to the change point, avoiding 'leakage' of the transformed basis functions of the constrained part into the unconstrained part.

Sum-to-zero ('centering') constraint is applied to the LSCOP-splines after imposing the shape-constrained model matrix transformation. Linear extrapolation is (only) used for prediction that requires extrapolation, i.e. prediction outside the range of the interior knots.

References

Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559

Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences

Examples

Run this code

## Local monotone increasing LSCOP-spline example...
## simulating data...
require(scam)
set.seed(4)
n <- 200
x <- sort(runif(n)*6)
f <- 4*tanh(2*x-5)-5*exp(-(x-4)^2) ## increasing up untill xt=2.978
xc <- 2.978
y <- f+rnorm(n)*.7
old.par <- par(mfrow=c(1,2))
plot(x,y,cex=.5)
lines(x,f,lty=2); abline(v=xc,lty=2,col="red")
## fit model ...
b <- scam(y~s(x,bs="lmpi",xt=list(xc=xc),k=15))
summary(b)
lines(x,fitted(b),col=2,lwd=2)
plot(b,shade=TRUE)
par(old.par)

## unconstrained model to compare with...
g <- scam(y~s(x))
plot(x,y,cex=.5)
lines(x,f,lty=2); abline(v=xc,lty=2,col=2)
lines(x,fitted(g),col=4,lwd=2)
lines(x,fitted(b),col=2,lwd=1)

## example with bs="lipl"...
set.seed(5)
n <- 200
x <- sort(runif(n)*4-1)
f <- 1.5*(1-exp(-2*x))
xc <- 1.5
y <- f+rnorm(n)*1.5
b <- scam(y~s(x,bs="lipl",xt=list(xc=xc)))
g <- scam(y~s(x,bs="ps")) ## unconstrained
old.par <- par(mfrow=c(2,2))
plot(x,y,cex=.5)
lines(x,f,lty=2); abline(v=xc,lty=2,col="red")
lines(x,fitted(b),col=2,lwd=2)
plot(b,shade=TRUE)
plot(g,shade=TRUE)
par(old.par)

Run the code above in your browser using DataLab