\chapter{Modeling the spatial and temporal dependence in fMRI data}
\label{ch:STM}

%======================================
\section{Introduction}
%======================================
\label{s:ch2_intro}

%------ Background, general info ----------------
Functional magnetic resonance imaging (fMRI) studies yield large data sets that contain temporal correlations from repeated scanning
(within and between scanning sessions) and complex spatial correlations.
Ignoring correlations affects the precision of estimates of model parameters and consequently may lead to inaccurate
statistical tests. Spatio-temporal modeling may mitigate these shortcomings by incorporating more physiologically
plausible assumptions and by borrowing strength across related measures of neural activity.
%
One session of a typical fMRI neuroactivation
study acquires 3-D scans every 2-3 seconds while the subject performs different experimental tasks.
fMRI studies may also involve  multiple sessions (e.g corresponding to pre- and post-treatment periods).
Typical analyses investigating task-related changes in measured brain activity use a two-stage
procedure in which the first stage involves subject-specific and voxel-specific models relating neural processing
to experimental tasks. The second stage specifies voxel-specific models for group (or population) level parameters.
%------------------------------------------------------
% Existing approaches/current state-of-the-art
%-------------------------------------------------------
The stage I analysis accounts for temporal correlations between serial scans within one session using
variants of autoregressive models \citep{bullmore:06a,purdon:01,friston:02}.
Despite addressing scan-to-scan correlations within a session, the stage I
estimates for multiple sessions may also exhibit correlations, typically not accounted for
at the second stage. Our focus is on modeling these latter correlations.

Some attempts have been made at modeling the repeated measures correlations and spatial correlations in fMRI
data, though typically not in the same statistical model.
\cite{worsley:02} introduce a random effects analysis for combining sessions and the method removes trends in
effects over time.
\citet{bowman:03} fit repeated measures covariance structures within a linear model
to address between-session correlations in positron emission tomography (PET) neuroimaging data.
Neither of these two approaches addresses the spatial correlations between localized brain activity
measurements, and the method by \citet{bowman:03} involves long computations for PET data, making it less suitable
for higher dimensional fMRI data.
%----- Spatial work -------
\cite{bowman:05} presents a two-stage model %for the estimation and testing of localized brain activity
in which the second stage accounts for spatial dependencies within brain regions or networks, and
\citet{bowman:08} give a more flexible Bayesian model to capture correlations both within and between
brain regions. Neither of these approaches, however, accounts for temporal or repeated measures
associations between the multiple experimental effects for a subject.
%
%-- simultaneous spatio-temporal modeling --------
There have been a few attempts to model both temporal and spatial correlations simultaneously, but these approaches
are generally hampered by costly computation. %, making them not practical for routine applications for fMRI data.
The computations generally involve inverting large covariance matrices, and separable models (over space
and time) are often considered \citep{benali:97,hartvig:02}. %,gelfand:01}.
\cite{gossl:01} use a Bayesian approach to fit a semi-parametric spatial and temporal model for fMRI data.
Their approach requires extensive computations, limiting its applicability in practice.
%
\cite{woolrich:04} propose a Bayesian framework to model the noise via a non-separable
space-time, simultaneous autoregressive model. This approach is very time consuming for the fMRI
datasets, taking roughly 6 hours for processing a single slice of 3-D fMRI data.

%
We propose a two-stage model that accounts for both spatial and temporal correlations in fMRI data, and our model leads
to fast parameter estimation. In the second stage, we construct a simultaneous autoregressive model to
capture spatio-temporal correlations between the multiple (session)
effects at a given location and between pairs of voxels within defined anatomical regions.
We use maximum likelihood (ML) methods to estimate parameters from our spatio-temporal model.
We overcome computational challenges involved with estimation by deriving an algorithm that
simplifies the calculations of inverses and determinants of large matrices, leading to fast estimation of the model parameters.
Our spatio-temporal model provides a unified
framework for both voxel-level and region-level inferences.
Using our proposed model, we analyze fMRI data from a study of
inhibitory control in cocaine addicts to evaluate the effects of behavioral therapy on
neural processing related to inhibitory control.
To further delineate benefits of our spatio-temporal model, we conduct two simulation studies: one to validate
the accuracy of our estimation methods and the other to evaluate the relative efficiency of our proposed
spatio-temporal model compared to the general linear model $\bigl($GLM (not to be confused with {\it generalized linear models} as in \cite{mcculagh:89}\,)$\bigr)$.

%======================================
\section{Experimental Data}
%======================================
\label{s:experimental_data}
To illustrate the use of our proposed model and to give motivation for its
development, we apply it to data from an fMRI study evaluating the impact of cocaine addiction and
treatment-related abstinence on neural responses to motor inhibition tasks.
The sample included 12 cocaine addicts enrolled in an intensive outpatient behavioral treatment program and 15 healthy controls,
matched by age, sex (all males), race, handedness, education and early life adversity.
All the subjects were scanned while performing an inhibitory control task used to model
a characteristic hallmark of drug addiction: an impairment in the ability to inhibit behaviors (e.g. drug-seeking behaviors).
The cocaine addicts were scanned in two separate sessions, before and after treatment, and similarly control subjects had baseline and follow-up scans.

{\it Experimental conditions:} The inhibitory control task, referred to as the {\it STOP-signal} task,
was designed to evaluate the ability to cancel a prepotent motor response \citep{aron:06}.
In this task, subjects are presented with visual {\it GO-stimuli}, consisting of uppercase alphabetical letters
appearing on the screen for 0.5 seconds
with an inter-stimulus interval of 2.3 seconds. They are instructed to respond to this signal by pressing a button as quickly as possible.
An auditory STOP-signal (also lasting 0.5 seconds) appears randomly in 16\% of the trials.
The occurrence of the STOP-signal following a GO-stimulus is an indicator to refrain from pushing the button (the prepotent response).
Hence, a successful performance requires the inhibition of a prepotent behavior.
The aim of our analysis is to identify brain locations where there are inhibitory control-related increases in brain activity
following treatment for cocaine addiction.


%======================================
\section{Methods}
%======================================
\label{s:ch2_methods}
Data from a single scan in our study are represented as a $53 \times 63 \times 46$
3-D rectangular lattice, comprised of a large number of voxels,
indexed by a mapping of coordinates $(x,y,z)$ to $v=1,\ldots,V$ ($V> 150,000$ in our study).
Each voxel contains intensity information corresponding to measures of localized brain activity
(BOLD fMRI responses).
Since the number of intracranial voxels is too large to estimate a global
correlation matrix including all voxel pairs,
we consider an approach that partitions the voxels into mutually exclusive subregions
 based on known neuroanatomical divisions.
