%=======================================================================================================
\chapter{Evaluating Functional Autocorrelation within Spatially Distributed Neural Processing Networks}
\label{ch:autocorr}
%=======================================================================================================
%======================
\section{Introduction}
%======================
Data-driven statistical approaches, such as cluster analysis or independent component analysis, applied to in vivo functional neuroimaging data help to identify neural processing networks that exhibit similar task-related or resting-state patterns of activity. Ideally, the measured brain activity for voxels within such networks should exhibit high autocorrelation.  An important limitation is that the algorithms do not typically quantify or statistically test the strength or nature of the within-network relatedness between voxels.

There have been numerous applications of clustering analysis to functional imaging data \citep{bowman:04b}.
Despite the successful application of clustering methods to neuroimaging data, statistical procedures for testing the validity of the resulting cluster structure are still lacking.  Although a clustering solution suggests the existence of networks of functionally related voxels, the clusters exhibit unknown levels of autocorrelation, which is a desirable characteristic of behaviorally relevant neural circuits.  Typical cluster analyses neither quantify the strength or nature of the relatedness between within-cluster voxels using an easily interpretable measure nor apply formal tests of hypotheses to evaluate the statistical significance of the level of autocorrelation.  These limitations are particularly noteworthy since clustering procedures always produce a solution, even in the absence of true functional connections in the data.  Therefore, a valuable addition to cluster analyses applied in functional neuroimaging studies is to measure the degree of functional autocorrelation within a clustering solution and to evaluate the statistical significance of the observed associations.

%
%Some clustering criteria, and stopping rules for hierarchical clustering methods, quantify measures based on within-cluster similarity, e.g. proportional to the change in within-cluster variation, possibly relative to between-cluster variation.  However, these quantities are limited
%because they are relative measures, e.g comparing two consecutive levels of a clustering hierarchy.  Also, these measures often lack easy interpretability and are not amenable to formal hypothesis testing.  The probability distributions of clustering criteria are difficult to derive theoretically \citep{hartigan:77,hartigan:78} and thus do not provide a basis for formal hypothesis testing using a conventional Neyman-Pearson framework.  Furthermore, even if one successfully derives the distribution of a particular clustering criterion, the result may be somewhat limited in practice in cases where an alternative clustering algorithm provides better performance.
%%Rather than focusing on the distribution of a single clustering criterion, we take a more general approach of evaluating statistical significance of %the functional autocorrelation present in the final clustering solution, regardless of the clustering algorithm employed to obtain the solution.

