Recently, there is an emerging interest for applications of tensor factorization in big-data analytics and machine learning. Tensor factorization can extract latent features and perform dimension reduction that can facilitate discoveries of new mechanisms hidden in the data. The nonnegative tensor factorization extracts latent features that are naturally sparse and are parts of the data, which makes them easily interpretable. This easy interpretability places the nonnegative factorization as a uniquely suitable method for exploratory data analysis and unsupervised learning. The standard Canonical Polyadic Decomposition (CPD) algorithm for tensor factorization experiences difficulties when applied to tensors with rank deficient factors. For example, the rank deficiency, or linear dependence in the factors, cannot be easily reproduced by the standard PARAFAC because the presence of noise in the real-world data can force the algorithm to extract linearly independent factors. Methods for low-rank approximation and extraction of latent features in the rank deficient case, such as PARALIND family of models, have been successfully developed for general tensors. In this paper, we propose a similar approach for factorization of nonnegative tensors with rank deficiency. Firstly, we determine the minimal nonnegative cones containing the initial tensor and by using a nonnegative Tucker decomposition determine its nonnegative multirank. Secondly, by a nonnegative Tucker decomposition we derived the core tensor and factor matrices corresponding to this multirank. Thirdly, we apply a nonnegative CPD to the derived core tensor to avoid the problematic rank deficiency. Finally, we combine both factorizations to obtain the final CPD factors and demonstrate our approach in several synthetic and real-world examples.