| smooth.construct.symm.smooth.spec {sparseFLMM} | R Documentation |
Symmetric bivariate smooths constructor
Description
The symm class is a smooth class that is appropriate for symmetric bivariate smooths, e.g. of covariance functions,
using tensor-product smooths in a gam formula. A constraint matrix is constructed
(see make_summation_matrix) to impose
a (skew-)symmetry constraint on the (cyclic) spline coefficients,
which considerably reduces the number of coefficients that have to be estimated.
Usage
## S3 method for class 'symm.smooth.spec'
smooth.construct(object, data, knots)
Arguments
object |
is a smooth specification object or a smooth object. |
data |
a data frame, model frame or list containing the values of the (named) covariates at which the smooth term is to be evaluated. |
knots |
an optional data frame supplying any knot locations to be supplied for basis construction. |
Details
By default a symmetric bivariate B-spline smooth g is specified,
in the sense that g(s, t) = g(t, s). By setting
s(..., bs = "symm", xt = list(skew = TRUE)), a skew-symmetric (or anti-smmetric)
smooth with g(s, t) = -g(t, s) can be specified instead.
In both cases, the smooth can also be constraint to be cyclic
with the property g(s, t) = g(s + c, t) = g(s, t + c)
for some fixed constant c via specifying xt = list(cyclic = TRUE).
Note that this does not correspond to specifying a tensor-product smooth from
cyclic marginal B-splines as given by the cp-smooth.
In the cyclic case, it is recommended to explicitly specify the range of the domain
of the smooth via the knots argument, as this determines the period and
often deviates from the observed range.
The underlying procedure is the following: First, the marginal spline design matrices and the corresponding marginal difference penalties are built. Second, the tensor product of the marginal design matrices is built and the bivariate penalty matrix is set up. Third, the constraint matrix is applied to the tensor product design matrix and to the penalty matrix.
Value
An object of class "symm.smooth". See smooth.construct for the elements it will contain.
Author(s)
Jona Cederbaum, Almond Stoecker
References
Cederbaum, Scheipl, Greven (2016): Fast symmetric additive covariance smoothing. Submitted on arXiv.
See Also
smooth.construct and smoothCon for details on constructors
Examples
require(sparseFLMM)
# (skew-)symmetric smooths ---------------------------------------
# generate random surface
dat1 <- data.frame(arg1 = 1:50)
dat2 <- expand.grid(arg1 = 1:50, arg2 = 1:50)
Bskew <- Predict.matrix(
smooth.construct(
s(arg1, arg2, bs = "symm", xt = list(skew = TRUE)),
data = dat2, knots = NULL ),
data = dat2 )
Bsymm <- Predict.matrix(
smooth.construct(
s(arg1, arg2, bs = "symm", xt = list(skew = FALSE)),
data = dat2, knots = NULL ),
data = dat2 )
set.seed(934811)
dat2$yskew <- c(Bskew %*% rnorm(ncol(Bskew)))
dat2$ysymm <- c(Bsymm %*% rnorm(ncol(Bsymm)))
# fit sum of skew-symmetric and symmetric parts with corresponding smooths
modpa <- gam( I(yskew + ysymm) ~ s(arg1, arg2, bs = "symm", xt = list(skew = TRUE)) +
s(arg1, arg2, bs = "symm", xt = list(skew = FALSE)), data = dat2)
# predict surfaces
preds <- predict(modpa, type = "terms")
dat1 <- as.list(dat1)
dat1$arg2 <- dat1$arg1
dat1$predskew <- matrix(preds[,1], nrow = length(dat1$arg1))
dat1$predsymm <- matrix(preds[,2], nrow = length(dat1$arg1))
cols <- hcl.colors(12, "RdBu")
opar <- par(mfcol = c(2,2))
# symm part (intercept missing)
with(dat1, image(arg1, arg2, predsymm, asp = 1,
main = "Symmetric part of y",
col = cols))
with(dat1, image(arg1, arg2, asp = 1,
main = "Fit via symm.smooth",
matrix(dat2$ysymm, nrow = length(arg1)),
col = cols))
# skew-symm part
with(dat1, image(arg1, arg2, predskew, asp = 1,
main = "Skew-symmetric part of y",
col = cols))
with(dat1, image(arg1, arg2, asp = 1,
main = "Fit via symm.smooth",
matrix(dat2$yskew, nrow = length(arg1)),
col = cols))
par(opar)
stopifnot(all.equal(dat1$predskew, - t(dat1$predskew)))
stopifnot(all.equal(dat1$predsymm, t(dat1$predsymm)))
# cyclic (skew-)symmetric splines ---------------------------------------
# fit the above example with cyclic smooths
modpac <- gam( I(yskew + ysymm) ~ s(arg1, arg2, bs = "symm",
xt = list(skew = TRUE, cyclic = TRUE)) +
s(arg1, arg2, bs = "symm", xt = list(skew = FALSE, cyclic = TRUE)),
knots = list(arg1 = c(1, 50), arg2 = c(1,50)),
# specify arg range to specify 'wavelength'!
data = dat2)
plot(modpac, asp = 1, se = FALSE, pages = 1)
predsc <- predict(modpac, type = "terms")
dat1$predskewc <- matrix(predsc[,1], nrow = length(dat1$arg1))
dat1$predsymmc <- matrix(predsc[,2], nrow = length(dat1$arg1))
# check cyclic margins
opar <- par(mfrow = c(1,2))
with(dat1, matplot(arg1, predsymmc[, c(1,10, 40)], t = "l",
main = "symmetric smooth"))
abline(h = dat1$predsymmc[1, c(1,10, 40)], col = "darkgrey")
abline(v = c(1,50), col = "darkgrey")
with(dat1, matplot(arg1, predskewc[, c(1,10, 40)], t = "l",
main = "skew-symmetric smooth"))
abline(h = dat1$predskewc[1, c(1,10, 40)], col = "darkgrey")
abline(v = c(1,50), col = "darkgrey")
par(opar)
# 1D point symmetric B-splines --------------------------------------------
# generate toy data
dat <- data.frame( x = 1:100 )
ps_obj <- with(dat, s(x, bs = "ps"))
B <- Predict.matrix(smooth.construct(ps_obj, dat, NULL), dat)
set.seed(3904)
dat$y <- B %*% rnorm(ncol(B))
plot(dat, t = "l")
# fit skew-symmetric spline
mod0 <- gam( y ~ s(x, bs = "symm", xt = list(skew = TRUE)),
knots = list(x = c(0,100)), # specify x range to determine inversion point
dat = dat )
lines(dat$x, predict(mod0), col = "cornflowerblue", lty = "dashed")
# or a symmetric spline to first part only
mod1 <- gam( y ~ s(x, bs = "symm"),
knots = list(x=c(0,50)),
dat = dat[1:50, ])
lines(dat[1:50, ]$x, predict(mod1), col = "darkred", lty = "dashed")