svcm: Fitting space-varying coefficient models
In svcm: 2d and 3d Space-Varying Coefficient Models

Description Usage Arguments Details Value Warnings Background Note References See Also Examples

View source: R/svcm.R

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

1
2
3

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

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

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

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

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

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.

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.

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.

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

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