Similarly, even though ICA has proved promising, there is a need to study the properties of ICA as applied to fMRI data. One major issue in application of ICA is that the reliability of the estimated independent components (ICs) is unknown. The results from ICA may vary considerably with different algorithm starting points, sampling of subjects, preprocessing steps or type of decomposition algorithm \citep{calhoun:03}. Therefore, validation of estimated ICs has become important for correct interpretation of ICA results.
%%We will address the issue of evaluating a global functional autocorrelation of a group ICA result.
%\cite{esposito:02} examined the  `infomax' and the fixed-point based algorithms and found (based both on simulated and real data) that both algorithms perform well in the sense of identifying independent components that predict true or predicted patterns, but neither algorithm dominates the other.
%Log-likelihood \citep{mckeown:98} ``was used as a measure of how significantly the estimated independent sources fit the statistical structure of the data and receiver operating characteristics (ROC) and linear correlation analyses were used to evaluate the algorithms' accuracy of estimating the spatial layout and the temporal dynamics of simulated and real activations \ldots'' \cite{himberg:04} evaluate statistical reproducibility and
%consistency of ICA results.
There has been previous work related to evaluating the results of ICA by testing mutual independence between the
extracted source signals \citep{murata:01,shimizu:01,chiu:03,stogbauer:04,wu:09}.
There have also been several measures proposed to evaluate the fit of the estimated ICs to the data \citep{mckeown:98,esposito:02,himberg:04}
and some attempts to investigate the reliability of the ICA estimates \citep{himberg:03,himberg:04}.
Still, it may remain unclear how much autocorrelation is present in the identified components.

We introduce a global index to evaluate the validity of neural processing networks identified by functional clustering or by ICA.  Specifically, we propose the use of Moran's $I$ statistic \citep{cliff:73} to quantify and test the autocorrelation present within computed functional networks.  We tailor our definition of the global index to neuroimaging applications through the use of network-based neighborhoods, and we compute network-specific contributions to the overall autocorrelation.  Consequently, our methods target autocorrelation exhibited by voxels within defined neural processing networks, rather than between voxels globally throughout the entire brain.  Thus, we make no assumptions (or statements) about statistical dependence between voxels on a global level.  Hypothesis testing results from our framework augment the descriptive findings of data-driven analysis by determining the presence, direction, and statistical significance of functional autocorrelation within neural processing networks.
For applications of the proposed methods to resting-state fMRI data, we present an explicit temporal component.
We illustrate the use of our methodology with data from two studies: a PET study of regional cerebral blood flow (rCBF) correlates of parametrically manipulated working memory among individuals with schizophrenia and an fMRI resting-state study of depression.  Moreover, we provide empirical support for the use of the proposed methods using a bootstrap analysis and a simulation study.

%===========================
\section{Experimental Data}
%===========================
\noindent
{\it PET Data on Working Memory in Schizophrenia Patients.}
We illustrate cluster analysis using PET data from a study of schizophrenia, but the cluster analysis is also applicable to fMRI data.
Our data are from a parametric study of working memory in an executive function task in individuals with schizophrenia.  We use data from $N=16$ subjects, each having a total of 8 PET scans obtained using the blood flow tracer $[^{15}$O]H$_2$. After aligning and re-slicing the PET images for each subject \citep{woods:98a}, we spatially normalize each subject's scans to a population-representative PET atlas
\citep{woods:98b} centered in Talairach stereotaxic coordinates \citep{talairach:88}. We refrain from spatial smoothing of the PET data during preprocessing to avoid artificially inducing correlations between nearby voxels.  We analyze data from 9,919 voxels, consisting primarily of gray-matter voxels that exhibit at least a 1.6\% change throughout the study. The study design encompasses 2 replicate sessions of 4 scans, where the scans in each session represent four experimental conditions that parametrically vary working memory load-the ability to manipulate information in short--term memory.  Specifically, the conditions are defined as digit shadowing (minimal working memory load) and low, moderate, and high working memory loads in the form of serial addition tasks.  In the digit shadowing condition, the subjects simply repeat a number that they receive by an auditory presentation, which imposes a negligible load on working memory.  In the other three experimental conditions, the subjects receive auditory presentations of a series of positive integers and are instructed to provide the sum of the current number and the preceding number, which requires that they store the previous number in short-term memory and mentally suppress the previous sum.  The range of the integers in the working memory load presentations distinguishes the conditions.  Both numbers in the low load condition fall between 1 and 3, both numbers in the moderate load condition are between 1 and 5, and the high load condition includes integers between 1 and 9.  On average, the sums for the three active load conditions during our experiment are 4.23 (low), 7.22 (moderate), and 10.09 (high).

{\it FMRI Data from a Study on Depression.}
We use data from a second study intended to characterize the impact of childhood abuse/early trauma, and its developmental course, on brain structure and connectivity. For our purpose, we analyze fMRI resting-state scans of seven women with a history of major depression. While in the scanner, subjects had their eyes open and were looking at a visual fixation cross. Functional data were acquired on a Siemens 3T HRRT (high resolution research tomography) scanner, in 3D mode. During a 7.5 min acquisition, a series of 210 scans were acquired with TR=2sec, 20 axial slices, 3.4 $\times$ 3.4 $\times$ 4 mm$^3$ resolution. The functional runs were collected with a Z-saga sequence to avoid orbitofrontal signal ablation.
A group ICA was applied to investigate the resting-state functional connectivity of the depressed patients. There is substantial interest in determining associations between the brain activity characteristics of different regions while subjects are in a resting state. For example, a set of regions has been consistently identified in these investigations and has been labeled as the {\it default mode network} (DMN).

The first step in the analysis (preprocessing and the group ICA) of the fMRI data was done in FSL (FMRIB's Software Library, {\tt www.fmrib.ox.ac.uk/fsl}) software package. The preprocessing and the first level, individual subject  analysis of the PET data was done in SPM (Statistical Parametric Mapping, Welcome Department of Cognitive Neurology, London, UK, {\tt www.fil.ion.ucl.ac.uk/spm/}). The rest of the analyses were carried out in Matlab.

\section{Methods}

\noindent{\it Notation and terminology.} The framework we present allows estimation (and testing) for the presence of autocorrelation separately for each experimental condition $p$ ($p=1,\ldots,P)$  or scan $s$ ($s = 1,\ldots,S$). Data from a single scan are represented as a 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$.  Each voxel contains intensity information corresponding to measures of localized brain activity for fMRI data, or rCBF in case of PET data.
%As will be described below, ICA spatially partitions the brain into IC's, and similarly cluster analysis partitions the brain into
%functional clusters, where both may be regarded as internally homogeneous {\it neural processing networks}.

%=======================
\subsection{Clustering}
%=======================
Summary statistics, such as means or regression coefficients, are often inputs for clustering procedures applied to task-related fMRI or PET data \citep{simon:04}.  We develop notation for this case, but the methodology extends to settings that cluster the original response data directly, e.g. in resting-state fMRI studies.  We denote rCBF, as a proxy for brain activity, by $\bY_k(v)=\bigl( Y_{k1}(v),\ldots,Y_{kS}(v) \bigr)'$ , representing localized blood flow measurements from all $S$ scans for subject $k$, $k=1,\ldots,K$
(at voxel $v$).  For clustering, we typically consider statistics $\mathbf{T}(v)=\bigl(T_1(v),\ldots,T_P(v) \bigr)'=f\bigl(\bY_1(v),\ldots,\bY_K(v)\bigr)$, $P\leq S$, that summarize data from all individuals, e.g. where $T_p(v)$ is the mean or estimated effect associated with experimental condition $p$.  By selecting the $p$th element of the summary vector from every voxel, we also define the vector
$\mathbf{T}_p=\bigl(T_p(1),\ldots, T_p(V)\bigr)'$, which will facilitate our upcoming discussion of functional autocorrelation.
	
Clustering procedures classify the $V$ voxels in an image into $G$ groups, with each cluster consisting of $V_g$ voxels, where $g=1,\ldots, G$ and
$V=\sum_{g=1}^G V_g$. Most methods classify voxels by measuring the distance (dissimilarity) between the activity time courses for every pair of voxels and combining voxels with small distances.  For example, one popular measure of distance between the activity in voxels $i$ and $j$ is the Euclidean distance given by
%-----------------------------
\begin{equation}\label{dist}
d\bigl(\bT(v_i),\bT(v_j)\bigr)\!=\! \bigl[\bigl(\bT(v_i)-\bT(v_j)\bigr)' \bigl(\bT(v_i)-\bT(v_j)\bigr)\bigr]^{1/2}.
\end{equation}
Descriptions of particular clustering algorithms are available in \cite{rencher:02}, and \cite{bowman:04b} highlight several of these algorithms for neuroimaging applications.
%----------------------------

Some clustering criteria and stopping rules for hierarchical clustering methods quantify measures based on within-cluster similarity, e.g. proportional to the change in within-cluster variation, relative to between-cluster variation.  However, these quantities are limited
because they are relative measures, e.g comparing two consecutive levels of a clustering hierarchy.  Also, these measures often lack easy interpretability and are not amenable to formal hypothesis testing.  The probability distributions of clustering criteria are difficult to derive theoretically \citep{hartigan:77,hartigan:78} and thus do not provide a basis for formal hypothesis testing using a conventional Neyman-Pearson framework.  Furthermore, even if one successfully derives the distribution of a particular clustering criterion, the result may be somewhat limited in practice in cases where an alternative clustering algorithm provides better performance.
We evaluate statistical significance of the functional autocorrelation present in the final clustering solution, regardless of the clustering
algorithm employed to obtain the solution and without relying on specific distributional assumptions of the data.

%====================
\subsubsection{ICA}
%====================
\noindent{\it Classical ICA.}
Let $\bY$ denote an $S \times V$ matrix of the observed voxel time courses, where $S$ is the total number of scans, and $V$ is the number of voxels included in the analysis; $\mathbf{C}$ is an $N \times V$ random matrix whose rows $C_i$ are to be filled with the unknown realizations of the $N$ spatial components (images, reshaped in 1-D vectors), and $\mathbf{A}$ is an $S \times N$ mixing matrix, whose columns contain the associated time-courses of the $N$ components. The spatial ICA problem for fMRI time series can be formulated as an estimation of the following model:
%-----------------------------
\begin{equation}\label{ICA1}
\bY =\bA \bC \, .
\end{equation}
%-----------------------------

There are no assumptions about the mixing matrix $\bA$, while the constraint on the spatial processes $C_i$ is that they are (ideally) mutually statistically independent. The amount of statistical dependence within a fixed number of spatial components can be quantified by means of their mutual information, an important function in information theory \citep{common:94}. Thus, the ICA decomposition of $\bY$ can be defined (up to a multiplicative constant and to the sign) as an invertible transformation: $\bC = \bU \bY$, where the matrix $\bU$ (so-called unmixing matrix) is determined such that the mutual information of the target components $C_i$ is minimized (i.e., such that $C_i$'s are ``as independent as possible''). The matrix $\bA$ is the (pseudo)inverse of $\bU$. In \cite{mckeown:98}, the sources are estimated by iteratively optimizing the unmixing matrix $\bU$ so that $\bC = \bU \bY$ contains mutually independent rows, using the `infomax' algorithm.

