Machine learning (ML) has been extensively applied in brain imaging studies to aid the diagnosis of psychiatric disorders and the selection of potential biomarkers. Due to the high dimensionality of imaging data and heterogeneous subtypes of psychiatric disorders, the reproducibility of ML results in brain imaging studies has drawn increasing attention. The reproducibility in brain imaging has been primarily examined in terms of prediction accuracy. However, achieving high prediction accuracy and discovering relevant features are two separate but related goals. An important yet under-investigated problem is the reproducibility of feature selection in brain imaging studies. We propose a new metric to quantify the reproducibility of neuroimaging feature selection via bootstrapping. We estimate the reproducibility index (R-index) for each feature as the reciprocal coefficient of variation of absolute mean difference across a larger number of bootstrap samples. We then integrate the R-index in regularized classification models as penalty weight. Reproducible features with a larger R-index are assigned smaller penalty weights and thus are more likely to be selected by our proposed models. Both simulated and multimodal neuroimaging data are used to examine the performance of our proposed models. Results show that our proposed R-index models are effective in separating informative features from noise features. Additionally, the proposed models yield similar or higher prediction accuracy than the standard regularized classification models while further reducing coefficient estimation error. Improvements achieved by the proposed models are essential to advance our understanding of the selected brain imaging features as well as their associations with psychiatric disorders.