We revisit a multidimensional varying-coefficient model (VCM), by allowing regressor coefficients to vary smoothly in more than one dimension, thereby extending the VCM of Hastie and Tibshirani. The motivating example is 3-dimensional, involving a special type of nuclear magnetic resonance measurement technique that is being used to estimate the diffusion tensor at each point in the human brain. We aim to improve the current state of the art, which is to apply a multiple regression model for each voxel separately using information from six or more volume images. We present a model, based on P-spline tensor products, to introduce spatial smoothness of the estimated diffusion tensor. Since the regression design matrix is space-invariant, a 4-dimensional tensor product model results, allowing more efficient computation with penalized array regression

Eilers, Paul H. C.

Heim, Susanne

Marx, Brian D.

English

Universität München: Elektronischen Publikationen

Eilers, Heim, Marx:
Varying Coefficient Tensor Models for Brain Imaging
Sonderforschungsbereich 386, Paper 436 (2005)
Online unter: http://epub.ub.uni-muenchen.de/
Projektpartner
Varying Coefficient Tensor Models for Brain
Imaging
Paul H.C. Eilers, Susanne Heim, and Brian D. Marx
1 Department of Medical Statistics, Leiden Univerity Medical Center, 2300 RA,
Leiden, The Netherlands
2 Susanne Heim, Institute of Statistics, Ludwig-Maximilians-University, D-80539
Munich, Germany
3 Brian D. Marx, Department of Experimental Statistics, Louisiana State Uni-
versity, Baton Rouge, LA 70803 USA
Abstract:
We revisit a multidimensional varying-coefficient model (VCM), by allowing re-
gressor coefficients to vary smoothly in more than one dimension, thereby extend-
ing the VCM of Hastie and Tibshirani. The motivating example is 3-dimensional,
involving a special type of nuclear magnetic resonance measurement technique
that is being used to estimate the diffusion tensor at each point in the human
brain. We aim to improve the current state of the art, which is to apply a multi-
ple regression model for each voxel separately using information from six or more
volume images. We present a model, based on P-spline tensor products, to intro-
duce spatial smoothness of the estimated diffusion tensor. Since the regression
design matrix is space-invariant, a 4-dimensional tensor product model results,
allowing more efficient computation with penalized array regression.
Keywords: Brain imaging; P-splines; Varying coefficient model
1 Introduction
Varying-coefficient models (VCM) were introduced by Hastie and Tib-
shirani (1993) to handle situations in which regression coefficients vary
smoothly (interact) with another variable. Unlike standard regression, where
coefficients are assumed to be constant, VCMs produce coefficients that are
smooth curves, e.g. reflecting slow change over time or space. Hastie and
Tibshirani used iterative backfitting to estimate these curves, but Eilers and
Marx (2002) showed how to use P-splines to eliminate backfitting, in effect
estimating all varying coefficient terms simultaneously. In this paper, we
show how to use multidimensional P-splines (tensor products of B-splines,
combined with roughness penalties) to construct VCMs with coefficients
that smoothly change in 3-dimensional space.
The motivating application, i.e. estimation of the the diffusion tensor in
brain images from data obtained at a set of magnetic gradients, as well
2 Varying Coefficient Tensor Models
as a proposal for applying varying coefficient tensor models (VCTM) are
described in Heim et al. (2004). Presently linear regression is applied in
each image voxel separately; the six unique elements of the diffusion tensor
correspond to six coefficients in the regression. Noise can be reduced, and
possibly also acquisition time, if the coefficients can be forced to change
smoothly in space.
To allow enough detail in the model, relatively large numbers of B-spline
knots are needed. Also the number of the data cubes (M = 128×128×24) is
large. We propose to implement an efficient method for smoothing on grids
with tensor products that was developed by Eilers, Currie, and Durbàn
(2005). Such an approach avoids the computation of large Kronecker prod-
ucts of B-spline bases, and with a slight modification, this algorithm turns
the 3-dimensional VCM into smoothing with 4-dimensional tensor prod-
ucts. In the next section we briefly describe the physical background. The
VCM model and penalized estimation are described in Section 3. Array
regression is presented in Section 4. An illustrative example, and a variety
of details are discussed in Sections 5 and 6, respectively.
2 Diffusion tensor (DT) imaging
In general, the data basis consists of a set of diffusion weighted images
{Sr : r = 1, . . . , R} and a non–weighted reference image S0 at each of
M voxels of a specified 3-dimensional volume. The voxels are indexed in
space by s = (s1, s2, s3). Each response image is acquired at a different
uniform magnetic gradient. Diffusion is recorded along each gradient gr =
(g1r, g2r, g3r)
′, r = 1, . . . , R, as attenuation of the radio frequency signals
which are emitted by the atomic nuclei of the randomly travelling water
molecules in a magnetic resonance tomograph. The resulting covariance
ellipsoid of the local diffusivity is characterized by the voxel-wise diffusion
tensor
D(s) =


