Ranking and averaging independent component analysis by reproducibility (RAICAR)

Abstract

Independent component analysis (ICA) is a data-driven approach that has exhibited great utility for functional magnetic resonance imaging (fMRI). Standard ICA implementations, however, do not provide the number and relative importance of the resulting components. In addition, ICA algorithms utilizing gradient-based optimization give decompositions that are dependent on initialization values, which can lead to dramatically different results. In this work, a new method, RAICAR (Ranking and Averaging Independent Component Analysis by Reproducibility), is introduced to address these issues for spatial ICA applied to fMRI. RAICAR utilizes repeated ICA realizations and relies on the reproducibility between them to rank and select components. Different realizations are aligned based on correlations, leading to aligned components. Each component is ranked and thresholded based on between-realization correlations. Furthermore, different realizations of each aligned component are selectively averaged to generate the final estimate of the given component. Reliability and accuracy of this method are demonstrated with both simulated and experimental fMRI data.

Figure 1

Each submatrix in the cross‐realization correlation matrix (CRCM), R _ij, is the spatial cross‐correlation matrix between realizations i and j. The circled dot represents the global maximum in the CRCM. After finding this maximum, the mth row in each submatrix R _ai and the nth column in each submatrix R _ib are searched for a corresponding submatrix‐specific maximum. These maxima are indicated by dots with their positions given in parentheses.

Figure 2

Simulated spatial sources and their mixing time courses. A: Spatial map of the numbered sources. All the sources are equal in area. B: The corresponding mixing time courses of the sources. The bottom panel shows the added global baseline.

Figure 3

Reproducibility ranking of the simulated data. A: Bimodal distribution of the correlation coefficients. The majority of the correlation coefficients lie in the lower range of 0–0.60, while the remaining fall in the upper range of 0.80–1.00. These two ranges are separated by a broad valley (roughly 0.60–0.80). The solid red curve is the smoothed histogram used to determine the threshold. The arrows indicate the three SCC thresholds used. B: The reproducibility index plots generated using the three SCC thresholds. The half‐maximum cut‐offs are shown with horizontal green lines, indicating six components in each case. The orders of the components derived with all three thresholds are also the same.

Figure 4

Detectability of RAICAR and individual ICA realizations using the sixth simulated source (with lowest SNR). ROC curves of the individual ICA realizations (black) and RAICAR results (red) are shown. The light red region shows the spread of 10 RAICAR repetitions and the red curve shows their mean. The individual ICA realizations exhibit variable results, while the repeated RAICAR results are all virtually identical and outperform the majority of individual ICA realizations.

Figure 5

Reproducibility rankings obtained from simulation 2. A: resting‐state data only. B: resting‐state data with low CNR sources. C: resting‐state data with high CNR sources. The increased CNR level shifts the components towards the left (increasing their rank).

Figure 6

Reproducibility ranking of the event‐related, delayed motor data in one subject. Columns correspond to different sets of ICA realizations. The top row shows the histograms of the SCCs for the three sets of ICA realizations. The correlation coefficients are distributed in two modes, one near zero and the other near one. The bottom three rows are ordered reproducibility index plots for three different thresholds. It can be seen that the number of components passing the cut‐off (shown as horizontal gray lines) do not vary significantly with the SCC thresholds and different sets of ICA realizations.

Figure 7

Several components extracted by RAICAR in the delayed motor dataset. Component maps are displayed on the transparent glass brain with corresponding time courses shown in blue; the red and black curves illustrate the ideal response of “Cue” and “Go” task, respectively. Components 1 and 2 are due to cardiac noise; Component 10 is task‐related activation in the sensory and eye field regions; Component 12 shows the activations corresponding to the visual motion control; Component 13 is related to head motion; Component 14 corresponds to the “Cue” task; Component 18 shows the activations in prefrontal cortex; Component 19 is activated by the “Go” task.

Figure 8

The reproducibility rank obtained from the constant force grip dataset. 17 components were above the reproducibility cut‐off.

Figure 9

The 17 components for the constant force grip dataset. The component maps are shown on the transparent glass brain with the corresponding time courses shown below. The ideal task response is shown (red) for components that are highly correlated with it. The first component may reflect the default mode network; Components 2, 4, and 5 arise from task‐related activations, which include functional areas for motor control and execution; Component 11 seems to be due to task‐related artifacts at the base of the brain. Other components are either artifacts or unexplained brain activation.

Figure 10

Comparison of the component maps extracted by RAICAR and a randomly selected individual ICA realization. The top row of each component shows the RAICAR map and the bottom row shows the individual ICA map. For higher ranked components, both results tend to be similar. While for lower ranked components the individual ICA maps tend to be noisier than the corresponding RAICAR maps.

Figure 11

Convergence of RAICAR. A: The positions of the components as a function of K, generated from simulation 1. The positions of the components do not change when K is larger than 20. B: The area under the ROC curves as a function of K, generated from the lowest SNR source in Simulation 1. The estimation error of the component maps does not change substantially when K is larger than 20. C, D, and E: The average variance of the component maps as a function of the number of ICA realizations (K), generated from the six sources in simulation 1, the 21 sources in the delayed motor dataset, and the 17 sources in the constant force grip dataset. For all three datasets, these curves rapidly approach asymptotic levels. When K = 30 (as used for our reported results), the average variance is no more than 0.005.

Figure 12

Comparison of the estimated number of components for RAICAR and ICASSO. The comparison was conducted using the sources from Simulation 1 with different time series lengths. RAICAR estimates the correct number of sources at all lengths, while ICASSO results vary with data length. When the length exceeds 440, ICASSO estimates could not be obtained due to large memory requirements.