Specifically, we define anatomical regions using a parcellation of the brain based on Brodmann areas (BA)
\citep{brodmann:09} and additionally separate BA's in the left and right hemispheres.
To provide broader coverage, we also add several subcortical regions (left and right caudate, thalamus and putamen) to our parcellation,
which were obtained from the Automated Anatomical Labeling (AAL) map based on \citep{tzourio:02}.
Other such parcellations exist, including the Freesurfer software \citep{fischl:02,fischl:04} and
LONI (Laboratory of Neuro Imaging, UCLA) Probabilistic Brain Atlas \citep{shattuck:08}.
Building our model based on well-established anatomical parcellations, rather than data-driven functional clusters
aids interpretability.
%--------------------------------------
\subsection{Statistical Model}
%--------------------------------------
Following the conventional two-stage modeling approach for fMRI data, we fit
a GLM at the first stage for each individual's vector of serial BOLD responses,
separately for each voxel.
The GLM regresses each voxel's BOLD responses for subject $i$
($i=1,\ldots,K$) on within-subject design variables (e.g. session or stimulus
indicators) and on covariates that are not of intrinsic interest, such as high-pass
filtering variables.
The model accounts for short-range scan-to-scan correlations within a session using
a first-order autoregressive process with white noise %\citet{bullmore:06a}, \citet{purdon:01}, or
\citep{purdon:01, friston:02}.
The regression coefficient, $B_{igs}^{(p)}(v)$, for voxel $v$ represents a
summary measure of an individual's (mean) neural activity associated with session
or stimulus $p$.
We add a subscript $g$ to denote the neuroanatomic region to which voxel $v$ belongs
($g=1,\ldots,G$), with the region consisting of $V_g$ voxels in total.

