Dynamic functional connectivity analysis is an effective way to capture the networks that are functionally associated and continuously changing over the scanning period. However, these methods mostly analyze the dynamic associations across the activation patterns of the spatial networks while assuming that the spatial networks are stationary. Hence, a model that allows for the variability in both domains and reduces the assumptions imposed on the data provides an effective way for extracting spatiotemporal networks. Independent vector analysis (IVA) is a joint blind source separation technique that allows for estimation of spatial and temporal features while successfully preserving variability. However, its performance is affected for higher number of datasets. Hence, we develop an effective two-stage method to extract time-varying spatial and temporal features using IVA, mitigating the problems with higher number of datasets while preserving the variability across subjects and time. The first stage is used to extract reference signals using group-independent component analysis (GICA) that are used in a parameter-tuned constrained IVA framework to estimate time-varying representations of these signals by preserving the variability through tuning the constraint parameter. This approach effectively captures variability across time from a large-scale resting-state fMRI data acquired from healthy controls and patients with schizophrenia and identifies more functionally relevant connections that are significantly different among healthy controls and patients with schizophrenia, compared with the widely used GICA method alone.