Non‐Negative Matrix Factorization (NMF) is a widely used dimension reduction method that factorizes a non‐negative data matrix into two lower dimensional non‐negative matrices: one is the basis or feature matrix which consists of the variables and the other is the coefficients matrix which is the projections of data points to the new basis. The features can be interpreted as sub‐structures of the data. The number of sub‐structures in the feature matrix is also called the rank. This parameter controls the model complexity and is the only tuning parameter for the NMF model. An appropriate rank will extract the key latent features while minimizing the noise from the original data. However due to the large amount of optimization error always present in the NMF computation, the rank selection has been a difficult problem. We develop a novel rank selection method based on hypothesis testing, using a deconvolved bootstrap distribution to assess the significance level accurately. Through simulations, we compare our method with a rank selection method based on hypothesis testing using bootstrap distribution without deconvolution and a method based on cross‐validation; we demonstrate that our method is not only accurate at estimating the true ranks for NMF, especially when the features are hard to distinguish, but also efficient at computation. When applied to real microbiome data (eg, OTU data and functional metagenomic data), our method also shows the ability to extract interpretable subcommunities in the data.