D1(s) D4(s) D5(s)
D4(s) D2(s) D6(s)
D5(s) D6(s) D3(s)

 .
Neuronal fibers in white matter are densely packed and highly ordered such
that the water molecules therein preferentially pass along the biophysiolog-
ical structures instead of perpendicular to them. Conversely, we determine
the principal diffusion direction with the dominant eigenvector of the local
diffusion tensor and identify it with the local fiber orientation in space.
The relation of the measured signal loss Sr(s) in voxel s and applied gra-
dient gr is given by the Stejskal-Tanner equation
Sr(s) = S0(s) exp {−bgr
′D(s)gr} , r = 1, . . . , R . (1)
Here, b is a scalar comprising several acquisition parameters such as gradi-
ent strength and duration.
Eilers, Marx, and Heim 3
To fully determine the six free parameters of the DT in (1), one has to
apply gradients in r ≥ 6 independent directions. Usually, measurements
are repeated for the same gradient set or the number of different gradients
are chosen to be over determined in order to mitigate the effects of noise.
Following Papadakis et al. (1999), (1) can be reformulated, for voxel s, as
yr(s) = −
1
b
log
(
Sr
S0
(s)
)
= xr
′β(s) , (2)
with the unknown elements of the diffusion tensor
β(s) = (D1, D2, D3, D4, D5, D6)
′(s) and a nonlinear mapping of the gradi-
ent to a six-dimensional vector xr = (g
2
1r, g
2
2r, g
2
3r, 2g1rg2r, 2g1rg3r, 2g2rg3r)
′.
In (2) we recognize a linear regression problem, in which the 6-dimensional
vector β(s) has to estimated for each voxel, using 6 or more images obtained
at different magnetic gradients. Note that xr is the same for all voxels, and
thus we face a repeated measures design.
3 The varying-coefficient model
The current state-of-the-art method of estimating the diffusion matrix from
DT imaging experiments is through voxel-wise multiple regressions. Yet the
voxel spatial structure encourages an assumption of smooth coefficients. A
major benefit of smoothed coefficients is that the DT can be approximated
at any arbitrary position within the recorded volume. Consequently, a fiber
tracking algorithm is allowed to operate on a markedly finer grid, which
should result in more precise neuronal tract reconstruction. Hence, the βt
(t = 1, . . . , T = 6) could be preferably modelled non-parametrically via
penalized tensor product B-splines or multidimensional P-splines.
We first present a direct VCM presentation. Consider
βt(s) =
KLH
∑
v=1
Tv(s)γtv = T(s)γt, (3)
where T(s) denotes the 3-dimensional tensor product of 1-dimensional B-
splines, evaluated at voxel s, γt is the vector of unknown amplitudes, and
K × L × H is determined by the (generous and regularly gridded) knot
partition and degree of the basis function.
Using tensor coefficient expression in (3), we aim to find a practical solution
to the penalized objective
Q = ‖y −
T
∑
t=1
Υxt(s)T(s)γt‖
2 + Penalty, (4)
where Υxt(s) = diag{xt(s)}. The penalty term places difference penalties
on the rows, columns, and layers of tensor product coefficients, for each
4 Varying Coefficient Tensor Models
regressor t, such that
Penalty =
T
∑
t=1
K
∑
k=1
L
∑
l=1
H
∑
h=1
{λ1(∆
d
1γtklh)
2 + λ2(∆
d
2γtklh)
2 + λ3(∆
d
3γtklh)
2}.
(5)
Here ∆d1 (∆
d
2, ∆
d
3 ) denotes the d-th order differences across a row (down a
column, along a layer) of the K ×L×H matrix of tensor product B-spline
coefficients, Γt = [γklh]. In theory, an explicit solution to γ = (γ
′
1, . . . , γ
′
T )
′
can be found using
γ̂ = (U′U + P)
−1
U′y, (6)
where U = [Υx1(s)T, . . . ,ΥxT (s)T], and P consists of a carefully arranged
matrix representation (using block diagonal matrices of Kronecker prod-
ucts) of the differences found in (5). Given γ, then the varying coefficient
volumes can be built. For practical brain imaging applications, we may
need, e.g., K × L × H = 32 × 32 × 8 knots, thus T has approximately
3 × 109 elements: if each floating point takes 8 bytes, then T will use sev-
eral of Gb memory, which is beyond reach of current computers.
As the straight-forward approach to VCM runs into difficulties, we can alle-
viate computation by taking advantage of the repeated measures structure,
i.e. x does not depend on s or vary across voxels. Since the x are on a grid,
the mean can be expressed as a tensor product involving x. Hence we can
apply the fast, compact smoothing algorithm of Eilers, Currie, and Durbàn
(2005), which uses a relatively moderate amount of memory, making the
brain VCM model tractable again.
4 Array regression applied to brain imaging data
To simplify the presentation, we will consider a VCM that models only data
in one dimension, e.g. assuming smooth regression coefficients along a line
of voxels. We then will move into two dimensional smooth coefficients using
slice images. To start, we first present the current approach of voxel-wise
multiple regression. For each voxel, we model the mean E(yr) =
∑T
t=1 xrtβt
or E(y) = Xβ, where in our example r = 1, . . . , R and T = 6 given by the
magnetic field gradients. The matrix X has dimension R×T , where R ≥ T ,
with normal equations X ′Xβ̂ = X ′y.
With a line of I voxels, the y vectors can be placed into a matrix Y of
dimension R × I, as well as the coefficients into a matrix A of dimension
T×I. Thus Y = XA+E, where E is a matrix of random (normal errors). To
ensure smoothness on the rows of A, we project each row onto a (common)
B-spline basis, and hence A′ = BΓ, where B (I × K) of degree q is built
along the index 1, . . . , I, and Γ is of dimension K×T . Thus we have E(Y ) =
XΓ′B′, which shows directly the tensor structure of X and B. There is only
one set of penalties on Γ: one penalty for each row. Thus the penalty term
in (5) is greatly simplified, as it only involves sums over t and k.
Eilers, Marx, and Heim 5
With 2-dimensional brain slices, instead of voxels, we conventionally think
of pixels indexed by i = 1, . . . , I and j = 1, . . . , J . Standard matrix notation
now breaks down, as Y and Γ become 3-dimensional arrays. The R images
at each pixel are denoted by Y = [yrij ] r = 1, . . . , R and the gradient
information by X = [xrt], t = 1, . . . , T . If we arrange Y as a matrix with
R rows and IJ columns (each column corresponding to one pixel), we can
write Y = XA + E, with a T × IJ coefficient matrix A = [αtij ] and E
is a R × IJ matrix of random (normal) errors. In this format, pixel-wise
regression would simply lead to the normal equations X ′XÂ = X ′Y . We
again impose spatial smoothness of the elements of A by projecting them
onto (K × L) tensor products of B-splines. Let B = [bik] and B̃ = [b̃jl]
(k = 1, . . . ,K; l = 1, . . . , L) be B-spline bases, of degree q, along the two
axes of the image. Thus αtij =
∑
k
∑
l bik b̃jlγklt. For simplicity we use the
same tensor product B-spline basis for each t. This leads to
µrij = E(yrij) =
T
∑
t=1
xrt
K
∑
k=1
L
∑
l=1
bik b̃jlγtkl, (7)
in which, by virtue of space-invariant x, we recognize a 3-dimensional tensor
product model. Again, we ensure further smoothness through additional
difference penalties on the coefficients of the tensor products. The simplified
version of (5) is written as
Penalty =
T
∑
t=1
{λ1
K
∑
k=1
L
∑
l=1
(∆d1γtkl)
2 + λ2
K
∑
k=1
L
∑
l=1
(∆d2γtkl)
2}. (8)
Here ∆d1 (∆
d
2) denotes the d-th order differences along a row (down a col-
umn) of the T × K × L array of tensor coefficients, Γ = [γtkl].
Eilers et al. (2005) showed how to fit a multidimensional tensor product
model with penalties efficiently. Their method adapts to any number of di-
mensions, and thus can be applied to the original 3-dimensional DT varying
coefficient model. Lack of space does not allow us to go into the details of
3-dimensional extensions, but can be found in Eilers et al. (2005). In this
setting, both Y and Γ will become 4-dimensional arrays, and the B-spline
basis will come in a triple: B, B̃, B̆. The mean of Y can again be expressed
in a tensor form, E(Y ) = XΓ?B?, where Γ? and B? are both carefully
arranged matrices involving Kronecker products.
5 Illustrative example
The example consists of the 128×128 middle axial slice of a brain volume of
diffusion tensor images. First we consider 1-dimensional VCM models. Fig-
ure 1 (top panel) shows the response signals for the 75th row of the brain
slice and the varying coefficients for the diagonal tensor elements (mid-
dle panel) and the off-diagonal elements (bottom panel). Figure 2 shows
6 Varying Coefficient Tensor Models
 
