In many areas of science, like neuroscience, genomics and text mining, several important and challenging research questions imply the study of (subject) heterogeneity present in three-way data. In clinical neuroscience, for example, disclosing differences or heterogeneity between subjects in resting state networks (RSNs) underlying multi-subject fMRI data (i.e., time by voxel by subject three-way data) may advance the subtyping of psychiatric and mental diseases. Recently, the Clusterwise Independent Component Analysis (C-ICA) method was proposed that enables the disclosure of heterogeneity between subjects in RSNs that is present in multi-subject rs-fMRI data [1]. Up to now, however, no publicly available software exists that allows to fit C-ICA to empirical data at hand. The goal of this paper, therefore, is to present the CICA R package, which contains the necessary functions to estimate the C-ICA parameters and to interpret and visualize the analysis output. Further, the package also includes functions to select suitable initial values for the C-ICA model parameters and to determine the optimal number of clusters and components for a given empirical data set (i.e., model selection). The use of the main functions of the package is discussed and demonstrated with simulated data. Herewith, the necessary analytical choices that have to be made by the user (e.g., starting values) are explained and showed step by step. The rich functionality of the package is further illustrated by applying C-ICA to empirical rs-fMRI data from a group of Alzheimer patients and elderly control subjects and to multi-country stock market data. Finally, extensions regarding the C-ICA algorithm and procedures for model selection that could be implemented in future releases of the package are discussed.