Identification of spatiotemporal interactions within/between neuron populations is critical for detection and characterization of large-scale neuronal interactions underlying perception, cognition, and behavior. Univariate analysis has been employed successfully in many neuroimaging studies. However, univariate analysis does not explicitly test for interactions between distributed areas of activity and is not sensitive to distributed responses across the brain. Multivariate analysis can explicitly test for multiple statistical models, including the designed paradigm, and allows for spatial and temporal model detection. Here, we investigate multivariate analysis approaches that take into consideration the 4D (time and space) covariance structure of the data. Principal component analysis (PCA) and independent component analysis (ICA) are two popular multivariate approaches with distinct mathematical constraints. Common difficulties in using these two different decompositions include the following: classification of the revealed components (task-related signal versus noise), overall signal-to-noise sensitivity, and the relatively low computational efficiency (multivariate analysis requires the entire raw data set and more time for model identification analysis). Using both Monte Carlo simulations and empirical data, we derived and tested the generalized partial least squares (gPLS) framework, which can incorporate both PCA and ICA decompositions with computational efficiency. The gPLS method explicitly incorporates the experimental design to simplify the identification of characteristic spatiotemporal patterns. We performed parametric modeling studies of a blocked-design experiment under various conditions, including background noise distribution, sampling rate, and hemodynamic response delay. We used a randomized grouping approach to manipulate the degrees of freedom of PCA and ICA in gPLS to characterize both paradigm coherent and transient brain responses. Simulation data suggest that in the gPLS framework, PCA mostly outperforms ICA as measured by the receiver operating curves (ROCs) in SNR from 0.01 to 100, the hemodynamic response delays from 0 to 3 TR in fMRI, background noise models of Guassian, sub-Gaussian, and super-Gaussian distributions and the number of observations from 5, 10, to 20 in each block of a six-block experiment. Further, due to selective averaging, the gPLS method performs robustly in low signal-to-noise ratio (<1) experiments. We also tested PCA and ICA using PLS in a simulated event-related fMRI data to show their similar detection. Finally, we tested our gPLS approach on empirical fMRI motor data. Using the randomized grouping method, we are able to identify both transient responses and consistent paradigm/model coherent components in the 10-epoch block design motor fMRI experiment. Overall, studies of synthetic and empirical data suggest that PLS analysis, using PCA decomposition, provides a stable and powerful tool for exploration of fMRI/behavior data.