svcm: Fitting space-varying coefficient models

In svcm: 2d and 3d Space-Varying Coefficient Models

Description

'svcm' is used to fit a 2d or 3d space-varying coefficient model or to merely smooth the input data using penalized B-splines. The smoothing works either sequentially or multidimensionally involving tensor products. So far, only space-invariant regressors are allowed. Data must be on a regular grid without missings.

Usage

svcm(Y, X, vsize = c(1, 1, 1), knots = c(10, 10, 10),
     deg = c(1, 1, 1), opd = c(1, 1, 1), search = TRUE,
     lambda.init = rep(0.001, 3), method = "grid", type = "SEQ", ...)

Arguments

Y	array of observational data. Last dimension must correspond to the number of rows of X.
X	(r x p)-design matrix
vsize	numeric vector of the voxel size
knots	vector of the numbers of inner knots in x-, y- (and z-) direction
deg	vector of degrees of the basis functions in x-, y- (and z-) direction
opd	vector of the order of the difference penalties in x-, y- (and z-) direction
search	logical. If TRUE, the smoothing parameter will be optimized using method and GCV. If FALSE, lambda.init specifies the fixed smoothing parameter.
lambda.init	compulsory; initial value of global or dimension-specific smoothing parameter. See Details.
method	optimization method to be used. See Details.
type	character. "SEQ" (sequential) or "TP" (tensor product).

...	parameters to be passed to the optimization algorithm specified by method.

Details

The purpose of lambda.init is three-fold: First, the length determines the use of either global or dimension-specific penalties. Second, it serves as fixed smoothing parameter if search is deactivated. Third, it is used as initial value from the optim algorithm which runs in case of a multidimensional tuning parameter when no grid search is desired.

Unless method equals "grid", optimize is called in the case of a global tuning parameter requiring a specified interval to be passed. While optimize does not take a starting value explicitly, a startvalue can be passed to cleversearch, e.g. startvalue = lambda.init.

In the case of a dimension-specific tuning parameter, method "grid" evokes a full or greedy grid search (see cleversearch). Amongst others, simplex method "Nelder-Mead" or quasi-Newton "L-BFGS-B" with positivity constraint for the smoothing parameter are conceivable, too. For further specification see optim.

For simple smoothing of Y set X = matrix(1, 1, 1) and ascertain that the last dimension of Y matches dim(X)[1].

Value

A list with components:

fitted	fitted values as array of the same dimension as Y
effects	effects of dimension (n.x, n.y, p) resp. (n.x, n.y, n.z, p).
coeff	coefficients (amplitudes of the basis functions) of dimension (p, r.x, r.y) resp. (p, r.x, r.y, r.z) with r.x number of basis functions in x-direction.
knots	see knots.
deg	see deg.
opd	see opd.
vsize	see vsize.
type	character describing the SVCM variant used. See type.
call	the matched call.
opt	a list with components depending on search, i.e. on whether optimization was performed or not: time calculation/optimization time par initial value lambda.init or the best parameter found. value GCV value corresponding to par. GCVtab matrix of the search process with values of lambda and corresponding GCV value. ... further values returned by optim(), optimize()

Warnings

This model assumes the regressors to be space-invariant. Data must be on a regular grid without missings.

Background

In the general case of 2d mesh data, Dierckx (1982, 1993) demonstrates the equivalence of successive univariate smoothing with smoothing based on a full bivariate B-spline matrix. However, the equivalence does no longer hold if penalties are introduced. Dierckx proposes the so-called 'new smoothing spline' as approximation to the multidimensional penalized smoothing (type = "TP"). While Dierckx determines the penalty structure through the spline degree and the equidistance between adjacent knots, the present implementation (type = "SEQ") uses penalties of simple differences.

The calculation of GCV involves an inversion which is achieved using the recursion formula for band matrices in Hutchinson/de Hoog (1985). My collegue Thomas Kneib not only recommended this paper but also provided us with the basic.

Note

The observations in Y are assigned to the center of the respective grid unit sized vsize. Hence the basis functions are evaluated at these coordinates.

