Dimensionality estimation for optimal detection of functional networks in BOLD fMRI data

Abstract

Estimation of the intrinsic dimensionality of fMRI data is an important part of data analysis that helps to separate the signal of interest from noise. We have studied multiple methods of dimensionality estimation proposed in the literature and used these estimates to select a subset of principal components that was subsequently processed by linear discriminant analysis (LDA). Using simulated multivariate Gaussian data, we show that the dimensionality that optimizes signal detection (in terms of the receiver operating characteristic (ROC) metric) goes through a transition from many dimensions to a single dimension as a function of the signal-to-noise ratio. This transition happens when the loci of activation are organized into a spatial network and the variance of the networked, task-related signals is high enough for the signal to be easily detected in the data. We show that reproducibility of activation maps is a metric that captures this switch in intrinsic dimensionality. Except for reproducibility, all of the methods of dimensionality estimation we considered failed to capture this transition: optimization of Bayesian evidence, minimum description length, supervised and unsupervised LDA prediction, and Stein's unbiased risk estimator. This failure results in sub-optimal ROC performance of LDA in the presence of a spatially distributed network, and may have caused LDA to underperform in many of the reported comparisons in the literature. Using real fMRI data sets, including multi-subject group and within-subject longitudinal analysis we demonstrate the existence of these dimensionality transitions in real data.

Figure 1

The phantom in baseline (left) and activation (right) states. Noise is not displayed.

Figure 2

Scheme of the ROC calculation process.

Figure 3

Median dimensionality estimates in simulations, as calculated by various methods (see legend and text), shown as a function of the relative signal variance, V, defined as the variance of the amplitude of the Gaussian activation blobs relative to the variance of the independent background Gaussian noise added to each voxel. CNR is set to 0.3 for Figure 3A and to 1.0 for Figure 3B. The three panels from left to right in A and B show three levels of correlation, ρ, between Gaussian activation blob amplitudes. Range bars on the first (V=0.1) and last (V=1.6) data points reflect the 25%-75% interquartile distribution range across 500 simulation estimates.

Figure 3

Median dimensionality estimates in simulations, as calculated by various methods (see legend and text), shown as a function of the relative signal variance, V, defined as the variance of the amplitude of the Gaussian activation blobs relative to the variance of the independent background Gaussian noise added to each voxel. CNR is set to 0.3 for Figure 3A and to 1.0 for Figure 3B. The three panels from left to right in A and B show three levels of correlation, ρ, between Gaussian activation blob amplitudes. Range bars on the first (V=0.1) and last (V=1.6) data points reflect the 25%-75% interquartile distribution range across 500 simulation estimates.

Figure 4

Partial ROC area (corresponding to false positive frequency [0, 0.1]) as a function of the relative signal variance, V, calculated for different methods of analysis: linear discriminants (LD, on the principal component subspace, with subspace size selected by various methods reported in Fig. 3), univariate general linear model (GLM), penalized discriminant analysis (PDA) with a ridge penalty, and the first component from principal component analysis (PCA). CNR is set to 0.3 for Figure 4A and to 1 for Figure 4B. The three panels from left to right in A and B show three levels of correlation, ρ, between Gaussian activation blob amplitudes. Error bars reflect the 16 partial ROC areas across the centre voxels of the sixteen Gaussian activation blobs (see Fig. 1).

Figure 4

Partial ROC area (corresponding to false positive frequency [0, 0.1]) as a function of the relative signal variance, V, calculated for different methods of analysis: linear discriminants (LD, on the principal component subspace, with subspace size selected by various methods reported in Fig. 3), univariate general linear model (GLM), penalized discriminant analysis (PDA) with a ridge penalty, and the first component from principal component analysis (PCA). CNR is set to 0.3 for Figure 4A and to 1 for Figure 4B. The three panels from left to right in A and B show three levels of correlation, ρ, between Gaussian activation blob amplitudes. Error bars reflect the 16 partial ROC areas across the centre voxels of the sixteen Gaussian activation blobs (see Fig. 1).

Figure 5

Plot of the first 10 eigenvalues of the covariance matrix of a single data set, for CNR = 0.3 (top) and CNR = 1 (bottom). ρ is set to 0.5, and V varies from 0.1 to 1.6. Eigenvalues are averaged across 500 simulated data sets.

Figure 6

Box plot of median dimensionality estimates estimated for 20 simulated large data sets (N = 1500). CNR is 1, ρ is to 0.5, and V is 1.1.

Figure 7

Scatter plots of prediction accuracy vs. reproducibility of spatial maps, as calculated for a linear discriminant on the principal components subspace of the simulated data. Size of the subspace was varied from K=1 (smallest symbol) to K=40 (largest symbol) components; here we show results for K=1, 2, 4, 7, 10, 15, 22, 30, 40. Different trajectories correspond to different levels of relative signal variance, V, and the three panels correspond to three levels of the coupling between activation blob amplitudes, ρ. CNR was fixed at 1.

Figure 8

Asymptotic relationship between global signal-to-noise ratio (gSNR) and optimal dimensionality. Marker size indicates relative signal variance, V, from 0.1 (small) to 1.6 (large). Three colours encode three different levels of CNR, and spatial correlation is encoded by markers.

Figure 9

Prediction-reproducibility plots, calculated on real data: the aging study (A; 2 age groups and 2 behavioral tasks), and the stroke recovery study (B; 3 stroke patients). Marker size indicates K, size of the PC subspace, with smaller markers corresponding to lower K. The range of K was from 4 to 140 (A) and from 1 to 300 (B).

Figure 9

Prediction-reproducibility plots, calculated on real data: the aging study (A; 2 age groups and 2 behavioral tasks), and the stroke recovery study (B; 3 stroke patients). Marker size indicates K, size of the PC subspace, with smaller markers corresponding to lower K. The range of K was from 4 to 140 (A) and from 1 to 300 (B).

Figure 10

Optimal dimensionality and gSNR in real data: the aging study (A) and the stroke recovery study (B).

Figure 10

Optimal dimensionality and gSNR in real data: the aging study (A) and the stroke recovery study (B).