VGAM (version 0.9-0)

vsmooth.spline: Vector cubic smoothing spline

Description

Fits a vector cubic smoothing spline.

Usage

vsmooth.spline(x, y, w = NULL, df = rep(5, M), spar = NULL,
               all.knots = FALSE, iconstraint = diag(M),
               xconstraint = diag(M), 
               constraints = list("(Intercepts)" = diag(M), x = diag(M)), 
               var.arg = FALSE, scale.w = TRUE, nk = NULL,
               control.spar = list())

Arguments

x

A vector, matrix or a list. If a list, the x component is used. If a matrix, the first column is used. x may also be a complex vector, in which case the real part is used, and the imaginary part is used for the response.

y

A vector, matrix or a list. If a list, the y component is used. If a matrix, all but the first column is used. In this help file, M is the number of columns of y if there are no constraints on the functions.

w

The weight matrices or the number of observations. If the weight matrices, then this must be a n-row matrix with the elements in matrix-band form (see iam). If a vector, then these are the number of observations. By defa

df

Numerical vector containing the degrees of freedom for each component function (smooth). If necessary, the vector is recycled to have length equal to the number of component functions to be estimated (M if there are no constraints), which eq

spar

Numerical vector containing the non-negative smoothing parameters for each component function (smooth). If necessary, the vector is recycled to have length equal to the number of component functions to be estimated (M if there are no constra

all.knots

Logical. If TRUE then each distinct value of x will be a knot. By default, only a subset of the unique values of x are used; typically, the number of knots is O(n^0.25) for n large, but if

iconstraint

A M-row constraint matrix for the intercepts. It must be of full column rank. By default, the constraint matrix for the intercepts is the M by M identity matrix, meaning no constraints.

xconstraint

A M-row constraint matrix for x. It must be of full column rank. By default, the constraint matrix for the intercepts is the M by M identity matrix, meaning no constraints.

constraints

An alternative to specifying iconstraint and xconstraint, this is a list with two components corresponding to the intercept and x respectively. They must both be a M-row constraint matrix with full colum

var.arg

Logical: return the pointwise variances of the fit? Currently, this corresponds only to the nonlinear part of the fit, and may be wrong.

scale.w

Logical. By default, the weights w are scaled so that the diagonal elements have mean 1.

nk

Number of knots. If used, this argument overrides all.knots, and must lie between 6 and n+2 inclusive.

control.spar

See smooth.spline.

Value

• An object of class "vsmooth.spline" (see vsmooth.spline-class).

Details

The algorithm implemented is detailed in Yee (2000). It involves decomposing the component functions into a linear and nonlinear part, and using B-splines. The cost of the computation is O(n M^3).

The argument spar contains scaled smoothing parameters.

References

Yee, T. W. (2000) Vector Splines and Other Vector Smoothers. Pages 529--534. In: Bethlehem, J. G. and van der Heijde, P. G. M. Proceedings in Computational Statistics COMPSTAT 2000. Heidelberg: Physica-Verlag.

See Also

vsmooth.spline-class, plot.vsmooth.spline, predict.vsmooth.spline, iam, s, smooth.spline.

Examples

nn <- 20; x <- 2 + 5*(nn:1)/nn
x[2:4] <- x[5:7]      # Allow duplication
y1 <- sin(x) + rnorm(nn, sd = 0.13)
y2 <- cos(x) + rnorm(nn, sd = 0.13)
y3 <- 1 + sin(x) + rnorm(nn, sd = 0.13) # Run this for constraints
y <- cbind(y1, y2, y3)
ww <- cbind(rep(3, nn), 4, (1:nn)/nn)

(fit <- vsmooth.spline(x, y, w = ww, df = 5))
plot(fit) # The 1st and 3rd functions do not differ by a constant

mat <- matrix(c(1,0,1, 0,1,0), 3, 2)
(fit2 <- vsmooth.spline(x, y, w = ww, df = 5, iconstr = mat, xconstr = mat))
# The 1st and 3rd functions do differ by a constant:
mycols <- c("orange", "blue", "orange")
plot(fit2, lcol = mycols, pcol = mycols, las = 1)


p <- predict(fit, x = model.matrix(fit, type = "lm"), deriv = 0)
max(abs(fit@y - with(p, y))) # Should be zero

par(mfrow <- c(3, 1))
ux <- seq(1, 8, len = 100)
for(d in 1:3) {
    p <- predict(fit, x = ux, deriv = d)
with(p, matplot(x, y, type = "l", main = paste("deriv =", d),
                         lwd = 2, ylab = "", cex.axis = 1.5,
                         cex.lab = 1.5, cex.main = 1.5))
}