\noindent{\it Probabilistic ICA (PICA).}
The model of equation \eqref{ICA1}  does not include random noise. The PICA model, which extends model \eqref{ICA1}, is formulated as a generative linear latent variables model. The model ``is characterized by assuming that the $S$-variate vector of observations is generated from a set of $q$ statistically independent non-Gaussian sources via a linear instantaneous mixing process corrupted by additive noise $\boldsymbol{\eta}(t)$'' \citep{beckmann:04}:
%-----------------------------
$$
\mathbf{y}_v=\mathbf{A c}_v + \boldsymbol{\mu} + \boldsymbol{\eta}_v, \quad \forall v\in \{1,\ldots,V\}\, .
$$
%-----------------------------
where $\mathbf{y}_v$ denotes the $S$-dimensional column vector of individual measurements at voxel location $v$, $\mathbf{c}_v$ denotes the $q$-dimensional column vector  of non-Gaussian source signals contained in the data and $\boldsymbol{\eta}_v$ denotes Gaussian noise
$\boldsymbol{\eta}_v \sim N(0,\si^2\boldsymbol{\Si}_v)$.
The number of source processes is assumed to be smaller than the number of observations in time ($q < S$).
The model is closely related to factor analysis.
%Summarize briefly ML estimation of the mixing matrix and sigma2, model order selection, estimation of the unmixing matrix  and inference] In brief,
The individual component maps are obtained using a modified fixed-point iteration scheme to
optimize for non-Gaussian source estimates via maximizing the negentropy. The maps are then separately
transformed to spatial $Z$-scores. %In contrast to raw IC estimates, which only encode the estimated signal,
These $Z$-score maps depend on the amount of variability explained by the entire decomposition at each voxel
location relative to the residual noise. % similar to statistical parametric maps from a GLM analysis.
Next, Gaussian/Gamma mixture models are fitted to the individual $Z$-maps in order to infer voxel
locations that are significantly modulated by the associated time-course.
For more details on PICA, we refer the reader to \cite{beckmann:04}.