References

Dierckx P. (1982) A fast algorithm for smoothing data on a rectangular grid while using spline functions. SIAM Journal on Numerical Analysis 19(6), 1286-1304.

Dierckx P. (1993) Curve and surface fitting with splines. Oxford: Monographs on Numerical Analysis, Oxford University Press.

Heim S., Fahrmeir L., Eilers P. H. C. and Marx B. D. (2006) Space-Varying Coefficient Models for Brain Imaging. Ludwig-Maximilians-University, SFB 386, Discussion Paper 455.

Hutchinson M. F. and de Hoog F. R. (1985) Smoothing noisy data with spline functions. Journal of Numerical Mathematics 47, 99-106.

See Also

optimize for one-dimensional minimization,
optim here explicitly method "L-BGFS-B",
cleversearch for clever or full grid search.

Examples

## 2d DT-MRI data
data(brain2d)
X <- matrix(c(0.5, 0.5,   0,   0, 0.5, 0.5,
                0,   0, 0.5, 0.5, 0.5, 0.5,
              0.5, 0.5, 0.5, 0.5,   0,   0,
                0,   0,   0,   0,   1,  -1,
                1,  -1,   0,   0,   0,   0,
                0,   0,   1,  -1,   0,   0), nrow = 6)
##2d grid search for lambda; 60*50*6=18000 parameters (amplitudes) in total
M <- svcm(brain2d, X, knots=c(60, 50), deg=c(1, 1), vsize=c(1.875, 1.875),
          type="SEQ", lambda.init=rep(0.1, 2), search=TRUE,
          method="grid", lower=rep(-5, 2), upper=rep(0, 2), ngrid=10) 
str(M$opt)

## show results
zlim <- range(brain2d, M$fit)
old.par <- par(no.readonly=TRUE) 
par(pin=c(6, 1.2), mfrow=c(3, 1)) 
image(t(matrix(brain2d, nrow=dim(brain2d)[1])), axes=FALSE, zlim=zlim,
      col=grey.colors(256)) 
title("Observations: Six Diffusion Weighted Images")
image(t(matrix(M$fitted, nrow=dim(M$fit)[1])), axes=FALSE, zlim=zlim,
      col=grey.colors(256)) 
title("Fitted Values")
image(t(matrix(M$effects, nrow=dim(M$eff)[1])), axes=FALSE, 
      col=grey.colors(256))
title("Estimated Coefficient Surfaces: Six Elements of the Diffusion Tensor") 
par(old.par)


## 3d DT-MRI data; same regressors as in 2d; fixed global smoothing parameter
data(brain3d)
M3d <- svcm(brain3d, X, knots=c(5, 10, 5), deg=c(1, 1, 1), search=FALSE,
            vsize=c(1.875, 1.875, 4.0), type="TP", lambda.init=1)

## visualize results
zlim <- range(brain3d[,,,1], M3d$fit[,,,1])
old.par <- par(no.readonly = TRUE) 
par(pin=c(1.8, 5), mfrow=c(1, 3))
image(matrix(aperm(brain3d[,,,1], c(2,1,3)), nrow=dim(brain3d)[2]),
      axes=FALSE, zlim=zlim, col=grey.colors(256)) 
title("(a) Obs: 1st DWI")
image(matrix(aperm(M3d$fit[,,,1], c(2,1,3)), nrow=dim(brain3d)[2]),
      axes=FALSE, zlim=zlim, col=grey.colors(256)) 
title("(b) Fits of 1st DWI")
image(matrix(aperm(M3d$eff[,,,1], c(2,1,3)), nrow=dim(brain3d)[2]),
             axes=FALSE, col=grey.colors(256)) 
title("(c) Effects: 1st DT element")
title("Six axial slices of the 1st DWI-transform (a) and its fit (b);
      \n\n(c) corresponds to the first diffusion tensor component.",
       outer=TRUE, line=-5) 
par(old.par)

svcm documentation built on May 2, 2019, 1:29 p.m.