svcm (version 0.1.2)

svcm: Fitting 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.

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.

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].

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)