Fitting spacevarying coefficient models
Description
'svcm' is used to fit a 2d or 3d spacevarying coefficient model or to merely smooth the input data using penalized Bsplines. The smoothing works either sequentially or multidimensionally involving tensor products. So far, only spaceinvariant regressors are allowed. Data must be on a regular grid without missings.
Usage
1 2 3 
Arguments
Y 
array of observational data. Last dimension must correspond to
the number of rows of 
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 
lambda.init 
compulsory; initial value of global or dimensionspecific smoothing parameter. See Details. 
method 
optimization method to be used. See Details. 
type 
character. 
... 
parameters to be passed to the optimization algorithm
specified by 
Details
The purpose of lambda.init
is threefold: First, the length
determines the use of either global or dimensionspecific
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 dimensionspecific tuning parameter, method
"grid"
evokes a full or greedy grid search (see
cleversearch
). Amongst others, simplex method
"NelderMead"
or quasiNewton "LBFGSB"
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

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 xdirection. 
knots 
see 
deg 
see 
opd 
see 
vsize 
see 
type 
character describing the SVCM variant used. See 
call 
the matched call. 
opt 
a list with components depending on

Warnings
This model assumes the regressors to be spaceinvariant. 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 Bspline matrix. However, the equivalence
does no longer hold if penalties are introduced. Dierckx proposes the
socalled '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), 12861304.
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) SpaceVarying Coefficient Models for Brain Imaging. LudwigMaximiliansUniversity, 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, 99106.
See Also
optimize
for onedimensional minimization,
optim
here explicitly method "LBGFSB"
,
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 DTMRI 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 DTMRI 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 DWItransform (a) and its fit (b);
\n\n(c) corresponds to the first diffusion tensor component.",
outer=TRUE, line=5)
par(old.par)