At the second stage, we propose a spatio-temporal (ST) autoregressive model
to capture temporal correlations between the multiple sessions and
spatial correlations between pairs of intra-regional voxels.
We express the spatio-temporal model as follows:
%
\begin{equation}
\label{2nd_stage}
\begin{split}
B_{igs}^{(p)}(v) & = \mathbf{x}_i' \bb_{gs}^{(p)}(v)+\rho_{gs} \; \frac
{1}{V_g-1} \sum_{v^* \in {\cal N}_v } \Bigl(B_{igs}^{(p)}(v^*)
- \mathbf{x}_i' \bb_{gs}^{(p)}(v^*)\Bigr)
\\ &+ \xi_{gs} \; \frac {1}{q-1} \sum_{p^* \ne p}
\Bigl(B_{igs}^{(p^*)}(v) - \mathbf{x}_i' \bb_{gs}^{(p^*)}(v)\Bigr)
+e_{igs}(v)\, ,
\end{split}
\end{equation}
%
where $\mathbf{x}_i$ is a vector of between-subject design variables,
$\bb_{gs}^{(p)}(v)$ are the group-level parameters,
$e_{igs}(v) \sim N(0,\phi_{gs}^2)$ ($\phi_{gs}^2$ is the
between-subject variation), $\rho_{gs}$ is a spatial
dependence parameter, $\xi_{gs}$ reflects the temporal
dependence between sessions, ${\cal N}_v$ is the set
of all voxels within the same neuroanatomic region as $v$,
and $s$ denotes a subject group (e.g., patients, controls).
%
When $\rho_{gs}=\xi_{gs}=0$, model \eqref{2nd_stage} resembles a typical stage II
GLM from a random effects analysis of brain imaging data, with constant
variances within each region.

%--------------------------------------
\subsection{Parametric Covariance Model}
%--------------------------------------
\label{ss:covar_model}
Our model involves three parameters defining
the variance-covariance structure: the spatial ($\rho$) and temporal ($\xi$) dependence parameters
and the variance parameter ($\phi^2$).
In this section, we present details about the assumed covariance/correlation model.
%
Let $\bB_{igs} = \Bigl( \bB_{igs}^{(1)'}, \ldots ,\bB_{igs}^{(q)'}
\Bigr)'$, % denote the vector of individualized coefficients,
where each $\bB_{igs}^{(p)}$  contains
individualized parameters for all voxels in region $g$, and similarly
$\bb_{gs} = \Bigl( \bb_{gs}^{(1)'}, \ldots ,\bb_{gs}^{(q)'}\Bigr)'$.
Model \eqref{2nd_stage} implies
% (see \eqref{e:matrix_form}, Section \ref{ss:matrix_form})
that $\bB_{igs}$  has the following
multivariate normal distribution $\bB_{igs} \sim {\rm MVN}\Bigl(
\bb_{gs}, \mathbf{\Omega}_{gs}  \Bigr)$ (see Appendix \ref{ss:STM_matrix_form}),
where $\mathbf{\Omega}_{gs} =\phi_{gs}^2 \mathbf{\Psi}_{gs} \mathbf{\Psi}_{gs}$,
and $\mathbf{\Psi}_{gs}$ is a $q V_g \times q V_g$ matrix with
%
%--------- Correlation matrix -------------
\begin{equation}
\mathbf{\Psi}_{gs}^{-1}=\left[%
\begin{array}{cccc}
  \I_Vg-\rho_{gs} \Wg & -\frac{\xi_{gs}}{q-1}\I_Vg & \ldots & -\frac{\xi_{gs}}{q-1}\I_Vg \\
  -\frac{\xi_{gs}}{q-1}\I_Vg & \I_Vg-\rho_{gs} \Wg & \ldots & -\frac{\xi_{gs}}{q-1}\I_Vg \\
  \vdots & \vdots & \ddots & \vdots \\
  -\frac{\xi_{gs}}{q-1}\I_Vg & -\frac{\xi_{gs}}{q-1}\I_Vg & \ldots & \I_Vg-\rho_{gs} \Wg \\
\end{array}
\right]_{qV_g \times qV_g}
\label{Psi_inv}
\end{equation}
($\I_Vg$ and $\J_Vg$ denote identity and unity matrices of size $V_g$).
In \eqref{Psi_inv}, $\Wg=[1/(V_g-1)](\bJ_{V_g}-\bI_{V_g})$ specifies the neighborhood structure
consisting of all voxels that fall in the same anatomical region.
The blocks along the main diagonal of $\mathbf{\Psi}_{gs}^{-1}$ allow
for spatial correlations between measures of task-related brain activity for
voxels in region $g$,  i.e. between elements of $\bB_{igs}^{(p)}$.
%
The off-diagonal blocks of $\mathbf{\Psi}_{gs}^{-1}$ capture
correlations between the summary measures of brain activity
associated with various scanning sessions.
%
Our model assumes an exchangeable covariance structure between voxels in the same
neuroanatomic region.
Although the complexity of human brain function is likely to render departures from this assumption,
the exchangeable structure provides an improvement over the often used independence assumption and
seems reasonable for statistical modeling purposes based on descriptive empirical results (not shown here;
see Appendix \ref{ss:exploratory_analysis}).

\vspace{1ex}
{\it Calculating the matrix $\bPsi_{gs}$.} Calculating $\bPsi_{gs}$ involves inverting a matrix of
size $qV_g \times qV_g$, which can become unwieldy for large regions, and
this calculation is performed iteratively during estimation of the model.
We derive an algorithm to facilitate calculations of inverses and determinants of our large
highly structured covariance matrices, thereby enabling estimation of our model for fMRI
applications.
%We represent $\invPsi_{gs}$ as a partitioned matrix with the diagonal
%blocks having the same compound symmetric form, and the off-diagonal blocks being
%proportional to the identity matrix $\bI_{V_g}$.
By recursive calculation on $q$, the number of blocks in $\invPsi_{gs}$, we represent $\bPsi_{gs}$ as
follows

\begin{equation}
\bPsi_{gs} = \begin{bmatrix}
           d\I_Vg + f\J_Vg & u\I_Vg + z\J_Vg & \cdots & u\I_Vg + z\J_Vg  \\
           u\I_Vg + z\J_Vg & d\I_Vg + f\J_Vg  & \cdots & u\I_Vg + z\J_Vg \\
           \vdots & \vdots & \vdots & \ddots \\
           u\I_Vg + z\J_Vg & u\I_Vg + z\J_Vg & \cdots & d\I_Vg + f\J_Vg
         \end{bmatrix}_{q V_g\times q V_g} \, .
\label{f:Psi}
\end{equation}

\vspace{2ex}
\noindent
Hence, $\bPsi_{gs}$ is determined by four functions: $d(\rho_{gs},\xi_{gs},V_g,q)$,
$f(\rho_{gs},\xi_{gs},V_g,q)$, \\
$u(\rho_{gs},\xi_{gs},V_g,q)$, and $z(\rho_{gs},\xi_{gs},V_g,q)$.
The explicit formulas for $d$, $f$, $u$ and $z$ are quite lengthy and are given in the Appendix \ref{ss:coefficients_of_psi}.
We employ \eqref{f:Psi} to circumvent issues of long computing times and excessive memory consumption.

Estimated spatial and temporal parameters are not directly interpretable as
spatial and temporal correlations. The spatial and temporal correlations,
denoted $\eta_\rho$ and $\eta_\xi$, respectively, require calculation of
the covariance matrix $\mathbf{\Omega}_{gs}$, i.e. of
$\mathbf{\Psi}_{gs} \mathbf{\Psi}_{gs}$.
Leveraging the special structure of $\bPsi_{gs}$, we obtain $\mathbf{\Psi}_{gs} \mathbf{\Psi}_{gs}$
as a partitioned matrix, with equal diagonal blocks $A \bI_{V_g} + B \bJ_{V_g}$, and
off-diagonal blocks $C \bI_{V_g} + D  \bJ_{V_g}$,
where
\begin{align}\label{PsiPsi}
A & = d^2+(q-1)u^2 & B & = 2df+V_g f^2+(q-1)(2uz+V_g z^2) \notag \\
C & = 2du+(q-2)u^2 & D & = 2(fu+dz+V_g fz)+ (q-2) (2uz+V_g z^2)\, .
\end{align}
The expressions for the spatial correlation, $\eta_\rho$, and the temporal correlation, $\eta_\xi$, are then given by
$\eta_\rho=B/(A+B)$ and $\eta_\xi=(C+D)/(A+B)$.
%--------------------------------------
\subsection{Estimation}
%--------------------------------------
\label{ss:estimation}
We perform estimation of the second stage spatio-temporal model using ML methods.
The log-likelihood function for model \eqref{2nd_stage} is given by
\begin{eqnarray}
\label{log-likelihood}
{\cal L} (\boldsymbol{\phi},\boldsymbol{\rho},\boldsymbol{\xi},\bb
| \bB) \propto \sum_{s=1}^2 \sum_{i \in K_s} \sum_{g=1}^G \biggl\{  -\tfrac{q
V_g}{2} \ln (\phi_{gs}^2) + \ln \bigl(\text{abs}|\mathbf{\Psi}_{gs}
^{-1}| \bigr) \notag \\                                    % NOTE: use \notag
- \tfrac{1}{2\phi_{gs}^2}  (\mathbf{B}_{igs} - \mathbf{X}_i
\beta_{gs})' \mathbf{\Psi}_{gs}^{-1} \mathbf{\Psi}_{gs}^{-1}
(\mathbf{B}_{igs} - \mathbf{X}_i \beta_{gs}) \biggr\}\, ,
\end{eqnarray}
where $K_s$ denotes the number of subjects in subgroup $s$ (cocaine addicts and controls).
%
The ML estimator of the mean parameter vector
$\bb_{gs}$ is given by
$
\hat{\bb}_{gs}=\frac{1}{K_s}\sum\limits_{i=1}^{K_s}
\boldsymbol{B}_{igs}, \;  g=1,\ldots,G
$.
$\hat{\bb}_{gs}$ is unbiased (as is $\hat{\bb}_{gs}$ from a GLM analysis)
and does not depend on the covariance parameters.
Estimation of the covariance parameters, therefore, proceeds using
the partially maximized likelihood function
${\cal L} (\boldsymbol{\phi},\boldsymbol{\rho},\boldsymbol{\xi},\hat{\bb}
| \bB)$.
We use the Fisher scoring algorithm to estimate the covariance parameters
$\boldsymbol{\rho}_s=(\rho_{1s},\ldots,\rho_{Gs})'$, $\boldsymbol{\xi}_s=(\xi_{1s},\ldots,\xi_{Gs})'$, and
$\boldsymbol{\phi}_g=(\phi_{1s},\ldots,\phi_{Gs})'$.
%Fisher scoring is attractive here since it simplifies computations by replacing the Hessian
%matrix by its expected value and since it is more robust to poor starting values than the
%Newton-Raphon's algorithm \citep{jennrich:86}.

\vspace{1ex}
{\it Note on the parameter space.} % ($\rho$ and $\xi$)}
To apply our model, the matrix $\mathbf{\Psi}_{gs}^{-1}$
must be nonsingular, or equivalently $\mathbf{\Omega}_{gs}$ must be positive definite.
In the context of our fMRI data,  with $V_g>20$, for all $g$, and $q=2$,
the resulting explicit parameter constraints are as follows:
%
\begin{align*}\label{constraints}
\rho_{gs} + (V_g -1)(1+\xi_{gs}) & \ne 0 & \rho_{gs} + (V_g -1)(1-\xi_{gs}) & \ne 0 \\
\rho_{gs} - \xi_{gs} & \ne 1 & \rho_{gs} + \xi_{gs} & \ne 1\, .
\end{align*}
%
We monitor these boundary constraints during our iterative estimation procedure.

Estimation involves iterative calculations of the determinant and the inverse of a
$q V_g \times q V_g$ matrix (up to roughly 5,800$\times$5,800 in our data application).
Employing the simplifications discussed in Section \ref{ss:covar_model}
enables fast and efficient computations.
Appendix \ref{ss:score_fcts} contains details regarding the score functions necessary for estimation.

%--------------------------------------
\subsection{Inferences}
%--------------------------------------
\label{ss:inferences}
Inferences that are commonly sought in functional neuroimaging studies target the identification of differences in neural processing between experimental tasks, sessions (e.g. treatment periods), or groups.  One may pursue these inferences at an extremely localized (voxel) level or at a more spatially coarse regional level.  In either case, estimation and inferences in this setting produce maps of distributed brain activity and the corresponding thresholded maps showing statistically significant (or highly probable, in a Bayesian context) differences.
We seek inferences about linear functions of $\bb_{gs}^{(p)}$, where for the inhibitory control study in cocaine addicts,
$p=1$ (baseline or pre-treatment) or $p=2$ (follow-up or post-treatment) and $s=1$ (cocaine addicts) or $s=2$ (controls).

{\it Voxel-level inferences.}
Following estimation using our spatio-temporal model, we obtain $t$-statistic images from the voxel-specific contrast estimates,
then threshold the $t$-statistics to determine voxels exhibiting statistically significant changes (or differences) in
brain activity. To adjust for typical changes observed with repeated scanning sessions in the inhibitory control study,
 we identify voxels for which the changes in brain activity
following treatment are larger in cocaine addicts than the corresponding changes between follow-up and baseline activity in
control subjects.
%
Specifically, we estimate and test hypotheses about
$\boldsymbol{\th}_g^* = \mathbf{C}^*_g \boldsymbol{\be}_g $, where
$\mathbf{C}_g^* = [-\bI_{V_g} \,\, \bI_{V_g} \,\, \bI_{V_g} \,\, -\bI_{V_g}]_{(V_g \times 4V_g)}$ and
$\bb_g=[\bb_{g1}^{(1)'}, \bb_{g1}^{(2)'}, \bb_{g2}^{(1)'}, \bb_{g2}^{(2)'}]'_{(4V_g \times 1)}$.
Each element of $\boldsymbol{\th}_g^*$ represents a voxel-specific parameter.
%
We construct Wald-type statistics for hypothesis testing, with
$\text{Var}(\hat{\boldsymbol{\th}}_g^*) = \mathbf{C}_g^* \, \text{Var}(\hat{\boldsymbol{\be}}_g )\, \mathbf{C}_g^{* '}$
and $df=2(n_p+n_c-2)$, where $n_p$ is the number of addicts, and $n_c$ is the number of controls.
%
We apply the stringent threshold of $\al=0.005$, which is one approach that has
a strong precedent in the neuroimaging literature, but other thresholding approaches, such as false discovery rate \citep{benjamini:95}
and random field theory \citep{friston:95a} are also available.

{\it Region-level inferences.}
Our spatio-temporal model also enables analyses targeting
an entire anatomical brain region corresponding to the underlying
anatomical parcellation. %(in our case Brodmann areas).
%
For inferences in region $g$, define $\th_g=\mathbf{C}_g \bb_g$, with
$\mathbf{C}_g=\tfrac{1}{V_g} [-\mathbf{1}_{V_g}' \,\, \mathbf{1}_{V_g}' \,\, \mathbf{1}_{V_g}' \,\, -\mathbf{1}_{V_g}']_{(1 \times 4 V_g)}$,
and construct the $t$-statistic using
$\text{Var} (\hat{\th}_g) = \mathbf{C}_g \, \text{Var} (\hat{\bb}_g) \, \mathbf{C}_g' $
and degrees of freedom  $df_r = 2V_g(n_p+n_c-2)$.
Combining both sessions for subgroups $s$, i.e.,
$\bb_{gs}=[\bb_{gs}^{(1)'}, \bb_{gs}^{(2)'}]'$, the variance of $\hat{\bb}_{gs}$ is given by
$\text{Var}(\hat{\bb}_{gs}) = \frac{1}{n_p+n_c} \phi_{gs}^2 \bPsi_{gs} \bPsi_{gs}$.
A notable advantage of modeling spatial correlations using our model is that the resulting regional
level analyses account for spatial dependencies between intra-regional voxel pairs, rather than implicitly assuming independence.
%

%{\it Inferences about $\rho_{gs}$ and $\xi_{gs}$. }. Inferences about spatial and temporal estimates are of interest.
{\it Statistical significance of spatial and temporal correlations.}
It is often informative to examine the magnitudes of spatial and temporal correlations as discussed in
Section 3.2. % and in some cases to
% make inferences about these correlations.
Given the complexity of the estimators of these quantities,
it is difficult to derive analytical expressions for the variances of these estimators.
If one has inherent interest in testing hypotheses about these parameters, then standard bootstrapping methods are
applicable. However, this approach may be computationally demanding. %, for example taking 100's % or even 1,000's
%of hours (possibly reduced by using multiple processors on a network).
Given that testing hypotheses about correlations was not of chief interest in the inhibitory control study,
we present estimates of correlations mainly for descriptive purposes.
%We provide bootstrap standard errors of the estimated correlations.
We estimate standard errors of the correlations using 30 bootstrap samples drawn (with replacement) from
the subjects in our data, separately for controls and patients. We apply our model and estimate the spatial and temporal
correlations for each sample, and then calculate the bootstrap standard errors as the standard
deviations of the 30 estimated values for each parameter, which should provide reasonable estimates
of variability \citep{efron:98}.

%====================================================================
\section{Application to Inhibitory Control in Cocaine Addicts Study}
%====================================================================
\label{s:application}
We used our spatio-temporal model to analyze fMRI data from the inhibitory control study in cocaine addicts.
Spatial and temporal parameters were estimated separately for cocaine addicts and controls.
The aim of our analysis was to evaluate whether the inhibitory control-related changes
in measured brain activity following treatment for the addicts were larger than the corresponding
between-session changes in the control subjects.

All of the image pre-processing (slice timing, realigning, normalizing, and smoothing %(3mm,6mm, or 9mm FWHM Gaussian kernel),
[6mm (2 voxels) full width at half maximum (FWHM) Gaussian kernel],
as well as the first level, single subject analyses, were carried out in SPM5\footnote{SPM5 is
is a MATLAB software package implementing Statistical Parametric Mapping for neuroimaging data
available for download from the Wellcome Trust Centre for Neuroimaging web page
{\tt http://www.fil.ion.ucl.ac.uk/spm/software/spm5/}.}.
To account for inter-subject neuroanatomic differences that may persist after spatial normalization,
we chose to apply very focal spatial smoothing, rather than the more spatially expansive smoothing often applied
in practice. Our smoothing should only have a small effect on subsequent spatial modeling, since most of the
Brodmann areas (within which we model spatial correlations) are quite large, relative to the
small size of the smoothing kernel.
The second level analysis, implementing our proposed spatio-temporal model,
was carried out using MATLAB.
Our analysis includes 38 BA's in each hemisphere, after excluding those with fewer than 20 voxels,
and we supplement the BA's with six subcortical regions, resulting in a parcellation consisting of 82 brain regions.
%We used an anatomical parcellation consisting of 82 cortical regions
%(after excluding those with fewer than 20 voxels) that was based on Brodmann areas, additionally separating regions falling in
%the left and the right hemispheres and supplemented with some subcortical regions.
We perform exploratory analyses, computing crude estimates of spatial and temporal correlations (Appendix \ref{ss:exploratory_analysis}), which
provide support for the existence of correlations from both sources as well as for the exchangeability assumption in our covariance model.


%----------------------------------
\subsection{Voxel-level inferences}
%----------------------------------
\label{ss:voxel_lev_inf}
Figure \ref{f:voxel_lev_inf}(a) shows voxels that reveal statistically significant increases in inhibitory control-related
brain activity following treatment at $\alpha=0.005$.
The axial slices correspond to 5 mm, 14 mm, 44 mm, and 65 mm above the anterior–-posterior commissural plane.
%
The significant areas include the right frontopolar cortex (BA 10), left middle temporal gyrus (BA 21), and
retrosubicular area (BA 48) in the 5 mm slice;  left and right thalamus and the right inferior frontal cortex (BA 45) at 14 mm;
left visual association cortex V3 (BA 19) and somatosensory association cortex (BA 7) and the left angular gyrus (BA 39) at 44 mm;
and the right pre-supplementary motor area (pre-SMA, BA 6) at 65 mm.
Other regions (not shown) exhibiting significant treatment-related
increases in brain activity include the right inferior prefrontal gyrus (BA 47),
left and right fusiforum gyrus (BA 37), and left primary auditory cortex (BA 41).

Figure \ref{f:voxel_lev_inf}(b) shows results from the corresponding GLM--based analysis
commonly used in the neuroimaging literature. Many of the areas
showing statistical significance agree with the results from our approach (e.g. Brodmann areas 6, 7, 10, and 19
% (e.g. Brodmann areas 10, 11, 19, 37 and 47)
emerged in both analyses), but there are some differences. For example,
our method identifies increased treatment-related brain activity in BA 48 and BA's 21 and 22 (5 mm),
which do not emerge from the GLM analysis. Also, both methods identify areas in the thalamus, but our method detects
a large right thalamic activation not revealed by the GLM. Similarly, both methods identify
voxels in BA's 6 (pre-SMA) and 39, but our model yields more spatially extensive activations.
Overall, our ST model produces more statistically significant voxels than the GLM.
However, neither method produces significant voxels when applying more conservative family-wise error or false
discovery rate multiple testing procedures.
The differences between the two models become more apparent at larger significance levels (e.g. $\al=0.01$, not shown),
suggesting that with a larger sample size and stronger activation signal, greater differences are likely to emerge.
We provide interpretive remarks about our results in Section \ref{ss:implications}.

\begin{figure}[!h]
%\begin{center}
\begin{tabular*}{\textwidth}{@{}c@{\extracolsep{\fill}}c@{\extracolsep{\fill}}c@{\extracolsep{\fill}}}
\hline
 \underline{Model} & Voxel-level inferences & Region-level \\
        &                        & inferences \\
ST  &
\includegraphics[scale=0.52]{Figures/STM_vox_lev_new_w_labels.eps} &
\, \includegraphics[scale=0.52]{Figures/STM_reg_lev_mask_label.eps} \\
 & (a) & (c) \\ [2ex]
GLM &
\includegraphics[scale=0.52]{Figures/GLM_vox_lev_new_w_labels.eps} &
\, \includegraphics[scale=0.52]{Figures/GLM_reg_lev_inf_0p05.eps} \\
 & (b) & (d)
\end{tabular*}
\caption[Voxel-level and region-level inferences for two models: the proposed spatio-temporal model and the GLM]
{Voxel-level and region-level inferences for two models: the proposed spatio-temporal (ST) model (top) and the
 GLM (bottom). (a): voxels that achieve statistical significance in the ST model analysis at $\al=0.005$ in favor
of a one-sided alternative hypothesis of increased activity following treatment for cocaine addiction.
Significant increases occur in the right frontopolar area (BA 10), left middle and superior temporal gyri (BA's 21/22),
right retrosubicular area (BA 48), right inferior frontal cortex (BA 45), left and right thalamus (Thal),
left somatosensory and visual association cortices (BA's 7 and 19), left angular gyrus (BA 39), and right BA 6 (pre-SMA);
(b): voxel-level results from a GLM analysis;
%
(c): region, left BA 39 (angular gyrus), that achieves significance in our ST model analysis at $\al=0.05$;
(d): corresponding region-level results from the GLM.
Slice labels denote the distance, in mm, above the anterior-posterior commissural plane.}

\label{f:voxel_lev_inf}
%\end{center}
\end{figure}


%------------------------------------
\subsection{Region-level inferences}
%------------------------------------
\label{ss:reg_lev_inf}
Our region-level analysis of treatment-emergent changes in neural processing,
using methods described in Section \ref{ss:inferences},
revealed only one region (left BA 39, angular gyrus)
achieving statistical significance at $\al=0.05$, uncorrected (see Figure \ref{f:voxel_lev_inf}(c)).
The comparative region-level GLM analysis (that implicitly assumes spatial and temporal independence)
did not yield any statistically significant differences.
Neither method reveals significant differences at more stringent thresholds such as
$\al=0.005$ (uncorrected).

%--------------------------------------------------------------
\subsection{Spatial and temporal correlations}
%--------------------------------------------------------------
\label{ss:correlations}
Figure \ref{f:correlations} shows the model-based estimates of spatial and temporal correlations for
both the cocaine addicts and control subjects in one axial slice of the brain.
The color scale indicates the strength of correlations, with brighter shades indicating correlations of
larger magnitude. For temporal correlations, positive values are shown in one color (top colorbar) and negative values
are shown in a separate color (bottom colorbar). For clearer distinction, we add stripes to
regions with negative correlations.
All of the correlations in Figure \ref{f:correlations} would have the value zero in a GLM analysis, signaling the
need to account for them in our model.

%===== if using EPS figures =========
\begin{figure}
\begin{center}
\begin{tabular}{ccc}
& Cocaine addicts & Controls \\
Spatial & \includegraphics[scale=1]{Figures/Pats_spatial_sl24_w_labels.eps} &
\includegraphics[scale=1]{Figures/Ctrls_spatial_sl24_w_labels.eps} \\
& (a) & (b) \\ [2ex]
Temporal & \includegraphics[scale=1]{Figures/Pats_temporal_sl20_w_labels_striped.eps} &
\includegraphics[scale=1]{Figures/Ctrls_temporal_sl20_w_labels_striped.eps} \\
& (c) & (d)
\end{tabular}
\end{center}
\caption[Model-based estimates of the spatial and temporal correlations for the inhibitory control study]
{Model-based estimates of the spatial and temporal correlations for the inhibitory control study
(a) addicts: spatial correlations, (b) controls: spatial correlations, (c) addicts: temporal correlations,
(d) controls: temporal correlations. For (a) and (b), the darker to lighter intensities in the
colorbar represent lower to higher correlations. For (c) and (d), the top colorbar represents positive
correlations, again with lighter intensities reflecting higher correlations, and the bottom
colorbar represents negative correlations (stripes added to maps) with brighter intensities
indicating stronger negative correlations.
%
The temporal correlations for control subjects range from -0.28 to 0.64 and are in general larger
than those for the cocaine addicts, which range from -0.24 to 0.45. The regional spatial
correlations are not consistently higher for either group, ranging from 0.08 to 0.54
among cocaine addicts and from 0.04 to 0.65 for healthy controls.}
\label{f:correlations}
\end{figure}



\begin{table}
\centering
 % \def\~{\hphantom{0}}
% \begin{minipage}{175mm}
\caption[Model-based estimates of spatial and temporal correlations
and estimated bootstrap standard errors]
{Model-based estimates of spatial and temporal correlations
and estimated bootstrap standard errors (in parentheses).}
\label{t:table1}
\begin{tabular*}{\textwidth}{@{}l@{\extracolsep{\fill}}l@{\extracolsep{\fill}}l@{\extracolsep{\fill}}l@{\extracolsep{\fill}}l@{\extracolsep{\fill}}}
\hline \hline
 Group & Region &  Temporal               & Region & Spatial \\
       &        & correlations (std. err.) &        & correlations (std. err.) \\
\hline\hline
\underline{Addicts} &  & & & \\ [2pt]
 & BA 43 (L) & 0.45 (0.19) &  BA 43 (L) & 0.54 (0.10) \\
 & BA 22 (R) & 0.39 (0.15) &  BA 29 (L) & 0.50 (0.12) \\
 & BA 42 (R) & 0.37 (0.18) & BA 3 (L) & 0.47 (0.10)\\
 % & BA 42 (R) & 0.37 (0.18) & BA 29 (L)$^{*_1}$ & 0.50 (0.12)\\
 & BA 3 (L) & 0.37 (0.17) &  caudate (R) & 0.47  (0.10) \\
 & BA 9 (R) & 0.34 (0.28) &  BA 43 (R) & 0.43 (0.07) \\ \hdashline %\hline
 & BA 35 (R) & -0.24 (0.06) &  BA 48 (R) & 0.08  (0.07) \\
\hline\hline
\underline{Controls} &  & & & \\ [2pt]
 & BA 22 (L) & 0.64 (0.07) & BA 26 (L) & 0.65 (0.05)\\
 & BA 42 (R) & 0.63 (0.09) & BA 26 (R) & 0.60 (0.08) \\
 & BA 22 (R) & 0.60 (0.09) & BA 5 (R) & 0.50 (0.12) \\
 & BA 21 (L) & 0.58 (0.10) & putamen (L,R) & 0.50 (0.08, 0.06)\\
 & BA 43 (R) & 0.46 (0.10) & BA 38 (L) & 0.47 (0.10) \\ \hdashline
 & BA 26 (R) & -0.28 (0.11) & BA 48 (R) & 0.04 (0.02)\\
\hline\hline
\end{tabular*}
%\end{minipage}
%\vspace*{-6pt}
%{\small $^{*_1}$  Associated temporal dependence parameter ($\xi$)
%achieved statistical significance at level $\al=0.05$. }
\end{table}

Table \ref{t:table1} gives the model-based estimates, along with bootstrap standard errors (in
parentheses), of both temporal and spatial correlations for select regions in the inhibitory control study.
The associated spatial ($\rho$) and temporal ($\xi$) dependence
parameters achieve statistical significance for all regions at $\al=0.01$, except for BA 29 (L), which \
is significant at $\al=0.05$.
%
The selected regions in Table 1 are those that have correlations with the largest magnitude.
The temporal correlations for cocaine addicts range from -0.24 to 0.45 and tend to be larger in control subjects
for whom they range from -0.28 to 0.64. The regional spatial correlations are not consistently higher for either group,
ranging from 0.08 to 0.54 among cocaine addicts and from 0.04 to 0.65 for healthy controls.

The brain regions with highest temporal correlations differ between the two groups. For the controls, the
highest temporal correlations ($\sim 0.6$) are in right BA 42, partially covering the auditory association cortex,
and in left BA 21, left BA 22, and right BA 22, regions in the temporal lobe which are also
revealed by our voxel-level analysis.
For the cocaine addicts, the highest temporal correlations are somewhat smaller ($\sim 0.4$) and appear in left
BA 43 and in right BA 22.
Spatial correlations for the control group are highest ($\geq 0.6$) in left and right BA 26, a retrosplenial region in the cingulate cortex.
For the cocaine addicts, the highest estimated spatial correlations ($\sim$ 0.5) are in left BA 29, also in the retrosplenial region,
and in left BA 43.
For both groups, the region with the lowest estimated spatial correlation ($\leq 0.08$) is the
right BA 48, which is the second largest region (2,875 voxels). This suggests that it is the least spatially
homogeneous region and that subdividing large BA's may be warranted.
%The complete lists of the spatial and temporal correlations for each of the 82 regions,
%both for the cocaine addicts and the healthy controls, appear in the Tables A1 and A2 (Appendix A), respectively.

\subsection{Implications of Results}
\label{ss:implications}
%---- OUR FINDINGS -----------------
The study on inhibitory control in cocaine addicts that is the subject of this novel image analysis approach is,
to the best of our knowledge, the first to examine the changes in brain activity associated with
an intensive addiction behavioral therapy.
%
Inhibiting prepotent or automatic responses is critical to the organization of successful goal-directed behaviors.
 \cite{aron:07a} state that it is ``... likely that the simple Stop signal task taps into a control circuit \ldots [and]  variation (or damage) to key nodes in this control circuitry (or to their connections) could produce important individual differences, for example in liability toward and recovery from addiction (\cite{garavan:07}.''
  \citet{aron:07} demonstrated that the ability to stop motor responses depends critically on a network of structures including the right inferior frontal cortex and the subthalamic nucleus (STN), both of which are connected to the pre-SMA.

Our results agree, in several aspects, with the findings from these and other previous neuroimaging studies.
%Clint:
Greater activation of the right orbitofrontal cortex is necessary for behavioral inhibition in impulsive individuals
\citep{horn:03}. Additionally, the right ventral prefrontal cortex has been selectively implicated in the neuropathophysiology of drug addiction \citep{goldstein:02,bolla:03,volkow:05}.
%\citep{goldstein:01, goldstein:02,bolla:03,volkow:05}.
% The number of studies showed that the ability to stop motor responses depends critically on the right inferior frontal cortex (IFC) (\citet{aron:07}, % \cite{floden:06}).
Voxel-level analysis of our data found that the major brain effect of treatment is represented by an enabling of the
right inferior prefrontal gyrus (BA 47) and right orbitofrontal cortex (BA 11) responses to a demand for inhibition of habitual motor responses.
Our findings are promising since they suggest a positive effect of the behavioral therapy on the neural processing deficits associated with
cocaine-addiction, though we do not strictly regard these effects to be causal.

Our analysis was not able to explore voxels in the STN because this brainstem area is not represented in the Brodmann area parcellation.
This limitation is not an inherent shortcoming of our model, and one can easily include additional regions such as the STN in the analysis, if desired.
The pre-SMA is identified as a significant treatment-related area in our voxel-level analysis.
The caudal pre-SMA is critical to controlled action selection \citep{nachev:07} and is functionally compromised in cocaine addicts
\citep{kaufman:03}. Enhancing the neural response of the pre-SMA to a demand for a shift from habitual responses to controlled response inhibition represents a plausible neural correlate of drug refusal skills acquired in behavioral therapies targeting relapse prevention.
Our region level analysis identifies the angular gyrus (BA 39), which is important in visuospatial attention
\citep{cattaneo:08}.

In areas exhibiting significant voxel-level post-treatment increases in task-related brain activity, we observe distinct patterns
of intra-regional correlations between patients and controls.
Based on the estimates of the spatial correlations, the controls exhibit higher spatial coherence in, for example,
right BA's 11 $\bigl($controls = 0.24(0.11), addicts = 0.14(0.04)$\bigr)$  and 47 $\bigl($controls = 0.28(0.13), addicts = 0.18(0.10)$\bigr)$,
which are frontal areas involved in executive functioning, planning and decision making, as well as in both left and right thalamus
$\bigl($e.g. in right thalamus, controls = 0.32(0.09), addicts = 0.17(0.07)$\bigr)$, which play a central role in the flow of information to the cortex \citep{sherman:02}.
Cocaine addicts showed more coherent functioning in left BA's 19 $\bigl($addicts = 0.16(0.02), controls = 0.09(0.07)$\bigr)$ and 39
$\bigl($addicts = 0.31(0.10), controls = 0.16(0.10)$\bigr)$, which are involved in visual attention. The latter finding is
particularly interesting since BA 39 was identified in our region-level analysis.

%======================================
\section{Simulations}
%======================================
\label{s:simulations}
We performed two simulation studies.
%{\it Simulation 1.}
The goal of the first simulation study was to evaluate the accuracy of our estimation procedure.
% {\it Description of the simulated data.}
We selected regions (Brodmann areas) with voxel sizes ranging from 100 to 2000.
To generate data for the simulations, we used our experimental data from the inhibitory control
fMRI study to specify {\it true} values for $\bb_{gs}=(\bb_{gs}^{(1)}, \bb_{gs}^{(2)})$ (where the superscripts
denote sessions 1 and 2), $\rho_{gs}$, $\xi_{gs}$ and
$\phi_{gs}$. We drew 200 samples from the multivariate normal probability distribution  ${\rm MVN}(
\bb_{gs}, \phi_{gs}^2 \mathbf{\Psi}_{gs} \mathbf{\Psi}_{gs})$, %with $\bPsi_{gs} =(\bI -\rho_{gs} \Wgs^* - \xi_{gs} \Mgs^*)^{-1}$,
%as defined previously,
and we applied our spatio-temporal model to estimate the true parameters $\bb_{gs}$, $\rho_{gs}$, $\xi_{gs}$, and $\phi_{gs}$.

%{\it Simulation 2.}
In the second simulation study, we evaluated the relative efficiency of our spatio-temporal model compared to the GLM,
both for voxel-level and region level estimators of secondary parameters of interest
($\th_g$ and $\boldsymbol{\th}_g^*$).
We selected two regions, region A with 100 randomly selected voxels
from BA 24, and region B with 1,000 randomly selected voxels from BA 19. We set the true parameters
for $\rho_{gs}$, $\xi_{gs}$, $\phi_{gs}$ and drew 200 samples from ${\rm MVN}(\bb_{gs}, \phi_{gs}^2 \mathbf{\Psi}_{gs} \mathbf{\Psi}_{gs})$,
where $\bb_{gs}, \phi_{gs}^2$ and $\mathbf{\Psi}_{gs}$ were again obtained from our experimental data.
In addition to estimating the true parameters, we estimated the secondary parameters of interest
$\boldsymbol{\theta}_g^*$ (voxel-level) and $\theta_g$ (region-level), defined to address the stated objective of our study
(see Section \ref{ss:inferences}). We calculated the estimated variances for $\hat{\boldsymbol{\theta}}_g^*$ and $\hat{\theta}_g$,
both applying the GLM and our spatio-temporal model, and compared the voxel-level and region-level relative efficiencies
(spatio-temporal model vs. GLM), averaged over 200 samples in each case.

{\it Simulation results.}
The first simulation study revealed extremely small biases in the estimator of the mean model parameters $\bb_{gs}$ (not shown).
Table \ref{t:table2} shows the simulation results for  the covariance parameters for the
generated data.
The following regions were selected: BA 24 (L), BA 32 (R), BA 19 (L) and BA 48 (R).
The differences between the median values of the estimated parameters and the true values, from 200 generated samples
are shown for each of the regions (1)-(4).
%
%\vspace{2ex}
%\begin{table}
% \centering
% % \def\~{\hphantom{0}}
% \caption[Estimation accuracy of the spatio-temporal model]{Estimation accuracy of the spatio-temporal model obtained from 200 simulated data sets for the cocaine addicts and the healthy control subjects, for four selected regions. The difference between the median of the estimated values from the generated samples and the true value is shown.}
%% Starting value for $\phi^2$ was estimated from the simulated data.}
%\label{t:table2}
%\begin{tabular*}{\textwidth}{@{}l@{\extracolsep{\fill}}c@{\extracolsep{\fill}}c@{\extracolsep{\fill}}c@{\extracolsep{\fill}}c@{\extracolsep{\fill}}c@{\extracolsep{\fill}}c@{\extracolsep{\fill}}}
%% {|c|c|c|c|c|c|c|c|c|}
%\hline
%Group &Estimated  & True  & \multicolumn{4}{c}{Region$^*$ (no. voxels)} \\
%          \cline{4-7} \\ [-30pt]
%& parameter & value & 1 (100) & 2 (500) & 3 (1000) & 4 (2000)  \\
%\hline
%\underline{Addicts} & & & & & & \\
%& $\rho$   & \bf{0.5} & -0.0097	& -0.0001  & -0.0015  & -0.0097\\
%& $ \xi$   & \bf{0.2} & 0.0001	& 0.0001  & -0.0001	& 0.0001\\
%& $\phi^2$ & \bf{1.8} & -0.1592	& -0.1537  & -0.1483 & -0.1501 \\
%\hline
%\underline{Controls} & & & & & & \\
%& $\rho$   & \bf{0.5} & -0.0075	& -0.0011 &	-0.0034	& 0.0009\\
%& $ \xi$   & \bf{0.2} & 0.0012	& 0.0009 & -0.0005	& 0.0000\\
%& $\phi^2$ & \bf{1.8} & -0.1137	& -0.1209	& -0.1196	& -0.1191\\
%\end{tabular*}
%\vspace*{-6pt}
%{\small $^*$ 1=BA 24(L), 2=BA 32(R), 3=BA 19(L), 4=BA 48(R). Corresponding region sizes are 100, 500, 1000, 2000 voxels, respectively.
%Starting value for $\rho$: 0.7, starting value for $\xi$: 0.1.}
%\end{table}


\begin{table}
 \centering
 % \def\~{\hphantom{0}}
 \caption[Estimation accuracy of the spatio-temporal model]
 {Estimation accuracy of the spatio-temporal model obtained from 200 simulated data sets for the cocaine addicts
and the healthy control subjects, for four selected regions. The difference between the median of the estimated values from the generated samples and the true value is shown.}

\vspace{2ex}
% Starting value for $\phi^2$ was estimated from the simulated data.}
\label{t:table2}
{\small
\begin{tabular*}{\textwidth}{@{}l@{\extracolsep{\fill}}c@{\extracolsep{\fill}}c@{\extracolsep{\fill}}c@{\extracolsep{\fill}}c@{\extracolsep{\fill}}c@{\extracolsep{\fill}}c@{\extracolsep{\fill}}}
% {|c|c|c|c|c|c|c|c|c|}
\hline \hline
\multirow{2}{*}{Group} & Est.  & True  & \multicolumn{4}{c}{Region$^*$ (no. voxels)} \\
          \cline{4-7} \\ [-2ex]
& param. & value & 1 (100) & 2 (500) & 3 (1000) & 4 (2000)  \\
\hline
\multicolumn{3}{c}{} & \multicolumn{4}{c}{Difference (s.e.) }\\ [-2ex]
\multirow{3}{*}{\underline{Addicts}} & & & & & & \\ [-1ex]
& $\rho$   & \bf{0.5} & -0.0097 (0.0621) & -0.0001 (0.0646) & -0.0015 (0.0581) & -0.0097 (0.0650)\\
& $ \xi$   & \bf{0.2} &  0.0001 (0.0146) &  0.0001 (0.0063) & -0.0001 (0.0044) &  0.0001 (0.0032)\\
& $\phi^2$ & \bf{1.8} & -0.1592 (0.0490) & -0.1537 (0.0213) & -0.1483 (0.0148) & -0.1501 (0.0120)\\
\hline
\multirow{3}{*}{\underline{Controls}} & & & & & & \\ [-1ex]
& $\rho$   & \bf{0.5} & -0.0075 (0.0531) & -0.0011 (0.0547) & -0.0034 (0.0574) &  0.0009 (0.0517)\\
& $ \xi$   & \bf{0.2} &  0.0012	(0.0111) &  0.0009 (0.0056) & -0.0005 (0.0039) &  0.0000 (0.0028)\\
& $\phi^2$ & \bf{1.8} & -0.1137	(0.0461) & -0.1209 (0.0191) & -0.1196 (0.0140) & -0.1191 (0.0107)\\
\hline \hline
\end{tabular*}
}

{\small $^*$ 1=BA 24(L), 2=BA 32(R), 3=BA 19(L), 4=BA 48(R). Corresponding region sizes are 100, 500, 1000, 2000 voxels, respectively.
Starting value for $\rho$: 0.7, starting value for $\xi$: 0.1.}
\end{table}

The second simulation study revealed similar results for regions A and B, and we only present the region B results for brevity.
To summarize the voxel-level relative efficiency results, we plot the distribution of relative efficiencies across 1,000
voxels from region B (Figure \ref{f:voxel_lev_rel_eff}). We see that the GLM-based voxel-level variance estimates are on
average higher than the corresponding estimates obtained from our spatio-temporal model.
The simulation-based estimate of the region-level relative efficiency is 1.72 for region B.
This demonstrates that, on the region level, our model also leads to more precise estimates compared to GLM, and therefore
may lead increased statistical power.

%The estimation for 100 generated samples (for one group) took 1 hours and 5 minutes.
%----- Conclusion from the simulations results --------
We found that estimates for both $\rho_{gs}$ and $\xi_{gs}$ (spatial and temporal
parameters) are quite accurate, while the between subject variability $\phi^2_{gs}$
was slightly underestimated (roughly 6\%), which likely stems from the well-known
downward bias of ML estimates of variance components \citep{laird:82}.
Restricted maximum likelihood (REML) estimation can be used for our model
and would presumably mitigate the observed bias. We favor ML estimation in our context, however, because it
substantially facilitates computations.

%----------------------
% Figure 3
%----------------------
\begin{figure}[!h]
\begin{center}
\includegraphics[scale=0.7]{Figures/Voxel_Rel_Eff_regB_new.eps}
\caption[The distribution of relative efficiency values from 1,000 voxels
randomly selected from left BA 19]{The distribution of relative efficiency values from 1,000 voxels
randomly selected from left BA 19. Each voxel-level quantity is calculated
as the mean relative efficiency across 200 simulated data sets. The relative efficiency
was calculated based on the estimates from our model, relative to GLM.}
\label{f:voxel_lev_rel_eff}
\end{center}
\end{figure}

%======================================
\section{Discussion}
%======================================
\label{s:discussion}
We propose a novel spatio-temporal modeling framework for functional neuroimaging data, overcoming unsupported
assumptions of independence between the multiple summary statistics for each subject (e.g. from different
sessions or stimuli) and between different brain locations
within neuroanatomic regions.
We used Brodmann area templates, which are popular in the neuroimaging community, to define neuroanatomic regions,
but other maps may be applicable.
Our proposed framework is targeted to fMRI data; however, the general ideas extend easily to PET brain imaging data measuring
distributed patterns of regional cerebral blood flow.
Our approach provides a unified framework for voxel-level and region-level
inferences. An important contribution of our work is that we derive efficient computational solutions, e.g.
enabling the inversion of large covariance matrices, to facilitate implementation.
Our simulation studies demonstrated that our method is quite accurate and that the estimated
standard errors are on average smaller compared to GLM, leading to better precision of estimates of model
parameters, and more accurate statistical tests by borrowing the strength across related measures of neural activity.

One limitation of our model is that it does not account for correlations {\it between} regions.
Doing so within our current framework would lead to substantially increased, perhaps prohibitive, computations.
\citet{bowman:08} establish a Bayesian framework that models between-region correlations, in addition
to within-region correlations, but the number of regions included in the analysis is consequently constrained
by the sample size. One advantage of our proposed spatio-temporal model is that the number of regions included
is unconstrained.
Our current maximum likelihood estimation procedure requires close monitoring to ensure convergence,
confinement within the parameter space, and avoidance of the other specified  boundary conditions.
% possibly comment on starting value ..

In summary, our proposed spatio-temporal model provides an appealing, computationally efficient alternative to standard GLM-based methods for analyzing fMRI data.  Our model is based on assumptions that are more neurophysiologically plausible, capturing correlations between different brain locations and between estimates of neural activity at different scanning sessions.  These correlations lead to interpretive advantages over the GLM, e.g. revealing information about the degree of coherence in brain activity within defined neuroanatomic regions.  Our simulation studies demonstrate that our model estimates are quite accurate and that the standard errors associated with estimates of our mean model parameters are on average smaller than those from a GLM.  This increase in efficiency will often lead to more powerful statistical tests and the detection of more statistically significant voxels.