\noindent{\it Group ICA}.
In \cite{beckmann:05}, the single-session PICA has been extended to higher dimensions allowing for a model-free analysis of multi-subject or multi-session fMRI data. For group ICA, either tensorial ICA (where the data is decomposed into spatial maps, time courses and subject/session modes), or a simpler temporal concatenation approach (a single 2D ICA is run on the concatenated data matrix obtained by stacking all 2D data matrices of every single data set on top of each other) are possible. The latter approach is recommended to use when the interest is in finding a common spatial pattern, but we cannot assume that the associated temporal response is consistent between subjects, such as in the analysis of data acquired without stimulation (e.g. resting-state data). We hence apply this approach for IC analysis of our fMRI data set.
The technique is derived from parallel factor analysis \citep{harshman:70,harshman:94}.
%[Briefly summarize temporal and tensorial approach if enough space?].
\footnote{The research by Beckmann and Smith described above has been implemented as MELODIC (Multivariate Exploratory Linear Optimized Decomposition into Independent Components - a stand alone C++ program). It is freely available as part of FSL (FMRIB's Software Library).}

We define $\bY(v)=\bigl(\bY_1(v)',\ldots,\bY_K(v)' \bigr)'$, representing concatenated
measures of localized brain activity at voxel $v$, across all subjects. In our upcoming fMRI example, $\bY(v)$ is a $(210*7) \times 1$ vector.
We also define the mean vector of brain activity, across all voxels in all ICs as
$\overline{\bT}=\tfrac 1V \sum_{v=1}^V\bY(v)$.


\subsection{Functional Autocorrelation Statistic}
Functional autocorrelation measures the extent to which voxels within neural processing networks
exhibit similar patterns of brain activity. For any selected element of the summary vector
(indexed by $p$), we measure the functional autocorrelation of a clustering, or ICA solution using Moran's
$I$ \citep{moran:50}, defined as
%-----------------------------
\begin{equation}\label{moransI}
I_p=\frac{V}{S_0}\frac{\bZ_p' \bW \bZ_p}{\bZ_p' \bZ_p} \,.
\end{equation}
%-----------------------------
%Our notation differs somewhat when calculating \eqref{moransI} for a clustering result from that for an ICA results.
In the case of clustering, $\bZ_p=[\bT_p-(\overline{T}_p \b1_V)]$, with $\b1_V$ denoting a $(V\times 1)$ vector of ones
and $\overline{T}_p$ denoting the mean summary element across all voxels
(i.e., $\overline{T}_p=\frac 1V \sum_{v=1}^V T_p(v)$);
$S_0=\sum_{g=1}^G (V_g-1)V_g$ and $\bW$ is the symmetric spatial
connectivity matrix representing the clustering solution. The matrix $\bW$ defines all pairs of voxels
that are in the same cluster, i.e., that exhibit functional connectivity, by including a nonzero value in
the appropriate off-diagonal element. Specifically, we define $\bW$ such that the diagonal elements
are $W_{vv}=0$ , and the $(u,v)$th off-diagonal element is either $W_{vu}=1$, if voxels $u$ and $v$ are in the
same cluster, or $W_{vu}=0$, otherwise.

When applied to group ICA results, we define $\bZ=[\bigl(\bY(1)', \ldots, \bY(V)'\bigr)'-
(\mathbf{1}_V \otimes \overline{\bT})]$, where $\otimes$ denotes the Kronecker product, or the direct product
operator. Here we drop the index $p$ for convenience, since in
our fMRI data example there is only one experimental condition.
$\bW$ is now a block matrix: $\bW=\bW^1 \otimes \bI_{(K*S)\times(K*S)}$, where $\bW^1$ is a $V \times V$ connectivity matrix
(i.e. $W^1_{vv}=0$ and $W^1_{vu}=1$, if voxels $u$ and $v$ are in the same IC and 0 otherwise).
$V$ and $S_0$ are the same as defined above.
Connectivity matrices in other areas of statistical application, e.g. geostatistics, typically
specify connections between locations that are physically adjacent, e.g. neighboring cities. We
adapt Moran's $I$ by tailoring our definition of adjacency to reflect the distributed neural
processing networks, which typically contain voxels that are not all spatially contiguous.
Since our application defines the proximity of voxels based on measures of
brain function, rather than physical or spatial distance, we refer to Moran's $I$ as a {\it functional}
autocorrelation index. However, the phrase {\it spatial} autocorrelation is also appropriate since we
view functional networks as spatially dissociable regions on a brain map.

The massive amount of data collected in functional neuroimaging studies yields high
dimensional spatial connectivity matrices, which may lead to computational difficulties.
In the case of clustering (and similarly for ICA), we can
re-express $I_p$, using the following, perhaps more computationally convenient, expression
%-----------------------------
\begin{equation}\label{moransI_2}
I_p=\frac{V}{S_0}\frac{\sum\limits_{v=1}^V \sum\limits_{u=1}^V W_{vu}
[T_p(v)-\overline{T}_p][T_p(u)-\overline{T}_p]}{\sum\limits_{v=1}^V[T_p(v)-\overline{T}_p]^2} \, .
\end{equation}
%-----------------------------
A similar formula can be obtained when \eqref{moransI} is applied to an ICA solution from fMRI data.
In \eqref{moransI_2}, the numerator of Moran's $I$ calculates the products of the mean-centered summary statistics for
within-cluster voxel pairs and computes the sum of all such products. The denominator gives a
measure of variation about the overall mean. These interpretations reveal the conceptual
similarity of Moran's $I$ to the Pearson correlation coefficient, with Moran's $I$ representing a
spatially weighted version of Pearson's measure. Moran's $I$ generally ranges between -1 and 1,
but the precise limits depend on the eigenvalues of $\bW$.
A large positive value of $I_p$ indicates the presence of functional clusters in which there is high
within-cluster similarity and a negative value indicates clusters exhibiting dissimilarity between
voxels. Randomly assigning voxels to clusters typically results in uncorrelated rCBF between
the voxels within each group.
%
%Similarly, when applied to the ICA result, formula  \eqref{moransI} can be re-expressed as:
%\begin{equation}\label{moransI_3}
%I=\frac{V}{S_0}\frac{\sum\limits_{v=1}^V \sum\limits_{u=1}^V W_{vu} [\bY(v)-\bar{\bT}]'[\bY(u)-\bar{\bT}]}
%{\sum\limits_{v=1}^V[\bY(v)-\bar{\bT}]'[\bY(v)-\bar{\bT}]}
%\end{equation}

\subsection{Hypothesis Testing}
The randomization scheme, characterized by the random allocation of voxels to neural processing networks, serves
as the basis for establishing the null hypothesis and related distributional properties for
evaluating the statistical significance of any observed value of Moran's $I$. To conduct a test to
determine if statistically significant functional autocorrelation is present in a clustering solution,
compute the test statistic
%-----------------------------
\begin{equation}\label{test_stat}
Z_p=\frac{I_p-E(I_p)}{\sqrt{V(I_p)}} \, .
\end{equation}
%-----------------------------
where $E(I_p)$ and $V(I_p)$ represent the theoretical mean and variance of Moran's $I$ under
randomization, respectively (we provide expressions below). Next, calculate
Pr($|Z_p|>z_{\al/2}$), the probability that the magnitude of the observed statistic $Z_p$ exceeds the critical value from a standard normal distribution corresponding to type-I error at level $\al$.
The test rejects the null hypothesis of uncorrelated networks, e.g. no true cluster structure, when the absolute value of $Z_p$ is sufficiently large. Therefore, a statistically significant test result indicates that the observed functional autocorrelation exhibited by a clustering or ICA solution is large relative to the expected value of the autocorrelation under randomization \citep{goodchild:86}.

Calculating the test statistic $Z_p$ requires expressions for the expected value and the variance of Moran's $I$ under the random assignment of voxels to functional networks.  The mean and the variance of Moran's $I$ under randomization are given by following formulas \citep{moran:48a,moran:50,cliff:73,cliff:81}:
%-----------------------------
\begin{equation}\label{mean}
\text{E}(I_p)=-\frac{1}{V-1}\, ,
\end{equation}
%-----------------------------
%and the variance is
%-----------------------------
\begin{equation}\label{variance}
\begin{split}
\text{V}(I_p) &=\frac{V\bigl[(V^2-3V+3)S_1-VS_2+3S_0^2 \bigr]}{(V-1)(V-2)(V-3)S_0^2} \\
& - \frac{m_p\bigl[(V(V-1)S_1-2VS_2+6S_0^2 \bigr]}{(V-1)(V-2)(V-3)S_0^2} - \text{E}^2(I_p)
\end{split}
\end{equation}
%-----------------------------
%
where $S_1=2S_0$, $S_2=\sum_{g=1}^G 4V_g(V_g-1)^2$, and
%-----------------------------
\begin{equation}\label{m_p}
m_p=\frac{V^{-1}\sum_{v=1}^V \bigl[T_p(v)-\overline{T}_p\bigr]^4}{\bigl\{V^{-1} \sum_{v=1}^V \bigl[ T_p(v)-\overline{T}_p \bigr]^2 \bigr\}^2}\,.
\end{equation}
%-----------------------------
%(Formula \eqref{m_p} is written here for case of clustering; it is easy to derive the corresponding
%formula for the ICA case.)

Using Monte Carlo simulations based on our PET data, we empirically validate the theoretical properties of Moran's $I$ under randomization to protect against potential bias, e.g. caused by intrinsic spatial autocorrelation in PET data (see the Results section, Bootstrap analysis).  Substituting the computed quantities \eqref{mean}--\eqref{m_p} into equation \eqref{test_stat} gives the observed value of the test statistic upon which to base hypothesis testing and statistical inferences regarding functional autocorrelation.

\section{Results}\label{s:results_topic2}

\noindent
{\it Working Memory Data.}
We begin by fitting the working memory data using a general linear model with expected value parameters representing the four experimental conditions (digit shadowing and low, moderate, and high working memory loads) and a covariate adjustment for global cerebral blood flow (gCBF).  We obtain ordinary least-squares estimates of the regression parameters and perform clustering using the vector of summary statistics
$\mathbf{T}(v)=\hat{\boldsymbol{\be}}(v)$ for each voxel.  We perform hierarchical clustering using Ward's method \citep{rencher:02} and determine the number of clusters using the cubic clustering criterion (CCC) \citep{searle:83}.  Generally, CCC compares the square of the observed correlation coefficient ($R^2$) and an approximation of its expected value under specified regularity conditions. A plot of CCC ranging from 1 to 40 clusters (Figure~\ref{CCC_plot}) clearly identified $G=29$ as the number of clusters present in the data.  Therefore, the final clustering solution contains 29 clusters, and we depict the cluster map in Figure~\ref{clust_images}(a).


\begin{figure}
\begin{tabular}{cc}
(a) Cluster Map & (b) Cluster 29 \\
\includegraphics[scale=1]{Figures/Fig2a.ps} & \includegraphics[scale=1]{Figures/Fig2b.ps}\\
(c) Cluster 23 & (d) Cluster 20 \\
\includegraphics[scale=1]{Figures/Fig2c.ps} & \includegraphics[scale=1]{Figures/Fig2d.ps}
\end{tabular}
\caption[Cluster maps from Ward's method displayed on axial brain slices ranging from -36mm to +52mm relative to the anterior/posterior commissure plane (z).]{Cluster maps from Ward's method displayed on axial brain slices ranging from -36mm to +52mm relative to the anterior/posterior commissure plane (z). (a) View of the 29 clusters, with each cluster containing voxels that exhibit similar responses across varying loads on working memory.  (b) Cluster 29 includes voxels in the inferior frontal gyrus (BA 47), thalamus, anterior and posterior cingulate (BA 30;32), and the lingual gyrus (BA 17;18).  (c) Cluster 23 consists of voxels in the precentral and middle frontal gyri (BA 9), in the middle, medial, and superior frontal gyri (BA 6), and along the supramarginal gyrus (BA 40).  (d) Cluster 20 spans portions of the right middle frontal gyrus (BA 10), the left medial frontal gyrus (BA 11), and the anterior cingulate (BA 25).}
\label{clust_images}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.8]{Figures/CCC_plot.eps}
\caption[Plotted values of the cubic clustering criterion (CCC)]{Plotted values of the cubic clustering criterion (CCC).  The figure identifies 29 clusters present in the data.}
\label{CCC_plot}
\end{figure}

The cluster sizes range widely with the smallest cluster (1) containing only 8 voxels and the
largest cluster (17) containing 773 voxels.  The average number of voxels in a cluster is 342 with a standard deviation of 165 voxels.  We calculated the mean normalized rCBF values for all voxels within a cluster (averaged across levels of working memory load), and arranged them in increasing order.  Cluster 1 contains voxels with outlying observations and has a mean of only 14.0.  The most active cluster (29) has a mean value of 69.8.  Ten of the clusters identified by Ward's method exhibit brain activity that is lower than the global average of 50, while the other clusters show increased activity relative to the global average value.

Typical cluster analyses cease after establishing a cluster solution and perhaps quantifying associated descriptive statistics.  While the descriptive cluster map provides some insights on functional associations in the brain, it does not quantify the similarity of voxels within the 29 clusters.  Our methodology allows us to compute a global measure of the functional autocorrelation using Moran's $I$ statistic and to evaluate the validity of the clustering solution by examining the statistical significance of the autocorrelation.  Overall, the clusters exhibit very strong positive functional autocorrelation for the digit shadowing condition (0.80), the low working memory load condition (0.81), the moderate load condition (0.81), and the high load condition (0.81).  The large positive values of the functional autocorrelation index suggest that there is generally strong coherence between the voxel profiles within the defined neural processing clusters.  Furthermore, these functional autocorrelations are highly statistically significant with $p$-values all less than 0.0001, confirming the validity of the clustering solution relative to the expected associations under spatial randomness.
%	
%Table \ref{autocorr_PET} lists the global measures of the functional autocorrelation for all four conditions in the working memory PET data.
%\begin{table}
%\centering
%\begin{tabular}{c|cccc}
%\hline\hline \\ [-1ex]
%Condition & digit shadowing & low & moderate & high \\ [1ex]
%\hline \\ [-1ex]
%Functional autocorrelation & 0.80 & 0.81 & 0.81 & 0.81\\ [1ex]
%\hline\hline
%\end{tabular}
%\caption{Global measures of the functional autocorrelation for the 4 experimental conditions in the PET data on working memory in schizophrenia patients.}
%\label{autocorr_PET}
%\end{table}
%Overall, the clusters exhibit very strong positive functional autocorrelation for all four experimental conditions.  The large positive values of the functional autocorrelation index suggest that there is generally strong coherence between the voxel profiles within the defined neural processing clusters.  Furthermore, these functional autocorrelations are highly statistically significant with $p$-values all less than 0.0001, confirming the validity of the clustering solution relative to the expected associations under spatial randomness.

Our analysis identifies clusters containing voxels that, on average, reveal activity profiles resembling particular polynomial curves as a function of varying working memory load levels. The representative mean profiles for three clusters appear in Figure \ref{profiles}; we display the associated standard errors using vertical bars.  The plots reveal potentially important functional characteristics of the relationship of cluster-level activity to working memory.  The nonlinear trends apparent in the plots provide insights on the nature of the stimulus-response relationship in the cluster and may assist in guiding subsequent statistical modeling, particularly in studies of parametrically varying experimental conditions.
\begin{figure}
\begin{tabular}{ccc}
\includegraphics[scale=0.4]{Figures/cl20.eps} &
\includegraphics[scale=0.4]{Figures/cl23.eps} &
\includegraphics[scale=0.4]{Figures/cl29.eps}\\
(a) & (b) & (c)
\end{tabular}
\caption[Mean rCBF profiles, along with standard errors, for 3 representative clusters across levels of working memory load ranging from digit shadowing (DS) to the high load condition.]{Mean rCBF profiles, along with standard errors, for 3 representative clusters across levels of working memory load ranging from digit shadowing (DS) to the high load.  Many clusters (not all shown here) exhibit nonlinear trends over the parametrically varied memory loads.}
\label{profiles}
\end{figure}

Cluster 29 is the most active cluster, i.e. has the highest normalized rCBF, associated with the neural processing underlying working memory.  The voxels within this cluster reveal a roughly quadratic trend in rCBF, exhibiting an increase in measured brain activity from digit shadowing through the moderate working memory load, followed by slight attenuation in the activity at the high load (see Figure~\ref{profiles}(c)).
Axial slices of the cluster map in Figure~\ref{clust_images}(b) show that this neural processing cluster consists of voxels localized bilaterally in the inferior frontal gyrus (Brodmann area (BA) 47) \citep{brodmann:94}, extending into the left superior temporal gyrus (BA 22); thalamus; anterior (BA 32) and posterior (BA 30 and 31) cingulate cortex; and the lingual gyrus (BA 17 and 18).  The quadratic trend that these voxels display suggests that, in individuals with schizophrenia, the activity in this cluster increases to accommodate more substantial memory loads, but then diminishes when the working memory load exceeds capacity.  There was not a precipitous fall in the performance of the serial addition task imposing a high working memory load, rather a near linear decrease in accuracy from the digit shadowing to high load conditions.  Performance at high load levels may be subsumed by other clusters (e.g. cluster 20).

Cluster 20 defines a network of functionally related voxels spanning portions of the prefrontal cortex including the right middle frontal gyrus (BA 10), left ventromedial frontal gyrus (BA 11), and the subcallosal cortex (BA 25), all visible in axial slices $z=-12$ mm through -4 mm of Figure
~\ref{clust_images}(d).  As indicated by Figure~\ref{profiles}(a), these voxels exhibit an approximate quadratic pattern of task-related activity with declining activity from digit shadowing to moderate working memory load and a slight increase at the high load level.

The full set of individual cluster maps identify functional associations related to the neural activity involved in performing parametrically varying working memory tasks.  Here, our detailed inspection of two clusters reveals distinct profiles of task-related brain activity and reveals associations between brain regions that have been previously linked to working memory tasks.  Numerous studies have detected the involvement of both dorsal lateral and ventral lateral regions of the prefrontal cortex in working memory functions \citep{desposito:01}. In addition, several studies have implicated distributed activations in the posterior parietal cortex (BA 40/7), anterior cingulate cortex (BA 32), and the thalamus, as well as premotor areas (BA 6 and 8) \citep{desposito:98,perlstein:03,picard:96}.
%The descriptive results of a cluster analysis further reveal regions with homogeneous task-related activity, possibly coactions, and prompt hypotheses about the neural coding of varying task demands, specific functional connections between regions, and the effect of medications or other treatments on such codes and connections.

The global autocorrelation measures augment the descriptive findings of classification procedures by evaluating the strength and direction of the associations within the defined neural processing clusters.  Local measures of autocorrelation do not extend readily to neuroimaging classification applications because they would require definitions of within-cluster adjacency or functional connections.  To provide a crude measure of the contribution that each cluster has on the overall Moran's $I$ autocorrelation statistic, we compute the percent contribution of each
cluster.
\begin{figure}
\centering
\includegraphics[scale=1]{Figures/Fig5.eps}
\caption{The cluster-specific contributions to the overall autocorrelation index.}
\label{contrib}
\end{figure}
%
Figure~\ref{contrib} displays the results of this local measure.  The plot reveals that several clusters with large or small mean rCBF values contribute more heavily to the global autocorrelation measure.  Generally, a cluster's contribution to the overall level of global autocorrelation gives an indication of the degree of similarity of rCBF profiles from voxels within the cluster and of the spatial extent of the cluster.  Both a high degree of within-cluster homogeneity in working memory-related rCBF values and a large cluster size lead to greater influence on the overall level of functional autocorrelation.  Clusters 3, 4, 5, 23, and 29 have the largest impact on the overall functional coherence, attesting to the high degree of functional autocorrelation in these clusters.


\vspace{2ex}
{\it Depression Data.}
 We first perform a group ICA analysis of the fMRI data using PICA \citep{beckmann:04} as implemented in the FSL tool MELODIC (Multivariate Exploratory Linear Decomposition into Independent Components) Version 3.09. The following data pre-processing was applied to the input data: masking of non-brain voxels; voxel-wise de-meaning of the data; normalisation of the voxel-wise variance.  Pre-processed data were whitened and projected into a 33-dimensional subspace using Principal Component Analysis. The whitened observations were decomposed into sets of vectors which describe signal variation across the temporal domain (time-courses), the session/subject domain and across the spatial
domain (maps) by optimizing for non-Gaussian spatial source distributions using a fixed-point iteration technique \citep{hyvarinen:99}. Estimated component maps were divided by the standard deviation of the residual noise and thresholded by fitting a mixture model to the histogram of intensity values \citep{beckmann:04}. The optimal number of components determined by the PICA algorithm \citep{beckmann:04} is 33.
Many brain voxels were allocated to more than one IC, and some voxels were not assigned to any of the ICs. Out of a total of 19,257 intracranial
voxels, 17,197 were assigned to at least one of the ICs. We assign each of the latter voxels to a single IC, based on the largest (in magnitude)  $Z$-score.
%
\begin{figure}
\centering
(a) \\
\includegraphics[scale=0.45]{Figures/IC4.eps}\\
\includegraphics[scale=0.6]{Figures/t4.eps}\\

(b)\\
\includegraphics[scale=0.45]{Figures/IC6.eps}\\  % IC17 - another possibility
\includegraphics[scale=0.6]{Figures/t6.eps}\\

(c)\\
\includegraphics[scale=0.45]{Figures/IC29.eps}\\
\includegraphics[scale=0.6]{Figures/t29.eps}\\
\caption[ICA results with optimal number of components (33).]{ICA results with optimal number of components (33);
(a) IC 4 (with the corresponding time course): the DMN; (b) IC 6,
(c) IC 29.}
\label{ICs}
\end{figure}

\begin{table}
\centering
\begin{tabular}{rrrc}
\hline \hline
\multirow{2}{*} {IC} & \% of global & \multirow{2}{*}{\# voxels} & \multirow{2}{*}{BA's}\\
   & autocorrelation & & \\
\hline\\
 4 & 12.41 & 843 & 7, 9, 10, 11, 23, 30, 39, 40 \\
 5 & 6.88 & 1032 & putamen, 3, 4, 6, 48 \\
 6 & 7.24 & 968 & thalamus(l), 6, 9, 10, 46 \\
 17 & 10.65 & 1464 & 6, 9, 24, 32, 46\\
 27 & 0 & 45 &  21(r), 22(r), thalamus (r)\\ %4(r),
 29 & 4.24 & 827 & 40(r), 44(r), 45(r), 47\\
 \hline \hline
\end{tabular}
\caption[Network-specific contributions to the overall autocorrelation index
of selected ICs (as identified by ICA) of the fMRI data
from the study on depression.]
{Network-specific contributions to the overall autocorrelation index
of selected ICs (as identified by ICA) of the fMRI data from the study on depression.
{\footnotesize (l and r denote the `left' and `right' hemisphere).}}
\label{table1}
\end{table}


For this data set, the functional autocorrelation index is 0.17. Even though it is modest compared to the one for PET data, it is
statistically significant with a $p$-value less than 0.0001.
%The corresponding crude local measures of autocorrelation for the 33 ICs are all relatively small (between 0 and 0.021).
There are only 4 ICs for which the crude estimates of local autocorrelation are $>0.01$, and Figure~\ref{ICs} shows three of them.
IC 4 represents the DMN, and IC 6 contains voxels in several brain areas (e.g. dorsolateral prefrontal cortex (BA 9), anterior
prefrontal cortex (BA 10), and left thalamus) relevant in depression \citep{mayberg:97, goldapple:04}).
Table~\ref{table1} lists the IC-specific contribution to the global autocorrelation index, number of voxels and Brodmann
regions in which the IC's voxels lie for several ICs.
%Similarly to what we found for PET data, a large IC size leads to greater influence on
%global autocorrelation measure. However, some ICs (e.g. IC 29, see Table~\ref{table1}) with relatively large number of voxels have a small
%contribution to the overall level of functional autocorrelation.

Since for this data set, we calculate the global Moran's $I$ statistic based on the voxels' time courses (not the summary statistic vectors), the statistic $I$ is small, mainly due to the fact that the denominator in \eqref{moransI} increases very fast since the number of time points is large (210). For this reason, we apply another approach to evaluate the global autocorrelation. We calculate the global Moran's $I$ separately for each time point, this way obtaining a $(210*7 \times 1)$ vector of global indices of autocorrelation, for our ICA solution. The obtained values are of a larger scale. They range from 0.03 to 0.53, with the peak values within subjects reaching at least 0.3.
The corresponding IC-specific contributions to the overall autocorrelation index agree with the previous results reported in Table~\ref{table1}.
Namely, the ICs with the largest range of `local' Moran's $I$ statistics are still IC 4 (range of $I$'s [0, 0.19]), IC  5 (range of $I$'s
[0, 0.14]), IC 6 (range of $I$'s: [0, 0.15]), and IC 17 (range of $I$'s [0, 0.26]). Therefore, both approaches bring us to similar conclusions about which of the 33 ICs are most functionally and spatially coherent neural networks. There are several possible reasons why, for this data set, the estimated global autocorrelation index is relatively small. First, the sample size is small for a group analysis (only 7 subjects). This is  resting-state data, so the signal to noise ratio is likely to be lower than for task data (such as our PET data set). Also, even though the ICA
analysis should ideally result in functionally correlated ICs, the approach is based on conceptually different criteria which results in ICs that are as statistically independent as possible.

%Since, for this data set, we calculate the global Moran's $I$ statistic based on the voxels' time courses and not the summary statistic vectors, the statistic $I$ is small. This is mainly due to the fact that the denominator in \eqref{moransI} increases very fast (compared to the numerator) because  of the large number of time points ($S=210$). This motivated us to investigate another approach for evaluating the global autocorrelation within functional neural networks. We calculate the global Moran's $I$ separately for each time point, this way obtaining a $((210*7) \times 1)$ vector of global indices of autocorrelation. The plots of the resulting statistics are given in figure~\ref{voxel_Is}.  The obtained values range from 0.03 to 0.53, with the peak values within subjects reaching at least 0.3.
%
%
%\begin{figure}
%\centering
%\includegraphics[scale=0.8]{Figures/M_I_profiles.eps}
%\caption[Global $I$ statistics calculated separately for each time point, and for each subject]{Global $I$ statistics calculated separately for each time point, and for each subject. The black line represents the mean profile across the 7 subjects.}
%\label{voxel_Is}
%\end{figure}
%The corresponding IC-specific contributions to the overall autocorrelation index agree with the previous result reported in Table~\ref{table1}.
%Namely, the IC's with the largest range of `local' Moran's $I$ statistics are still IC \#4 (range of $I$'s: [0, 0.19]), IC \# 5 (range of $I$'s:
%[0, 0.14]), IC \#6 (range of $I$'s: [0, 0.15]), and IC \# 17 (range of $I$'s: [0, 0.26]). Hence, both approaches bring us to a similar conclusion about which of the 33 IC's are the most functionally and spatially coherent neural networks. There are several possible reasons why, for this data set, the estimated global autocorrelation index is relatively small. First, the sample size is small (only 7 subjects). Since this is a resting-state data, the signal to noise ratio is moderate compared to the task data (such as our PET data set). Also, even though the ICA
%analysis should result in functionally correlated IC's, the approach is based on conceptually different criteria from functional clustering. While clustering tries to assign functionally correlated voxels to same clusters, the IC's are determined so that they are as statistically independent from each other as possible.

%\section{Simulations}
%label{ss:simulations}
\noindent
{\it Bootstrap Analysis}
%
\begin{table*}
\centering
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}c|rrrr}%{c|rrrr}
\hline \hline
Mean& \multicolumn{4}{c}{Working memory load} \\
Moran's $I$ & digit shadowing (DS) & low & moderate & high \\
\hline
Theoretical & -0.000101 & -0.000101 & -0.000101 & -0.000101 \\
Bootstrap & -0.00033 & -0.00025 & - 0.00029 & -0.00030\\
%Bias & $-0.36\times 10^{-5} $ & $-0.22\times 10^{-5} $ &
%$-0.51\times 10^{-5} $ & $-0.45\times 10^{-5} $ \\
Bias & $-2.29\times 10^{-4} $ & $-1.49\times 10^{-4} $ &
$-1.89\times 10^{-4} $ & $-1.99\times 10^{-4} $ \\
Variability & $ 0.62\times 10^{-3}$ & $0.62\times 10^{-3}$ & $0.63\times 10^{-3}$ & $0.63\times 10^{-3}$\\
\hline \hline
\end{tabular*}
\caption[Summary of the Bootstrap analysis designed to evaluate a potential bias of Moran's $I$ ]{Summary of the Bootstrap analysis designed to evaluate a potential bias of Moran's $I$ (e.g. caused by intrinsic spatial autocorrelation) in PET data.}
\label{t:bias_PET}
\end{table*}
We do not perform spatial smoothing of the PET data prior to conducting the cluster analysis to avoid artificially inflating the correlations between the rCBF measurements in nearby voxels.  A related issue is whether the intrinsic spatial correlations in PET data, e.g. stemming from the acquisition and reconstruction processes, may bias Moran's $I$ toward higher values.  If true, then the upward bias of Moran's $I$ would call for modifications to the hypothesis testing framework outlined previously.  We conduct a Monte Carlo simulation study to evaluate the expected value and the variability of Moran's $I$ under the null hypothesis of random cluster assignments.

%Bootstrapping is a well-established procedure for empirically constructing the distribution of an estimated quantity - in our case, Moran's $I$ statistic \citep{efron:86}.
We conduct a bootstrap analysis that includes 500 samples drawn, with replacement, from the working memory data.  For each sample, we compute summary statistics representing the four working memory load conditions, adjusted for gCBF.  To establish the null distribution of Moran's $I$, we randomly allocate voxels to clusters in each bootstrap sample, with the cluster sizes matching those obtained in our original analysis, and then we calculate Moran's $I$ for the resulting randomized clustering solutions.  The 500 computed statistics represent the empirical distribution of Moran's $I$ under the null hypothesis.  The bootstrap distribution allows a comparison of empirical and theortical properties of Moran's $I$ including both its mean and variance.


Table \ref{t:bias_PET} summarizes the results of the bootstrap analysis. Theoretically, the mean of Moran's $I$ is equal to -0.000101.
%Based on the bootstrap null distribution of Moran's $I$ for each of the four working memory load conditions, the mean values of Moran's $I$ are -0.00033, -0.00025, -0.00029, and -0.00030  for working memory loads ranging from digit shadow (DS) to the high load, respectively.
The similarity of the theoretical and empirical mean values provides strong evidence of negligible bias in Moran's $I$ for our PET application.  We find the empirical distributions of the bootstrap biases, reflecting estimates of Moran's $I$ obtained from the boostrap samples (under randomization) minus the theoretical value of the mean.  Averaging over all the bootstrap samples gives the bootstrap estimate of bias, for each condition.  Although slightly negative,
% ($-0.36*10^{-5}\pm 0.62 * 10^{-3}$, $-0.22*10^{-5} \pm 0.62* 10^{-3}$, $-0.51*10^{-5}\pm 0.63* 10^{-3}$, and $-0.45*10^{-5}\pm 0.63* 10^{-3}$),
the biases are all extremely small.  Similarly, the bootstrap variance estimates exhibit close correspondence to the theoretical values.  The magnitudes of all differences between the empirical and theoretical variances are less than $0.4224\times 10^{-6}$.  The accuracy of the theoretical mean and variance expressions relative to their empirical counterparts obtained from the bootstrap estimates supports the use of the hypothesis testing procedure outlined in the Methods section.
%\begin{figure*}
%\begin{tabular}{cc}
%\includegraphics[scale=0.6]{Figs/Figure6a.eps} & \includegraphics[scale=0.6]{Figs/Figure6b.eps}\\
%\includegraphics[scale=0.8]{Figs/Bootstrap_mean_Intermed.eps} & \includegraphics[scale=0.6]{Figs/Figure6d.eps}
%\end{tabular}
%\caption{Estimated bias of Moran's $I$ statistic, under the null hypothesis of randomization, based on a bootstrap analysis consisting of 500 samples.  The bootstrap analysis indicates that Moran's $I$ exhibits negligible bias for the PET working memory data under the null hypothesis, supporting the use of the inferential procedures described in the paper.}
%\label{bias}
%\end{figure*}


\vspace{2ex}
\noindent
{\it Simulated Data Example.}
We analyze simulated data based on the working memory PET study to illustrate further the utility of Moran's $I$
for detecting the presence of functional autocorrelation in clustering solutions.
We simulate data using a mixture model that takes a weighted average of the global mean and the
mean of the voxel-specific summary values within the same neural processing cluster, with weights defined by a spatial dependence (correlation) parameter $\rho$.  Specifically, we generate data for a particular voxel within cluster $g$ using
\begin{equation}
\hat{T}_v=(1-\rho)\mu + \rho \overline{T}_{N_v}+\ve_v ,
\label{simul_model}
\end{equation}
where $\ve_v\sim N(0,\si_v^2)$, $\mu$ is the overall mean, e.g. 50; and $\overline{T}_{N_v}$ is the mean of the summary statistics from the
$V_g-1$ voxels in the same cluster as voxel $v$, i.e., from the functional {\it neighborhood} of $v$, denoted by $N_v$.
We consider six simulated data sets with the correlation parameter varying across the set of values $\rho \in \{0, 0.2, 0.4, 0.6, 0.8, 1\}$.
We use the 29 clusters computed from the working memory data to define the cluster structure and the corresponding cluster sizes $V_g$.
We utilize estimates of the variance parameters and the summary statistic vector from the working memory data.

The error terms from different voxels in model \eqref{simul_model} are independent.  However, the model yields simulated values that depend on the activity levels from other voxels within the same neural processing cluster, when $\rho>0$.
The activity in a given voxel is uncorrelated with the activity of the other voxels within the same cluster when $\rho=0$.
%As the correlation between the brain activity at a particular voxel and the activity in other voxels within the same cluster approaches its maximum
%($\rho \rightarrow 1$), the simulated value approaches the mean of the functionally related voxels.  The simulation model
%\eqref{simul_model} relates to a simultaneous autoregressive model, for which the Neyman-Pearson test of the null hypothesis $\rho=0$ against the alternative $\rho=\rho_0>0$ reduces to a multiple of Moran's $I$, for normally distributed data \citep{ripley:81}.

We compute Moran's $I$ for the clustering solution corresponding to each simulated data set, and Figure~\ref{autocorr} displays the relationship between Moran's $I$ and the correlation parameter for the simulated data. The plot illustrates that Moran's $I$ is effective for detecting and quantifying functional autocorrelation in clustered data.  The value of Moran's $I$ is zero in the absence of functional autocorrelation.  In practice, the importance of revealing the absence of functional autocorrelations lies in the fact that a simple cluster analysis would yield misleading results by providing a clustering solution when no true underlying structure is present in the data.  Moran's $I$ increases as a function of $\rho$, reinforcing that the global index captures the functional autocorrelation inherent within the computed clusters.  The simulated data do not achieve perfect correlation, even when
$\rho=1$, primarily due to the addition of random noise specified by our simulation model.  Consequently, Moran's $I$ is approximately equal to 0.8 when the correlation parameter for the simulated data equals 1.
%
\begin{figure}
\centering
\includegraphics[scale=0.9]{Figures/Fig6.eps}
\caption{Plot of Moran's $I$ versus the correlation parameter $\rho$ for the simulated data.
Data with increased functional autocorrelation yield larger values of Moran's $I$.}
\label{autocorr}
\end{figure}

%++++= HERE ====
\section{Discussion}
%Here we have appo. of validity of the clustering solution. What we propose hopes to establish
%some baseline criteria to asses that validity. However, in certain situations, one may establish a stronger criterion for validity of the proposed method.
%(``straw man's comparison'')
Data-driven approaches, such as cluster analysis and ICA, are extremely useful tools for describing the functional structure and patterns present in neuroimaging data.  Cluster analyses help dissociate groups of voxels that exhibit similar patterns of brain activity over time or across experimental conditions.  ICA is very useful for detecting resting-state neural networks.
In this work, we extend the interpretations given by those two methods by introducing Moran's $I$ for statistical estimation and hypothesis testing of the functional autocorrelation present in neural processing networks as identified by clustering or ICA.  Applying these procedures helps to evaluate the validity of the typological descriptions given by a clustering or an ICA solution.  Computing Moran's $I$ is fast and easy to program using standard statistical software such as R or SAS.  We use MATLAB to calculate the functional autocorrelation of neural processing networks of (functional clusters or ICA components) from PET working memory data, fMRI depression data, and simulated data.
%Following the cluster analysis, computing Moran's $I$ only took approximately 15 seconds to run on a Pentium 4, 3.2 gigahertz desktop computer with 1 %gigabyte of random access memory.

We propose a definition of Moran's $I$ that estimates the functional autocorrelation separately for each element of the summary statistic vector, e.g. pertaining to a specified experimental condition such as high working memory load. For studies that cluster time courses measured under numerous experimental stimuli, it is perhaps more beneficial to aggregate the spatial index across the components of the summary statistic vector.  One can easily modify Moran's $I$ to define a more crude measure that ``averages'' across all elements of the summary vector.  When computing the functional network-specific contributions to the global autocorrelation index, the network sizes impact the corresponding contributions.  The data must provide extremely strong evidence that a small neural processing network (e.g. a small cluster) contributes substantially to the global autocorrelation for the network to weigh heavily in the global calculation.  The number of networks may also affect estimates of the functional autocorrelation.  Assigning voxels into a large number of groups will allow the solution to reflect strong resemblances among within-cluster voxels, whereas designating fewer groups will often lead to less coherence within clusters.

We demonstrate the use of Moran's $I$ using PET neuroimaging data analyzed using cluster analysis and fMRI resting-state data analyzed using ICA.  In addition, Moran's $I$ may combine with results from related descriptive procedures such as principal component analysis \citep{friston:93}.  Structural equation modeling is a useful approach to examine the direct and indirect effects that one brain region has on another \citep{mcintosh:94}, but it requires the specification of anatomical models and substantially simplifies spatial representations from a large number of voxels to a small number of regions.  Both cluster analyses and ICA generally serve as exploratory tools, but often the ultimate objectives of neuroimaging studies require the use of inferential or confirmatory statistical procedures.  In this light, conducting such analyses, followed by an assessment of the functional autocorrelation, helps provide insights about characteristics of neuroimaging data that may prove valuable for subsequent modeling, estimation, and hypothesis testing.

%
%\vspace{2ex}
%{\bf Future Work.}
%There are several areas where we can expand our research in investigating the measures of the level of autocorrelation within functionally defined neural networks. One of our goals is to evaluate the accuracy of the normal approximation formulas under randomization assumption \citep{cliff:73} derived for neuroimaging data.	In particular, we think of the spatial correlations present in neuroimaging data. Inherent spatial correlations can be present in the data due to neurophysiology, data preprocessing, scanning acquisition etc. There are couple of different scenarios in our evaluation of the appropriateness of these approximations: (i) to consider the data at the voxel level (randomization assumption: randomly select voxels and assign them into clusters) and (ii) consider regional based clustering (with much smaller number of regions, compared to the number of intracranial voxels).
%%than selected throughout the brain.
%Because the number of voxels or regions play the central role in these calculations, we may observe smaller or larger degrees of bias when we consider these two settings.
%Thus far, we have started investigating the first setting (voxel-level analysis) in our bootstrap analysis presented above.
%We want to investigate both settings because we feel the answers may vary.
%%	What gives the possibility: the large number of voxels really drives the story
%%	What shall we do if we find the differences?
%In case we find that the answers vary, we hope to propose an alternative hypothesis testing framework (possibly non-parametric).
%%In conclusion, our future work will be concentrated on the investigation into new testing framework (pending these investigations into appropriateness of the approximation
%
%Furthermore, in our simulation study based on the PET data, the normal approximation seems to be reasonable well, so we felt that there would be little gain by developing alternative hypothesis testing procedure. However, it is an open question whether or not the regional level analysis will need to do some adjustment for inherent spatial correlations in the data.

