# 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

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

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

## 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: timecalculation/optimization time parinitial value `lambda.init` or the best parameter found. valueGCV value corresponding to `par`. GCVtabmatrix 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

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52``` ```## 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))), axes=FALSE, zlim=zlim, col=grey.colors(256)) title("Observations: Six Diffusion Weighted Images") image(t(matrix(M\$fitted, nrow=dim(M\$fit))), axes=FALSE, zlim=zlim, col=grey.colors(256)) title("Fitted Values") image(t(matrix(M\$effects, nrow=dim(M\$eff))), 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)), 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)), 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)), 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.