S
ig
na
l
40 60 80 100
0.
0
0.
00
2
0.
00
4
Signals at six gradients
 
S
ig
na
l
40 60 80 100
-0
.0
01
0.
00
1
Smooth regression: diagonal tensor products
 
S
ig
na
l
40 60 80 100
-0
.0
00
6
0.
0
0.
00
06
Smooth regression: off-diagonal tensor products 
FIGURE 1. Signals at six gradients and smoothed DT coefficients.
the corresponding coefficients using pixel-wise multiple regression. Figure
3 (left, top and bottom) presents example 32×32 images from the first and
third gradients. The center panels present the diagonal D1(s) (center) and
the off-diagonal D4(s) element pixel-wise. The right panels presents the
corresponding DT elements computed with VCM models and interpolated
to a 128 × 128 resolution; thus demonstrating a benefit from an analytic
P-spline approximation.
6 Details
In all cases, cubic B-splines, and a second order difference penalty, with
the penalty parameter set to 0.01, were used. We used a knot at every
second pixel for both the 1D and 2D cases (in the latter we analyzed only
a 32 × 32 image). Unfortunately, the demands on computer memory and
computation time load increase sharply with knot density. For 3D images
of the original size it will essentially be impossible to explicitly form and
solve the penalized normal equation. We would have to deal with about
50 by 50 by 10 by 6 knots, or 1.5 105 equations. We will research direct
iterative methods to circumvent theses difficulties. Our approach naturally
accommodates observation weights, so that missing data does not present
a problem, and adjustments can be made for voxels with uncertainty, e.g.
Eilers, Marx, and Heim 7
40 60 80 100
0.
0
0.
00
1
0.
00
3
Pixel-wise regression: diagonal tensor elements
40 60 80 100
-0
.0
00
6
0.
0
0.
00
04
Pixel-wise regression: off-diagonal tensor elements
FIGURE 2. Estimated pixel-wise DT coefficients.
at the boundaries. Further as parameter estimation is grounded in classical
methods, model summaries, such as: effective dimension and information
criteria, and delete-one diagnostics, such as cross-validation (CV) are eas-
ily obtained. Accessible CV can also be a practical matter when searching
for the “optimal” penalty tuning parameters, e.g. when isotropy is relaxed.
Anisotropic penalization can be achieved by placing different λs on the
indexing axes of Γ, leading to a higher dimensional grid-search. We are
currently researching methods to optimally choose the penalty parameters
(λ), including an E-M algorithm that alternates between: a) estimating
the residual variance and the var(γ̂t) to update penalty parameters, and
b) estimating Γ with these updated penalty parameters. Since we are ul-
timately interested in plotting the smooth coefficients, we are also investi-
gating choosing λ based on mean square error of the coefficient estimates,
as an alternative prediction oriented criteria. Further, our penalized VCM
approach can be transplanted into a generalized linear model framework,
e.g. smoothing coefficients associated with a binary or threshold response.
References
Basser, P.J and Jones, D.K. (2002). Diffusion-tensor MRI: theory, experi-
mental design and data analysis – a technical review. NMR Biomed,
15, 456-467.
8 Varying Coefficient Tensor Models
Image at gradient 1
Image at gradient 3
Smooth interpolated DT element 1,1
Smooth interpolated DT element 1,2
Pixel−wise DT element 1,1
Pixel−wise DT element 1,2
FIGURE 3. Two 32 by 32 pixel example images (left, top and bottom), two DT
elements, as computed by pixel-wise regression (middle, top and bottom) and the
same DT elements as computed with VCM model, and interpolated to a 128 by
128 grid (right, top and bottom).
Eilers, P.H.C, Currie, I., and Durbàn, M. (2005). Low memory, high speed
smoothing on large multidimensional grids. Computational Statistics
and Data Analysis. In Press.
Eilers, P.H.C. and Marx, B.D. (2002). Generalized linear additive smooth
structures. J. of Computat. and Graphical Statistics, 11(4), 758-783.
Hastie, T. and Tibshirani, R. (1993). Varying-Coefficient Models (with dis-
cussion). Journal of the Royal Statistical Society B 55: 757-796.
Heim, S. Hahn, K., Auer, D.P., Fahrmeir, L. (2004). Proceedings of the 19th
IWSM, Florence, Italy.
Papadakis, N.G., Xing, D. Huang, CL-H, Hall, L., and Carpenter, T.A. (1999).
A comparative study of acquisition schemes for diffusion tensor imag-
ing using MRI. J Magn Reson, 137, 67-82.


