smooth.construct.symm.smooth.spec: Symmetric bivariate smooths constructor

In sparseFLMM: Functional Linear Mixed Models for Irregularly or Sparsely Sampled Data

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")

sparseFLMM documentation built on June 19, 2021, 5:06 p.m.