Varying Coefficient Tensor Models for Brain Imaging

Open Access LMU

Eilers, Heim, Marx:Varying Coefficient Tensor Models for Brain ImagingSonderforschungsbereich 386, Paper 436 (2005)Online unter: http://epub.ub.uni-muenchen.de/ProjektpartnerVarying Coefficient Tensor Models for BrainImagingPaul H.C. Eilers, Susanne Heim, and Brian D. Marx1 Department of Medical Statistics, Leiden Univerity Medical Center, 2300 RA,Leiden, The Netherlands2 Susanne Heim, Institute of Statistics, Ludwig-Maximilians-University, D-80539Munich, Germany3 Brian D. Marx, Department of Experimental Statistics, Louisiana State Uni-versity, Baton Rouge, LA 70803 USAAbstract:We revisit a multidimensional varying-coefficient model (VCM), by allowing re-gressor coefficients to vary smoothly in more than one dimension, thereby extend-ing the VCM of Hastie and Tibshirani. The motivating example is 3-dimensional,involving a special type of nuclear magnetic resonance measurement techniquethat is being used to estimate the diffusion tensor at each point in the humanbrain. We aim to improve the current state of the art, which is to apply a multi-ple regression model for each voxel separately using information from six or morevolume images. We present a model, based on P-spline tensor products, to intro-duce spatial smoothness of the estimated diffusion tensor. Since the regressiondesign matrix is space-invariant, a 4-dimensional tensor product model results,allowing more efficient computation with penalized array regression.Keywords: Brain imaging; P-splines; Varying coefficient model1 IntroductionVarying-coefficient models (VCM) were introduced by Hastie and Tib-shirani (1993) to handle situations in which regression coefficients varysmoothly (interact) with another variable. Unlike standard regression, wherecoefficients are assumed to be constant, VCMs produce coefficients that aresmooth curves, e.g. reflecting slow change over time or space. Hastie andTibshirani used iterative backfitting to estimate these curves, but Eilers andMarx (2002) showed how to use P-splines to eliminate backfitting, in effectestimating all varying coefficient terms simultaneously. In this paper, weshow how to use multidimensional P-splines (tensor products of B-splines,combined with roughness penalties) to construct VCMs with coefficientsthat smoothly change in 3-dimensional space.The motivating application, i.e. estimation of the the diffusion tensor inbrain images from data obtained at a set of magnetic gradients, as well2 Varying Coefficient Tensor Modelsas a proposal for applying varying coefficient tensor models (VCTM) aredescribed in Heim et al. (2004). Presently linear regression is applied ineach image voxel separately; the six unique elements of the diffusion tensorcorrespond to six coefficients in the regression. Noise can be reduced, andpossibly also acquisition time, if the coefficients can be forced to changesmoothly in space.To allow enough detail in the model, relatively large numbers of B-splineknots are needed. Also the number of the data cubes (M = 128×128×24) islarge. We propose to implement an efficient method for smoothing on gridswith tensor products that was developed by Eilers, Currie, and Durba`n(2005). Such an approach avoids the computation of large Kronecker prod-ucts of B-spline bases, and with a slight modification, this algorithm turnsthe 3-dimensional VCM into smoothing with 4-dimensional tensor prod-ucts. In the next section we briefly describe the physical background. TheVCM model and penalized estimation are described in Section 3. Arrayregression is presented in Section 4. An illustrative example, and a varietyof details are discussed in Sections 5 and 6, respectively.2 Diffusion tensor (DT) imagingIn general, the data basis consists of a set of diffusion weighted images{Sr : r = 1, . . . , R} and a non–weighted reference image S0 at each ofM voxels of a specified 3-dimensional volume. The voxels are indexed inspace by s = (s1, s2, s3). Each response image is acquired at a differentuniform magnetic gradient. Diffusion is recorded along each gradient gr =(g1r, g2r, g3r)′, r = 1, . . . , R, as attenuation of the radio frequency signalswhich are emitted by the atomic nuclei of the randomly travelling watermolecules in a magnetic resonance tomograph. The resulting covarianceellipsoid of the local diffusivity is characterized by the voxel-wise diffusiontensorD(s) =D1(s) D4(s) D5(s)D4(s) D2(s) D6(s)D5(s) D6(s) D3(s) .Neuronal fibers in white matter are densely packed and highly ordered suchthat the water molecules therein preferentially pass along the biophysiolog-ical structures instead of perpendicular to them. Conversely, we determinethe principal diffusion direction with the dominant eigenvector of the localdiffusion tensor and identify it with the local fiber orientation in space.The relation of the measured signal loss Sr(s) in voxel s and applied gra-dient gr is given by the Stejskal-Tanner equationSr(s) = S0(s) exp {−bgr′D(s)gr} , r = 1, . . . , R . (1)Here, b is a scalar comprising several acquisition parameters such as gradi-ent strength and duration.Eilers, Marx, and Heim 3To fully determine the six free parameters of the DT in (1), one has toapply gradients in r ≥ 6 independent directions. Usually, measurementsare repeated for the same gradient set or the number of different gradientsare chosen to be over determined in order to mitigate the effects of noise.Following Papadakis et al. (1999), (1) can be reformulated, for voxel s, asyr(s) = −1blog(SrS0(s))= xr′β(s) , (2)with the unknown elements of the diffusion tensorβ(s) = (D1, D2, D3, D4, D5, D6)′(s) and a nonlinear mapping of the gradi-ent to a six-dimensional vector xr = (g21r, g22r, g23r, 2g1rg2r, 2g1rg3r, 2g2rg3r)′.In (2) we recognize a linear regression problem, in which the 6-dimensionalvector β(s) has to estimated for each voxel, using 6 or more images obtainedat different magnetic gradients. Note that xr is the same for all voxels, andthus we face a repeated measures design.3 The varying-coefficient modelThe current state-of-the-art method of estimating the diffusion matrix fromDT imaging experiments is through voxel-wise multiple regressions. Yet thevoxel spatial structure encourages an assumption of smooth coefficients. Amajor benefit of smoothed coefficients is that the DT can be approximatedat any arbitrary position within the recorded volume. Consequently, a fibertracking algorithm is allowed to operate on a markedly finer grid, whichshould result in more precise neuronal tract reconstruction. Hence, the βt(t = 1, . . . , T = 6) could be preferably modelled non-parametrically viapenalized tensor product B-splines or multidimensional P-splines.We first present a direct VCM presentation. Considerβt(s) =KLH∑v=1Tv(s)γtv = T(s)γt, (3)where T(s) denotes the 3-dimensional tensor product of 1-dimensional B-splines, evaluated at voxel s, γt is the vector of unknown amplitudes, andK × L × H is determined by the (generous and regularly gridded) knotpartition and degree of the basis function.Using tensor coefficient expression in (3), we aim to find a practical solutionto the penalized objectiveQ = ‖y −T∑t=1Υxt(s)T(s)γt‖2 + Penalty, (4)where Υxt(s) = diag{xt(s)}. The penalty term places difference penaltieson the rows, columns, and layers of tensor product coefficients, for each4 Varying Coefficient Tensor Modelsregressor t, such thatPenalty =T∑t=1K∑k=1L∑l=1H∑h=1{λ1(∆d1γtklh)2 + λ2(∆d2γtklh)2 + λ3(∆d3γtklh)2}.(5)Here ∆d1 (∆d2, ∆d3 ) denotes the d-th order differences across a row (down acolumn, along a layer) of the K ×L×H matrix of tensor product B-splinecoefficients, Γt = [γklh]. In theory, an explicit solution to γ = (γ′1, . . . , γ′T )′can be found usingγˆ = (U′U + P)−1U′y, (6)where U = [Υx1(s)T, . . . ,ΥxT (s)T], and P consists of a carefully arrangedmatrix representation (using block diagonal matrices of Kronecker prod-ucts) of the differences found in (5). Given γ, then the varying coefficientvolumes can be built. For practical brain imaging applications, we mayneed, e.g., K × L × H = 32 × 32 × 8 knots, thus T has approximately3× 109 elements: if each floating point takes 8 bytes, then T will use sev-eral of Gb memory, which is beyond reach of current computers.As the straight-forward approach to VCM runs into difficulties, we can alle-viate computation by taking advantage of the repeated measures structure,i.e. x does not depend on s or vary across voxels. Since the x are on a grid,the mean can be expressed as a tensor product involving x. Hence we canapply the fast, compact smoothing algorithm of Eilers, Currie, and Durba`n(2005), which uses a relatively moderate amount of memory, making thebrain VCM model tractable again.4 Array regression applied to brain imaging dataTo simplify the presentation, we will consider a VCM that models only datain one dimension, e.g. assuming smooth regression coefficients along a lineof voxels. We then will move into two dimensional smooth coefficients usingslice images. To start, we first present the current approach of voxel-wisemultiple regression. For each voxel, we model the mean E(yr) =∑Tt=1 xrtβtor E(y) = Xβ, where in our example r = 1, . . . , R and T = 6 given by themagnetic field gradients. The matrix X has dimension R×T , where R ≥ T ,with normal equations X ′Xβˆ = X ′y.With a line of I voxels, the y vectors can be placed into a matrix Y ofdimension R × I, as well as the coefficients into a matrix A of dimensionT×I. Thus Y = XA+E, where E is a matrix of random (normal errors). Toensure smoothness on the rows of A, we project each row onto a (common)B-spline basis, and hence A′ = BΓ, where B (I × K) of degree q is builtalong the index 1, . . . , I, and Γ is of dimension K×T . Thus we have E(Y ) =XΓ′B′, which shows directly the tensor structure of X and B. There is onlyone set of penalties on Γ: one penalty for each row. Thus the penalty termin (5) is greatly simplified, as it only involves sums over t and k.Eilers, Marx, and Heim 5With 2-dimensional brain slices, instead of voxels, we conventionally thinkof pixels indexed by i = 1, . . . , I and j = 1, . . . , J . Standard matrix notationnow breaks down, as Y and Γ become 3-dimensional arrays. The R imagesat each pixel are denoted by Y = [yrij ] r = 1, . . . , R and the gradientinformation by X = [xrt], t = 1, . . . , T . If we arrange Y as a matrix withR rows and IJ columns (each column corresponding to one pixel), we canwrite Y = XA + E, with a T × IJ coefficient matrix A = [αtij ] and Eis a R × IJ matrix of random (normal) errors. In this format, pixel-wiseregression would simply lead to the normal equations X ′XAˆ = X ′Y . Weagain impose spatial smoothness of the elements of A by projecting themonto (K × L) tensor products of B-splines. Let B = [bik] and B˜ = [b˜jl](k = 1, . . . ,K; l = 1, . . . , L) be B-spline bases, of degree q, along the twoaxes of the image. Thus αtij =∑k∑l bik b˜jlγklt. For simplicity we use thesame tensor product B-spline basis for each t. This leads toµrij = E(yrij) =T∑t=1xrtK∑k=1L∑l=1bik b˜jlγtkl, (7)in which, by virtue of space-invariant x, we recognize a 3-dimensional tensorproduct model. Again, we ensure further smoothness through additionaldifference penalties on the coefficients of the tensor products. The simplifiedversion of (5) is written asPenalty =T∑t=1{λ1K∑k=1L∑l=1(∆d1γtkl)2 + λ2K∑k=1L∑l=1(∆d2γtkl)2}. (8)Here ∆d1 (∆d2) denotes the d-th order differences along a row (down a col-umn) of the T ×K × L array of tensor coefficients, Γ = [γtkl].Eilers et al. (2005) showed how to fit a multidimensional tensor productmodel with penalties efficiently. Their method adapts to any number of di-mensions, and thus can be applied to the original 3-dimensional DT varyingcoefficient model. Lack of space does not allow us to go into the details of3-dimensional extensions, but can be found in Eilers et al. (2005). In thissetting, both Y and Γ will become 4-dimensional arrays, and the B-splinebasis will come in a triple: B, B˜, B˘. The mean of Y can again be expressedin a tensor form, E(Y ) = XΓ?B?, where Γ? and B? are both carefullyarranged matrices involving Kronecker products.5 Illustrative exampleThe example consists of the 128×128 middle axial slice of a brain volume ofdiffusion tensor images. First we consider 1-dimensional VCM models. Fig-ure 1 (top panel) shows the response signals for the 75th row of the brainslice and the varying coefficients for the diagonal tensor elements (mid-dle panel) and the off-diagonal elements (bottom panel). Figure 2 shows6 Varying Coefficient Tensor Models Signal40 60 80 1000.00.0020.004Signals at six gradients Signal40 60 80 100-0.0010.001Smooth regression: diagonal tensor products Signal40 60 80 100-0.00060.00.0006Smooth regression: off-diagonal tensor products FIGURE 1. Signals at six gradients and smoothed DT coefficients.the corresponding coefficients using pixel-wise multiple regression. Figure3 (left, top and bottom) presents example 32×32 images from the first andthird gradients. The center panels present the diagonal D1(s) (center) andthe off-diagonal D4(s) element pixel-wise. The right panels presents thecorresponding DT elements computed with VCM models and interpolatedto a 128 × 128 resolution; thus demonstrating a benefit from an analyticP-spline approximation.6 DetailsIn all cases, cubic B-splines, and a second order difference penalty, withthe penalty parameter set to 0.01, were used. We used a knot at everysecond pixel for both the 1D and 2D cases (in the latter we analyzed onlya 32 × 32 image). Unfortunately, the demands on computer memory andcomputation time load increase sharply with knot density. For 3D imagesof the original size it will essentially be impossible to explicitly form andsolve the penalized normal equation. We would have to deal with about50 by 50 by 10 by 6 knots, or 1.5 105 equations. We will research directiterative methods to circumvent theses difficulties. Our approach naturallyaccommodates observation weights, so that missing data does not presenta problem, and adjustments can be made for voxels with uncertainty, e.g.Eilers, Marx, and Heim 740 60 80 1000.00.0010.003Pixel-wise regression: diagonal tensor elements40 60 80 100-0.00060.00.0004Pixel-wise regression: off-diagonal tensor elementsFIGURE 2. Estimated pixel-wise DT coefficients.at the boundaries. Further as parameter estimation is grounded in classicalmethods, model summaries, such as: effective dimension and informationcriteria, and delete-one diagnostics, such as cross-validation (CV) are eas-ily obtained. Accessible CV can also be a practical matter when searchingfor the “optimal” penalty tuning parameters, e.g. when isotropy is relaxed.Anisotropic penalization can be achieved by placing different λs on theindexing axes of Γ, leading to a higher dimensional grid-search. We arecurrently researching methods to optimally choose the penalty parameters(λ), including an E-M algorithm that alternates between: a) estimatingthe residual variance and the var(γˆt) to update penalty parameters, andb) estimating Γ with these updated penalty parameters. Since we are ul-timately interested in plotting the smooth coefficients, we are also investi-gating choosing λ based on mean square error of the coefficient estimates,as an alternative prediction oriented criteria. Further, our penalized VCMapproach can be transplanted into a generalized linear model framework,e.g. smoothing coefficients associated with a binary or threshold response.ReferencesBasser, P.J and Jones, D.K. (2002). Diffusion-tensor MRI: theory, experi-mental design and data analysis – a technical review. NMR Biomed,15, 456-467.8 Varying Coefficient Tensor ModelsImage at gradient 1Image at gradient 3Smooth interpolated DT element 1,1Smooth interpolated DT element 1,2Pixel−wise DT element 1,1Pixel−wise DT element 1,2FIGURE 3. Two 32 by 32 pixel example images (left, top and bottom), two DTelements, as computed by pixel-wise regression (middle, top and bottom) and thesame DT elements as computed with VCM model, and interpolated to a 128 by128 grid (right, top and bottom).Eilers, P.H.C, Currie, I., and Durba`n, M. (2005). Low memory, high speedsmoothing on large multidimensional grids. Computational Statisticsand Data Analysis. In Press.Eilers, P.H.C. and Marx, B.D. (2002). Generalized linear additive smoothstructures. J. of Computat. and Graphical Statistics, 11(4), 758-783.Hastie, T. and Tibshirani, R. (1993). Varying-Coefficient Models (with dis-cussion). Journal of the Royal Statistical Society B 55: 757-796.Heim, S. Hahn, K., Auer, D.P., Fahrmeir, L. (2004). Proceedings of the 19thIWSM, Florence, Italy.Papadakis, N.G., Xing, D. Huang, CL-H, Hall, L., and Carpenter, T.A. (1999).A comparative study of acquisition schemes for diffusion tensor imag-ing using MRI. J Magn Reson, 137, 67-82.

A comparative study of acquisition schemes for diusion tensor imaging using MRI.

Diusion-tensor MRI: theory, experimental design and data analysis { a technical review.

Generalized linear additive smooth structures.

Low memory, high speed smoothing on large multidimensional grids. Computational Statistics and Data Analysis.

Varying-Coecient Models (with discussion).

EconStor (ZBW Kiel)

Varying coefficient tensor models for brain imaging

https://epub.ub.uni-muenchen.de/1805/1/paper_436.pdf

Varying Coefficient Tensor Models for Brain Imaging

Abstract

Similar works

Full text

Available Versions

Universität München: Elektronischen Publikationen

Open Access LMU

EconStor (ZBW Kiel)